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
52namespace Intrepid2{
53
54
55 template <typename DT, 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,DT> 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,DT>("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 DT::memory_space(), sideNodeMapHost);
173 Kokkos::deep_copy(sideNodeMap_, sideNodeMapHost);
174 }
175
176 template <typename DT, 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,DT> tempPointViewType;
209 typedef Kokkos::DynRankView<PT,Kokkos::LayoutStride,DT> 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);
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<DT>::setJacobian(subcvJacobian, // C, P, D, D
249 sidePoint, // P, D
250 subcvCoordsNode, // C, N, D
251 subcvCellTopo_);
252
253 CellTools<DT>::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*,DT> mapViewType;
259 const auto sideNodeMap = Kokkos::subview(sideNodeMap_, node, Kokkos::ALL());
260
261 //typedef typename ExecSpace<typename PointViewType::execution_space,DT>::ExecSpaceType ExecSpaceType;
262
263 const auto loopSize = numCells;
264 Kokkos::RangePolicy<typename DT::execution_space,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<DT>::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 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 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.
static void setJacobian(Kokkos::DynRankView< jacobianValueType, jacobianProperties... > jacobian, const Kokkos::DynRankView< pointValueType, pointProperties... > points, const WorksetType worksetCell, const Teuchos::RCP< HGradBasisType > basis, const int startCell=0, const int endCell=-1)
Computes the Jacobian matrix DF of the reference-to-physical frame map F.
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.
CubatureControlVolumeBoundary(const shards::CellTopology cellTopology, const ordinal_type sideIndex)
virtual void getCubature(PointViewType cubPoints, weightViewType cubWeights, PointViewType cellCoords) const override
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
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...
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...