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