49 #ifndef __INTREPID2_CUBATURE_CONTROLVOLUME_BOUNDARY_DEF_HPP__
50 #define __INTREPID2_CUBATURE_CONTROLVOLUME_BOUNDARY_DEF_HPP__
55 template <
typename SpT,
typename PT,
typename WT>
58 const ordinal_type sideIndex) {
61 primaryCellTopo_ = cellTopology;
64 const ordinal_type spaceDim = primaryCellTopo_.getDimension();
67 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
70 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
77 sideIndex_ = sideIndex;
80 const ordinal_type sideDim = spaceDim - 1;
82 const ordinal_type numPrimarySides = primaryCellTopo_.getSubcellCount(sideDim);
83 const ordinal_type numPrimarySideNodes = primaryCellTopo_.getNodeCount(sideDim, sideIndex_);
85 const ordinal_type numSubcvSideNodes = subcvCellTopo_.getNodeCount(sideDim, 0);
87 boundarySidesHost_ = Kokkos::View<ordinal_type**,Kokkos::HostSpace>(
"CubatureControlVolumeBoundary::boundarySidesHost",
88 numPrimarySides, numSubcvSideNodes);
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;
100 case shards::Hexahedron<8>::key: {
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;
110 boundarySidesHost_(4,0) = 4;
111 boundarySidesHost_(4,1) = 4;
112 boundarySidesHost_(4,2) = 4;
113 boundarySidesHost_(4,3) = 4;
116 boundarySidesHost_(5,0) = 5;
117 boundarySidesHost_(5,1) = 5;
118 boundarySidesHost_(5,2) = 5;
119 boundarySidesHost_(5,3) = 5;
122 case shards::Tetrahedron<4>::key: {
123 boundarySidesHost_(0,0) = 0;
124 boundarySidesHost_(0,1) = 3;
125 boundarySidesHost_(0,2) = 0;
127 boundarySidesHost_(1,0) = 0;
128 boundarySidesHost_(1,1) = 3;
129 boundarySidesHost_(1,2) = 3;
131 boundarySidesHost_(2,0) = 3;
132 boundarySidesHost_(2,1) = 4;
133 boundarySidesHost_(2,2) = 0;
135 boundarySidesHost_(3,0) = 4;
136 boundarySidesHost_(3,1) = 4;
137 boundarySidesHost_(3,2) = 4;
141 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
142 ">>> ERROR (CubatureControlVolumeBoundary: invalid cell topology.");
146 Kokkos::DynRankView<PT,SpT> sideCenterLocal(
"CubatureControlVolumeBoundary::sideCenterLocal",
149 sidePoints_ = Kokkos::DynRankView<PT,SpT>(
"CubatureControlVolumeBoundary::sidePoints",
150 numPrimarySideNodes, spaceDim);
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());
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);
172 sideNodeMap_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), sideNodeMapHost);
173 Kokkos::deep_copy(sideNodeMap_, sideNodeMapHost);
176 template <
typename SpT,
typename PT,
typename WT>
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).");
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");
195 const ordinal_type spaceDim = cellCoords.extent(2);
196 const ordinal_type sideDim = spaceDim - 1;
197 const size_type numPrimarySideNodes = primaryCellTopo_.getNodeCount(sideDim, sideIndex_);
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");
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.");
208 typedef Kokkos::DynRankView<PT,SpT> tempPointViewType;
209 typedef Kokkos::DynRankView<PT,Kokkos::LayoutStride,SpT> tempPointStrideViewType;
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;
217 const ordinal_type numNodesPerSubcv = subcvCellTopo_.getNodeCount();
218 tempPointViewType subcvCoords(
"CubatureControlVolumeBoundary::subcvCoords",
219 numCells, numNodesPerCell, numNodesPerSubcv, spaceDim);
225 const ordinal_type numSubcvPoints = 1;
227 tempPointViewType subcvJacobian(
"CubatureControlVolumeBoundary::subcvJacobian",
228 numCells, numSubcvPoints, spaceDim, spaceDim);
230 tempPointViewType subcvJacobianDet(
"CubatureControlVolumeBoundary::subcvJacobianDet",
231 numCells, numSubcvPoints);
233 tempPointViewType weights(
"CubatureControlVolumeBoundary::subcvWeights",
235 Kokkos::deep_copy(weights, spaceDim == 2 ? 2.0 : 4.0);
237 tempPointViewType scratch(
"CubatureControlVolumeBoundary::scratch",
238 numCells*numSubcvPoints*spaceDim*spaceDim);
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());
245 const auto idx = primaryCellTopo_.getNodeMap(sideDim, sideIndex_, node);
246 auto subcvCoordsNode = Kokkos::subdynrankview(subcvCoords, Kokkos::ALL(), idx, Kokkos::ALL(), Kokkos::ALL());
256 auto cubPointsNode = Kokkos::subdynrankview(cubPoints, Kokkos::ALL(), node, Kokkos::ALL());
258 typedef Kokkos::View<ordinal_type*,SpT> mapViewType;
259 const auto sideNodeMap = Kokkos::subview(sideNodeMap_, node, Kokkos::ALL());
261 typedef typename ExecSpace<typename pointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
263 const auto loopSize = numCells;
264 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
268 Kokkos::parallel_for( policy, FunctorType(cubPointsNode,
273 const auto sideOrd = boundarySidesHost_(sideIndex_, node);
276 auto cubWeightsNode = Kokkos::subdynrankview(cubWeights, Kokkos::ALL(), sideRange);
279 std::cout <<
"subcv jacobian rank = " << subcvJacobian.rank() << std::endl;
CubatureControlVolumeBoundary(const shards::CellTopology cellTopology, const ordinal_type sideIndex)
virtual void getCubature(pointViewType cubPoints, weightViewType cubWeights, pointViewType cellCoords) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).