MOAB: Mesh Oriented datABase  (version 5.5.0)
geomObject.cpp
Go to the documentation of this file.
1 #include <cmath>
2 #include <vector>
3 #include <stdexcept>
4 #include <cassert>
5 
6 class geomObject
7 {
8  public:
9  virtual ~geomObject() {}
10  virtual void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const = 0;
11  virtual double compute_projecterror( int dim, double* oldcoords ) const
12  {
13  double* newcoords = new double[dim]();
14  project_points2geom( dim, oldcoords, newcoords, NULL );
15  double ans = 0;
16 
17  for( int i = 0; i < dim; ++i )
18  {
19  ans += ( newcoords[i] - oldcoords[i] ) * ( newcoords[i] - oldcoords[i] );
20  }
21 
22  delete[] newcoords;
23  return sqrt( ans );
24  }
25  virtual void compute_projecterror( int dim,
26  int nverts,
27  double* oldcoords,
28  double& l1err,
29  double& l2err,
30  double& linferr ) const
31  {
32  l1err = l2err = linferr = 0;
33 
34  for( int i = 0; i < nverts; ++i )
35  {
36  double err = compute_projecterror( dim, oldcoords + i * dim );
37  l1err += err;
38  l2err += err * err;
39  linferr = std::max( linferr, err );
40  }
41 
42  l1err /= nverts;
43  l2err = sqrt( l2err / nverts );
44  }
45  inline double Twonorm( int dim, double* vec ) const
46  {
47  double ans = 0;
48 
49  for( int i = 0; i < dim; ++i )
50  {
51  ans += vec[i] * vec[i];
52  }
53 
54  return sqrt( ans );
55  }
56 };
57 
58 class sphere : public geomObject
59 {
60  public:
61  sphere( double x = 0, double y = 0, double z = 0, double r = 1 )
62  : centerx( x ), centery( y ), centerz( z ), radius( r )
63  {
64  assert( r > 0 );
65  }
66  virtual ~sphere() {}
67  void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const
68  {
69  if( oldcoords == NULL || newcoords == NULL )
70  {
71  throw std::invalid_argument( "NULL pointer" );
72  }
73 
74  std::vector< double > direction;
75  direction.push_back( oldcoords[0] - centerx );
76  direction.push_back( oldcoords[1] - centery );
77 
78  if( dim == 3 )
79  {
80  direction.push_back( oldcoords[2] - centerz );
81  }
82 
83  double len = geomObject::Twonorm( dim, &( direction[0] ) );
84  assert( len > 0 );
85 
86  for( int i = 0; i < dim; ++i )
87  {
88  direction[i] /= len;
89  }
90 
91  newcoords[0] = centerx + direction[0] * radius;
92  newcoords[1] = centery + direction[1] * radius;
93 
94  if( dim == 2 )
95  {
96  if( derivs )
97  {
98  derivs[0] = -direction[1];
99  derivs[1] = direction[0];
100  }
101 
102  return;
103  }
104  else if( dim == 3 )
105  {
106  newcoords[2] = centerz + direction[2] * radius;
107 
108  if( derivs )
109  {
110  derivs[0] = direction[0];
111  derivs[1] = direction[1];
112  derivs[2] = direction[2];
113  }
114  }
115  else
116  {
117  throw std::invalid_argument( "dim must be 2 or 3" );
118  }
119  }
120 
121  private:
122  double centerx, centery, centerz, radius;
123 };
124 
125 class torus : public geomObject
126 {
127  public:
128  torus( double x = 0, double y = 0, double z = 0, double aa = 0.3, double cc = 1.0 )
129  : centerx( x ), centery( y ), centerz( z ), a( aa ), c( cc )
130  {
131  assert( aa > 0 && cc > 0 );
132  }
133  virtual ~torus() {}
134  void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const
135  {
136  assert( dim == 3 && oldcoords && newcoords );
137  double transfer[3] = { oldcoords[0] - centerx, oldcoords[1] - centery, oldcoords[2] - centerz };
138  double twodnrm = geomObject::Twonorm( 2, transfer );
139  double tubecenter[3] = { c * transfer[0] / twodnrm + centerx, c * transfer[1] / twodnrm + centery, centerz };
140  double direction[3] = { oldcoords[0] - tubecenter[0], oldcoords[1] - tubecenter[1],
141  oldcoords[2] - tubecenter[2] };
142  double len = geomObject::Twonorm( 3, direction );
143  assert( len > 0 );
144  direction[0] /= len;
145  direction[1] /= len;
146  direction[2] /= len;
147 
148  if( derivs )
149  {
150  derivs[0] = direction[0];
151  derivs[1] = direction[1];
152  derivs[2] = direction[2];
153  }
154 
155  newcoords[0] = tubecenter[0] + a * direction[0];
156  newcoords[1] = tubecenter[1] + a * direction[1];
157  newcoords[2] = tubecenter[2] + a * direction[2];
158  }
159 
160  private:
161  double centerx, centery, centerz, a, c;
162 };