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 
52 namespace Intrepid2 {
53 
54  // first create method
55  template<typename SpT, typename PT, typename WT>
56  Teuchos::RCP<Cubature<SpT,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<SpT,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<SpT,PT,WT>(degree[0], polytype));
71  else
72  r_val = Teuchos::rcp(new CubatureDirectLineGauss<SpT,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<SpT,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<SpT,PT,WT>(degree[0], polytype);
87  const auto y_line = ( degree[1] == degree[0] ? x_line : CubaturePolylib<SpT,PT,WT>(degree[1], polytype) );
88  r_val = Teuchos::rcp(new CubatureTensor<SpT,PT,WT>( x_line, y_line ));
89  } else {
90  const auto x_line = CubatureDirectLineGauss<SpT,PT,WT>(degree[0]);
91  const auto y_line = ( degree[1] == degree[0] ? x_line : CubatureDirectLineGauss<SpT,PT,WT>(degree[1]) );
92  r_val = Teuchos::rcp(new CubatureTensor<SpT,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<SpT,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<SpT,PT,WT>(degree[0], polytype);
107  const auto y_line = ( degree[1] == degree[0] ? x_line : CubaturePolylib<SpT,PT,WT>(degree[1], polytype) );
108  const auto z_line = ( degree[2] == degree[0] ? x_line :
109  degree[2] == degree[1] ? y_line : CubaturePolylib<SpT,PT,WT>(degree[2], polytype) );
110 
111  r_val = Teuchos::rcp(new CubatureTensor<SpT,PT,WT>( x_line, y_line, z_line ));
112  } else {
113  const auto x_line = CubatureDirectLineGauss<SpT,PT,WT>(degree[0]);
114  const auto y_line = ( degree[1] == degree[0] ? x_line : CubatureDirectLineGauss<SpT,PT,WT>(degree[1]) );
115  const auto z_line = ( degree[2] == degree[0] ? x_line :
116  degree[2] == degree[1] ? y_line : CubatureDirectLineGauss<SpT,PT,WT>(degree[2]) );
117 
118  r_val = Teuchos::rcp(new CubatureTensor<SpT,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<SpT,PT,WT>(degree[0]);
127  const auto z_line = CubatureDirectLineGauss<SpT,PT,WT>(degree[1]);
128  r_val = Teuchos::rcp(new CubatureTensor<SpT,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<SpT,PT,WT>(degree[0]);
139  const auto z_line = CubatureDirectLineGaussJacobi20<SpT,PT,WT>(degree[0]);
140  r_val = Teuchos::rcp(new CubatureTensorPyr<SpT,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 SpT, typename PT, typename WT>
161  Teuchos::RCP<Cubature<SpT,PT,WT> >
163  create( const shards::CellTopology cellTopology,
164  const std::vector<ordinal_type> &degree,
165  const EPolyType polytype) {
166  return create<SpT,PT,WT>(cellTopology.getBaseKey(), degree, polytype);
167  }
168 
169 
170  template<typename SpT, typename PT, typename WT>
171  Teuchos::RCP<Cubature<SpT,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<SpT,PT,WT>(topologyKey, degreeArray, polytype);
179  }
180 
181  template<typename SpT, typename PT, typename WT>
182  Teuchos::RCP<Cubature<SpT,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<SpT,PT,WT>(cellTopology.getBaseKey(), degreeArray, polytype);
190  }
191 
192 
193  // template<typename SpT>
194  // template<typename cellVertexValueType, class ...cellVertexProperties>
195  // Teuchos::RCP<Cubature<SpT,PT,WT> >
196  // DefaultCubatureFactory::
197  // create<SpT,PT,WT>(const shards::CellTopology& cellTopology,
198  // const Kokkos::DynRankView<cellVertexValueType,cellVertexProperties> cellVertices,
199  // ordinal_type degree){
200  // return Teuchos::rcp(new CubaturePolygon<SpT,PT,WT>(cellTopology,cellVertices, degree));
201  // }
202 
203 } // namespace Intrepid2
204 
205 #endif
static Teuchos::RCP< Cubature< ExecSpaceType, pointValueType, weightValueType > > create(unsigned topologyKey, const std::vector< ordinal_type > &degree, const EPolyType polytype=POLYTYPE_MAX)
Factory method.