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