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