Intrepid
Intrepid_PointToolsDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid 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 Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49namespace Intrepid {
50
51 template<class Scalar, class ArrayType>
52 void PointTools::getLattice(ArrayType &pts ,
53 const shards::CellTopology& cellType ,
54 const int order ,
55 const int offset ,
56 const EPointType pointType )
57 {
58 switch (pointType) {
59 case POINTTYPE_EQUISPACED:
60 getEquispacedLattice<Scalar,ArrayType>( cellType , pts , order , offset );
61 break;
62 case POINTTYPE_WARPBLEND:
63
64 getWarpBlendLattice<Scalar,ArrayType>( cellType , pts , order , offset );
65 break;
66 default:
67 TEUCHOS_TEST_FOR_EXCEPTION( true ,
68 std::invalid_argument ,
69 "PointTools::getLattice: invalid EPointType" );
70 }
71 return;
72 }
73
74 template<class Scalar, class ArrayType>
75 void PointTools::getGaussPoints( ArrayType &pts ,
76 const int order )
77 {
78
79 Scalar *z = new Scalar[order+1];
80 Scalar *w = new Scalar[order+1];
81
82 IntrepidPolylib::zwgj( z , w , order + 1 , 0.0 , 0.0 );
83 for (int i=0;i<order+1;i++) {
84 pts(i,0) = z[i];
85 }
86
87 delete []z;
88 delete []w;
89 }
90
91
92 template<class Scalar, class ArrayType>
93 void PointTools::getEquispacedLattice(const shards::CellTopology& cellType ,
94 ArrayType &points ,
95 const int order ,
96 const int offset )
97
98 {
99 switch (cellType.getKey()) {
100 case shards::Tetrahedron<4>::key:
101 case shards::Tetrahedron<8>::key:
102 case shards::Tetrahedron<10>::key:
103 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) )
104 || points.dimension(1) != 3 ,
105 std::invalid_argument ,
106 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
107 getEquispacedLatticeTetrahedron<Scalar,ArrayType>( points , order , offset );
108 break;
109 case shards::Triangle<3>::key:
110 case shards::Triangle<4>::key:
111 case shards::Triangle<6>::key:
112 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) )
113 || points.dimension(1) != 2 ,
114 std::invalid_argument ,
115 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
116 getEquispacedLatticeTriangle<Scalar,ArrayType>( points , order , offset );
117 break;
118 case shards::Line<2>::key:
119 case shards::Line<3>::key:
120 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) )
121 || points.dimension(1) != 1 ,
122 std::invalid_argument ,
123 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
124 getEquispacedLatticeLine<Scalar,ArrayType>( points , order , offset );
125 break;
126 default:
127 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
128 ">>> ERROR (Intrepid::PointTools::getEquispacedLattice): Illegal cell type" );
129 }
130
131 }
132
133 template<class Scalar, class ArrayType>
134 void PointTools::getWarpBlendLattice( const shards::CellTopology& cellType ,
135 ArrayType &points ,
136 const int order ,
137 const int offset )
138
139 {
140
141 switch (cellType.getKey()) {
142 case shards::Tetrahedron<4>::key:
143 case shards::Tetrahedron<8>::key:
144 case shards::Tetrahedron<10>::key:
145 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) )
146 || points.dimension(1) != 3 ,
147 std::invalid_argument ,
148 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
149
150 getWarpBlendLatticeTetrahedron<Scalar,ArrayType>( points , order , offset );
151 break;
152 case shards::Triangle<3>::key:
153 case shards::Triangle<4>::key:
154 case shards::Triangle<6>::key:
155 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) )
156 || points.dimension(1) != 2 ,
157 std::invalid_argument ,
158 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
159
160 getWarpBlendLatticeTriangle<Scalar,ArrayType>( points , order , offset );
161 break;
162 case shards::Line<2>::key:
163 case shards::Line<3>::key:
164 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) )
165 || points.dimension(1) != 1 ,
166 std::invalid_argument ,
167 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
168
169 getWarpBlendLatticeLine<Scalar,ArrayType>( points , order , offset );
170 break;
171 default:
172 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
173 ">>> ERROR (Intrepid::PointTools::getWarpBlendLattice): Illegal cell type" );
174 }
175
176 }
177
178 template<class Scalar, class ArrayType>
179 void PointTools::getEquispacedLatticeLine( ArrayType &points ,
180 const int order ,
181 const int offset )
182 {
183 TEUCHOS_TEST_FOR_EXCEPTION( order < 0 ,
184 std::invalid_argument ,
185 ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeLine): order must be positive" );
186 if (order == 0) {
187 points(0,0) = 0.0;
188 }
189 else {
190 const Scalar h = 2.0 / order;
191
192 for (int i=offset;i<=order-offset;i++) {
193 points(i-offset,0) = -1.0 + h * (Scalar) i;
194 }
195 }
196
197 return;
198 }
199
200 template<class Scalar, class ArrayType>
202 const int order ,
203 const int offset )
204 {
205 TEUCHOS_TEST_FOR_EXCEPTION( order <= 0 ,
206 std::invalid_argument ,
207 ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeLine): order must be positive" );
208
209 const Scalar h = 1.0 / order;
210 int cur = 0;
211
212 for (int i=offset;i<=order-offset;i++) {
213 for (int j=offset;j<=order-i-offset;j++) {
214 points(cur,0) = (Scalar)0.0 + (Scalar) j * h ;
215 points(cur,1) = (Scalar)0.0 + (Scalar) i * h;
216 cur++;
217 }
218 }
219
220 return;
221 }
222
223 template<class Scalar, class ArrayType>
225 const int order ,
226 const int offset )
227 {
228 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
229 std::invalid_argument ,
230 ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeTetrahedron): order must be positive" );
231
232 const Scalar h = 1.0 / order;
233 int cur = 0;
234
235 for (int i=offset;i<=order-offset;i++) {
236 for (int j=offset;j<=order-i-offset;j++) {
237 for (int k=offset;k<=order-i-j-offset;k++) {
238 points(cur,0) = (Scalar) k * h;
239 points(cur,1) = (Scalar) j * h;
240 points(cur,2) = (Scalar) i * h;
241 cur++;
242 }
243 }
244 }
245
246 return;
247 }
248
249 template<class Scalar, class ArrayType>
250 void PointTools::getWarpBlendLatticeLine( ArrayType &points ,
251 const int order ,
252 const int offset )
253 {
254 Scalar *z = new Scalar[order+1];
255 Scalar *w = new Scalar[order+1];
256
257 // order is order of polynomial degree. The Gauss-Lobatto points are accurate
258 // up to degree 2*i-1
259 IntrepidPolylib::zwglj( z , w , order+1 , 0.0 , 0.0 );
260
261 for (int i=offset;i<=order-offset;i++) {
262 points(i-offset,0) = z[i];
263 }
264
265 delete []z;
266 delete []w;
267
268 return;
269 }
270
271 template<class Scalar, class ArrayType>
272 void PointTools::warpFactor( const int order ,
273 const ArrayType &xnodes ,
274 const ArrayType &xout ,
275 ArrayType &warp)
276 {
277 TEUCHOS_TEST_FOR_EXCEPTION( ( warp.dimension(0) != xout.dimension(0) ) ,
278 std::invalid_argument ,
279 ">>> ERROR (PointTools::warpFactor): xout and warp must be same size." );
280
281 warp.initialize();
282
283 FieldContainer<Scalar> d( xout.dimension(0) );
284 d.initialize();
285
286 FieldContainer<Scalar> xeq( order + 1 ,1);
287 PointTools::getEquispacedLatticeLine<Scalar,ArrayType>( xeq , order , 0 );
288 xeq.resize( order + 1 );
289
290 TEUCHOS_TEST_FOR_EXCEPTION( ( xeq.dimension(0) != xnodes.dimension(0) ) ,
291 std::invalid_argument ,
292 ">>> ERROR (PointTools::warpFactor): xeq and xnodes must be same size." );
293
294 for (int i=0;i<=order;i++) {
295
296 for (int k=0;k<d.dimension(0);k++) {
297 d(k) = xnodes(i) - xeq(i);
298 }
299
300 for (int j=1;j<order;j++) {
301 if (i != j) {
302 for (int k=0;k<d.dimension(0);k++) {
303 d(k) = d(k) * ( (xout(k)-xeq(j)) / (xeq(i)-xeq(j)) );
304 }
305 }
306 }
307
308 // deflate end roots
309 if ( i != 0 ) {
310 for (int k=0;k<d.dimension(0);k++) {
311 d(k) = -d(k) / (xeq(i) - xeq(0));
312 }
313 }
314
315 if (i != order ) {
316 for (int k=0;k<d.dimension(0);k++) {
317 d(k) = d(k) / (xeq(i) - xeq(order));
318 }
319 }
320
321 for (int k=0;k<d.dimension(0);k++) {
322 warp(k) += d(k);
323 }
324
325 }
326
327
328 return;
329 }
330
331 template<class Scalar, class ArrayType>
333 const int order ,
334 const int offset )
335 {
336 /* get Gauss-Lobatto points */
337
338 Intrepid::FieldContainer<Scalar> gaussX( order + 1 , 1 );
339
340 PointTools::getWarpBlendLatticeLine<Scalar,FieldContainer<Scalar> >( gaussX , order , 0 );
341
342 gaussX.resize(gaussX.dimension(0));
343
344 Scalar alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999,
345 1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
346
347 Scalar alpha;
348
349 if (order >= 1 && order < 16) {
350 alpha = alpopt[order-1];
351 }
352 else {
353 alpha = 5.0 / 3.0;
354 }
355
356 const int p = order; /* switch to Warburton's notation */
357 int N = (p+1)*(p+2)/2;
358
359 /* equidistributed nodes on equilateral triangle */
365
366 int sk = 0;
367 for (int n=1;n<=p+1;n++) {
368 for (int m=1;m<=p+2-n;m++) {
369 L1(sk) = (n-1) / (Scalar)p;
370 L3(sk) = (m-1) / (Scalar)p;
371 L2(sk) = 1.0 - L1(sk) - L3(sk);
372 sk++;
373 }
374 }
375
376 for (int n=0;n<N;n++) {
377 X(n) = -L2(n) + L3(n);
378 Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772;
379 }
380
381 /* get blending function for each node at each edge */
385
386 for (int n=0;n<N;n++) {
387 blend1(n) = 4.0 * L2(n) * L3(n);
388 blend2(n) = 4.0 * L1(n) * L3(n);
389 blend3(n) = 4.0 * L1(n) * L2(n);
390 }
391
392 /* get difference of each barycentric coordinate */
396
397 for (int k=0;k<N;k++) {
398 L3mL2(k) = L3(k)-L2(k);
399 L1mL3(k) = L1(k)-L3(k);
400 L2mL1(k) = L2(k)-L1(k);
401 }
402
403 FieldContainer<Scalar> warpfactor1(N);
404 FieldContainer<Scalar> warpfactor2(N);
405 FieldContainer<Scalar> warpfactor3(N);
406
407 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L3mL2 , warpfactor1 );
408 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L1mL3 , warpfactor2 );
409 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L2mL1 , warpfactor3 );
410
411 FieldContainer<Scalar> warp1(N);
412 FieldContainer<Scalar> warp2(N);
413 FieldContainer<Scalar> warp3(N);
414
415 for (int k=0;k<N;k++) {
416 warp1(k) = blend1(k) * warpfactor1(k) *
417 ( 1.0 + alpha * alpha * L1(k) * L1(k) );
418 warp2(k) = blend2(k) * warpfactor2(k) *
419 ( 1.0 + alpha * alpha * L2(k) * L2(k) );
420 warp3(k) = blend3(k) * warpfactor3(k) *
421 ( 1.0 + alpha * alpha * L3(k) * L3(k) );
422 }
423
424 for (int k=0;k<N;k++) {
425 X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k);
426 Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k);
427 }
428
429 FieldContainer<Scalar> warXY(N,2);
430
431 for (int k=0;k<N;k++) {
432 warXY(k,0) = X(k);
433 warXY(k,1) = Y(k);
434 }
435
436
437 // finally, convert the warp-blend points to the correct triangle
438 FieldContainer<Scalar> warburtonVerts(1,3,2);
439 warburtonVerts(0,0,0) = -1.0;
440 warburtonVerts(0,0,1) = -1.0/sqrt(3.0);
441 warburtonVerts(0,1,0) = 1.0;
442 warburtonVerts(0,1,1) = -1.0/sqrt(3.0);
443 warburtonVerts(0,2,0) = 0.0;
444 warburtonVerts(0,2,1) = 2.0/sqrt(3.0);
445
446 FieldContainer<Scalar> refPts(N,2);
447
449 warXY ,
450 warburtonVerts ,
451 shards::getCellTopologyData< shards::Triangle<3> >(),
452 0 );
453
454 // now write from refPts into points, taking care of offset
455 int noffcur = 0; // index into refPts
456 int offcur = 0; // index int points
457 for (int i=0;i<=order;i++) {
458 for (int j=0;j<=order-i;j++) {
459 if ( (i >= offset) && (i <= order-offset) &&
460 (j >= offset) && (j <= order-i-offset) ) {
461 points(offcur,0) = refPts(noffcur,0);
462 points(offcur,1) = refPts(noffcur,1);
463 offcur++;
464 }
465 noffcur++;
466 }
467 }
468
469 return;
470 }
471
472
473 template<class Scalar, class ArrayType>
474 void PointTools::warpShiftFace3D( const int order ,
475 const Scalar pval ,
476 const ArrayType &L1,
477 const ArrayType &L2,
478 const ArrayType &L3,
479 const ArrayType &L4,
480 ArrayType &dxy)
481 {
482 evalshift<Scalar,ArrayType>(order,pval,L2,L3,L4,dxy);
483 return;
484 }
485
486 template<class Scalar, class ArrayType>
487 void PointTools::evalshift( const int order ,
488 const Scalar pval ,
489 const ArrayType &L1 ,
490 const ArrayType &L2 ,
491 const ArrayType &L3 ,
492 ArrayType &dxy )
493 {
494 // get Gauss-Lobatto-nodes
495 FieldContainer<Scalar> gaussX(order+1,1);
496 PointTools::getWarpBlendLatticeLine<Scalar,FieldContainer<Scalar> >( gaussX , order , 0 );
497 gaussX.resize(order+1);
498 const int N = L1.dimension(0);
499
500 // Warburton code reverses them
501 for (int k=0;k<=order;k++) {
502 gaussX(k) *= -1.0;
503 }
504
505 // blending function at each node for each edge
506 FieldContainer<Scalar> blend1(N);
507 FieldContainer<Scalar> blend2(N);
508 FieldContainer<Scalar> blend3(N);
509
510 for (int i=0;i<N;i++) {
511 blend1(i) = L2(i) * L3(i);
512 blend2(i) = L1(i) * L3(i);
513 blend3(i) = L1(i) * L2(i);
514 }
515
516 // amount of warp for each node for each edge
517 FieldContainer<Scalar> warpfactor1(N);
518 FieldContainer<Scalar> warpfactor2(N);
519 FieldContainer<Scalar> warpfactor3(N);
520
521 // differences of barycentric coordinates
522 FieldContainer<Scalar> L3mL2(N);
523 FieldContainer<Scalar> L1mL3(N);
524 FieldContainer<Scalar> L2mL1(N);
525
526 for (int k=0;k<N;k++) {
527 L3mL2(k) = L3(k)-L2(k);
528 L1mL3(k) = L1(k)-L3(k);
529 L2mL1(k) = L2(k)-L1(k);
530 }
531
532 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor1 , order , gaussX , L3mL2 );
533 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor2 , order , gaussX , L1mL3 );
534 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor3 , order , gaussX , L2mL1 );
535
536 for (int k=0;k<N;k++) {
537 warpfactor1(k) *= 4.0;
538 warpfactor2(k) *= 4.0;
539 warpfactor3(k) *= 4.0;
540 }
541
542 FieldContainer<Scalar> warp1(N);
543 FieldContainer<Scalar> warp2(N);
544 FieldContainer<Scalar> warp3(N);
545
546 for (int k=0;k<N;k++) {
547 warp1(k) = blend1(k) * warpfactor1(k) *
548 ( 1.0 + pval * pval * L1(k) * L1(k) );
549 warp2(k) = blend2(k) * warpfactor2(k) *
550 ( 1.0 + pval * pval * L2(k) * L2(k) );
551 warp3(k) = blend3(k) * warpfactor3(k) *
552 ( 1.0 + pval * pval * L3(k) * L3(k) );
553 }
554
555 for (int k=0;k<N;k++) {
556 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);
557 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);
558 }
559
560 return;
561
562 }
563
564 /* one-d edge warping function */
565 template<class Scalar, class ArrayType>
566 void PointTools::evalwarp( ArrayType &warp ,
567 const int order ,
568 const ArrayType &xnodes ,
569 const ArrayType &xout )
570 {
571 FieldContainer<Scalar> xeq(order+1);
572 FieldContainer<Scalar> d(xout.dimension(0));
573
574 d.initialize();
575
576 for (int i=0;i<=order;i++) {
577 xeq(i) = -1.0 + 2.0 * ( order - i ) / order;
578 }
579
580
581
582 for (int i=0;i<=order;i++) {
583 d.initialize( xnodes(i) - xeq(i) );
584 for (int j=1;j<order;j++) {
585 if (i!=j) {
586 for (int k=0;k<d.dimension(0);k++) {
587 d(k) = d(k) * (xout(k)-xeq(j))/(xeq(i)-xeq(j));
588 }
589 }
590 }
591 if (i!=0) {
592 for (int k=0;k<d.dimension(0);k++) {
593 d(k) = -d(k)/(xeq(i)-xeq(0));
594 }
595 }
596 if (i!=order) {
597 for (int k=0;k<d.dimension(0);k++) {
598 d(k) = d(k)/(xeq(i)-xeq(order));
599 }
600 }
601
602 for (int k=0;k<d.dimension(0);k++) {
603 warp(k) += d(k);
604 }
605 }
606
607 return;
608 }
609
610
611 template<class Scalar, class ArrayType>
613 const int order ,
614 const int offset )
615 {
616 Scalar alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603,
617 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655};
618 Scalar alpha;
619
620 if (order <= 15) {
621 alpha = alphastore[order-1];
622 }
623 else {
624 alpha = 1.0;
625 }
626
627 const int N = (order+1)*(order+2)*(order+3)/6;
628 Scalar tol = 1.e-10;
629
630 FieldContainer<Scalar> shift(N,3);
631 shift.initialize();
632
633 /* create 3d equidistributed nodes on Warburton tet */
634 FieldContainer<Scalar> equipoints(N,3);
635 int sk = 0;
636 for (int n=0;n<=order;n++) {
637 for (int m=0;m<=order-n;m++) {
638 for (int q=0;q<=order-n-m;q++) {
639 equipoints(sk,0) = -1.0 + (q * 2.0 ) / order;
640 equipoints(sk,1) = -1.0 + (m * 2.0 ) / order;
641 equipoints(sk,2) = -1.0 + (n * 2.0 ) / order;
642 sk++;
643 }
644 }
645 }
646
647
648 /* construct barycentric coordinates for equispaced lattice */
649 FieldContainer<Scalar> L1(N);
650 FieldContainer<Scalar> L2(N);
651 FieldContainer<Scalar> L3(N);
652 FieldContainer<Scalar> L4(N);
653 for (int i=0;i<N;i++) {
654 L1(i) = (1.0 + equipoints(i,2)) / 2.0;
655 L2(i) = (1.0 + equipoints(i,1)) / 2.0;
656 L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0;
657 L4(i) = (1.0 + equipoints(i,0)) / 2.0;
658 }
659
660
661 /* vertices of Warburton tet */
662 FieldContainer<Scalar> warVerts(4,3);
663 warVerts(0,0) = -1.0;
664 warVerts(0,1) = -1.0/sqrt(3.0);
665 warVerts(0,2) = -1.0/sqrt(6.0);
666 warVerts(1,0) = 1.0;
667 warVerts(1,1) = -1.0/sqrt(3.0);
668 warVerts(1,2) = -1.0/sqrt(6.0);
669 warVerts(2,0) = 0.0;
670 warVerts(2,1) = 2.0 / sqrt(3.0);
671 warVerts(2,2) = -1.0/sqrt(6.0);
672 warVerts(3,0) = 0.0;
673 warVerts(3,1) = 0.0;
674 warVerts(3,2) = 3.0 / sqrt(6.0);
675
676
677 /* tangents to faces */
678 FieldContainer<Scalar> t1(4,3);
679 FieldContainer<Scalar> t2(4,3);
680 for (int i=0;i<3;i++) {
681 t1(0,i) = warVerts(1,i) - warVerts(0,i);
682 t1(1,i) = warVerts(1,i) - warVerts(0,i);
683 t1(2,i) = warVerts(2,i) - warVerts(1,i);
684 t1(3,i) = warVerts(2,i) - warVerts(0,i);
685 t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
686 t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
687 t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) );
688 t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) );
689 }
690
691 /* normalize tangents */
692 for (int n=0;n<4;n++) {
693 /* Compute norm of t1(n) and t2(n) */
694 Scalar normt1n = 0.0;
695 Scalar normt2n = 0.0;
696 for (int i=0;i<3;i++) {
697 normt1n += (t1(n,i) * t1(n,i));
698 normt2n += (t2(n,i) * t2(n,i));
699 }
700 normt1n = sqrt(normt1n);
701 normt2n = sqrt(normt2n);
702 /* normalize each tangent now */
703 for (int i=0;i<3;i++) {
704 t1(n,i) /= normt1n;
705 t2(n,i) /= normt2n;
706 }
707 }
708
709 /* undeformed coordinates */
710 FieldContainer<Scalar> XYZ(N,3);
711 for (int i=0;i<N;i++) {
712 for (int j=0;j<3;j++) {
713 XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j);
714 }
715 }
716
717 for (int face=1;face<=4;face++) {
718 FieldContainer<Scalar> La, Lb, Lc, Ld;
719 FieldContainer<Scalar> warp(N,2);
720 FieldContainer<Scalar> blend(N);
721 FieldContainer<Scalar> denom(N);
722 switch (face) {
723 case 1:
724 La = L1; Lb = L2; Lc = L3; Ld = L4; break;
725 case 2:
726 La = L2; Lb = L1; Lc = L3; Ld = L4; break;
727 case 3:
728 La = L3; Lb = L1; Lc = L4; Ld = L2; break;
729 case 4:
730 La = L4; Lb = L1; Lc = L3; Ld = L2; break;
731 }
732
733 /* get warp tangential to face */
734 warpShiftFace3D<Scalar,ArrayType>(order,alpha,La,Lb,Lc,Ld,warp);
735
736 for (int k=0;k<N;k++) {
737 blend(k) = Lb(k) * Lc(k) * Ld(k);
738 }
739
740 for (int k=0;k<N;k++) {
741 denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k));
742 }
743
744 for (int k=0;k<N;k++) {
745 if (denom(k) > tol) {
746 blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k);
747 }
748 }
749
750
751 // compute warp and blend
752 for (int k=0;k<N;k++) {
753 for (int j=0;j<3;j++) {
754 shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j)
755 + blend(k) * warp(k,1) * t2(face-1,j);
756 }
757 }
758
759 for (int k=0;k<N;k++) {
760 if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) {
761 for (int j=0;j<3;j++) {
762 shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j);
763 }
764 }
765 }
766
767 }
768
769 FieldContainer<Scalar> updatedPoints(N,3);
770 for (int k=0;k<N;k++) {
771 for (int j=0;j<3;j++) {
772 updatedPoints(k,j) = XYZ(k,j) + shift(k,j);
773 }
774 }
775
776 warVerts.resize( 1 , 4 , 3 );
777
778 // now we convert to Pavel's reference triangle!
779 FieldContainer<Scalar> refPts(N,3);
780 CellTools<Scalar>::mapToReferenceFrame( refPts ,updatedPoints ,
781 warVerts ,
782 shards::getCellTopologyData<shards::Tetrahedron<4> >() ,
783 0 );
784
785 // now write from refPts into points, taking offset into account
786 int noffcur = 0;
787 int offcur = 0;
788 for (int i=0;i<=order;i++) {
789 for (int j=0;j<=order-i;j++) {
790 for (int k=0;k<=order-i-j;k++) {
791 if ( (i >= offset) && (i <= order-offset) &&
792 (j >= offset) && (j <= order-i-offset) &&
793 (k >= offset) && (k <= order-i-j-offset) ) {
794 points(offcur,0) = refPts(noffcur,0);
795 points(offcur,1) = refPts(noffcur,1);
796 points(offcur,2) = refPts(noffcur,2);
797 offcur++;
798 }
799 noffcur++;
800 }
801 }
802 }
803
804
805
806 }
807
808
809} // face Intrepid
static void mapToReferenceFrame(ArrayRefPoint &refPoints, const ArrayPhysPoint &physPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes , the inverse of the reference-to-physical frame map using a default initial guess.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
static void zwglj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
static void zwgj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Jacobi zeros and weights.
static void getWarpBlendLatticeTetrahedron(ArrayType &points, const int order, const int offset=0)
Returns Warburton's warp-blend points of a given order on the reference tetrahedron....
static void warpFactor(const int order, const ArrayType &xnodes, const ArrayType &xout, ArrayType &warp)
interpolates Warburton's warp function on the line
static void getLattice(ArrayType &pts, const shards::CellTopology &cellType, const int order, const int offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference simplex (currently disabled for other ce...
static void getWarpBlendLatticeTriangle(ArrayType &points, const int order, const int offset=0)
Returns Warburton's warp-blend points of a given order on the reference triangle. The output array is...
static void getGaussPoints(ArrayType &pts, const int order)
static void getEquispacedLatticeTriangle(ArrayType &points, const int order, const int offset=0)
Computes an equispaced lattice of a given order on the reference triangle. The output array is (P,...
static void getEquispacedLattice(const shards::CellTopology &cellType, ArrayType &points, const int order, const int offset=0)
Computes an equispaced lattice of a given order on a reference simplex (currently disabled for other ...
static void evalshift(const int order, const Scalar pval, const ArrayType &L1, const ArrayType &L2, const ArrayType &L3, ArrayType &dxy)
Used internally to evaluate the point shift for warp-blend points on faces of tets.
static void evalwarp(ArrayType &warp, const int order, const ArrayType &xnodes, const ArrayType &xout)
Used internally to compute the warp on edges of a triangle in warp-blend points.
static void getWarpBlendLatticeLine(ArrayType &points, const int order, const int offset=0)
Returns the Gauss-Lobatto points of a given order on the reference line [-1,1]. The output array is (...
static void getWarpBlendLattice(const shards::CellTopology &cellType, ArrayType &points, const int order, const int offset=0)
Computes a warped lattice (ala Warburton's warp-blend points of a given order on a reference simplex ...
static void getEquispacedLatticeLine(ArrayType &points, const int order, const int offset=0)
Computes an equispaced lattice of a given order on the reference line [-1,1]. The output array is (P,...
static void getEquispacedLatticeTetrahedron(ArrayType &points, const int order, const int offset=0)
Computes an equispaced lattice of a given order on the reference tetrahedron. The output array is (P,...
static int getLatticeSize(const shards::CellTopology &cellType, const int order, const int offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...
static void warpShiftFace3D(const int order, const Scalar pval, const ArrayType &L1, const ArrayType &L2, const ArrayType &L3, const ArrayType &L4, ArrayType &dxy)
This is used internally to compute the tetrahedral warp-blend points one each face.