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