MOAB: Mesh Oriented datABase  (version 5.4.1)
geomObject.cpp
Go to the documentation of this file.
00001 #include <cmath>
00002 #include <vector>
00003 #include <stdexcept>
00004 #include <cassert>
00005 
00006 class geomObject
00007 {
00008   public:
00009     virtual ~geomObject() {}
00010     virtual void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const = 0;
00011     virtual double compute_projecterror( int dim, double* oldcoords ) const
00012     {
00013         double* newcoords = new double[dim]();
00014         project_points2geom( dim, oldcoords, newcoords, NULL );
00015         double ans = 0;
00016 
00017         for( int i = 0; i < dim; ++i )
00018         {
00019             ans += ( newcoords[i] - oldcoords[i] ) * ( newcoords[i] - oldcoords[i] );
00020         }
00021 
00022         delete[] newcoords;
00023         return sqrt( ans );
00024     }
00025     virtual void compute_projecterror( int dim,
00026                                        int nverts,
00027                                        double* oldcoords,
00028                                        double& l1err,
00029                                        double& l2err,
00030                                        double& linferr ) const
00031     {
00032         l1err = l2err = linferr = 0;
00033 
00034         for( int i = 0; i < nverts; ++i )
00035         {
00036             double err = compute_projecterror( dim, oldcoords + i * dim );
00037             l1err += err;
00038             l2err += err * err;
00039             linferr = std::max( linferr, err );
00040         }
00041 
00042         l1err /= nverts;
00043         l2err = sqrt( l2err / nverts );
00044     }
00045     inline double Twonorm( int dim, double* vec ) const
00046     {
00047         double ans = 0;
00048 
00049         for( int i = 0; i < dim; ++i )
00050         {
00051             ans += vec[i] * vec[i];
00052         }
00053 
00054         return sqrt( ans );
00055     }
00056 };
00057 
00058 class sphere : public geomObject
00059 {
00060   public:
00061     sphere( double x = 0, double y = 0, double z = 0, double r = 1 )
00062         : centerx( x ), centery( y ), centerz( z ), radius( r )
00063     {
00064         assert( r > 0 );
00065     }
00066     virtual ~sphere() {}
00067     void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const
00068     {
00069         if( oldcoords == NULL || newcoords == NULL )
00070         {
00071             throw std::invalid_argument( "NULL pointer" );
00072         }
00073 
00074         std::vector< double > direction;
00075         direction.push_back( oldcoords[0] - centerx );
00076         direction.push_back( oldcoords[1] - centery );
00077 
00078         if( dim == 3 )
00079         {
00080             direction.push_back( oldcoords[2] - centerz );
00081         }
00082 
00083         double len = geomObject::Twonorm( dim, &( direction[0] ) );
00084         assert( len > 0 );
00085 
00086         for( int i = 0; i < dim; ++i )
00087         {
00088             direction[i] /= len;
00089         }
00090 
00091         newcoords[0] = centerx + direction[0] * radius;
00092         newcoords[1] = centery + direction[1] * radius;
00093 
00094         if( dim == 2 )
00095         {
00096             if( derivs )
00097             {
00098                 derivs[0] = -direction[1];
00099                 derivs[1] = direction[0];
00100             }
00101 
00102             return;
00103         }
00104         else if( dim == 3 )
00105         {
00106             newcoords[2] = centerz + direction[2] * radius;
00107 
00108             if( derivs )
00109             {
00110                 derivs[0] = direction[0];
00111                 derivs[1] = direction[1];
00112                 derivs[2] = direction[2];
00113             }
00114         }
00115         else
00116         {
00117             throw std::invalid_argument( "dim must be 2 or 3" );
00118         }
00119     }
00120 
00121   private:
00122     double centerx, centery, centerz, radius;
00123 };
00124 
00125 class torus : public geomObject
00126 {
00127   public:
00128     torus( double x = 0, double y = 0, double z = 0, double aa = 0.3, double cc = 1.0 )
00129         : centerx( x ), centery( y ), centerz( z ), a( aa ), c( cc )
00130     {
00131         assert( aa > 0 && cc > 0 );
00132     }
00133     virtual ~torus() {}
00134     void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const
00135     {
00136         assert( dim == 3 && oldcoords && newcoords );
00137         double transfer[3]   = { oldcoords[0] - centerx, oldcoords[1] - centery, oldcoords[2] - centerz };
00138         double twodnrm       = geomObject::Twonorm( 2, transfer );
00139         double tubecenter[3] = { c * transfer[0] / twodnrm + centerx, c * transfer[1] / twodnrm + centery, centerz };
00140         double direction[3]  = { oldcoords[0] - tubecenter[0], oldcoords[1] - tubecenter[1],
00141                                 oldcoords[2] - tubecenter[2] };
00142         double len           = geomObject::Twonorm( 3, direction );
00143         assert( len > 0 );
00144         direction[0] /= len;
00145         direction[1] /= len;
00146         direction[2] /= len;
00147 
00148         if( derivs )
00149         {
00150             derivs[0] = direction[0];
00151             derivs[1] = direction[1];
00152             derivs[2] = direction[2];
00153         }
00154 
00155         newcoords[0] = tubecenter[0] + a * direction[0];
00156         newcoords[1] = tubecenter[1] + a * direction[1];
00157         newcoords[2] = tubecenter[2] + a * direction[2];
00158     }
00159 
00160   private:
00161     double centerx, centery, centerz, a, c;
00162 };
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines