49 #ifndef __INTREPID2_CUBATURE_CONTROLVOLUME_SIDE_DEF_HPP__
50 #define __INTREPID2_CUBATURE_CONTROLVOLUME_SIDE_DEF_HPP__
54 template <
typename SpT,
typename PT,
typename WT>
59 primaryCellTopo_ = cellTopology;
62 const auto spaceDim = primaryCellTopo_.getDimension();
65 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
68 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
76 const auto maxNumNodesPerSide = 10;
77 const auto numSideNodeMaps = (spaceDim == 2 ? 1 : 2);
79 const ordinal_type sideOrd[2] = { 1, 5 };
80 Kokkos::View<ordinal_type**,Kokkos::HostSpace> sideNodeMapHost(
"CubatureControlVolumeSide::sideNodeMapHost",
81 numSideNodeMaps, maxNumNodesPerSide);
83 const auto sideDim = spaceDim - 1;
84 for (ordinal_type i=0;i<numSideNodeMaps;++i) {
85 const auto side = sideOrd[i];
86 sideNodeMapHost(i,0) = subcvCellTopo_.getNodeCount(sideDim, side);
87 for (ordinal_type j=0;j<sideNodeMapHost(i,0);++j)
88 sideNodeMapHost(i,j+1) = subcvCellTopo_.getNodeMap(sideDim, side, j);
90 sideNodeMap_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), sideNodeMapHost);
91 Kokkos::deep_copy(sideNodeMap_, sideNodeMapHost);
93 Kokkos::DynRankView<PT,SpT> sideCenterLocal(
"CubatureControlVolumeSide::sideCenterLocal",
97 sidePoints_ = Kokkos::DynRankView<PT,SpT>(
"CubatureControlVolumeSide::sidePoints", numSideNodeMaps, spaceDim);
98 for (ordinal_type i=0;i<numSideNodeMaps;++i) {
99 const auto sideRange = Kokkos::pair<ordinal_type,ordinal_type>(i, i+1);
100 auto sidePoint = Kokkos::subdynrankview(sidePoints_, sideRange, Kokkos::ALL());
109 template <
typename SpT,
typename PT,
typename WT>
113 weightViewType cubWeights,
114 pointViewType cellCoords )
const {
115 #ifdef HAVE_INTREPID2_DEBUG
116 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.rank() != 3, std::invalid_argument,
117 ">>> ERROR (CubatureControlVolumeSide): cubPoints must have rank 3 (C,P,D).");
118 INTREPID2_TEST_FOR_EXCEPTION( cubWeights.rank() != 3, std::invalid_argument,
119 ">>> ERROR (CubatureControlVolumeSide): cubWeights must have rank 2 (C,P,D).");
120 INTREPID2_TEST_FOR_EXCEPTION( cellCoords.rank() != 3, std::invalid_argument,
121 ">>> ERROR (CubatureControlVolumeSide): cellCoords must have rank 3 of (C,P,D).");
123 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(0) != cellCoords.extent(0) ||
124 cubPoints.extent(0) != cubWeights.extent(0), std::invalid_argument,
125 ">>> ERROR (CubatureControlVolumeSide): cubPoints, cubWeights and cellCoords dimension(0) are not consistent, numCells");
127 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(2) != cellCoords.extent(2) ||
128 cubPoints.extent(2) != cubWeights.extent(2) ||
129 static_cast<ordinal_type
>(cubPoints.extent(2)) != getDimension(), std::invalid_argument,
130 ">>> ERROR (CubatureControlVolumeSide): cubPoints, cellCoords, this->getDimension() are not consistent, spaceDim.");
132 typedef Kokkos::DynRankView<PT,SpT> tempPointViewType;
135 const auto numCells = cellCoords.extent(0);
136 const auto numNodesPerCell = cellCoords.extent(1);
137 const auto spaceDim = cellCoords.extent(2);
139 const auto numNodesPerSubcv = subcvCellTopo_.getNodeCount();
140 tempPointViewType subcvCoords(
"CubatureControlVolumeSide::subcvCoords",
141 numCells, numNodesPerCell, numNodesPerSubcv, spaceDim);
146 const auto numSideNodeMaps = (spaceDim == 2 ? 1 : 2);
147 const ordinal_type sideOrd[2] = { 1, 5 };
149 Kokkos::pair<ordinal_type,ordinal_type> nodeRangePerSide[2];
152 switch (primaryCellTopo_.getKey()) {
153 case shards::Triangle<3>::key:
154 case shards::Quadrilateral<4>::key:
155 nodeRangePerSide[0].first = 0;
156 nodeRangePerSide[0].second = nodeRangePerSide[0].first + numNodesPerCell;
158 case shards::Hexahedron<8>::key:
159 nodeRangePerSide[0].first = 0;
160 nodeRangePerSide[0].second = nodeRangePerSide[0].first + numNodesPerCell;
161 nodeRangePerSide[1].first = numNodesPerCell;
162 nodeRangePerSide[1].second = nodeRangePerSide[1].first + 4;
164 case shards::Tetrahedron<4>::key:
165 nodeRangePerSide[0].first = 0;
166 nodeRangePerSide[0].second = nodeRangePerSide[0].first + numNodesPerCell;
167 nodeRangePerSide[1].first = 3;
168 nodeRangePerSide[1].second = nodeRangePerSide[1].first + 3;
172 for (ordinal_type i=0;i<numSideNodeMaps;++i) {
173 const auto numSubcvPoints = 1;
174 const auto numNodesPerThisSide = nodeRangePerSide[i].second - nodeRangePerSide[i].first;
175 tempPointViewType subcvJacobian(
"CubatureControlVolume::subcvJacobian",
176 numCells, numNodesPerThisSide, numSubcvPoints, spaceDim, spaceDim);
178 tempPointViewType subcvSideNormal(
"CubatureControlVolume::subcvSideNormal",
179 numCells, numNodesPerThisSide, numSubcvPoints, spaceDim);
182 const auto sideRange = Kokkos::pair<ordinal_type,ordinal_type>(i, i+1);
183 const auto sidePoint = Kokkos::subdynrankview(sidePoints_, sideRange, Kokkos::ALL());
185 for (
auto node=0;node<numNodesPerThisSide;++node) {
186 auto subcvJacobianNode = Kokkos::subdynrankview(subcvJacobian, Kokkos::ALL(), node, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
187 auto subcvCoordsNode = Kokkos::subdynrankview(subcvCoords, Kokkos::ALL(), node, Kokkos::ALL(), Kokkos::ALL());
188 auto subcvSideNormalNode = Kokkos::subdynrankview(subcvSideNormal, Kokkos::ALL(), node, Kokkos::ALL(), Kokkos::ALL());
201 typedef Kokkos::View<ordinal_type*,SpT> mapViewType;
202 const auto sideNodeMap = Kokkos::subview(sideNodeMap_, i, Kokkos::ALL());
204 typedef typename ExecSpace<typename pointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
206 const auto loopSize = numCells;
207 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
211 auto cubPointsThisSide = Kokkos::subdynrankview(cubPoints, Kokkos::ALL(), nodeRangePerSide[i], Kokkos::ALL());
212 auto cubWeightsThisSide = Kokkos::subdynrankview(cubWeights, Kokkos::ALL(), nodeRangePerSide[i], Kokkos::ALL());
214 Kokkos::parallel_for( policy, FunctorType(cubPointsThisSide,
virtual void getCubature(pointViewType cubPoints, weightViewType cubWeights, pointViewType cellCoords) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
CubatureControlVolumeSide(const shards::CellTopology cellTopology)