MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include <iostream> 00002 #include <limits> 00003 #include <assert.h> 00004 00005 #include "ElemUtil.hpp" 00006 #include "moab/BoundBox.hpp" 00007 00008 namespace moab 00009 { 00010 namespace ElemUtil 00011 { 00012 00013 bool debug = false; 00014 00015 /**\brief Class representing a 3-D mapping function (e.g. shape function for volume element) */ 00016 class VolMap 00017 { 00018 public: 00019 /**\brief Return \f$\vec \xi\f$ corresponding to logical center of element */ 00020 virtual CartVect center_xi() const = 0; 00021 /**\brief Evaluate mapping function (calculate \f$\vec x = F($\vec \xi)\f$ )*/ 00022 virtual CartVect evaluate( const CartVect& xi ) const = 0; 00023 /**\brief Evaluate Jacobian of mapping function */ 00024 virtual Matrix3 jacobian( const CartVect& xi ) const = 0; 00025 /**\brief Evaluate inverse of mapping function (calculate \f$\vec \xi = F^-1($\vec x)\f$ )*/ 00026 bool solve_inverse( const CartVect& x, CartVect& xi, double tol ) const; 00027 }; 00028 00029 bool VolMap::solve_inverse( const CartVect& x, CartVect& xi, double tol ) const 00030 { 00031 const double error_tol_sqr = tol * tol; 00032 double det; 00033 xi = center_xi(); 00034 CartVect delta = evaluate( xi ) - x; 00035 Matrix3 J; 00036 while( delta % delta > error_tol_sqr ) 00037 { 00038 J = jacobian( xi ); 00039 det = J.determinant(); 00040 if( det < std::numeric_limits< double >::epsilon() ) return false; 00041 xi -= J.inverse() * delta; 00042 delta = evaluate( xi ) - x; 00043 } 00044 return true; 00045 } 00046 00047 /**\brief Shape function for trilinear hexahedron */ 00048 class LinearHexMap : public VolMap 00049 { 00050 public: 00051 LinearHexMap( const CartVect* corner_coords ) : corners( corner_coords ) {} 00052 virtual CartVect center_xi() const; 00053 virtual CartVect evaluate( const CartVect& xi ) const; 00054 virtual double evaluate_scalar_field( const CartVect& xi, const double* f_vals ) const; 00055 virtual Matrix3 jacobian( const CartVect& xi ) const; 00056 00057 private: 00058 const CartVect* corners; 00059 static const double corner_xi[8][3]; 00060 }; 00061 00062 const double LinearHexMap::corner_xi[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 }, 00063 { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } }; 00064 CartVect LinearHexMap::center_xi() const 00065 { 00066 return CartVect( 0.0 ); 00067 } 00068 00069 CartVect LinearHexMap::evaluate( const CartVect& xi ) const 00070 { 00071 CartVect x( 0.0 ); 00072 for( unsigned i = 0; i < 8; ++i ) 00073 { 00074 const double N_i = 00075 ( 1 + xi[0] * corner_xi[i][0] ) * ( 1 + xi[1] * corner_xi[i][1] ) * ( 1 + xi[2] * corner_xi[i][2] ); 00076 x += N_i * corners[i]; 00077 } 00078 x *= 0.125; 00079 return x; 00080 } 00081 00082 double LinearHexMap::evaluate_scalar_field( const CartVect& xi, const double* f_vals ) const 00083 { 00084 double f( 0.0 ); 00085 for( unsigned i = 0; i < 8; ++i ) 00086 { 00087 const double N_i = 00088 ( 1 + xi[0] * corner_xi[i][0] ) * ( 1 + xi[1] * corner_xi[i][1] ) * ( 1 + xi[2] * corner_xi[i][2] ); 00089 f += N_i * f_vals[i]; 00090 } 00091 f *= 0.125; 00092 return f; 00093 } 00094 00095 Matrix3 LinearHexMap::jacobian( const CartVect& xi ) const 00096 { 00097 Matrix3 J( 0.0 ); 00098 for( unsigned i = 0; i < 8; ++i ) 00099 { 00100 const double xi_p = 1 + xi[0] * corner_xi[i][0]; 00101 const double eta_p = 1 + xi[1] * corner_xi[i][1]; 00102 const double zeta_p = 1 + xi[2] * corner_xi[i][2]; 00103 const double dNi_dxi = corner_xi[i][0] * eta_p * zeta_p; 00104 const double dNi_deta = corner_xi[i][1] * xi_p * zeta_p; 00105 const double dNi_dzeta = corner_xi[i][2] * xi_p * eta_p; 00106 J( 0, 0 ) += dNi_dxi * corners[i][0]; 00107 J( 1, 0 ) += dNi_dxi * corners[i][1]; 00108 J( 2, 0 ) += dNi_dxi * corners[i][2]; 00109 J( 0, 1 ) += dNi_deta * corners[i][0]; 00110 J( 1, 1 ) += dNi_deta * corners[i][1]; 00111 J( 2, 1 ) += dNi_deta * corners[i][2]; 00112 J( 0, 2 ) += dNi_dzeta * corners[i][0]; 00113 J( 1, 2 ) += dNi_dzeta * corners[i][1]; 00114 J( 2, 2 ) += dNi_dzeta * corners[i][2]; 00115 } 00116 return J *= 0.125; 00117 } 00118 00119 bool nat_coords_trilinear_hex( const CartVect* corner_coords, const CartVect& x, CartVect& xi, double tol ) 00120 { 00121 return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol ); 00122 } 00123 // 00124 // nat_coords_trilinear_hex2 00125 // Duplicate functionality of nat_coords_trilinear_hex using hex_findpt 00126 // 00127 void nat_coords_trilinear_hex2( const CartVect hex[8], const CartVect& xyz, CartVect& ncoords, double etol ) 00128 00129 { 00130 const int ndim = 3; 00131 const int nverts = 8; 00132 const int vertMap[nverts] = { 0, 1, 3, 2, 4, 5, 7, 6 }; // Map from nat to lex ordering 00133 00134 const int n = 2; // linear 00135 realType coords[ndim * nverts]; // buffer 00136 00137 realType* xm[ndim]; 00138 for( int i = 0; i < ndim; i++ ) 00139 xm[i] = coords + i * nverts; 00140 00141 // stuff hex into coords 00142 for( int i = 0; i < nverts; i++ ) 00143 { 00144 realType vcoord[ndim]; 00145 hex[i].get( vcoord ); 00146 00147 for( int d = 0; d < ndim; d++ ) 00148 coords[d * nverts + vertMap[i]] = vcoord[d]; 00149 } 00150 00151 double dist = 0.0; 00152 ElemUtil::hex_findpt( xm, n, xyz, ncoords, dist ); 00153 if( 3 * MOAB_POLY_EPS < dist ) 00154 { 00155 // outside element, set extremal values to something outside range 00156 for( int j = 0; j < 3; j++ ) 00157 { 00158 if( ncoords[j] < ( -1.0 - etol ) || ncoords[j] > ( 1.0 + etol ) ) ncoords[j] *= 10; 00159 } 00160 } 00161 } 00162 bool point_in_trilinear_hex( const CartVect* hex, const CartVect& xyz, double etol ) 00163 { 00164 CartVect xi; 00165 return nat_coords_trilinear_hex( hex, xyz, xi, etol ) && ( fabs( xi[0] - 1.0 ) < etol ) && 00166 ( fabs( xi[1] - 1.0 ) < etol ) && ( fabs( xi[2] - 1.0 ) < etol ); 00167 } 00168 00169 bool point_in_trilinear_hex( const CartVect* hex, const CartVect& xyz, const CartVect& box_min, 00170 const CartVect& box_max, double etol ) 00171 { 00172 // all values scaled by 2 (eliminates 3 flops) 00173 const CartVect mid = box_max + box_min; 00174 const CartVect dim = box_max - box_min; 00175 const CartVect pt = 2 * xyz - mid; 00176 return ( fabs( pt[0] - dim[0] ) < etol ) && ( fabs( pt[1] - dim[1] ) < etol ) && 00177 ( fabs( pt[2] - dim[2] ) < etol ) && point_in_trilinear_hex( hex, xyz, etol ); 00178 } 00179 00180 // Wrapper to James Lottes' findpt routines 00181 // hex_findpt 00182 // Find the parametric coordinates of an xyz point inside 00183 // a 3d hex spectral element with n nodes per dimension 00184 // xm: coordinates fields, value of x,y,z for each of then n*n*n gauss-lobatto nodes. Nodes are 00185 // in lexicographical order (x is fastest-changing) n: number of nodes per dimension -- n=2 for 00186 // a linear element xyz: input, point to find rst: output: parametric coords of xyz inside the 00187 // element. If xyz is outside the element, rst will be the coords of the closest point dist: 00188 // output: distance between xyz and the point with parametric coords rst 00189 /*extern "C"{ 00190 #include "types.h" 00191 #include "poly.h" 00192 #include "tensor.h" 00193 #include "findpt.h" 00194 #include "extrafindpt.h" 00195 #include "errmem.h" 00196 }*/ 00197 00198 void hex_findpt( realType* xm[3], int n, CartVect xyz, CartVect& rst, double& dist ) 00199 { 00200 00201 // compute stuff that only depends on the order -- could be cached 00202 realType* z[3]; 00203 lagrange_data ld[3]; 00204 opt_data_3 data; 00205 00206 // triplicates 00207 for( int d = 0; d < 3; d++ ) 00208 { 00209 z[d] = tmalloc( realType, n ); 00210 lobatto_nodes( z[d], n ); 00211 lagrange_setup( &ld[d], z[d], n ); 00212 } 00213 00214 opt_alloc_3( &data, ld ); 00215 00216 // find nearest point 00217 realType x_star[3]; 00218 xyz.get( x_star ); 00219 00220 realType r[3] = { 0, 0, 0 }; // initial guess for parametric coords 00221 unsigned c = opt_no_constraints_3; 00222 dist = opt_findpt_3( &data, (const realType**)xm, x_star, r, &c ); 00223 // c tells us if we landed inside the element or exactly on a face, edge, or node 00224 00225 // copy parametric coords back 00226 rst = r; 00227 00228 // Clean-up (move to destructor if we decide to cache) 00229 opt_free_3( &data ); 00230 for( int d = 0; d < 3; ++d ) 00231 lagrange_free( &ld[d] ); 00232 for( int d = 0; d < 3; ++d ) 00233 free( z[d] ); 00234 } 00235 00236 // hex_eval 00237 // Evaluate a field in a 3d hex spectral element with n nodes per dimension, at some given 00238 // parametric coordinates field: field values for each of then n*n*n gauss-lobatto nodes. Nodes 00239 // are in lexicographical order (x is fastest-changing) n: number of nodes per dimension -- n=2 00240 // for a linear element rst: input: parametric coords of the point where we want to evaluate the 00241 // field value: output: value of field at rst 00242 00243 void hex_eval( realType* field, int n, CartVect rstCartVec, double& value ) 00244 { 00245 int d; 00246 realType rst[3]; 00247 rstCartVec.get( rst ); 00248 00249 // can cache stuff below 00250 lagrange_data ld[3]; 00251 realType* z[3]; 00252 for( d = 0; d < 3; ++d ) 00253 { 00254 z[d] = tmalloc( realType, n ); 00255 lobatto_nodes( z[d], n ); 00256 lagrange_setup( &ld[d], z[d], n ); 00257 } 00258 00259 // cut and paste -- see findpt.c 00260 const unsigned nf = n * n, ne = n, nw = 2 * n * n + 3 * n; 00261 realType* od_work = tmalloc( realType, 6 * nf + 9 * ne + nw ); 00262 00263 // piece that we shouldn't want to cache 00264 for( d = 0; d < 3; d++ ) 00265 { 00266 lagrange_0( &ld[d], rst[d] ); 00267 } 00268 00269 value = tensor_i3( ld[0].J, ld[0].n, ld[1].J, ld[1].n, ld[2].J, ld[2].n, field, od_work ); 00270 00271 // all this could be cached 00272 for( d = 0; d < 3; d++ ) 00273 { 00274 free( z[d] ); 00275 lagrange_free( &ld[d] ); 00276 } 00277 free( od_work ); 00278 } 00279 00280 // Gaussian quadrature points for a trilinear hex element. 00281 // Five 2d arrays are defined. 00282 // One for the single gaussian point solution, 2 point solution, 00283 // 3 point solution, 4 point solution and 5 point solution. 00284 // There are 2 columns, one for Weights and the other for Locations 00285 // Weight Location 00286 00287 const double gauss_1[1][2] = { { 2.0, 0.0 } }; 00288 00289 const double gauss_2[2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } }; 00290 00291 const double gauss_3[3][2] = { { 0.5555555555, -0.7745966692 }, 00292 { 0.8888888888, 0.0 }, 00293 { 0.5555555555, 0.7745966692 } }; 00294 00295 const double gauss_4[4][2] = { { 0.3478548451, -0.8611363116 }, 00296 { 0.6521451549, -0.3399810436 }, 00297 { 0.6521451549, 0.3399810436 }, 00298 { 0.3478548451, 0.8611363116 } }; 00299 00300 const double gauss_5[5][2] = { { 0.2369268851, -0.9061798459 }, 00301 { 0.4786286705, -0.5384693101 }, 00302 { 0.5688888889, 0.0 }, 00303 { 0.4786286705, 0.5384693101 }, 00304 { 0.2369268851, 0.9061798459 } }; 00305 00306 // Function to integrate the field defined by field_fn function 00307 // over the volume of the trilinear hex defined by the hex_corners 00308 00309 bool integrate_trilinear_hex( const CartVect* hex_corners, double* corner_fields, double& field_val, int num_pts ) 00310 { 00311 // Create the LinearHexMap object using the hex_corners array of CartVects 00312 LinearHexMap hex( hex_corners ); 00313 00314 // Use the correct table of points and locations based on the num_pts parameter 00315 const double( *g_pts )[2] = 0; 00316 switch( num_pts ) 00317 { 00318 case 1: 00319 g_pts = gauss_1; 00320 break; 00321 00322 case 2: 00323 g_pts = gauss_2; 00324 break; 00325 00326 case 3: 00327 g_pts = gauss_3; 00328 break; 00329 00330 case 4: 00331 g_pts = gauss_4; 00332 break; 00333 00334 case 5: 00335 g_pts = gauss_5; 00336 break; 00337 00338 default: // value out of range 00339 return false; 00340 } 00341 00342 // Test code - print Gaussian Quadrature data 00343 if( debug ) 00344 { 00345 for( int r = 0; r < num_pts; r++ ) 00346 for( int c = 0; c < 2; c++ ) 00347 std::cout << "g_pts[" << r << "][" << c << "]=" << g_pts[r][c] << std::endl; 00348 } 00349 // End Test code 00350 00351 double soln = 0.0; 00352 00353 for( int i = 0; i < num_pts; i++ ) 00354 { // Loop for xi direction 00355 double w_i = g_pts[i][0]; 00356 double xi_i = g_pts[i][1]; 00357 for( int j = 0; j < num_pts; j++ ) 00358 { // Loop for eta direction 00359 double w_j = g_pts[j][0]; 00360 double eta_j = g_pts[j][1]; 00361 for( int k = 0; k < num_pts; k++ ) 00362 { // Loop for zeta direction 00363 double w_k = g_pts[k][0]; 00364 double zeta_k = g_pts[k][1]; 00365 00366 // Calculate the "realType" space point given the "normal" point 00367 CartVect normal_pt( xi_i, eta_j, zeta_k ); 00368 00369 // Calculate the value of F(x(xi,eta,zeta),y(xi,eta,zeta),z(xi,eta,zeta) 00370 double field = hex.evaluate_scalar_field( normal_pt, corner_fields ); 00371 00372 // Calculate the Jacobian for this "normal" point and its determinant 00373 Matrix3 J = hex.jacobian( normal_pt ); 00374 double det = J.determinant(); 00375 00376 // Calculate integral and update the solution 00377 soln = soln + ( w_i * w_j * w_k * field * det ); 00378 } 00379 } 00380 } 00381 00382 // Set the output parameter 00383 field_val = soln; 00384 00385 return true; 00386 } 00387 00388 } // namespace ElemUtil 00389 00390 namespace Element 00391 { 00392 00393 Map::~Map() {} 00394 00395 inline const std::vector< CartVect >& Map::get_vertices() 00396 { 00397 return this->vertex; 00398 } 00399 // 00400 void Map::set_vertices( const std::vector< CartVect >& v ) 00401 { 00402 if( v.size() != this->vertex.size() ) { throw ArgError(); } 00403 this->vertex = v; 00404 } 00405 00406 bool Map::inside_box( const CartVect& xi, double& tol ) const 00407 { 00408 // bail out early, before doing an expensive NR iteration 00409 // compute box 00410 BoundBox box( this->vertex ); 00411 return box.contains_point( xi.array(), tol ); 00412 } 00413 00414 // 00415 CartVect Map::ievaluate( const CartVect& x, double tol, const CartVect& x0 ) const 00416 { 00417 // TODO: should differentiate between epsilons used for 00418 // Newton Raphson iteration, and epsilons used for curved boundary geometry errors 00419 // right now, fix the tolerance used for NR 00420 tol = 1.0e-10; 00421 const double error_tol_sqr = tol * tol; 00422 double det; 00423 CartVect xi = x0; 00424 CartVect delta = evaluate( xi ) - x; 00425 Matrix3 J; 00426 00427 int iters = 0; 00428 while( delta % delta > error_tol_sqr ) 00429 { 00430 if( ++iters > 10 ) throw Map::EvaluationError( x, vertex ); 00431 00432 J = jacobian( xi ); 00433 det = J.determinant(); 00434 if( det < std::numeric_limits< double >::epsilon() ) throw Map::EvaluationError( x, vertex ); 00435 xi -= J.inverse() * delta; 00436 delta = evaluate( xi ) - x; 00437 } 00438 return xi; 00439 } // Map::ievaluate() 00440 00441 SphericalQuad::SphericalQuad( const std::vector< CartVect >& vertices ) : LinearQuad( vertices ) 00442 { 00443 // project the vertices to the plane tangent at first vertex 00444 v1 = vertex[0]; // member data 00445 double v1v1 = v1 % v1; // this is 1, in general, for unit sphere meshes 00446 for( int j = 1; j < 4; j++ ) 00447 { 00448 // first, bring all vertices in the gnomonic plane 00449 // the new vertex will intersect the plane at vnew 00450 // so that (vnew-v1)%v1 is 0 ( vnew is in the tangent plane, i.e. normal to v1 ) 00451 // pos is the old position of the vertex, and it is in general on the sphere 00452 // vnew = alfa*pos; (alfa*pos-v1)%v1 = 0 <=> alfa*(pos%v1)=v1%v1 <=> alfa = 00453 // v1v1/(pos%v1) 00454 // <=> vnew = ( v1v1/(pos%v1) )*pos 00455 CartVect vnew = v1v1 / ( vertex[j] % v1 ) * vertex[j]; 00456 vertex[j] = vnew; 00457 } 00458 // will compute a transf matrix, such that a new point will be transformed with 00459 // newpos = transf * (vnew-v1), and it will be a point in the 2d plane 00460 // the transformation matrix will be oriented in such a way that orientation will be 00461 // positive 00462 CartVect vx = vertex[1] - v1; // this will become Ox axis 00463 // z axis will be along v1, in such a way that orientation of the quad is positive 00464 // look at the first 2 edges 00465 CartVect vz = vx * ( vertex[2] - vertex[1] ); 00466 vz = vz / vz.length(); 00467 00468 vx = vx / vx.length(); 00469 00470 CartVect vy = vz * vx; 00471 transf = Matrix3( vx[0], vx[1], vx[2], vy[0], vy[1], vy[2], vz[0], vz[1], vz[2] ); 00472 vertex[0] = CartVect( 0. ); 00473 for( int j = 1; j < 4; j++ ) 00474 vertex[j] = transf * ( vertex[j] - v1 ); 00475 } 00476 00477 CartVect SphericalQuad::ievaluate( const CartVect& x, double tol, const CartVect& x0 ) const 00478 { 00479 // project to the plane tangent at first vertex (gnomonic projection) 00480 double v1v1 = v1 % v1; 00481 CartVect vnew = v1v1 / ( x % v1 ) * x; // so that (vnew-v1)%v1 is 0 00482 vnew = transf * ( vnew - v1 ); 00483 // det will be positive now 00484 return Map::ievaluate( vnew, tol, x0 ); 00485 } 00486 00487 bool SphericalQuad::inside_box( const CartVect& pos, double& tol ) const 00488 { 00489 // project to the plane tangent at first vertex 00490 // CartVect v1=vertex[0]; 00491 double v1v1 = v1 % v1; 00492 CartVect vnew = v1v1 / ( pos % v1 ) * pos; // so that (x-v1)%v1 is 0 00493 vnew = transf * ( vnew - v1 ); 00494 return Map::inside_box( vnew, tol ); 00495 } 00496 00497 const double LinearTri::corner[3][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 } }; 00498 00499 LinearTri::LinearTri() : Map( 0 ), det_T( 0.0 ), det_T_inverse( 0.0 ) {} // LinearTri::LinearTri() 00500 00501 LinearTri::~LinearTri() {} 00502 00503 void LinearTri::set_vertices( const std::vector< CartVect >& v ) 00504 { 00505 this->Map::set_vertices( v ); 00506 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, 00507 v[1][2] - v[0][2], v[2][2] - v[0][2], 1 ); 00508 this->T_inverse = this->T.inverse(); 00509 this->det_T = this->T.determinant(); 00510 this->det_T_inverse = ( this->det_T < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / this->det_T ); 00511 } // LinearTri::set_vertices() 00512 00513 bool LinearTri::inside_nat_space( const CartVect& xi, double& tol ) const 00514 { 00515 // linear tri space is a triangle with vertices (0,0,0), (1,0,0), (0,1,0) 00516 // first check if outside bigger box, then below the line x+y=1 00517 return ( xi[0] >= -tol ) && ( xi[1] >= -tol ) && ( xi[2] >= -tol ) && ( xi[2] <= tol ) && 00518 ( xi[0] + xi[1] < 1.0 + tol ); 00519 } 00520 CartVect LinearTri::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const 00521 { 00522 return this->T_inverse * ( x - this->vertex[0] ); 00523 } // LinearTri::ievaluate 00524 00525 double LinearTri::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const 00526 { 00527 double f0 = field_vertex_value[0]; 00528 double f = f0; 00529 for( unsigned i = 1; i < 3; ++i ) 00530 { 00531 f += ( field_vertex_value[i] - f0 ) * xi[i - 1]; 00532 } 00533 return f; 00534 } // LinearTri::evaluate_scalar_field() 00535 00536 double LinearTri::integrate_scalar_field( const double* field_vertex_values ) const 00537 { 00538 double I( 0.0 ); 00539 for( unsigned int i = 0; i < 3; ++i ) 00540 { 00541 I += field_vertex_values[i]; 00542 } 00543 I *= this->det_T / 6.0; // TODO 00544 return I; 00545 } // LinearTri::integrate_scalar_field() 00546 00547 SphericalTri::SphericalTri( const std::vector< CartVect >& vertices ) 00548 { 00549 vertex.resize( vertices.size() ); 00550 vertex = vertices; 00551 // project the vertices to the plane tangent at first vertex 00552 v1 = vertex[0]; // member data 00553 double v1v1 = v1 % v1; // this is 1, in general, for unit sphere meshes 00554 for( int j = 1; j < 3; j++ ) 00555 { 00556 // first, bring all vertices in the gnomonic plane 00557 // the new vertex will intersect the plane at vnew 00558 // so that (vnew-v1)%v1 is 0 ( vnew is in the tangent plane, i.e. normal to v1 ) 00559 // pos is the old position of the vertex, and it is in general on the sphere 00560 // vnew = alfa*pos; (alfa*pos-v1)%v1 = 0 <=> alfa*(pos%v1)=v1%v1 <=> alfa = 00561 // v1v1/(pos%v1) 00562 // <=> vnew = ( v1v1/(pos%v1) )*pos 00563 CartVect vnew = v1v1 / ( vertex[j] % v1 ) * vertex[j]; 00564 vertex[j] = vnew; 00565 } 00566 // will compute a transf matrix, such that a new point will be transformed with 00567 // newpos = transf * (vnew-v1), and it will be a point in the 2d plane 00568 // the transformation matrix will be oriented in such a way that orientation will be 00569 // positive 00570 CartVect vx = vertex[1] - v1; // this will become Ox axis 00571 // z axis will be along v1, in such a way that orientation of the quad is positive 00572 // look at the first 2 edges 00573 CartVect vz = vx * ( vertex[2] - vertex[1] ); 00574 vz = vz / vz.length(); 00575 00576 vx = vx / vx.length(); 00577 00578 CartVect vy = vz * vx; 00579 transf = Matrix3( vx[0], vx[1], vx[2], vy[0], vy[1], vy[2], vz[0], vz[1], vz[2] ); 00580 vertex[0] = CartVect( 0. ); 00581 for( int j = 1; j < 3; j++ ) 00582 vertex[j] = transf * ( vertex[j] - v1 ); 00583 00584 LinearTri::set_vertices( vertex ); 00585 } 00586 00587 CartVect SphericalTri::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const 00588 { 00589 // project to the plane tangent at first vertex (gnomonic projection) 00590 double v1v1 = v1 % v1; 00591 CartVect vnew = v1v1 / ( x % v1 ) * x; // so that (vnew-v1)%v1 is 0 00592 vnew = transf * ( vnew - v1 ); 00593 // det will be positive now 00594 return LinearTri::ievaluate( vnew ); 00595 } 00596 00597 bool SphericalTri::inside_box( const CartVect& pos, double& tol ) const 00598 { 00599 // project to the plane tangent at first vertex 00600 // CartVect v1=vertex[0]; 00601 double v1v1 = v1 % v1; 00602 CartVect vnew = v1v1 / ( pos % v1 ) * pos; // so that (x-v1)%v1 is 0 00603 vnew = transf * ( vnew - v1 ); 00604 return Map::inside_box( vnew, tol ); 00605 } 00606 00607 // filescope for static member data that is cached 00608 const double LinearEdge::corner[2][3] = { { -1, 0, 0 }, { 1, 0, 0 } }; 00609 00610 LinearEdge::LinearEdge() : Map( 0 ) {} // LinearEdge::LinearEdge() 00611 00612 /* For each point, its weight and location are stored as an array. 00613 Hence, the inner dimension is 2, the outer dimension is gauss_count. 00614 We use a one-point Gaussian quadrature, since it integrates linear functions exactly. 00615 */ 00616 const double LinearEdge::gauss[1][2] = { { 2.0, 0.0 } }; 00617 00618 CartVect LinearEdge::evaluate( const CartVect& xi ) const 00619 { 00620 CartVect x( 0.0 ); 00621 for( unsigned i = 0; i < LinearEdge::corner_count; ++i ) 00622 { 00623 const double N_i = ( 1.0 + xi[0] * corner[i][0] ); 00624 x += N_i * this->vertex[i]; 00625 } 00626 x /= LinearEdge::corner_count; 00627 return x; 00628 } // LinearEdge::evaluate 00629 00630 Matrix3 LinearEdge::jacobian( const CartVect& xi ) const 00631 { 00632 Matrix3 J( 0.0 ); 00633 for( unsigned i = 0; i < LinearEdge::corner_count; ++i ) 00634 { 00635 const double xi_p = 1.0 + xi[0] * corner[i][0]; 00636 const double dNi_dxi = corner[i][0] * xi_p; 00637 J( 0, 0 ) += dNi_dxi * vertex[i][0]; 00638 } 00639 J( 1, 1 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */ 00640 J( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */ 00641 J /= LinearEdge::corner_count; 00642 return J; 00643 } // LinearEdge::jacobian() 00644 00645 double LinearEdge::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const 00646 { 00647 double f( 0.0 ); 00648 for( unsigned i = 0; i < LinearEdge::corner_count; ++i ) 00649 { 00650 const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1.0 + xi[1] * corner[i][1] ); 00651 f += N_i * field_vertex_value[i]; 00652 } 00653 f /= LinearEdge::corner_count; 00654 return f; 00655 } // LinearEdge::evaluate_scalar_field() 00656 00657 double LinearEdge::integrate_scalar_field( const double* field_vertex_values ) const 00658 { 00659 double I( 0.0 ); 00660 for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 ) 00661 { 00662 double x1 = this->gauss[j1][1]; 00663 double w1 = this->gauss[j1][0]; 00664 CartVect x( x1, 0.0, 0.0 ); 00665 I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * this->det_jacobian( x ); 00666 } 00667 return I; 00668 } // LinearEdge::integrate_scalar_field() 00669 00670 bool LinearEdge::inside_nat_space( const CartVect& xi, double& tol ) const 00671 { 00672 // just look at the box+tol here 00673 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ); 00674 } 00675 00676 const double LinearHex::corner[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 }, 00677 { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } }; 00678 00679 LinearHex::LinearHex() : Map( 0 ) {} // LinearHex::LinearHex() 00680 00681 LinearHex::~LinearHex() {} 00682 /* For each point, its weight and location are stored as an array. 00683 Hence, the inner dimension is 2, the outer dimension is gauss_count. 00684 We use a one-point Gaussian quadrature, since it integrates linear functions exactly. 00685 */ 00686 // const double LinearHex::gauss[1][2] = { { 2.0, 0.0 } }; 00687 const double LinearHex::gauss[2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } }; 00688 // const double LinearHex::gauss[4][2] = { { 0.3478548451, -0.8611363116 }, 00689 // { 0.6521451549, -0.3399810436 }, 00690 // { 0.6521451549, 0.3399810436 }, 00691 // { 0.3478548451, 0.8611363116 } }; 00692 00693 CartVect LinearHex::evaluate( const CartVect& xi ) const 00694 { 00695 CartVect x( 0.0 ); 00696 for( unsigned i = 0; i < 8; ++i ) 00697 { 00698 const double N_i = 00699 ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] ) * ( 1 + xi[2] * corner[i][2] ); 00700 x += N_i * this->vertex[i]; 00701 } 00702 x *= 0.125; 00703 return x; 00704 } // LinearHex::evaluate 00705 00706 Matrix3 LinearHex::jacobian( const CartVect& xi ) const 00707 { 00708 Matrix3 J( 0.0 ); 00709 for( unsigned i = 0; i < 8; ++i ) 00710 { 00711 const double xi_p = 1 + xi[0] * corner[i][0]; 00712 const double eta_p = 1 + xi[1] * corner[i][1]; 00713 const double zeta_p = 1 + xi[2] * corner[i][2]; 00714 const double dNi_dxi = corner[i][0] * eta_p * zeta_p; 00715 const double dNi_deta = corner[i][1] * xi_p * zeta_p; 00716 const double dNi_dzeta = corner[i][2] * xi_p * eta_p; 00717 J( 0, 0 ) += dNi_dxi * vertex[i][0]; 00718 J( 1, 0 ) += dNi_dxi * vertex[i][1]; 00719 J( 2, 0 ) += dNi_dxi * vertex[i][2]; 00720 J( 0, 1 ) += dNi_deta * vertex[i][0]; 00721 J( 1, 1 ) += dNi_deta * vertex[i][1]; 00722 J( 2, 1 ) += dNi_deta * vertex[i][2]; 00723 J( 0, 2 ) += dNi_dzeta * vertex[i][0]; 00724 J( 1, 2 ) += dNi_dzeta * vertex[i][1]; 00725 J( 2, 2 ) += dNi_dzeta * vertex[i][2]; 00726 } 00727 return J *= 0.125; 00728 } // LinearHex::jacobian() 00729 00730 double LinearHex::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const 00731 { 00732 double f( 0.0 ); 00733 for( unsigned i = 0; i < 8; ++i ) 00734 { 00735 const double N_i = 00736 ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] ) * ( 1 + xi[2] * corner[i][2] ); 00737 f += N_i * field_vertex_value[i]; 00738 } 00739 f *= 0.125; 00740 return f; 00741 } // LinearHex::evaluate_scalar_field() 00742 00743 double LinearHex::integrate_scalar_field( const double* field_vertex_values ) const 00744 { 00745 double I( 0.0 ); 00746 for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 ) 00747 { 00748 double x1 = this->gauss[j1][1]; 00749 double w1 = this->gauss[j1][0]; 00750 for( unsigned int j2 = 0; j2 < this->gauss_count; ++j2 ) 00751 { 00752 double x2 = this->gauss[j2][1]; 00753 double w2 = this->gauss[j2][0]; 00754 for( unsigned int j3 = 0; j3 < this->gauss_count; ++j3 ) 00755 { 00756 double x3 = this->gauss[j3][1]; 00757 double w3 = this->gauss[j3][0]; 00758 CartVect x( x1, x2, x3 ); 00759 I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * w2 * w3 * this->det_jacobian( x ); 00760 } 00761 } 00762 } 00763 return I; 00764 } // LinearHex::integrate_scalar_field() 00765 00766 bool LinearHex::inside_nat_space( const CartVect& xi, double& tol ) const 00767 { 00768 // just look at the box+tol here 00769 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) && 00770 ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol ); 00771 } 00772 00773 // those are not just the corners, but for simplicity, keep this name 00774 // 00775 const int QuadraticHex::corner[27][3] = { 00776 { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, // corner nodes: 0-7 00777 { -1, 1, -1 }, // mid-edge nodes: 8-19 00778 { -1, -1, 1 }, // center-face nodes 20-25 center node 26 00779 { 1, -1, 1 }, // 00780 { 1, 1, 1 }, { -1, 1, 1 }, // 4 ----- 19 ----- 7 00781 { 0, -1, -1 }, // . | . | 00782 { 1, 0, -1 }, // 16 25 18 | 00783 { 0, 1, -1 }, // . | . | 00784 { -1, 0, -1 }, // 5 ----- 17 ----- 6 | 00785 { -1, -1, 0 }, // | 12 | 23 15 00786 { 1, -1, 0 }, // | | | 00787 { 1, 1, 0 }, // | 20 | 26 | 22 | 00788 { -1, 1, 0 }, // | | | 00789 { 0, -1, 1 }, // 13 21 | 14 | 00790 { 1, 0, 1 }, // | 0 ----- 11 ----- 3 00791 { 0, 1, 1 }, // | . | . 00792 { -1, 0, 1 }, // | 8 24 | 10 00793 { 0, -1, 0 }, // | . | . 00794 { 1, 0, 0 }, // 1 ----- 9 ----- 2 00795 { 0, 1, 0 }, // 00796 { -1, 0, 0 }, { 0, 0, -1 }, { 0, 0, 1 }, { 0, 0, 0 } 00797 }; 00798 // QuadraticHex::QuadraticHex(const std::vector<CartVect>& vertices) : Map(vertices){}; 00799 QuadraticHex::QuadraticHex() : Map( 0 ) {} 00800 00801 QuadraticHex::~QuadraticHex() {} 00802 double SH( const int i, const double xi ) 00803 { 00804 switch( i ) 00805 { 00806 case -1: 00807 return ( xi * xi - xi ) / 2; 00808 case 0: 00809 return 1 - xi * xi; 00810 case 1: 00811 return ( xi * xi + xi ) / 2; 00812 default: 00813 return 0.; 00814 } 00815 } 00816 double DSH( const int i, const double xi ) 00817 { 00818 switch( i ) 00819 { 00820 case -1: 00821 return xi - 0.5; 00822 case 0: 00823 return -2 * xi; 00824 case 1: 00825 return xi + 0.5; 00826 default: 00827 return 0.; 00828 } 00829 } 00830 00831 CartVect QuadraticHex::evaluate( const CartVect& xi ) const 00832 { 00833 00834 CartVect x( 0.0 ); 00835 for( int i = 0; i < 27; i++ ) 00836 { 00837 const double sh = SH( corner[i][0], xi[0] ) * SH( corner[i][1], xi[1] ) * SH( corner[i][2], xi[2] ); 00838 x += sh * vertex[i]; 00839 } 00840 00841 return x; 00842 } 00843 00844 bool QuadraticHex::inside_nat_space( const CartVect& xi, double& tol ) const 00845 { // just look at the box+tol here 00846 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) && 00847 ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol ); 00848 } 00849 00850 Matrix3 QuadraticHex::jacobian( const CartVect& xi ) const 00851 { 00852 Matrix3 J( 0.0 ); 00853 for( int i = 0; i < 27; i++ ) 00854 { 00855 const double sh[3] = { SH( corner[i][0], xi[0] ), SH( corner[i][1], xi[1] ), SH( corner[i][2], xi[2] ) }; 00856 const double dsh[3] = { DSH( corner[i][0], xi[0] ), DSH( corner[i][1], xi[1] ), 00857 DSH( corner[i][2], xi[2] ) }; 00858 00859 for( int j = 0; j < 3; j++ ) 00860 { 00861 J( j, 0 ) += dsh[0] * sh[1] * sh[2] * vertex[i][j]; // dxj/dr first column 00862 J( j, 1 ) += sh[0] * dsh[1] * sh[2] * vertex[i][j]; // dxj/ds 00863 J( j, 2 ) += sh[0] * sh[1] * dsh[2] * vertex[i][j]; // dxj/dt 00864 } 00865 } 00866 00867 return J; 00868 } 00869 double QuadraticHex::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const 00870 { 00871 double x = 0.0; 00872 for( int i = 0; i < 27; i++ ) 00873 { 00874 const double sh = SH( corner[i][0], xi[0] ) * SH( corner[i][1], xi[1] ) * SH( corner[i][2], xi[2] ); 00875 x += sh * field_vertex_values[i]; 00876 } 00877 00878 return x; 00879 } 00880 double QuadraticHex::integrate_scalar_field( const double* /*field_vertex_values*/ ) const 00881 { 00882 return 0.; // TODO: gaussian integration , probably 2x2x2 00883 } 00884 00885 const double LinearTet::corner[4][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } }; 00886 00887 LinearTet::LinearTet() : Map( 0 ), det_T( 0.0 ), det_T_inverse( 0.0 ) {} // LinearTet::LinearTet() 00888 00889 LinearTet::~LinearTet() {} 00890 00891 void LinearTet::set_vertices( const std::vector< CartVect >& v ) 00892 { 00893 this->Map::set_vertices( v ); 00894 this->T = 00895 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], 00896 v[3][1] - v[0][1], v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] ); 00897 this->T_inverse = this->T.inverse(); 00898 this->det_T = this->T.determinant(); 00899 this->det_T_inverse = ( this->det_T < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / this->det_T ); 00900 } // LinearTet::set_vertices() 00901 00902 double LinearTet::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const 00903 { 00904 double f0 = field_vertex_value[0]; 00905 double f = f0; 00906 for( unsigned i = 1; i < 4; ++i ) 00907 { 00908 f += ( field_vertex_value[i] - f0 ) * xi[i - 1]; 00909 } 00910 return f; 00911 } // LinearTet::evaluate_scalar_field() 00912 00913 CartVect LinearTet::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const 00914 { 00915 return this->T_inverse * ( x - this->vertex[0] ); 00916 } // LinearTet::ievaluate 00917 00918 double LinearTet::integrate_scalar_field( const double* field_vertex_values ) const 00919 { 00920 double I( 0.0 ); 00921 for( unsigned int i = 0; i < 4; ++i ) 00922 { 00923 I += field_vertex_values[i]; 00924 } 00925 I *= this->det_T / 24.0; 00926 return I; 00927 } // LinearTet::integrate_scalar_field() 00928 bool LinearTet::inside_nat_space( const CartVect& xi, double& tol ) const 00929 { 00930 // linear tet space is a tetra with vertices (0,0,0), (1,0,0), (0,1,0), (0, 0, 1) 00931 // first check if outside bigger box, then below the plane x+y+z=1 00932 return ( xi[0] >= -tol ) && ( xi[1] >= -tol ) && ( xi[2] >= -tol ) && ( xi[0] + xi[1] + xi[2] < 1.0 + tol ); 00933 } 00934 // SpectralHex 00935 00936 // filescope for static member data that is cached 00937 int SpectralHex::_n; 00938 realType* SpectralHex::_z[3]; 00939 lagrange_data SpectralHex::_ld[3]; 00940 opt_data_3 SpectralHex::_data; 00941 realType* SpectralHex::_odwork; 00942 00943 bool SpectralHex::_init = false; 00944 00945 SpectralHex::SpectralHex() : Map( 0 ) 00946 { 00947 _xyz[0] = _xyz[1] = _xyz[2] = NULL; 00948 } 00949 // the preferred constructor takes pointers to GL blocked positions 00950 SpectralHex::SpectralHex( int order, double* x, double* y, double* z ) : Map( 0 ) 00951 { 00952 Init( order ); 00953 _xyz[0] = x; 00954 _xyz[1] = y; 00955 _xyz[2] = z; 00956 } 00957 SpectralHex::SpectralHex( int order ) : Map( 0 ) 00958 { 00959 Init( order ); 00960 _xyz[0] = _xyz[1] = _xyz[2] = NULL; 00961 } 00962 SpectralHex::~SpectralHex() 00963 { 00964 if( _init ) freedata(); 00965 _init = false; 00966 } 00967 void SpectralHex::Init( int order ) 00968 { 00969 if( _init && _n == order ) return; 00970 if( _init && _n != order ) 00971 { 00972 // TODO: free data cached 00973 freedata(); 00974 } 00975 // compute stuff that depends only on order 00976 _init = true; 00977 _n = order; 00978 // triplicates! n is the same in all directions !!! 00979 for( int d = 0; d < 3; d++ ) 00980 { 00981 _z[d] = tmalloc( realType, _n ); 00982 lobatto_nodes( _z[d], _n ); 00983 lagrange_setup( &_ld[d], _z[d], _n ); 00984 } 00985 opt_alloc_3( &_data, _ld ); 00986 00987 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n; 00988 _odwork = tmalloc( realType, 6 * nf + 9 * ne + nw ); 00989 } 00990 void SpectralHex::freedata() 00991 { 00992 for( int d = 0; d < 3; d++ ) 00993 { 00994 free( _z[d] ); 00995 lagrange_free( &_ld[d] ); 00996 } 00997 opt_free_3( &_data ); 00998 free( _odwork ); 00999 } 01000 01001 void SpectralHex::set_gl_points( double* x, double* y, double* z ) 01002 { 01003 _xyz[0] = x; 01004 _xyz[1] = y; 01005 _xyz[2] = z; 01006 } 01007 CartVect SpectralHex::evaluate( const CartVect& xi ) const 01008 { 01009 // piece that we shouldn't want to cache 01010 int d = 0; 01011 for( d = 0; d < 3; d++ ) 01012 { 01013 lagrange_0( &_ld[d], xi[d] ); 01014 } 01015 CartVect result; 01016 for( d = 0; d < 3; d++ ) 01017 { 01018 result[d] = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, 01019 _xyz[d], // this is the "field" 01020 _odwork ); 01021 } 01022 return result; 01023 } 01024 // replicate the functionality of hex_findpt 01025 CartVect SpectralHex::ievaluate( CartVect const& xyz, double tol, const CartVect& x0 ) const 01026 { 01027 // find nearest point 01028 realType x_star[3]; 01029 xyz.get( x_star ); 01030 01031 realType r[3] = { 0, 0, 0 }; // initial guess for parametric coords 01032 x0.get( r ); 01033 unsigned c = opt_no_constraints_3; 01034 realType dist = opt_findpt_3( &_data, (const realType**)_xyz, x_star, r, &c ); 01035 // if it did not converge, get out with throw... 01036 if( dist > 10 * tol ) // outside the element 01037 { 01038 std::vector< CartVect > dummy; 01039 throw Map::EvaluationError( xyz, dummy ); 01040 } 01041 // c tells us if we landed inside the element or exactly on a face, edge, or node 01042 // also, dist shows the distance to the computed point. 01043 // copy parametric coords back 01044 return CartVect( r ); 01045 } 01046 Matrix3 SpectralHex::jacobian( const CartVect& xi ) const 01047 { 01048 realType x_i[3]; 01049 xi.get( x_i ); 01050 // set the positions of GL nodes, before evaluations 01051 _data.elx[0] = _xyz[0]; 01052 _data.elx[1] = _xyz[1]; 01053 _data.elx[2] = _xyz[2]; 01054 opt_vol_set_intp_3( &_data, x_i ); 01055 Matrix3 J( 0. ); 01056 // it is organized differently 01057 J( 0, 0 ) = _data.jac[0]; // dx/dr 01058 J( 0, 1 ) = _data.jac[1]; // dx/ds 01059 J( 0, 2 ) = _data.jac[2]; // dx/dt 01060 J( 1, 0 ) = _data.jac[3]; // dy/dr 01061 J( 1, 1 ) = _data.jac[4]; // dy/ds 01062 J( 1, 2 ) = _data.jac[5]; // dy/dt 01063 J( 2, 0 ) = _data.jac[6]; // dz/dr 01064 J( 2, 1 ) = _data.jac[7]; // dz/ds 01065 J( 2, 2 ) = _data.jac[8]; // dz/dt 01066 return J; 01067 } 01068 double SpectralHex::evaluate_scalar_field( const CartVect& xi, const double* field ) const 01069 { 01070 // piece that we shouldn't want to cache 01071 int d; 01072 for( d = 0; d < 3; d++ ) 01073 { 01074 lagrange_0( &_ld[d], xi[d] ); 01075 } 01076 01077 double value = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork ); 01078 return value; 01079 } 01080 double SpectralHex::integrate_scalar_field( const double* field_vertex_values ) const 01081 { 01082 // set the position of GL points 01083 // set the positions of GL nodes, before evaluations 01084 _data.elx[0] = _xyz[0]; 01085 _data.elx[1] = _xyz[1]; 01086 _data.elx[2] = _xyz[2]; 01087 double xi[3]; 01088 // triple loop; the most inner loop is in r direction, then s, then t 01089 double integral = 0.; 01090 // double volume = 0; 01091 int index = 0; // used fr the inner loop 01092 for( int k = 0; k < _n; k++ ) 01093 { 01094 xi[2] = _ld[2].z[k]; 01095 // double wk= _ld[2].w[k]; 01096 for( int j = 0; j < _n; j++ ) 01097 { 01098 xi[1] = _ld[1].z[j]; 01099 // double wj= _ld[1].w[j]; 01100 for( int i = 0; i < _n; i++ ) 01101 { 01102 xi[0] = _ld[0].z[i]; 01103 // double wi= _ld[0].w[i]; 01104 opt_vol_set_intp_3( &_data, xi ); 01105 double wk = _ld[2].J[k]; 01106 double wj = _ld[1].J[j]; 01107 double wi = _ld[0].J[i]; 01108 Matrix3 J( 0. ); 01109 // it is organized differently 01110 J( 0, 0 ) = _data.jac[0]; // dx/dr 01111 J( 0, 1 ) = _data.jac[1]; // dx/ds 01112 J( 0, 2 ) = _data.jac[2]; // dx/dt 01113 J( 1, 0 ) = _data.jac[3]; // dy/dr 01114 J( 1, 1 ) = _data.jac[4]; // dy/ds 01115 J( 1, 2 ) = _data.jac[5]; // dy/dt 01116 J( 2, 0 ) = _data.jac[6]; // dz/dr 01117 J( 2, 1 ) = _data.jac[7]; // dz/ds 01118 J( 2, 2 ) = _data.jac[8]; // dz/dt 01119 double bm = wk * wj * wi * J.determinant(); 01120 integral += bm * field_vertex_values[index++]; 01121 // volume +=bm; 01122 } 01123 } 01124 } 01125 // std::cout << "volume: " << volume << "\n"; 01126 return integral; 01127 } 01128 // this is the same as a linear hex, although we should not need it 01129 bool SpectralHex::inside_nat_space( const CartVect& xi, double& tol ) const 01130 { 01131 // just look at the box+tol here 01132 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) && 01133 ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol ); 01134 } 01135 01136 // SpectralHex 01137 01138 // filescope for static member data that is cached 01139 const double LinearQuad::corner[4][3] = { { -1, -1, 0 }, { 1, -1, 0 }, { 1, 1, 0 }, { -1, 1, 0 } }; 01140 01141 LinearQuad::LinearQuad() : Map( 0 ) {} // LinearQuad::LinearQuad() 01142 01143 LinearQuad::~LinearQuad() {} 01144 01145 /* For each point, its weight and location are stored as an array. 01146 Hence, the inner dimension is 2, the outer dimension is gauss_count. 01147 We use a one-point Gaussian quadrature, since it integrates linear functions exactly. 01148 */ 01149 const double LinearQuad::gauss[1][2] = { { 2.0, 0.0 } }; 01150 01151 CartVect LinearQuad::evaluate( const CartVect& xi ) const 01152 { 01153 CartVect x( 0.0 ); 01154 for( unsigned i = 0; i < LinearQuad::corner_count; ++i ) 01155 { 01156 const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] ); 01157 x += N_i * this->vertex[i]; 01158 } 01159 x /= LinearQuad::corner_count; 01160 return x; 01161 } // LinearQuad::evaluate 01162 01163 Matrix3 LinearQuad::jacobian( const CartVect& xi ) const 01164 { 01165 // this basically ignores the z component: xi[2] or vertex[][2] 01166 Matrix3 J( 0.0 ); 01167 for( unsigned i = 0; i < LinearQuad::corner_count; ++i ) 01168 { 01169 const double xi_p = 1 + xi[0] * corner[i][0]; 01170 const double eta_p = 1 + xi[1] * corner[i][1]; 01171 const double dNi_dxi = corner[i][0] * eta_p; 01172 const double dNi_deta = corner[i][1] * xi_p; 01173 J( 0, 0 ) += dNi_dxi * vertex[i][0]; 01174 J( 1, 0 ) += dNi_dxi * vertex[i][1]; 01175 J( 0, 1 ) += dNi_deta * vertex[i][0]; 01176 J( 1, 1 ) += dNi_deta * vertex[i][1]; 01177 } 01178 J( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */ 01179 J /= LinearQuad::corner_count; 01180 return J; 01181 } // LinearQuad::jacobian() 01182 01183 double LinearQuad::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const 01184 { 01185 double f( 0.0 ); 01186 for( unsigned i = 0; i < LinearQuad::corner_count; ++i ) 01187 { 01188 const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] ); 01189 f += N_i * field_vertex_value[i]; 01190 } 01191 f /= LinearQuad::corner_count; 01192 return f; 01193 } // LinearQuad::evaluate_scalar_field() 01194 01195 double LinearQuad::integrate_scalar_field( const double* field_vertex_values ) const 01196 { 01197 double I( 0.0 ); 01198 for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 ) 01199 { 01200 double x1 = this->gauss[j1][1]; 01201 double w1 = this->gauss[j1][0]; 01202 for( unsigned int j2 = 0; j2 < this->gauss_count; ++j2 ) 01203 { 01204 double x2 = this->gauss[j2][1]; 01205 double w2 = this->gauss[j2][0]; 01206 CartVect x( x1, x2, 0.0 ); 01207 I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * w2 * this->det_jacobian( x ); 01208 } 01209 } 01210 return I; 01211 } // LinearQuad::integrate_scalar_field() 01212 01213 bool LinearQuad::inside_nat_space( const CartVect& xi, double& tol ) const 01214 { 01215 // just look at the box+tol here 01216 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ); 01217 } 01218 01219 // filescope for static member data that is cached 01220 int SpectralQuad::_n; 01221 realType* SpectralQuad::_z[2]; 01222 lagrange_data SpectralQuad::_ld[2]; 01223 opt_data_2 SpectralQuad::_data; 01224 realType* SpectralQuad::_odwork; 01225 realType* SpectralQuad::_glpoints; 01226 bool SpectralQuad::_init = false; 01227 01228 SpectralQuad::SpectralQuad() : Map( 0 ) 01229 { 01230 _xyz[0] = _xyz[1] = _xyz[2] = NULL; 01231 } 01232 // the preferred constructor takes pointers to GL blocked positions 01233 SpectralQuad::SpectralQuad( int order, double* x, double* y, double* z ) : Map( 0 ) 01234 { 01235 Init( order ); 01236 _xyz[0] = x; 01237 _xyz[1] = y; 01238 _xyz[2] = z; 01239 } 01240 SpectralQuad::SpectralQuad( int order ) : Map( 4 ) 01241 { 01242 Init( order ); 01243 _xyz[0] = _xyz[1] = _xyz[2] = NULL; 01244 } 01245 SpectralQuad::~SpectralQuad() 01246 { 01247 if( _init ) freedata(); 01248 _init = false; 01249 } 01250 void SpectralQuad::Init( int order ) 01251 { 01252 if( _init && _n == order ) return; 01253 if( _init && _n != order ) 01254 { 01255 // TODO: free data cached 01256 freedata(); 01257 } 01258 // compute stuff that depends only on order 01259 _init = true; 01260 _n = order; 01261 // duplicates! n is the same in all directions !!! 01262 for( int d = 0; d < 2; d++ ) 01263 { 01264 _z[d] = tmalloc( realType, _n ); 01265 lobatto_nodes( _z[d], _n ); 01266 lagrange_setup( &_ld[d], _z[d], _n ); 01267 } 01268 opt_alloc_2( &_data, _ld ); 01269 01270 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n; 01271 _odwork = tmalloc( realType, 6 * nf + 9 * ne + nw ); 01272 _glpoints = tmalloc( realType, 3 * nf ); 01273 } 01274 01275 void SpectralQuad::freedata() 01276 { 01277 for( int d = 0; d < 2; d++ ) 01278 { 01279 free( _z[d] ); 01280 lagrange_free( &_ld[d] ); 01281 } 01282 opt_free_2( &_data ); 01283 free( _odwork ); 01284 free( _glpoints ); 01285 } 01286 01287 void SpectralQuad::set_gl_points( double* x, double* y, double* z ) 01288 { 01289 _xyz[0] = x; 01290 _xyz[1] = y; 01291 _xyz[2] = z; 01292 } 01293 CartVect SpectralQuad::evaluate( const CartVect& xi ) const 01294 { 01295 // piece that we shouldn't want to cache 01296 int d = 0; 01297 for( d = 0; d < 2; d++ ) 01298 { 01299 lagrange_0( &_ld[d], xi[d] ); 01300 } 01301 CartVect result; 01302 for( d = 0; d < 3; d++ ) 01303 { 01304 result[d] = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _xyz[d], _odwork ); 01305 } 01306 return result; 01307 } 01308 // replicate the functionality of hex_findpt 01309 CartVect SpectralQuad::ievaluate( CartVect const& xyz, double /*tol*/, const CartVect& /*x0*/ ) const 01310 { 01311 // find nearest point 01312 realType x_star[3]; 01313 xyz.get( x_star ); 01314 01315 realType r[2] = { 0, 0 }; // initial guess for parametric coords 01316 unsigned c = opt_no_constraints_3; 01317 realType dist = opt_findpt_2( &_data, (const realType**)_xyz, x_star, r, &c ); 01318 // if it did not converge, get out with throw... 01319 if( dist > 0.9e+30 ) 01320 { 01321 std::vector< CartVect > dummy; 01322 throw Map::EvaluationError( xyz, dummy ); 01323 } 01324 01325 // c tells us if we landed inside the element or exactly on a face, edge, or node 01326 // also, dist shows the distance to the computed point. 01327 // copy parametric coords back 01328 return CartVect( r[0], r[1], 0. ); 01329 } 01330 01331 Matrix3 SpectralQuad::jacobian( const CartVect& /*xi*/ ) const 01332 { 01333 // not implemented 01334 Matrix3 J( 0. ); 01335 return J; 01336 } 01337 01338 double SpectralQuad::evaluate_scalar_field( const CartVect& xi, const double* field ) const 01339 { 01340 // piece that we shouldn't want to cache 01341 int d; 01342 for( d = 0; d < 2; d++ ) 01343 { 01344 lagrange_0( &_ld[d], xi[d] ); 01345 } 01346 01347 double value = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, field, _odwork ); 01348 return value; 01349 } 01350 double SpectralQuad::integrate_scalar_field( const double* /*field_vertex_values*/ ) const 01351 { 01352 return 0.; // not implemented 01353 } 01354 // this is the same as a linear hex, although we should not need it 01355 bool SpectralQuad::inside_nat_space( const CartVect& xi, double& tol ) const 01356 { 01357 // just look at the box+tol here 01358 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ); 01359 } 01360 // something we don't do for spectral hex; we do it here because 01361 // we do not store the position of gl points in a tag yet 01362 void SpectralQuad::compute_gl_positions() 01363 { 01364 // will need to use shape functions on a simple linear quad to compute gl points 01365 // so we know the position of gl points in parametric space, we will just compute those 01366 // from the 3d vertex position (corner nodes of the quad), using simple mapping 01367 assert( this->vertex.size() == 4 ); 01368 static double corner_xi[4][2] = { { -1., -1. }, { 1., -1. }, { 1., 1. }, { -1., 1. } }; 01369 // we will use the cached lobatto nodes in parametric space _z[d] (the same in both 01370 // directions) 01371 int indexGL = 0; 01372 int n2 = _n * _n; 01373 for( int i = 0; i < _n; i++ ) 01374 { 01375 double eta = _z[0][i]; 01376 for( int j = 0; j < _n; j++ ) 01377 { 01378 double csi = _z[1][j]; // we could really use the same _z[0] array of lobatto nodes 01379 CartVect pos( 0.0 ); 01380 for( int k = 0; k < 4; k++ ) 01381 { 01382 const double N_k = ( 1 + csi * corner_xi[k][0] ) * ( 1 + eta * corner_xi[k][1] ); 01383 pos += N_k * vertex[k]; 01384 } 01385 pos *= 0.25; // these are x, y, z of gl points; reorder them 01386 _glpoints[indexGL] = pos[0]; // x 01387 _glpoints[indexGL + n2] = pos[1]; // y 01388 _glpoints[indexGL + 2 * n2] = pos[2]; // z 01389 indexGL++; 01390 } 01391 } 01392 // now, we can set the _xyz pointers to internal memory allocated to these! 01393 _xyz[0] = &( _glpoints[0] ); 01394 _xyz[1] = &( _glpoints[n2] ); 01395 _xyz[2] = &( _glpoints[2 * n2] ); 01396 } 01397 void SpectralQuad::get_gl_points( double*& x, double*& y, double*& z, int& psize ) 01398 { 01399 x = (double*)_xyz[0]; 01400 y = (double*)_xyz[1]; 01401 z = (double*)_xyz[2]; 01402 psize = _n * _n; 01403 } 01404 } // namespace Element 01405 01406 } // namespace moab