Intrepid2
Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.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_TRI_h
50#define Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h
51
52#include <Kokkos_DynRankView.hpp>
53
54#include <Intrepid2_config.h>
55
56#include "Intrepid2_Basis.hpp"
59#include "Intrepid2_Utils.hpp"
60
61namespace Intrepid2
62{
68 template<class DeviceType, class OutputScalar, class PointScalar,
69 class OutputFieldType, class InputPointsType>
71 {
72 using ExecutionSpace = typename DeviceType::execution_space;
73 using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
74 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
75 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
76
77 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
78 using TeamMember = typename TeamPolicy::member_type;
79
80 EOperator opType_;
81
82 OutputFieldType output_; // F,P
83 InputPointsType inputPoints_; // P,D
84
85 int polyOrder_;
86 bool defineVertexFunctions_;
87 int numFields_, numPoints_;
88
89 size_t fad_size_output_;
90
91 static const int numVertices = 3;
92 static const int numEdges = 3;
93 const int edge_start_[numEdges] = {0,1,0}; // edge i is from edge_start_[i] to edge_end_[i]
94 const int edge_end_[numEdges] = {1,2,2}; // edge i is from edge_start_[i] to edge_end_[i]
95
96 Hierarchical_HGRAD_TRI_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints,
97 int polyOrder, bool defineVertexFunctions)
98 : opType_(opType), output_(output), inputPoints_(inputPoints),
99 polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
100 fad_size_output_(getScalarDimensionForView(output))
101 {
102 numFields_ = output.extent_int(0);
103 numPoints_ = output.extent_int(1);
104 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
105 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)/2, std::invalid_argument, "output field size does not match basis cardinality");
106 }
107
108 KOKKOS_INLINE_FUNCTION
109 void operator()( const TeamMember & teamMember ) const
110 {
111 auto pointOrdinal = teamMember.league_rank();
112 OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
113 if (fad_size_output_ > 0) {
114 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
115 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
116 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
117 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
118 }
119 else {
120 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
121 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
122 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
123 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
124 }
125
126 const auto & x = inputPoints_(pointOrdinal,0);
127 const auto & y = inputPoints_(pointOrdinal,1);
128
129 // write as barycentric coordinates:
130 const PointScalar lambda[3] = {1. - x - y, x, y};
131 const PointScalar lambda_dx[3] = {-1., 1., 0.};
132 const PointScalar lambda_dy[3] = {-1., 0., 1.};
133
134 const int num1DEdgeFunctions = polyOrder_ - 1;
135
136 switch (opType_)
137 {
138 case OPERATOR_VALUE:
139 {
140 // vertex functions come first, according to vertex ordering: (0,0), (1,0), (0,1)
141 for (int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
142 {
143 output_(vertexOrdinal,pointOrdinal) = lambda[vertexOrdinal];
144 }
145 if (!defineVertexFunctions_)
146 {
147 // "DG" basis case
148 // here, we overwrite the first vertex function with 1:
149 output_(0,pointOrdinal) = 1.0;
150 }
151
152 // edge functions
153 int fieldOrdinalOffset = 3;
154 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
155 {
156 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
157 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
158
159 Polynomials::shiftedScaledIntegratedLegendreValues(edge_field_values_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
160 for (int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
161 {
162 // the first two integrated legendre functions are essentially the vertex functions; hence the +2 on on the RHS here:
163 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal) = edge_field_values_at_point(edgeFunctionOrdinal+2);
164 }
165 fieldOrdinalOffset += num1DEdgeFunctions;
166 }
167
168 // face functions
169 {
170 // these functions multiply the edge functions from the 01 edge by integrated Jacobi functions, appropriately scaled
171 const double jacobiScaling = 1.0; // s0 + s1 + s2
172
173 const int max_ij_sum = polyOrder_;
174 const int min_i = 2;
175 const int min_j = 1;
176 const int min_ij_sum = min_i + min_j;
177 for (int ij_sum = min_ij_sum; ij_sum <= max_ij_sum; ij_sum++)
178 {
179 for (int i=min_i; i<=ij_sum-min_j; i++)
180 {
181 const int j = ij_sum - i;
182 const int edgeBasisOrdinal = i+numVertices-2; // i+1: where the value of the edge function is stored in output_
183 const auto & edgeValue = output_(edgeBasisOrdinal,pointOrdinal);
184 const double alpha = i*2.0;
185
186 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_values_at_point, alpha, polyOrder_-2, lambda[2], jacobiScaling);
187 const auto & jacobiValue = jacobi_values_at_point(j);
188 output_(fieldOrdinalOffset,pointOrdinal) = edgeValue * jacobiValue;
189 fieldOrdinalOffset++;
190 }
191 }
192 }
193 } // end OPERATOR_VALUE
194 break;
195 case OPERATOR_GRAD:
196 case OPERATOR_D1:
197 {
198 // vertex functions
199 if (defineVertexFunctions_)
200 {
201 // standard, "CG" basis case
202 // first vertex function is 1-x-y
203 output_(0,pointOrdinal,0) = -1.0;
204 output_(0,pointOrdinal,1) = -1.0;
205 }
206 else
207 {
208 // "DG" basis case
209 // here, the first "vertex" function is 1, so the derivative is 0:
210 output_(0,pointOrdinal,0) = 0.0;
211 output_(0,pointOrdinal,1) = 0.0;
212 }
213 // second vertex function is x
214 output_(1,pointOrdinal,0) = 1.0;
215 output_(1,pointOrdinal,1) = 0.0;
216 // third vertex function is y
217 output_(2,pointOrdinal,0) = 0.0;
218 output_(2,pointOrdinal,1) = 1.0;
219
220 // edge functions
221 int fieldOrdinalOffset = 3;
222 /*
223 Per Fuentes et al. (see Appendix E.1, E.2), the edge functions, defined for i ≥ 2, are
224 [L_i](s0,s1) = L_i(s1; s0+s1)
225 and have gradients:
226 grad [L_i](s0,s1) = [P_{i-1}](s0,s1) grad s1 + [R_{i-1}](s0,s1) grad (s0 + s1)
227 where
228 [R_{i-1}](s0,s1) = R_{i-1}(s1; s0+s1) = d/dt L_{i}(s0; s0+s1)
229 The P_i we have implemented in shiftedScaledLegendreValues, while d/dt L_{i+1} is
230 implemented in shiftedScaledIntegratedLegendreValues_dt.
231 */
232 // rename the scratch memory to match our usage here:
233 auto & P_i_minus_1 = edge_field_values_at_point;
234 auto & L_i_dt = jacobi_values_at_point;
235 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
236 {
237 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
238 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
239
240 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
241 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
242 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
243 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
244
245 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
246 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
247 for (int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
248 {
249 // the first two (integrated) Legendre functions are essentially the vertex functions; hence the +2 here:
250 const int i = edgeFunctionOrdinal+2;
251 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = P_i_minus_1(i-1) * s1_dx + L_i_dt(i) * (s1_dx + s0_dx);
252 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = P_i_minus_1(i-1) * s1_dy + L_i_dt(i) * (s1_dy + s0_dy);
253 }
254 fieldOrdinalOffset += num1DEdgeFunctions;
255 }
256
257 /*
258 Fuentes et al give the face functions as phi_{ij}, with gradient:
259 grad phi_{ij}(s0,s1,s2) = [L^{2i}_j](s0+s1,s2) grad [L_i](s0,s1) + [L_i](s0,s1) grad [L^{2i}_j](s0+s1,s2)
260 where:
261 - grad [L_i](s0,s1) is the edge function gradient we computed above
262 - [L_i](s0,s1) is the edge function which we have implemented above (in OPERATOR_VALUE)
263 - L^{2i}_j is a Jacobi polynomial with:
264 [L^{2i}_j](s0,s1) = L^{2i}_j(s1;s0+s1)
265 and the gradient for j ≥ 1 is
266 grad [L^{2i}_j](s0,s1) = [P^{2i}_{j-1}](s0,s1) grad s1 + [R^{2i}_{j-1}(s0,s1)] grad (s0 + s1)
267 Here,
268 [P^{2i}_{j-1}](s0,s1) = P^{2i}_{j-1}(s1,s0+s1)
269 and
270 [R^{2i}_{j-1}(s0,s1)] = d/dt L^{2i}_j(s1,s0+s1)
271 We have implemented P^{alpha}_{j} as shiftedScaledJacobiValues,
272 and d/dt L^{alpha}_{j} as shiftedScaledIntegratedJacobiValues_dt.
273 */
274 // rename the scratch memory to match our usage here:
275 auto & P_2i_j_minus_1 = edge_field_values_at_point;
276 auto & L_2i_j_dt = jacobi_values_at_point;
277 auto & L_i = other_values_at_point;
278 auto & L_2i_j = other_values2_at_point;
279 {
280 // face functions multiply the edge functions from the 01 edge by integrated Jacobi functions, appropriately scaled
281 const double jacobiScaling = 1.0; // s0 + s1 + s2
282
283 const int max_ij_sum = polyOrder_;
284 const int min_i = 2;
285 const int min_j = 1;
286 const int min_ij_sum = min_i + min_j;
287 for (int ij_sum = min_ij_sum; ij_sum <= max_ij_sum; ij_sum++)
288 {
289 for (int i=min_i; i<=ij_sum-min_j; i++)
290 {
291 const int j = ij_sum - i;
292 // the edge function here is for edge 01, in the first set of edge functions.
293 const int edgeBasisOrdinal = i+numVertices-2; // i+1: where the value of the edge function is stored in output_
294 const auto & grad_L_i_dx = output_(edgeBasisOrdinal,pointOrdinal,0);
295 const auto & grad_L_i_dy = output_(edgeBasisOrdinal,pointOrdinal,1);
296
297 const double alpha = i*2.0;
298
299 Polynomials::shiftedScaledIntegratedLegendreValues (L_i, polyOrder_, lambda[1], lambda[0]+lambda[1]);
300 Polynomials::shiftedScaledIntegratedJacobiValues_dt(L_2i_j_dt, alpha, polyOrder_, lambda[2], jacobiScaling);
301 Polynomials::shiftedScaledIntegratedJacobiValues ( L_2i_j, alpha, polyOrder_, lambda[2], jacobiScaling);
302 Polynomials::shiftedScaledJacobiValues(P_2i_j_minus_1, alpha, polyOrder_-1, lambda[2], jacobiScaling);
303
304 const auto & s0_dx = lambda_dx[0];
305 const auto & s0_dy = lambda_dy[0];
306 const auto & s1_dx = lambda_dx[1];
307 const auto & s1_dy = lambda_dy[1];
308 const auto & s2_dx = lambda_dx[2];
309 const auto & s2_dy = lambda_dy[2];
310
311 const OutputScalar basisValue_dx = L_2i_j(j) * grad_L_i_dx + L_i(i) * (P_2i_j_minus_1(j-1) * s2_dx + L_2i_j_dt(j) * (s0_dx + s1_dx + s2_dx));
312 const OutputScalar basisValue_dy = L_2i_j(j) * grad_L_i_dy + L_i(i) * (P_2i_j_minus_1(j-1) * s2_dy + L_2i_j_dt(j) * (s0_dy + s1_dy + s2_dy));
313
314 output_(fieldOrdinalOffset,pointOrdinal,0) = basisValue_dx;
315 output_(fieldOrdinalOffset,pointOrdinal,1) = basisValue_dy;
316 fieldOrdinalOffset++;
317 }
318 }
319 }
320 }
321 break;
322 case OPERATOR_D2:
323 case OPERATOR_D3:
324 case OPERATOR_D4:
325 case OPERATOR_D5:
326 case OPERATOR_D6:
327 case OPERATOR_D7:
328 case OPERATOR_D8:
329 case OPERATOR_D9:
330 case OPERATOR_D10:
331 INTREPID2_TEST_FOR_ABORT( true,
332 ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
333 default:
334 // unsupported operator type
335 device_assert(false);
336 }
337 }
338
339 // Provide the shared memory capacity.
340 // This function takes the team_size as an argument,
341 // which allows team_size-dependent allocations.
342 size_t team_shmem_size (int team_size) const
343 {
344 // we will use shared memory to create a fast buffer for basis computations
345 size_t shmem_size = 0;
346 if (fad_size_output_ > 0)
347 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
348 else
349 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);
350
351 return shmem_size;
352 }
353 };
354
371 template<typename DeviceType,
372 typename OutputScalar = double,
373 typename PointScalar = double,
374 bool defineVertexFunctions = true> // if defineVertexFunctions is true, first three basis functions are 1-x-y, x, and y. Otherwise, they are 1, x, and y.
376 : public Basis<DeviceType,OutputScalar,PointScalar>
377 {
378 public:
381
384
385 using typename BasisBase::OutputViewType;
386 using typename BasisBase::PointViewType;
387 using typename BasisBase::ScalarViewType;
388
389 using typename BasisBase::ExecutionSpace;
390
391 protected:
392 int polyOrder_; // the maximum order of the polynomial
393 EPointType pointType_;
394 public:
405 IntegratedLegendreBasis_HGRAD_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
406 :
407 polyOrder_(polyOrder),
408 pointType_(pointType)
409 {
410 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
411
412 this->basisCardinality_ = ((polyOrder+2) * (polyOrder+1)) / 2;
413 this->basisDegree_ = polyOrder;
414 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
415 this->basisType_ = BASIS_FEM_HIERARCHICAL;
416 this->basisCoordinates_ = COORDINATES_CARTESIAN;
417 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
418
419 const int degreeLength = 1;
420 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);
421 this->fieldOrdinalH1PolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);
422
423 int fieldOrdinalOffset = 0;
424 // **** vertex functions **** //
425 const int numVertices = this->basisCellTopology_.getVertexCount();
426 const int numFunctionsPerVertex = 1;
427 const int numVertexFunctions = numVertices * numFunctionsPerVertex;
428 for (int i=0; i<numVertexFunctions; i++)
429 {
430 // for H(grad) on triangle, if defineVertexFunctions is false, first three basis members are linear
431 // if not, then the only difference is that the first member is constant
432 this->fieldOrdinalPolynomialDegree_ (i,0) = 1;
433 this->fieldOrdinalH1PolynomialDegree_(i,0) = 1;
434 }
435 if (!defineVertexFunctions)
436 {
437 this->fieldOrdinalPolynomialDegree_ (0,0) = 0;
438 this->fieldOrdinalH1PolynomialDegree_(0,0) = 0;
439 }
440 fieldOrdinalOffset += numVertexFunctions;
441
442 // **** edge functions **** //
443 const int numFunctionsPerEdge = polyOrder - 1; // bubble functions: all but the vertices
444 const int numEdges = this->basisCellTopology_.getEdgeCount();
445 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
446 {
447 for (int i=0; i<numFunctionsPerEdge; i++)
448 {
449 this->fieldOrdinalPolynomialDegree_(i+fieldOrdinalOffset,0) = i+2; // vertex functions are 1st order; edge functions start at order 2
450 this->fieldOrdinalH1PolynomialDegree_(i+fieldOrdinalOffset,0) = i+2; // vertex functions are 1st order; edge functions start at order 2
451 }
452 fieldOrdinalOffset += numFunctionsPerEdge;
453 }
454
455 // **** face functions **** //
456 const int max_ij_sum = polyOrder;
457 const int min_i = 2;
458 const int min_j = 1;
459 const int min_ij_sum = min_i + min_j;
460 for (int ij_sum = min_ij_sum; ij_sum <= max_ij_sum; ij_sum++)
461 {
462 for (int i=min_i; i<=ij_sum-min_j; i++)
463 {
464 const int j = ij_sum - i;
465 this->fieldOrdinalPolynomialDegree_(fieldOrdinalOffset,0) = i+j;
466 this->fieldOrdinalH1PolynomialDegree_(fieldOrdinalOffset,0) = i+j;
467 fieldOrdinalOffset++;
468 }
469 }
470 const int numFaces = 1;
471 const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2;
472 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
473
474 // initialize tags
475 {
476 const auto & cardinality = this->basisCardinality_;
477
478 // Basis-dependent initializations
479 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
480 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
481 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
482 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
483
484 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
485 const int vertexDim = 0, edgeDim = 1, faceDim = 2;
486
487 if (defineVertexFunctions) {
488 {
489 int tagNumber = 0;
490 for (int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
491 {
492 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerVertex; functionOrdinal++)
493 {
494 tagView(tagNumber*tagSize+0) = vertexDim; // vertex dimension
495 tagView(tagNumber*tagSize+1) = vertexOrdinal; // vertex id
496 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
497 tagView(tagNumber*tagSize+3) = numFunctionsPerVertex; // total number of dofs in this vertex
498 tagNumber++;
499 }
500 }
501 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
502 {
503 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
504 {
505 tagView(tagNumber*tagSize+0) = edgeDim; // edge dimension
506 tagView(tagNumber*tagSize+1) = edgeOrdinal; // edge id
507 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
508 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge; // total number of dofs on this edge
509 tagNumber++;
510 }
511 }
512 for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
513 {
514 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
515 {
516 tagView(tagNumber*tagSize+0) = faceDim; // face dimension
517 tagView(tagNumber*tagSize+1) = faceOrdinal; // face id
518 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
519 tagView(tagNumber*tagSize+3) = numFunctionsPerFace; // total number of dofs on this face
520 tagNumber++;
521 }
522 }
523 }
524 } else {
525 for (ordinal_type i=0;i<cardinality;++i) {
526 tagView(i*tagSize+0) = faceDim; // face dimension
527 tagView(i*tagSize+1) = 0; // face id
528 tagView(i*tagSize+2) = i; // local dof id
529 tagView(i*tagSize+3) = cardinality; // total number of dofs on this face
530 }
531 }
532
533 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
534 // tags are constructed on host
536 this->ordinalToTag_,
537 tagView,
538 this->basisCardinality_,
539 tagSize,
540 posScDim,
541 posScOrd,
542 posDfOrd);
543 }
544 }
545
550 const char* getName() const override {
551 return "Intrepid2_IntegratedLegendreBasis_HGRAD_TRI";
552 }
553
556 virtual bool requireOrientation() const override {
557 return (this->getDegree() > 2);
558 }
559
560 // since the getValues() below only overrides the FEM variant, we specify that
561 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
562 // (It's an error to use the FVD variant on this basis.)
564
583 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
584 const EOperator operatorType = OPERATOR_VALUE ) const override
585 {
586 auto numPoints = inputPoints.extent_int(0);
587
589
590 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
591
592 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
593 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
594 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
595 const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
596
597 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
598 Kokkos::parallel_for("Hierarchical_HGRAD_TRI_Functor", policy, functor);
599 }
600
610 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
611 if(subCellDim == 1) {
612 return Teuchos::rcp(new
614 (this->basisDegree_));
615 }
616 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
617 }
618
624 getHostBasis() const override {
625 using HostDeviceType = typename Kokkos::HostSpace::device_type;
627 return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
628 }
629 };
630} // end namespace Intrepid2
631
632#endif /* Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h */
Header file for the abstract base class Intrepid2::Basis.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
KOKKOS_INLINE_FUNCTION void device_assert(bool val)
H(grad) basis on the line based on integrated Legendre polynomials.
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.
ordinal_type getDegree() const
Returns the degree of 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.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line: e...
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...
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
virtual bool requireOrientation() const override
True if orientation is required.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
IntegratedLegendreBasis_HGRAD_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Functor for computing values for the IntegratedLegendreBasis_HGRAD_TRI class.