Intrepid2
Intrepid2_CubatureControlVolumeSideDef.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_SIDE_DEF_HPP__
50#define __INTREPID2_CUBATURE_CONTROLVOLUME_SIDE_DEF_HPP__
51
52namespace Intrepid2 {
53
54 template <typename DT, typename PT, typename WT>
56 CubatureControlVolumeSide(const shards::CellTopology cellTopology) {
57
58 // define primary cell topology with given one
59 primaryCellTopo_ = cellTopology;
60
61 // subcell is defined either hex or quad according to dimension
62 const auto spaceDim = primaryCellTopo_.getDimension();
63 switch (spaceDim) {
64 case 2:
65 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
66 break;
67 case 3:
68 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
69 break;
70 }
71
72 // computation order is always one;
73 degree_ = 1;
74
75 // precompute reference side points that are repeatedly used in get cubature
76 const auto maxNumNodesPerSide = 10;
77 const auto numSideNodeMaps = (spaceDim == 2 ? 1 : 2);
78
79 const ordinal_type sideOrd[2] = { 1, 5 };
80 Kokkos::View<ordinal_type**,Kokkos::HostSpace> sideNodeMapHost("CubatureControlVolumeSide::sideNodeMapHost",
81 numSideNodeMaps, maxNumNodesPerSide);
82
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);
89 }
90 sideNodeMap_ = Kokkos::create_mirror_view(typename DT::memory_space(), sideNodeMapHost);
91 Kokkos::deep_copy(sideNodeMap_, sideNodeMapHost);
92
93 Kokkos::DynRankView<PT,DT> sideCenterLocal("CubatureControlVolumeSide::sideCenterLocal",
94 1, sideDim);
95
96 // map to reference subcell function relies on uvm; some utility functions in cell tools still need uvm
97 sidePoints_ = Kokkos::DynRankView<PT,DT>("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());
102 sideCenterLocal,
103 sideDim,
104 sideOrd[i],
105 subcvCellTopo_);
106 }
107 }
108
109 template <typename DT, typename PT, typename WT>
110 void
112 getCubature( PointViewType cubPoints,
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).");
122
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");
126
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.");
131#endif
132 typedef Kokkos::DynRankView<PT,DT> tempPointViewType;
133
134 // get array dimensions
135 const auto numCells = cellCoords.extent(0);
136 const auto numNodesPerCell = cellCoords.extent(1);
137 const auto spaceDim = cellCoords.extent(2);
138
139 const auto numNodesPerSubcv = subcvCellTopo_.getNodeCount();
140 tempPointViewType subcvCoords("CubatureControlVolumeSide::subcvCoords",
141 numCells, numNodesPerCell, numNodesPerSubcv, spaceDim);
143 cellCoords,
144 primaryCellTopo_);
145
146 const auto numSideNodeMaps = (spaceDim == 2 ? 1 : 2);
147 const ordinal_type sideOrd[2] = { 1, 5 };
148
149 Kokkos::pair<ordinal_type,ordinal_type> nodeRangePerSide[2];
150
151 // the second rage is cell specific to handle remained sides
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;
157 break;
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;
163 break;
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;
169 break;
170 }
171
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);
177
178 tempPointViewType subcvSideNormal("CubatureControlVolume::subcvSideNormal",
179 numCells, numNodesPerThisSide, numSubcvPoints, spaceDim);
180
181 // numNodesPerCell is maximum 8; this repeated run is necessary because of cell tools input consideration
182 const auto sideRange = Kokkos::pair<ordinal_type,ordinal_type>(i, i+1);
183 const auto sidePoint = Kokkos::subdynrankview(sidePoints_, sideRange, Kokkos::ALL());
184
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());
189
190 CellTools<DT>::setJacobian(subcvJacobianNode, // C, P, D, D
191 sidePoint, // P, D
192 subcvCoordsNode, // C, N, D
193 subcvCellTopo_);
194
195 CellTools<DT>::getPhysicalSideNormals(subcvSideNormalNode, // C, P, D
196 subcvJacobianNode,
197 sideOrd[i],
198 subcvCellTopo_); // C, P, D, D
199 }
200
201 typedef Kokkos::View<ordinal_type*,DT> mapViewType;
202 const auto sideNodeMap = Kokkos::subview(sideNodeMap_, i, Kokkos::ALL());
203
204 //typedef typename ExecSpace<typename PointViewType::execution_space,DT>::ExecSpaceType ExecSpaceType;
205
206 const auto loopSize = numCells;
207 Kokkos::RangePolicy<typename DT::execution_space,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
208
210
211 auto cubPointsThisSide = Kokkos::subdynrankview(cubPoints, Kokkos::ALL(), nodeRangePerSide[i], Kokkos::ALL());
212 auto cubWeightsThisSide = Kokkos::subdynrankview(cubWeights, Kokkos::ALL(), nodeRangePerSide[i], Kokkos::ALL());
213
214 Kokkos::parallel_for( policy, FunctorType(cubPointsThisSide,
215 cubWeightsThisSide,
216 subcvCoords,
217 subcvSideNormal,
218 sideNodeMap) );
219 }
220 }
221
222}
223
224#endif
225
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.
static void getPhysicalSideNormals(Kokkos::DynRankView< sideNormalValueType, sideNormalProperties... > sideNormals, const Kokkos::DynRankView< worksetJacobianValueType, worksetJacobianProperties... > worksetJacobians, const ordinal_type worksetSideOrd, const shards::CellTopology parentCell)
Computes non-normalized normal vectors to physical sides in a side workset .
CubatureControlVolumeSide(const shards::CellTopology cellTopology)
virtual void getCubature(PointViewType cubPoints, weightViewType cubWeights, PointViewType cellCoords) const override
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).