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