Intrepid2
Intrepid2_DerivedBasis_HGRAD_HEX.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),
38 // Mauro Perego (mperego@sandia.gov), or
39 // Nate Roberts (nvrober@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
57 #ifndef Intrepid2_DerivedBasis_HGRAD_HEX_h
58 #define Intrepid2_DerivedBasis_HGRAD_HEX_h
59 
61 
63 
64 namespace Intrepid2
65 {
66  // TODO: make this a subclass of TensorBasis3 instead, following what we've done for H(curl) and H(div)
67  template<class HGRAD_LINE>
69  : public Basis_TensorBasis<typename HGRAD_LINE::BasisBase>
70  {
71  public:
72  using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
73  using OutputValueType = typename HGRAD_LINE::OutputValueType;
74  using PointValueType = typename HGRAD_LINE::PointValueType;
75 
76  using OutputViewType = typename HGRAD_LINE::OutputViewType;
77  using PointViewType = typename HGRAD_LINE::PointViewType ;
78  using ScalarViewType = typename HGRAD_LINE::ScalarViewType;
79 
80  using LineBasis = HGRAD_LINE;
82  using BasisBase = typename HGRAD_LINE::BasisBase;
84 
85  std::string name_;
86  ordinal_type order_x_;
87  ordinal_type order_y_;
88  ordinal_type order_z_;
89  EPointType pointType_;
90 
97  Basis_Derived_HGRAD_HEX(int polyOrder_x, int polyOrder_y, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
98  :
99  TensorBasis(Teuchos::rcp( new QuadBasis(polyOrder_x,polyOrder_y, pointType)),
100  Teuchos::rcp( new LineBasis(polyOrder_z, pointType)))
101  {
102  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
103 
104  std::ostringstream basisName;
105  basisName << "HGRAD_HEX (" << this->TensorBasis::getName() << ")";
106  name_ = basisName.str();
107 
108  order_x_ = polyOrder_x;
109  order_y_ = polyOrder_y;
110  order_z_ = polyOrder_z;
111  pointType_ = pointType;
112  }
113 
114 
119  Basis_Derived_HGRAD_HEX(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) :
120  Basis_Derived_HGRAD_HEX(polyOrder, polyOrder, polyOrder, pointType) {}
121 
122 
125  virtual bool requireOrientation() const override {
126  return (this->getDofCount(1,0) > 1); //if it has more than 1 DOF per edge, than it needs orientations
127  }
128 
131  virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator operatorType) const override
132  {
133  const EOperator VALUE = Intrepid2::OPERATOR_VALUE;
134  const EOperator GRAD = Intrepid2::OPERATOR_GRAD;
135 
136  if (operatorType == VALUE)
137  {
138  return OperatorTensorDecomposition(VALUE,VALUE);
139  }
140  else if (operatorType == GRAD)
141  {
142  // to evaluate gradient, we need both OP_VALUE and OP_GRAD (thanks to product rule)
143  // for quad x line, we will put derivative * value in first component, and value * derivative in second
144 
145  std::vector< std::vector<EOperator> > ops;
146  ops.push_back(std::vector<EOperator>{GRAD, VALUE});
147  ops.push_back(std::vector<EOperator>{VALUE, GRAD});
148 
149  std::vector<double> weights(ops.size(), 1.0);
150 
151  return OperatorTensorDecomposition(ops, weights);
152  }
153  else
154  {
155  INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
156  }
157  }
158 
159  using BasisBase::getValues;
160 
168  virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
169  const PointViewType inputPoints1, const PointViewType inputPoints2,
170  bool tensorPoints) const override
171  {
172  Intrepid2::EOperator op1, op2;
173  if (operatorType == Intrepid2::OPERATOR_VALUE)
174  {
175  op1 = Intrepid2::OPERATOR_VALUE;
176  op2 = Intrepid2::OPERATOR_VALUE;
177 
178  this->TensorBasis::getValues(outputValues,
179  inputPoints1, op1,
180  inputPoints2, op2, tensorPoints);
181  }
182  else if (operatorType == Intrepid2::OPERATOR_GRAD)
183  {
184  // to evaluate gradient, we actually need both OP_VALUE and OP_GRAD (thanks to product rule)
185  // for 1D line x line, we will put derivative * value in first component, and value * derivative in second
186 
187  // outputValues1 and outputValues2 are computed by basis1 and basis2 -- these are tensorial components
188  // outputValuesComponent1 is a slice of the final output container (similarly, outputValuesComponent2)
189  // when the component basis is 1D, it expects not to have a "dimension" component in the output container
190  // the int argument in the dimension component creates a subview that skips the dimension component; the std::pair argument retains it
191  auto outputValuesComponent1 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::make_pair(0,2));
192  auto outputValuesComponent2 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),2);
193 
194  // compute first component -- derivative happens in x and y, and value taken in z
195  op1 = Intrepid2::OPERATOR_GRAD;
196  op2 = Intrepid2::OPERATOR_VALUE;
197 
198  this->TensorBasis::getValues(outputValuesComponent1,
199  inputPoints1, op1,
200  inputPoints2, op2, tensorPoints);
201 
202  // second component -- value in x and y, derivative in z
203  op1 = Intrepid2::OPERATOR_VALUE;
204  op2 = Intrepid2::OPERATOR_GRAD;
205 
206  this->TensorBasis::getValues(outputValuesComponent2,
207  inputPoints1, op1,
208  inputPoints2, op2, tensorPoints);
209  }
210  else
211  {
212  INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
213  }
214  }
215 
220  virtual
221  const char*
222  getName() const override {
223  return name_.c_str();
224  }
225 
235  Teuchos::RCP<BasisBase>
236  getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
237  if(subCellDim == 1) {
238  switch(subCellOrd) {
239  case 0:
240  case 2:
241  case 4:
242  case 6:
243  return Teuchos::rcp( new LineBasis(order_x_, pointType_) );
244  case 1:
245  case 3:
246  case 5:
247  case 7:
248  return Teuchos::rcp( new LineBasis(order_y_, pointType_) );
249  case 8:
250  case 9:
251  case 10:
252  case 11:
253  return Teuchos::rcp( new LineBasis(order_z_, pointType_) );
254  }
255  } else if(subCellDim == 2) {
256  switch(subCellOrd) {
257  case 0:
258  return Teuchos::rcp( new QuadBasis(order_x_, order_z_, pointType_) );
259  case 1:
260  return Teuchos::rcp( new QuadBasis(order_y_,order_z_, pointType_) );
261  case 2:
262  return Teuchos::rcp( new QuadBasis(order_x_, order_z_, pointType_) );
263  case 3:
264  return Teuchos::rcp( new QuadBasis(order_z_, order_y_, pointType_) );
265  case 4:
266  return Teuchos::rcp( new QuadBasis(order_y_, order_x_, pointType_) );
267  case 5:
268  return Teuchos::rcp( new QuadBasis(order_x_, order_y_, pointType_) );
269  }
270  }
271 
272  INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
273  }
274 
280  getHostBasis() const override {
282 
283  auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, order_z_, pointType_));
284 
285  return hostBasis;
286  }
287  };
288 } // end namespace Intrepid2
289 
290 #endif /* Intrepid2_DerivedBasis_HGRAD_HEX_h */
Implementation of bases that are tensor products of two or three component bases. ...
Teuchos::RCP< BasisBase > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
For a multi-component tensor basis, specifies the operators to be applied to the components to produc...
BasisPtr< typename Kokkos::HostSpace::device_type, OutputType, PointType > HostBasisPtr
Pointer to a Basis whose device type is on the host (Kokkos::HostSpace::device_type), allowing host access to input and output views, and ensuring host execution of basis evaluation.
Implementation of H(grad) basis on the quadrilateral that is templated on H(grad) on the line...
virtual bool requireOrientation() const override
True if orientation is required.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
virtual const char * getName() const override
Returns basis name.
EPointType
Enumeration of types of point distributions in Intrepid.
virtual void getValues(OutputViewType outputValues, const EOperator operatorType, const PointViewType inputPoints1, const PointViewType inputPoints2, bool tensorPoints) const override
multi-component getValues() method (required/called by TensorBasis)
virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator operatorType) const override
Returns a simple decomposition of the specified operator: what operator(s) should be applied to basis...
virtual void getValues(BasisValues< OutputValueType, DeviceType > outputValues, const TensorPoints< PointValueType, DeviceType > inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell, using point and output value containers that allow pre...
virtual const char * getName() const override
Returns basis name.
Basis_Derived_HGRAD_HEX(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Basis defined as the tensor product of two component bases.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
Basis_Derived_HGRAD_HEX(int polyOrder_x, int polyOrder_y, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.