Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
Classes | Functions | Variables
moab::ElemUtil Namespace Reference

Classes

class  VolMap
 Class representing a 3-D mapping function (e.g. shape function for volume element) More...
 
class  LinearHexMap
 Shape function for trilinear hexahedron. More...
 

Functions

bool nat_coords_trilinear_hex (const CartVect *corner_coords, const CartVect &x, CartVect &xi, double tol)
 
void nat_coords_trilinear_hex2 (const CartVect hex[8], const CartVect &xyz, CartVect &ncoords, double etol)
 
bool point_in_trilinear_hex (const CartVect *hex, const CartVect &xyz, double etol)
 
bool point_in_trilinear_hex (const CartVect *hex, const CartVect &xyz, const CartVect &box_min, const CartVect &box_max, double etol)
 
void hex_findpt (realType *xm[3], int n, CartVect xyz, CartVect &rst, double &dist)
 
void hex_eval (realType *field, int n, CartVect rstCartVec, double &value)
 
bool integrate_trilinear_hex (const CartVect *hex_corners, double *corner_fields, double &field_val, int num_pts)
 
void nat_coords_trilinear_hex2 (const CartVect *hex_corners, const CartVect &x, CartVect &xi, double til)
 

Variables

bool debug = false
 
const double gauss_1 [1][2] = { { 2.0, 0.0 } }
 
const double gauss_2 [2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } }
 
const double gauss_3 [3][2]
 
const double gauss_4 [4][2]
 
const double gauss_5 [5][2]
 

Function Documentation

◆ hex_eval()

void moab::ElemUtil::hex_eval ( realType *  field,
int  n,
CartVect  rstCartVec,
double &  value 
)

Definition at line 246 of file ElemUtil.cpp.

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  }

References moab::CartVect::get(), lagrange_0(), lagrange_free(), lagrange_setup(), lobatto_nodes(), tensor_i3(), and tmalloc.

◆ hex_findpt()

void moab::ElemUtil::hex_findpt ( realType *  xm[3],
int  n,
CartVect  xyz,
CartVect &  rst,
double &  dist 
)

Definition at line 201 of file ElemUtil.cpp.

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  }

References moab::CartVect::get(), lagrange_free(), lagrange_setup(), lobatto_nodes(), opt_alloc_3(), opt_findpt_3(), opt_free_3(), opt_no_constraints_3, and tmalloc.

Referenced by nat_coords_trilinear_hex2().

◆ integrate_trilinear_hex()

bool moab::ElemUtil::integrate_trilinear_hex ( const CartVect *  hex_corners,
double *  corner_fields,
double &  field_val,
int  num_pts 
)

Definition at line 312 of file ElemUtil.cpp.

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  }

References debug, moab::Matrix3::determinant(), moab::ElemUtil::LinearHexMap::evaluate_scalar_field(), gauss_1, gauss_2, gauss_3, gauss_4, gauss_5, hex_corners, and moab::ElemUtil::LinearHexMap::jacobian().

◆ nat_coords_trilinear_hex()

bool moab::ElemUtil::nat_coords_trilinear_hex ( const CartVect *  corner_coords,
const CartVect &  x,
CartVect &  xi,
double  tol 
)

Definition at line 119 of file ElemUtil.cpp.

120  {
121  return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
122  }

References moab::ElemUtil::VolMap::solve_inverse().

Referenced by point_in_trilinear_hex().

◆ nat_coords_trilinear_hex2() [1/2]

void moab::ElemUtil::nat_coords_trilinear_hex2 ( const CartVect *  hex_corners,
const CartVect &  x,
CartVect &  xi,
double  til 
)

◆ nat_coords_trilinear_hex2() [2/2]

void moab::ElemUtil::nat_coords_trilinear_hex2 ( const CartVect  hex[8],
const CartVect &  xyz,
CartVect &  ncoords,
double  etol 
)

Definition at line 127 of file ElemUtil.cpp.

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  }

References moab::CartVect::get(), hex_findpt(), and MOAB_POLY_EPS.

◆ point_in_trilinear_hex() [1/2]

bool moab::ElemUtil::point_in_trilinear_hex ( const CartVect *  hex,
const CartVect &  xyz,
const CartVect &  box_min,
const CartVect &  box_max,
double  etol 
)

Definition at line 169 of file ElemUtil.cpp.

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  }

References box_max(), box_min(), and point_in_trilinear_hex().

◆ point_in_trilinear_hex() [2/2]

bool moab::ElemUtil::point_in_trilinear_hex ( const CartVect *  hex,
const CartVect &  xyz,
double  etol 
)

Definition at line 162 of file ElemUtil.cpp.

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  }

References nat_coords_trilinear_hex().

Referenced by main(), and point_in_trilinear_hex().

Variable Documentation

◆ debug

bool moab::ElemUtil::debug = false

Definition at line 13 of file ElemUtil.cpp.

Referenced by integrate_trilinear_hex().

◆ gauss_1

const double moab::ElemUtil::gauss_1[1][2] = { { 2.0, 0.0 } }

Definition at line 290 of file ElemUtil.cpp.

Referenced by integrate_trilinear_hex().

◆ gauss_2

const double moab::ElemUtil::gauss_2[2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } }

Definition at line 292 of file ElemUtil.cpp.

Referenced by integrate_trilinear_hex().

◆ gauss_3

const double moab::ElemUtil::gauss_3[3][2]
Initial value:
= { { 0.5555555555, -0.7745966692 },
{ 0.8888888888, 0.0 },
{ 0.5555555555, 0.7745966692 } }

Definition at line 294 of file ElemUtil.cpp.

Referenced by integrate_trilinear_hex().

◆ gauss_4

const double moab::ElemUtil::gauss_4[4][2]
Initial value:
= { { 0.3478548451, -0.8611363116 },
{ 0.6521451549, -0.3399810436 },
{ 0.6521451549, 0.3399810436 },
{ 0.3478548451, 0.8611363116 } }

Definition at line 298 of file ElemUtil.cpp.

Referenced by integrate_trilinear_hex().

◆ gauss_5

const double moab::ElemUtil::gauss_5[5][2]
Initial value:
= { { 0.2369268851, -0.9061798459 },
{ 0.4786286705, -0.5384693101 },
{ 0.5688888889, 0.0 },
{ 0.4786286705, 0.5384693101 },
{ 0.2369268851, 0.9061798459 } }

Definition at line 303 of file ElemUtil.cpp.

Referenced by integrate_trilinear_hex().