Overture  Version 25
Inverse.h
Go to the documentation of this file.
1 #ifndef INVERSE_H
2 #define INVERSE_H "Inverse.h"
3 
4 #include "Mapping.h"
5 #include "BoundingBox.h" // define BoundingBox and BoundingBoxStack
6 
7 class Mapping;
8 class IntersectionMapping;
9 class GenericDataBase;
10 
11 //===========================================================================
12 // Class to define an Approximate Global Inverse
13 //
14 //===========================================================================
15 class ApproximateGlobalInverse
16 {
17  friend class Mapping;
18  friend class IntersectionMapping;
19  public:
20  static const real bogus; // Bogus value to indicate no convergence
21 
22  protected:
23 
24  Mapping *map;
25  int domainDimension;
26  int rangeDimension;
27  int uninitialized;
28  bool useRobustApproximateInverse;
29 
30  int gridDefined;
31  RealArray & grid; // holds grid points, reference to map->grid
32  IntegerArray dimension; // holds dimensions of grid
33  IntegerArray indexRange; // holds index range for edges of the grid, may depend on singularities
34 
35  void constructGrid();
36  int base,bound;
37  int base0,bound0;
38 
39  RealArray boundingBox; // grid is contained in this box
40  BoundingBox boundingBoxTree[2][3]; // root of tree of boxes for each side
41 
42  BoundingBox *serialBoundingBox; // in parallel these are bounding boxes for each local grid array.
43 
44  // BoundingBox box,box1,box2;
45  // BoundingBoxStack boxStack;
46 
47  real boundingBoxExtensionFactor; // relative amount to increase the bounding box each direction.
48  real stencilWalkBoundingBoxExtensionFactor; // stencil walk need to converge on a larger region
49  bool findBestGuess; // always find the best guess (ignore bounding boxes);
50 
51  RealArray xOrigin;
52  RealArray xTangent;
53 
54  Index Axes;
55  Index xAxes;
56 
57  // void setGrid( const RealArray & grid, const IntegerArray & gridIndexRange );
58 
59  void getPeriodicImages( const RealArray & x, RealArray & xI, int & nI,
60  const int & periodicityOfSpace, const RealArray & periodicityVector );
61 
62 
63  public: // ** for now ***
64 
65  void intersectLine( const RealArray & x, int & nI, RealArray & xI,
66  const RealArray & vector, const RealArray & xOrigin, const RealArray & xTangent );
67  void intersectPlane( const RealArray & x, int & nI, RealArray & xI,
68  const RealArray & vector, const RealArray & xOrigin, const RealArray & xTangent );
69  void intersectCube( const RealArray & x, int & nI, RealArray & xI,
70  const RealArray & vector, const RealArray & xOrigin, const RealArray & xTangent );
71 
72  void initializeBoundingBoxTrees();
73  void binarySearchOverBoundary( real x[3], real & minimumDistance, int iv[3], int side=-1, int axis=-1 );
74 /* ---
75  void robustBinarySearchOverBoundary( real x[3],
76  real & minimumDistance,
77  int iv[3],
78  int side,
79  int axis );
80 --- */
81  int insideGrid( int side, int axis, real x[], int iv[], real & dot );
82 
83  int distanceToCell( real x[], int iv[], real & signedDistance, const real minimumDistance );
84 
85  void findNearestGridPoint( const int base, const int bound, RealArray & x, RealArray & r );
86  int findNearestCell(real x[3], int iv[3], real & minimumDistance );
87 
88  void initializeStencilWalk();
89  void countCrossingsWithPolygon(const RealArray & x,
90  IntegerArray & crossings,
91  const int & side=Start,
92  const int & axis=axis1,
93  RealArray & xCross = Overture::nullRealArray(),
94  const IntegerArray & mask = Overture::nullIntArray(),
95  const unsigned int & maskBit = UINT_MAX, // (UINT_MAX : all bits on)
96  const int & maskRatio1 =1 ,
97  const int & maskRatio2 =1 ,
98  const int & maskRatio3 =1 );
99 
100  public:
101  ApproximateGlobalInverse( Mapping & map );
102  virtual ~ApproximateGlobalInverse();
103  virtual void inverse( const RealArray & x, RealArray & r, RealArray & rx,
104  MappingWorkSpace & workSpace, MappingParameters & params );
105 
106  void initialize(); // initialize if not already done so
107  void reinitialize(); // this will force a re-initialize the inverse
108  const RealArray & getGrid() const; // return the grid used for the inverse
109  const RealArray & getBoundingBox() const;
110  const BoundingBox & getBoundingBoxTree(int side, int axis) const;
111 
112  real getParameter( const MappingParameters::realParameter & param ) const;
113  int getParameter( const MappingParameters::intParameter & param ) const;
114  virtual void setParameter( const MappingParameters::realParameter & param, const real & value );
115  virtual void setParameter( const MappingParameters::intParameter & param, const int & value );
116  virtual void useRobustInverse(const bool trueOrFalse=TRUE );
117  virtual bool usingRobustInverse() const;
118 
119  // return size of this object
120  virtual real sizeOf(FILE *file = NULL ) const;
121 
122  virtual int get( const GenericDataBase & dir, const aString & name); // get from a database file
123  virtual int put( GenericDataBase & dir, const aString & name) const; // put to a database file
124 
125  public:
126  // these are used for timings and statistics
127  static real timeForApproximateInverse,
128  timeForFindNearestGridPoint,
129  timeForBinarySearchOverBoundary,
130  timeForBinarySearchOnLeaves;
131 
132  static int numberOfStencilWalks,
133  numberOfStencilSearches,
134  numberOfBinarySearches,
135  numberOfBoxesChecked,
136  numberOfBoundingBoxes;
137 
138  static void printStatistics();
139 
140  // these next arrays are for the stencil walk
141  static int numberOfStencilDir2D[9]; // (3,3);
142  static int stencilDir2D1[8*3*3]; // (8,3,3)
143  static int stencilDir2D2[8*3*3];
144  static int stencilDir2D3[8*3*3];
145  static int numberOfStencilDir3D[27]; // (3,3,3);
146  static int stencilDir3D1[27*3*3*3]; // (27,3,3,3)
147  static int stencilDir3D2[27*3*3*3];
148  static int stencilDir3D3[27*3*3*3];
149 
150 };
151 
152 
153 //===========================================================================
154 // Class to define an Exact Local Inverse
155 //
156 //===========================================================================
157 const int maximumNumberOfRecursionLevels=5;
158 
159 class ExactLocalInverse
160 {
161  private:
162  Mapping *map;
163  int domainDimension;
164  int rangeDimension;
165  RealArray periodVector;
166  int base,bound;
167  Index Axes;
168  Index xAxes;
169  int uninitialized;
170  bool useRobustExactLocalInverse;
171 
172  real nonConvergenceValue; // value given to inverse when there is no convergence
173  real newtonToleranceFactor; // convergence tolerance is this times the machine epsilon
174  real newtonDivergenceValue; // newton is deemed to have diverged if the r value is this much outside [0,1]
175  real newtonL2Factor; // extra factor used in inverting the closest point to a curve or surface
176 
177  // Work Arrays for Newton:
178  // RealArray y,yr,r2,yr2; // these arrays cannot be static if the inverse is called recursively
179  // bool workArraysTooSmall;
180 
181  void initialize();
182  inline void periodicShift( RealArray & r, const Index & I );
183  void underdeterminedLS(const RealArray & xt,
184  const RealArray & tx, // *** should not be const -- do for IBM compiler
185  const RealArray & dy,
186  const RealArray & dr, // *** should not be const -- do for IBM compiler
187  real & det );
188 
189  inline void invert(RealArray & yr, RealArray & dy, RealArray & det,
190  RealArray & ry, RealArray & dr, IntegerArray & status);
191  inline void invertL2(RealArray & yr, RealArray & dy, RealArray & det, RealArray & yr2, RealArray & yrr,
192  RealArray & ry, RealArray & dr, IntegerArray & status );
193  void minimizeByBisection(RealArray & r, RealArray & x, RealArray & dr, IntegerArray & status, real & eps );
194 
195  public:
196  ExactLocalInverse( Mapping & map );
197  virtual ~ExactLocalInverse();
198  virtual void inverse( const RealArray & x1, RealArray & r1, RealArray & rx1,
199  MappingWorkSpace & workSpace, const int computeGlobalInverse=FALSE );
200 
201  void reinitialize();
202 
203  real getParameter( const MappingParameters::realParameter & param ) const;
204  virtual void setParameter( const MappingParameters::realParameter & param, const real & value );
205 
206  virtual void useRobustInverse(const bool trueOrFalse=true );
207 
208  // return size of this object
209  virtual real sizeOf(FILE *file = NULL ) const;
210 
211  virtual int get( const GenericDataBase & dir, const aString & name); // get from a database file
212  virtual int put( GenericDataBase & dir, const aString & name) const; // put to a database file
213 
214  public:
215  static real timeForExactInverse;
216  static int numberOfNewtonInversions,
217  numberOfNewtonSteps;
218 
219  protected:
220  bool mappingHasACoordinateSingularity;
221 
222  int compressConvergedPoints(Index & I,
223  RealArray & x,
224  RealArray & r,
225  RealArray & ry,
226  RealArray & det,
227  IntegerArray & status,
228  const RealArray & x1,
229  RealArray & r1,
230  RealArray & rx1,
231  MappingWorkSpace & workSpace,
232  const int computeGlobalInverse );
233 
234 };
235 
236 
237 
238 
239 #endif // "Inverse.h"
Generated on Fri Jan 4 2013 10:17:54 for Overture by   doxygen 1.8.3