Intrepid2
Intrepid2_HierarchicalBasis_HCURL_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_HierarchicalBasis_HCURL_TRI_h
50#define Intrepid2_HierarchicalBasis_HCURL_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 int numFields_, numPoints_;
87
88 size_t fad_size_output_;
89
90 static const int numVertices = 3;
91 static const int numEdges = 3;
92 static const int numFaceFamilies = 2;
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 const int face_family_start_[numFaceFamilies] = {0,1};
96 const int face_family_middle_[numFaceFamilies] = {1,2};
97 const int face_family_end_ [numFaceFamilies] = {2,0};
98
99 Hierarchical_HCURL_TRI_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints, int polyOrder)
100 : opType_(opType), output_(output), inputPoints_(inputPoints),
101 polyOrder_(polyOrder),
102 fad_size_output_(getScalarDimensionForView(output))
103 {
104 numFields_ = output.extent_int(0);
105 numPoints_ = output.extent_int(1);
106 const int expectedCardinality = 3 * polyOrder_ + polyOrder_ * (polyOrder-1);
107
108 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
109 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != expectedCardinality, std::invalid_argument, "output field size does not match basis cardinality");
110 }
111
112 KOKKOS_INLINE_FUNCTION
113 void operator()( const TeamMember & teamMember ) const
114 {
115 auto pointOrdinal = teamMember.league_rank();
116 OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
117 if (fad_size_output_ > 0) {
118 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
119 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
120 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
121 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
122 }
123 else {
124 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
125 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
126 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
127 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
128 }
129
130 const auto & x = inputPoints_(pointOrdinal,0);
131 const auto & y = inputPoints_(pointOrdinal,1);
132
133 // write as barycentric coordinates:
134 const PointScalar lambda[3] = {1. - x - y, x, y};
135 const PointScalar lambda_dx[3] = {-1., 1., 0.};
136 const PointScalar lambda_dy[3] = {-1., 0., 1.};
137
138 const int num1DEdgeFunctions = polyOrder_; // per edge
139
140 switch (opType_)
141 {
142 case OPERATOR_VALUE:
143 {
144 // edge functions
145 int fieldOrdinalOffset = 0;
146 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
147 {
148 const auto & s0 = lambda [edge_start_[edgeOrdinal]];
149 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
150 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
151
152 const auto & s1 = lambda [ edge_end_[edgeOrdinal]];
153 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
154 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
155
156 Polynomials::shiftedScaledLegendreValues(edge_field_values_at_point, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
157 for (int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
158 {
159 const auto & legendreValue = edge_field_values_at_point(edgeFunctionOrdinal);
160 const PointScalar xWeight = s0 * s1_dx - s1 * s0_dx;
161 const PointScalar yWeight = s0 * s1_dy - s1 * s0_dy;
162 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = legendreValue * xWeight;
163 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = legendreValue * yWeight;
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_ - 1;
174
175 // following ESEAS, we interleave the face families. This groups all the face dofs of a given degree together.
176 const int faceFieldOrdinalOffset = fieldOrdinalOffset;
177 for (int familyOrdinal=1; familyOrdinal<=2; familyOrdinal++)
178 {
179 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
180 const auto &s2 = lambda[ face_family_end_[familyOrdinal-1]];
181 for (int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
182 {
183 for (int i=0; i<ij_sum; i++)
184 {
185 const int j = ij_sum - i; // j >= 1
186 // family 1 involves edge functions from edge (0,1) (edgeOrdinal 0); family 2 involves functions from edge (1,2) (edgeOrdinal 1)
187 const int edgeBasisOrdinal = i + (familyOrdinal-1)*num1DEdgeFunctions;
188 const auto & edgeValue_x = output_(edgeBasisOrdinal,pointOrdinal,0);
189 const auto & edgeValue_y = output_(edgeBasisOrdinal,pointOrdinal,1);
190 const double alpha = i*2.0 + 1;
191
192 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_values_at_point, alpha, polyOrder_-1, s2, jacobiScaling);
193 const auto & jacobiValue = jacobi_values_at_point(j);
194 output_(fieldOrdinal,pointOrdinal,0) = edgeValue_x * jacobiValue;
195 output_(fieldOrdinal,pointOrdinal,1) = edgeValue_y * jacobiValue;
196
197 fieldOrdinal += 2; // 2 because there are two face families, and we interleave them.
198 }
199 }
200 }
201 }
202 } // end OPERATOR_VALUE
203 break;
204 case OPERATOR_CURL:
205 {
206 // edge functions
207 int fieldOrdinalOffset = 0;
208 /*
209 Per Fuentes et al. (see Appendix E.1, E.2), the curls of the edge functions, are
210 (i+2) * [P_i](s0,s1) * (grad s0 \times grad s1)
211 The P_i we have implemented in shiftedScaledLegendreValues.
212 */
213 // rename the scratch memory to match our usage here:
214 auto & P_i = edge_field_values_at_point;
215 auto & L_2ip1_j = jacobi_values_at_point;
216 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
217 {
218 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
219 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
220
221 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
222 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
223 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
224 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
225
226 const OutputScalar grad_s0_cross_grad_s1 = s0_dx * s1_dy - s1_dx * s0_dy;
227
228 Polynomials::shiftedScaledLegendreValues(P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
229 for (int i=0; i<num1DEdgeFunctions; i++)
230 {
231 output_(i+fieldOrdinalOffset,pointOrdinal) = (i+2) * P_i(i) * grad_s0_cross_grad_s1;
232 }
233 fieldOrdinalOffset += num1DEdgeFunctions;
234 }
235
236 /*
237 Fuentes et al give the face functions as E^f_{ij}, with curl:
238 [L^{2i+1}_j](s0+s1,s2) curl(E^E_i(s0,s1)) + grad[L^(2i+1)_j](s0+s1,s2) \times E^E_i(s0,s1)
239 where:
240 - E^E_i is the ith edge function on the edge s0 to s1
241 - L^{2i+1}_j is an shifted, scaled integrated Jacobi polynomial.
242 - For family 1, s0s1s2 = 012
243 - For family 2, s0s1s2 = 120
244 - Note that grad[L^(2i+1)_j](s0+s1,s2) is computed as [P^{2i+1}_{j-1}](s0+s1,s2) (grad s2) + [R^{2i+1}_{j-1}] grad (s0+s1+s2),
245 but for triangles (s0+s1+s2) is always 1, so that the grad (s0+s1+s2) is 0.
246 - Here,
247 [P^{2i+1}_{j-1}](s0,s1) = P^{2i+1}_{j-1}(s1,s0+s1)
248 and
249 [R^{2i+1}_{j-1}(s0,s1)] = d/dt L^{2i+1}_j(s1,s0+s1)
250 We have implemented P^{alpha}_{j} as shiftedScaledJacobiValues,
251 and d/dt L^{alpha}_{j} as shiftedScaledIntegratedJacobiValues_dt.
252 */
253 // rename the scratch memory to match our usage here:
254 auto & P_2ip1_j = other_values_at_point;
255
256 // following ESEAS, we interleave the face families. This groups all the face dofs of a given degree together.
257 const int faceFieldOrdinalOffset = fieldOrdinalOffset;
258 for (int familyOrdinal=1; familyOrdinal<=2; familyOrdinal++)
259 {
260 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
261
262 const auto &s0_index = face_family_start_ [familyOrdinal-1];
263 const auto &s1_index = face_family_middle_[familyOrdinal-1];
264 const auto &s2_index = face_family_end_ [familyOrdinal-1];
265 const auto &s0 = lambda[s0_index];
266 const auto &s1 = lambda[s1_index];
267 const auto &s2 = lambda[s2_index];
268 const double jacobiScaling = 1.0; // s0 + s1 + s2
269
270 const auto & s0_dx = lambda_dx[s0_index];
271 const auto & s0_dy = lambda_dy[s0_index];
272 const auto & s1_dx = lambda_dx[s1_index];
273 const auto & s1_dy = lambda_dy[s1_index];
274 const auto & s2_dx = lambda_dx[s2_index];
275 const auto & s2_dy = lambda_dy[s2_index];
276
277 const OutputScalar grad_s0_cross_grad_s1 = s0_dx * s1_dy - s1_dx * s0_dy;
278
279 Polynomials::shiftedScaledLegendreValues (P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
280 // [L^{2i+1}_j](s0+s1,s2) curl(E^E_i(s0,s1)) + grad[L^(2i+1)_j](s0+s1,s2) \times E^E_i(s0,s1)
281 // grad[L^(2i+1)_j](s0+s1,s2) \times E^E_i(s0,s1)
282// - Note that grad[L^(2i+1)_j](s0+s1,s2) is computed as [P^{2i+1}_{j-1}](s0+s1,s2) (grad s2) + [R^{2i+1}_{j-1}] grad (s0+s1+s2),
283 const PointScalar xEdgeWeight = s0 * s1_dx - s1 * s0_dx;
284 const PointScalar yEdgeWeight = s0 * s1_dy - s1 * s0_dy;
285 OutputScalar grad_s2_cross_xy_edgeWeight = s2_dx * yEdgeWeight - xEdgeWeight * s2_dy;
286
287 const int max_ij_sum = polyOrder_ - 1;
288 for (int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
289 {
290 for (int i=0; i<ij_sum; i++)
291 {
292 const int j = ij_sum - i; // j >= 1
293 const OutputScalar edgeCurl = (i+2.) * P_i(i) * grad_s0_cross_grad_s1;
294
295 const double alpha = i*2.0 + 1;
296
297 Polynomials::shiftedScaledJacobiValues(P_2ip1_j, alpha, polyOrder_-1, PointScalar(s2), jacobiScaling);
298 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ip1_j, alpha, polyOrder_-1, s2, jacobiScaling);
299
300 const PointScalar & edgeValue = P_i(i);
301 output_(fieldOrdinal,pointOrdinal) = L_2ip1_j(j) * edgeCurl + P_2ip1_j(j-1) * edgeValue * grad_s2_cross_xy_edgeWeight;
302
303 fieldOrdinal += 2; // 2 because there are two face families, and we interleave them.
304 }
305 }
306 }
307 }
308 break;
309 case OPERATOR_GRAD:
310 case OPERATOR_D1:
311 case OPERATOR_D2:
312 case OPERATOR_D3:
313 case OPERATOR_D4:
314 case OPERATOR_D5:
315 case OPERATOR_D6:
316 case OPERATOR_D7:
317 case OPERATOR_D8:
318 case OPERATOR_D9:
319 case OPERATOR_D10:
320 INTREPID2_TEST_FOR_ABORT( true,
321 ">>> ERROR: (Intrepid2::Hierarchical_HCURL_TRI_Functor) Unsupported differential operator");
322 default:
323 // unsupported operator type
324 device_assert(false);
325 }
326 }
327
328 // Provide the shared memory capacity.
329 // This function takes the team_size as an argument,
330 // which allows team_size-dependent allocations.
331 size_t team_shmem_size (int team_size) const
332 {
333 // we will use shared memory to create a fast buffer for basis computations
334 size_t shmem_size = 0;
335 if (fad_size_output_ > 0)
336 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
337 else
338 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);
339
340 return shmem_size;
341 }
342 };
343
358 template<typename DeviceType,
359 typename OutputScalar = double,
360 typename PointScalar = double,
361 bool useCGBasis = true> // if useCGBasis is false, all basis functions will be associated with the interior
363 : public Basis<DeviceType,OutputScalar,PointScalar>
364 {
365 public:
367
369
372
373 using typename BasisBase::OutputViewType;
374 using typename BasisBase::PointViewType;
375 using typename BasisBase::ScalarViewType;
376
377 using typename BasisBase::ExecutionSpace;
378
379 protected:
380 int polyOrder_; // the maximum order of the polynomial
381 public:
386 HierarchicalBasis_HCURL_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
387 :
388 polyOrder_(polyOrder)
389 {
390 const int numEdgeFunctions = polyOrder * 3;
391 const int numFaceFunctions = polyOrder * (polyOrder-1); // two families, each with p*(p-1)/2 functions
392 this->basisCardinality_ = numEdgeFunctions + numFaceFunctions;
393 this->basisDegree_ = polyOrder;
394 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
395 this->basisType_ = BASIS_FEM_HIERARCHICAL;
396 this->basisCoordinates_ = COORDINATES_CARTESIAN;
397 this->functionSpace_ = FUNCTION_SPACE_HCURL;
398
399 const int degreeLength = 1;
400 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Hierarchical H(curl) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);
401 this->fieldOrdinalH1PolynomialDegree_ = OrdinalTypeArray2DHost("Hierarchical H(curl) triangle polynomial H^1 degree lookup", this->basisCardinality_, degreeLength);
402
403 int fieldOrdinalOffset = 0;
404 // **** vertex functions **** //
405 // no vertex functions in H(curl)
406
407 // **** edge functions **** //
408 const int numFunctionsPerEdge = polyOrder; // p functions associated with each edge
409 const int numEdges = this->basisCellTopology_.getEdgeCount();
410 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
411 {
412 for (int i=0; i<numFunctionsPerEdge; i++)
413 {
414 this->fieldOrdinalPolynomialDegree_(i+fieldOrdinalOffset,0) = i+1; // the multiplicands involving the gradients of the vertex functions are first degree polynomials; hence the +1 (the remaining multiplicands are order i = 0,…,p-1).
415 this->fieldOrdinalH1PolynomialDegree_(i+fieldOrdinalOffset,0) = i+2;
416 }
417 fieldOrdinalOffset += numFunctionsPerEdge;
418 }
419 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != numEdgeFunctions, std::invalid_argument, "Internal error: basis enumeration is incorrect");
420
421 // **** face functions **** //
422 const int max_ij_sum = polyOrder-1;
423 const int faceFieldOrdinalOffset = fieldOrdinalOffset;
424 for (int faceFamilyOrdinal=1; faceFamilyOrdinal<=2; faceFamilyOrdinal++)
425 {
426 // following ESEAS, we interleave the face families. This groups all the face dofs of a given degree together.
427 int fieldOrdinal = faceFieldOrdinalOffset + faceFamilyOrdinal - 1;
428 for (int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
429 {
430 for (int i=0; i<ij_sum; i++)
431 {
432 this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = ij_sum+1;
433 this->fieldOrdinalH1PolynomialDegree_(fieldOrdinal,0) = ij_sum+2;
434 fieldOrdinal += 2; // 2 because there are two face families, and we interleave them.
435 }
436 }
437 fieldOrdinalOffset = fieldOrdinal - 1; // due to the interleaving increment, we've gone two past the last face ordinal. Set offset to be one past.
438 }
439
440 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
441
442 // initialize tags
443 {
444 const auto & cardinality = this->basisCardinality_;
445
446 // Basis-dependent initializations
447 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
448 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
449 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
450 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
451
452 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
453 const ordinal_type edgeDim = 1, faceDim = 2;
454
455 if (useCGBasis) {
456 {
457 int tagNumber = 0;
458 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
459 {
460 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
461 {
462 tagView(tagNumber*tagSize+0) = edgeDim; // edge dimension
463 tagView(tagNumber*tagSize+1) = edgeOrdinal; // edge id
464 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
465 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge; // total number of dofs on this edge
466 tagNumber++;
467 }
468 }
469 const int numFunctionsPerFace = numFaceFunctions; // just one face in the triangle
470 const int numFaces = 1;
471 for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
472 {
473 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
474 {
475 tagView(tagNumber*tagSize+0) = faceDim; // face dimension
476 tagView(tagNumber*tagSize+1) = faceOrdinal; // face id
477 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
478 tagView(tagNumber*tagSize+3) = numFunctionsPerFace; // total number of dofs on this face
479 tagNumber++;
480 }
481 }
482 }
483 }
484 else
485 {
486 // DG basis: all functions are associated with interior
487 for (ordinal_type i=0;i<cardinality;++i) {
488 tagView(i*tagSize+0) = faceDim; // face dimension
489 tagView(i*tagSize+1) = 0; // face id
490 tagView(i*tagSize+2) = i; // local dof id
491 tagView(i*tagSize+3) = cardinality; // total number of dofs on this face
492 }
493 }
494
495 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
496 // tags are constructed on host
498 this->ordinalToTag_,
499 tagView,
500 this->basisCardinality_,
501 tagSize,
502 posScDim,
503 posScOrd,
504 posDfOrd);
505 }
506 }
507
512 const char* getName() const override {
513 return "Intrepid2_HierarchicalBasis_HCURL_TRI";
514 }
515
518 virtual bool requireOrientation() const override {
519 return (this->getDegree() > 2);
520 }
521
522 // since the getValues() below only overrides the FEM variant, we specify that
523 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
524 // (It's an error to use the FVD variant on this basis.)
526
545 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
546 const EOperator operatorType = OPERATOR_VALUE ) const override
547 {
548 auto numPoints = inputPoints.extent_int(0);
549
551
552 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);
553
554 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
555 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
556 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
557 const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
558
559 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
560 Kokkos::parallel_for("Hierarchical_HCURL_TRI_Functor", policy, functor);
561 }
562
572 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
573 if(subCellDim == 1) {
574 return Teuchos::rcp(new
576 (this->basisDegree_-1));
577 }
578 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
579 }
580
586 getHostBasis() const override {
587 using HostDeviceType = typename Kokkos::HostSpace::device_type;
589 return Teuchos::rcp( new HostBasisType(polyOrder_) );
590 }
591 };
592} // end namespace Intrepid2
593
594#endif /* Intrepid2_HierarchicalBasis_HCURL_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.
For mathematical details of the construction, see:
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 void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
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...
const char * getName() const override
Returns basis name.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
HierarchicalBasis_HCURL_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual bool requireOrientation() const override
True if orientation is required.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line.
Functor for computing values for the HierarchicalBasis_HCURL_TRI class.