00001 #ifdef CH_LANG_CC
00002
00003
00004
00005
00006
00007
00008
00009 #endif
00010
00011
00012
00013 #ifndef _POLYGEOM_H_
00014 #define _POLYGEOM_H_
00015
00016 #include "REAL.H"
00017 #include "RealVect.H"
00018 #include "Tuple.H"
00019 #include "Vector.H"
00020
00021 #include "VolIndex.H"
00022
00023 #include "NamespaceHeader.H"
00024
00025 class EBISBox;
00026
00027
00028
00029
00030
00031
00032
00033
00034 class PolyGeom
00035 {
00036 public:
00037
00038
00039
00040
00041
00042 static void invertMatrix(Real a_AInverse[SpaceDim][SpaceDim],
00043 const Real a_A[SpaceDim][SpaceDim], bool a_test);
00044
00045
00046
00047
00048
00049
00050 static void getMinor(Real a_Aminor[SpaceDim-1][SpaceDim-1],
00051 const Real a_A[SpaceDim ][SpaceDim ],
00052 int a_row, int acol);
00053
00054
00055 static void
00056 getJacobianAndInverse(Real a_Jacobian[SpaceDim][SpaceDim],
00057 Real a_Jinverse[SpaceDim][SpaceDim],
00058 RealVect& a_normal,
00059 RealVect a_tangents[SpaceDim-1])
00060 {
00061
00062
00063 for (int icol = 0; icol < SpaceDim; icol++)
00064 {
00065 a_Jacobian[0][icol] = a_normal[icol];
00066 for (int irow = 1; irow < SpaceDim; irow++)
00067 {
00068 a_Jacobian[irow][icol] = a_tangents[irow-1][icol];
00069 }
00070 }
00071
00072
00073 PolyGeom::invertMatrix(a_Jinverse, a_Jacobian, true);
00074 }
00075
00076
00077
00078 static void
00079 getTangentVectors( RealVect a_tangen[CH_SPACEDIM-1],
00080 const RealVect& a_normal)
00081 {
00082 #if CH_SPACEDIM==2
00083 a_tangen[0][1] = a_normal[0];
00084 a_tangen[0][0] = -a_normal[1];
00085 #else
00086
00087 RealVect notNormal = a_normal + RealVect::Unit;
00088 Real dotProd = PolyGeom::dot(notNormal, a_normal);
00089
00090
00091 if (Abs(dotProd) < 1.0e-6)
00092 {
00093 for (int idir = 0; idir < SpaceDim; idir++)
00094 {
00095 notNormal[idir] += Real(idir+1);
00096 }
00097 dotProd = PolyGeom::dot(notNormal, a_normal);
00098 CH_assert(Abs(dotProd) > 1.0e-6);
00099 }
00100 Real sumSquare;
00101 PolyGeom::unifyVector(notNormal, sumSquare);
00102
00103 a_tangen[0] = PolyGeom::cross(notNormal, a_normal);
00104
00105 a_tangen[1] = PolyGeom::cross(a_tangen[0], a_normal);
00106 #endif
00107
00108
00109 for (int itan = 0; itan < CH_SPACEDIM-1; itan++)
00110 {
00111 Real sumSquare;
00112 PolyGeom::unifyVector(a_tangen[itan], sumSquare);
00113 }
00114 }
00115
00116
00117 static Real determinantSD(const Real a_A[SpaceDim][SpaceDim]);
00118
00119
00120 static Real determinantSDM1(const Real a_A[SpaceDim-1][SpaceDim-1]);
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 static void
00132 pointToLine(RealVect& a_closestPt,
00133 RealVect& a_normal,
00134 const RealVect& a_point,
00135 const RealVect& a_pointOnLine,
00136 const RealVect& a_direction);
00137
00138
00139 static Real
00140 distancePointToPlane(const RealVect& a_point,
00141 const RealVect& a_normal,
00142 const RealVect& a_pointOnLine);
00143
00144
00145
00146
00147 static void setTolerance(const Real& a_tolerance);
00148
00149
00150
00151
00152 static void setVectDx(const RealVect& a_vectDx);
00153
00154
00155
00156
00157 static void setVolumeTolerance(const Real& a_tolerance);
00158
00159
00160
00161
00162 static void setAreaTolerance(const Real& a_tolerance);
00163
00164
00165
00166
00167 static void setLengthTolerance(const Real& a_tolerance);
00168
00169
00170
00171
00172
00173
00174 static RealVect cross(const RealVect& a_xvec1,
00175 const RealVect& a_xvec0);
00176
00177
00178
00179
00180
00181
00182 static Real dot(const RealVect& a_xvec1,
00183 const RealVect& a_xvec0);
00184
00185
00186
00187
00188 static const Real& getTolerance();
00189
00190
00191
00192
00193 static const RealVect& getVectDx();
00194
00195
00196
00197
00198 static const Real& getVolumeTolerance();
00199
00200
00201
00202
00203 static const Real& getAreaTolerance();
00204
00205
00206
00207
00208 static const Real& getLengthTolerance();
00209
00210
00211
00212
00213
00214
00215
00216 static RealVect normal(const VolIndex& a_vof,
00217 const EBISBox& a_ebisBox,
00218 const Real& a_bndryArea);
00219
00220
00221
00222
00223 static Real bndryArea(const VolIndex& a_vof,
00224 const EBISBox& a_ebisBox);
00225
00226
00227
00228
00229
00230
00231
00232
00233 static Real computeAlpha(const Real& a_volFrac,
00234 const RealVect& a_normal);
00235
00236
00237
00238
00239
00240
00241 static Real computeVolume(const Real& a_alpha,
00242 const RealVect& a_normal);
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253 static void computeNormalAndAlpha(Real& a_alpha,
00254 RealVect& a_normal,
00255 const int& a_upDir,
00256 const Tuple<RealVect, CH_SPACEDIM>& a_poly);
00257
00258
00259 static Real
00260 matrixSolveComp(const Vector<Vector<Real> >& a_A,
00261 const Vector<Real>& a_rhs,
00262 const int& a_icomp);
00263
00264
00265 static Real
00266 determinant(const Vector< Vector<Real> >& a_A);
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 static Real sAndZVolume(const Real& a_alpha,
00278 const RealVect& a_normal);
00279
00280
00281 PolyGeom()
00282 {;}
00283
00284 ~PolyGeom()
00285 {;}
00286
00287
00288
00289
00290
00291
00292 static void sortVector(RealVect& vect, IntVect& ivmap);
00293
00294
00295 static void posifyVector(RealVect& vect, IntVect& signs);
00296
00297
00298 static void unifyVector(RealVect& normal, Real& sumSquare);
00299
00300
00301 static Tuple<int, CH_SPACEDIM-1> computeTanDirs(int upDir);
00302
00303 protected:
00304
00305 static Real twoDFunc(const Real& arg);
00306 static Real threeDFunc(const Real& arg);
00307
00308 static Real computeAnyVolume(const Real& a_alpha,
00309 const Real& a_norm0,
00310 const Real& a_norm1,
00311 const Real& a_norm2);
00312
00313 static RealVect s_vectDx;
00314
00315 static Real s_tolerance;
00316
00317 static Real s_lengthtolerance;
00318
00319 static Real s_areatolerance;
00320
00321 static RealVect tetCentroid(const RealVect& normal,
00322 const Real& alpha);
00323
00324 static Real tetVolume(const RealVect& normal,
00325 const Real& alpha);
00326
00327
00328 private:
00329 };
00330
00331 #include "NamespaceFooter.H"
00332 #endif