00001 #ifndef _MY_LightOrthoMesh_
00002 #define _MY_LightOrthoMesh_
00003 #include "globals.h"
00004 #include <assert.h>
00005 #include "sfunctions.h"
00006 #include "wellinfo.h"
00007 #include <vector>
00008 #include "delindexlst.h"
00009 #include "opointer.h"
00010
00011 #ifndef _GLOBALS__
00012
00013 class Point3D
00014 {
00015 private:
00016 double data[3];
00017 public:
00018 Point3D(double d1,double d2,double d3)
00019 {
00020 data[0]=d1;
00021 data[1]=d2;
00022 data[2]=d3;
00023 }
00024 Point3D(){}
00025
00026 ~Point3D(){}
00027 double& operator[](int index)
00028 {
00029 assert(index < 3);
00030 return data[index];
00031 }
00032 double operator[](int index) const
00033 {
00034 assert(index < 3);
00035 return data[index];
00036 }
00037
00038
00039 };
00040
00041 enum CellDirection3D {LEFT_CELL,RIGHT_CELL,FRONT_CELL,BACK_CELL,BOTTOM_CELL,UP_CELL};
00042 enum VertexDirection3D{VERTEX_000,VERTEX_100,VERTEX_010,VERTEX_110,VERTEX_001,VERTEX_101,VERTEX_011,VERTEX_111};
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 enum FaceDirection3D{LEFT_FACE=0,RIGHT_FACE=1,FRONT_FACE=2,BACK_FACE=3,BOTTOM_FACE=4,UP_FACE=5};
00057 #endif
00058
00059 class OrthoCellAccessor;
00060 class OrthoCellAccessorWithHoles;
00061 class OrthoFaceAccessor;
00062 class OrthoFaceAccessorWithHoles;
00063 class OrthoVerticeAccessor;
00064 class OrthoVerticeAccessorWithHoles;
00065 class CashOrthoFaceAccessor;
00066 class FaceInfoCash;
00067 class UnitTests;
00068
00069
00070 template<class Acc>
00071 class IteratorAccessor
00072 {
00073 protected:
00074 Acc m_data;
00075
00076 public:
00077 IteratorAccessor(){}
00078 IteratorAccessor(const Acc &acc):m_data(acc){}
00079 Acc& operator*(){return m_data;}
00080 bool operator!=(IteratorAccessor<Acc> it) {return (*it) != m_data;}
00081 Acc* operator->(){return &m_data;}
00082 const Acc* operator->() const {return &m_data;}
00083 void operator++(int) {m_data++;}
00084
00085 };
00086
00087
00095 template<class Acc>
00096 class InnerIteratorAccessor : public IteratorAccessor<Acc>
00097 {
00098 public:
00099
00100 InnerIteratorAccessor(Acc& acc):IteratorAccessor<Acc>(acc){}
00101 void operator++(int) {this->m_data.advance_inner();}
00102 };
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00159 class OrthoMesh
00160 {
00161 friend class OrthoCellAccessor;
00162 friend class OrthoCellAccessorWithHoles;
00163 friend class OrthoFaceAccessorWithHoles;
00164 friend class OrthoFaceAccessor;
00165 friend class OrthoVerticeAccessor;
00166 friend class OrthoVerticeAccessorWithHoles;
00167 friend class FaceInfoCash;
00168 friend class CashOrthoFaceAccessor;
00169 friend class UnitTests;
00170
00171
00172 private:
00173 VecWellInfo m_wells;
00174 const static unsigned _VDirTbl[][3];
00175 Point3D _p1,_p2;
00176 Point3D _DX,_Barycenter0;
00177 double _dx,_dy,_dz;
00178 Point3D _fxP0,_fyP0,_fzP0;
00179 unsigned _nElemX,_nElemY,_nElemZ;
00180 unsigned _nElemXY,_nElemYZ,_nElemXZ;
00181 unsigned _numRawCells;
00182 unsigned _numInnerRawCells;
00183 unsigned _numRawVertices;
00184 unsigned _numRawFaces;
00185 unsigned _numCells;
00186 unsigned _numVertices;
00187 unsigned _numFaces;
00188 unsigned _numHGhostCells;
00189 unsigned _numRawGhostCells;
00190 unsigned _ZStride;
00191 unsigned _YStride;
00192 unsigned _ZVStride;
00193 unsigned _YVStride;
00194 unsigned _lastI,_lastJ,_lastK;
00195 unsigned _fxYStride,_fxZStride,_fyYStride,_fyZStride,_fzYStride,_fzZStride;
00196 unsigned _fyOff,_fzOff;
00197 unsigned _numFacesAtBoundary;
00198 unsigned _LastI,_LastJ,_LastK;
00199 unsigned _GhostLeftOff,_GhostRightOff,_GhostBottomOff,_GhostUpOff,_GhostFrontOff,_GhostBackOff;
00200 unsigned _HGhostLeftOff,_HGhostRightOff,_HGhostBottomOff,_HGhostUpOff,_HGhostFrontOff,_HGhostBackOff;
00201 unsigned _HGhostWellOff;
00202
00203 double _faceArea[3];
00204 double _faceAreaPerCellVol[3];
00205 double _cellVolume;
00206 double _grid_tol;
00207
00208
00209 Index m_innerFacesChangeNormalIndex[3];
00210 Index m_innerFacesChangeKOffset[3];
00211
00212
00213 OrthoCellAccessor *_pBeginCell,*_pEndCell;
00214 OPointer<OrthoCellAccessor> _pBeginInnerCell,_pEndInnerCell;
00215 OrthoFaceAccessor *_pBeginFace,*_pEndFace,*_pBeginInnerFace;
00216 OrthoVerticeAccessor *_pBeginVert,*_pEndVert;
00217 OrthoVerticeAccessorWithHoles *_pBeginVertH,*_pEndVertH;
00218 OrthoCellAccessorWithHoles *_pBeginCellH,*_pEndCellH;
00219 OrthoFaceAccessorWithHoles *_pBeginFaceH,*_pEndFaceH;
00220
00221 DelIndexLst _delFaceLst;
00222 DelIndexLst _delVertLst;
00223 DelIndexLst _delCellLst;
00224
00225
00226
00227
00228 std::vector<FaceInfoCash> m_vFCash;
00229 typedef std::vector<FaceInfoCash>::const_iterator FaceCellsTbl_It;
00230 private:
00231 void initOrthoMesh(Point3D &p1,Point3D &p2,unsigned nElemX,unsigned nElemY,unsigned nElemZ);
00232
00233
00234
00235 public:
00236 enum NORMAL_AXIS {NORMAL_X=0,NORMAL_Y=1,NORMAL_Z=2};
00237 const static unsigned FACES_PER_CELL;
00238 const static unsigned VERTICES_PER_CELL;
00239 const static unsigned INVALID_INDEX;
00240 typedef IteratorAccessor<OrthoCellAccessorWithHoles> Cell_It;
00241 typedef IteratorAccessor<OrthoFaceAccessorWithHoles> Face_It;
00242 typedef IteratorAccessor<OrthoCellAccessor> Raw_Cell_It;
00243 typedef IteratorAccessor<OrthoFaceAccessor> Raw_Face_It;
00244 typedef IteratorAccessor<OrthoVerticeAccessor> Raw_Vertex_It;
00245 typedef IteratorAccessor<OrthoVerticeAccessorWithHoles> Vertex_It;
00246 typedef IteratorAccessor<CashOrthoFaceAccessor> Cash_Face_It;
00247
00248 typedef InnerIteratorAccessor<OrthoCellAccessor> Inner_Raw_Cell_It;
00249 typedef InnerIteratorAccessor<OrthoFaceAccessor> Inner_Raw_Face_It;
00250
00251
00252 protected:
00253 unsigned getCellIndexFromPosition(unsigned i,unsigned j,unsigned k) const;
00254 void getCellPositionFromIndex(unsigned index,unsigned &i,unsigned &j,unsigned &k) const;
00255 void getVerticePositionFromIndex(unsigned index,unsigned &i,unsigned &j,unsigned &k) const;
00256 unsigned getVerticeIndexFromPosition(unsigned i,unsigned j,unsigned k) const;
00257 void getFacePositionFromIndex(unsigned index,unsigned &i,unsigned &j,unsigned &k,NORMAL_AXIS &normal) const;
00258 Index getFaceIndexFromPosition(unsigned i,unsigned j,unsigned k,NORMAL_AXIS normal) const;
00259 void getBarycenterFromCell(unsigned i,unsigned j,unsigned k,Point3D &p) const;
00260 void getVerticeFromCell(VertexDirection3D dir, unsigned i,unsigned j,unsigned k,Point3D &p) const;
00261 unsigned getVerticeIndexFromCell(VertexDirection3D dir,unsigned i,unsigned j,unsigned k) const;
00262 unsigned getFaceIndexFromCell(FaceDirection3D dir,unsigned i,unsigned j,unsigned k) const;
00263 unsigned getFaceIndexFromCell(FaceDirection3D dirEnum,unsigned i,unsigned j,unsigned k,DelIndexLst::Iterator &itX,DelIndexLst::Iterator &itY,DelIndexLst::Iterator &itZ) const;
00264 void getFaceInfoFromCell(FaceDirection3D dirEnum,unsigned ci,unsigned cj,unsigned ck,
00265 unsigned *fi,unsigned *fj,unsigned *fk,unsigned *findex,NORMAL_AXIS *fnormal) const;
00266
00267 bool isCellAtBoundary(unsigned i,unsigned j,unsigned k) const;
00268 static void advanceIJK(unsigned &i,unsigned &j,unsigned &k,unsigned minI,unsigned maxI,unsigned minJ,unsigned maxJ);
00269 Face_It get_face_from_raw_index(unsigned raw_index) const;
00270
00271
00272 private:
00273
00274
00275 Index m_FaceStridesTbl[3][2];
00276
00277
00278 public:
00279 OrthoMesh(Point3D p1,Point3D p2,unsigned nElemX,unsigned nElemY,unsigned nElemZ);
00280 OrthoMesh(Point3D &p1,Point3D &p2,unsigned nElemX,unsigned nElemY,unsigned nElemZ,VecWellInfo &wells);
00281 void putWells(VecWellInfo &wells);
00282 void buildCashInfo();
00283 ~OrthoMesh();
00284 unsigned numRawCells() const {return _numRawCells;}
00285 unsigned numRawFaces() const {return _numRawFaces;}
00286 unsigned numRawVertices() const {return _numRawVertices;}
00287 unsigned numCells() const {return _numCells;}
00288 unsigned numFaces() const {return _numFaces;}
00289 unsigned numVertices() const {return _numVertices;}
00290 unsigned numGhostCells() const {return _numHGhostCells;}
00291 unsigned numFacesAtBoundary() const {return _numFacesAtBoundary;}
00292 unsigned numElemX() const {return _nElemX;}
00293 unsigned numElemY() const {return _nElemY;}
00294 unsigned numElemZ() const {return _nElemZ;}
00295 Point3D numElems() const {return Point3D(_nElemX,_nElemY,_nElemZ);}
00296 Point3D getP() const {return _p1;}
00297 Point3D getQ() const {return _p2;}
00298 void print(std::string fileName="");
00299 static bool isValid(unsigned index){return index != INVALID_INDEX;}
00300
00301 Point3D getDX() const {return _DX;}
00302 double get_dx() const {return _dx;}
00303 double get_dy() const {return _dy;}
00304 double get_dz() const {return _dz;}
00305 VecWellInfo getWells(){return m_wells;}
00306 Face_It begin_face() const;
00307 Face_It get_face(unsigned index);
00308 Cash_Face_It begin_cash_face();
00309 Cash_Face_It end_cash_face() const;
00310
00311 Inner_Raw_Face_It begin_inner_raw_face();
00312 Inner_Raw_Face_It end_inner_raw_face();
00313
00314
00315 Face_It end_face() const;
00316 Cell_It begin_cell() const;
00317 Cell_It get_cell(unsigned index) const;
00318 Cell_It end_cell() const;
00319 Vertex_It begin_vertice() const;
00320 Vertex_It end_vertice() const;
00321 Vertex_It get_vertice(unsigned index) const;
00322 Raw_Cell_It begin_raw_cell() const;
00323 Inner_Raw_Cell_It begin_inner_raw_cell() const;
00324 Inner_Raw_Cell_It end_inner_raw_cell() const;
00325
00326 Raw_Cell_It end_raw_cell() const;
00327 Raw_Face_It begin_raw_face() const;
00328 Raw_Face_It end_raw_face() const;
00329
00330
00331 Raw_Vertex_It begin_raw_vertice() const;
00332 Raw_Vertex_It end_raw_vertice() const;
00333
00334 Vertex_It getNearestVertice(Point3D p);
00335
00336 static unsigned invalidIndex(){return INVALID_INDEX;}
00337 bool nonZeroInwardNormalComponenentIsPositive(FaceDirection3D face_no) const;
00338 void projectCentralValuesAtVertices(const VecDouble &cValues,VecDouble &vValues) const;
00339 void projectVelocitiesInFacesToCells(Matrix &vel,const VecDouble &vNormalFaceVel);
00340 void projectVelocitiesInFacesToVertices(Matrix &vel,const VecDouble &vNormalFaceVel);
00341 void getMaxAbsVelocitiesInFacesByComponents(VecDouble &vel, const VecDouble &vNormalFaceVel);
00342
00343 double cellVolume() const {return _cellVolume;}
00344 double cellMeasure() const {return cellVolume();}
00345 double meshVolume() const {return numCells()*cellVolume();}
00346
00347 void setCentralValuesFromFunction(Function3D &f,VecDouble &cValues,unsigned component=0);
00348 void setCentralValuesFromFunction(Function3D &f,Matrix &MValuess);
00349 void setFacesValuesFromFunction(Function3D &f,VecDouble &fValues,unsigned deg=0,bool bOnlyPrescribedBoundary=false);
00350 void setFacesValuesFromFunction(Function3D &f,Matrix &MValues,const bool bOnlyAtBoundary=false);
00351
00352 double getIntegralAtCells(const VecDouble &cValues) const;
00353 double getIntegralAtVertices(const VecDouble &vValues) const;
00354
00355
00356
00357 void getFacesInSlab(double Const,VecIndex &fIndices,unsigned dim);
00358 void getFacesInYZSlab(double X,VecIndex &fIndices);
00359
00360
00361 void getCellsIndicesInDomain(Function3D &f,VecIndex &vec);
00362 void tagFacesInDomain(Function3D &f,VecTag &tags,char value);
00363 void getFacesIndicesInDomain(Function3D &f, VecIndex &vec);
00364 void printCells(VecIndex &vec,std::ostream &out);
00365 void printCell(Cell_It cell,std::ostream &out);
00366 double getGridTolerance(){return _grid_tol;}
00367 unsigned getFacesInYOffset() {return _fyOff;}
00368 unsigned getFacesInZOffset() {return _fzOff;}
00369
00370 Cell_It getCellAt(Point3D &p);
00371
00372
00373 static double getOutwardNormalOrientation(FaceDirection3D faceDir);
00374
00375
00376 };
00377
00378
00379 #include "orthocellaccessor.h"
00380 #include "orthofaceaccessor.h"
00381 #include "orthoverticeaccessor.h"
00382 #include "orthoverticeaccessorwithholes.h"
00383 #include "orthocellaccessorwithholes.h"
00384 #include "orthofaceaccessorwithholes.h"
00385 #include "faceinfocash.h"
00386 #include "cashorthofaceaccessor.h"
00387 #include "delindexlst.h"
00388
00389
00390
00391 #endif
00392