Intrepid2
Intrepid2_PointToolsDef.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), or
38// Mauro Perego (mperego@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
48#ifndef __INTREPID2_POINTTOOLS_DEF_HPP__
49#define __INTREPID2_POINTTOOLS_DEF_HPP__
50
51#if defined(_MSC_VER) || defined(_WIN32) && defined(__ICL)
52// M_PI, M_SQRT2, etc. are hidden in MSVC by #ifdef _USE_MATH_DEFINES
53 #ifndef _USE_MATH_DEFINES
54 #define _USE_MATH_DEFINES
55 #endif
56 #include <math.h>
57#endif
58
59namespace Intrepid2 {
60
61// -----------------------------------------------------------------------------------------
62// Front interface
63// -----------------------------------------------------------------------------------------
64
65inline
66ordinal_type
68getLatticeSize( const shards::CellTopology cellType,
69 const ordinal_type order,
70 const ordinal_type offset ) {
71#ifdef HAVE_INTREPID2_DEBUG
72 INTREPID2_TEST_FOR_EXCEPTION( order < 0 || offset < 0,
73 std::invalid_argument ,
74 ">>> ERROR (PointTools::getLatticeSize): order and offset must be positive values." );
75#endif
76 ordinal_type r_val = 0;
77 switch (cellType.getBaseKey()) {
78 case shards::Tetrahedron<>::key: {
79 const auto effectiveOrder = order - 4 * offset;
80 r_val = (effectiveOrder < 0 ? 0 :(effectiveOrder+1)*(effectiveOrder+2)*(effectiveOrder+3)/6);
81 break;
82 }
83 case shards::Triangle<>::key: {
84 const auto effectiveOrder = order - 3 * offset;
85 r_val = (effectiveOrder < 0 ? 0 : (effectiveOrder+1)*(effectiveOrder+2)/2);
86 break;
87 }
88 case shards::Line<>::key: {
89 const auto effectiveOrder = order - 2 * offset;
90 r_val = (effectiveOrder < 0 ? 0 : (effectiveOrder+1));
91 break;
92 }
93 case shards::Quadrilateral<>::key: {
94 const auto effectiveOrder = order - 2 * offset;
95 r_val = std::pow(effectiveOrder < 0 ? 0 : (effectiveOrder+1),2);
96 break;
97 }
98 case shards::Hexahedron<>::key: {
99 const auto effectiveOrder = order - 2 * offset;
100 r_val = std::pow(effectiveOrder < 0 ? 0 : (effectiveOrder+1),3);
101 break;
102 }
103 case shards::Pyramid<>::key: {
104 const auto effectiveOrder = order - 2 * offset;
105 r_val = (effectiveOrder < 0 ? 0 : (effectiveOrder+1)*(effectiveOrder+2)*(2*effectiveOrder+3)/6);
106 break;
107 }
108 default: {
109 INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
110 ">>> ERROR (Intrepid2::PointTools::getLatticeSize): the specified cell type is not supported." );
111 }
112 }
113 return r_val;
114}
115
116template<typename pointValueType, class ...pointProperties>
117void
119getLattice( Kokkos::DynRankView<pointValueType,pointProperties...> points,
120 const shards::CellTopology cell,
121 const ordinal_type order,
122 const ordinal_type offset,
123 const EPointType pointType ) {
124#ifdef HAVE_INTREPID2_DEBUG
125 INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 2,
126 std::invalid_argument ,
127 ">>> ERROR (PointTools::getLattice): points rank must be 2." );
128 INTREPID2_TEST_FOR_EXCEPTION( order < 0 || offset < 0,
129 std::invalid_argument ,
130 ">>> ERROR (PointTools::getLattice): order and offset must be positive values." );
131
132 const size_type latticeSize = getLatticeSize( cell, order, offset );
133 const size_type spaceDim = cell.getDimension();
134
135 INTREPID2_TEST_FOR_EXCEPTION( points.extent(0) != latticeSize ||
136 points.extent(1) != spaceDim,
137 std::invalid_argument ,
138 ">>> ERROR (PointTools::getLattice): dimension does not match to lattice size." );
139#endif
140
141 switch (cell.getBaseKey()) {
142 case shards::Tetrahedron<>::key: getLatticeTetrahedron( points, order, offset, pointType ); break;
143 case shards::Pyramid<>::key: getLatticePyramid ( points, order, offset, pointType ); break;
144 case shards::Triangle<>::key: getLatticeTriangle ( points, order, offset, pointType ); break;
145 case shards::Line<>::key: getLatticeLine ( points, order, offset, pointType ); break;
146 case shards::Quadrilateral<>::key: {
147 auto hostPoints = Kokkos::create_mirror_view(points);
148 shards::CellTopology line(shards::getCellTopologyData<shards::Line<2> >());
149 const ordinal_type numPoints = getLatticeSize( line, order, offset );
150 auto linePoints = getMatchingViewWithLabel(hostPoints, "linePoints", numPoints, 1);
151 getLatticeLine( linePoints, order, offset, pointType );
152 ordinal_type idx=0;
153 for (ordinal_type j=0; j<numPoints; ++j) {
154 for (ordinal_type i=0; i<numPoints; ++i, ++idx) {
155 hostPoints(idx,0) = linePoints(i,0);
156 hostPoints(idx,1) = linePoints(j,0);
157 }
158 }
159 Kokkos::deep_copy(points,hostPoints);
160 }
161 break;
162 case shards::Hexahedron<>::key: {
163 auto hostPoints = Kokkos::create_mirror_view(points);
164 shards::CellTopology line(shards::getCellTopologyData<shards::Line<2> >());
165 const ordinal_type numPoints = getLatticeSize( line, order, offset );
166 auto linePoints = getMatchingViewWithLabel(hostPoints, "linePoints", numPoints, 1);
167 getLatticeLine( linePoints, order, offset, pointType );
168 ordinal_type idx=0;
169 for (ordinal_type k=0; k<numPoints; ++k) {
170 for (ordinal_type j=0; j<numPoints; ++j) {
171 for (ordinal_type i=0; i<numPoints; ++i, ++idx) {
172 hostPoints(idx,0) = linePoints(i,0);
173 hostPoints(idx,1) = linePoints(j,0);
174 hostPoints(idx,2) = linePoints(k,0);
175 }
176 }
177 }
178 Kokkos::deep_copy(points,hostPoints);
179 }
180 break;
181 default: {
182 INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
183 ">>> ERROR (Intrepid2::PointTools::getLattice): the specified cell type is not supported." );
184 }
185 }
186}
187
188template<typename pointValueType, class ...pointProperties>
189void
191getGaussPoints( Kokkos::DynRankView<pointValueType,pointProperties...> points,
192 const ordinal_type order ) {
193#ifdef HAVE_INTREPID2_DEBUG
194 INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 2,
195 std::invalid_argument ,
196 ">>> ERROR (PointTools::getGaussPoints): points rank must be 1." );
197 INTREPID2_TEST_FOR_EXCEPTION( order < 0,
198 std::invalid_argument ,
199 ">>> ERROR (PointTools::getGaussPoints): order must be positive value." );
200#endif
201 const ordinal_type np = order + 1;
202 const double alpha = 0.0, beta = 0.0;
203
204 // until view and dynrankview inter-operatible, we use views in a consistent way
205 Kokkos::View<pointValueType*,Kokkos::HostSpace>
206 zHost("PointTools::getGaussPoints::z", np),
207 wHost("PointTools::getGaussPoints::w", np);
208
209 // sequential means that the code is decorated with KOKKOS_INLINE_FUNCTION
210 // and it does not invoke parallel for inside (cheap operation), which means
211 // that gpu memory is not accessible unless this is called inside of functor.
212 Polylib::Serial::Cubature<POLYTYPE_GAUSS>::getValues(zHost, wHost, np, alpha, beta);
213
214 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
215 auto pts = Kokkos::subview( points, range_type(0,np), 0 );
216 // should be fixed after view and dynrankview are inter-operatible
217 auto z = Kokkos::DynRankView<pointValueType,Kokkos::HostSpace>(zHost.data(), np);
218 Kokkos::deep_copy(pts, z);
219}
220
221// -----------------------------------------------------------------------------------------
222// Internal implementation
223// -----------------------------------------------------------------------------------------
224
225template<typename pointValueType, class ...pointProperties>
226void
228getLatticeLine( Kokkos::DynRankView<pointValueType,pointProperties...> points,
229 const ordinal_type order,
230 const ordinal_type offset,
231 const EPointType pointType ) {
232 switch (pointType) {
233 case POINTTYPE_EQUISPACED: getEquispacedLatticeLine( points, order, offset ); break;
234 case POINTTYPE_WARPBLEND: getWarpBlendLatticeLine( points, order, offset ); break;
235 default: {
236 INTREPID2_TEST_FOR_EXCEPTION( true ,
237 std::invalid_argument ,
238 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
239 }
240 }
241}
242
243template<typename pointValueType, class ...pointProperties>
244void
246getLatticeTriangle( Kokkos::DynRankView<pointValueType,pointProperties...> points,
247 const ordinal_type order,
248 const ordinal_type offset,
249 const EPointType pointType ) {
250 switch (pointType) {
251 case POINTTYPE_EQUISPACED: getEquispacedLatticeTriangle( points, order, offset ); break;
252 case POINTTYPE_WARPBLEND: getWarpBlendLatticeTriangle ( points, order, offset ); break;
253 default: {
254 INTREPID2_TEST_FOR_EXCEPTION( true ,
255 std::invalid_argument ,
256 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
257 }
258 }
259}
260
261template<typename pointValueType, class ...pointProperties>
262void
264getLatticeTetrahedron( Kokkos::DynRankView<pointValueType,pointProperties...> points,
265 const ordinal_type order,
266 const ordinal_type offset,
267 const EPointType pointType ) {
268 switch (pointType) {
269 case POINTTYPE_EQUISPACED: getEquispacedLatticeTetrahedron( points, order, offset ); break;
270 case POINTTYPE_WARPBLEND: getWarpBlendLatticeTetrahedron ( points, order, offset ); break;
271 default: {
272 INTREPID2_TEST_FOR_EXCEPTION( true ,
273 std::invalid_argument ,
274 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
275 }
276 }
277}
278
279template<typename pointValueType, class ...pointProperties>
280void
282getLatticePyramid( Kokkos::DynRankView<pointValueType,pointProperties...> points,
283 const ordinal_type order,
284 const ordinal_type offset,
285 const EPointType pointType )
286{
287 switch (pointType) {
288 case POINTTYPE_EQUISPACED: getEquispacedLatticePyramid( points, order, offset ); break;
289 default: {
290 INTREPID2_TEST_FOR_EXCEPTION( true ,
291 std::invalid_argument ,
292 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
293 }
294 }
295}
296
297// -----------------------------------------------------------------------------------------
298
299template<typename pointValueType, class ...pointProperties>
300void
302getEquispacedLatticeLine( Kokkos::DynRankView<pointValueType,pointProperties...> points,
303 const ordinal_type order,
304 const ordinal_type offset ) {
305 auto pointsHost = Kokkos::create_mirror_view(points);
306
307 if (order == 0)
308 pointsHost(0,0) = 0.0;
309 else {
310 const pointValueType h = 2.0 / order;
311 const ordinal_type ibeg = offset, iend = order-offset+1;
312 for (ordinal_type i=ibeg;i<iend;++i)
313 pointsHost(i-ibeg, 0) = -1.0 + h * (pointValueType) i;
314 }
315
316 Kokkos::deep_copy(points, pointsHost);
317}
318
319template<typename pointValueType, class ...pointProperties>
320void
322getWarpBlendLatticeLine( Kokkos::DynRankView<pointValueType,pointProperties...> points,
323 const ordinal_type order,
324 const ordinal_type offset ) {
325 // order is order of polynomial degree. The Gauss-Lobatto points are accurate
326 // up to degree 2*i-1
327 const ordinal_type np = order + 1;
328 const double alpha = 0.0, beta = 0.0;
329 const ordinal_type s = np - 2*offset;
330
331 if (s > 0) {
332 // until view and dynrankview inter-operatible, we use views in a consistent way
333 Kokkos::View<pointValueType*,Kokkos::HostSpace>
334 zHost("PointTools::getGaussPoints::z", np),
335 wHost("PointTools::getGaussPoints::w", np);
336
337 // sequential means that the code is decorated with KOKKOS_INLINE_FUNCTION
338 // and it does not invoke parallel for inside (cheap operation), which means
339 // that gpu memory is not accessible unless this is called inside of functor.
340 Polylib::Serial::Cubature<POLYTYPE_GAUSS_LOBATTO>::getValues(zHost, wHost, np, alpha, beta);
341
342 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
343 auto pts = Kokkos::subview( points, range_type(0, s), 0 );
344
345 // this should be fixed after view and dynrankview is interoperatable
346 auto z = Kokkos::DynRankView<pointValueType,Kokkos::HostSpace>(zHost.data() + offset, np-offset);
347
348 const auto common_range = range_type(0, std::min(pts.extent(0), z.extent(0)));
349 Kokkos::deep_copy(Kokkos::subview(pts, common_range),
350 Kokkos::subview(z, common_range));
351 }
352}
353
354template<typename pointValueType, class ...pointProperties>
355void
357getEquispacedLatticeTriangle( Kokkos::DynRankView<pointValueType,pointProperties...> points,
358 const ordinal_type order,
359 const ordinal_type offset ) {
360 TEUCHOS_TEST_FOR_EXCEPTION( order <= 0 ,
361 std::invalid_argument ,
362 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTriangle): order must be positive" );
363
364 auto pointsHost = Kokkos::create_mirror_view(points);
365
366 const pointValueType h = 1.0 / order;
367 ordinal_type cur = 0;
368
369 for (ordinal_type i=offset;i<=order-offset;i++) {
370 for (ordinal_type j=offset;j<=order-i-offset;j++) {
371 pointsHost(cur,0) = j * h ;
372 pointsHost(cur,1) = i * h;
373 cur++;
374 }
375 }
376 Kokkos::deep_copy(points, pointsHost);
377}
378
379template<typename pointValueType, class ...pointProperties>
380void
382getEquispacedLatticeTetrahedron( Kokkos::DynRankView<pointValueType,pointProperties...> points,
383 const ordinal_type order ,
384 const ordinal_type offset ) {
385 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
386 std::invalid_argument ,
387 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTetrahedron): order must be positive" );
388
389 auto pointsHost = Kokkos::create_mirror_view(points);
390
391 const pointValueType h = 1.0 / order;
392 ordinal_type cur = 0;
393
394 for (ordinal_type i=offset;i<=order-offset;i++) {
395 for (ordinal_type j=offset;j<=order-i-offset;j++) {
396 for (ordinal_type k=offset;k<=order-i-j-offset;k++) {
397 pointsHost(cur,0) = k * h;
398 pointsHost(cur,1) = j * h;
399 pointsHost(cur,2) = i * h;
400 cur++;
401 }
402 }
403 }
404 Kokkos::deep_copy(points, pointsHost);
405}
406
407template<typename pointValueType, class ...pointProperties>
408void
410getEquispacedLatticePyramid( Kokkos::DynRankView<pointValueType,pointProperties...> points,
411 const ordinal_type order ,
412 const ordinal_type offset ) {
413 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
414 std::invalid_argument ,
415 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticePyramid): order must be positive" );
416
417 auto pointsHost = Kokkos::create_mirror_view(points);
418
419 const pointValueType h = 1.0 / order;
420 ordinal_type cur = 0;
421
422 for (ordinal_type i=offset;i<=order-offset;i++) { // z dimension (goes from 0 to 1)
423 pointValueType z = i * h;
424 for (ordinal_type j=offset;j<=order-i-offset;j++) { // y (goes from -(1-z) to 1-z)
425 for (ordinal_type k=offset;k<=order-i-offset;k++) { // x (goes from -(1-z) to 1-z)
426 pointsHost(cur,0) = (2 * k * h - 1.) * (1-z);
427 pointsHost(cur,1) = (2 * j * h - 1.) * (1-z);
428 pointsHost(cur,2) = z;
429 cur++;
430 }
431 }
432 }
433 Kokkos::deep_copy(points, pointsHost);
434}
435
436template<typename pointValueType, class ...pointProperties>
437void
439warpFactor( Kokkos::DynRankView<pointValueType,pointProperties...> warp,
440 const ordinal_type order,
441 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes ,
442 const Kokkos::DynRankView<pointValueType,pointProperties...> xout
443 )
444 {
445 TEUCHOS_TEST_FOR_EXCEPTION( ( warp.extent(0) != xout.extent(0) ) ,
446 std::invalid_argument ,
447 ">>> ERROR (PointTools::warpFactor): xout and warp must be same size." );
448
449 Kokkos::deep_copy(warp, pointValueType(0.0));
450
451 ordinal_type xout_dim0 = xout.extent(0);
452 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> d("d", xout_dim0 );
453
454 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> xeq_("xeq", order + 1 ,1);
455 PointTools::getEquispacedLatticeLine( xeq_ , order , 0 );
456 const auto xeq = Kokkos::subview(xeq_, Kokkos::ALL(),0);
457
458 TEUCHOS_TEST_FOR_EXCEPTION( ( xeq.extent(0) != xnodes.extent(0) ) ,
459 std::invalid_argument ,
460 ">>> ERROR (PointTools::warpFactor): xeq and xnodes must be same size." );
461
462
463 for (ordinal_type i=0;i<=order;i++) {
464
465 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
466
467 for (ordinal_type j=1;j<order;j++) {
468 if (i != j) {
469 for (ordinal_type k=0;k<xout_dim0;k++) {
470 d(k) = d(k) * ( (xout(k)-xeq(j)) / (xeq(i)-xeq(j)) );
471 }
472 }
473 }
474
475 // deflate end roots
476 if ( i != 0 ) {
477 for (ordinal_type k=0;k<xout_dim0;k++) {
478 d(k) = -d(k) / (xeq(i) - xeq(0));
479 }
480 }
481
482 if (i != order ) {
483 for (ordinal_type k=0;k<xout_dim0;k++) {
484 d(k) = d(k) / (xeq(i) - xeq(order));
485 }
486 }
487
488 for (ordinal_type k=0;k<xout_dim0;k++) {
489 warp(k) += d(k);
490 }
491
492 }
493 }
494
495template<typename pointValueType, class ...pointProperties>
496void
498getWarpBlendLatticeTriangle( Kokkos::DynRankView<pointValueType,pointProperties...> points,
499 const ordinal_type order ,
500 const ordinal_type offset)
501 {
502 /* get Gauss-Lobatto points */
503
504 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> gaussX("gaussX", order + 1 , 1 );
505
506 PointTools::getWarpBlendLatticeLine( gaussX , order , 0 );
507 //auto gaussX = Kokkos::subdynrankview(gaussX_, Kokkos::ALL(),0);
508
509 // gaussX.resize(gaussX.extent(0));
510
511 pointValueType alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999,
512 1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
513
514 pointValueType alpha;
515
516 if (order >= 1 && order < 16) {
517 alpha = alpopt[order-1];
518 }
519 else {
520 alpha = 5.0 / 3.0;
521 }
522
523 const ordinal_type p = order; /* switch to Warburton's notation */
524 ordinal_type N = (p+1)*(p+2)/2;
525
526 /* equidistributed nodes on equilateral triangle */
527 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1("L1", N );
528 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2("L2", N );
529 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3("L3", N );
530 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> X("X", N);
531 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> Y("Y", N);
532
533 ordinal_type sk = 0;
534 for (ordinal_type n=1;n<=p+1;n++) {
535 for (ordinal_type m=1;m<=p+2-n;m++) {
536 L1(sk) = (n-1) / (pointValueType)p;
537 L3(sk) = (m-1) / (pointValueType)p;
538 L2(sk) = 1.0 - L1(sk) - L3(sk);
539 sk++;
540 }
541 }
542
543 for (ordinal_type n=0;n<N;n++) {
544 X(n) = -L2(n) + L3(n);
545 Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772;
546 }
547
548 /* get blending function for each node at each edge */
549 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend1("blend1", N);
550 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend2("blend2", N);
551 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend3("blend3", N);
552
553 for (ordinal_type n=0;n<N;n++) {
554 blend1(n) = 4.0 * L2(n) * L3(n);
555 blend2(n) = 4.0 * L1(n) * L3(n);
556 blend3(n) = 4.0 * L1(n) * L2(n);
557 }
558
559 /* get difference of each barycentric coordinate */
560 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3mL2("L3mL2", N);
561 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1mL3("L1mL3", N);
562 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2mL1("L2mL1", N);
563
564 for (ordinal_type k=0;k<N;k++) {
565 L3mL2(k) = L3(k)-L2(k);
566 L1mL3(k) = L1(k)-L3(k);
567 L2mL1(k) = L2(k)-L1(k);
568 }
569
570 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor1("warpfactor1", N);
571 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor2("warpfactor2", N);
572 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor3("warpfactor3", N);
573
574 warpFactor( warpfactor1, order , gaussX , L3mL2 );
575 warpFactor( warpfactor2, order , gaussX , L1mL3 );
576 warpFactor( warpfactor3, order , gaussX , L2mL1 );
577
578 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp1("warp1", N);
579 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp2("warp2", N);
580 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp3("warp3", N);
581
582 for (ordinal_type k=0;k<N;k++) {
583 warp1(k) = blend1(k) * warpfactor1(k) *
584 ( 1.0 + alpha * alpha * L1(k) * L1(k) );
585 warp2(k) = blend2(k) * warpfactor2(k) *
586 ( 1.0 + alpha * alpha * L2(k) * L2(k) );
587 warp3(k) = blend3(k) * warpfactor3(k) *
588 ( 1.0 + alpha * alpha * L3(k) * L3(k) );
589 }
590
591 for (ordinal_type k=0;k<N;k++) {
592 X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k);
593 Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k);
594 }
595
596 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warXY("warXY", 1, N,2);
597
598 for (ordinal_type k=0;k<N;k++) {
599 warXY(0, k,0) = X(k);
600 warXY(0, k,1) = Y(k);
601 }
602
603
604 // finally, convert the warp-blend points to the correct triangle
605 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warburtonVerts("warburtonVerts", 1,3,2);
606 warburtonVerts(0,0,0) = -1.0;
607 warburtonVerts(0,0,1) = -1.0/std::sqrt(3.0);
608 warburtonVerts(0,1,0) = 1.0;
609 warburtonVerts(0,1,1) = -1.0/std::sqrt(3.0);
610 warburtonVerts(0,2,0) = 0.0;
611 warburtonVerts(0,2,1) = 2.0/std::sqrt(3.0);
612
613 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> refPts("refPts", 1, N,2);
614
616 warXY ,
617 warburtonVerts ,
618 shards::getCellTopologyData< shards::Triangle<3> >()
619 );
620
621
622 auto pointsHost = Kokkos::create_mirror_view(points);
623 // now write from refPts into points, taking care of offset
624 ordinal_type noffcur = 0; // index into refPts
625 ordinal_type offcur = 0; // index ordinal_type points
626 for (ordinal_type i=0;i<=order;i++) {
627 for (ordinal_type j=0;j<=order-i;j++) {
628 if ( (i >= offset) && (i <= order-offset) &&
629 (j >= offset) && (j <= order-i-offset) ) {
630 pointsHost(offcur,0) = refPts(0, noffcur,0);
631 pointsHost(offcur,1) = refPts(0, noffcur,1);
632 offcur++;
633 }
634 noffcur++;
635 }
636 }
637
638 Kokkos::deep_copy(points, pointsHost);
639
640 }
641
642
643template<typename pointValueType, class ...pointProperties>
644void
646warpShiftFace3D( Kokkos::DynRankView<pointValueType,pointProperties...> dxy,
647 const ordinal_type order ,
648 const pointValueType pval ,
649 const Kokkos::DynRankView<pointValueType,pointProperties...> /* L1 */,
650 const Kokkos::DynRankView<pointValueType,pointProperties...> L2,
651 const Kokkos::DynRankView<pointValueType,pointProperties...> L3,
652 const Kokkos::DynRankView<pointValueType,pointProperties...> L4
653 )
654 {
655 evalshift(dxy,order,pval,L2,L3,L4);
656 return;
657 }
658
659template<typename pointValueType, class ...pointProperties>
660void
662evalshift( Kokkos::DynRankView<pointValueType,pointProperties...> dxy,
663 const ordinal_type order ,
664 const pointValueType pval ,
665 const Kokkos::DynRankView<pointValueType,pointProperties...> L1 ,
666 const Kokkos::DynRankView<pointValueType,pointProperties...> L2 ,
667 const Kokkos::DynRankView<pointValueType,pointProperties...> L3
668 )
669 {
670 // get Gauss-Lobatto-nodes
671 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> gaussX("gaussX",order+1,1);
672 PointTools::getWarpBlendLatticeLine( gaussX , order , 0 );
673 //gaussX.resize(order+1);
674 const ordinal_type N = L1.extent(0);
675
676 // Warburton code reverses them
677 for (ordinal_type k=0;k<=order;k++) {
678 gaussX(k,0) *= -1.0;
679 }
680
681 // blending function at each node for each edge
682 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend1("blend1",N);
683 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend2("blend2",N);
684 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend3("blend3",N);
685
686 for (ordinal_type i=0;i<N;i++) {
687 blend1(i) = L2(i) * L3(i);
688 blend2(i) = L1(i) * L3(i);
689 blend3(i) = L1(i) * L2(i);
690 }
691
692 // amount of warp for each node for each edge
693 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor1("warpfactor1",N);
694 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor2("warpfactor2",N);
695 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor3("warpfactor3",N);
696
697 // differences of barycentric coordinates
698 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3mL2("L3mL2",N);
699 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1mL3("L1mL3",N);
700 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2mL1("L2mL1",N);
701
702 for (ordinal_type k=0;k<N;k++) {
703 L3mL2(k) = L3(k)-L2(k);
704 L1mL3(k) = L1(k)-L3(k);
705 L2mL1(k) = L2(k)-L1(k);
706 }
707
708 evalwarp( warpfactor1 , order , gaussX , L3mL2 );
709 evalwarp( warpfactor2 , order , gaussX , L1mL3 );
710 evalwarp( warpfactor3 , order , gaussX , L2mL1 );
711
712 for (ordinal_type k=0;k<N;k++) {
713 warpfactor1(k) *= 4.0;
714 warpfactor2(k) *= 4.0;
715 warpfactor3(k) *= 4.0;
716 }
717
718 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp1("warp1",N);
719 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp2("warp2",N);
720 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp3("warp3",N);
721
722 for (ordinal_type k=0;k<N;k++) {
723 warp1(k) = blend1(k) * warpfactor1(k) *
724 ( 1.0 + pval * pval * L1(k) * L1(k) );
725 warp2(k) = blend2(k) * warpfactor2(k) *
726 ( 1.0 + pval * pval * L2(k) * L2(k) );
727 warp3(k) = blend3(k) * warpfactor3(k) *
728 ( 1.0 + pval * pval * L3(k) * L3(k) );
729 }
730
731 for (ordinal_type k=0;k<N;k++) {
732 dxy(k,0) = 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos( 4.0*M_PI/3.0 ) * warp3(k);
733 dxy(k,1) = 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4.0*M_PI/3.0 ) * warp3(k);
734 }
735 }
736
737 /* one-d edge warping function */
738template<typename pointValueType, class ...pointProperties>
739void
741evalwarp(Kokkos::DynRankView<pointValueType,pointProperties...> warp ,
742 const ordinal_type order ,
743 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes,
744 const Kokkos::DynRankView<pointValueType,pointProperties...> xout )
745 {
746 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> xeq("xeq",order+1);
747
748 ordinal_type xout_dim0 = xout.extent(0);
749 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> d("d",xout_dim0);
750
751 //Kokkos::deep_copy(d, 0.0);
752
753 for (ordinal_type i=0;i<=order;i++) {
754 xeq(i) = -1.0 + 2.0 * ( order - i ) / order;
755 }
756
757 for (ordinal_type i=0;i<=order;i++) {
758 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
759 for (ordinal_type j=1;j<order;j++) {
760 if (i!=j) {
761 for (ordinal_type k=0;k<xout_dim0;k++) {
762 d(k) = d(k) * (xout(k)-xeq(j))/(xeq(i)-xeq(j));
763 }
764 }
765 }
766 if (i!=0) {
767 for (ordinal_type k=0;k<xout_dim0;k++) {
768 d(k) = -d(k)/(xeq(i)-xeq(0));
769 }
770 }
771 if (i!=order) {
772 for (ordinal_type k=0;k<xout_dim0;k++) {
773 d(k) = d(k)/(xeq(i)-xeq(order));
774 }
775 }
776
777 for (ordinal_type k=0;k<xout_dim0;k++) {
778 warp(k) += d(k);
779 }
780 }
781 }
782
783
784template<typename pointValueType, class ...pointProperties>
785void
787getWarpBlendLatticeTetrahedron(Kokkos::DynRankView<pointValueType,pointProperties...> points,
788 const ordinal_type order ,
789 const ordinal_type offset )
790 {
791 pointValueType alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603,
792 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655};
793 pointValueType alpha;
794
795 if (order <= 15) {
796 alpha = alphastore[std::max(order-1,0)];
797 }
798 else {
799 alpha = 1.0;
800 }
801
802 const ordinal_type N = (order+1)*(order+2)*(order+3)/6;
803 pointValueType tol = 1.e-10;
804
805 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> shift("shift",N,3);
806 Kokkos::deep_copy(shift,0.0);
807
808 /* create 3d equidistributed nodes on Warburton tet */
809 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> equipoints("equipoints", N,3);
810 ordinal_type sk = 0;
811 for (ordinal_type n=0;n<=order;n++) {
812 for (ordinal_type m=0;m<=order-n;m++) {
813 for (ordinal_type q=0;q<=order-n-m;q++) {
814 equipoints(sk,0) = -1.0 + (q * 2.0 ) / order;
815 equipoints(sk,1) = -1.0 + (m * 2.0 ) / order;
816 equipoints(sk,2) = -1.0 + (n * 2.0 ) / order;
817 sk++;
818 }
819 }
820 }
821
822
823 /* construct barycentric coordinates for equispaced lattice */
824 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1("L1",N);
825 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2("L2",N);
826 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3("L3",N);
827 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L4("L4",N);
828 for (ordinal_type i=0;i<N;i++) {
829 L1(i) = (1.0 + equipoints(i,2)) / 2.0;
830 L2(i) = (1.0 + equipoints(i,1)) / 2.0;
831 L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0;
832 L4(i) = (1.0 + equipoints(i,0)) / 2.0;
833 }
834
835
836 /* vertices of Warburton tet */
837 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warVerts_("warVerts",1,4,3);
838 auto warVerts = Kokkos::subview(warVerts_,0,Kokkos::ALL(),Kokkos::ALL());
839 warVerts(0,0) = -1.0;
840 warVerts(0,1) = -1.0/sqrt(3.0);
841 warVerts(0,2) = -1.0/sqrt(6.0);
842 warVerts(1,0) = 1.0;
843 warVerts(1,1) = -1.0/sqrt(3.0);
844 warVerts(1,2) = -1.0/sqrt(6.0);
845 warVerts(2,0) = 0.0;
846 warVerts(2,1) = 2.0 / sqrt(3.0);
847 warVerts(2,2) = -1.0/sqrt(6.0);
848 warVerts(3,0) = 0.0;
849 warVerts(3,1) = 0.0;
850 warVerts(3,2) = 3.0 / sqrt(6.0);
851
852
853 /* tangents to faces */
854 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t1("t1",4,3);
855 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t2("t2",4,3);
856 for (ordinal_type i=0;i<3;i++) {
857 t1(0,i) = warVerts(1,i) - warVerts(0,i);
858 t1(1,i) = warVerts(1,i) - warVerts(0,i);
859 t1(2,i) = warVerts(2,i) - warVerts(1,i);
860 t1(3,i) = warVerts(2,i) - warVerts(0,i);
861 t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
862 t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
863 t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) );
864 t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) );
865 }
866
867 /* normalize tangents */
868 for (ordinal_type n=0;n<4;n++) {
869 /* Compute norm of t1(n) and t2(n) */
870 pointValueType normt1n = 0.0;
871 pointValueType normt2n = 0.0;
872 for (ordinal_type i=0;i<3;i++) {
873 normt1n += (t1(n,i) * t1(n,i));
874 normt2n += (t2(n,i) * t2(n,i));
875 }
876 normt1n = sqrt(normt1n);
877 normt2n = sqrt(normt2n);
878 /* normalize each tangent now */
879 for (ordinal_type i=0;i<3;i++) {
880 t1(n,i) /= normt1n;
881 t2(n,i) /= normt2n;
882 }
883 }
884
885 /* undeformed coordinates */
886 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> XYZ("XYZ",N,3);
887 for (ordinal_type i=0;i<N;i++) {
888 for (ordinal_type j=0;j<3;j++) {
889 XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j);
890 }
891 }
892
893 for (ordinal_type face=1;face<=4;face++) {
894 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> La, Lb, Lc, Ld;
895 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp("warp",N,2);
896 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend("blend",N);
897 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> denom("denom",N);
898 switch (face) {
899 case 1:
900 La = L1; Lb = L2; Lc = L3; Ld = L4; break;
901 case 2:
902 La = L2; Lb = L1; Lc = L3; Ld = L4; break;
903 case 3:
904 La = L3; Lb = L1; Lc = L4; Ld = L2; break;
905 case 4:
906 La = L4; Lb = L1; Lc = L3; Ld = L2; break;
907 }
908
909 /* get warp tangential to face */
910 warpShiftFace3D(warp,order,alpha,La,Lb,Lc,Ld);
911
912 for (ordinal_type k=0;k<N;k++) {
913 blend(k) = Lb(k) * Lc(k) * Ld(k);
914 }
915
916 for (ordinal_type k=0;k<N;k++) {
917 denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k));
918 }
919
920 for (ordinal_type k=0;k<N;k++) {
921 if (denom(k) > tol) {
922 blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k);
923 }
924 }
925
926
927 // compute warp and blend
928 for (ordinal_type k=0;k<N;k++) {
929 for (ordinal_type j=0;j<3;j++) {
930 shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j)
931 + blend(k) * warp(k,1) * t2(face-1,j);
932 }
933 }
934
935 for (ordinal_type k=0;k<N;k++) {
936 if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) {
937 for (ordinal_type j=0;j<3;j++) {
938 shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j);
939 }
940 }
941 }
942
943 }
944
945 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> updatedPoints("updatedPoints",1,N,3);
946 for (ordinal_type k=0;k<N;k++) {
947 for (ordinal_type j=0;j<3;j++) {
948 updatedPoints(0,k,j) = XYZ(k,j) + shift(k,j);
949 }
950 }
951
952 //warVerts.resize( 1 , 4 , 3 );
953
954 // now we convert to Pavel's reference triangle!
955 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> refPts("refPts",1,N,3);
957 warVerts_ ,
958 shards::getCellTopologyData<shards::Tetrahedron<4> >()
959 );
960
961 auto pointsHost = Kokkos::create_mirror_view(points);
962 // now write from refPts into points, taking offset into account
963 ordinal_type noffcur = 0;
964 ordinal_type offcur = 0;
965 for (ordinal_type i=0;i<=order;i++) {
966 for (ordinal_type j=0;j<=order-i;j++) {
967 for (ordinal_type k=0;k<=order-i-j;k++) {
968 if ( (i >= offset) && (i <= order-offset) &&
969 (j >= offset) && (j <= order-i-offset) &&
970 (k >= offset) && (k <= order-i-j-offset) ) {
971 pointsHost(offcur,0) = refPts(0,noffcur,0);
972 pointsHost(offcur,1) = refPts(0,noffcur,1);
973 pointsHost(offcur,2) = refPts(0,noffcur,2);
974 offcur++;
975 }
976 noffcur++;
977 }
978 }
979 }
980
981 Kokkos::deep_copy(points, pointsHost);
982 }
983
984
985} // face Intrepid
986#endif
Kokkos::DynRankView< typename ViewType::value_type, typename DeduceLayout< ViewType >::result_layout, typename ViewType::device_type > getMatchingViewWithLabel(const ViewType &view, const std::string &label, DimArgs... dims)
Creates and returns a view that matches the provided view in Kokkos Layout.
static void mapToReferenceFrame(Kokkos::DynRankView< refPointValueType, refPointProperties... > refPoints, const Kokkos::DynRankView< physPointValueType, physPointProperties... > physPoints, const Kokkos::DynRankView< worksetCellValueType, worksetCellProperties... > worksetCell, const shards::CellTopology cellTopo)
Computes , the inverse of the reference-to-physical frame map using a default initial guess.
static void getWarpBlendLatticeLine(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0)
Returns the Gauss-Lobatto points of a given order on the reference line [-1,1]. The output array is (...
static void getEquispacedLatticeTriangle(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0)
Computes an equispaced lattice of a given order on the reference triangle. The output array is (P,...
static void evalshift(Kokkos::DynRankView< pointValueType, pointProperties... > dxy, const ordinal_type order, const pointValueType pval, const Kokkos::DynRankView< pointValueType, pointProperties... > L1, const Kokkos::DynRankView< pointValueType, pointProperties... > L2, const Kokkos::DynRankView< pointValueType, pointProperties... > L3)
Used internally to evaluate the point shift for warp-blend points on faces of tets.
static void getLatticePyramid(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference pyramid. The output array is (P,...
static void getLatticeTetrahedron(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference tetrahedron. The output array is (P,...
static void warpFactor(Kokkos::DynRankView< pointValueType, pointProperties... > warp, const ordinal_type order, const Kokkos::DynRankView< pointValueType, pointProperties... > xnodes, const Kokkos::DynRankView< pointValueType, pointProperties... > xout)
interpolates Warburton's warp function on the line
static ordinal_type getLatticeSize(const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...
static void getWarpBlendLatticeTetrahedron(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0)
Returns Warburton's warp-blend points of a given order on the reference tetrahedron....
static void getLatticeTriangle(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference triangle. The output array is (P,...
static void evalwarp(Kokkos::DynRankView< pointValueType, pointProperties... > warp, const ordinal_type order, const Kokkos::DynRankView< pointValueType, pointProperties... > xnodes, const Kokkos::DynRankView< pointValueType, pointProperties... > xout)
Used internally to compute the warp on edges of a triangle in warp-blend points.
static void getGaussPoints(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order)
static void getLatticeLine(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference line. The output array is (P,...
static void getWarpBlendLatticeTriangle(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0)
Returns Warburton's warp-blend points of a given order on the reference triangle. The output array is...
static void getEquispacedLatticeTetrahedron(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0)
Computes an equispaced lattice of a given order on the reference tetrahedron. The output array is (P,...
static void getEquispacedLatticePyramid(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0)
Computes an equispaced lattice of a given order on the reference pyramid. The output array is (P,...
static void warpShiftFace3D(Kokkos::DynRankView< pointValueType, pointProperties... > dxy, const ordinal_type order, const pointValueType pval, const Kokkos::DynRankView< pointValueType, pointProperties... > L1, const Kokkos::DynRankView< pointValueType, pointProperties... > L2, const Kokkos::DynRankView< pointValueType, pointProperties... > L3, const Kokkos::DynRankView< pointValueType, pointProperties... > L4)
This is used internally to compute the tetrahedral warp-blend points one each face.
static void getLattice(Kokkos::DynRankView< pointValueType, pointProperties... > points, const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference simplex, quadrilateral or hexahedron (cu...
static void getEquispacedLatticeLine(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order, const ordinal_type offset=0)
Computes an equispaced lattice of a given order on the reference line [-1,1]. The output array is (P,...