moab
ElemUtil.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines