Intrepid2
Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.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
49#ifndef Intrepid2_IntegratedLegendreBasis_HGRAD_LINE_h
50#define Intrepid2_IntegratedLegendreBasis_HGRAD_LINE_h
51
52#include <Kokkos_DynRankView.hpp>
53
54#include <Intrepid2_config.h>
55
56#include "Intrepid2_Basis.hpp"
58#include "Intrepid2_Utils.hpp"
59
60namespace Intrepid2
61{
67 template<class DeviceType, class OutputScalar, class PointScalar,
68 class OutputFieldType, class InputPointsType>
70 {
71 using ExecutionSpace = typename DeviceType::execution_space;
72 using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
73 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
74 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
75
76 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
77 using TeamMember = typename TeamPolicy::member_type;
78
79 EOperator opType_; // OPERATOR_VALUE or OPERATOR_GRAD
80
81 OutputFieldType output_; // F,P
82 InputPointsType inputPoints_; // P,D
83
84 int polyOrder_;
85 bool defineVertexFunctions_;
86 int numFields_, numPoints_;
87
88 size_t fad_size_output_;
89
90 Hierarchical_HGRAD_LINE_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints,
91 int polyOrder, bool defineVertexFunctions)
92 : opType_(opType), output_(output), inputPoints_(inputPoints),
93 polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
94 fad_size_output_(getScalarDimensionForView(output))
95 {
96 numFields_ = output.extent_int(0);
97 numPoints_ = output.extent_int(1);
98 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
99 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != polyOrder_+1, std::invalid_argument, "output field size does not match basis cardinality");
100 }
101
102 KOKKOS_INLINE_FUNCTION
103 void operator()( const TeamMember & teamMember ) const
104 {
105 auto pointOrdinal = teamMember.league_rank();
106 OutputScratchView field_values_at_point;
107 if (fad_size_output_ > 0) {
108 field_values_at_point = OutputScratchView(teamMember.team_shmem(), numFields_, fad_size_output_);
109 }
110 else {
111 field_values_at_point = OutputScratchView(teamMember.team_shmem(), numFields_);
112 }
113
114 const auto & input_x = inputPoints_(pointOrdinal,0);
115 const bool taking_derivative = (opType_ != OPERATOR_VALUE);
116 const bool callingShiftedScaledLegendre = (opType_ == OPERATOR_VALUE) || (opType_ == OPERATOR_GRAD) || (opType_ == OPERATOR_D1);
117
118 // shiftedScaledIntegratedLegendreValues{_dx} expects x in [0,1]
119 const PointScalar x = callingShiftedScaledLegendre ? PointScalar((input_x + 1.0)/2.0) : PointScalar(input_x);
120 const double legendreScaling = 1.0;
121 const double outputScaling = taking_derivative ? 0.5 : 1.0; // output scaling -- 0.5 if we take derivatives, 1.0 otherwise
122
123 switch (opType_)
124 {
125 case OPERATOR_VALUE:
126 // field values are integrated Legendre polynomials, except for the first and second field,
127 // which may be 1 and x or x and 1-x, depending on whether the vertex compatibility flag is set.
128 Polynomials::shiftedScaledIntegratedLegendreValues(field_values_at_point, polyOrder_, x, legendreScaling);
129
130 // note that because shiftedScaledIntegratedLegendreValues determines field values recursively, there is not much
131 // opportunity at that level for further parallelism
132
133 if (defineVertexFunctions_)
134 {
135 field_values_at_point(0) = 1. - x;
136 field_values_at_point(1) = x;
137 }
138 break;
139 case OPERATOR_GRAD:
140 case OPERATOR_D1:
141 // field values are Legendre polynomials, except for the first and second field,
142 // which may be 0 and 1 or -1 and 1, depending on whether the vertex compatibility flag is set.
143 Polynomials::shiftedScaledIntegratedLegendreValues_dx(field_values_at_point, polyOrder_, x, legendreScaling);
144
145 // note that because shiftedScaledIntegratedLegendreValues_dx determines field values recursively, there is not much
146 // opportunity at that level for further parallelism
147
148 if (defineVertexFunctions_)
149 {
150 field_values_at_point(0) = -1.0; // derivative of 1-x
151 field_values_at_point(1) = 1.0; // derivative of x
152 }
153 break;
154 case OPERATOR_D2:
155 case OPERATOR_D3:
156 case OPERATOR_D4:
157 case OPERATOR_D5:
158 case OPERATOR_D6:
159 case OPERATOR_D7:
160 case OPERATOR_D8:
161 case OPERATOR_D9:
162 case OPERATOR_D10:
163 {
164 auto derivativeOrder = getOperatorOrder(opType_) - 1;
165 Polynomials::legendreDerivativeValues(field_values_at_point, polyOrder_, x, derivativeOrder);
166
167 // L_i is defined in terms of an integral of P_(i-1), so we need to shift the values by 1
168 if (numFields_ >= 3)
169 {
170 OutputScalar Pn_minus_one = field_values_at_point(1);
171 for (int fieldOrdinal=2; fieldOrdinal<numFields_; fieldOrdinal++)
172 {
173 OutputScalar Pn = field_values_at_point(fieldOrdinal);
174 field_values_at_point(fieldOrdinal) = Pn_minus_one;
175 Pn_minus_one = Pn;
176 }
177 }
178 if (numFields_ >= 1) field_values_at_point(0) = 0.0;
179 if (numFields_ >= 2) field_values_at_point(1) = 0.0;
180 // legendreDerivativeValues works on [-1,1], so no per-derivative scaling is necessary
181 // however, there is a factor of 0.5 that comes from the scaling of the Legendre polynomials prior to integration
182 // in the shiftedScaledIntegratedLegendreValues -- the first derivative of our integrated polynomials is 0.5 times the Legendre polynomial
183 break;
184 }
185 default:
186 // unsupported operator type
187 device_assert(false);
188 }
189
190 // copy the values into the output container
191 for (int fieldOrdinal=0; fieldOrdinal<numFields_; fieldOrdinal++)
192 {
193 // access() allows us to write one line that applies both to gradient (for which outputValues has rank 3, but third rank has only one entry) and to value (rank 2)
194 output_.access(fieldOrdinal,pointOrdinal,0) = outputScaling * field_values_at_point(fieldOrdinal);
195 }
196 }
197
198 // Provide the shared memory capacity.
199 // This function takes the team_size as an argument,
200 // which allows team_size-dependent allocations.
201 size_t team_shmem_size (int team_size) const
202 {
203 // we want to use shared memory to create a fast buffer that we can use for basis computations
204 size_t shmem_size = 0;
205 if (fad_size_output_ > 0)
206 shmem_size += OutputScratchView::shmem_size(numFields_, fad_size_output_);
207 else
208 shmem_size += OutputScratchView::shmem_size(numFields_);
209
210 return shmem_size;
211 }
212 };
213
231 template<typename DeviceType,
232 typename OutputScalar = double,
233 typename PointScalar = double,
234 bool defineVertexFunctions = true, // if defineVertexFunctions is true, first and second basis functions are x and 1-x. Otherwise, they are 1 and x.
235 bool useMinusOneToOneReferenceElement = true> // if useMinusOneToOneReferenceElement is true, basis is define on [-1,1]. Otherwise, [0,1].
237 : public Basis<DeviceType,OutputScalar,PointScalar>
238 {
239 public:
242
245
246 using typename BasisBase::OutputViewType;
247 using typename BasisBase::PointViewType;
248 using typename BasisBase::ScalarViewType;
249
250 using typename BasisBase::ExecutionSpace;
251
252 protected:
253 int polyOrder_; // the maximum order of the polynomial
254 bool defineVertexFunctions_; // if true, first and second basis functions are x and 1-x. Otherwise, they are 1 and x.
255 EPointType pointType_;
256 public:
267 IntegratedLegendreBasis_HGRAD_LINE(int polyOrder, EPointType pointType=POINTTYPE_DEFAULT)
268 :
269 polyOrder_(polyOrder),
270 pointType_(pointType)
271 {
272 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
273
274 this->basisCardinality_ = polyOrder+1;
275 this->basisDegree_ = polyOrder;
276 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
277 this->basisType_ = BASIS_FEM_HIERARCHICAL;
278 this->basisCoordinates_ = COORDINATES_CARTESIAN;
279 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
280
281 const int degreeLength = 1;
282 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) line polynomial degree lookup", this->basisCardinality_, degreeLength);
283 this->fieldOrdinalH1PolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) line polynomial H^1 degree lookup", this->basisCardinality_, degreeLength);
284
285 for (int i=0; i<this->basisCardinality_; i++)
286 {
287 // for H(grad) line, if defineVertexFunctions is false, first basis member is constant, second is first-degree, etc.
288 // if defineVertexFunctions is true, then the only difference is that the entry is also degree 1
289 this->fieldOrdinalPolynomialDegree_ (i,0) = i;
290 this->fieldOrdinalH1PolynomialDegree_(i,0) = i;
291 }
292 if (defineVertexFunctions)
293 {
294 this->fieldOrdinalPolynomialDegree_ (0,0) = 1;
295 this->fieldOrdinalH1PolynomialDegree_(0,0) = 1;
296 }
297
298 // initialize tags
299 {
300 const auto & cardinality = this->basisCardinality_;
301
302 // Basis-dependent initializations
303 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
304 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
305 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
306 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
307
308 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
309
310 if (defineVertexFunctions) {
311 {
312 const ordinal_type v0 = 0;
313 tagView(v0*tagSize+0) = 0; // vertex dof
314 tagView(v0*tagSize+1) = 0; // vertex id
315 tagView(v0*tagSize+2) = 0; // local dof id
316 tagView(v0*tagSize+3) = 1; // total number of dofs in this vertex
317
318 const ordinal_type v1 = 1;
319 tagView(v1*tagSize+0) = 0; // vertex dof
320 tagView(v1*tagSize+1) = 1; // vertex id
321 tagView(v1*tagSize+2) = 0; // local dof id
322 tagView(v1*tagSize+3) = 1; // total number of dofs in this vertex
323
324 const ordinal_type iend = cardinality - 2;
325 for (ordinal_type i=0;i<iend;++i) {
326 const auto e = i + 2;
327 tagView(e*tagSize+0) = 1; // edge dof
328 tagView(e*tagSize+1) = 0; // edge id
329 tagView(e*tagSize+2) = i; // local dof id
330 tagView(e*tagSize+3) = iend; // total number of dofs in this edge
331 }
332 }
333 } else {
334 for (ordinal_type i=0;i<cardinality;++i) {
335 tagView(i*tagSize+0) = 1; // edge dof
336 tagView(i*tagSize+1) = 0; // edge id
337 tagView(i*tagSize+2) = i; // local dof id
338 tagView(i*tagSize+3) = cardinality; // total number of dofs in this edge
339 }
340 }
341
342 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
343 // tags are constructed on host
345 this->ordinalToTag_,
346 tagView,
347 this->basisCardinality_,
348 tagSize,
349 posScDim,
350 posScOrd,
351 posDfOrd);
352 }
353 }
354
359 const char* getName() const override {
360 return "Intrepid2_IntegratedLegendreBasis_HGRAD_LINE";
361 }
362
365 virtual bool requireOrientation() const override {
366 return false;
367 }
368
369 // since the getValues() below only overrides the FEM variant, we specify that
370 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
371 // (It's an error to use the FVD variant on this basis.)
373
392 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
393 const EOperator operatorType = OPERATOR_VALUE ) const override
394 {
395 auto numPoints = inputPoints.extent_int(0);
396
398
399 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
400
401 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
402 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
403 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
404 const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
405
406 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
407 Kokkos::parallel_for("Hierarchical_HGRAD_LINE_Functor", policy, functor);
408 }
409
415 getHostBasis() const override {
416 using HostDeviceType = typename Kokkos::HostSpace::device_type;
418 return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
419 }
420 };
421} // end namespace Intrepid2
422
423#endif /* Intrepid2_IntegratedLegendreBasis_HGRAD_LINE_h */
Header file for the abstract base class Intrepid2::Basis.
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
KOKKOS_INLINE_FUNCTION void device_assert(bool val)
Free functions, callable from device code, that implement various polynomials useful in basis definit...
Header function for Intrepid2::Util class and other utility functions.
KOKKOS_INLINE_FUNCTION constexpr unsigned getScalarDimensionForView(const ViewType &view)
Returns the size of the Scalar dimension for the View. This is 0 for non-AD types....
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
EBasis basisType_
Type of the basis.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
void setOrdinalTagData(OrdinalTypeView3D &tagToOrdinal, OrdinalTypeView2D &ordinalToTag, const OrdinalTypeView1D tags, const ordinal_type basisCard, const ordinal_type tagSize, const ordinal_type posScDim, const ordinal_type posScOrd, const ordinal_type posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
OrdinalTypeArray2DHost fieldOrdinalH1PolynomialDegree_
H^1 polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
EFunctionSpace functionSpace_
The function space in which the basis is defined.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
virtual BasisPtr< typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
IntegratedLegendreBasis_HGRAD_LINE(int polyOrder, EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
virtual bool requireOrientation() const override
True if orientation is required.
Functor for computing values for the IntegratedLegendreBasis_HGRAD_LINE class.