Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ElemUtil.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <limits>
3 #include <cassert>
4 
5 #include "ElemUtil.hpp"
6 #include "moab/BoundBox.hpp"
7 
8 namespace moab
9 {
10 namespace ElemUtil
11 {
12 
13  bool debug = false;
14 
15  /**\brief Class representing a 3-D mapping function (e.g. shape function for volume element) */
16  class VolMap
17  {
18  public:
19  /**\brief Return \f$\vec \xi\f$ corresponding to logical center of element */
20  virtual CartVect center_xi() const = 0;
21  /**\brief Evaluate mapping function (calculate \f$\vec x = F($\vec \xi)\f$ )*/
22  virtual CartVect evaluate( const CartVect& xi ) const = 0;
23  /**\brief Evaluate Jacobian of mapping function */
24  virtual Matrix3 jacobian( const CartVect& xi ) const = 0;
25  /**\brief Evaluate inverse of mapping function (calculate \f$\vec \xi = F^-1($\vec x)\f$ )*/
26  bool solve_inverse( const CartVect& x, CartVect& xi, double tol ) const;
27  };
28 
29  bool VolMap::solve_inverse( const CartVect& x, CartVect& xi, double tol ) const
30  {
31  const double error_tol_sqr = tol * tol;
32  double det;
33  xi = center_xi();
34  CartVect delta = evaluate( xi ) - x;
35  Matrix3 J;
36  while( delta % delta > error_tol_sqr )
37  {
38  J = jacobian( xi );
39  det = J.determinant();
40  if( det < std::numeric_limits< double >::epsilon() ) return false;
41  xi -= J.inverse() * delta;
42  delta = evaluate( xi ) - x;
43  }
44  return true;
45  }
46 
47  /**\brief Shape function for trilinear hexahedron */
48  class LinearHexMap : public VolMap
49  {
50  public:
51  LinearHexMap( const CartVect* corner_coords ) : corners( corner_coords ) {}
52  virtual CartVect center_xi() const;
53  virtual CartVect evaluate( const CartVect& xi ) const;
54  virtual double evaluate_scalar_field( const CartVect& xi, const double* f_vals ) const;
55  virtual Matrix3 jacobian( const CartVect& xi ) const;
56 
57  private:
58  const CartVect* corners;
59  static const double corner_xi[8][3];
60  };
61 
62  const double LinearHexMap::corner_xi[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
63  { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } };
65  {
66  return CartVect( 0.0 );
67  }
68 
70  {
71  CartVect x( 0.0 );
72  for( unsigned i = 0; i < 8; ++i )
73  {
74  const double N_i =
75  ( 1 + xi[0] * corner_xi[i][0] ) * ( 1 + xi[1] * corner_xi[i][1] ) * ( 1 + xi[2] * corner_xi[i][2] );
76  x += N_i * corners[i];
77  }
78  x *= 0.125;
79  return x;
80  }
81 
82  double LinearHexMap::evaluate_scalar_field( const CartVect& xi, const double* f_vals ) const
83  {
84  double f( 0.0 );
85  for( unsigned i = 0; i < 8; ++i )
86  {
87  const double N_i =
88  ( 1 + xi[0] * corner_xi[i][0] ) * ( 1 + xi[1] * corner_xi[i][1] ) * ( 1 + xi[2] * corner_xi[i][2] );
89  f += N_i * f_vals[i];
90  }
91  f *= 0.125;
92  return f;
93  }
94 
96  {
97  Matrix3 J( 0.0 );
98  for( unsigned i = 0; i < 8; ++i )
99  {
100  const double xi_p = 1 + xi[0] * corner_xi[i][0];
101  const double eta_p = 1 + xi[1] * corner_xi[i][1];
102  const double zeta_p = 1 + xi[2] * corner_xi[i][2];
103  const double dNi_dxi = corner_xi[i][0] * eta_p * zeta_p;
104  const double dNi_deta = corner_xi[i][1] * xi_p * zeta_p;
105  const double dNi_dzeta = corner_xi[i][2] * xi_p * eta_p;
106  J( 0, 0 ) += dNi_dxi * corners[i][0];
107  J( 1, 0 ) += dNi_dxi * corners[i][1];
108  J( 2, 0 ) += dNi_dxi * corners[i][2];
109  J( 0, 1 ) += dNi_deta * corners[i][0];
110  J( 1, 1 ) += dNi_deta * corners[i][1];
111  J( 2, 1 ) += dNi_deta * corners[i][2];
112  J( 0, 2 ) += dNi_dzeta * corners[i][0];
113  J( 1, 2 ) += dNi_dzeta * corners[i][1];
114  J( 2, 2 ) += dNi_dzeta * corners[i][2];
115  }
116  return J *= 0.125;
117  }
118 
119  bool nat_coords_trilinear_hex( const CartVect* corner_coords, const CartVect& x, CartVect& xi, double tol )
120  {
121  return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
122  }
123  //
124  // nat_coords_trilinear_hex2
125  // Duplicate functionality of nat_coords_trilinear_hex using hex_findpt
126  //
127  void nat_coords_trilinear_hex2( const CartVect hex[8], const CartVect& xyz, CartVect& ncoords, double etol )
128 
129  {
130  const int ndim = 3;
131  const int nverts = 8;
132  const int vertMap[nverts] = { 0, 1, 3, 2, 4, 5, 7, 6 }; // Map from nat to lex ordering
133 
134  const int n = 2; // linear
135  realType coords[ndim * nverts]; // buffer
136 
137  realType* xm[ndim];
138  for( int i = 0; i < ndim; i++ )
139  xm[i] = coords + i * nverts;
140 
141  // stuff hex into coords
142  for( int i = 0; i < nverts; i++ )
143  {
144  realType vcoord[ndim];
145  hex[i].get( vcoord );
146 
147  for( int d = 0; d < ndim; d++ )
148  coords[d * nverts + vertMap[i]] = vcoord[d];
149  }
150 
151  double dist = 0.0;
152  ElemUtil::hex_findpt( xm, n, xyz, ncoords, dist );
153  if( 3 * MOAB_POLY_EPS < dist )
154  {
155  // outside element, set extremal values to something outside range
156  for( int j = 0; j < 3; j++ )
157  {
158  if( ncoords[j] < ( -1.0 - etol ) || ncoords[j] > ( 1.0 + etol ) ) ncoords[j] *= 10;
159  }
160  }
161  }
162  bool point_in_trilinear_hex( const CartVect* hex, const CartVect& xyz, double etol )
163  {
164  CartVect xi;
165  return nat_coords_trilinear_hex( hex, xyz, xi, etol ) && ( fabs( xi[0] - 1.0 ) < etol ) &&
166  ( fabs( xi[1] - 1.0 ) < etol ) && ( fabs( xi[2] - 1.0 ) < etol );
167  }
168 
170  const CartVect& xyz,
171  const CartVect& box_min,
172  const CartVect& box_max,
173  double etol )
174  {
175  // all values scaled by 2 (eliminates 3 flops)
176  const CartVect mid = box_max + box_min;
177  const CartVect dim = box_max - box_min;
178  const CartVect pt = 2 * xyz - mid;
179  return ( fabs( pt[0] - dim[0] ) < etol ) && ( fabs( pt[1] - dim[1] ) < etol ) &&
180  ( fabs( pt[2] - dim[2] ) < etol ) && point_in_trilinear_hex( hex, xyz, etol );
181  }
182 
183  // Wrapper to James Lottes' findpt routines
184  // hex_findpt
185  // Find the parametric coordinates of an xyz point inside
186  // a 3d hex spectral element with n nodes per dimension
187  // xm: coordinates fields, value of x,y,z for each of then n*n*n gauss-lobatto nodes. Nodes are
188  // in lexicographical order (x is fastest-changing) n: number of nodes per dimension -- n=2 for
189  // a linear element xyz: input, point to find rst: output: parametric coords of xyz inside the
190  // element. If xyz is outside the element, rst will be the coords of the closest point dist:
191  // output: distance between xyz and the point with parametric coords rst
192  /*extern "C"{
193  #include "types.h"
194  #include "poly.h"
195  #include "tensor.h"
196  #include "findpt.h"
197  #include "extrafindpt.h"
198  #include "errmem.h"
199  }*/
200 
201  void hex_findpt( realType* xm[3], int n, CartVect xyz, CartVect& rst, double& dist )
202  {
203 
204  // compute stuff that only depends on the order -- could be cached
205  realType* z[3];
206  lagrange_data ld[3];
207  opt_data_3 data;
208 
209  // triplicates
210  for( int d = 0; d < 3; d++ )
211  {
212  z[d] = tmalloc( realType, n );
213  lobatto_nodes( z[d], n );
214  lagrange_setup( &ld[d], z[d], n );
215  }
216 
217  opt_alloc_3( &data, ld );
218 
219  // find nearest point
220  realType x_star[3];
221  xyz.get( x_star );
222 
223  realType r[3] = { 0, 0, 0 }; // initial guess for parametric coords
224  unsigned c = opt_no_constraints_3;
225  dist = opt_findpt_3( &data, (const realType**)xm, x_star, r, &c );
226  // c tells us if we landed inside the element or exactly on a face, edge, or node
227 
228  // copy parametric coords back
229  rst = r;
230 
231  // Clean-up (move to destructor if we decide to cache)
232  opt_free_3( &data );
233  for( int d = 0; d < 3; ++d )
234  lagrange_free( &ld[d] );
235  for( int d = 0; d < 3; ++d )
236  free( z[d] );
237  }
238 
239  // hex_eval
240  // Evaluate a field in a 3d hex spectral element with n nodes per dimension, at some given
241  // parametric coordinates field: field values for each of then n*n*n gauss-lobatto nodes. Nodes
242  // are in lexicographical order (x is fastest-changing) n: number of nodes per dimension -- n=2
243  // for a linear element rst: input: parametric coords of the point where we want to evaluate the
244  // field value: output: value of field at rst
245 
246  void hex_eval( realType* field, int n, CartVect rstCartVec, double& value )
247  {
248  int d;
249  realType rst[3];
250  rstCartVec.get( rst );
251 
252  // can cache stuff below
253  lagrange_data ld[3];
254  realType* z[3];
255  for( d = 0; d < 3; ++d )
256  {
257  z[d] = tmalloc( realType, n );
258  lobatto_nodes( z[d], n );
259  lagrange_setup( &ld[d], z[d], n );
260  }
261 
262  // cut and paste -- see findpt.c
263  const unsigned nf = n * n, ne = n, nw = 2 * n * n + 3 * n;
264  realType* od_work = tmalloc( realType, 6 * nf + 9 * ne + nw );
265 
266  // piece that we shouldn't want to cache
267  for( d = 0; d < 3; d++ )
268  {
269  lagrange_0( &ld[d], rst[d] );
270  }
271 
272  value = tensor_i3( ld[0].J, ld[0].n, ld[1].J, ld[1].n, ld[2].J, ld[2].n, field, od_work );
273 
274  // all this could be cached
275  for( d = 0; d < 3; d++ )
276  {
277  free( z[d] );
278  lagrange_free( &ld[d] );
279  }
280  free( od_work );
281  }
282 
283  // Gaussian quadrature points for a trilinear hex element.
284  // Five 2d arrays are defined.
285  // One for the single gaussian point solution, 2 point solution,
286  // 3 point solution, 4 point solution and 5 point solution.
287  // There are 2 columns, one for Weights and the other for Locations
288  // Weight Location
289 
290  const double gauss_1[1][2] = { { 2.0, 0.0 } };
291 
292  const double gauss_2[2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } };
293 
294  const double gauss_3[3][2] = { { 0.5555555555, -0.7745966692 },
295  { 0.8888888888, 0.0 },
296  { 0.5555555555, 0.7745966692 } };
297 
298  const double gauss_4[4][2] = { { 0.3478548451, -0.8611363116 },
299  { 0.6521451549, -0.3399810436 },
300  { 0.6521451549, 0.3399810436 },
301  { 0.3478548451, 0.8611363116 } };
302 
303  const double gauss_5[5][2] = { { 0.2369268851, -0.9061798459 },
304  { 0.4786286705, -0.5384693101 },
305  { 0.5688888889, 0.0 },
306  { 0.4786286705, 0.5384693101 },
307  { 0.2369268851, 0.9061798459 } };
308 
309  // Function to integrate the field defined by field_fn function
310  // over the volume of the trilinear hex defined by the hex_corners
311 
312  bool integrate_trilinear_hex( const CartVect* hex_corners, double* corner_fields, double& field_val, int num_pts )
313  {
314  // Create the LinearHexMap object using the hex_corners array of CartVects
315  LinearHexMap hex( hex_corners );
316 
317  // Use the correct table of points and locations based on the num_pts parameter
318  const double( *g_pts )[2] = 0;
319  switch( num_pts )
320  {
321  case 1:
322  g_pts = gauss_1;
323  break;
324 
325  case 2:
326  g_pts = gauss_2;
327  break;
328 
329  case 3:
330  g_pts = gauss_3;
331  break;
332 
333  case 4:
334  g_pts = gauss_4;
335  break;
336 
337  case 5:
338  g_pts = gauss_5;
339  break;
340 
341  default: // value out of range
342  return false;
343  }
344 
345  // Test code - print Gaussian Quadrature data
346  if( debug )
347  {
348  for( int r = 0; r < num_pts; r++ )
349  for( int c = 0; c < 2; c++ )
350  std::cout << "g_pts[" << r << "][" << c << "]=" << g_pts[r][c] << std::endl;
351  }
352  // End Test code
353 
354  double soln = 0.0;
355 
356  for( int i = 0; i < num_pts; i++ )
357  { // Loop for xi direction
358  double w_i = g_pts[i][0];
359  double xi_i = g_pts[i][1];
360  for( int j = 0; j < num_pts; j++ )
361  { // Loop for eta direction
362  double w_j = g_pts[j][0];
363  double eta_j = g_pts[j][1];
364  for( int k = 0; k < num_pts; k++ )
365  { // Loop for zeta direction
366  double w_k = g_pts[k][0];
367  double zeta_k = g_pts[k][1];
368 
369  // Calculate the "realType" space point given the "normal" point
370  CartVect normal_pt( xi_i, eta_j, zeta_k );
371 
372  // Calculate the value of F(x(xi,eta,zeta),y(xi,eta,zeta),z(xi,eta,zeta)
373  double field = hex.evaluate_scalar_field( normal_pt, corner_fields );
374 
375  // Calculate the Jacobian for this "normal" point and its determinant
376  Matrix3 J = hex.jacobian( normal_pt );
377  double det = J.determinant();
378 
379  // Calculate integral and update the solution
380  soln = soln + ( w_i * w_j * w_k * field * det );
381  }
382  }
383  }
384 
385  // Set the output parameter
386  field_val = soln;
387 
388  return true;
389  }
390 
391 } // namespace ElemUtil
392 
393 namespace Element
394 {
395 
397 
398  inline const std::vector< CartVect >& Map::get_vertices()
399  {
400  return this->vertex;
401  }
402  //
403  void Map::set_vertices( const std::vector< CartVect >& v )
404  {
405  if( v.size() != this->vertex.size() )
406  {
407  throw ArgError();
408  }
409  this->vertex = v;
410  }
411 
412  bool Map::inside_box( const CartVect& xi, double& tol ) const
413  {
414  // bail out early, before doing an expensive NR iteration
415  // compute box
416  BoundBox box( this->vertex );
417  return box.contains_point( xi.array(), tol );
418  }
419 
420  //
421  CartVect Map::ievaluate( const CartVect& x, double tol, const CartVect& x0 ) const
422  {
423  // TODO: should differentiate between epsilons used for
424  // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
425  // right now, fix the tolerance used for NR
426  tol = 1.0e-10;
427  const double error_tol_sqr = tol * tol;
428  double det;
429  CartVect xi = x0;
430  CartVect delta = evaluate( xi ) - x;
431  Matrix3 J;
432 
433  int iters = 0;
434  while( delta % delta > error_tol_sqr )
435  {
436  if( ++iters > 10 ) throw Map::EvaluationError( x, vertex );
437 
438  J = jacobian( xi );
439  det = J.determinant();
440  if( det < std::numeric_limits< double >::epsilon() ) throw Map::EvaluationError( x, vertex );
441  xi -= J.inverse() * delta;
442  delta = evaluate( xi ) - x;
443  }
444  return xi;
445  } // Map::ievaluate()
446 
447  SphericalQuad::SphericalQuad( const std::vector< CartVect >& vertices ) : LinearQuad( vertices )
448  {
449  // project the vertices to the plane tangent at first vertex
450  v1 = vertex[0]; // member data
451  double v1v1 = v1 % v1; // this is 1, in general, for unit sphere meshes
452  for( int j = 1; j < 4; j++ )
453  {
454  // first, bring all vertices in the gnomonic plane
455  // the new vertex will intersect the plane at vnew
456  // so that (vnew-v1)%v1 is 0 ( vnew is in the tangent plane, i.e. normal to v1 )
457  // pos is the old position of the vertex, and it is in general on the sphere
458  // vnew = alfa*pos; (alfa*pos-v1)%v1 = 0 <=> alfa*(pos%v1)=v1%v1 <=> alfa =
459  // v1v1/(pos%v1)
460  // <=> vnew = ( v1v1/(pos%v1) )*pos
461  CartVect vnew = v1v1 / ( vertex[j] % v1 ) * vertex[j];
462  vertex[j] = vnew;
463  }
464  // will compute a transf matrix, such that a new point will be transformed with
465  // newpos = transf * (vnew-v1), and it will be a point in the 2d plane
466  // the transformation matrix will be oriented in such a way that orientation will be
467  // positive
468  CartVect vx = vertex[1] - v1; // this will become Ox axis
469  // z axis will be along v1, in such a way that orientation of the quad is positive
470  // look at the first 2 edges
471  CartVect vz = vx * ( vertex[2] - vertex[1] );
472  vz = vz / vz.length();
473 
474  vx = vx / vx.length();
475 
476  CartVect vy = vz * vx;
477  transf = Matrix3( vx[0], vx[1], vx[2], vy[0], vy[1], vy[2], vz[0], vz[1], vz[2] );
478  vertex[0] = CartVect( 0. );
479  for( int j = 1; j < 4; j++ )
480  vertex[j] = transf * ( vertex[j] - v1 );
481  }
482 
483  CartVect SphericalQuad::ievaluate( const CartVect& x, double tol, const CartVect& x0 ) const
484  {
485  // project to the plane tangent at first vertex (gnomonic projection)
486  double v1v1 = v1 % v1;
487  CartVect vnew = v1v1 / ( x % v1 ) * x; // so that (vnew-v1)%v1 is 0
488  vnew = transf * ( vnew - v1 );
489  // det will be positive now
490  return Map::ievaluate( vnew, tol, x0 );
491  }
492 
493  bool SphericalQuad::inside_box( const CartVect& pos, double& tol ) const
494  {
495  // project to the plane tangent at first vertex
496  // CartVect v1=vertex[0];
497  double v1v1 = v1 % v1;
498  CartVect vnew = v1v1 / ( pos % v1 ) * pos; // so that (x-v1)%v1 is 0
499  vnew = transf * ( vnew - v1 );
500  return Map::inside_box( vnew, tol );
501  }
502 
503  const double LinearTri::corner[3][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 } };
504 
505  LinearTri::LinearTri() : Map( 0 ), det_T( 0.0 ), det_T_inverse( 0.0 ) {} // LinearTri::LinearTri()
506 
508 
509  void LinearTri::set_vertices( const std::vector< CartVect >& v )
510  {
511  this->Map::set_vertices( v );
512  this->T = Matrix3( v[1][0] - v[0][0], v[2][0] - v[0][0], 0, v[1][1] - v[0][1], v[2][1] - v[0][1], 0,
513  v[1][2] - v[0][2], v[2][2] - v[0][2], 1 );
514  this->T_inverse = this->T.inverse();
515  this->det_T = this->T.determinant();
516  this->det_T_inverse = ( this->det_T < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / this->det_T );
517  } // LinearTri::set_vertices()
518 
519  bool LinearTri::inside_nat_space( const CartVect& xi, double& tol ) const
520  {
521  // linear tri space is a triangle with vertices (0,0,0), (1,0,0), (0,1,0)
522  // first check if outside bigger box, then below the line x+y=1
523  return ( xi[0] >= -tol ) && ( xi[1] >= -tol ) && ( xi[2] >= -tol ) && ( xi[2] <= tol ) &&
524  ( xi[0] + xi[1] < 1.0 + tol );
525  }
526  CartVect LinearTri::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const
527  {
528  return this->T_inverse * ( x - this->vertex[0] );
529  } // LinearTri::ievaluate
530 
531  double LinearTri::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
532  {
533  double f0 = field_vertex_value[0];
534  double f = f0;
535  for( unsigned i = 1; i < 3; ++i )
536  {
537  f += ( field_vertex_value[i] - f0 ) * xi[i - 1];
538  }
539  return f;
540  } // LinearTri::evaluate_scalar_field()
541 
542  double LinearTri::integrate_scalar_field( const double* field_vertex_values ) const
543  {
544  double I( 0.0 );
545  for( unsigned int i = 0; i < 3; ++i )
546  {
547  I += field_vertex_values[i];
548  }
549  I *= this->det_T / 6.0; // TODO
550  return I;
551  } // LinearTri::integrate_scalar_field()
552 
553  SphericalTri::SphericalTri( const std::vector< CartVect >& vertices )
554  {
555  vertex.resize( vertices.size() );
556  vertex = vertices;
557  // project the vertices to the plane tangent at first vertex
558  v1 = vertex[0]; // member data
559  double v1v1 = v1 % v1; // this is 1, in general, for unit sphere meshes
560  for( int j = 1; j < 3; j++ )
561  {
562  // first, bring all vertices in the gnomonic plane
563  // the new vertex will intersect the plane at vnew
564  // so that (vnew-v1)%v1 is 0 ( vnew is in the tangent plane, i.e. normal to v1 )
565  // pos is the old position of the vertex, and it is in general on the sphere
566  // vnew = alfa*pos; (alfa*pos-v1)%v1 = 0 <=> alfa*(pos%v1)=v1%v1 <=> alfa =
567  // v1v1/(pos%v1)
568  // <=> vnew = ( v1v1/(pos%v1) )*pos
569  CartVect vnew = v1v1 / ( vertex[j] % v1 ) * vertex[j];
570  vertex[j] = vnew;
571  }
572  // will compute a transf matrix, such that a new point will be transformed with
573  // newpos = transf * (vnew-v1), and it will be a point in the 2d plane
574  // the transformation matrix will be oriented in such a way that orientation will be
575  // positive
576  CartVect vx = vertex[1] - v1; // this will become Ox axis
577  // z axis will be along v1, in such a way that orientation of the quad is positive
578  // look at the first 2 edges
579  CartVect vz = vx * ( vertex[2] - vertex[1] );
580  vz = vz / vz.length();
581 
582  vx = vx / vx.length();
583 
584  CartVect vy = vz * vx;
585  transf = Matrix3( vx[0], vx[1], vx[2], vy[0], vy[1], vy[2], vz[0], vz[1], vz[2] );
586  vertex[0] = CartVect( 0. );
587  for( int j = 1; j < 3; j++ )
588  vertex[j] = transf * ( vertex[j] - v1 );
589 
591  }
592 
593  CartVect SphericalTri::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const
594  {
595  // project to the plane tangent at first vertex (gnomonic projection)
596  double v1v1 = v1 % v1;
597  CartVect vnew = v1v1 / ( x % v1 ) * x; // so that (vnew-v1)%v1 is 0
598  vnew = transf * ( vnew - v1 );
599  // det will be positive now
600  return LinearTri::ievaluate( vnew );
601  }
602 
603  bool SphericalTri::inside_box( const CartVect& pos, double& tol ) const
604  {
605  // project to the plane tangent at first vertex
606  // CartVect v1=vertex[0];
607  double v1v1 = v1 % v1;
608  CartVect vnew = v1v1 / ( pos % v1 ) * pos; // so that (x-v1)%v1 is 0
609  vnew = transf * ( vnew - v1 );
610  return Map::inside_box( vnew, tol );
611  }
612 
613  // filescope for static member data that is cached
614  const double LinearEdge::corner[2][3] = { { -1, 0, 0 }, { 1, 0, 0 } };
615 
616  LinearEdge::LinearEdge() : Map( 0 ) {} // LinearEdge::LinearEdge()
617 
618  /* For each point, its weight and location are stored as an array.
619  Hence, the inner dimension is 2, the outer dimension is gauss_count.
620  We use a one-point Gaussian quadrature, since it integrates linear functions exactly.
621  */
622  const double LinearEdge::gauss[1][2] = { { 2.0, 0.0 } };
623 
625  {
626  CartVect x( 0.0 );
627  for( unsigned i = 0; i < LinearEdge::corner_count; ++i )
628  {
629  const double N_i = ( 1.0 + xi[0] * corner[i][0] );
630  x += N_i * this->vertex[i];
631  }
633  return x;
634  } // LinearEdge::evaluate
635 
637  {
638  Matrix3 J( 0.0 );
639  for( unsigned i = 0; i < LinearEdge::corner_count; ++i )
640  {
641  const double xi_p = 1.0 + xi[0] * corner[i][0];
642  const double dNi_dxi = corner[i][0] * xi_p;
643  J( 0, 0 ) += dNi_dxi * vertex[i][0];
644  }
645  J( 1, 1 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */
646  J( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */
648  return J;
649  } // LinearEdge::jacobian()
650 
651  double LinearEdge::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
652  {
653  double f( 0.0 );
654  for( unsigned i = 0; i < LinearEdge::corner_count; ++i )
655  {
656  const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1.0 + xi[1] * corner[i][1] );
657  f += N_i * field_vertex_value[i];
658  }
660  return f;
661  } // LinearEdge::evaluate_scalar_field()
662 
663  double LinearEdge::integrate_scalar_field( const double* field_vertex_values ) const
664  {
665  double I( 0.0 );
666  for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 )
667  {
668  double x1 = this->gauss[j1][1];
669  double w1 = this->gauss[j1][0];
670  CartVect x( x1, 0.0, 0.0 );
671  I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * this->det_jacobian( x );
672  }
673  return I;
674  } // LinearEdge::integrate_scalar_field()
675 
676  bool LinearEdge::inside_nat_space( const CartVect& xi, double& tol ) const
677  {
678  // just look at the box+tol here
679  return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol );
680  }
681 
682  const double LinearHex::corner[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
683  { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } };
684 
685  LinearHex::LinearHex() : Map( 0 ) {} // LinearHex::LinearHex()
686 
688  /* For each point, its weight and location are stored as an array.
689  Hence, the inner dimension is 2, the outer dimension is gauss_count.
690  We use a one-point Gaussian quadrature, since it integrates linear functions exactly.
691  */
692  // const double LinearHex::gauss[1][2] = { { 2.0, 0.0 } };
693  const double LinearHex::gauss[2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } };
694  // const double LinearHex::gauss[4][2] = { { 0.3478548451, -0.8611363116 },
695  // { 0.6521451549, -0.3399810436 },
696  // { 0.6521451549, 0.3399810436 },
697  // { 0.3478548451, 0.8611363116 } };
698 
700  {
701  CartVect x( 0.0 );
702  for( unsigned i = 0; i < 8; ++i )
703  {
704  const double N_i =
705  ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] ) * ( 1 + xi[2] * corner[i][2] );
706  x += N_i * this->vertex[i];
707  }
708  x *= 0.125;
709  return x;
710  } // LinearHex::evaluate
711 
713  {
714  Matrix3 J( 0.0 );
715  for( unsigned i = 0; i < 8; ++i )
716  {
717  const double xi_p = 1 + xi[0] * corner[i][0];
718  const double eta_p = 1 + xi[1] * corner[i][1];
719  const double zeta_p = 1 + xi[2] * corner[i][2];
720  const double dNi_dxi = corner[i][0] * eta_p * zeta_p;
721  const double dNi_deta = corner[i][1] * xi_p * zeta_p;
722  const double dNi_dzeta = corner[i][2] * xi_p * eta_p;
723  J( 0, 0 ) += dNi_dxi * vertex[i][0];
724  J( 1, 0 ) += dNi_dxi * vertex[i][1];
725  J( 2, 0 ) += dNi_dxi * vertex[i][2];
726  J( 0, 1 ) += dNi_deta * vertex[i][0];
727  J( 1, 1 ) += dNi_deta * vertex[i][1];
728  J( 2, 1 ) += dNi_deta * vertex[i][2];
729  J( 0, 2 ) += dNi_dzeta * vertex[i][0];
730  J( 1, 2 ) += dNi_dzeta * vertex[i][1];
731  J( 2, 2 ) += dNi_dzeta * vertex[i][2];
732  }
733  return J *= 0.125;
734  } // LinearHex::jacobian()
735 
736  double LinearHex::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
737  {
738  double f( 0.0 );
739  for( unsigned i = 0; i < 8; ++i )
740  {
741  const double N_i =
742  ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] ) * ( 1 + xi[2] * corner[i][2] );
743  f += N_i * field_vertex_value[i];
744  }
745  f *= 0.125;
746  return f;
747  } // LinearHex::evaluate_scalar_field()
748 
749  double LinearHex::integrate_scalar_field( const double* field_vertex_values ) const
750  {
751  double I( 0.0 );
752  for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 )
753  {
754  double x1 = this->gauss[j1][1];
755  double w1 = this->gauss[j1][0];
756  for( unsigned int j2 = 0; j2 < this->gauss_count; ++j2 )
757  {
758  double x2 = this->gauss[j2][1];
759  double w2 = this->gauss[j2][0];
760  for( unsigned int j3 = 0; j3 < this->gauss_count; ++j3 )
761  {
762  double x3 = this->gauss[j3][1];
763  double w3 = this->gauss[j3][0];
764  CartVect x( x1, x2, x3 );
765  I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * w2 * w3 * this->det_jacobian( x );
766  }
767  }
768  }
769  return I;
770  } // LinearHex::integrate_scalar_field()
771 
772  bool LinearHex::inside_nat_space( const CartVect& xi, double& tol ) const
773  {
774  // just look at the box+tol here
775  return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) &&
776  ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol );
777  }
778 
779  // those are not just the corners, but for simplicity, keep this name
780  //
781  const int QuadraticHex::corner[27][3] = {
782  { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, // corner nodes: 0-7
783  { -1, 1, -1 }, // mid-edge nodes: 8-19
784  { -1, -1, 1 }, // center-face nodes 20-25 center node 26
785  { 1, -1, 1 }, //
786  { 1, 1, 1 }, { -1, 1, 1 }, // 4 ----- 19 ----- 7
787  { 0, -1, -1 }, // . | . |
788  { 1, 0, -1 }, // 16 25 18 |
789  { 0, 1, -1 }, // . | . |
790  { -1, 0, -1 }, // 5 ----- 17 ----- 6 |
791  { -1, -1, 0 }, // | 12 | 23 15
792  { 1, -1, 0 }, // | | |
793  { 1, 1, 0 }, // | 20 | 26 | 22 |
794  { -1, 1, 0 }, // | | |
795  { 0, -1, 1 }, // 13 21 | 14 |
796  { 1, 0, 1 }, // | 0 ----- 11 ----- 3
797  { 0, 1, 1 }, // | . | .
798  { -1, 0, 1 }, // | 8 24 | 10
799  { 0, -1, 0 }, // | . | .
800  { 1, 0, 0 }, // 1 ----- 9 ----- 2
801  { 0, 1, 0 }, //
802  { -1, 0, 0 }, { 0, 0, -1 }, { 0, 0, 1 }, { 0, 0, 0 } };
803  // QuadraticHex::QuadraticHex(const std::vector<CartVect>& vertices) : Map(vertices){};
805 
807  double SH( const int i, const double xi )
808  {
809  switch( i )
810  {
811  case -1:
812  return ( xi * xi - xi ) / 2;
813  case 0:
814  return 1 - xi * xi;
815  case 1:
816  return ( xi * xi + xi ) / 2;
817  default:
818  return 0.;
819  }
820  }
821  double DSH( const int i, const double xi )
822  {
823  switch( i )
824  {
825  case -1:
826  return xi - 0.5;
827  case 0:
828  return -2 * xi;
829  case 1:
830  return xi + 0.5;
831  default:
832  return 0.;
833  }
834  }
835 
837  {
838 
839  CartVect x( 0.0 );
840  for( int i = 0; i < 27; i++ )
841  {
842  const double sh = SH( corner[i][0], xi[0] ) * SH( corner[i][1], xi[1] ) * SH( corner[i][2], xi[2] );
843  x += sh * vertex[i];
844  }
845 
846  return x;
847  }
848 
849  bool QuadraticHex::inside_nat_space( const CartVect& xi, double& tol ) const
850  { // just look at the box+tol here
851  return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) &&
852  ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol );
853  }
854 
856  {
857  Matrix3 J( 0.0 );
858  for( int i = 0; i < 27; i++ )
859  {
860  const double sh[3] = { SH( corner[i][0], xi[0] ), SH( corner[i][1], xi[1] ), SH( corner[i][2], xi[2] ) };
861  const double dsh[3] = { DSH( corner[i][0], xi[0] ), DSH( corner[i][1], xi[1] ),
862  DSH( corner[i][2], xi[2] ) };
863 
864  for( int j = 0; j < 3; j++ )
865  {
866  J( j, 0 ) += dsh[0] * sh[1] * sh[2] * vertex[i][j]; // dxj/dr first column
867  J( j, 1 ) += sh[0] * dsh[1] * sh[2] * vertex[i][j]; // dxj/ds
868  J( j, 2 ) += sh[0] * sh[1] * dsh[2] * vertex[i][j]; // dxj/dt
869  }
870  }
871 
872  return J;
873  }
874  double QuadraticHex::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const
875  {
876  double x = 0.0;
877  for( int i = 0; i < 27; i++ )
878  {
879  const double sh = SH( corner[i][0], xi[0] ) * SH( corner[i][1], xi[1] ) * SH( corner[i][2], xi[2] );
880  x += sh * field_vertex_values[i];
881  }
882 
883  return x;
884  }
885  double QuadraticHex::integrate_scalar_field( const double* /*field_vertex_values*/ ) const
886  {
887  return 0.; // TODO: gaussian integration , probably 2x2x2
888  }
889 
890  const double LinearTet::corner[4][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
891 
892  LinearTet::LinearTet() : Map( 0 ), det_T( 0.0 ), det_T_inverse( 0.0 ) {} // LinearTet::LinearTet()
893 
895 
896  void LinearTet::set_vertices( const std::vector< CartVect >& v )
897  {
898  this->Map::set_vertices( v );
899  this->T =
900  Matrix3( v[1][0] - v[0][0], v[2][0] - v[0][0], v[3][0] - v[0][0], v[1][1] - v[0][1], v[2][1] - v[0][1],
901  v[3][1] - v[0][1], v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] );
902  this->T_inverse = this->T.inverse();
903  this->det_T = this->T.determinant();
904  this->det_T_inverse = ( this->det_T < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / this->det_T );
905  } // LinearTet::set_vertices()
906 
907  double LinearTet::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
908  {
909  double f0 = field_vertex_value[0];
910  double f = f0;
911  for( unsigned i = 1; i < 4; ++i )
912  {
913  f += ( field_vertex_value[i] - f0 ) * xi[i - 1];
914  }
915  return f;
916  } // LinearTet::evaluate_scalar_field()
917 
918  CartVect LinearTet::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const
919  {
920  return this->T_inverse * ( x - this->vertex[0] );
921  } // LinearTet::ievaluate
922 
923  double LinearTet::integrate_scalar_field( const double* field_vertex_values ) const
924  {
925  double I( 0.0 );
926  for( unsigned int i = 0; i < 4; ++i )
927  {
928  I += field_vertex_values[i];
929  }
930  I *= this->det_T / 24.0;
931  return I;
932  } // LinearTet::integrate_scalar_field()
933  bool LinearTet::inside_nat_space( const CartVect& xi, double& tol ) const
934  {
935  // linear tet space is a tetra with vertices (0,0,0), (1,0,0), (0,1,0), (0, 0, 1)
936  // first check if outside bigger box, then below the plane x+y+z=1
937  return ( xi[0] >= -tol ) && ( xi[1] >= -tol ) && ( xi[2] >= -tol ) && ( xi[0] + xi[1] + xi[2] < 1.0 + tol );
938  }
939  // SpectralHex
940 
941  // filescope for static member data that is cached
942  int SpectralHex::_n;
947 
948  bool SpectralHex::_init = false;
949 
951  {
952  _xyz[0] = _xyz[1] = _xyz[2] = NULL;
953  }
954  // the preferred constructor takes pointers to GL blocked positions
955  SpectralHex::SpectralHex( int order, double* x, double* y, double* z ) : Map( 0 )
956  {
957  Init( order );
958  _xyz[0] = x;
959  _xyz[1] = y;
960  _xyz[2] = z;
961  }
962  SpectralHex::SpectralHex( int order ) : Map( 0 )
963  {
964  Init( order );
965  _xyz[0] = _xyz[1] = _xyz[2] = NULL;
966  }
968  {
969  if( _init ) freedata();
970  _init = false;
971  }
972  void SpectralHex::Init( int order )
973  {
974  if( _init && _n == order ) return;
975  if( _init && _n != order )
976  {
977  // TODO: free data cached
978  freedata();
979  }
980  // compute stuff that depends only on order
981  _init = true;
982  _n = order;
983  // triplicates! n is the same in all directions !!!
984  for( int d = 0; d < 3; d++ )
985  {
986  _z[d] = tmalloc( realType, _n );
987  lobatto_nodes( _z[d], _n );
988  lagrange_setup( &_ld[d], _z[d], _n );
989  }
990  opt_alloc_3( &_data, _ld );
991 
992  unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
993  _odwork = tmalloc( realType, 6 * nf + 9 * ne + nw );
994  }
995  void SpectralHex::freedata()
996  {
997  for( int d = 0; d < 3; d++ )
998  {
999  free( _z[d] );
1000  lagrange_free( &_ld[d] );
1001  }
1002  opt_free_3( &_data );
1003  free( _odwork );
1004  }
1005 
1006  void SpectralHex::set_gl_points( double* x, double* y, double* z )
1007  {
1008  _xyz[0] = x;
1009  _xyz[1] = y;
1010  _xyz[2] = z;
1011  }
1012  CartVect SpectralHex::evaluate( const CartVect& xi ) const
1013  {
1014  // piece that we shouldn't want to cache
1015  int d = 0;
1016  for( d = 0; d < 3; d++ )
1017  {
1018  lagrange_0( &_ld[d], xi[d] );
1019  }
1020  CartVect result;
1021  for( d = 0; d < 3; d++ )
1022  {
1023  result[d] = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n,
1024  _xyz[d], // this is the "field"
1025  _odwork );
1026  }
1027  return result;
1028  }
1029  // replicate the functionality of hex_findpt
1030  CartVect SpectralHex::ievaluate( CartVect const& xyz, double tol, const CartVect& x0 ) const
1031  {
1032  // find nearest point
1033  realType x_star[3];
1034  xyz.get( x_star );
1035 
1036  realType r[3] = { 0, 0, 0 }; // initial guess for parametric coords
1037  x0.get( r );
1038  unsigned c = opt_no_constraints_3;
1039  realType dist = opt_findpt_3( &_data, (const realType**)_xyz, x_star, r, &c );
1040  // if it did not converge, get out with throw...
1041  if( dist > 10 * tol ) // outside the element
1042  {
1043  std::vector< CartVect > dummy;
1044  throw Map::EvaluationError( xyz, dummy );
1045  }
1046  // c tells us if we landed inside the element or exactly on a face, edge, or node
1047  // also, dist shows the distance to the computed point.
1048  // copy parametric coords back
1049  return CartVect( r );
1050  }
1051  Matrix3 SpectralHex::jacobian( const CartVect& xi ) const
1052  {
1053  realType x_i[3];
1054  xi.get( x_i );
1055  // set the positions of GL nodes, before evaluations
1056  _data.elx[0] = _xyz[0];
1057  _data.elx[1] = _xyz[1];
1058  _data.elx[2] = _xyz[2];
1059  opt_vol_set_intp_3( &_data, x_i );
1060  Matrix3 J( 0. );
1061  // it is organized differently
1062  J( 0, 0 ) = _data.jac[0]; // dx/dr
1063  J( 0, 1 ) = _data.jac[1]; // dx/ds
1064  J( 0, 2 ) = _data.jac[2]; // dx/dt
1065  J( 1, 0 ) = _data.jac[3]; // dy/dr
1066  J( 1, 1 ) = _data.jac[4]; // dy/ds
1067  J( 1, 2 ) = _data.jac[5]; // dy/dt
1068  J( 2, 0 ) = _data.jac[6]; // dz/dr
1069  J( 2, 1 ) = _data.jac[7]; // dz/ds
1070  J( 2, 2 ) = _data.jac[8]; // dz/dt
1071  return J;
1072  }
1073  double SpectralHex::evaluate_scalar_field( const CartVect& xi, const double* field ) const
1074  {
1075  // piece that we shouldn't want to cache
1076  int d;
1077  for( d = 0; d < 3; d++ )
1078  {
1079  lagrange_0( &_ld[d], xi[d] );
1080  }
1081 
1082  double value = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork );
1083  return value;
1084  }
1085  double SpectralHex::integrate_scalar_field( const double* field_vertex_values ) const
1086  {
1087  // set the position of GL points
1088  // set the positions of GL nodes, before evaluations
1089  _data.elx[0] = _xyz[0];
1090  _data.elx[1] = _xyz[1];
1091  _data.elx[2] = _xyz[2];
1092  double xi[3];
1093  // triple loop; the most inner loop is in r direction, then s, then t
1094  double integral = 0.;
1095  // double volume = 0;
1096  int index = 0; // used fr the inner loop
1097  for( int k = 0; k < _n; k++ )
1098  {
1099  xi[2] = _ld[2].z[k];
1100  // double wk= _ld[2].w[k];
1101  for( int j = 0; j < _n; j++ )
1102  {
1103  xi[1] = _ld[1].z[j];
1104  // double wj= _ld[1].w[j];
1105  for( int i = 0; i < _n; i++ )
1106  {
1107  xi[0] = _ld[0].z[i];
1108  // double wi= _ld[0].w[i];
1109  opt_vol_set_intp_3( &_data, xi );
1110  double wk = _ld[2].J[k];
1111  double wj = _ld[1].J[j];
1112  double wi = _ld[0].J[i];
1113  Matrix3 J( 0. );
1114  // it is organized differently
1115  J( 0, 0 ) = _data.jac[0]; // dx/dr
1116  J( 0, 1 ) = _data.jac[1]; // dx/ds
1117  J( 0, 2 ) = _data.jac[2]; // dx/dt
1118  J( 1, 0 ) = _data.jac[3]; // dy/dr
1119  J( 1, 1 ) = _data.jac[4]; // dy/ds
1120  J( 1, 2 ) = _data.jac[5]; // dy/dt
1121  J( 2, 0 ) = _data.jac[6]; // dz/dr
1122  J( 2, 1 ) = _data.jac[7]; // dz/ds
1123  J( 2, 2 ) = _data.jac[8]; // dz/dt
1124  double bm = wk * wj * wi * J.determinant();
1125  integral += bm * field_vertex_values[index++];
1126  // volume +=bm;
1127  }
1128  }
1129  }
1130  // std::cout << "volume: " << volume << "\n";
1131  return integral;
1132  }
1133  // this is the same as a linear hex, although we should not need it
1134  bool SpectralHex::inside_nat_space( const CartVect& xi, double& tol ) const
1135  {
1136  // just look at the box+tol here
1137  return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) &&
1138  ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol );
1139  }
1140 
1141  // SpectralHex
1142 
1143  // filescope for static member data that is cached
1144  const double LinearQuad::corner[4][3] = { { -1, -1, 0 }, { 1, -1, 0 }, { 1, 1, 0 }, { -1, 1, 0 } };
1145 
1146  LinearQuad::LinearQuad() : Map( 0 ) {} // LinearQuad::LinearQuad()
1147 
1149 
1150  /* For each point, its weight and location are stored as an array.
1151  Hence, the inner dimension is 2, the outer dimension is gauss_count.
1152  We use a one-point Gaussian quadrature, since it integrates linear functions exactly.
1153  */
1154  const double LinearQuad::gauss[1][2] = { { 2.0, 0.0 } };
1155 
1157  {
1158  CartVect x( 0.0 );
1159  for( unsigned i = 0; i < LinearQuad::corner_count; ++i )
1160  {
1161  const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] );
1162  x += N_i * this->vertex[i];
1163  }
1165  return x;
1166  } // LinearQuad::evaluate
1167 
1169  {
1170  // this basically ignores the z component: xi[2] or vertex[][2]
1171  Matrix3 J( 0.0 );
1172  for( unsigned i = 0; i < LinearQuad::corner_count; ++i )
1173  {
1174  const double xi_p = 1 + xi[0] * corner[i][0];
1175  const double eta_p = 1 + xi[1] * corner[i][1];
1176  const double dNi_dxi = corner[i][0] * eta_p;
1177  const double dNi_deta = corner[i][1] * xi_p;
1178  J( 0, 0 ) += dNi_dxi * vertex[i][0];
1179  J( 1, 0 ) += dNi_dxi * vertex[i][1];
1180  J( 0, 1 ) += dNi_deta * vertex[i][0];
1181  J( 1, 1 ) += dNi_deta * vertex[i][1];
1182  }
1183  J( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */
1185  return J;
1186  } // LinearQuad::jacobian()
1187 
1188  double LinearQuad::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
1189  {
1190  double f( 0.0 );
1191  for( unsigned i = 0; i < LinearQuad::corner_count; ++i )
1192  {
1193  const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] );
1194  f += N_i * field_vertex_value[i];
1195  }
1197  return f;
1198  } // LinearQuad::evaluate_scalar_field()
1199 
1200  double LinearQuad::integrate_scalar_field( const double* field_vertex_values ) const
1201  {
1202  double I( 0.0 );
1203  for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 )
1204  {
1205  double x1 = this->gauss[j1][1];
1206  double w1 = this->gauss[j1][0];
1207  for( unsigned int j2 = 0; j2 < this->gauss_count; ++j2 )
1208  {
1209  double x2 = this->gauss[j2][1];
1210  double w2 = this->gauss[j2][0];
1211  CartVect x( x1, x2, 0.0 );
1212  I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * w2 * this->det_jacobian( x );
1213  }
1214  }
1215  return I;
1216  } // LinearQuad::integrate_scalar_field()
1217 
1218  bool LinearQuad::inside_nat_space( const CartVect& xi, double& tol ) const
1219  {
1220  // just look at the box+tol here
1221  return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol );
1222  }
1223 
1224  // filescope for static member data that is cached
1225  int SpectralQuad::_n;
1231  bool SpectralQuad::_init = false;
1232 
1234  {
1235  _xyz[0] = _xyz[1] = _xyz[2] = NULL;
1236  }
1237  // the preferred constructor takes pointers to GL blocked positions
1238  SpectralQuad::SpectralQuad( int order, double* x, double* y, double* z ) : Map( 0 )
1239  {
1240  Init( order );
1241  _xyz[0] = x;
1242  _xyz[1] = y;
1243  _xyz[2] = z;
1244  }
1245  SpectralQuad::SpectralQuad( int order ) : Map( 4 )
1246  {
1247  Init( order );
1248  _xyz[0] = _xyz[1] = _xyz[2] = NULL;
1249  }
1251  {
1252  if( _init ) freedata();
1253  _init = false;
1254  }
1255  void SpectralQuad::Init( int order )
1256  {
1257  if( _init && _n == order ) return;
1258  if( _init && _n != order )
1259  {
1260  // TODO: free data cached
1261  freedata();
1262  }
1263  // compute stuff that depends only on order
1264  _init = true;
1265  _n = order;
1266  // duplicates! n is the same in all directions !!!
1267  for( int d = 0; d < 2; d++ )
1268  {
1269  _z[d] = tmalloc( realType, _n );
1270  lobatto_nodes( _z[d], _n );
1271  lagrange_setup( &_ld[d], _z[d], _n );
1272  }
1273  opt_alloc_2( &_data, _ld );
1274 
1275  unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
1276  _odwork = tmalloc( realType, 6 * nf + 9 * ne + nw );
1277  _glpoints = tmalloc( realType, 3 * nf );
1278  }
1279 
1280  void SpectralQuad::freedata()
1281  {
1282  for( int d = 0; d < 2; d++ )
1283  {
1284  free( _z[d] );
1285  lagrange_free( &_ld[d] );
1286  }
1287  opt_free_2( &_data );
1288  free( _odwork );
1289  free( _glpoints );
1290  }
1291 
1292  void SpectralQuad::set_gl_points( double* x, double* y, double* z )
1293  {
1294  _xyz[0] = x;
1295  _xyz[1] = y;
1296  _xyz[2] = z;
1297  }
1299  {
1300  // piece that we shouldn't want to cache
1301  int d = 0;
1302  for( d = 0; d < 2; d++ )
1303  {
1304  lagrange_0( &_ld[d], xi[d] );
1305  }
1306  CartVect result;
1307  for( d = 0; d < 3; d++ )
1308  {
1309  result[d] = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _xyz[d], _odwork );
1310  }
1311  return result;
1312  }
1313  // replicate the functionality of hex_findpt
1314  CartVect SpectralQuad::ievaluate( CartVect const& xyz, double /*tol*/, const CartVect& /*x0*/ ) const
1315  {
1316  // find nearest point
1317  realType x_star[3];
1318  xyz.get( x_star );
1319 
1320  realType r[2] = { 0, 0 }; // initial guess for parametric coords
1321  unsigned c = opt_no_constraints_3;
1322  realType dist = opt_findpt_2( &_data, (const realType**)_xyz, x_star, r, &c );
1323  // if it did not converge, get out with throw...
1324  if( dist > 0.9e+30 )
1325  {
1326  std::vector< CartVect > dummy;
1327  throw Map::EvaluationError( xyz, dummy );
1328  }
1329 
1330  // c tells us if we landed inside the element or exactly on a face, edge, or node
1331  // also, dist shows the distance to the computed point.
1332  // copy parametric coords back
1333  return CartVect( r[0], r[1], 0. );
1334  }
1335 
1336  Matrix3 SpectralQuad::jacobian( const CartVect& /*xi*/ ) const
1337  {
1338  // not implemented
1339  Matrix3 J( 0. );
1340  return J;
1341  }
1342 
1343  double SpectralQuad::evaluate_scalar_field( const CartVect& xi, const double* field ) const
1344  {
1345  // piece that we shouldn't want to cache
1346  int d;
1347  for( d = 0; d < 2; d++ )
1348  {
1349  lagrange_0( &_ld[d], xi[d] );
1350  }
1351 
1352  double value = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, field, _odwork );
1353  return value;
1354  }
1355  double SpectralQuad::integrate_scalar_field( const double* /*field_vertex_values*/ ) const
1356  {
1357  return 0.; // not implemented
1358  }
1359  // this is the same as a linear hex, although we should not need it
1360  bool SpectralQuad::inside_nat_space( const CartVect& xi, double& tol ) const
1361  {
1362  // just look at the box+tol here
1363  return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol );
1364  }
1365  // something we don't do for spectral hex; we do it here because
1366  // we do not store the position of gl points in a tag yet
1368  {
1369  // will need to use shape functions on a simple linear quad to compute gl points
1370  // so we know the position of gl points in parametric space, we will just compute those
1371  // from the 3d vertex position (corner nodes of the quad), using simple mapping
1372  assert( this->vertex.size() == 4 );
1373  static double corner_xi[4][2] = { { -1., -1. }, { 1., -1. }, { 1., 1. }, { -1., 1. } };
1374  // we will use the cached lobatto nodes in parametric space _z[d] (the same in both
1375  // directions)
1376  int indexGL = 0;
1377  int n2 = _n * _n;
1378  for( int i = 0; i < _n; i++ )
1379  {
1380  double eta = _z[0][i];
1381  for( int j = 0; j < _n; j++ )
1382  {
1383  double csi = _z[1][j]; // we could really use the same _z[0] array of lobatto nodes
1384  CartVect pos( 0.0 );
1385  for( int k = 0; k < 4; k++ )
1386  {
1387  const double N_k = ( 1 + csi * corner_xi[k][0] ) * ( 1 + eta * corner_xi[k][1] );
1388  pos += N_k * vertex[k];
1389  }
1390  pos *= 0.25; // these are x, y, z of gl points; reorder them
1391  _glpoints[indexGL] = pos[0]; // x
1392  _glpoints[indexGL + n2] = pos[1]; // y
1393  _glpoints[indexGL + 2 * n2] = pos[2]; // z
1394  indexGL++;
1395  }
1396  }
1397  // now, we can set the _xyz pointers to internal memory allocated to these!
1398  _xyz[0] = &( _glpoints[0] );
1399  _xyz[1] = &( _glpoints[n2] );
1400  _xyz[2] = &( _glpoints[2 * n2] );
1401  }
1402  void SpectralQuad::get_gl_points( double*& x, double*& y, double*& z, int& psize )
1403  {
1404  x = (double*)_xyz[0];
1405  y = (double*)_xyz[1];
1406  z = (double*)_xyz[2];
1407  psize = _n * _n;
1408  }
1409 } // namespace Element
1410 
1411 } // namespace moab