Intrepid2
Intrepid2_DerivedBasis_HDIV_QUAD.hpp
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_HIV_QUAD_h
54#define Intrepid2_DerivedBasis_HIV_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 LineGradBasis = HGRAD_LINE;
80 using LineHVolBasis = HVOL_LINE;
81
82 using BasisBase = typename HGRAD_LINE::BasisBase;
83
85 public:
91 Basis_Derived_HDIV_Family1_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType)
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_HDIV;
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 DIV = Intrepid2::OPERATOR_DIV;
107
108 // ESEAS implements H(div) as rotated H(curl), which involves weighting family1 with -1.
109 // We follow that here to simplify verification tests that involve ESEAS.
110 const double weight = -1.0;
111 if (operatorType == VALUE)
112 {
113 std::vector< std::vector<EOperator> > ops(2);
114 ops[0] = std::vector<EOperator>{};
115 ops[1] = std::vector<EOperator>{VALUE,VALUE};
116 std::vector<double> weights {0.0,weight};
117 return OperatorTensorDecomposition(ops, weights);
118 }
119 else if (operatorType == DIV)
120 {
121 // family 1 is nonzero in the y component, so the div is (VALUE,GRAD)
122 std::vector< std::vector<EOperator> > ops(1); // scalar value
123 ops[0] = std::vector<EOperator>{VALUE,GRAD};
124 std::vector<double> weights {weight};
125 return OperatorTensorDecomposition(ops,weights);
126 }
127 else
128 {
129 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
130 }
131 }
132
134
142 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
143 const PointViewType inputPoints1, const PointViewType inputPoints2,
144 bool tensorPoints) const override
145 {
146 // ESEAS implements H(div) as rotated H(curl), which involves weighting family1 with -1.
147 // We follow that here to simplify verification tests that involve ESEAS.
148 const double weight = -1.0;
149
150 Intrepid2::EOperator op1, op2;
151 if (operatorType == Intrepid2::OPERATOR_VALUE)
152 {
153 op1 = Intrepid2::OPERATOR_VALUE;
154 op2 = Intrepid2::OPERATOR_VALUE;
155
156 // family 1 goes in the y component; 0 in the x component
157 auto outputValuesComponent1 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),0);
158 auto outputValuesComponent2 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),1);
159
160 this->TensorBasis::getValues(outputValuesComponent2,
161 inputPoints1, op1,
162 inputPoints2, op2, tensorPoints, weight);
163 // place 0 in the y component
164 Kokkos::deep_copy(outputValuesComponent1,0.0);
165 }
166 else if (operatorType == Intrepid2::OPERATOR_DIV)
167 {
168 // family 1 gets a d/dy applied to the second (nonzero) vector component
169 // since this is H(VOL)(x) * H(GRAD)(y), this amounts to taking the derivative in the second tensorial component
170 op1 = Intrepid2::OPERATOR_VALUE;
171 op2 = Intrepid2::OPERATOR_GRAD;
172
173 this->TensorBasis::getValues(outputValues,
174 inputPoints1, op1,
175 inputPoints2, op2, tensorPoints, weight);
176 }
177 else
178 {
179 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
180 }
181 }
182
194 virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
195 auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
196 auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
197 Kokkos::deep_copy(dofCoeffs1,0.0);
198 this->TensorBasis::getDofCoeffs(dofCoeffs2);
199 //multiply by weight = -1.0
200 Kokkos::parallel_for( Kokkos::RangePolicy<ExecutionSpace>(0, dofCoeffs2.extent(0)),
201 KOKKOS_LAMBDA (const int i){ dofCoeffs2(i) *= -1.0; });
202 }
203
204 };
205
206 template<class HGRAD_LINE, class HVOL_LINE>
208 : public Basis_TensorBasis<typename HGRAD_LINE::BasisBase>
209 {
210 using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
211 using OutputValueType = typename HGRAD_LINE::OutputValueType;
212 using PointValueType = typename HGRAD_LINE::PointValueType;
213
214 using OutputViewType = typename HGRAD_LINE::OutputViewType;
215 using PointViewType = typename HGRAD_LINE::PointViewType ;
216 using ScalarViewType = typename HGRAD_LINE::ScalarViewType;
217
218 using LineGradBasis = HGRAD_LINE;
219 using LineHVolBasis = HVOL_LINE;
220
221 using BasisBase = typename HGRAD_LINE::BasisBase;
222
224 public:
230 Basis_Derived_HDIV_Family2_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType = POINTTYPE_DEFAULT)
231 :
232 TensorBasis(Teuchos::rcp(new LineGradBasis(polyOrder_x,pointType)),
233 Teuchos::rcp(new LineHVolBasis(polyOrder_y-1,pointType)))
234 {
235 this->functionSpace_ = FUNCTION_SPACE_HDIV;
236 this->setShardsTopologyAndTags();
237 }
238
241 virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const override
242 {
243 const EOperator VALUE = Intrepid2::OPERATOR_VALUE;
244 const EOperator GRAD = Intrepid2::OPERATOR_GRAD;
245 const EOperator DIV = Intrepid2::OPERATOR_DIV;
246
247 const double weight = 1.0; // family 2 (x component nonzero)
248 if (operatorType == VALUE)
249 {
250 std::vector< std::vector<EOperator> > ops(2);
251 ops[0] = std::vector<EOperator>{VALUE,VALUE};
252 ops[1] = std::vector<EOperator>{};
253 std::vector<double> weights {weight, 0.0};
254 return OperatorTensorDecomposition(ops, weights);
255 }
256 else if (operatorType == DIV)
257 {
258 // family 2 is nonzero in the x component, so the div is (GRAD,VALUE)
259 std::vector< std::vector<EOperator> > ops(1); // scalar value
260 ops[0] = std::vector<EOperator>{GRAD,VALUE};
261 std::vector<double> weights {weight};
262 return OperatorTensorDecomposition(ops,weights);
263 }
264 else
265 {
266 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
267 }
268 }
269
271
279 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
280 const PointViewType inputPoints1, const PointViewType inputPoints2,
281 bool tensorPoints) const override
282 {
283 Intrepid2::EOperator op1, op2;
284 if (operatorType == Intrepid2::OPERATOR_VALUE)
285 {
286 op1 = Intrepid2::OPERATOR_VALUE;
287 op2 = Intrepid2::OPERATOR_VALUE;
288
289 // family 2 goes in the x component; 0 in the x component
290 auto outputValuesComponent1 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),0);
291 auto outputValuesComponent2 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),1);
292
293 this->TensorBasis::getValues(outputValuesComponent1,
294 inputPoints1, op1,
295 inputPoints2, op2, tensorPoints);
296 // place 0 in the y component
297 Kokkos::deep_copy(outputValuesComponent2, 0.0);
298 }
299 else if (operatorType == Intrepid2::OPERATOR_DIV)
300 {
301 // family 2 gets a d/dx applied to the first (nonzero) vector component
302 // since this is H(GRAD)(x) * H(VOL)(y), this amounts to taking the derivative in the first tensorial component
303 op1 = Intrepid2::OPERATOR_GRAD;
304 op2 = Intrepid2::OPERATOR_VALUE;
305
306 this->TensorBasis::getValues(outputValues,
307 inputPoints1, op1,
308 inputPoints2, op2, tensorPoints);
309 }
310 else
311 {
312 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
313 }
314 }
315
327 virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
328 auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
329 auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
330 this->TensorBasis::getDofCoeffs(dofCoeffs1);
331 Kokkos::deep_copy(dofCoeffs2,0.0);
332 }
333
334 };
335
336 template<class HGRAD_LINE, class HVOL_LINE>
338 : public Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>
339 {
342 using DirectSumBasis = Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>;
343
344 public:
345 using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
346 using OutputValueType = typename HGRAD_LINE::OutputValueType;
347 using PointValueType = typename HGRAD_LINE::PointValueType;
348
349 using BasisBase = typename HGRAD_LINE::BasisBase;
350
351 protected:
352 std::string name_;
353 ordinal_type order_x_;
354 ordinal_type order_y_;
355 EPointType pointType_;
356 public:
362 Basis_Derived_HDIV_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT)
363 :
364 DirectSumBasis(Teuchos::rcp(new Family1(polyOrder_x, polyOrder_y, pointType)),
365 Teuchos::rcp(new Family2(polyOrder_x, polyOrder_y, pointType)))
366 {
367 this->functionSpace_ = FUNCTION_SPACE_HDIV;
368
369 std::ostringstream basisName;
370 basisName << "HDIV_QUAD (" << this->DirectSumBasis::getName() << ")";
371 name_ = basisName.str();
372
373 order_x_ = polyOrder_x;
374 order_y_ = polyOrder_y;
375 pointType_ = pointType;
376 }
377
382 Basis_Derived_HDIV_QUAD(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) : Basis_Derived_HDIV_QUAD(polyOrder, polyOrder, pointType) {}
383
386 virtual bool requireOrientation() const override {
387 return true;
388 }
389
394 virtual
395 const char*
396 getName() const override {
397 return name_.c_str();
398 }
399
410 Teuchos::RCP<BasisBase>
411 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
412 if(subCellDim == 1) {
413 switch(subCellOrd) {
414 case 0:
415 case 2:
416 return Teuchos::rcp( new HVOL_LINE(order_x_-1, pointType_) );
417 case 1:
418 case 3:
419 return Teuchos::rcp( new HVOL_LINE(order_y_-1, pointType_) );
420 }
421 }
422
423 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
424 }
425
431 getHostBasis() const override {
433
434 auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, pointType_));
435
436 return hostBasis;
437 }
438 };
439} // end namespace Intrepid2
440
441#endif /* Intrepid2_DerivedBasis_HIV_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.
Basis_Derived_HDIV_Family1_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType)
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 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_HDIV_Family2_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...
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 bool requireOrientation() const override
True if orientation is required.
Teuchos::RCP< BasisBase > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
Basis_Derived_HDIV_QUAD(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual const char * getName() const override
Returns basis name.
Basis_Derived_HDIV_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...