Intrepid2
Intrepid2_DerivedBasis_HCURL_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 Mauro Perego (mperego@sandia.gov) or
38// Nate Roberts (nvrober@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
70#ifndef Intrepid2_DerivedBasis_HCURL_WEDGE_h
71#define Intrepid2_DerivedBasis_HCURL_WEDGE_h
72
73#include <Kokkos_DynRankView.hpp>
74
76#include "Intrepid2_Sacado.hpp"
77
80
81namespace Intrepid2
82{
83 template<class HCURL_TRI, class HGRAD_LINE>
85 : public Basis_TensorBasis<typename HGRAD_LINE::BasisBase>
86 {
87 public:
88 using OutputViewType = typename HGRAD_LINE::OutputViewType;
89 using PointViewType = typename HGRAD_LINE::PointViewType ;
90 using ScalarViewType = typename HGRAD_LINE::ScalarViewType;
91
92 using BasisBase = typename HGRAD_LINE::BasisBase;
93
94 using DeviceType = typename BasisBase::DeviceType;
95 using ExecutionSpace = typename BasisBase::ExecutionSpace;
96 using OutputValueType = typename BasisBase::OutputValueType;
97 using PointValueType = typename BasisBase::PointValueType;
98
99 using TriCurlBasis = HCURL_TRI;
100 using LineGradBasis = HGRAD_LINE;
101
103 public:
109 Basis_Derived_HCURL_Family1_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType = POINTTYPE_DEFAULT)
110 :
111 TensorBasis(Teuchos::rcp( new TriCurlBasis(polyOrder_xy,pointType)),
112 Teuchos::rcp( new LineGradBasis(polyOrder_z,pointType)))
113 {
114 this->functionSpace_ = FUNCTION_SPACE_HCURL;
115 this->setShardsTopologyAndTags();
116 }
117
119
122 virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const override
123 {
124 const EOperator & VALUE = OPERATOR_VALUE;
125 const EOperator & GRAD = OPERATOR_GRAD;
126 const EOperator & CURL = OPERATOR_CURL;
127 if (operatorType == VALUE)
128 {
129 // family 1 goes in x,y components
130 std::vector< std::vector<EOperator> > ops(2);
131 ops[0] = std::vector<EOperator>{VALUE,VALUE}; // occupies x,y
132 ops[1] = std::vector<EOperator>{}; // zero z
133 std::vector<double> weights {1.0, 0.0};
134 return OperatorTensorDecomposition(ops, weights);
135 }
136 else if (operatorType == CURL)
137 {
138 // curl of (f_x(x,y) g(z), f_y(x,y) g(z), 0), where f is in H(curl,tri), g in H(grad,line)
139 // x,y components of curl: rot(f) dg/dz, where rot(f) is a 90-degree rotation of f.
140 // z component of curl: curl(f) g, where curl(f) is the 2D curl, a scalar.
141 std::vector< std::vector<EOperator> > ops(2);
142 ops[0] = std::vector<EOperator>{VALUE,GRAD}; // occupies the x,y components
143 ops[1] = std::vector<EOperator>{CURL,VALUE}; // z component
144 std::vector<double> weights {1.0, 1.0};
145 OperatorTensorDecomposition opDecomposition(ops, weights);
146 opDecomposition.setRotateXYNinetyDegrees(true);
147 return opDecomposition;
148 }
149 else
150 {
151 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
152 }
153 }
154
162 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
163 const PointViewType inputPoints1, const PointViewType inputPoints2,
164 bool tensorPoints) const override
165 {
166 EOperator op1, op2;
167 if (operatorType == OPERATOR_VALUE)
168 {
169 op1 = OPERATOR_VALUE;
170 op2 = OPERATOR_VALUE;
171
172 // family 1 values go in the x,y components; 0 in the z component
173 auto outputValuesComponent12 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::pair<int,int>{0,2});
174 auto outputValuesComponent3 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),2);
175
176 // place 0 in the z component
177 Kokkos::deep_copy(outputValuesComponent3, 0.0);
178 this->TensorBasis::getValues(outputValuesComponent12,
179 inputPoints1, op1,
180 inputPoints2, op2, tensorPoints);
181
182 }
183 else if (operatorType == OPERATOR_CURL)
184 {
185 // curl of (f_x(x,y) g(z), f_y(x,y) g(z), 0), where f is in H(curl,tri), g in H(grad,line)
186 // x,y components of curl: rot(f) dg/dz, where rot(f) is a 90-degree rotation of f.
187 // z component of curl: curl(f) g, where curl(f) is the 2D curl, a scalar.
188
189 auto outputValuesComponent12 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::pair<int,int>{0,2});
190 auto outputValuesComponent3 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),2);
191
192 op1 = OPERATOR_VALUE;
193 op2 = OPERATOR_GRAD;
194
195 this->TensorBasis::getValues(outputValuesComponent12,
196 inputPoints1, op1,
197 inputPoints2, op2, tensorPoints);
198
199 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{outputValuesComponent12.extent_int(0),outputValuesComponent12.extent_int(1)});
200 Kokkos::parallel_for("wedge family 1 curl: rotateXYNinetyDegrees CW", policy,
201 KOKKOS_LAMBDA (const int &fieldOrdinal, const int &pointOrdinal) {
202 const auto f_x = outputValuesComponent12(fieldOrdinal,pointOrdinal,0); // copy
203 const auto &f_y = outputValuesComponent12(fieldOrdinal,pointOrdinal,1); // reference
204 outputValuesComponent12(fieldOrdinal,pointOrdinal,0) = -f_y;
205 outputValuesComponent12(fieldOrdinal,pointOrdinal,1) = f_x;
206 });
207
208 op1 = OPERATOR_CURL;
209 op2 = OPERATOR_VALUE;
210 this->TensorBasis::getValues(outputValuesComponent3,
211 inputPoints1, op1,
212 inputPoints2, op2, tensorPoints);
213 }
214 else
215 {
216 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
217 }
218 }
219
231 virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
232 auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
233 auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
234 this->TensorBasis::getDofCoeffs(dofCoeffs1);
235 Kokkos::deep_copy(dofCoeffs2,0.0);
236 }
237 };
238
239 template<class HGRAD_TRI, class HVOL_LINE>
241 : public Basis_TensorBasis<typename HVOL_LINE::BasisBase>
242 {
243
244 public:
245 using OutputViewType = typename HVOL_LINE::OutputViewType;
246 using PointViewType = typename HVOL_LINE::PointViewType ;
247 using ScalarViewType = typename HVOL_LINE::ScalarViewType;
248
249 using TriGradBasis = HGRAD_TRI;
250 using LineHVolBasis = HVOL_LINE;
251
252 using BasisBase = typename HVOL_LINE::BasisBase;
254
255 using DeviceType = typename BasisBase::DeviceType;
256 using ExecutionSpace = typename BasisBase::ExecutionSpace;
257 using OutputValueType = typename BasisBase::OutputValueType;
258 using PointValueType = typename BasisBase::PointValueType;
259
265 Basis_Derived_HCURL_Family2_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType = POINTTYPE_DEFAULT)
266 :
267 TensorBasis(Teuchos::rcp( new TriGradBasis(polyOrder_xy, pointType) ),
268 Teuchos::rcp( new LineHVolBasis(polyOrder_z-1, pointType) ))
269 {
270 this->functionSpace_ = FUNCTION_SPACE_HCURL;
271 this->setShardsTopologyAndTags();
272 }
273
276 OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const override
277 {
278 if (operatorType == OPERATOR_CURL)
279 {
280 // curl of (0,0,f) is (df/dy, -df/dx, 0)
281
282 // this is a rotation of gradient of the triangle H^1 basis, times the H(vol) line value
283 std::vector< std::vector<EOperator> > ops(2);
284 ops[0] = std::vector<EOperator>{OPERATOR_GRAD,OPERATOR_VALUE}; // occupies the x,y components
285 ops[1] = std::vector<EOperator>{};
286 std::vector<double> weights {-1.0, 0.0}; // -1 because the rotation goes from (df/dx,df/dy) --> (-df/dy,df/dx), and we want (df/dy,-df/dx)
287 OperatorTensorDecomposition opDecomposition(ops, weights);
288 opDecomposition.setRotateXYNinetyDegrees(true);
289 return opDecomposition;
290 }
291 else if (OPERATOR_VALUE == operatorType)
292 {
293 // family 2 goes in z component
294 std::vector< std::vector<EOperator> > ops(2);
295 ops[0] = std::vector<EOperator>{}; // because family 1 identifies this as spanning (x,y), empty op here will also span (x,y)
296 ops[1] = std::vector<EOperator>{OPERATOR_VALUE,OPERATOR_VALUE}; // z component
297 std::vector<double> weights {0.0, 1.0};
298 return OperatorTensorDecomposition(ops, weights);
299 }
300 else
301 {
302 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type");
303 }
304 }
305
307
315 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
316 const PointViewType inputPoints1, const PointViewType inputPoints2,
317 bool tensorPoints) const override
318 {
319 EOperator op1, op2;
320 if (operatorType == OPERATOR_VALUE)
321 {
322 op1 = OPERATOR_VALUE;
323 op2 = OPERATOR_VALUE;
324
325 // family 2 values go in z component, 0 in (x,y)
326 auto outputValuesComponent12 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::pair<int,int>{0,2});
327 auto outputValuesComponent3 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),2);
328
329 // place 0 in the x,y components
330 Kokkos::deep_copy(outputValuesComponent12, 0.0);
331 this->TensorBasis::getValues(outputValuesComponent3,
332 inputPoints1, op1,
333 inputPoints2, op2, tensorPoints);
334
335 }
336 else if (operatorType == OPERATOR_CURL)
337 {
338 // curl of (0,0,f) is (df/dy, -df/dx, 0)
339
340 auto outputValuesComponent12 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::pair<int,int>{0,2});
341 auto outputValuesComponent3 = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),2);
342
343 op1 = OPERATOR_GRAD;
344 op2 = OPERATOR_VALUE;
345
346 this->TensorBasis::getValues(outputValuesComponent12,
347 inputPoints1, op1,
348 inputPoints2, op2, tensorPoints);
349
350 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{outputValuesComponent12.extent_int(0),outputValuesComponent12.extent_int(1)});
351 Kokkos::parallel_for("wedge family 2 curl: rotateXYNinetyDegrees CCW", policy,
352 KOKKOS_LAMBDA (const int &fieldOrdinal, const int &pointOrdinal) {
353 const auto f_x = outputValuesComponent12(fieldOrdinal,pointOrdinal,0); // copy
354 const auto &f_y = outputValuesComponent12(fieldOrdinal,pointOrdinal,1); // reference
355 outputValuesComponent12(fieldOrdinal,pointOrdinal,0) = f_y;
356 outputValuesComponent12(fieldOrdinal,pointOrdinal,1) = -f_x;
357 });
358
359 Kokkos::deep_copy(outputValuesComponent3, 0.0);
360 }
361 else
362 {
363 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported");
364 }
365 }
366
378 virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
379 auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
380 auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
381 Kokkos::deep_copy(dofCoeffs1,0.0);
382 this->TensorBasis::getDofCoeffs(dofCoeffs2);
383 }
384 };
385
386 template<class HGRAD_TRI, class HCURL_TRI, class HGRAD_LINE, class HVOL_LINE>
388 : public Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>
389 {
392 using DirectSumBasis = Basis_DirectSumBasis <typename HGRAD_LINE::BasisBase>;
393 public:
394 using BasisBase = typename HGRAD_LINE::BasisBase;
395
396 protected:
397 std::string name_;
398 ordinal_type order_xy_;
399 ordinal_type order_z_;
400 EPointType pointType_;
401
402 public:
403 using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace;
404 using OutputValueType = typename HGRAD_LINE::OutputValueType;
405 using PointValueType = typename HGRAD_LINE::PointValueType;
406
412 Basis_Derived_HCURL_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
413 :
414 DirectSumBasis(Teuchos::rcp( new Family1(polyOrder_xy, polyOrder_z, pointType) ),
415 Teuchos::rcp( new Family2(polyOrder_xy, polyOrder_z, pointType) ))
416 {
417 this->functionSpace_ = FUNCTION_SPACE_HCURL;
418
419 std::ostringstream basisName;
420 basisName << "HCURL_WEDGE (" << this->DirectSumBasis::getName() << ")";
421 name_ = basisName.str();
422
423 order_xy_ = polyOrder_xy;
424 order_z_ = polyOrder_z;
425 pointType_ = pointType;
426 }
427
432 Basis_Derived_HCURL_WEDGE(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) : Basis_Derived_HCURL_WEDGE(polyOrder, polyOrder, pointType) {}
433
436 virtual bool requireOrientation() const override
437 {
438 return true;
439 }
440
445 virtual
446 const char*
447 getName() const override {
448 return name_.c_str();
449 }
450
456 getHostBasis() const override {
458
459 auto hostBasis = Teuchos::rcp(new HostBasis(order_xy_, order_z_, pointType_));
460
461 return hostBasis;
462 }
463 };
464} // end namespace Intrepid2
465
466#endif /* Intrepid2_DerivedBasis_HCURL_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.
Basis_Derived_HCURL_Family1_WEDGE(int polyOrder_xy, int polyOrder_z, 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 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.
Basis_Derived_HCURL_Family2_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual void getDofCoeffs(ScalarViewType dofCoeffs) const override
Fills in coefficients of degrees of freedom for Lagrangian basis on the reference cell.
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(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.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
Basis_Derived_HCURL_WEDGE(int polyOrder_xy, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual const char * getName() const override
Returns basis name.
Basis_Derived_HCURL_WEDGE(int polyOrder, 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...