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_DynRankView.hpp>
57
59#include "Intrepid2_Sacado.hpp"
60
63
64namespace Intrepid2
65{
66 template<class HGRAD_LINE, class HVOL_LINE>
68 : public Basis_TensorBasis<typename HGRAD_LINE::BasisBase>
69 {
70 public:
71 using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
72 using OutputValueType = typename HGRAD_LINE::OutputValueType;
73 using PointValueType = typename HGRAD_LINE::PointValueType;
74
75 using OutputViewType = typename HGRAD_LINE::OutputViewType;
76 using PointViewType = typename HGRAD_LINE::PointViewType ;
77 using ScalarViewType = typename HGRAD_LINE::ScalarViewType;
78
79 using BasisBase = typename HGRAD_LINE::BasisBase;
80
81 using LineGradBasis = HGRAD_LINE;
82 using LineHVolBasis = HVOL_LINE;
83
85 public:
91 Basis_Derived_HCURL_Family1_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType = POINTTYPE_DEFAULT)
92 :
93 TensorBasis(Teuchos::rcp( new LineHVolBasis(polyOrder_x-1,pointType)),
94 Teuchos::rcp( new LineGradBasis(polyOrder_y,pointType)))
95 {
96 this->functionSpace_ = FUNCTION_SPACE_HCURL;
97 this->setShardsTopologyAndTags();
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 this->setShardsTopologyAndTags();
227 }
228
231 virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const override
232 {
233 const EOperator VALUE = Intrepid2::OPERATOR_VALUE;
234 const EOperator GRAD = Intrepid2::OPERATOR_GRAD;
235 const EOperator CURL = Intrepid2::OPERATOR_CURL;
236 if (operatorType == VALUE)
237 {
238 // family 2 goes in y component
239 std::vector< std::vector<EOperator> > ops(2);
240 ops[0] = std::vector<EOperator>{};
241 ops[1] = std::vector<EOperator>{VALUE,VALUE};
242 std::vector<double> weights {0.0, 1.0};
243 return OperatorTensorDecomposition(ops, weights);
244 }
245 else if (operatorType == CURL)
246 {
247 // family 2 gets a d/dx applied to the second (nonzero) vector component
248 // since this is H(GRAD)(x) * H(VOL)(y), this amounts to taking the derivative in the first tensorial component
249 return OperatorTensorDecomposition(GRAD,VALUE,1.0);
250 }
251 else
252 {
253 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
254 }
255 }
256
258
266 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
267 const PointViewType inputPoints1, const PointViewType inputPoints2,
268 bool tensorPoints) const override
269 {
270 Intrepid2::EOperator op1, op2;
271 if (operatorType == Intrepid2::OPERATOR_VALUE)
272 {
273 op1 = Intrepid2::OPERATOR_VALUE;
274 op2 = Intrepid2::OPERATOR_VALUE;
275
276 // family 2 goes in the y component; 0 in the x component
277 auto outputValuesComponent1 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),0);
278 auto outputValuesComponent2 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),1);
279
280 // place 0 in the x component
281 Kokkos::deep_copy(outputValuesComponent1, 0.0);
282 this->TensorBasis::getValues(outputValuesComponent2,
283 inputPoints1, op1,
284 inputPoints2, op2, tensorPoints);
285
286 }
287 else if (operatorType == Intrepid2::OPERATOR_CURL)
288 {
289 // family 2 gets a d/dx applied to the second (nonzero) vector component
290 // since this is H(GRAD)(x) * H(VOL)(y), this amounts to taking the derivative in the first tensorial component
291 op1 = Intrepid2::OPERATOR_GRAD;
292 op2 = Intrepid2::OPERATOR_VALUE;
293
294 this->TensorBasis::getValues(outputValues,
295 inputPoints1, op1,
296 inputPoints2, op2, tensorPoints);
297 }
298 else
299 {
300 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
301 }
302 }
303
315 virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
316 auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
317 auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
318 Kokkos::deep_copy(dofCoeffs1,0.0);
319 this->TensorBasis::getDofCoeffs(dofCoeffs2);
320 }
321 };
322
323 template<class HGRAD_LINE, class HVOL_LINE>
325 : public Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>
326 {
329 using DirectSumBasis = Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>;
330 public:
331 using BasisBase = typename HGRAD_LINE::BasisBase;
332
333 protected:
334 std::string name_;
335 ordinal_type order_x_;
336 ordinal_type order_y_;
337 EPointType pointType_;
338
339 public:
340 using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
341 using OutputValueType = typename HGRAD_LINE::OutputValueType;
342 using PointValueType = typename HGRAD_LINE::PointValueType;
343
349 Basis_Derived_HCURL_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
350 :
351 DirectSumBasis(Teuchos::rcp( new Family1(polyOrder_x, polyOrder_y, pointType) ),
352 Teuchos::rcp( new Family2(polyOrder_x, polyOrder_y, pointType) ))
353 {
354 this->functionSpace_ = FUNCTION_SPACE_HCURL;
355
356 std::ostringstream basisName;
357 basisName << "HCURL_QUAD (" << this->DirectSumBasis::getName() << ")";
358 name_ = basisName.str();
359
360 order_x_ = polyOrder_x;
361 order_y_ = polyOrder_y;
362 pointType_ = pointType;
363 }
364
369 Basis_Derived_HCURL_QUAD(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) : Basis_Derived_HCURL_QUAD(polyOrder, polyOrder, pointType) {}
370
373 virtual bool requireOrientation() const override
374 {
375 return true;
376 }
377
382 virtual
383 const char*
384 getName() const override {
385 return name_.c_str();
386 }
387
398 Teuchos::RCP<BasisBase>
399 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
400 if(subCellDim == 1) {
401 switch(subCellOrd) {
402 case 0:
403 case 2:
404 return Teuchos::rcp( new HVOL_LINE(order_x_-1, pointType_) );
405 case 1:
406 case 3:
407 return Teuchos::rcp( new HVOL_LINE(order_y_-1, pointType_) );
408 }
409 }
410
411 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
412 }
413
419 getHostBasis() const override {
421
422 auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, pointType_));
423
424 return hostBasis;
425 }
426 };
427} // end namespace Intrepid2
428
429#endif /* Intrepid2_DerivedBasis_HCURL_QUAD_h */
BasisPtr< typename Kokkos::HostSpace::device_type, OutputType, PointType > HostBasisPtr
Pointer to a Basis whose device type is on the host (Kokkos::HostSpace::device_type),...
Implementation of a basis that is the direct sum of two other bases.
Free functions, callable from device code, that implement various polynomials useful in basis definit...
Header file to include all Sacado headers that are required if using Intrepid2 with Sacado types.
Implementation of bases that are tensor products of two or three component bases.
virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const override
Returns a simple decomposition of the specified operator: what operator(s) should be applied to basis...
Basis_Derived_HCURL_Family1_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
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 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.
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...
Basis_Derived_HCURL_Family2_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual const char * getName() const override
Returns basis name.
Teuchos::RCP< BasisBase > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
Basis_Derived_HCURL_QUAD(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual bool requireOrientation() const override
True if orientation is required.
Basis_Derived_HCURL_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
A basis that is the direct sum of two other bases.
virtual const char * getName() const override
Returns basis name.
Basis defined as the tensor product of two component bases.
virtual void getDofCoeffs(typename BasisBase::ScalarViewType dofCoeffs) const override
Fills in coefficients of degrees of freedom on the reference cell.
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...
For a multi-component tensor basis, specifies the operators to be applied to the components to produc...