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