Chombo + EB + MF  3.2
PolyGeom.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
8  */
9 #endif
10
11 // ANAG, LBNL
12
13 #ifndef _POLYGEOM_H_
14 #define _POLYGEOM_H_
15
16 #include "REAL.H"
17 #include "RealVect.H"
18 #include "Tuple.H"
19 #include "Vector.H"
20
21 #include "VolIndex.H"
22
24
25 class EBISBox;
26 ///
27 /**
28  PolyGeom encapsulates the functions
29  to generate centroids, normals, etc given
30  the geometric information.
31  It is meant to be used simply as a static class.
32  Its only data member (a static) is its tolerance.
33 */
34 class PolyGeom
35 {
36 public:
37
38  ///
39  /**
40  return the inverse of the input matrix. done using cramers rule
41  */
42  static void invertMatrix(Real a_AInverse[SpaceDim][SpaceDim],
43  const Real a_A[SpaceDim][SpaceDim], bool a_test);
44
45  ///
46  /**
47  // calculate the cofactor of element (row,col)
48  //this is the matrix which excludes row and column row, col
49  */
50  static void getMinor(Real a_Aminor[SpaceDim-1][SpaceDim-1],
51  const Real a_A[SpaceDim ][SpaceDim ],
52  int a_row, int acol);
53
54  /****/
55  static void
56  getJacobianAndInverse(Real a_Jacobian[SpaceDim][SpaceDim],
57  Real a_Jinverse[SpaceDim][SpaceDim],
58  RealVect& a_normal,
59  RealVect a_tangents[SpaceDim-1])
60  {
61  //make the first row = the normal
62  //make next row the tangent vectors
63  for (int icol = 0; icol < SpaceDim; icol++)
64  {
65  a_Jacobian[0][icol] = a_normal[icol];
66  for (int irow = 1; irow < SpaceDim; irow++)
67  {
68  a_Jacobian[irow][icol] = a_tangents[irow-1][icol];
69  }
70  }
71  //compute the inverse of the matrix. the true is turn on
72  //the test function.
73  PolyGeom::invertMatrix(a_Jinverse, a_Jacobian, true);
74  }
75
76 //given a normal vector, get spacedim-1 tangential vectors that are
77 //unit vectors and
78  static void
80  const RealVect& a_normal)
81  {
82 #if CH_SPACEDIM==2
83  a_tangen[0][1] = a_normal[0];
84  a_tangen[0][0] = -a_normal[1];
85 #else
86  //get some vector that is not parallel to a_normal
87  RealVect notNormal = a_normal + RealVect::Unit;
88  Real dotProd = PolyGeom::dot(notNormal, a_normal);
89  //there is one bizzare case where I have just scaled up a_normal
90  //in this case, add something diff to each component.
91  if (Abs(dotProd) < 1.0e-6)
92  {
93  for (int idir = 0; idir < SpaceDim; idir++)
94  {
95  notNormal[idir] += Real(idir+1);
96  }
97  dotProd = PolyGeom::dot(notNormal, a_normal);
98  CH_assert(Abs(dotProd) > 1.0e-6);
99  }
100  Real sumSquare;
101  PolyGeom::unifyVector(notNormal, sumSquare);
102  //make the first tangential vector = cross(notnormal, normal);
103  a_tangen[0] = PolyGeom::cross(notNormal, a_normal);
104  //make the second one the cross between that and normal
105  a_tangen[1] = PolyGeom::cross(a_tangen[0], a_normal);
106 #endif
107
108  //make them all unit vectors
109  for (int itan = 0; itan < CH_SPACEDIM-1; itan++)
110  {
111  Real sumSquare;
112  PolyGeom::unifyVector(a_tangen[itan], sumSquare);
113  }
114  }
115
116  ///
117  static Real determinantSD(const Real a_A[SpaceDim][SpaceDim]);
118
119  ///
120  static Real determinantSDM1(const Real a_A[SpaceDim-1][SpaceDim-1]);
121
122  ///
123  /**
124  find equation of line
125  between a line formed by a_pointOnLine = p0 and a_direction=v
126  l(t) = p0 + v t
127  and a point in space a_point.
128  returns closest point on line and the normal between the
129  closest point and the point in space
130  */
131  static void
132  pointToLine(RealVect& a_closestPt,
133  RealVect& a_normal,
134  const RealVect& a_point,
135  const RealVect& a_pointOnLine,
136  const RealVect& a_direction);
137
138  ///return distance from point to a plane
139  static Real
140  distancePointToPlane(const RealVect& a_point,
141  const RealVect& a_normal,
142  const RealVect& a_pointOnLine);
143
144  ///
145  /**
146  */
147  static void setTolerance(const Real& a_tolerance);
148
149  ///
150  /**
151  */
152  static void setVectDx(const RealVect& a_vectDx);
153
154  ///
155  /**
156  */
157  static void setVolumeTolerance(const Real& a_tolerance);
158
159  ///
160  /**
161  */
162  static void setAreaTolerance(const Real& a_tolerance);
163
164  ///
165  /**
166  */
167  static void setLengthTolerance(const Real& a_tolerance);
168
169  ///
170  /**
171  compute the cross product between xvec0 and xvec1
172  (returns xvec1 x xvec0)
173  */
174  static RealVect cross(const RealVect& a_xvec1,
175  const RealVect& a_xvec0);
176
177  ///
178  /**
179  compute the dot product between xvec0 and xvec1
180  (returns xvec1 dot xvec0)
181  */
182  static Real dot(const RealVect& a_xvec1,
183  const RealVect& a_xvec0);
184
185  ///
186  /**
187  */
188  static const Real& getTolerance();
189
190  ///
191  /**
192  */
193  static const RealVect& getVectDx();
194
195  ///
196  /**
197  */
198  static const Real& getVolumeTolerance();
199
200  ///
201  /**
202  */
203  static const Real& getAreaTolerance();
204
205  ///
206  /**
207  */
208  static const Real& getLengthTolerance();
209
210  ///
211  /**
212  Return the normal to the boundary at the
213  input VoF. If said normal is undefined,
214  returns the zero vector.
215  */
216  static RealVect normal(const VolIndex& a_vof,
217  const EBISBox& a_ebisBox,
218  const Real& a_bndryArea);
219
220  ///
221  /**
222  */
223  static Real bndryArea(const VolIndex& a_vof,
224  const EBISBox& a_ebisBox);
225
226  ///
227  /**
228  Return the inhomogeneous component of the plane
229  eqution ni xi = alpha.
230  Volume fraction must be positive and <= 1.0
231  and the normals can be in any order.
232  */
233  static Real computeAlpha(const Real& a_volFrac,
234  const RealVect& a_normal);
235
236  ///
237  /**
238  Compute the volume fraction in the cell given
239  alpha and the normal.
240  */
241  static Real computeVolume(const Real& a_alpha,
242  const RealVect& a_normal);
243
244  ///
245  /**
246  Get the normal and alpha given the points in a polygon.
247  Solves the linear system \\
248  ni xi - alpha = 0\\
249  \sum ni*ni = 1\\
250  for ni and alpha.
251  Always returns nromal[updir] >= 0
252  */
253  static void computeNormalAndAlpha(Real& a_alpha,
254  RealVect& a_normal,
255  const int& a_upDir,
256  const Tuple<RealVect, CH_SPACEDIM>& a_poly);
257
258  // solve ax = b (gives icomp component of x using cramers rule).
259  static Real
260  matrixSolveComp(const Vector<Vector<Real> >& a_A,
261  const Vector<Real>& a_rhs,
262  const int& a_icomp);
263
264  //return the determinant of a
265  static Real
266  determinant(const Vector< Vector<Real> >& a_A);
267
268  /**
269  andzvolume computes the volume under
270  the curve (covered volume).
271  Does so using Scardovelli
272  and Zaleski's direct formulae.
273  The normal has to be all positive numbers.
274  The normals also must be ordered from lowest
275  to highest.
276  */
277  static Real sAndZVolume(const Real& a_alpha,
278  const RealVect& a_normal);
279
280  ///
282  {;}
283  ///
285  {;}
286
287  ///
288  /**
289  Sort vector from low to high. All comps must be positive.
290  Return reverse mapping
291  */
292  static void sortVector(RealVect& vect, IntVect& ivmap);
293
294  ///make vector all pos and return the signs
295  static void posifyVector(RealVect& vect, IntVect& signs);
296
297  ///make vector into a unit vector and return the sum of squares
298  static void unifyVector(RealVect& normal, Real& sumSquare);
299
300  ///
301  static Tuple<int, CH_SPACEDIM-1> computeTanDirs(int upDir);
302
303 protected:
304
305  static Real twoDFunc(const Real& arg);
306  static Real threeDFunc(const Real& arg);
307
308  static Real computeAnyVolume(const Real& a_alpha,
309  const Real& a_norm0,
310  const Real& a_norm1,
311  const Real& a_norm2);
312
314
316
318
320
321  static RealVect tetCentroid(const RealVect& normal,
322  const Real& alpha);
323
324  static Real tetVolume(const RealVect& normal,
325  const Real& alpha);
326
327
328 private:
329 };
330
331 #include "NamespaceFooter.H"
332 #endif
static Real determinant(const Vector< Vector< Real > > &a_A)
static Real matrixSolveComp(const Vector< Vector< Real > > &a_A, const Vector< Real > &a_rhs, const int &a_icomp)
static void getMinor(Real a_Aminor[SpaceDim-1][SpaceDim-1], const Real a_A[SpaceDim][SpaceDim], int a_row, int acol)
#define CH_SPACEDIM
Definition: SPACE.H:51
#define CH_assert(cond)
Definition: CHArray.H:37
static RealVect cross(const RealVect &a_xvec1, const RealVect &a_xvec0)
static Real bndryArea(const VolIndex &a_vof, const EBISBox &a_ebisBox)
static void setVectDx(const RealVect &a_vectDx)
one dimensional dynamic array
Definition: Vector.H:53
static Real s_areatolerance
Definition: PolyGeom.H:319
static const RealVect & getVectDx()
static RealVect tetCentroid(const RealVect &normal, const Real &alpha)
Definition: EBISBox.H:46
PolyGeom()
Definition: PolyGeom.H:281
static void getJacobianAndInverse(Real a_Jacobian[SpaceDim][SpaceDim], Real a_Jinverse[SpaceDim][SpaceDim], RealVect &a_normal, RealVect a_tangents[SpaceDim-1])
Definition: PolyGeom.H:56
static void getTangentVectors(RealVect a_tangen[CH_SPACEDIM-1], const RealVect &a_normal)
Definition: PolyGeom.H:79
static void invertMatrix(Real a_AInverse[SpaceDim][SpaceDim], const Real a_A[SpaceDim][SpaceDim], bool a_test)
static RealVect normal(const VolIndex &a_vof, const EBISBox &a_ebisBox, const Real &a_bndryArea)
~PolyGeom()
Definition: PolyGeom.H:284
static Real s_lengthtolerance
Definition: PolyGeom.H:317
static Real computeVolume(const Real &a_alpha, const RealVect &a_normal)
const int SpaceDim
Definition: SPACE.H:38
static const RealVect Unit
Definition: RealVect.H:427
static const Real & getTolerance()
static void setLengthTolerance(const Real &a_tolerance)
Definition: PolyGeom.H:34
static Real s_tolerance
Definition: PolyGeom.H:315
static void sortVector(RealVect &vect, IntVect &ivmap)
static Real determinantSD(const Real a_A[SpaceDim][SpaceDim])
static Real threeDFunc(const Real &arg)
double Real
Definition: REAL.H:33
static Tuple< int, CH_SPACEDIM-1 > computeTanDirs(int upDir)
static Real computeAnyVolume(const Real &a_alpha, const Real &a_norm0, const Real &a_norm1, const Real &a_norm2)
static Real dot(const RealVect &a_xvec1, const RealVect &a_xvec0)
Ordered Tuples for Types T.
Definition: Tuple.H:30
T Abs(const T &a_a)
Definition: Misc.H:53
static const Real & getVolumeTolerance()
static void setAreaTolerance(const Real &a_tolerance)
A Real vector in SpaceDim-dimensional space.
Definition: RealVect.H:41
static const Real & getLengthTolerance()
static void computeNormalAndAlpha(Real &a_alpha, RealVect &a_normal, const int &a_upDir, const Tuple< RealVect, CH_SPACEDIM > &a_poly)
static Real twoDFunc(const Real &arg)
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
Volume of Fluid Index.
Definition: VolIndex.H:31
static void setTolerance(const Real &a_tolerance)
static Real distancePointToPlane(const RealVect &a_point, const RealVect &a_normal, const RealVect &a_pointOnLine)
return distance from point to a plane
static RealVect s_vectDx
Definition: PolyGeom.H:313
static void posifyVector(RealVect &vect, IntVect &signs)
make vector all pos and return the signs
static Real determinantSDM1(const Real a_A[SpaceDim-1][SpaceDim-1])
static const Real & getAreaTolerance()
static Real tetVolume(const RealVect &normal, const Real &alpha)
static void pointToLine(RealVect &a_closestPt, RealVect &a_normal, const RealVect &a_point, const RealVect &a_pointOnLine, const RealVect &a_direction)
static void unifyVector(RealVect &normal, Real &sumSquare)
make vector into a unit vector and return the sum of squares
static Real computeAlpha(const Real &a_volFrac, const RealVect &a_normal)
static Real sAndZVolume(const Real &a_alpha, const RealVect &a_normal)
static void setVolumeTolerance(const Real &a_tolerance)