Intrepid2
Intrepid2_DerivedBasis_HCURL_QUAD.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 
53 #ifndef Intrepid2_DerivedBasis_HCURL_QUAD_h
54 #define Intrepid2_DerivedBasis_HCURL_QUAD_h
55 
56 #include <Kokkos_View.hpp>
57 #include <Kokkos_DynRankView.hpp>
58 
60 #include "Intrepid2_Sacado.hpp"
61 
64 
65 namespace Intrepid2
66 {
67  template<class HGRAD_LINE, class HVOL_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 BasisBase = typename HGRAD_LINE::BasisBase;
81 
82  using LineGradBasis = HGRAD_LINE;
83  using LineHVolBasis = HVOL_LINE;
84 
86  public:
92  Basis_Derived_HCURL_Family1_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType = POINTTYPE_DEFAULT)
93  :
94  TensorBasis(Teuchos::rcp( new LineHVolBasis(polyOrder_x-1,pointType)),
95  Teuchos::rcp( new LineGradBasis(polyOrder_y,pointType)))
96  {
97  this->functionSpace_ = FUNCTION_SPACE_HCURL;
98  }
99 
102  virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator operatorType) const override
103  {
104  const EOperator VALUE = Intrepid2::OPERATOR_VALUE;
105  const EOperator GRAD = Intrepid2::OPERATOR_GRAD;
106  const EOperator CURL = Intrepid2::OPERATOR_CURL;
107  if (operatorType == VALUE)
108  {
109  // family 1 goes in x component
110  std::vector< std::vector<EOperator> > ops(2);
111  ops[0] = std::vector<EOperator>{VALUE,VALUE};
112  ops[1] = std::vector<EOperator>{};
113  std::vector<double> weights {1.0, 0.0};
114  return OperatorTensorDecomposition(ops, weights);
115  }
116  else if (operatorType == CURL)
117  {
118  // family 1 gets a -d/dy applied to the first (nonzero) vector component
119  // since this is H(VOL)(x) * H(GRAD)(y), this amounts to taking the derivative in the second tensorial component
120  return OperatorTensorDecomposition(VALUE,GRAD,-1.0);
121  }
122  else
123  {
124  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
125  }
126  }
127 
129 
137  virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
138  const PointViewType inputPoints1, const PointViewType inputPoints2,
139  bool tensorPoints) const override
140  {
141  Intrepid2::EOperator op1, op2;
142  if (operatorType == Intrepid2::OPERATOR_VALUE)
143  {
144  op1 = Intrepid2::OPERATOR_VALUE;
145  op2 = Intrepid2::OPERATOR_VALUE;
146 
147  // family 1 goes in the x component; 0 in the y component
148  OutputViewType outputValuesComponent1 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),0);
149  OutputViewType outputValuesComponent2 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),1);
150 
151  this->TensorBasis::getValues(outputValuesComponent1,
152  inputPoints1, op1,
153  inputPoints2, op2, tensorPoints);
154  // place 0 in the y component
155  Kokkos::deep_copy(outputValuesComponent2,0);
156  }
157  else if (operatorType == Intrepid2::OPERATOR_CURL)
158  {
159  // family 1 gets a -d/dy applied to the first (nonzero) vector component
160  // since this is H(VOL)(x) * H(GRAD)(y), this amounts to taking the derivative in the second tensorial component
161  op1 = Intrepid2::OPERATOR_VALUE;
162  op2 = Intrepid2::OPERATOR_GRAD;
163 
164  double weight = -1.0; // the minus sign in front of d/dy
165  this->TensorBasis::getValues(outputValues,
166  inputPoints1, op1,
167  inputPoints2, op2, tensorPoints, weight);
168  }
169  else
170  {
171  INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
172  }
173  }
174 
186  virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
187  auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
188  auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
189  this->TensorBasis::getDofCoeffs(dofCoeffs1);
190  Kokkos::deep_copy(dofCoeffs2,0.0);
191  }
192  };
193 
194  template<class HGRAD_LINE, class HVOL_LINE>
196  : public Basis_TensorBasis<typename HGRAD_LINE::BasisBase>
197  {
198 
199  public:
200  using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
201  using OutputValueType = typename HGRAD_LINE::OutputValueType;
202  using PointValueType = typename HGRAD_LINE::PointValueType;
203 
204  using OutputViewType = typename HGRAD_LINE::OutputViewType;
205  using PointViewType = typename HGRAD_LINE::PointViewType ;
206  using ScalarViewType = typename HGRAD_LINE::ScalarViewType;
207 
208  using LineGradBasis = HGRAD_LINE;
209  using LineHVolBasis = HVOL_LINE;
210 
211  using BasisBase = typename HGRAD_LINE::BasisBase;
212 
214 
220  Basis_Derived_HCURL_Family2_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType = POINTTYPE_DEFAULT)
221  :
222  TensorBasis(Teuchos::rcp( new LineGradBasis(polyOrder_x,pointType) ),
223  Teuchos::rcp( new LineHVolBasis(polyOrder_y-1,pointType) ))
224  {
225  this->functionSpace_ = FUNCTION_SPACE_HCURL;
226  }
227 
230  virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator operatorType) const override
231  {
232  const EOperator VALUE = Intrepid2::OPERATOR_VALUE;
233  const EOperator GRAD = Intrepid2::OPERATOR_GRAD;
234  const EOperator CURL = Intrepid2::OPERATOR_CURL;
235  if (operatorType == VALUE)
236  {
237  // family 2 goes in y component
238  std::vector< std::vector<EOperator> > ops(2);
239  ops[0] = std::vector<EOperator>{};
240  ops[1] = std::vector<EOperator>{VALUE,VALUE};
241  std::vector<double> weights {0.0, 1.0};
242  return OperatorTensorDecomposition(ops, weights);
243  }
244  else if (operatorType == CURL)
245  {
246  // family 2 gets a d/dx applied to the second (nonzero) vector component
247  // since this is H(GRAD)(x) * H(VOL)(y), this amounts to taking the derivative in the first tensorial component
248  return OperatorTensorDecomposition(GRAD,VALUE,1.0);
249  }
250  else
251  {
252  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
253  }
254  }
255 
257 
265  virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
266  const PointViewType inputPoints1, const PointViewType inputPoints2,
267  bool tensorPoints) const override
268  {
269  Intrepid2::EOperator op1, op2;
270  if (operatorType == Intrepid2::OPERATOR_VALUE)
271  {
272  op1 = Intrepid2::OPERATOR_VALUE;
273  op2 = Intrepid2::OPERATOR_VALUE;
274 
275  // family 2 goes in the y component; 0 in the x component
276  auto outputValuesComponent1 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),0);
277  auto outputValuesComponent2 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),1);
278 
279  // place 0 in the x component
280  Kokkos::deep_copy(outputValuesComponent1, 0.0);
281  this->TensorBasis::getValues(outputValuesComponent2,
282  inputPoints1, op1,
283  inputPoints2, op2, tensorPoints);
284 
285  }
286  else if (operatorType == Intrepid2::OPERATOR_CURL)
287  {
288  // family 2 gets a d/dx applied to the second (nonzero) vector component
289  // since this is H(GRAD)(x) * H(VOL)(y), this amounts to taking the derivative in the first tensorial component
290  op1 = Intrepid2::OPERATOR_GRAD;
291  op2 = Intrepid2::OPERATOR_VALUE;
292 
293  this->TensorBasis::getValues(outputValues,
294  inputPoints1, op1,
295  inputPoints2, op2, tensorPoints);
296  }
297  else
298  {
299  INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
300  }
301  }
302 
314  virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
315  auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
316  auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
317  Kokkos::deep_copy(dofCoeffs1,0.0);
318  this->TensorBasis::getDofCoeffs(dofCoeffs2);
319  }
320  };
321 
322  template<class HGRAD_LINE, class HVOL_LINE>
324  : public Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>
325  {
329 
330  using BasisBase = typename HGRAD_LINE::BasisBase;
331 
332  protected:
333  std::string name_;
334  ordinal_type order_x_;
335  ordinal_type order_y_;
336  EPointType pointType_;
337 
338  public:
339  using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
340  using OutputValueType = typename HGRAD_LINE::OutputValueType;
341  using PointValueType = typename HGRAD_LINE::PointValueType;
342 
348  Basis_Derived_HCURL_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
349  :
350  DirectSumBasis(Teuchos::rcp( new Family1(polyOrder_x, polyOrder_y, pointType) ),
351  Teuchos::rcp( new Family2(polyOrder_x, polyOrder_y, pointType) ))
352  {
353  this->functionSpace_ = FUNCTION_SPACE_HCURL;
354 
355  std::ostringstream basisName;
356  basisName << "HCURL_QUAD (" << this->DirectSumBasis::getName() << ")";
357  name_ = basisName.str();
358 
359  order_x_ = polyOrder_x;
360  order_y_ = polyOrder_y;
361  pointType_ = pointType;
362  }
363 
368  Basis_Derived_HCURL_QUAD(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) : Basis_Derived_HCURL_QUAD(polyOrder, polyOrder, pointType) {}
369 
372  virtual bool requireOrientation() const override
373  {
374  return (this->getDofCount(1,0) > 0); //if it has edge DOFs, than it needs orientations
375  }
376 
381  virtual
382  const char*
383  getName() const override {
384  return name_.c_str();
385  }
386 
397  Teuchos::RCP<BasisBase>
398  getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
399  if(subCellDim == 1) {
400  switch(subCellOrd) {
401  case 0:
402  case 2:
403  return Teuchos::rcp( new HVOL_LINE(order_x_-1, pointType_) );
404  case 1:
405  case 3:
406  return Teuchos::rcp( new HVOL_LINE(order_y_-1, pointType_) );
407  }
408  }
409 
410  INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
411  }
412 
418  getHostBasis() const override {
420 
421  auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, pointType_));
422 
423  return hostBasis;
424  }
425  };
426 } // end namespace Intrepid2
427 
428 #endif /* Intrepid2_DerivedBasis_HCURL_QUAD_h */
Basis_Derived_HCURL_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator operatorType) const override
Returns a simple decomposition of the specified operator: what operator(s) should be applied to basis...
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.
virtual void getDofCoeffs(ScalarViewType dofCoeffs) const override
Fills in coefficients of degrees of freedom for Lagrangian basis on the reference cell...
virtual void getDofCoeffs(ScalarViewType dofCoeffs) const override
Fills in coefficients of degrees of freedom for Lagrangian basis on the reference cell...
Free functions, callable from device code, that implement various polynomials useful in basis definit...
virtual bool requireOrientation() const override
True if orientation is required.
For a multi-component tensor basis, specifies the operators to be applied to the components to produc...
Implementation of a basis that is the direct sum of two other bases.
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.
Basis_Derived_HCURL_Family1_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual void getDofCoeffs(typename BasisBase::ScalarViewType dofCoeffs) const override
Fills in coefficients of degrees of freedom on the reference cell.
A basis that is the direct sum of two other bases.
virtual const char * getName() const override
Returns basis name.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
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)
EPointType
Enumeration of types of point distributions in Intrepid.
Header file to include all Sacado headers that are required if using Intrepid2 with Sacado types...
virtual const char * getName() const override
Returns basis name.
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 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)
Basis defined as the tensor product of two component bases.
Basis_Derived_HCURL_Family2_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Basis_Derived_HCURL_QUAD(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.