ElemUtil.cpp

Go to the documentation of this file.

#include <iostream>
#include <limits>
#include <cassert>
 
#include "ElemUtil.hpp"
#include "moab/BoundBox.hpp"
 
namespace moab
{
namespace ElemUtil
{
 
 bool debug = false;
 
 /**\brief Class representing a 3-D mapping function (e.g. shape function for volume element) */
 class VolMap
 {
 public:
 /**\brief Return \f$\vec \xi\f$ corresponding to logical center of element */
 virtual CartVect center_xi() const = 0;
 /**\brief Evaluate mapping function (calculate \f$\vec x = F($\vec \xi)\f$ )*/
 virtual CartVect evaluate( const CartVect& xi ) const = 0;
 /**\brief Evaluate Jacobian of mapping function */
 virtual Matrix3 jacobian( const CartVect& xi ) const = 0;
 /**\brief Evaluate inverse of mapping function (calculate \f$\vec \xi = F^-1($\vec x)\f$ )*/
 bool solve_inverse( const CartVect& x, CartVect& xi, double tol ) const;
 };
 
 bool VolMap::solve_inverse( const CartVect& x, CartVect& xi, double tol ) const
 {
 const double error_tol_sqr = tol * tol;
 double det;
 xi = center_xi();
 CartVect delta = evaluate( xi ) - x;
 Matrix3 J;
 while( delta % delta > error_tol_sqr )
 {
 J = jacobian( xi );
 det = J.determinant();
 if( det < std::numeric_limits< double >::epsilon() ) return false;
 xi -= J.inverse() * delta;
 delta = evaluate( xi ) - x;
 }
 return true;
 }
 
 /**\brief Shape function for trilinear hexahedron */
 class LinearHexMap : public VolMap
 {
 public:
 LinearHexMap( const CartVect* corner_coords ) : corners( corner_coords ) {}
 virtual CartVect center_xi() const;
 virtual CartVect evaluate( const CartVect& xi ) const;
 virtual double evaluate_scalar_field( const CartVect& xi, const double* f_vals ) const;
 virtual Matrix3 jacobian( const CartVect& xi ) const;
 
 private:
 const CartVect* corners;
 static const double corner_xi[8][3];
 };
 
 const double LinearHexMap::corner_xi[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
 { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } };
 CartVect LinearHexMap::center_xi() const
 {
 return CartVect( 0.0 );
 }
 
 CartVect LinearHexMap::evaluate( const CartVect& xi ) const
 {
 CartVect x( 0.0 );
 for( unsigned i = 0; i < 8; ++i )
 {
 const double N_i =
 ( 1 + xi[0] * corner_xi[i][0] ) * ( 1 + xi[1] * corner_xi[i][1] ) * ( 1 + xi[2] * corner_xi[i][2] );
 x += N_i * corners[i];
 }
 x *= 0.125;
 return x;
 }
 
 double LinearHexMap::evaluate_scalar_field( const CartVect& xi, const double* f_vals ) const
 {
 double f( 0.0 );
 for( unsigned i = 0; i < 8; ++i )
 {
 const double N_i =
 ( 1 + xi[0] * corner_xi[i][0] ) * ( 1 + xi[1] * corner_xi[i][1] ) * ( 1 + xi[2] * corner_xi[i][2] );
 f += N_i * f_vals[i];
 }
 f *= 0.125;
 return f;
 }
 
 Matrix3 LinearHexMap::jacobian( const CartVect& xi ) const
 {
 Matrix3 J( 0.0 );
 for( unsigned i = 0; i < 8; ++i )
 {
 const double xi_p = 1 + xi[0] * corner_xi[i][0];
 const double eta_p = 1 + xi[1] * corner_xi[i][1];
 const double zeta_p = 1 + xi[2] * corner_xi[i][2];
 const double dNi_dxi = corner_xi[i][0] * eta_p * zeta_p;
 const double dNi_deta = corner_xi[i][1] * xi_p * zeta_p;
 const double dNi_dzeta = corner_xi[i][2] * xi_p * eta_p;
 J( 0, 0 ) += dNi_dxi * corners[i][0];
 J( 1, 0 ) += dNi_dxi * corners[i][1];
 J( 2, 0 ) += dNi_dxi * corners[i][2];
 J( 0, 1 ) += dNi_deta * corners[i][0];
 J( 1, 1 ) += dNi_deta * corners[i][1];
 J( 2, 1 ) += dNi_deta * corners[i][2];
 J( 0, 2 ) += dNi_dzeta * corners[i][0];
 J( 1, 2 ) += dNi_dzeta * corners[i][1];
 J( 2, 2 ) += dNi_dzeta * corners[i][2];
 }
 return J *= 0.125;
 }
 
 bool nat_coords_trilinear_hex( const CartVect* corner_coords, const CartVect& x, CartVect& xi, double tol )
 {
 return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
 }
 //
 // nat_coords_trilinear_hex2
 // Duplicate functionality of nat_coords_trilinear_hex using hex_findpt
 //
 void nat_coords_trilinear_hex2( const CartVect hex[8], const CartVect& xyz, CartVect& ncoords, double etol )
 
 {
 const int ndim = 3;
 const int nverts = 8;
 const int vertMap[nverts] = { 0, 1, 3, 2, 4, 5, 7, 6 }; // Map from nat to lex ordering
 
 const int n = 2; // linear
 realType coords[ndim * nverts]; // buffer
 
 realType* xm[ndim];
 for( int i = 0; i < ndim; i++ )
 xm[i] = coords + i * nverts;
 
 // stuff hex into coords
 for( int i = 0; i < nverts; i++ )
 {
 realType vcoord[ndim];
 hex[i].get( vcoord );
 
 for( int d = 0; d < ndim; d++ )
 coords[d * nverts + vertMap[i]] = vcoord[d];
 }
 
 double dist = 0.0;
 ElemUtil::hex_findpt( xm, n, xyz, ncoords, dist );
 if( 3 * MOAB_POLY_EPS < dist )
 {
 // outside element, set extremal values to something outside range
 for( int j = 0; j < 3; j++ )
 {
 if( ncoords[j] < ( -1.0 - etol ) || ncoords[j] > ( 1.0 + etol ) ) ncoords[j] *= 10;
 }
 }
 }
 bool point_in_trilinear_hex( const CartVect* hex, const CartVect& xyz, double etol )
 {
 CartVect xi;
 return nat_coords_trilinear_hex( hex, xyz, xi, etol ) && ( fabs( xi[0] - 1.0 ) < etol ) &&
 ( fabs( xi[1] - 1.0 ) < etol ) && ( fabs( xi[2] - 1.0 ) < etol );
 }
 
 bool point_in_trilinear_hex( const CartVect* hex,
 const CartVect& xyz,
 const CartVect& box_min,
 const CartVect& box_max,
 double etol )
 {
 // all values scaled by 2 (eliminates 3 flops)
 const CartVect mid = box_max + box_min;
 const CartVect dim = box_max - box_min;
 const CartVect pt = 2 * xyz - mid;
 return ( fabs( pt[0] - dim[0] ) < etol ) && ( fabs( pt[1] - dim[1] ) < etol ) &&
 ( fabs( pt[2] - dim[2] ) < etol ) && point_in_trilinear_hex( hex, xyz, etol );
 }
 
 // Wrapper to James Lottes' findpt routines
 // hex_findpt
 // Find the parametric coordinates of an xyz point inside
 // a 3d hex spectral element with n nodes per dimension
 // xm: coordinates fields, value of x,y,z for each of then n*n*n gauss-lobatto nodes. Nodes are
 // in lexicographical order (x is fastest-changing) n: number of nodes per dimension -- n=2 for
 // a linear element xyz: input, point to find rst: output: parametric coords of xyz inside the
 // element. If xyz is outside the element, rst will be the coords of the closest point dist:
 // output: distance between xyz and the point with parametric coords rst
 /*extern "C"{
 #include "types.h"
 #include "poly.h"
 #include "tensor.h"
 #include "findpt.h"
 #include "extrafindpt.h"
 #include "errmem.h"
 }*/
 
 void hex_findpt( realType* xm[3], int n, CartVect xyz, CartVect& rst, double& dist )
 {
 
 // compute stuff that only depends on the order -- could be cached
 realType* z[3];
 lagrange_data ld[3];
 opt_data_3 data;
 
 // triplicates
 for( int d = 0; d < 3; d++ )
 {
 z[d] = tmalloc( realType, n );
 lobatto_nodes( z[d], n );
 lagrange_setup( &ld[d], z[d], n );
 }
 
 opt_alloc_3( &data, ld );
 
 // find nearest point
 realType x_star[3];
 xyz.get( x_star );
 
 realType r[3] = { 0, 0, 0 }; // initial guess for parametric coords
 unsigned c = opt_no_constraints_3;
 dist = opt_findpt_3( &data, (const realType**)xm, x_star, r, &c );
 // c tells us if we landed inside the element or exactly on a face, edge, or node
 
 // copy parametric coords back
 rst = r;
 
 // Clean-up (move to destructor if we decide to cache)
 opt_free_3( &data );
 for( int d = 0; d < 3; ++d )
 lagrange_free( &ld[d] );
 for( int d = 0; d < 3; ++d )
 free( z[d] );
 }
 
 // hex_eval
 // Evaluate a field in a 3d hex spectral element with n nodes per dimension, at some given
 // parametric coordinates field: field values for each of then n*n*n gauss-lobatto nodes. Nodes
 // are in lexicographical order (x is fastest-changing) n: number of nodes per dimension -- n=2
 // for a linear element rst: input: parametric coords of the point where we want to evaluate the
 // field value: output: value of field at rst
 
 void hex_eval( realType* field, int n, CartVect rstCartVec, double& value )
 {
 int d;
 realType rst[3];
 rstCartVec.get( rst );
 
 // can cache stuff below
 lagrange_data ld[3];
 realType* z[3];
 for( d = 0; d < 3; ++d )
 {
 z[d] = tmalloc( realType, n );
 lobatto_nodes( z[d], n );
 lagrange_setup( &ld[d], z[d], n );
 }
 
 // cut and paste -- see findpt.c
 const unsigned nf = n * n, ne = n, nw = 2 * n * n + 3 * n;
 realType* od_work = tmalloc( realType, 6 * nf + 9 * ne + nw );
 
 // piece that we shouldn't want to cache
 for( d = 0; d < 3; d++ )
 {
 lagrange_0( &ld[d], rst[d] );
 }
 
 value = tensor_i3( ld[0].J, ld[0].n, ld[1].J, ld[1].n, ld[2].J, ld[2].n, field, od_work );
 
 // all this could be cached
 for( d = 0; d < 3; d++ )
 {
 free( z[d] );
 lagrange_free( &ld[d] );
 }
 free( od_work );
 }
 
 // Gaussian quadrature points for a trilinear hex element.
 // Five 2d arrays are defined.
 // One for the single gaussian point solution, 2 point solution,
 // 3 point solution, 4 point solution and 5 point solution.
 // There are 2 columns, one for Weights and the other for Locations
 // Weight Location
 
 const double gauss_1[1][2] = { { 2.0, 0.0 } };
 
 const double gauss_2[2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } };
 
 const double gauss_3[3][2] = { { 0.5555555555, -0.7745966692 },
 { 0.8888888888, 0.0 },
 { 0.5555555555, 0.7745966692 } };
 
 const double gauss_4[4][2] = { { 0.3478548451, -0.8611363116 },
 { 0.6521451549, -0.3399810436 },
 { 0.6521451549, 0.3399810436 },
 { 0.3478548451, 0.8611363116 } };
 
 const double gauss_5[5][2] = { { 0.2369268851, -0.9061798459 },
 { 0.4786286705, -0.5384693101 },
 { 0.5688888889, 0.0 },
 { 0.4786286705, 0.5384693101 },
 { 0.2369268851, 0.9061798459 } };
 
 // Function to integrate the field defined by field_fn function
 // over the volume of the trilinear hex defined by the hex_corners
 
 bool integrate_trilinear_hex( const CartVect* hex_corners, double* corner_fields, double& field_val, int num_pts )
 {
 // Create the LinearHexMap object using the hex_corners array of CartVects
 LinearHexMap hex( hex_corners );
 
 // Use the correct table of points and locations based on the num_pts parameter
 const double( *g_pts )[2] = 0;
 switch( num_pts )
 {
 case 1:
 g_pts = gauss_1;
 break;
 
 case 2:
 g_pts = gauss_2;
 break;
 
 case 3:
 g_pts = gauss_3;
 break;
 
 case 4:
 g_pts = gauss_4;
 break;
 
 case 5:
 g_pts = gauss_5;
 break;
 
 default: // value out of range
 return false;
 }
 
 // Test code - print Gaussian Quadrature data
 if( debug )
 {
 for( int r = 0; r < num_pts; r++ )
 for( int c = 0; c < 2; c++ )
 std::cout << "g_pts[" << r << "][" << c << "]=" << g_pts[r][c] << std::endl;
 }
 // End Test code
 
 double soln = 0.0;
 
 for( int i = 0; i < num_pts; i++ )
 { // Loop for xi direction
 double w_i = g_pts[i][0];
 double xi_i = g_pts[i][1];
 for( int j = 0; j < num_pts; j++ )
 { // Loop for eta direction
 double w_j = g_pts[j][0];
 double eta_j = g_pts[j][1];
 for( int k = 0; k < num_pts; k++ )
 { // Loop for zeta direction
 double w_k = g_pts[k][0];
 double zeta_k = g_pts[k][1];
 
 // Calculate the "realType" space point given the "normal" point
 CartVect normal_pt( xi_i, eta_j, zeta_k );
 
 // Calculate the value of F(x(xi,eta,zeta),y(xi,eta,zeta),z(xi,eta,zeta)
 double field = hex.evaluate_scalar_field( normal_pt, corner_fields );
 
 // Calculate the Jacobian for this "normal" point and its determinant
 Matrix3 J = hex.jacobian( normal_pt );
 double det = J.determinant();
 
 // Calculate integral and update the solution
 soln = soln + ( w_i * w_j * w_k * field * det );
 }
 }
 }
 
 // Set the output parameter
 field_val = soln;
 
 return true;
 }
 
} // namespace ElemUtil
 
namespace Element
{
 
 Map::~Map() {}
 
 inline const std::vector< CartVect >& Map::get_vertices()
 {
 return this->vertex;
 }
 //
 void Map::set_vertices( const std::vector< CartVect >& v )
 {
 if( v.size() != this->vertex.size() )
 {
 throw ArgError();
 }
 this->vertex = v;
 }
 
 bool Map::inside_box( const CartVect& xi, double& tol ) const
 {
 // bail out early, before doing an expensive NR iteration
 // compute box
 BoundBox box( this->vertex );
 return box.contains_point( xi.array(), tol );
 }
 
 //
 CartVect Map::ievaluate( const CartVect& x, double tol, const CartVect& x0 ) const
 {
 // TODO: should differentiate between epsilons used for
 // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
 // right now, fix the tolerance used for NR
 tol = 1.0e-10;
 const double error_tol_sqr = tol * tol;
 double det;
 CartVect xi = x0;
 CartVect delta = evaluate( xi ) - x;
 Matrix3 J;
 
 int iters = 0;
 while( delta % delta > error_tol_sqr )
 {
 if( ++iters > 10 ) throw Map::EvaluationError( x, vertex );
 
 J = jacobian( xi );
 det = J.determinant();
 if( det < std::numeric_limits< double >::epsilon() ) throw Map::EvaluationError( x, vertex );
 xi -= J.inverse() * delta;
 delta = evaluate( xi ) - x;
 }
 return xi;
 } // Map::ievaluate()
 
 SphericalQuad::SphericalQuad( const std::vector< CartVect >& vertices ) : LinearQuad( vertices )
 {
 // project the vertices to the plane tangent at first vertex
 v1 = vertex[0]; // member data
 double v1v1 = v1 % v1; // this is 1, in general, for unit sphere meshes
 for( int j = 1; j < 4; j++ )
 {
 // first, bring all vertices in the gnomonic plane
 // the new vertex will intersect the plane at vnew
 // so that (vnew-v1)%v1 is 0 ( vnew is in the tangent plane, i.e. normal to v1 )
 // pos is the old position of the vertex, and it is in general on the sphere
 // vnew = alfa*pos; (alfa*pos-v1)%v1 = 0 <=> alfa*(pos%v1)=v1%v1 <=> alfa =
 // v1v1/(pos%v1)
 // <=> vnew = ( v1v1/(pos%v1) )*pos
 CartVect vnew = v1v1 / ( vertex[j] % v1 ) * vertex[j];
 vertex[j] = vnew;
 }
 // will compute a transf matrix, such that a new point will be transformed with
 // newpos = transf * (vnew-v1), and it will be a point in the 2d plane
 // the transformation matrix will be oriented in such a way that orientation will be
 // positive
 CartVect vx = vertex[1] - v1; // this will become Ox axis
 // z axis will be along v1, in such a way that orientation of the quad is positive
 // look at the first 2 edges
 CartVect vz = vx * ( vertex[2] - vertex[1] );
 vz = vz / vz.length();
 
 vx = vx / vx.length();
 
 CartVect vy = vz * vx;
 transf = Matrix3( vx[0], vx[1], vx[2], vy[0], vy[1], vy[2], vz[0], vz[1], vz[2] );
 vertex[0] = CartVect( 0. );
 for( int j = 1; j < 4; j++ )
 vertex[j] = transf * ( vertex[j] - v1 );
 }
 
 CartVect SphericalQuad::ievaluate( const CartVect& x, double tol, const CartVect& x0 ) const
 {
 // project to the plane tangent at first vertex (gnomonic projection)
 double v1v1 = v1 % v1;
 CartVect vnew = v1v1 / ( x % v1 ) * x; // so that (vnew-v1)%v1 is 0
 vnew = transf * ( vnew - v1 );
 // det will be positive now
 return Map::ievaluate( vnew, tol, x0 );
 }
 
 bool SphericalQuad::inside_box( const CartVect& pos, double& tol ) const
 {
 // project to the plane tangent at first vertex
 // CartVect v1=vertex[0];
 double v1v1 = v1 % v1;
 CartVect vnew = v1v1 / ( pos % v1 ) * pos; // so that (x-v1)%v1 is 0
 vnew = transf * ( vnew - v1 );
 return Map::inside_box( vnew, tol );
 }
 
 const double LinearTri::corner[3][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 } };
 
 LinearTri::LinearTri() : Map( 0 ), det_T( 0.0 ), det_T_inverse( 0.0 ) {} // LinearTri::LinearTri()
 
 LinearTri::~LinearTri() {}
 
 void LinearTri::set_vertices( const std::vector< CartVect >& v )
 {
 this->Map::set_vertices( v );
 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,
 v[1][2] - v[0][2], v[2][2] - v[0][2], 1 );
 this->T_inverse = this->T.inverse();
 this->det_T = this->T.determinant();
 this->det_T_inverse = ( this->det_T < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / this->det_T );
 } // LinearTri::set_vertices()
 
 bool LinearTri::inside_nat_space( const CartVect& xi, double& tol ) const
 {
 // linear tri space is a triangle with vertices (0,0,0), (1,0,0), (0,1,0)
 // first check if outside bigger box, then below the line x+y=1
 return ( xi[0] >= -tol ) && ( xi[1] >= -tol ) && ( xi[2] >= -tol ) && ( xi[2] <= tol ) &&
 ( xi[0] + xi[1] < 1.0 + tol );
 }
 CartVect LinearTri::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const
 {
 return this->T_inverse * ( x - this->vertex[0] );
 } // LinearTri::ievaluate
 
 double LinearTri::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
 {
 double f0 = field_vertex_value[0];
 double f = f0;
 for( unsigned i = 1; i < 3; ++i )
 {
 f += ( field_vertex_value[i] - f0 ) * xi[i - 1];
 }
 return f;
 } // LinearTri::evaluate_scalar_field()
 
 double LinearTri::integrate_scalar_field( const double* field_vertex_values ) const
 {
 double I( 0.0 );
 for( unsigned int i = 0; i < 3; ++i )
 {
 I += field_vertex_values[i];
 }
 I *= this->det_T / 6.0; // TODO
 return I;
 } // LinearTri::integrate_scalar_field()
 
 SphericalTri::SphericalTri( const std::vector< CartVect >& vertices )
 {
 vertex.resize( vertices.size() );
 vertex = vertices;
 // project the vertices to the plane tangent at first vertex
 v1 = vertex[0]; // member data
 double v1v1 = v1 % v1; // this is 1, in general, for unit sphere meshes
 for( int j = 1; j < 3; j++ )
 {
 // first, bring all vertices in the gnomonic plane
 // the new vertex will intersect the plane at vnew
 // so that (vnew-v1)%v1 is 0 ( vnew is in the tangent plane, i.e. normal to v1 )
 // pos is the old position of the vertex, and it is in general on the sphere
 // vnew = alfa*pos; (alfa*pos-v1)%v1 = 0 <=> alfa*(pos%v1)=v1%v1 <=> alfa =
 // v1v1/(pos%v1)
 // <=> vnew = ( v1v1/(pos%v1) )*pos
 CartVect vnew = v1v1 / ( vertex[j] % v1 ) * vertex[j];
 vertex[j] = vnew;
 }
 // will compute a transf matrix, such that a new point will be transformed with
 // newpos = transf * (vnew-v1), and it will be a point in the 2d plane
 // the transformation matrix will be oriented in such a way that orientation will be
 // positive
 CartVect vx = vertex[1] - v1; // this will become Ox axis
 // z axis will be along v1, in such a way that orientation of the quad is positive
 // look at the first 2 edges
 CartVect vz = vx * ( vertex[2] - vertex[1] );
 vz = vz / vz.length();
 
 vx = vx / vx.length();
 
 CartVect vy = vz * vx;
 transf = Matrix3( vx[0], vx[1], vx[2], vy[0], vy[1], vy[2], vz[0], vz[1], vz[2] );
 vertex[0] = CartVect( 0. );
 for( int j = 1; j < 3; j++ )
 vertex[j] = transf * ( vertex[j] - v1 );
 
 LinearTri::set_vertices( vertex );
 }
 
 CartVect SphericalTri::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const
 {
 // project to the plane tangent at first vertex (gnomonic projection)
 double v1v1 = v1 % v1;
 CartVect vnew = v1v1 / ( x % v1 ) * x; // so that (vnew-v1)%v1 is 0
 vnew = transf * ( vnew - v1 );
 // det will be positive now
 return LinearTri::ievaluate( vnew );
 }
 
 bool SphericalTri::inside_box( const CartVect& pos, double& tol ) const
 {
 // project to the plane tangent at first vertex
 // CartVect v1=vertex[0];
 double v1v1 = v1 % v1;
 CartVect vnew = v1v1 / ( pos % v1 ) * pos; // so that (x-v1)%v1 is 0
 vnew = transf * ( vnew - v1 );
 return Map::inside_box( vnew, tol );
 }
 
 // filescope for static member data that is cached
 const double LinearEdge::corner[2][3] = { { -1, 0, 0 }, { 1, 0, 0 } };
 
 LinearEdge::LinearEdge() : Map( 0 ) {} // LinearEdge::LinearEdge()
 
 /* For each point, its weight and location are stored as an array.
 Hence, the inner dimension is 2, the outer dimension is gauss_count.
 We use a one-point Gaussian quadrature, since it integrates linear functions exactly.
 */
 const double LinearEdge::gauss[1][2] = { { 2.0, 0.0 } };
 
 CartVect LinearEdge::evaluate( const CartVect& xi ) const
 {
 CartVect x( 0.0 );
 for( unsigned i = 0; i < LinearEdge::corner_count; ++i )
 {
 const double N_i = ( 1.0 + xi[0] * corner[i][0] );
 x += N_i * this->vertex[i];
 }
 x /= LinearEdge::corner_count;
 return x;
 } // LinearEdge::evaluate
 
 Matrix3 LinearEdge::jacobian( const CartVect& xi ) const
 {
 Matrix3 J( 0.0 );
 for( unsigned i = 0; i < LinearEdge::corner_count; ++i )
 {
 const double xi_p = 1.0 + xi[0] * corner[i][0];
 const double dNi_dxi = corner[i][0] * xi_p;
 J( 0, 0 ) += dNi_dxi * vertex[i][0];
 }
 J( 1, 1 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */
 J( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */
 J /= LinearEdge::corner_count;
 return J;
 } // LinearEdge::jacobian()
 
 double LinearEdge::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
 {
 double f( 0.0 );
 for( unsigned i = 0; i < LinearEdge::corner_count; ++i )
 {
 const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1.0 + xi[1] * corner[i][1] );
 f += N_i * field_vertex_value[i];
 }
 f /= LinearEdge::corner_count;
 return f;
 } // LinearEdge::evaluate_scalar_field()
 
 double LinearEdge::integrate_scalar_field( const double* field_vertex_values ) const
 {
 double I( 0.0 );
 for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 )
 {
 double x1 = this->gauss[j1][1];
 double w1 = this->gauss[j1][0];
 CartVect x( x1, 0.0, 0.0 );
 I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * this->det_jacobian( x );
 }
 return I;
 } // LinearEdge::integrate_scalar_field()
 
 bool LinearEdge::inside_nat_space( const CartVect& xi, double& tol ) const
 {
 // just look at the box+tol here
 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol );
 }
 
 const double LinearHex::corner[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
 { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } };
 
 LinearHex::LinearHex() : Map( 0 ) {} // LinearHex::LinearHex()
 
 LinearHex::~LinearHex() {}
 /* For each point, its weight and location are stored as an array.
 Hence, the inner dimension is 2, the outer dimension is gauss_count.
 We use a one-point Gaussian quadrature, since it integrates linear functions exactly.
 */
 // const double LinearHex::gauss[1][2] = { { 2.0, 0.0 } };
 const double LinearHex::gauss[2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } };
 // const double LinearHex::gauss[4][2] = { { 0.3478548451, -0.8611363116 },
 // { 0.6521451549, -0.3399810436 },
 // { 0.6521451549, 0.3399810436 },
 // { 0.3478548451, 0.8611363116 } };
 
 CartVect LinearHex::evaluate( const CartVect& xi ) const
 {
 CartVect x( 0.0 );
 for( unsigned i = 0; i < 8; ++i )
 {
 const double N_i =
 ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] ) * ( 1 + xi[2] * corner[i][2] );
 x += N_i * this->vertex[i];
 }
 x *= 0.125;
 return x;
 } // LinearHex::evaluate
 
 Matrix3 LinearHex::jacobian( const CartVect& xi ) const
 {
 Matrix3 J( 0.0 );
 for( unsigned i = 0; i < 8; ++i )
 {
 const double xi_p = 1 + xi[0] * corner[i][0];
 const double eta_p = 1 + xi[1] * corner[i][1];
 const double zeta_p = 1 + xi[2] * corner[i][2];
 const double dNi_dxi = corner[i][0] * eta_p * zeta_p;
 const double dNi_deta = corner[i][1] * xi_p * zeta_p;
 const double dNi_dzeta = corner[i][2] * xi_p * eta_p;
 J( 0, 0 ) += dNi_dxi * vertex[i][0];
 J( 1, 0 ) += dNi_dxi * vertex[i][1];
 J( 2, 0 ) += dNi_dxi * vertex[i][2];
 J( 0, 1 ) += dNi_deta * vertex[i][0];
 J( 1, 1 ) += dNi_deta * vertex[i][1];
 J( 2, 1 ) += dNi_deta * vertex[i][2];
 J( 0, 2 ) += dNi_dzeta * vertex[i][0];
 J( 1, 2 ) += dNi_dzeta * vertex[i][1];
 J( 2, 2 ) += dNi_dzeta * vertex[i][2];
 }
 return J *= 0.125;
 } // LinearHex::jacobian()
 
 double LinearHex::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
 {
 double f( 0.0 );
 for( unsigned i = 0; i < 8; ++i )
 {
 const double N_i =
 ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] ) * ( 1 + xi[2] * corner[i][2] );
 f += N_i * field_vertex_value[i];
 }
 f *= 0.125;
 return f;
 } // LinearHex::evaluate_scalar_field()
 
 double LinearHex::integrate_scalar_field( const double* field_vertex_values ) const
 {
 double I( 0.0 );
 for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 )
 {
 double x1 = this->gauss[j1][1];
 double w1 = this->gauss[j1][0];
 for( unsigned int j2 = 0; j2 < this->gauss_count; ++j2 )
 {
 double x2 = this->gauss[j2][1];
 double w2 = this->gauss[j2][0];
 for( unsigned int j3 = 0; j3 < this->gauss_count; ++j3 )
 {
 double x3 = this->gauss[j3][1];
 double w3 = this->gauss[j3][0];
 CartVect x( x1, x2, x3 );
 I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * w2 * w3 * this->det_jacobian( x );
 }
 }
 }
 return I;
 } // LinearHex::integrate_scalar_field()
 
 bool LinearHex::inside_nat_space( const CartVect& xi, double& tol ) const
 {
 // just look at the box+tol here
 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) &&
 ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol );
 }
 
 // those are not just the corners, but for simplicity, keep this name
 //
 const int QuadraticHex::corner[27][3] = {
 { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, // corner nodes: 0-7
 { -1, 1, -1 }, // mid-edge nodes: 8-19
 { -1, -1, 1 }, // center-face nodes 20-25 center node 26
 { 1, -1, 1 }, //
 { 1, 1, 1 }, { -1, 1, 1 }, // 4 ----- 19 ----- 7
 { 0, -1, -1 }, // . | . |
 { 1, 0, -1 }, // 16 25 18 |
 { 0, 1, -1 }, // . | . |
 { -1, 0, -1 }, // 5 ----- 17 ----- 6 |
 { -1, -1, 0 }, // | 12 | 23 15
 { 1, -1, 0 }, // | | |
 { 1, 1, 0 }, // | 20 | 26 | 22 |
 { -1, 1, 0 }, // | | |
 { 0, -1, 1 }, // 13 21 | 14 |
 { 1, 0, 1 }, // | 0 ----- 11 ----- 3
 { 0, 1, 1 }, // | . | .
 { -1, 0, 1 }, // | 8 24 | 10
 { 0, -1, 0 }, // | . | .
 { 1, 0, 0 }, // 1 ----- 9 ----- 2
 { 0, 1, 0 }, //
 { -1, 0, 0 }, { 0, 0, -1 }, { 0, 0, 1 }, { 0, 0, 0 } };
 // QuadraticHex::QuadraticHex(const std::vector<CartVect>& vertices) : Map(vertices){};
 QuadraticHex::QuadraticHex() : Map( 0 ) {}
 
 QuadraticHex::~QuadraticHex() {}
 double SH( const int i, const double xi )
 {
 switch( i )
 {
 case -1:
 return ( xi * xi - xi ) / 2;
 case 0:
 return 1 - xi * xi;
 case 1:
 return ( xi * xi + xi ) / 2;
 default:
 return 0.;
 }
 }
 double DSH( const int i, const double xi )
 {
 switch( i )
 {
 case -1:
 return xi - 0.5;
 case 0:
 return -2 * xi;
 case 1:
 return xi + 0.5;
 default:
 return 0.;
 }
 }
 
 CartVect QuadraticHex::evaluate( const CartVect& xi ) const
 {
 
 CartVect x( 0.0 );
 for( int i = 0; i < 27; i++ )
 {
 const double sh = SH( corner[i][0], xi[0] ) * SH( corner[i][1], xi[1] ) * SH( corner[i][2], xi[2] );
 x += sh * vertex[i];
 }
 
 return x;
 }
 
 bool QuadraticHex::inside_nat_space( const CartVect& xi, double& tol ) const
 { // just look at the box+tol here
 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) &&
 ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol );
 }
 
 Matrix3 QuadraticHex::jacobian( const CartVect& xi ) const
 {
 Matrix3 J( 0.0 );
 for( int i = 0; i < 27; i++ )
 {
 const double sh[3] = { SH( corner[i][0], xi[0] ), SH( corner[i][1], xi[1] ), SH( corner[i][2], xi[2] ) };
 const double dsh[3] = { DSH( corner[i][0], xi[0] ), DSH( corner[i][1], xi[1] ),
 DSH( corner[i][2], xi[2] ) };
 
 for( int j = 0; j < 3; j++ )
 {
 J( j, 0 ) += dsh[0] * sh[1] * sh[2] * vertex[i][j]; // dxj/dr first column
 J( j, 1 ) += sh[0] * dsh[1] * sh[2] * vertex[i][j]; // dxj/ds
 J( j, 2 ) += sh[0] * sh[1] * dsh[2] * vertex[i][j]; // dxj/dt
 }
 }
 
 return J;
 }
 double QuadraticHex::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const
 {
 double x = 0.0;
 for( int i = 0; i < 27; i++ )
 {
 const double sh = SH( corner[i][0], xi[0] ) * SH( corner[i][1], xi[1] ) * SH( corner[i][2], xi[2] );
 x += sh * field_vertex_values[i];
 }
 
 return x;
 }
 double QuadraticHex::integrate_scalar_field( const double* /*field_vertex_values*/ ) const
 {
 return 0.; // TODO: gaussian integration , probably 2x2x2
 }
 
 const double LinearTet::corner[4][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
 
 LinearTet::LinearTet() : Map( 0 ), det_T( 0.0 ), det_T_inverse( 0.0 ) {} // LinearTet::LinearTet()
 
 LinearTet::~LinearTet() {}
 
 void LinearTet::set_vertices( const std::vector< CartVect >& v )
 {
 this->Map::set_vertices( v );
 this->T =
 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],
 v[3][1] - v[0][1], v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] );
 this->T_inverse = this->T.inverse();
 this->det_T = this->T.determinant();
 this->det_T_inverse = ( this->det_T < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / this->det_T );
 } // LinearTet::set_vertices()
 
 double LinearTet::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
 {
 double f0 = field_vertex_value[0];
 double f = f0;
 for( unsigned i = 1; i < 4; ++i )
 {
 f += ( field_vertex_value[i] - f0 ) * xi[i - 1];
 }
 return f;
 } // LinearTet::evaluate_scalar_field()
 
 CartVect LinearTet::ievaluate( const CartVect& x, double /*tol*/, const CartVect& /*x0*/ ) const
 {
 return this->T_inverse * ( x - this->vertex[0] );
 } // LinearTet::ievaluate
 
 double LinearTet::integrate_scalar_field( const double* field_vertex_values ) const
 {
 double I( 0.0 );
 for( unsigned int i = 0; i < 4; ++i )
 {
 I += field_vertex_values[i];
 }
 I *= this->det_T / 24.0;
 return I;
 } // LinearTet::integrate_scalar_field()
 bool LinearTet::inside_nat_space( const CartVect& xi, double& tol ) const
 {
 // linear tet space is a tetra with vertices (0,0,0), (1,0,0), (0,1,0), (0, 0, 1)
 // first check if outside bigger box, then below the plane x+y+z=1
 return ( xi[0] >= -tol ) && ( xi[1] >= -tol ) && ( xi[2] >= -tol ) && ( xi[0] + xi[1] + xi[2] < 1.0 + tol );
 }
 // SpectralHex
 
 // filescope for static member data that is cached
 int SpectralHex::_n;
 realType* SpectralHex::_z[3];
 lagrange_data SpectralHex::_ld[3];
 opt_data_3 SpectralHex::_data;
 realType* SpectralHex::_odwork;
 
 bool SpectralHex::_init = false;
 
 SpectralHex::SpectralHex() : Map( 0 )
 {
 _xyz[0] = _xyz[1] = _xyz[2] = NULL;
 }
 // the preferred constructor takes pointers to GL blocked positions
 SpectralHex::SpectralHex( int order, double* x, double* y, double* z ) : Map( 0 )
 {
 Init( order );
 _xyz[0] = x;
 _xyz[1] = y;
 _xyz[2] = z;
 }
 SpectralHex::SpectralHex( int order ) : Map( 0 )
 {
 Init( order );
 _xyz[0] = _xyz[1] = _xyz[2] = NULL;
 }
 SpectralHex::~SpectralHex()
 {
 if( _init ) freedata();
 _init = false;
 }
 void SpectralHex::Init( int order )
 {
 if( _init && _n == order ) return;
 if( _init && _n != order )
 {
 // TODO: free data cached
 freedata();
 }
 // compute stuff that depends only on order
 _init = true;
 _n = order;
 // triplicates! n is the same in all directions !!!
 for( int d = 0; d < 3; d++ )
 {
 _z[d] = tmalloc( realType, _n );
 lobatto_nodes( _z[d], _n );
 lagrange_setup( &_ld[d], _z[d], _n );
 }
 opt_alloc_3( &_data, _ld );
 
 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
 _odwork = tmalloc( realType, 6 * nf + 9 * ne + nw );
 }
 void SpectralHex::freedata()
 {
 for( int d = 0; d < 3; d++ )
 {
 free( _z[d] );
 lagrange_free( &_ld[d] );
 }
 opt_free_3( &_data );
 free( _odwork );
 }
 
 void SpectralHex::set_gl_points( double* x, double* y, double* z )
 {
 _xyz[0] = x;
 _xyz[1] = y;
 _xyz[2] = z;
 }
 CartVect SpectralHex::evaluate( const CartVect& xi ) const
 {
 // piece that we shouldn't want to cache
 int d = 0;
 for( d = 0; d < 3; d++ )
 {
 lagrange_0( &_ld[d], xi[d] );
 }
 CartVect result;
 for( d = 0; d < 3; d++ )
 {
 result[d] = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n,
 _xyz[d], // this is the "field"
 _odwork );
 }
 return result;
 }
 // replicate the functionality of hex_findpt
 CartVect SpectralHex::ievaluate( CartVect const& xyz, double tol, const CartVect& x0 ) const
 {
 // find nearest point
 realType x_star[3];
 xyz.get( x_star );
 
 realType r[3] = { 0, 0, 0 }; // initial guess for parametric coords
 x0.get( r );
 unsigned c = opt_no_constraints_3;
 realType dist = opt_findpt_3( &_data, (const realType**)_xyz, x_star, r, &c );
 // if it did not converge, get out with throw...
 if( dist > 10 * tol ) // outside the element
 {
 std::vector< CartVect > dummy;
 throw Map::EvaluationError( xyz, dummy );
 }
 // c tells us if we landed inside the element or exactly on a face, edge, or node
 // also, dist shows the distance to the computed point.
 // copy parametric coords back
 return CartVect( r );
 }
 Matrix3 SpectralHex::jacobian( const CartVect& xi ) const
 {
 realType x_i[3];
 xi.get( x_i );
 // set the positions of GL nodes, before evaluations
 _data.elx[0] = _xyz[0];
 _data.elx[1] = _xyz[1];
 _data.elx[2] = _xyz[2];
 opt_vol_set_intp_3( &_data, x_i );
 Matrix3 J( 0. );
 // it is organized differently
 J( 0, 0 ) = _data.jac[0]; // dx/dr
 J( 0, 1 ) = _data.jac[1]; // dx/ds
 J( 0, 2 ) = _data.jac[2]; // dx/dt
 J( 1, 0 ) = _data.jac[3]; // dy/dr
 J( 1, 1 ) = _data.jac[4]; // dy/ds
 J( 1, 2 ) = _data.jac[5]; // dy/dt
 J( 2, 0 ) = _data.jac[6]; // dz/dr
 J( 2, 1 ) = _data.jac[7]; // dz/ds
 J( 2, 2 ) = _data.jac[8]; // dz/dt
 return J;
 }
 double SpectralHex::evaluate_scalar_field( const CartVect& xi, const double* field ) const
 {
 // piece that we shouldn't want to cache
 int d;
 for( d = 0; d < 3; d++ )
 {
 lagrange_0( &_ld[d], xi[d] );
 }
 
 double value = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork );
 return value;
 }
 double SpectralHex::integrate_scalar_field( const double* field_vertex_values ) const
 {
 // set the position of GL points
 // set the positions of GL nodes, before evaluations
 _data.elx[0] = _xyz[0];
 _data.elx[1] = _xyz[1];
 _data.elx[2] = _xyz[2];
 double xi[3];
 // triple loop; the most inner loop is in r direction, then s, then t
 double integral = 0.;
 // double volume = 0;
 int index = 0; // used fr the inner loop
 for( int k = 0; k < _n; k++ )
 {
 xi[2] = _ld[2].z[k];
 // double wk= _ld[2].w[k];
 for( int j = 0; j < _n; j++ )
 {
 xi[1] = _ld[1].z[j];
 // double wj= _ld[1].w[j];
 for( int i = 0; i < _n; i++ )
 {
 xi[0] = _ld[0].z[i];
 // double wi= _ld[0].w[i];
 opt_vol_set_intp_3( &_data, xi );
 double wk = _ld[2].J[k];
 double wj = _ld[1].J[j];
 double wi = _ld[0].J[i];
 Matrix3 J( 0. );
 // it is organized differently
 J( 0, 0 ) = _data.jac[0]; // dx/dr
 J( 0, 1 ) = _data.jac[1]; // dx/ds
 J( 0, 2 ) = _data.jac[2]; // dx/dt
 J( 1, 0 ) = _data.jac[3]; // dy/dr
 J( 1, 1 ) = _data.jac[4]; // dy/ds
 J( 1, 2 ) = _data.jac[5]; // dy/dt
 J( 2, 0 ) = _data.jac[6]; // dz/dr
 J( 2, 1 ) = _data.jac[7]; // dz/ds
 J( 2, 2 ) = _data.jac[8]; // dz/dt
 double bm = wk * wj * wi * J.determinant();
 integral += bm * field_vertex_values[index++];
 // volume +=bm;
 }
 }
 }
 // std::cout << "volume: " << volume << "\n";
 return integral;
 }
 // this is the same as a linear hex, although we should not need it
 bool SpectralHex::inside_nat_space( const CartVect& xi, double& tol ) const
 {
 // just look at the box+tol here
 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol ) &&
 ( xi[2] >= -1. - tol ) && ( xi[2] <= 1. + tol );
 }
 
 // SpectralHex
 
 // filescope for static member data that is cached
 const double LinearQuad::corner[4][3] = { { -1, -1, 0 }, { 1, -1, 0 }, { 1, 1, 0 }, { -1, 1, 0 } };
 
 LinearQuad::LinearQuad() : Map( 0 ) {} // LinearQuad::LinearQuad()
 
 LinearQuad::~LinearQuad() {}
 
 /* For each point, its weight and location are stored as an array.
 Hence, the inner dimension is 2, the outer dimension is gauss_count.
 We use a one-point Gaussian quadrature, since it integrates linear functions exactly.
 */
 const double LinearQuad::gauss[1][2] = { { 2.0, 0.0 } };
 
 CartVect LinearQuad::evaluate( const CartVect& xi ) const
 {
 CartVect x( 0.0 );
 for( unsigned i = 0; i < LinearQuad::corner_count; ++i )
 {
 const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] );
 x += N_i * this->vertex[i];
 }
 x /= LinearQuad::corner_count;
 return x;
 } // LinearQuad::evaluate
 
 Matrix3 LinearQuad::jacobian( const CartVect& xi ) const
 {
 // this basically ignores the z component: xi[2] or vertex[][2]
 Matrix3 J( 0.0 );
 for( unsigned i = 0; i < LinearQuad::corner_count; ++i )
 {
 const double xi_p = 1 + xi[0] * corner[i][0];
 const double eta_p = 1 + xi[1] * corner[i][1];
 const double dNi_dxi = corner[i][0] * eta_p;
 const double dNi_deta = corner[i][1] * xi_p;
 J( 0, 0 ) += dNi_dxi * vertex[i][0];
 J( 1, 0 ) += dNi_dxi * vertex[i][1];
 J( 0, 1 ) += dNi_deta * vertex[i][0];
 J( 1, 1 ) += dNi_deta * vertex[i][1];
 }
 J( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */
 J /= LinearQuad::corner_count;
 return J;
 } // LinearQuad::jacobian()
 
 double LinearQuad::evaluate_scalar_field( const CartVect& xi, const double* field_vertex_value ) const
 {
 double f( 0.0 );
 for( unsigned i = 0; i < LinearQuad::corner_count; ++i )
 {
 const double N_i = ( 1 + xi[0] * corner[i][0] ) * ( 1 + xi[1] * corner[i][1] );
 f += N_i * field_vertex_value[i];
 }
 f /= LinearQuad::corner_count;
 return f;
 } // LinearQuad::evaluate_scalar_field()
 
 double LinearQuad::integrate_scalar_field( const double* field_vertex_values ) const
 {
 double I( 0.0 );
 for( unsigned int j1 = 0; j1 < this->gauss_count; ++j1 )
 {
 double x1 = this->gauss[j1][1];
 double w1 = this->gauss[j1][0];
 for( unsigned int j2 = 0; j2 < this->gauss_count; ++j2 )
 {
 double x2 = this->gauss[j2][1];
 double w2 = this->gauss[j2][0];
 CartVect x( x1, x2, 0.0 );
 I += this->evaluate_scalar_field( x, field_vertex_values ) * w1 * w2 * this->det_jacobian( x );
 }
 }
 return I;
 } // LinearQuad::integrate_scalar_field()
 
 bool LinearQuad::inside_nat_space( const CartVect& xi, double& tol ) const
 {
 // just look at the box+tol here
 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol );
 }
 
 // filescope for static member data that is cached
 int SpectralQuad::_n;
 realType* SpectralQuad::_z[2];
 lagrange_data SpectralQuad::_ld[2];
 opt_data_2 SpectralQuad::_data;
 realType* SpectralQuad::_odwork;
 realType* SpectralQuad::_glpoints;
 bool SpectralQuad::_init = false;
 
 SpectralQuad::SpectralQuad() : Map( 0 )
 {
 _xyz[0] = _xyz[1] = _xyz[2] = NULL;
 }
 // the preferred constructor takes pointers to GL blocked positions
 SpectralQuad::SpectralQuad( int order, double* x, double* y, double* z ) : Map( 0 )
 {
 Init( order );
 _xyz[0] = x;
 _xyz[1] = y;
 _xyz[2] = z;
 }
 SpectralQuad::SpectralQuad( int order ) : Map( 4 )
 {
 Init( order );
 _xyz[0] = _xyz[1] = _xyz[2] = NULL;
 }
 SpectralQuad::~SpectralQuad()
 {
 if( _init ) freedata();
 _init = false;
 }
 void SpectralQuad::Init( int order )
 {
 if( _init && _n == order ) return;
 if( _init && _n != order )
 {
 // TODO: free data cached
 freedata();
 }
 // compute stuff that depends only on order
 _init = true;
 _n = order;
 // duplicates! n is the same in all directions !!!
 for( int d = 0; d < 2; d++ )
 {
 _z[d] = tmalloc( realType, _n );
 lobatto_nodes( _z[d], _n );
 lagrange_setup( &_ld[d], _z[d], _n );
 }
 opt_alloc_2( &_data, _ld );
 
 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
 _odwork = tmalloc( realType, 6 * nf + 9 * ne + nw );
 _glpoints = tmalloc( realType, 3 * nf );
 }
 
 void SpectralQuad::freedata()
 {
 for( int d = 0; d < 2; d++ )
 {
 free( _z[d] );
 lagrange_free( &_ld[d] );
 }
 opt_free_2( &_data );
 free( _odwork );
 free( _glpoints );
 }
 
 void SpectralQuad::set_gl_points( double* x, double* y, double* z )
 {
 _xyz[0] = x;
 _xyz[1] = y;
 _xyz[2] = z;
 }
 CartVect SpectralQuad::evaluate( const CartVect& xi ) const
 {
 // piece that we shouldn't want to cache
 int d = 0;
 for( d = 0; d < 2; d++ )
 {
 lagrange_0( &_ld[d], xi[d] );
 }
 CartVect result;
 for( d = 0; d < 3; d++ )
 {
 result[d] = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _xyz[d], _odwork );
 }
 return result;
 }
 // replicate the functionality of hex_findpt
 CartVect SpectralQuad::ievaluate( CartVect const& xyz, double /*tol*/, const CartVect& /*x0*/ ) const
 {
 // find nearest point
 realType x_star[3];
 xyz.get( x_star );
 
 realType r[2] = { 0, 0 }; // initial guess for parametric coords
 unsigned c = opt_no_constraints_3;
 realType dist = opt_findpt_2( &_data, (const realType**)_xyz, x_star, r, &c );
 // if it did not converge, get out with throw...
 if( dist > 0.9e+30 )
 {
 std::vector< CartVect > dummy;
 throw Map::EvaluationError( xyz, dummy );
 }
 
 // c tells us if we landed inside the element or exactly on a face, edge, or node
 // also, dist shows the distance to the computed point.
 // copy parametric coords back
 return CartVect( r[0], r[1], 0. );
 }
 
 Matrix3 SpectralQuad::jacobian( const CartVect& /*xi*/ ) const
 {
 // not implemented
 Matrix3 J( 0. );
 return J;
 }
 
 double SpectralQuad::evaluate_scalar_field( const CartVect& xi, const double* field ) const
 {
 // piece that we shouldn't want to cache
 int d;
 for( d = 0; d < 2; d++ )
 {
 lagrange_0( &_ld[d], xi[d] );
 }
 
 double value = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, field, _odwork );
 return value;
 }
 double SpectralQuad::integrate_scalar_field( const double* /*field_vertex_values*/ ) const
 {
 return 0.; // not implemented
 }
 // this is the same as a linear hex, although we should not need it
 bool SpectralQuad::inside_nat_space( const CartVect& xi, double& tol ) const
 {
 // just look at the box+tol here
 return ( xi[0] >= -1. - tol ) && ( xi[0] <= 1. + tol ) && ( xi[1] >= -1. - tol ) && ( xi[1] <= 1. + tol );
 }
 // something we don't do for spectral hex; we do it here because
 // we do not store the position of gl points in a tag yet
 void SpectralQuad::compute_gl_positions()
 {
 // will need to use shape functions on a simple linear quad to compute gl points
 // so we know the position of gl points in parametric space, we will just compute those
 // from the 3d vertex position (corner nodes of the quad), using simple mapping
 assert( this->vertex.size() == 4 );
 static double corner_xi[4][2] = { { -1., -1. }, { 1., -1. }, { 1., 1. }, { -1., 1. } };
 // we will use the cached lobatto nodes in parametric space _z[d] (the same in both
 // directions)
 int indexGL = 0;
 int n2 = _n * _n;
 for( int i = 0; i < _n; i++ )
 {
 double eta = _z[0][i];
 for( int j = 0; j < _n; j++ )
 {
 double csi = _z[1][j]; // we could really use the same _z[0] array of lobatto nodes
 CartVect pos( 0.0 );
 for( int k = 0; k < 4; k++ )
 {
 const double N_k = ( 1 + csi * corner_xi[k][0] ) * ( 1 + eta * corner_xi[k][1] );
 pos += N_k * vertex[k];
 }
 pos *= 0.25; // these are x, y, z of gl points; reorder them
 _glpoints[indexGL] = pos[0]; // x
 _glpoints[indexGL + n2] = pos[1]; // y
 _glpoints[indexGL + 2 * n2] = pos[2]; // z
 indexGL++;
 }
 }
 // now, we can set the _xyz pointers to internal memory allocated to these!
 _xyz[0] = &( _glpoints[0] );
 _xyz[1] = &( _glpoints[n2] );
 _xyz[2] = &( _glpoints[2 * n2] );
 }
 void SpectralQuad::get_gl_points( double*& x, double*& y, double*& z, int& psize )
 {
 x = (double*)_xyz[0];
 y = (double*)_xyz[1];
 z = (double*)_xyz[2];
 psize = _n * _n;
 }
} // namespace Element
 
} // namespace moab