Intrepid2
Intrepid2_DerivedBasis_HDIV_WEDGE.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
64#ifndef Intrepid2_DerivedBasis_HDIV_WEDGE_h
65#define Intrepid2_DerivedBasis_HDIV_WEDGE_h
66
67#include <Kokkos_DynRankView.hpp>
68
70#include "Intrepid2_Sacado.hpp"
71
74
75namespace Intrepid2
76{
77 template<class HDIV_TRI, class HVOL_LINE>
79 : public Basis_TensorBasis<typename HVOL_LINE::BasisBase>
80 {
81
82 public:
83 using OutputViewType = typename HVOL_LINE::OutputViewType;
84 using PointViewType = typename HVOL_LINE::PointViewType ;
85 using ScalarViewType = typename HVOL_LINE::ScalarViewType;
86
87 using TriDivBasis = HDIV_TRI;
88 using LineHVolBasis = HVOL_LINE;
89
90 using BasisBase = typename HVOL_LINE::BasisBase;
92
93 using DeviceType = typename BasisBase::DeviceType;
94 using ExecutionSpace = typename BasisBase::ExecutionSpace;
95 using OutputValueType = typename BasisBase::OutputValueType;
96 using PointValueType = typename BasisBase::PointValueType;
97
103 Basis_Derived_HDIV_Family1_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType = POINTTYPE_DEFAULT)
104 :
105 TensorBasis(Teuchos::rcp( new TriDivBasis(polyOrder_xy, pointType) ),
106 Teuchos::rcp( new LineHVolBasis(polyOrder_z-1, pointType) ))
107 {
108 this->functionSpace_ = FUNCTION_SPACE_HDIV;
109 this->setShardsTopologyAndTags();
110 }
111
114 OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const override
115 {
116 if (operatorType == OPERATOR_DIV)
117 {
118 // div of (f(x,y)h(z),g(x,y)h(z),0) is (df/dx + dg/dy) * h
119
120 std::vector< std::vector<EOperator> > ops(1);
121 ops[0] = std::vector<EOperator>{OPERATOR_DIV,OPERATOR_VALUE};
122 std::vector<double> weights {1.0};
123 OperatorTensorDecomposition opDecomposition(ops, weights);
124 return opDecomposition;
125 }
126 else if (OPERATOR_VALUE == operatorType)
127 {
128 // family 1 goes in x,y components
129 std::vector< std::vector<EOperator> > ops(2);
130 ops[0] = std::vector<EOperator>{OPERATOR_VALUE,OPERATOR_VALUE}; // spans (x,y) due to H(div,tri) (first tensorial component)
131 ops[1] = std::vector<EOperator>{}; // empty z component
132 std::vector<double> weights {1.0,0.0};
133 return OperatorTensorDecomposition(ops, weights);
134 }
135 else
136 {
137 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
138 }
139 }
140
142
150 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
151 const PointViewType inputPoints1, const PointViewType inputPoints2,
152 bool tensorPoints) const override
153 {
154 EOperator op1, op2;
155 if (operatorType == OPERATOR_VALUE)
156 {
157 op1 = OPERATOR_VALUE;
158 op2 = OPERATOR_VALUE;
159
160 // family 1 values goes in (x,y) components, 0 in z
161 auto outputValuesComponent12 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::pair<int,int>{0,2});
162 auto outputValuesComponent3 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),2);
163
164 // place 0 in the z component
165 Kokkos::deep_copy(outputValuesComponent3, 0.0);
166 this->TensorBasis::getValues(outputValuesComponent12,
167 inputPoints1, op1,
168 inputPoints2, op2, tensorPoints);
169
170 }
171 else if (operatorType == OPERATOR_DIV)
172 {
173 // div of (f(x,y)h(z),g(x,y)h(z),0) is (df/dx + dg/dy) * h
174
175 op1 = OPERATOR_DIV;
176 op2 = OPERATOR_VALUE;
177
178 this->TensorBasis::getValues(outputValues,
179 inputPoints1, op1,
180 inputPoints2, op2, tensorPoints);
181 }
182 else
183 {
184 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
185 }
186 }
187
199 virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
200 auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
201 auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
202 Kokkos::deep_copy(dofCoeffs1,0.0);
203 this->TensorBasis::getDofCoeffs(dofCoeffs2);
204 }
205 };
206
207 template<class HVOL_TRI, class HGRAD_LINE>
209 : public Basis_TensorBasis<typename HGRAD_LINE::BasisBase>
210 {
211 public:
212 using OutputViewType = typename HGRAD_LINE::OutputViewType;
213 using PointViewType = typename HGRAD_LINE::PointViewType ;
214 using ScalarViewType = typename HGRAD_LINE::ScalarViewType;
215
216 using BasisBase = typename HGRAD_LINE::BasisBase;
217
218 using DeviceType = typename BasisBase::DeviceType;
219 using ExecutionSpace = typename BasisBase::ExecutionSpace;
220 using OutputValueType = typename BasisBase::OutputValueType;
221 using PointValueType = typename BasisBase::PointValueType;
222
223 using TriVolBasis = HVOL_TRI;
224 using LineGradBasis = HGRAD_LINE;
225
227 public:
233 Basis_Derived_HDIV_Family2_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType = POINTTYPE_DEFAULT)
234 :
235 TensorBasis(Teuchos::rcp( new TriVolBasis(polyOrder_xy-1,pointType)),
236 Teuchos::rcp( new LineGradBasis(polyOrder_z,pointType)))
237 {
238 this->functionSpace_ = FUNCTION_SPACE_HDIV;
239 this->setShardsTopologyAndTags();
240 }
241
243
246 virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const override
247 {
248 const EOperator & VALUE = OPERATOR_VALUE;
249 const EOperator & GRAD = OPERATOR_GRAD;
250 const EOperator & DIV = OPERATOR_DIV;
251 if (operatorType == VALUE)
252 {
253 // family 2 goes in z component
254 std::vector< std::vector<EOperator> > ops(2);
255 ops[0] = std::vector<EOperator>{}; // will span x,y because family 1's first component does
256 ops[1] = std::vector<EOperator>{VALUE,VALUE};
257 std::vector<double> weights {0.,1.0};
258 OperatorTensorDecomposition opDecomposition(ops, weights);
259 return opDecomposition;
260 }
261 else if (operatorType == DIV)
262 {
263 // div of (0, 0, f*g, where f=f(x,y), g=g(z), equals f * g'(z).
264 std::vector< std::vector<EOperator> > ops(1);
265 ops[0] = std::vector<EOperator>{VALUE,GRAD};
266 std::vector<double> weights {1.0};
267 OperatorTensorDecomposition opDecomposition(ops, weights);
268 return opDecomposition;
269 }
270 else
271 {
272 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
273 }
274 }
275
283 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
284 const PointViewType inputPoints1, const PointViewType inputPoints2,
285 bool tensorPoints) const override
286 {
287 EOperator op1, op2;
288 if (operatorType == OPERATOR_VALUE)
289 {
290 op1 = OPERATOR_VALUE;
291 op2 = OPERATOR_VALUE;
292
293 // family 2 values goes in z component
294 auto outputValuesComponent12 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::pair<int,int>{0,2});
295 auto outputValuesComponent3 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),2);
296
297 // place 0 in the x,y components
298 Kokkos::deep_copy(outputValuesComponent12, 0.0);
299 this->TensorBasis::getValues(outputValuesComponent3,
300 inputPoints1, op1,
301 inputPoints2, op2, tensorPoints);
302
303 }
304 else if (operatorType == OPERATOR_DIV)
305 {
306 // div of (0, 0, f*g, where f=f(x,y), g=g(z), equals f * g'(z).
307 op1 = OPERATOR_VALUE;
308 op2 = OPERATOR_GRAD;
309
310 this->TensorBasis::getValues(outputValues,
311 inputPoints1, op1,
312 inputPoints2, op2, tensorPoints);
313 }
314 else
315 {
316 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
317 }
318 }
319
331 virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
332 auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
333 auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
334 this->TensorBasis::getDofCoeffs(dofCoeffs1);
335 Kokkos::deep_copy(dofCoeffs2,0.0);
336 }
337 };
338
339
340 template<class HDIV_TRI, class HVOL_TRI, class HGRAD_LINE, class HVOL_LINE>
342 : public Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>
343 {
346 using DirectSumBasis = Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>;
347 public:
348 using BasisBase = typename HGRAD_LINE::BasisBase;
349
350 protected:
351 std::string name_;
352 ordinal_type order_xy_;
353 ordinal_type order_z_;
354 EPointType pointType_;
355
356 public:
357 using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
358 using OutputValueType = typename HGRAD_LINE::OutputValueType;
359 using PointValueType = typename HGRAD_LINE::PointValueType;
360
366 Basis_Derived_HDIV_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
367 :
368 DirectSumBasis(Teuchos::rcp( new Family1(polyOrder_xy, polyOrder_z, pointType) ),
369 Teuchos::rcp( new Family2(polyOrder_xy, polyOrder_z, pointType) ))
370 {
371 this->functionSpace_ = FUNCTION_SPACE_HDIV;
372
373 std::ostringstream basisName;
374 basisName << "HDIV_WEDGE (" << this->DirectSumBasis::getName() << ")";
375 name_ = basisName.str();
376
377 order_xy_ = polyOrder_xy;
378 order_z_ = polyOrder_z;
379 pointType_ = pointType;
380 }
381
386 Basis_Derived_HDIV_WEDGE(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) : Basis_Derived_HDIV_WEDGE(polyOrder, polyOrder, pointType) {}
387
390 virtual bool requireOrientation() const override
391 {
392 return true;
393 }
394
399 virtual
400 const char*
401 getName() const override {
402 return name_.c_str();
403 }
404
410 getHostBasis() const override {
412
413 auto hostBasis = Teuchos::rcp(new HostBasis(order_xy_, order_z_, pointType_));
414
415 return hostBasis;
416 }
417 };
418} // end namespace Intrepid2
419
420#endif /* Intrepid2_DerivedBasis_HDIV_WEDGE_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 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_Derived_HDIV_Family1_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
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 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_HDIV_Family2_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual bool requireOrientation() const override
True if orientation is required.
virtual const char * getName() const override
Returns basis name.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
Basis_Derived_HDIV_WEDGE(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Basis_Derived_HDIV_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
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...