Intrepid2
Intrepid2_CubatureControlVolumeBoundaryDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38 // Mauro Perego (mperego@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 
49 #ifndef __INTREPID2_CUBATURE_CONTROLVOLUME_BOUNDARY_DEF_HPP__
50 #define __INTREPID2_CUBATURE_CONTROLVOLUME_BOUNDARY_DEF_HPP__
51 
52 namespace Intrepid2{
53 
54 
55  template <typename SpT, typename PT, typename WT>
57  CubatureControlVolumeBoundary(const shards::CellTopology cellTopology,
58  const ordinal_type sideIndex) {
59 
60  // define primary cell topology with given one
61  primaryCellTopo_ = cellTopology;
62 
63  // subcell is defined either hex or quad according to dimension
64  const ordinal_type spaceDim = primaryCellTopo_.getDimension();
65  switch (spaceDim) {
66  case 2:
67  subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
68  break;
69  case 3:
70  subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
71  break;
72  }
73 
74  // computation order is always one;
75  degree_ = 1;
76 
77  sideIndex_ = sideIndex;
78 
79  // precompute subcontrol volume side index corresponding to primary cell side
80  const ordinal_type sideDim = spaceDim - 1;
81 
82  const ordinal_type numPrimarySides = primaryCellTopo_.getSubcellCount(sideDim);
83  const ordinal_type numPrimarySideNodes = primaryCellTopo_.getNodeCount(sideDim, sideIndex_);
84 
85  const ordinal_type numSubcvSideNodes = subcvCellTopo_.getNodeCount(sideDim, 0);
86 
87  boundarySidesHost_ = Kokkos::View<ordinal_type**,Kokkos::HostSpace>("CubatureControlVolumeBoundary::boundarySidesHost",
88  numPrimarySides, numSubcvSideNodes);
89 
90  // tabulate node map on boundary side
91  switch (primaryCellTopo_.getKey()) {
92  case shards::Triangle<3>::key:
93  case shards::Quadrilateral<4>::key: {
94  for (ordinal_type i=0;i<numPrimarySides;++i) {
95  boundarySidesHost_(i,0) = 0;
96  boundarySidesHost_(i,1) = 3;
97  }
98  break;
99  }
100  case shards::Hexahedron<8>::key: {
101  // sides 0-3
102  for (ordinal_type i=0;i<4;++i) {
103  boundarySidesHost_(i,0) = 0;
104  boundarySidesHost_(i,1) = 3;
105  boundarySidesHost_(i,2) = 3;
106  boundarySidesHost_(i,3) = 0;
107  }
108 
109  // side 4
110  boundarySidesHost_(4,0) = 4;
111  boundarySidesHost_(4,1) = 4;
112  boundarySidesHost_(4,2) = 4;
113  boundarySidesHost_(4,3) = 4;
114 
115  // side 5
116  boundarySidesHost_(5,0) = 5;
117  boundarySidesHost_(5,1) = 5;
118  boundarySidesHost_(5,2) = 5;
119  boundarySidesHost_(5,3) = 5;
120  break;
121  }
122  case shards::Tetrahedron<4>::key: {
123  boundarySidesHost_(0,0) = 0;
124  boundarySidesHost_(0,1) = 3;
125  boundarySidesHost_(0,2) = 0;
126 
127  boundarySidesHost_(1,0) = 0;
128  boundarySidesHost_(1,1) = 3;
129  boundarySidesHost_(1,2) = 3;
130 
131  boundarySidesHost_(2,0) = 3;
132  boundarySidesHost_(2,1) = 4;
133  boundarySidesHost_(2,2) = 0;
134 
135  boundarySidesHost_(3,0) = 4;
136  boundarySidesHost_(3,1) = 4;
137  boundarySidesHost_(3,2) = 4;
138  break;
139  }
140  default: {
141  INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
142  ">>> ERROR (CubatureControlVolumeBoundary: invalid cell topology.");
143  }
144  }
145 
146  Kokkos::DynRankView<PT,SpT> sideCenterLocal("CubatureControlVolumeBoundary::sideCenterLocal",
147  1, sideDim);
148  // map to reference subcell function relies on uvm; some utility functions in cell tools still need uvm
149  sidePoints_ = Kokkos::DynRankView<PT,SpT>("CubatureControlVolumeBoundary::sidePoints",
150  numPrimarySideNodes, spaceDim);
151 
152  for (ordinal_type i=0;i<numPrimarySideNodes;++i) {
153  const ordinal_type sideOrd = boundarySidesHost_(sideIndex_,i);
154  const auto sideRange = Kokkos::pair<ordinal_type,ordinal_type>(i, i+1);
155  const auto sidePoint = Kokkos::subdynrankview(sidePoints_, sideRange, Kokkos::ALL());
157  sideCenterLocal,
158  sideDim,
159  sideOrd,
160  subcvCellTopo_);
161  }
162 
163  const ordinal_type maxNumNodesPerSide = 10;
164  Kokkos::View<ordinal_type**,Kokkos::HostSpace> sideNodeMapHost("CubatureControlVolumeSide::sideNodeMapHost",
165  numPrimarySideNodes, maxNumNodesPerSide);
166  for (ordinal_type i=0;i<numPrimarySideNodes;++i) {
167  const ordinal_type sideOrd = boundarySidesHost_(sideIndex_,i);
168  sideNodeMapHost(i,0) = subcvCellTopo_.getNodeCount(sideDim, sideOrd);
169  for (ordinal_type j=0;j<sideNodeMapHost(i,0);++j)
170  sideNodeMapHost(i,j+1) = subcvCellTopo_.getNodeMap(sideDim, sideOrd, j);
171  }
172  sideNodeMap_ = Kokkos::create_mirror_view(typename SpT::memory_space(), sideNodeMapHost);
173  Kokkos::deep_copy(sideNodeMap_, sideNodeMapHost);
174  }
175 
176  template <typename SpT, typename PT, typename WT>
177  void
179  getCubature( PointViewType cubPoints,
180  weightViewType cubWeights,
181  PointViewType cellCoords ) const {
182 #ifdef HAVE_INTREPID2_DEBUG
183  INTREPID2_TEST_FOR_EXCEPTION( cubPoints.rank() != 3, std::invalid_argument,
184  ">>> ERROR (CubatureControlVolumeBoundary): cubPoints must have rank 3 (C,P,D).");
185  INTREPID2_TEST_FOR_EXCEPTION( cubWeights.rank() != 2, std::invalid_argument,
186  ">>> ERROR (CubatureControlVolumeBoundary): cubWeights must have rank 2 (C,P).");
187  INTREPID2_TEST_FOR_EXCEPTION( cellCoords.rank() != 3, std::invalid_argument,
188  ">>> ERROR (CubatureControlVolumeBoundary): cellCoords must have rank 3 of (C,P,D).");
189 
190  INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(0) != cellCoords.extent(0) ||
191  cubPoints.extent(0) != cubWeights.extent(0), std::invalid_argument,
192  ">>> ERROR (CubatureControlVolume): cubPoints, cubWeights and cellCoords dimension(0) are not consistent, numCells");
193 
194  {
195  const ordinal_type spaceDim = cellCoords.extent(2);
196  const ordinal_type sideDim = spaceDim - 1;
197  const size_type numPrimarySideNodes = primaryCellTopo_.getNodeCount(sideDim, sideIndex_);
198 
199  INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(1) != numPrimarySideNodes ||
200  cubWeights.extent(1) != numPrimarySideNodes, std::invalid_argument,
201  ">>> ERROR (CubatureControlVolume): cubPoints and cubWeights dimension(1) are not consistent, numPrimarySideNodes");
202  }
203  INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(2) != cellCoords.extent(2) ||
204  static_cast<ordinal_type>(cubPoints.extent(2)) != getDimension(), std::invalid_argument,
205  ">>> ERROR (CubatureControlVolume): cubPoints, cellCoords, this->getDimension() are not consistent, spaceDim.");
206 #endif
207 
208  typedef Kokkos::DynRankView<PT,SpT> tempPointViewType;
209  typedef Kokkos::DynRankView<PT,Kokkos::LayoutStride,SpT> tempPointStrideViewType;
210 
211  // get array dimensions
212  const ordinal_type numCells = cellCoords.extent(0);
213  const ordinal_type numNodesPerCell = cellCoords.extent(1);
214  const ordinal_type spaceDim = cellCoords.extent(2);
215  const ordinal_type sideDim = spaceDim - 1;
216 
217  const ordinal_type numNodesPerSubcv = subcvCellTopo_.getNodeCount();
218  tempPointViewType subcvCoords("CubatureControlVolumeBoundary::subcvCoords",
219  numCells, numNodesPerCell, numNodesPerSubcv, spaceDim);
220  CellTools<SpT>::getSubcvCoords(subcvCoords,
221  cellCoords,
222  primaryCellTopo_);
223 
224  //const auto numPrimarySides = primaryCellTopo_.getSubcellCount(sideDim);
225  const ordinal_type numSubcvPoints = 1;
226 
227  tempPointViewType subcvJacobian("CubatureControlVolumeBoundary::subcvJacobian",
228  numCells, numSubcvPoints, spaceDim, spaceDim);
229 
230  tempPointViewType subcvJacobianDet("CubatureControlVolumeBoundary::subcvJacobianDet",
231  numCells, numSubcvPoints);
232 
233  tempPointViewType weights("CubatureControlVolumeBoundary::subcvWeights",
234  numCells, 1);
235  Kokkos::deep_copy(weights, spaceDim == 2 ? 2.0 : 4.0);
236 
237  tempPointViewType scratch("CubatureControlVolumeBoundary::scratch",
238  numCells*numSubcvPoints*spaceDim*spaceDim);
239 
240  const ordinal_type numPrimarySideNodes = primaryCellTopo_.getNodeCount(sideDim, sideIndex_);
241  for (ordinal_type node=0;node<numPrimarySideNodes;++node) {
242  const auto sideRange = Kokkos::pair<ordinal_type,ordinal_type>(node, node+1);
243  const auto sidePoint = Kokkos::subdynrankview(sidePoints_, sideRange, Kokkos::ALL());
244 
245  const auto idx = primaryCellTopo_.getNodeMap(sideDim, sideIndex_, node);
246  auto subcvCoordsNode = Kokkos::subdynrankview(subcvCoords, Kokkos::ALL(), idx, Kokkos::ALL(), Kokkos::ALL());
247 
248  CellTools<SpT>::setJacobian(subcvJacobian, // C, P, D, D
249  sidePoint, // P, D
250  subcvCoordsNode, // C, N, D
251  subcvCellTopo_);
252 
253  CellTools<SpT>::setJacobianDet(subcvJacobianDet, // C, P
254  subcvJacobian); // C, P, D, D
255 
256  auto cubPointsNode = Kokkos::subdynrankview(cubPoints, Kokkos::ALL(), node, Kokkos::ALL());
257 
258  typedef Kokkos::View<ordinal_type*,SpT> mapViewType;
259  const auto sideNodeMap = Kokkos::subview(sideNodeMap_, node, Kokkos::ALL());
260 
261  typedef typename ExecSpace<typename PointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
262 
263  const auto loopSize = numCells;
264  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
265 
266  // compute centers
268  Kokkos::parallel_for( policy, FunctorType(cubPointsNode,
269  subcvCoordsNode,
270  sideNodeMap) );
271 
272  // compute weights
273  const auto sideOrd = boundarySidesHost_(sideIndex_, node);
274 
275  // cub weights node requires to have rank 2
276  auto cubWeightsNode = Kokkos::subdynrankview(cubWeights, Kokkos::ALL(), sideRange);
277  switch (spaceDim) {
278  case 2: {
279  FunctionSpaceTools<SpT>::computeEdgeMeasure(cubWeightsNode, // rank 2
280  subcvJacobian, // rank 4
281  weights, // rank 2
282  sideOrd,
283  subcvCellTopo_,
284  scratch);
285  break;
286  }
287  case 3: {
289  subcvJacobian,
290  weights,
291  sideOrd,
292  subcvCellTopo_,
293  scratch);
294  break;
295  }
296  }
297  }
298  }
299 }
300 
301 #endif
302 
static void mapToReferenceSubcell(Kokkos::DynRankView< refSubcellPointValueType, refSubcellPointProperties... > refSubcellPoints, const Kokkos::DynRankView< paramPointValueType, paramPointProperties... > paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
virtual void getCubature(PointViewType cubPoints, weightViewType cubWeights, PointViewType cellCoords) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
static void getSubcvCoords(Kokkos::DynRankView< subcvCoordValueType, subcvCoordProperties... > subcvCoords, const Kokkos::DynRankView< cellCoordValueType, cellCoordProperties... > cellCoords, const shards::CellTopology primaryCell)
Computes coordinates of sub-control volumes in each primary cell.
CubatureControlVolumeBoundary(const shards::CellTopology cellTopology, const ordinal_type sideIndex)
static void setJacobianDet(Kokkos::DynRankView< jacobianDetValueType, jacobianDetProperties... > jacobianDet, const Kokkos::DynRankView< jacobianValueType, jacobianProperties... > jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F...
static void computeEdgeMeasure(Kokkos::DynRankView< outputValValueType, outputValProperties... > outputVals, const Kokkos::DynRankView< inputJacValueType, inputJacProperties... > inputJac, const Kokkos::DynRankView< inputWeightValueType, inputWeightProperties... > inputWeights, const int whichEdge, const shards::CellTopology parentCell, const Kokkos::DynRankView< scratchValueType, scratchProperties... > scratch)
Returns the weighted integration measures outVals with dimensions (C,P) used for the computation of e...
static void setJacobian(Kokkos::DynRankView< jacobianValueType, jacobianProperties... > jacobian, const Kokkos::DynRankView< pointValueType, pointProperties... > points, const Kokkos::DynRankView< worksetCellValueType, worksetCellProperties... > worksetCell, const HGradBasisPtrType basis)
Computes the Jacobian matrix DF of the reference-to-physical frame map F.
static void computeFaceMeasure(Kokkos::DynRankView< outputValValueType, outputValProperties... > outputVals, const Kokkos::DynRankView< inputJacValueType, inputJacProperties... > inputJac, const Kokkos::DynRankView< inputWeightValueType, inputWeightPropertes... > inputWeights, const int whichFace, const shards::CellTopology parentCell, const Kokkos::DynRankView< scratchValueType, scratchProperties... > scratch)
Returns the weighted integration measures outputVals with dimensions (C,P) used for the computation o...