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