Intrepid2
Intrepid2_Cubature.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_HPP__
50#define __INTREPID2_CUBATURE_HPP__
51
52#include "Intrepid2_ConfigDefs.hpp"
55#include "Intrepid2_Types.hpp"
56#include "Intrepid2_Utils.hpp"
57
58namespace Intrepid2 {
59
60 /* \struct Intrepid2::CubatureTemplate
61 \brief Template for the cubature rules used by Intrepid. Cubature template consists of
62 cubature points and cubature weights. Intrepid provides a collection of cubature
63 templates for most standard cell topologies. The templates are defined in reference
64 coordinates using a standard reference cell for each canonical cell type. Cubature
65 points are always specified by a triple of (X,Y,Z) coordinates even if the cell
66 dimension is less than 3. The unused dimensions should be padded by zeroes.
67
68 For example, a set of Gauss rules on [-1,1] looks as the following array of CubatureTemplate structs:
69
70 \verbatim
71 cubature_rule[4] =
72 { // Collection of Gauss rules on [-1,1]
73 {
74 1, ----> number of points in the rule
75 {{0.0,0.0,0.0}}, ----> X,Y,Z coordinates of the cubature points
76 {0.5} ----> the cubature weight
77 },
78 {
79 2,
80 {{-sqrt(1.0/3.0),0.0,0.0},
81 {+sqrt(1.0/3.0),0.0,0.0}},
82 {1.0,1.0}
83 },
84 {
85 3,
86 {{-sqrt(3.0/5.0),0.0,0.0},
87 {0.0,0.0,0.0},
88 {+sqrt(3.0/5.0),0.0,0.0}},
89 {5.0/9.0, 8.0/9.0,5.0/9.0}
90 },
91 {
92 4,
93 {{-sqrt((3.0+4.0*sqrt(0.3))/7.0),0.0,0.0},
94 {-sqrt((3.0-4.0*sqrt(0.3))/7.0),0.0,0.0},
95 {+sqrt((3.0-4.0*sqrt(0.3))/7.0),0.0,0.0},
96 {+sqrt((3.0+4.0*sqrt(0.3))/7.0),0.0,0.0}},
97 //
98 {0.5-sqrt(10.0/3.0)/12.0,
99 0.5+sqrt(10.0/3.0)/12.0,
100 0.5+sqrt(10.0/3.0)/12.0,
101 0.5-sqrt(10.0/3.0)/12.0}
102 }
103 }; // end Gauss
104 \endverbatim
105
106 Also see data member documentation.
107 */
108
109
121 template<typename DeviceType = void,
122 typename pointValueType = double,
123 typename weightValueType = double>
124 class Cubature {
125 public:
126 using ExecSpaceType = typename DeviceType::execution_space;
127 using PointViewType = Kokkos::DynRankView<pointValueType,Kokkos::LayoutStride,DeviceType>;
128 using weightViewType = Kokkos::DynRankView<weightValueType,Kokkos::LayoutStride,DeviceType>;
129
130 using PointViewTypeAllocatable = Kokkos::DynRankView<pointValueType,DeviceType>; // uses default layout; allows us to allocate (in contrast to LayoutStride)
131 using WeightViewTypeAllocatable = Kokkos::DynRankView<weightValueType,DeviceType>; // uses default layout; allows us to allocate (in contrast to LayoutStride)
134
140 {
141 // default implementation has trivial tensor structure
142 PointViewTypeAllocatable allPointData("cubature points",this->getNumPoints(),this->getDimension());
143 Kokkos::Array< PointViewTypeAllocatable, 1> cubaturePointComponents;
144 cubaturePointComponents[0] = allPointData;
145 return TensorPointDataType(cubaturePointComponents);
146 }
147
153 {
154 // default implementation has trivial tensor structure
155 using WeightDataType = Data<weightValueType,DeviceType>;
156 WeightViewTypeAllocatable allWeightData("cubature weights",this->getNumPoints());
157 Kokkos::Array< WeightDataType, 1> cubatureWeightComponents;
158 cubatureWeightComponents[0] = WeightDataType(allWeightData);
159 return TensorWeightDataType(cubatureWeightComponents);
160 }
161
168 virtual
169 void
170 getCubature( PointViewType /* cubPoints */,
171 weightViewType /* cubWeights */ ) const {
172 INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
173 ">>> ERROR (Cubature::getCubature): this method should be overridden by derived classes.");
174 }
175
183 virtual
184 void
185 getCubature( PointViewType /* cubPoints */,
186 weightViewType /* cubWeights */,
187 PointViewType /* cellVertices */) const {
188 INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
189 ">>> ERROR (Cubature::getCubature): this method should be overridden by derived classes.");
190 }
191
197 virtual
198 void
199 getCubature( const TensorPointDataType & tensorCubPoints,
200 const TensorWeightDataType & tensorCubWeights) const {
201 // default implementation has trivial tensor product structure.
202
203 INTREPID2_TEST_FOR_EXCEPTION(1 != tensorCubPoints.numTensorComponents(), std::logic_error, "default implementation of getCubature only supports trivial tensor structure -- numTensorComponents() must be 1");
204 INTREPID2_TEST_FOR_EXCEPTION(1 != tensorCubWeights.numTensorComponents(), std::logic_error, "default implementation of getCubature only supports trivial tensor structure -- numTensorComponents() must be 1");
205
206 auto underlyingPointView = tensorCubPoints.getTensorComponent(0);
207
208 this->getCubature(underlyingPointView, tensorCubWeights.getTensorComponent(0).getUnderlyingView());
209 }
210
213 virtual
214 ordinal_type
215 getNumPoints() const {
216 INTREPID2_TEST_FOR_WARNING( true,
217 ">>> ERROR (Cubature::getNumPoints): this method should be overridden by derived classes.");
218 return 0;
219 }
220
221
224 virtual
225 ordinal_type
226 getDimension() const {
227 INTREPID2_TEST_FOR_WARNING( true,
228 ">>> ERROR (Cubature::getDimension): this method should be overridden by derived classes.");
229 return 0;
230 }
231
234 virtual
235 ordinal_type
236 getAccuracy() const {
237 INTREPID2_TEST_FOR_WARNING( true,
238 ">>> ERROR (Cubature::getDimension): this method should be overridden by derived classes.");
239 return 0;
240 }
241
244 virtual
245 const char*
246 getName() const {
247 return "Cubature";
248 }
249
250 Cubature() = default;
251 virtual ~Cubature() {}
252
253 };
254
255}// end namespace Intrepid2
256
257
388#endif
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
Contains definitions of custom data types in Intrepid2.
Header function for Intrepid2::Util class and other utility functions.
Defines the base class for cubature (integration) rules in Intrepid.
virtual ordinal_type getDimension() const
Returns dimension of the integration domain.
virtual const char * getName() const
Returns cubature name.
virtual void getCubature(PointViewType, weightViewType) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
virtual void getCubature(PointViewType, weightViewType, PointViewType) const
Returns cubature points and weights on physical cells (return arrays must be pre-sized/pre-allocated)...
virtual ordinal_type getNumPoints() const
Returns the number of cubature points.
virtual TensorPointDataType allocateCubaturePoints() const
Returns a points container appropriate for passing to getCubature().
virtual TensorWeightDataType allocateCubatureWeights() const
Returns a weight container appropriate for passing to getCubature().
virtual ordinal_type getAccuracy() const
Returns dimension of the integration domain.
virtual void getCubature(const TensorPointDataType &tensorCubPoints, const TensorWeightDataType &tensorCubWeights) const
Returns tensor cubature points and weights. For non-tensor cubatures, the tensor structures are trivi...
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
KOKKOS_INLINE_FUNCTION enable_if_t< rank==1, const Kokkos::View< typename RankExpander< DataScalar, rank >::value_type, DeviceType > & > getUnderlyingView() const
Returns the underlying view. Throws an exception if the underlying view is not rank 1.
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
KOKKOS_INLINE_FUNCTION const Data< Scalar, DeviceType > & getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Return the number of tensorial components.
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Returns the number of tensorial components.
KOKKOS_INLINE_FUNCTION ScalarView< PointScalar, DeviceType > getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.