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 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
111template<typename pointValueType, class ...pointProperties>
112void
114getLattice( 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
182template<typename pointValueType, class ...pointProperties>
183void
185getGaussPoints( 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
219template<typename pointValueType, class ...pointProperties>
220void
222getLatticeLine( 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
237template<typename pointValueType, class ...pointProperties>
238void
240getLatticeTriangle( 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
255template<typename pointValueType, class ...pointProperties>
256void
258getLatticeTetrahedron( 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
275template<typename pointValueType, class ...pointProperties>
276void
278getEquispacedLatticeLine( 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
295template<typename pointValueType, class ...pointProperties>
296void
298getWarpBlendLatticeLine( 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 const auto common_range = range_type(0, std::min(pts.extent(0), z.extent(0)));
325 Kokkos::deep_copy(Kokkos::subview(pts, common_range),
326 Kokkos::subview(z, common_range));
327 }
328}
329
330template<typename pointValueType, class ...pointProperties>
331void
333getEquispacedLatticeTriangle( Kokkos::DynRankView<pointValueType,pointProperties...> points,
334 const ordinal_type order,
335 const ordinal_type offset ) {
336 TEUCHOS_TEST_FOR_EXCEPTION( order <= 0 ,
337 std::invalid_argument ,
338 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTriangle): order must be positive" );
339
340 auto pointsHost = Kokkos::create_mirror_view(points);
341
342 const pointValueType h = 1.0 / order;
343 ordinal_type cur = 0;
344
345 for (ordinal_type i=offset;i<=order-offset;i++) {
346 for (ordinal_type j=offset;j<=order-i-offset;j++) {
347 pointsHost(cur,0) = j * h ;
348 pointsHost(cur,1) = i * h;
349 cur++;
350 }
351 }
352 Kokkos::deep_copy(points, pointsHost);
353}
354
355template<typename pointValueType, class ...pointProperties>
356void
358getEquispacedLatticeTetrahedron( Kokkos::DynRankView<pointValueType,pointProperties...> points,
359 const ordinal_type order ,
360 const ordinal_type offset ) {
361 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
362 std::invalid_argument ,
363 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTetrahedron): order must be positive" );
364
365 auto pointsHost = Kokkos::create_mirror_view(points);
366
367 const pointValueType h = 1.0 / order;
368 ordinal_type cur = 0;
369
370 for (ordinal_type i=offset;i<=order-offset;i++) {
371 for (ordinal_type j=offset;j<=order-i-offset;j++) {
372 for (ordinal_type k=offset;k<=order-i-j-offset;k++) {
373 pointsHost(cur,0) = k * h;
374 pointsHost(cur,1) = j * h;
375 pointsHost(cur,2) = i * h;
376 cur++;
377 }
378 }
379 }
380 Kokkos::deep_copy(points, pointsHost);
381}
382
383
384template<typename pointValueType, class ...pointProperties>
385void
387warpFactor( Kokkos::DynRankView<pointValueType,pointProperties...> warp,
388 const ordinal_type order,
389 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes ,
390 const Kokkos::DynRankView<pointValueType,pointProperties...> xout
391 )
392 {
393 TEUCHOS_TEST_FOR_EXCEPTION( ( warp.extent(0) != xout.extent(0) ) ,
394 std::invalid_argument ,
395 ">>> ERROR (PointTools::warpFactor): xout and warp must be same size." );
396
397 Kokkos::deep_copy(warp, pointValueType(0.0));
398
399 ordinal_type xout_dim0 = xout.extent(0);
400 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> d("d", xout_dim0 );
401
402 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> xeq_("xeq", order + 1 ,1);
403 PointTools::getEquispacedLatticeLine( xeq_ , order , 0 );
404 const auto xeq = Kokkos::subview(xeq_, Kokkos::ALL(),0);
405
406 TEUCHOS_TEST_FOR_EXCEPTION( ( xeq.extent(0) != xnodes.extent(0) ) ,
407 std::invalid_argument ,
408 ">>> ERROR (PointTools::warpFactor): xeq and xnodes must be same size." );
409
410
411 for (ordinal_type i=0;i<=order;i++) {
412
413 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
414
415 for (ordinal_type j=1;j<order;j++) {
416 if (i != j) {
417 for (ordinal_type k=0;k<xout_dim0;k++) {
418 d(k) = d(k) * ( (xout(k)-xeq(j)) / (xeq(i)-xeq(j)) );
419 }
420 }
421 }
422
423 // deflate end roots
424 if ( i != 0 ) {
425 for (ordinal_type k=0;k<xout_dim0;k++) {
426 d(k) = -d(k) / (xeq(i) - xeq(0));
427 }
428 }
429
430 if (i != order ) {
431 for (ordinal_type k=0;k<xout_dim0;k++) {
432 d(k) = d(k) / (xeq(i) - xeq(order));
433 }
434 }
435
436 for (ordinal_type k=0;k<xout_dim0;k++) {
437 warp(k) += d(k);
438 }
439
440 }
441 }
442
443template<typename pointValueType, class ...pointProperties>
444void
446getWarpBlendLatticeTriangle( Kokkos::DynRankView<pointValueType,pointProperties...> points,
447 const ordinal_type order ,
448 const ordinal_type offset)
449 {
450 /* get Gauss-Lobatto points */
451
452 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> gaussX("gaussX", order + 1 , 1 );
453
454 PointTools::getWarpBlendLatticeLine( gaussX , order , 0 );
455 //auto gaussX = Kokkos::subdynrankview(gaussX_, Kokkos::ALL(),0);
456
457 // gaussX.resize(gaussX.extent(0));
458
459 pointValueType alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999,
460 1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
461
462 pointValueType alpha;
463
464 if (order >= 1 && order < 16) {
465 alpha = alpopt[order-1];
466 }
467 else {
468 alpha = 5.0 / 3.0;
469 }
470
471 const ordinal_type p = order; /* switch to Warburton's notation */
472 ordinal_type N = (p+1)*(p+2)/2;
473
474 /* equidistributed nodes on equilateral triangle */
475 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1("L1", N );
476 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2("L2", N );
477 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3("L3", N );
478 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> X("X", N);
479 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> Y("Y", N);
480
481 ordinal_type sk = 0;
482 for (ordinal_type n=1;n<=p+1;n++) {
483 for (ordinal_type m=1;m<=p+2-n;m++) {
484 L1(sk) = (n-1) / (pointValueType)p;
485 L3(sk) = (m-1) / (pointValueType)p;
486 L2(sk) = 1.0 - L1(sk) - L3(sk);
487 sk++;
488 }
489 }
490
491 for (ordinal_type n=0;n<N;n++) {
492 X(n) = -L2(n) + L3(n);
493 Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772;
494 }
495
496 /* get blending function for each node at each edge */
497 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend1("blend1", N);
498 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend2("blend2", N);
499 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend3("blend3", N);
500
501 for (ordinal_type n=0;n<N;n++) {
502 blend1(n) = 4.0 * L2(n) * L3(n);
503 blend2(n) = 4.0 * L1(n) * L3(n);
504 blend3(n) = 4.0 * L1(n) * L2(n);
505 }
506
507 /* get difference of each barycentric coordinate */
508 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3mL2("L3mL2", N);
509 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1mL3("L1mL3", N);
510 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2mL1("L2mL1", N);
511
512 for (ordinal_type k=0;k<N;k++) {
513 L3mL2(k) = L3(k)-L2(k);
514 L1mL3(k) = L1(k)-L3(k);
515 L2mL1(k) = L2(k)-L1(k);
516 }
517
518 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor1("warpfactor1", N);
519 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor2("warpfactor2", N);
520 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor3("warpfactor3", N);
521
522 warpFactor( warpfactor1, order , gaussX , L3mL2 );
523 warpFactor( warpfactor2, order , gaussX , L1mL3 );
524 warpFactor( warpfactor3, order , gaussX , L2mL1 );
525
526 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp1("warp1", N);
527 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp2("warp2", N);
528 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp3("warp3", N);
529
530 for (ordinal_type k=0;k<N;k++) {
531 warp1(k) = blend1(k) * warpfactor1(k) *
532 ( 1.0 + alpha * alpha * L1(k) * L1(k) );
533 warp2(k) = blend2(k) * warpfactor2(k) *
534 ( 1.0 + alpha * alpha * L2(k) * L2(k) );
535 warp3(k) = blend3(k) * warpfactor3(k) *
536 ( 1.0 + alpha * alpha * L3(k) * L3(k) );
537 }
538
539 for (ordinal_type k=0;k<N;k++) {
540 X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k);
541 Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k);
542 }
543
544 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warXY("warXY", 1, N,2);
545
546 for (ordinal_type k=0;k<N;k++) {
547 warXY(0, k,0) = X(k);
548 warXY(0, k,1) = Y(k);
549 }
550
551
552 // finally, convert the warp-blend points to the correct triangle
553 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warburtonVerts("warburtonVerts", 1,3,2);
554 warburtonVerts(0,0,0) = -1.0;
555 warburtonVerts(0,0,1) = -1.0/std::sqrt(3.0);
556 warburtonVerts(0,1,0) = 1.0;
557 warburtonVerts(0,1,1) = -1.0/std::sqrt(3.0);
558 warburtonVerts(0,2,0) = 0.0;
559 warburtonVerts(0,2,1) = 2.0/std::sqrt(3.0);
560
561 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> refPts("refPts", 1, N,2);
562
564 warXY ,
565 warburtonVerts ,
566 shards::getCellTopologyData< shards::Triangle<3> >()
567 );
568
569
570 auto pointsHost = Kokkos::create_mirror_view(points);
571 // now write from refPts into points, taking care of offset
572 ordinal_type noffcur = 0; // index into refPts
573 ordinal_type offcur = 0; // index ordinal_type points
574 for (ordinal_type i=0;i<=order;i++) {
575 for (ordinal_type j=0;j<=order-i;j++) {
576 if ( (i >= offset) && (i <= order-offset) &&
577 (j >= offset) && (j <= order-i-offset) ) {
578 pointsHost(offcur,0) = refPts(0, noffcur,0);
579 pointsHost(offcur,1) = refPts(0, noffcur,1);
580 offcur++;
581 }
582 noffcur++;
583 }
584 }
585
586 Kokkos::deep_copy(points, pointsHost);
587
588 }
589
590
591template<typename pointValueType, class ...pointProperties>
592void
594warpShiftFace3D( Kokkos::DynRankView<pointValueType,pointProperties...> dxy,
595 const ordinal_type order ,
596 const pointValueType pval ,
597 const Kokkos::DynRankView<pointValueType,pointProperties...> /* L1 */,
598 const Kokkos::DynRankView<pointValueType,pointProperties...> L2,
599 const Kokkos::DynRankView<pointValueType,pointProperties...> L3,
600 const Kokkos::DynRankView<pointValueType,pointProperties...> L4
601 )
602 {
603 evalshift(dxy,order,pval,L2,L3,L4);
604 return;
605 }
606
607template<typename pointValueType, class ...pointProperties>
608void
610evalshift( Kokkos::DynRankView<pointValueType,pointProperties...> dxy,
611 const ordinal_type order ,
612 const pointValueType pval ,
613 const Kokkos::DynRankView<pointValueType,pointProperties...> L1 ,
614 const Kokkos::DynRankView<pointValueType,pointProperties...> L2 ,
615 const Kokkos::DynRankView<pointValueType,pointProperties...> L3
616 )
617 {
618 // get Gauss-Lobatto-nodes
619 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> gaussX("gaussX",order+1,1);
620 PointTools::getWarpBlendLatticeLine( gaussX , order , 0 );
621 //gaussX.resize(order+1);
622 const ordinal_type N = L1.extent(0);
623
624 // Warburton code reverses them
625 for (ordinal_type k=0;k<=order;k++) {
626 gaussX(k,0) *= -1.0;
627 }
628
629 // blending function at each node for each edge
630 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend1("blend1",N);
631 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend2("blend2",N);
632 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend3("blend3",N);
633
634 for (ordinal_type i=0;i<N;i++) {
635 blend1(i) = L2(i) * L3(i);
636 blend2(i) = L1(i) * L3(i);
637 blend3(i) = L1(i) * L2(i);
638 }
639
640 // amount of warp for each node for each edge
641 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor1("warpfactor1",N);
642 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor2("warpfactor2",N);
643 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor3("warpfactor3",N);
644
645 // differences of barycentric coordinates
646 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3mL2("L3mL2",N);
647 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1mL3("L1mL3",N);
648 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2mL1("L2mL1",N);
649
650 for (ordinal_type k=0;k<N;k++) {
651 L3mL2(k) = L3(k)-L2(k);
652 L1mL3(k) = L1(k)-L3(k);
653 L2mL1(k) = L2(k)-L1(k);
654 }
655
656 evalwarp( warpfactor1 , order , gaussX , L3mL2 );
657 evalwarp( warpfactor2 , order , gaussX , L1mL3 );
658 evalwarp( warpfactor3 , order , gaussX , L2mL1 );
659
660 for (ordinal_type k=0;k<N;k++) {
661 warpfactor1(k) *= 4.0;
662 warpfactor2(k) *= 4.0;
663 warpfactor3(k) *= 4.0;
664 }
665
666 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp1("warp1",N);
667 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp2("warp2",N);
668 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp3("warp3",N);
669
670 for (ordinal_type k=0;k<N;k++) {
671 warp1(k) = blend1(k) * warpfactor1(k) *
672 ( 1.0 + pval * pval * L1(k) * L1(k) );
673 warp2(k) = blend2(k) * warpfactor2(k) *
674 ( 1.0 + pval * pval * L2(k) * L2(k) );
675 warp3(k) = blend3(k) * warpfactor3(k) *
676 ( 1.0 + pval * pval * L3(k) * L3(k) );
677 }
678
679 for (ordinal_type k=0;k<N;k++) {
680 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);
681 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);
682 }
683 }
684
685 /* one-d edge warping function */
686template<typename pointValueType, class ...pointProperties>
687void
689evalwarp(Kokkos::DynRankView<pointValueType,pointProperties...> warp ,
690 const ordinal_type order ,
691 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes,
692 const Kokkos::DynRankView<pointValueType,pointProperties...> xout )
693 {
694 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> xeq("xeq",order+1);
695
696 ordinal_type xout_dim0 = xout.extent(0);
697 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> d("d",xout_dim0);
698
699 //Kokkos::deep_copy(d, 0.0);
700
701 for (ordinal_type i=0;i<=order;i++) {
702 xeq(i) = -1.0 + 2.0 * ( order - i ) / order;
703 }
704
705 for (ordinal_type i=0;i<=order;i++) {
706 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
707 for (ordinal_type j=1;j<order;j++) {
708 if (i!=j) {
709 for (ordinal_type k=0;k<xout_dim0;k++) {
710 d(k) = d(k) * (xout(k)-xeq(j))/(xeq(i)-xeq(j));
711 }
712 }
713 }
714 if (i!=0) {
715 for (ordinal_type k=0;k<xout_dim0;k++) {
716 d(k) = -d(k)/(xeq(i)-xeq(0));
717 }
718 }
719 if (i!=order) {
720 for (ordinal_type k=0;k<xout_dim0;k++) {
721 d(k) = d(k)/(xeq(i)-xeq(order));
722 }
723 }
724
725 for (ordinal_type k=0;k<xout_dim0;k++) {
726 warp(k) += d(k);
727 }
728 }
729 }
730
731
732template<typename pointValueType, class ...pointProperties>
733void
735getWarpBlendLatticeTetrahedron(Kokkos::DynRankView<pointValueType,pointProperties...> points,
736 const ordinal_type order ,
737 const ordinal_type offset )
738 {
739 pointValueType alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603,
740 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655};
741 pointValueType alpha;
742
743 if (order <= 15) {
744 alpha = alphastore[std::max(order-1,0)];
745 }
746 else {
747 alpha = 1.0;
748 }
749
750 const ordinal_type N = (order+1)*(order+2)*(order+3)/6;
751 pointValueType tol = 1.e-10;
752
753 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> shift("shift",N,3);
754 Kokkos::deep_copy(shift,0.0);
755
756 /* create 3d equidistributed nodes on Warburton tet */
757 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> equipoints("equipoints", N,3);
758 ordinal_type sk = 0;
759 for (ordinal_type n=0;n<=order;n++) {
760 for (ordinal_type m=0;m<=order-n;m++) {
761 for (ordinal_type q=0;q<=order-n-m;q++) {
762 equipoints(sk,0) = -1.0 + (q * 2.0 ) / order;
763 equipoints(sk,1) = -1.0 + (m * 2.0 ) / order;
764 equipoints(sk,2) = -1.0 + (n * 2.0 ) / order;
765 sk++;
766 }
767 }
768 }
769
770
771 /* construct barycentric coordinates for equispaced lattice */
772 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1("L1",N);
773 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2("L2",N);
774 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3("L3",N);
775 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L4("L4",N);
776 for (ordinal_type i=0;i<N;i++) {
777 L1(i) = (1.0 + equipoints(i,2)) / 2.0;
778 L2(i) = (1.0 + equipoints(i,1)) / 2.0;
779 L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0;
780 L4(i) = (1.0 + equipoints(i,0)) / 2.0;
781 }
782
783
784 /* vertices of Warburton tet */
785 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warVerts_("warVerts",1,4,3);
786 auto warVerts = Kokkos::subview(warVerts_,0,Kokkos::ALL(),Kokkos::ALL());
787 warVerts(0,0) = -1.0;
788 warVerts(0,1) = -1.0/sqrt(3.0);
789 warVerts(0,2) = -1.0/sqrt(6.0);
790 warVerts(1,0) = 1.0;
791 warVerts(1,1) = -1.0/sqrt(3.0);
792 warVerts(1,2) = -1.0/sqrt(6.0);
793 warVerts(2,0) = 0.0;
794 warVerts(2,1) = 2.0 / sqrt(3.0);
795 warVerts(2,2) = -1.0/sqrt(6.0);
796 warVerts(3,0) = 0.0;
797 warVerts(3,1) = 0.0;
798 warVerts(3,2) = 3.0 / sqrt(6.0);
799
800
801 /* tangents to faces */
802 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t1("t1",4,3);
803 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t2("t2",4,3);
804 for (ordinal_type i=0;i<3;i++) {
805 t1(0,i) = warVerts(1,i) - warVerts(0,i);
806 t1(1,i) = warVerts(1,i) - warVerts(0,i);
807 t1(2,i) = warVerts(2,i) - warVerts(1,i);
808 t1(3,i) = warVerts(2,i) - warVerts(0,i);
809 t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
810 t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
811 t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) );
812 t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) );
813 }
814
815 /* normalize tangents */
816 for (ordinal_type n=0;n<4;n++) {
817 /* Compute norm of t1(n) and t2(n) */
818 pointValueType normt1n = 0.0;
819 pointValueType normt2n = 0.0;
820 for (ordinal_type i=0;i<3;i++) {
821 normt1n += (t1(n,i) * t1(n,i));
822 normt2n += (t2(n,i) * t2(n,i));
823 }
824 normt1n = sqrt(normt1n);
825 normt2n = sqrt(normt2n);
826 /* normalize each tangent now */
827 for (ordinal_type i=0;i<3;i++) {
828 t1(n,i) /= normt1n;
829 t2(n,i) /= normt2n;
830 }
831 }
832
833 /* undeformed coordinates */
834 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> XYZ("XYZ",N,3);
835 for (ordinal_type i=0;i<N;i++) {
836 for (ordinal_type j=0;j<3;j++) {
837 XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j);
838 }
839 }
840
841 for (ordinal_type face=1;face<=4;face++) {
842 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> La, Lb, Lc, Ld;
843 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp("warp",N,2);
844 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend("blend",N);
845 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> denom("denom",N);
846 switch (face) {
847 case 1:
848 La = L1; Lb = L2; Lc = L3; Ld = L4; break;
849 case 2:
850 La = L2; Lb = L1; Lc = L3; Ld = L4; break;
851 case 3:
852 La = L3; Lb = L1; Lc = L4; Ld = L2; break;
853 case 4:
854 La = L4; Lb = L1; Lc = L3; Ld = L2; break;
855 }
856
857 /* get warp tangential to face */
858 warpShiftFace3D(warp,order,alpha,La,Lb,Lc,Ld);
859
860 for (ordinal_type k=0;k<N;k++) {
861 blend(k) = Lb(k) * Lc(k) * Ld(k);
862 }
863
864 for (ordinal_type k=0;k<N;k++) {
865 denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k));
866 }
867
868 for (ordinal_type k=0;k<N;k++) {
869 if (denom(k) > tol) {
870 blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k);
871 }
872 }
873
874
875 // compute warp and blend
876 for (ordinal_type k=0;k<N;k++) {
877 for (ordinal_type j=0;j<3;j++) {
878 shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j)
879 + blend(k) * warp(k,1) * t2(face-1,j);
880 }
881 }
882
883 for (ordinal_type k=0;k<N;k++) {
884 if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) {
885 for (ordinal_type j=0;j<3;j++) {
886 shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j);
887 }
888 }
889 }
890
891 }
892
893 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> updatedPoints("updatedPoints",1,N,3);
894 for (ordinal_type k=0;k<N;k++) {
895 for (ordinal_type j=0;j<3;j++) {
896 updatedPoints(0,k,j) = XYZ(k,j) + shift(k,j);
897 }
898 }
899
900 //warVerts.resize( 1 , 4 , 3 );
901
902 // now we convert to Pavel's reference triangle!
903 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> refPts("refPts",1,N,3);
905 warVerts_ ,
906 shards::getCellTopologyData<shards::Tetrahedron<4> >()
907 );
908
909 auto pointsHost = Kokkos::create_mirror_view(points);
910 // now write from refPts into points, taking offset into account
911 ordinal_type noffcur = 0;
912 ordinal_type offcur = 0;
913 for (ordinal_type i=0;i<=order;i++) {
914 for (ordinal_type j=0;j<=order-i;j++) {
915 for (ordinal_type k=0;k<=order-i-j;k++) {
916 if ( (i >= offset) && (i <= order-offset) &&
917 (j >= offset) && (j <= order-i-offset) &&
918 (k >= offset) && (k <= order-i-j-offset) ) {
919 pointsHost(offcur,0) = refPts(0,noffcur,0);
920 pointsHost(offcur,1) = refPts(0,noffcur,1);
921 pointsHost(offcur,2) = refPts(0,noffcur,2);
922 offcur++;
923 }
924 noffcur++;
925 }
926 }
927 }
928
929 Kokkos::deep_copy(points, pointsHost);
930 }
931
932
933} // face Intrepid
934#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 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 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,...
Gauss-Jacobi/Gauss-Radau-Jacobi/Gauss-Lobatto zeros and weights.