 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 } }; 64 CartVect LinearHexMap::center_xi() const 65 { 66 return CartVect( 0.0 ); 67 } 68 69 CartVect LinearHexMap::evaluate( const CartVect& xi ) const 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 95 Matrix3 LinearHexMap::jacobian( const CartVect& xi ) const 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 169 bool point_in_trilinear_hex( const CartVect* hex, 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 396 Map::~Map() {} 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 507 LinearTri::~LinearTri() {} 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 590 LinearTri::set_vertices( vertex ); 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 624 CartVect LinearEdge::evaluate( const CartVect& xi ) const 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 } 632 x /= LinearEdge::corner_count; 633 return x; 634 } // LinearEdge::evaluate 635 636 Matrix3 LinearEdge::jacobian( const CartVect& xi ) const 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 */ 647 J /= LinearEdge::corner_count; 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 } 659 f /= LinearEdge::corner_count; 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 687 LinearHex::~LinearHex() {} 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 699 CartVect LinearHex::evaluate( const CartVect& xi ) const 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 712 Matrix3 LinearHex::jacobian( const CartVect& xi ) const 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){}; 804 QuadraticHex::QuadraticHex() : Map( 0 ) {} 805 806 QuadraticHex::~QuadraticHex() {} 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 836 CartVect QuadraticHex::evaluate( const CartVect& xi ) const 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 855 Matrix3 QuadraticHex::jacobian( const CartVect& xi ) const 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 894 LinearTet::~LinearTet() {} 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; 943 realType* SpectralHex::_z[3]; 944 lagrange_data SpectralHex::_ld[3]; 945 opt_data_3 SpectralHex::_data; 946 realType* SpectralHex::_odwork; 947 948 bool SpectralHex::_init = false; 949 950 SpectralHex::SpectralHex() : Map( 0 ) 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 } 967 SpectralHex::~SpectralHex() 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 1148 LinearQuad::~LinearQuad() {} 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 1156 CartVect LinearQuad::evaluate( const CartVect& xi ) const 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 } 1164 x /= LinearQuad::corner_count; 1165 return x; 1166 } // LinearQuad::evaluate 1167 1168 Matrix3 LinearQuad::jacobian( const CartVect& xi ) const 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 */ 1184 J /= LinearQuad::corner_count; 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 } 1196 f /= LinearQuad::corner_count; 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; 1226 realType* SpectralQuad::_z[2]; 1227 lagrange_data SpectralQuad::_ld[2]; 1228 opt_data_2 SpectralQuad::_data; 1229 realType* SpectralQuad::_odwork; 1230 realType* SpectralQuad::_glpoints; 1231 bool SpectralQuad::_init = false; 1232 1233 SpectralQuad::SpectralQuad() : Map( 0 ) 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 } 1250 SpectralQuad::~SpectralQuad() 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 } 1298 CartVect SpectralQuad::evaluate( const CartVect& xi ) const 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 1367 void SpectralQuad::compute_gl_positions() 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
1.9.1.
