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 
49 namespace 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>
332  void PointTools::getWarpBlendLatticeTriangle( ArrayType &points ,
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 */
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
810 
811 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
812 #ifdef __GNUC__
813 #warning "The Intrepid package is deprecated"
814 #endif
815 #endif
816 
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...