Intrepid2
Intrepid2_CubatureControlVolumeDef.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_DEF_HPP__
50#define __INTREPID2_CUBATURE_CONTROLVOLUME_DEF_HPP__
51
52namespace Intrepid2 {
53
54 template <typename DT, typename PT, typename WT>
56 CubatureControlVolume(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 switch (primaryCellTopo_.getDimension()) {
63 case 2:
64 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
65 break;
66 case 3:
67 subcvCellTopo_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
68 break;
69 }
70
71 // computation order is always one;
72 degree_ = 1;
73
74 // create subcell cubature points and weights and cache them
75 const ordinal_type subcvDegree = 2;
76 auto subcvCubature = DefaultCubatureFactory::create<DT,PT,WT>(subcvCellTopo_, subcvDegree);
77
78 const ordinal_type numSubcvPoints = subcvCubature->getNumPoints();
79 const ordinal_type subcvDim = subcvCubature->getDimension();
80
81 subcvCubaturePoints_ = Kokkos::DynRankView<PT,DT>("CubatureControlVolume::subcvCubaturePoints_",
82 numSubcvPoints, subcvDim);
83 subcvCubatureWeights_ = Kokkos::DynRankView<WT,DT>("CubatureControlVolume::subcvCubatureWeights_",
84 numSubcvPoints);
85
86 subcvCubature->getCubature(subcvCubaturePoints_,
87 subcvCubatureWeights_);
88 }
89
90 template <typename DT, typename PT, typename WT>
91 void
93 getCubature( PointViewType cubPoints,
94 weightViewType cubWeights,
95 PointViewType cellCoords ) const {
96#ifdef HAVE_INTREPID2_DEBUG
97 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.rank() != 3, std::invalid_argument,
98 ">>> ERROR (CubatureControlVolume): cubPoints must have rank 3 (C,P,D).");
99 INTREPID2_TEST_FOR_EXCEPTION( cubWeights.rank() != 2, std::invalid_argument,
100 ">>> ERROR (CubatureControlVolume): cubWeights must have rank 2 (C,P).");
101 INTREPID2_TEST_FOR_EXCEPTION( cellCoords.rank() != 3, std::invalid_argument,
102 ">>> ERROR (CubatureControlVolume): cellCoords must have rank 3 of (C,P,D).");
103
104 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(0) != cellCoords.extent(0) ||
105 cubPoints.extent(0) != cubWeights.extent(0), std::invalid_argument,
106 ">>> ERROR (CubatureControlVolume): cubPoints, cubWeights and cellCoords dimension(0) are not consistent, numCells");
107
108 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(1) != cellCoords.extent(1) ||
109 cubPoints.extent(1) != cubWeights.extent(1), std::invalid_argument,
110 ">>> ERROR (CubatureControlVolume): cubPoints, cubWeights and cellCoords dimension(1) are not consistent, numNodesPerCell");
111
112 INTREPID2_TEST_FOR_EXCEPTION( cubPoints.extent(2) != cellCoords.extent(2) ||
113 static_cast<ordinal_type>(cubPoints.extent(2)) != getDimension(), std::invalid_argument,
114 ">>> ERROR (CubatureControlVolume): cubPoints, cellCoords, this->getDimension() are not consistent, spaceDim.");
115#endif
116 typedef Kokkos::DynRankView<PT,DT> tempPointViewType;
117
118 // get array dimensions
119 const ordinal_type numCells = cellCoords.extent(0);
120 const ordinal_type numNodesPerCell = cellCoords.extent(1);
121 const ordinal_type spaceDim = cellCoords.extent(2);
122
123 const ordinal_type numNodesPerSubcv = subcvCellTopo_.getNodeCount();
124 tempPointViewType subcvCoords("CubatureControlVolume::subcvCoords_",
125 numCells, numNodesPerCell, numNodesPerSubcv, spaceDim);
127 cellCoords,
128 primaryCellTopo_);
129
130 const ordinal_type numSubcvPoints = subcvCubaturePoints_.extent(0);
131 tempPointViewType subcvJacobian("CubatureControlVolume::subcvJacobian_",
132 numCells, numNodesPerCell, numSubcvPoints, spaceDim, spaceDim);
133
134 tempPointViewType subcvJacobianDet("CubatureControlVolume::subcvJacobDet_",
135 numCells, numNodesPerCell, numSubcvPoints);
136
137 // numNodesPerCell is maximum 8; this repeated run is necessary because of cell tools input consideration
138 for (ordinal_type node=0;node<numNodesPerCell;++node) {
139 auto subcvJacobianNode = Kokkos::subdynrankview(subcvJacobian, Kokkos::ALL(), node, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
140 auto subcvCoordsNode = Kokkos::subdynrankview(subcvCoords, Kokkos::ALL(), node, Kokkos::ALL(), Kokkos::ALL());
141 auto subcvJacobianDetNode = Kokkos::subdynrankview(subcvJacobianDet, Kokkos::ALL(), node, Kokkos::ALL());
142
143 CellTools<DT>::setJacobian(subcvJacobianNode, // C, P, D, D
144 subcvCubaturePoints_, // P, D
145 subcvCoordsNode, // C, N, D
146 subcvCellTopo_);
147
148 CellTools<DT>::setJacobianDet(subcvJacobianDetNode, // C, P
149 subcvJacobianNode); // C, P, D, D
150 }
151
152 //typedef typename ExecSpace<typename PointViewType::execution_space,DT>::ExecSpaceType ExecSpaceType;
153
154 const auto loopSize = numCells;
155 Kokkos::RangePolicy<typename DT::execution_space,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
156
158 Kokkos::parallel_for( policy, FunctorType(cubPoints,
159 cubWeights,
160 subcvCoords,
161 subcvCubatureWeights_,
162 subcvJacobianDet) );
163 }
164
165}
166
167#endif
168
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.
virtual void getCubature(PointViewType cubPoints, weightViewType cubWeights, PointViewType cellCoords) const override
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
CubatureControlVolume(const shards::CellTopology cellTopology)