Intrepid2
Intrepid2_CubatureControlVolumeBoundaryDef.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
16#ifndef __INTREPID2_CUBATURE_CONTROLVOLUME_BOUNDARY_DEF_HPP__
17#define __INTREPID2_CUBATURE_CONTROLVOLUME_BOUNDARY_DEF_HPP__
18
19namespace Intrepid2{
20
21
22 template <typename DT, typename PT, typename WT>
23 CubatureControlVolumeBoundary<DT,PT,WT>::
24 CubatureControlVolumeBoundary(const shards::CellTopology cellTopology,
25 const ordinal_type sideIndex) {
26
27 // define primary cell topology with given one
28 primaryCellTopo_ = cellTopology;
29
30 // subcell is defined either hex or quad according to dimension
31 const ordinal_type spaceDim = primaryCellTopo_.getDimension();
32 switch (spaceDim) {
33 case 2:
34 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
35 break;
36 case 3:
37 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
38 break;
39 }
40
41 // computation order is always one;
42 degree_ = 1;
43
44 sideIndex_ = sideIndex;
45
46 // precompute subcontrol volume side index corresponding to primary cell side
47 const ordinal_type sideDim = spaceDim - 1;
48
49 const ordinal_type numPrimarySides = primaryCellTopo_.getSubcellCount(sideDim);
50 const ordinal_type numPrimarySideNodes = primaryCellTopo_.getNodeCount(sideDim, sideIndex_);
51
52 const ordinal_type numSubcvSideNodes = subcvCellTopo_.getNodeCount(sideDim, 0);
53
54 boundarySidesHost_ = Kokkos::View<ordinal_type**,Kokkos::HostSpace>("CubatureControlVolumeBoundary::boundarySidesHost",
55 numPrimarySides, numSubcvSideNodes);
56
57 // tabulate node map on boundary side
58 switch (primaryCellTopo_.getKey()) {
59 case shards::Triangle<3>::key:
60 case shards::Quadrilateral<4>::key: {
61 for (ordinal_type i=0;i<numPrimarySides;++i) {
62 boundarySidesHost_(i,0) = 0;
63 boundarySidesHost_(i,1) = 3;
64 }
65 break;
66 }
67 case shards::Hexahedron<8>::key: {
68 // sides 0-3
69 for (ordinal_type i=0;i<4;++i) {
70 boundarySidesHost_(i,0) = 0;
71 boundarySidesHost_(i,1) = 3;
72 boundarySidesHost_(i,2) = 3;
73 boundarySidesHost_(i,3) = 0;
74 }
75
76 // side 4
77 boundarySidesHost_(4,0) = 4;
78 boundarySidesHost_(4,1) = 4;
79 boundarySidesHost_(4,2) = 4;
80 boundarySidesHost_(4,3) = 4;
81
82 // side 5
83 boundarySidesHost_(5,0) = 5;
84 boundarySidesHost_(5,1) = 5;
85 boundarySidesHost_(5,2) = 5;
86 boundarySidesHost_(5,3) = 5;
87 break;
88 }
89 case shards::Tetrahedron<4>::key: {
90 boundarySidesHost_(0,0) = 0;
91 boundarySidesHost_(0,1) = 3;
92 boundarySidesHost_(0,2) = 0;
93
94 boundarySidesHost_(1,0) = 0;
95 boundarySidesHost_(1,1) = 3;
96 boundarySidesHost_(1,2) = 3;
97
98 boundarySidesHost_(2,0) = 3;
99 boundarySidesHost_(2,1) = 4;
100 boundarySidesHost_(2,2) = 0;
101
102 boundarySidesHost_(3,0) = 4;
103 boundarySidesHost_(3,1) = 4;
104 boundarySidesHost_(3,2) = 4;
105 break;
106 }
107 default: {
108 INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
109 ">>> ERROR (CubatureControlVolumeBoundary: invalid cell topology.");
110 }
111 }
112
113 Kokkos::DynRankView<PT,DT> sideCenterLocal("CubatureControlVolumeBoundary::sideCenterLocal",
114 1, sideDim);
115 // map to reference subcell function relies on uvm; some utility functions in cell tools still need uvm
116 sidePoints_ = Kokkos::DynRankView<PT,DT>("CubatureControlVolumeBoundary::sidePoints",
117 numPrimarySideNodes, spaceDim);
118
119 for (ordinal_type i=0;i<numPrimarySideNodes;++i) {
120 const ordinal_type sideOrd = boundarySidesHost_(sideIndex_,i);
121 const auto sideRange = Kokkos::pair<ordinal_type,ordinal_type>(i, i+1);
122 const auto sidePoint = Kokkos::subdynrankview(sidePoints_, sideRange, Kokkos::ALL());
123 CellTools<DT>::mapToReferenceSubcell(sidePoint,
124 sideCenterLocal,
125 sideDim,
126 sideOrd,
127 subcvCellTopo_);
128 }
129
130 const ordinal_type maxNumNodesPerSide = 10;
131 Kokkos::View<ordinal_type**,Kokkos::HostSpace> sideNodeMapHost("CubatureControlVolumeSide::sideNodeMapHost",
132 numPrimarySideNodes, maxNumNodesPerSide);
133 for (ordinal_type i=0;i<numPrimarySideNodes;++i) {
134 const ordinal_type sideOrd = boundarySidesHost_(sideIndex_,i);
135 sideNodeMapHost(i,0) = subcvCellTopo_.getNodeCount(sideDim, sideOrd);
136 for (ordinal_type j=0;j<sideNodeMapHost(i,0);++j)
137 sideNodeMapHost(i,j+1) = subcvCellTopo_.getNodeMap(sideDim, sideOrd, j);
138 }
139 sideNodeMap_ = Kokkos::create_mirror_view(typename DT::memory_space(), sideNodeMapHost);
140 Kokkos::deep_copy(sideNodeMap_, sideNodeMapHost);
141 }
142
143 template <typename DT, typename PT, typename WT>
144 void
145 CubatureControlVolumeBoundary<DT,PT,WT>::
146 getCubature( PointViewType cubPoints,
147 weightViewType cubWeights,
148 PointViewType cellCoords ) const {
149#ifdef HAVE_INTREPID2_DEBUG
150 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.rank() != 3, std::invalid_argument,
151 ">>> ERROR (CubatureControlVolumeBoundary): cubPoints must have rank 3 (C,P,D).");
152 INTREPID2_TEST_FOR_EXCEPTION( cubWeights.rank() != 2, std::invalid_argument,
153 ">>> ERROR (CubatureControlVolumeBoundary): cubWeights must have rank 2 (C,P).");
154 INTREPID2_TEST_FOR_EXCEPTION( cellCoords.rank() != 3, std::invalid_argument,
155 ">>> ERROR (CubatureControlVolumeBoundary): cellCoords must have rank 3 of (C,P,D).");
156
157 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(0) != cellCoords.extent(0) ||
158 cubPoints.extent(0) != cubWeights.extent(0), std::invalid_argument,
159 ">>> ERROR (CubatureControlVolume): cubPoints, cubWeights and cellCoords dimension(0) are not consistent, numCells");
160
161 {
162 const ordinal_type spaceDim = cellCoords.extent(2);
163 const ordinal_type sideDim = spaceDim - 1;
164 const size_type numPrimarySideNodes = primaryCellTopo_.getNodeCount(sideDim, sideIndex_);
165
166 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(1) != numPrimarySideNodes ||
167 cubWeights.extent(1) != numPrimarySideNodes, std::invalid_argument,
168 ">>> ERROR (CubatureControlVolume): cubPoints and cubWeights dimension(1) are not consistent, numPrimarySideNodes");
169 }
170 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(2) != cellCoords.extent(2) ||
171 static_cast<ordinal_type>(cubPoints.extent(2)) != getDimension(), std::invalid_argument,
172 ">>> ERROR (CubatureControlVolume): cubPoints, cellCoords, this->getDimension() are not consistent, spaceDim.");
173#endif
174
175 typedef Kokkos::DynRankView<PT,DT> tempPointViewType;
176 typedef Kokkos::DynRankView<PT,Kokkos::LayoutStride,DT> tempPointStrideViewType;
177
178 // get array dimensions
179 const ordinal_type numCells = cellCoords.extent(0);
180 const ordinal_type numNodesPerCell = cellCoords.extent(1);
181 const ordinal_type spaceDim = cellCoords.extent(2);
182 const ordinal_type sideDim = spaceDim - 1;
183
184 const ordinal_type numNodesPerSubcv = subcvCellTopo_.getNodeCount();
185 tempPointViewType subcvCoords("CubatureControlVolumeBoundary::subcvCoords",
186 numCells, numNodesPerCell, numNodesPerSubcv, spaceDim);
187 CellTools<DT>::getSubcvCoords(subcvCoords,
188 cellCoords,
189 primaryCellTopo_);
190
191 //const auto numPrimarySides = primaryCellTopo_.getSubcellCount(sideDim);
192 const ordinal_type numSubcvPoints = 1;
193
194 tempPointViewType subcvJacobian("CubatureControlVolumeBoundary::subcvJacobian",
195 numCells, numSubcvPoints, spaceDim, spaceDim);
196
197 tempPointViewType subcvJacobianDet("CubatureControlVolumeBoundary::subcvJacobianDet",
198 numCells, numSubcvPoints);
199
200 tempPointViewType weights("CubatureControlVolumeBoundary::subcvWeights",
201 numCells, 1);
202 Kokkos::deep_copy(weights, spaceDim == 2 ? 2.0 : 4.0);
203
204 tempPointViewType scratch("CubatureControlVolumeBoundary::scratch",
205 numCells*numSubcvPoints*spaceDim*spaceDim);
206
207 const ordinal_type numPrimarySideNodes = primaryCellTopo_.getNodeCount(sideDim, sideIndex_);
208 for (ordinal_type node=0;node<numPrimarySideNodes;++node) {
209 const auto sideRange = Kokkos::pair<ordinal_type,ordinal_type>(node, node+1);
210 const auto sidePoint = Kokkos::subdynrankview(sidePoints_, sideRange, Kokkos::ALL());
211
212 const auto idx = primaryCellTopo_.getNodeMap(sideDim, sideIndex_, node);
213 auto subcvCoordsNode = Kokkos::subdynrankview(subcvCoords, Kokkos::ALL(), idx, Kokkos::ALL(), Kokkos::ALL());
214
215 CellTools<DT>::setJacobian(subcvJacobian, // C, P, D, D
216 sidePoint, // P, D
217 subcvCoordsNode, // C, N, D
218 subcvCellTopo_);
219
220 CellTools<DT>::setJacobianDet(subcvJacobianDet, // C, P
221 subcvJacobian); // C, P, D, D
222
223 auto cubPointsNode = Kokkos::subdynrankview(cubPoints, Kokkos::ALL(), node, Kokkos::ALL());
224
225 typedef Kokkos::View<ordinal_type*,DT> mapViewType;
226 const auto sideNodeMap = Kokkos::subview(sideNodeMap_, node, Kokkos::ALL());
227
228 //typedef typename ExecSpace<typename PointViewType::execution_space,DT>::ExecSpaceType ExecSpaceType;
229
230 const auto loopSize = numCells;
231 Kokkos::RangePolicy<typename DT::execution_space,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
232
233 // compute centers
234 typedef Functor<PointViewType,tempPointStrideViewType,mapViewType> FunctorType;
235 Kokkos::parallel_for( policy, FunctorType(cubPointsNode,
236 subcvCoordsNode,
237 sideNodeMap) );
238
239 // compute weights
240 const auto sideOrd = boundarySidesHost_(sideIndex_, node);
241
242 // cub weights node requires to have rank 2
243 auto cubWeightsNode = Kokkos::subdynrankview(cubWeights, Kokkos::ALL(), sideRange);
244 switch (spaceDim) {
245 case 2: {
246 FunctionSpaceTools<DT>::computeEdgeMeasure(cubWeightsNode, // rank 2
247 subcvJacobian, // rank 4
248 weights, // rank 2
249 sideOrd,
250 subcvCellTopo_,
251 scratch);
252 break;
253 }
254 case 3: {
255 FunctionSpaceTools<DT>::computeFaceMeasure(cubWeightsNode,
256 subcvJacobian,
257 weights,
258 sideOrd,
259 subcvCellTopo_,
260 scratch);
261 break;
262 }
263 }
264 }
265 }
266}
267
268#endif
269
Intrepid2::CellTools::setJacobianDet
static void setJacobianDet(JacobianDetViewType jacobianDet, const JacobianViewType jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
Definition Intrepid2_CellToolsDefJacobian.hpp:874
Intrepid2::CellTools::getSubcvCoords
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.
Definition Intrepid2_CellToolsDefControlVolume.hpp:340
Intrepid2::CellTools::mapToReferenceSubcell
static void mapToReferenceSubcell(refSubcellViewType refSubcellPoints, const paramPointViewType paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
Intrepid2::CellTools::setJacobian
static void setJacobian(JacobianViewType jacobian, const PointViewType 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.
Definition Intrepid2_CellToolsDefJacobian.hpp:799
Intrepid2::CubatureControlVolumeBoundary::CubatureControlVolumeBoundary
CubatureControlVolumeBoundary(const shards::CellTopology cellTopology, const ordinal_type sideIndex)
Definition Intrepid2_CubatureControlVolumeBoundaryDef.hpp:24
Intrepid2::CubatureControlVolumeBoundary::getCubature
virtual void getCubature(PointViewType cubPoints, weightViewType cubWeights, PointViewType cellCoords) const override
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Definition Intrepid2_CubatureControlVolumeBoundaryDef.hpp:146
Intrepid2::FunctionSpaceTools::computeFaceMeasure
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...
Intrepid2::FunctionSpaceTools::computeEdgeMeasure
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...
Definition Intrepid2_FunctionSpaceToolsDef.hpp:630
Generated by doxygen 1.9.8