Intrepid2
Intrepid2_DefaultCubatureFactoryDef.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_DEFAULT_CUBATURE_FACTORY_DEF_HPP__
50#define __INTREPID2_DEFAULT_CUBATURE_FACTORY_DEF_HPP__
51
52namespace Intrepid2 {
53
54 // first create method
55 template<typename DT, typename PT, typename WT>
56 Teuchos::RCP<Cubature<DT,PT,WT> >
58 create( unsigned topologyKey,
59 const std::vector<ordinal_type> &degree,
60 const EPolyType polytype ) {
61
62 // Create generic cubature.
63 Teuchos::RCP<Cubature<DT,PT,WT> > r_val;
64
65 switch (topologyKey) {
66 case shards::Line<>::key: {
67 INTREPID2_TEST_FOR_EXCEPTION( degree.size() < 1, std::invalid_argument,
68 ">>> ERROR (DefaultCubatureFactory): Provided degree array is of insufficient length.");
69 if (isValidPolyType(polytype))
70 r_val = Teuchos::rcp(new CubaturePolylib<DT,PT,WT>(degree[0], polytype));
71 else
72 r_val = Teuchos::rcp(new CubatureDirectLineGauss<DT,PT,WT>(degree[0]));
73 break;
74 }
75 case shards::Triangle<>::key: {
76 INTREPID2_TEST_FOR_EXCEPTION( degree.size() < 1, std::invalid_argument,
77 ">>> ERROR (DefaultCubatureFactory): Provided degree array is of insufficient length.");
78 r_val = Teuchos::rcp(new CubatureDirectTriDefault<DT,PT,WT>(degree[0]));
79 break;
80 }
81 case shards::Quadrilateral<>::key:
82 case shards::ShellQuadrilateral<>::key: {
83 INTREPID2_TEST_FOR_EXCEPTION( degree.size() < 2, std::invalid_argument,
84 ">>> ERROR (DefaultCubatureFactory): Provided degree array is of insufficient length.");
85 if (isValidPolyType(polytype)) {
86 const auto x_line = CubaturePolylib<DT,PT,WT>(degree[0], polytype);
87 const auto y_line = ( degree[1] == degree[0] ? x_line : CubaturePolylib<DT,PT,WT>(degree[1], polytype) );
88 r_val = Teuchos::rcp(new CubatureTensor<DT,PT,WT>( x_line, y_line ));
89 } else {
90 const auto x_line = CubatureDirectLineGauss<DT,PT,WT>(degree[0]);
91 const auto y_line = ( degree[1] == degree[0] ? x_line : CubatureDirectLineGauss<DT,PT,WT>(degree[1]) );
92 r_val = Teuchos::rcp(new CubatureTensor<DT,PT,WT>( x_line, y_line ));
93 }
94 break;
95 }
96 case shards::Tetrahedron<>::key: {
97 INTREPID2_TEST_FOR_EXCEPTION( degree.size() < 1, std::invalid_argument,
98 ">>> ERROR (DefaultCubatureFactory): Provided degree array is of insufficient length.");
99 r_val = Teuchos::rcp(new CubatureDirectTetDefault<DT,PT,WT>(degree[0]));
100 break;
101 }
102 case shards::Hexahedron<>::key: {
103 INTREPID2_TEST_FOR_EXCEPTION( degree.size() < 3, std::invalid_argument,
104 ">>> ERROR (DefaultCubatureFactory): Provided degree array is of insufficient length.");
105 if (isValidPolyType(polytype)) {
106 const auto x_line = CubaturePolylib<DT,PT,WT>(degree[0], polytype);
107 const auto y_line = ( degree[1] == degree[0] ? x_line : CubaturePolylib<DT,PT,WT>(degree[1], polytype) );
108 const auto z_line = ( degree[2] == degree[0] ? x_line :
109 degree[2] == degree[1] ? y_line : CubaturePolylib<DT,PT,WT>(degree[2], polytype) );
110
111 r_val = Teuchos::rcp(new CubatureTensor<DT,PT,WT>( x_line, y_line, z_line ));
112 } else {
113 const auto x_line = CubatureDirectLineGauss<DT,PT,WT>(degree[0]);
114 const auto y_line = ( degree[1] == degree[0] ? x_line : CubatureDirectLineGauss<DT,PT,WT>(degree[1]) );
115 const auto z_line = ( degree[2] == degree[0] ? x_line :
116 degree[2] == degree[1] ? y_line : CubatureDirectLineGauss<DT,PT,WT>(degree[2]) );
117
118 r_val = Teuchos::rcp(new CubatureTensor<DT,PT,WT>( x_line, y_line, z_line ));
119 }
120 break;
121 }
122 case shards::Wedge<>::key: {
123 INTREPID2_TEST_FOR_EXCEPTION( degree.size() < 2, std::invalid_argument,
124 ">>> ERROR (DefaultCubatureFactory): Provided degree array is of insufficient length.")
125 {
126 const auto xy_tri = CubatureDirectTriDefault<DT,PT,WT>(degree[0]);
127 const auto z_line = CubatureDirectLineGauss<DT,PT,WT>(degree[1]);
128 r_val = Teuchos::rcp(new CubatureTensor<DT,PT,WT>( xy_tri, z_line ));
129 }
130 break;
131 }
132 case shards::Pyramid<>::key: {
133 INTREPID2_TEST_FOR_EXCEPTION( degree.size() < 1, std::invalid_argument,
134 ">>> ERROR (DefaultCubatureFactory): Provided degree array is of insufficient length.");
135 {
136 // if direct gauss is used for pyramid
137 // need over-integration to account for the additional transformation
138 const auto xy_line = CubatureDirectLineGauss<DT,PT,WT>(degree[0]);
139 const auto z_line = CubatureDirectLineGaussJacobi20<DT,PT,WT>(degree[0]);
140 r_val = Teuchos::rcp(new CubatureTensorPyr<DT,PT,WT>( xy_line, xy_line, z_line ));
141 }
142 break;
143 }
144 default: {
145 INTREPID2_TEST_FOR_EXCEPTION( ( (topologyKey != shards::Line<>::key) &&
146 (topologyKey != shards::Triangle<>::key) &&
147 (topologyKey != shards::Quadrilateral<>::key) &&
148 (topologyKey != shards::ShellQuadrilateral<>::key) &&
149 (topologyKey != shards::Tetrahedron<>::key) &&
150 (topologyKey != shards::Hexahedron<>::key) &&
151 (topologyKey != shards::Pyramid<>::key) &&
152 (topologyKey != shards::Wedge<>::key) ),
153 std::invalid_argument,
154 ">>> ERROR (DefaultCubatureFactory): Invalid cell topology prevents cubature creation.");
155 }
156 }
157 return r_val;
158 }
159
160 template<typename DT, typename PT, typename WT>
161 Teuchos::RCP<Cubature<DT,PT,WT> >
163 create( const shards::CellTopology cellTopology,
164 const std::vector<ordinal_type> &degree,
165 const EPolyType polytype) {
166 return create<DT,PT,WT>(cellTopology.getBaseKey(), degree, polytype);
167 }
168
169
170 template<typename DT, typename PT, typename WT>
171 Teuchos::RCP<Cubature<DT,PT,WT> >
173 create( unsigned topologyKey,
174 const ordinal_type degree,
175 const EPolyType polytype ) {
176 // uniform order for 3 axes
177 const std::vector<ordinal_type> degreeArray(3, degree);
178 return create<DT,PT,WT>(topologyKey, degreeArray, polytype);
179 }
180
181 template<typename DT, typename PT, typename WT>
182 Teuchos::RCP<Cubature<DT,PT,WT> >
184 create( const shards::CellTopology cellTopology,
185 const ordinal_type degree,
186 const EPolyType polytype ) {
187 // uniform order for 3 axes
188 const std::vector<ordinal_type> degreeArray(3, degree);
189 return create<DT,PT,WT>(cellTopology.getBaseKey(), degreeArray, polytype);
190 }
191
192
193 // template<typename DT>
194 // template<typename cellVertexValueType, class ...cellVertexProperties>
195 // Teuchos::RCP<Cubature<DT,PT,WT> >
196 // DefaultCubatureFactory::
197 // create<DT,PT,WT>(const shards::CellTopology& cellTopology,
198 // const Kokkos::DynRankView<cellVertexValueType,cellVertexProperties> cellVertices,
199 // ordinal_type degree){
200 // return Teuchos::rcp(new CubaturePolygon<DT,PT,WT>(cellTopology,cellVertices, degree));
201 // }
202
203} // namespace Intrepid2
204
205#endif
KOKKOS_FORCEINLINE_FUNCTION bool isValidPolyType(const EPolyType polytype)
Verifies validity of a PolyType enum.
static Teuchos::RCP< Cubature< DeviceType, pointValueType, weightValueType > > create(unsigned topologyKey, const std::vector< ordinal_type > &degree, const EPolyType polytype=POLYTYPE_MAX)
Factory method.