CO2INJECTION: orthomesh.cpp Source File

00001 #include "orthomesh.h"
00002 #include <limits.h>
00003 #include "orthocellaccessor.h"
00004 #include "orthofaceaccessor.h"
00005 #include <stdio.h>
00006 #include "sfunctions.h"
00007 #include "numericmethods.h"
00008 #include <ostream>
00009 
00010 const unsigned OrthoMesh::_VDirTbl[][3] = {{0,0,0},
00011                                                 {1,0,0},
00012                                                 {0,1,0},
00013                                                 {1,1,0},
00014                                                 {0,0,1},
00015                                                 {1,0,1},
00016                                                 {0,1,1},
00017                                                 {1,1,1}};
00018 const unsigned OrthoMesh::FACES_PER_CELL=6;
00019 const unsigned OrthoMesh::VERTICES_PER_CELL=8;
00020                                        
00021 const unsigned OrthoMesh::INVALID_INDEX = UINT_MAX;
00022 
00023 
00024 void OrthoMesh::initOrthoMesh(Point3D &p1,Point3D &p2,unsigned nElemX,unsigned nElemY,unsigned nElemZ) 
00025  {
00026    _p1=p1;
00027    _p2=p2;
00028    _dx = _DX[0] = (_p2[0]-_p1[0])/((double) nElemX);
00029    _dy = _DX[1] = (_p2[1]-_p1[1])/((double) nElemY);
00030    _dz = _DX[2] = (_p2[2]-_p1[2])/((double) nElemZ);
00031    _cellVolume = _DX[0]*_DX[1]*_DX[2];
00032    
00033    _Barycenter0[0] =  _p1[0]+_DX[0]/2.0;
00034    _Barycenter0[1] =  _p1[1]+_DX[1]/2.0;
00035    _Barycenter0[2] =  _p1[2]+_DX[2]/2.0;
00036 
00037    _nElemX = nElemX;
00038    _nElemY = nElemY;
00039    _nElemZ = nElemZ;
00040 
00041    _LastI = nElemX-1;
00042    _LastJ = nElemY-1;
00043    _LastK = nElemZ-1;
00044 
00045    _nElemXY= nElemX*nElemY;
00046    _nElemXZ= nElemX*nElemZ;
00047    _nElemYZ= nElemY*nElemZ;
00048    
00049   
00050    
00051    _numRawCells = nElemX*nElemY*nElemZ;
00052    _numRawVertices=(nElemX+1)*(nElemY+1)*(nElemZ+1);
00053    _numInnerRawCells = ((nElemX > 2) ? nElemX-2 : 0)*((nElemY > 2) ? nElemY-2 : 0)*((nElemZ > 2) ? nElemZ-2 : 0);
00054 
00055 
00056    _ZStride = nElemX*nElemY;
00057    _YStride = nElemX;
00058 
00059 
00060 
00061    _ZVStride = (nElemX+1)*(nElemY+1);
00062    _YVStride = (nElemX+1);
00063 
00064    _fxYStride = nElemX+1;
00065    _fxZStride = (nElemX+1)*nElemY;
00066 
00067    _fyYStride = nElemX;
00068    _fyZStride = (nElemY+1)*nElemX;
00069 
00070    _fzYStride = nElemX;
00071    _fzZStride = nElemX*nElemY;
00072    
00073    _fyOff = (nElemX+1)*nElemY*nElemZ;
00074    _fzOff = _fyOff + nElemX*(nElemY+1)*nElemZ;
00075    _numRawFaces = (nElemX+1)*nElemY*nElemZ + nElemX*(nElemY+1)*nElemZ + nElemX*nElemY*(nElemZ+1);
00076 
00077    _numFaces    = _numRawFaces;
00078    _numVertices = _numRawVertices;
00079    _numCells    = _numRawCells;
00080    
00081    _numFacesAtBoundary = 2*(nElemX*nElemY + nElemX*nElemZ + nElemZ*nElemY); 
00082    
00083    _lastI = nElemX-1;
00084    _lastJ = nElemY-1;
00085    _lastK = nElemZ-1;
00086 
00087    _fxP0=_p1;
00088    _fxP0[0] += 0.0;
00089    _fxP0[1] += _DX[1]/2.0;
00090    _fxP0[2] += _DX[2]/2.0;
00091 
00092    _fyP0=_p1;
00093    _fyP0[0] += _DX[0]/2.0;
00094    _fyP0[1] += 0.0;
00095    _fyP0[2] += _DX[2]/2.0;
00096 
00097    _fzP0=_p1;
00098    _fzP0[0] += _DX[0]/2.0;
00099    _fzP0[1] += _DX[1]/2.0;
00100    _fzP0[2] += 0.0;
00101 
00102 
00103    _faceArea[NORMAL_X]=_DX[1]*_DX[2];
00104    _faceArea[NORMAL_Y]=_DX[0]*_DX[2];
00105    _faceArea[NORMAL_Z]=_DX[0]*_DX[1];
00106 
00107    _faceAreaPerCellVol[NORMAL_X]=1.0/_DX[0];
00108    _faceAreaPerCellVol[NORMAL_Y]=1.0/_DX[1];
00109    _faceAreaPerCellVol[NORMAL_Z]=1.0/_DX[2];
00110 
00111    //close the lists of deleted faces,cells and vertices
00112    _delFaceLst.close();
00113    _delCellLst.close();
00114    _delVertLst.close();
00115 
00116    _grid_tol = std::min(std::min(_DX[0],_DX[1]),_DX[2])/4.0;
00117    
00118 
00119 
00120    //Set the offsets to calculate indices of ghost cells. See OrthoFaceAccessor::ghostAdjCellIndices();
00121    _GhostLeftOff  = _numRawCells;
00122    _GhostRightOff = _GhostLeftOff  + _nElemYZ;
00123    _GhostFrontOff =_GhostRightOff  + _nElemYZ;
00124    _GhostBackOff  =_GhostFrontOff  + _nElemXZ;
00125    _GhostBottomOff= _GhostBackOff  + _nElemXZ;
00126    _GhostUpOff    = _GhostBottomOff+ _nElemXY;
00127 
00128    _numRawGhostCells = 2*(_nElemYZ + _nElemXY + _nElemXZ);
00129 
00130    _HGhostLeftOff   = _GhostLeftOff;
00131    _HGhostRightOff  = _GhostRightOff; 
00132    _HGhostFrontOff  = _GhostFrontOff;
00133    _HGhostBackOff   = _GhostBackOff;
00134    _HGhostBottomOff = _GhostBottomOff;
00135    _HGhostUpOff     = _GhostUpOff; 
00136 
00137    _numHGhostCells = _numRawGhostCells;
00138    
00139   
00140    
00141    //Set the iterators of the first and last elements
00142    //for the three geometries (vertices,cells,faces)
00143    _pBeginCell  = new OrthoCellAccessor(*this,0); 
00144    _pEndCell = new OrthoCellAccessor(*this,this->numRawCells()-1);
00145    (*_pEndCell)++;
00146    
00147    
00148    _pBeginFace  = new OrthoFaceAccessor(*this,0); 
00149    _pEndFace = new OrthoFaceAccessor(*this,this->numRawFaces()-1);
00150    (*_pEndFace)++;
00151 
00152    
00153 
00154 
00155    _pBeginVert  = new OrthoVerticeAccessor(*this,0); 
00156    _pEndVert = new OrthoVerticeAccessor(*this,this->numRawVertices()-1);
00157    (*_pEndVert)++;
00158 
00159 
00160 
00161 
00162    _pBeginCellH  = new OrthoCellAccessorWithHoles(*this,0); 
00163    _pEndCellH = new OrthoCellAccessorWithHoles(*this,this->numRawCells()-1);
00164    (*_pEndCellH)++;
00165 
00166    _pBeginFaceH  = new OrthoFaceAccessorWithHoles(*this,0); 
00167    _pEndFaceH = new OrthoFaceAccessorWithHoles(*this,this->numRawFaces()-1);
00168    (*_pEndFaceH)++;
00169 
00170 
00171    _pBeginVertH  = new OrthoVerticeAccessorWithHoles(*this,0); 
00172    _pEndVertH = new OrthoVerticeAccessorWithHoles(*this,this->numVertices()-1);
00173    (*_pEndVertH)++;
00174 
00175 
00176    //Set usefull tables. This tables are to avoid if then else conditions.
00177    //Some of theses tables are used by one or just two methods and the best
00178    //thing they are small and does not depende on the size of the mesh.
00179 
00180    //   m_faceIndexTable[NORMAL_X]
00181    m_FaceStridesTbl[NORMAL_X][0]=_fxYStride;   
00182    m_FaceStridesTbl[NORMAL_X][1]=_fxZStride;   
00183    m_FaceStridesTbl[NORMAL_Y][0]=_fyYStride;   
00184    m_FaceStridesTbl[NORMAL_Y][1]=_fyZStride;   
00185    m_FaceStridesTbl[NORMAL_Z][0]=_fzYStride;   
00186    m_FaceStridesTbl[NORMAL_Z][1]=_fzZStride;
00187 
00188 
00189    
00190 
00191    //This array is used in the OrthoFaceAccessor::advance_inner method 
00192    //m_innerFacesChangeKOffset is the  value of the offset you
00193    //need to add to the indice to go to the next J position.
00194    //m_innerFacesChangeKOffset is an array with offset for
00195    //each of the faces normal type (faces with normal_x, normal_y,normal_z)
00196    m_innerFacesChangeKOffset[NORMAL_X]=3;
00197    m_innerFacesChangeKOffset[NORMAL_Y]=2*(_nElemX)+1;
00198    m_innerFacesChangeKOffset[NORMAL_Z]=1;
00199 
00200 
00201    //This array is used in the OrthoFaceAccessor::advance_inner method 
00202    //It basically contains the indice of the next inner face when we
00203    //have an iterator at the last inner face with normal_x, normal_y
00204    //normal_z. So m_innerFacesChangeNormalOffset[NORMAL_X] is the indice
00205    //of the next inner face following the last inner face with normal_x.
00206    //That will be the indice of the face 1,1,1 with normal Y.
00207    //In the case of Z the next indice will be the invalid
00208    //cell with indice equal to _numRawFaces.
00209    m_innerFacesChangeNormalIndex[NORMAL_X]=_fyOff+_nElemX;
00210    m_innerFacesChangeNormalIndex[NORMAL_Y]=_fzOff+_nElemXY;
00211    m_innerFacesChangeNormalIndex[NORMAL_Z]=_numRawFaces;
00212    
00213 
00214    OrthoMesh::Raw_Cell_It cell = begin_raw_cell();
00215    while(cell != end_raw_cell() && cell->at_boundary())
00216      cell++;
00217    _pBeginInnerCell = new OrthoCellAccessor(*cell);
00218    _pEndInnerCell = new OrthoCellAccessor(*this,0,0,_LastK);
00219 
00220    OrthoMesh::Raw_Face_It face = begin_raw_face();
00221    OrthoMesh::Raw_Face_It endf = end_raw_face();
00222 
00223    while(face != end_raw_face() && face->at_boundary()) 
00224      face++;
00225    _pBeginInnerFace = new OrthoFaceAccessor(*face);
00226    
00227 
00228  }
00229 
00230 OrthoMesh::OrthoMesh(Point3D p1,Point3D p2,unsigned nElemX,unsigned nElemY,unsigned nElemZ)
00231 {
00232   initOrthoMesh(p1,p2,nElemX,nElemY,nElemZ);
00233 }
00234 
00235 OrthoMesh::OrthoMesh(Point3D &p1,Point3D &p2,unsigned nElemX,unsigned nElemY,unsigned nElemZ,VecWellInfo &wells)
00236 {
00237   initOrthoMesh(p1,p2,nElemX,nElemY,nElemZ);
00238   putWells(wells);
00239 }
00240 
00241 OrthoMesh::~OrthoMesh() 
00242 {
00243   delete _pBeginCell;
00244   delete _pEndCell;
00245   delete _pEndFace;
00246   delete _pBeginFace;
00247 }
00248 
00249 void OrthoMesh::getCellPositionFromIndex(unsigned index,unsigned &i,unsigned &j,unsigned &k) const 
00250 {
00251   assert(index < numRawCells()); 
00252   k = index/_ZStride;
00253   index = index%_ZStride;
00254   j = index/_YStride;
00255   i = index%_YStride;
00256 
00257 }
00258 
00259 unsigned OrthoMesh::getCellIndexFromPosition(unsigned i,unsigned j,unsigned k)  const
00260 {
00261   assert( i < _nElemX && j < _nElemY && k < _nElemZ); 
00262   return k*_ZStride + j*_YStride + i;
00263   
00264 }
00265 
00266 
00267 
00268 void OrthoMesh::getVerticePositionFromIndex(unsigned index,unsigned &i,unsigned &j,unsigned &k) const 
00269 {
00270   assert(index<numRawVertices());
00271   k=index/_ZVStride;
00272   index=index%_ZVStride;
00273   j=index/_YVStride;
00274   i=index%_YVStride;
00275   
00276 }
00277 
00278 unsigned OrthoMesh::getVerticeIndexFromPosition(unsigned i,unsigned j,unsigned k) const 
00279 {
00280   return i + j*_YVStride + k*_ZVStride;
00281 }
00282 
00283 
00284 
00285 unsigned OrthoMesh::getVerticeIndexFromCell(VertexDirection3D dir,unsigned i,unsigned j,unsigned k) const
00286 {
00287   const unsigned *vDir = &(_VDirTbl[dir][0]); 
00288   return (i+vDir[0]) + (j+vDir[1])*(_YVStride) + (k+vDir[2])*_ZVStride;
00289   
00290 }
00291 
00292 void OrthoMesh::getVerticeFromCell(VertexDirection3D dir, unsigned i,unsigned j,unsigned k,Point3D &p) const
00293 {
00294   const unsigned *vDir = &(_VDirTbl[dir][0]); 
00295   p[0]=(i+vDir[0])*_DX[0];
00296   p[1]=(j+vDir[1])*_DX[1];
00297   p[2]=(k+vDir[2])*_DX[2];
00298 }
00299 
00300 void OrthoMesh::getBarycenterFromCell(unsigned i,unsigned j,unsigned k,Point3D &p) const
00301 {
00302   p[0]= _Barycenter0[0] +  i*_DX[0];
00303   p[1]= _Barycenter0[1] +  j*_DX[1];
00304   p[2]= _Barycenter0[2] +  k*_DX[2];
00305 }
00306 
00307 
00308 
00309 
00310 
00311 
00312 
00313 
00314 
00315 
00316 
00317 
00342 unsigned OrthoMesh::getFaceIndexFromCell(FaceDirection3D dirEnum,unsigned i,unsigned j,unsigned k)  const
00343 {
00344   unsigned dir = (unsigned) dirEnum;
00345   unsigned result;
00346   if (dir < 2)
00347   {
00348     result = i + dir  + j*_fxYStride + k*_fxZStride;
00349   }
00350   else if (dir < 4)
00351   {
00352     result = i + (j + dir - 2)*_fyYStride + k*_fyZStride;
00353     result += _fyOff;
00354   }
00355   else
00356   {
00357     result = i + j*_fzYStride + (k +dir-4)*_fzZStride;
00358     result+= _fzOff;
00359   }
00360   return result;
00361 }
00362 
00387 void OrthoMesh::getFaceInfoFromCell(FaceDirection3D dirEnum,unsigned ci,unsigned cj,unsigned ck,
00388                                     unsigned *fi,unsigned *fj,unsigned *fk,unsigned *findex,NORMAL_AXIS *fnormal)  const
00389 {
00390   unsigned dir = (unsigned) dirEnum;
00391   if (dir < 2)
00392   {
00393     *findex = ci + dir  + cj*_fxYStride + ck*_fxZStride;
00394     *fi=ci+dir;
00395     *fj=cj;
00396     *fk=ck;
00397     *fnormal=NORMAL_X;
00398   }
00399   else if (dir < 4)
00400   {
00401     *findex = ci + (cj + dir - 2)*_fyYStride + ck*_fyZStride;
00402     *findex += _fyOff;
00403     *fi=ci;
00404     *fj=cj+dir-2;
00405     *fk=ck;
00406     *fnormal=NORMAL_Y;
00407   }
00408   else
00409   {
00410     *findex = ci + cj*_fzYStride + (ck +dir-4)*_fzZStride;
00411     *findex+= _fzOff;
00412     *fi=ci;
00413     *fj=cj;
00414     *fk=ck+dir-4;
00415     *fnormal=NORMAL_Z;
00416   }
00417 }
00418 
00419 
00420 
00437 unsigned OrthoMesh::getFaceIndexFromCell(FaceDirection3D dirEnum,unsigned i,unsigned j,unsigned k,
00438                                          DelIndexLst::Iterator &itX,DelIndexLst::Iterator &itY,DelIndexLst::Iterator &itZ)  const
00439 {
00440   unsigned dir = (unsigned) dirEnum;
00441   if (dir < 2)
00442   {
00443     return  _delFaceLst.adjust(i + dir  + j*_fxYStride + k*_fxZStride,itX);
00444 
00445   }
00446   else if (dir < 4)
00447   {
00448     return  _delFaceLst.adjust(i + (j + dir - 2)*_fyYStride + k*_fyZStride + _fyOff,itY);
00449   }
00450   else
00451   {
00452     return _delFaceLst.adjust( i + j*_fzYStride + (k +dir-4)*_fzZStride + _fzOff,itZ);
00453   }
00454 }
00455 
00456 
00457 
00458 
00459 bool OrthoMesh::isCellAtBoundary(unsigned i,unsigned j,unsigned k) const
00460 {
00461   return (i==0 || j ==0 || k==0 || i == _lastI || j == _lastJ || k == _lastK);
00462 }
00463 
00464 
00465 void OrthoMesh::getFacePositionFromIndex(unsigned index,unsigned &i,unsigned &j,unsigned &k,NORMAL_AXIS &normal) const
00466 {
00467   assert(index < numRawFaces());
00468   unsigned rm;
00469   if (index < _fyOff)
00470   {
00471     k  = index/_fxZStride;
00472     rm = index%_fxZStride; 
00473     j = rm/_fxYStride;
00474     i = rm%_fxYStride;
00475     normal = NORMAL_X;
00476   }
00477   else if (index < _fzOff)
00478   {
00479     index-=_fyOff;
00480     k  = index/_fyZStride;
00481     rm = index%_fyZStride; 
00482     j = rm/_fyYStride;
00483     i = rm%_fyYStride;
00484     normal = NORMAL_Y;
00485   }
00486   else
00487   {
00488     index-=_fzOff;
00489     k  = index/_fzZStride;
00490     rm = index%_fzZStride; 
00491     j = rm/_fzYStride;
00492     i = rm%_fzYStride;
00493     normal = NORMAL_Z;
00494   }
00495   
00496 }
00497 
00498 Index OrthoMesh::getFaceIndexFromPosition(unsigned i,unsigned j,unsigned k,NORMAL_AXIS normal) const 
00499 {
00500   return i+m_FaceStridesTbl[normal][0]*j + m_FaceStridesTbl[normal][1]*k;
00501   
00502   
00503 }
00504 
00505 
00506 /*
00507 void OrthoMesh::getBarycenterFromFace(OrthoFaceAccessor &face,Point3D &p) 
00508 {
00509   switch (face._normal)
00510   {
00511   case NORMAL_X:
00512     p[0]=_fxP0[0] + face._i*_DX[0];
00513     p[1]=_fxP0[1] + face._j*_DX[1];
00514     p[2]=_fxP0[2] + face._k*_DX[2];
00515     break;
00516   case NORMAL_Y:
00517     p[0]=_fyP0[0] + face._i*_DX[0];
00518     p[1]=_fyP0[1] + face._j*_DX[1];
00519     p[2]=_fyP0[2] + face._k*_DX[2];
00520     break;
00521   case NORMAL_Z:
00522     p[0]=_fzP0[0] + face._i*_DX[0];
00523     p[1]=_fzP0[1] + face._j*_DX[1];
00524     p[2]=_fzP0[2] + face._k*_DX[2];
00525     break;
00526   }
00527 }
00528 
00529 
00530 
00531 
00532 
00533 
00534 
00535 
00536 
00537 
00538 
00539 */
00552 void OrthoMesh::advanceIJK(unsigned &i,unsigned &j,unsigned &k,const unsigned minI,const unsigned maxI,const unsigned minJ,const unsigned maxJ)  
00553 {
00554   if (i != maxI)
00555   {
00556     i++;
00557   }
00558   else
00559   {
00560     i=minI;
00561     if (j != maxJ)
00562     {
00563       j++;
00564     }
00565     else{
00566       j=minJ;
00567       k++;
00568     }
00569   }
00570 }
00571 
00572 
00573 
00574 
00575 OrthoMesh::Cell_It OrthoMesh::begin_cell() const 
00576 {
00577   return Cell_It(*_pBeginCellH);
00578 }
00579 
00580 
00581 OrthoMesh::Inner_Raw_Cell_It OrthoMesh::begin_inner_raw_cell() const 
00582 {
00583   assert(_pBeginInnerCell);
00584   return Inner_Raw_Cell_It(*_pBeginInnerCell);
00585 }
00586 
00587 OrthoMesh::Inner_Raw_Cell_It OrthoMesh::end_inner_raw_cell() const 
00588 {
00589   assert(_pEndInnerCell);
00590   return Inner_Raw_Cell_It(*_pEndInnerCell);
00591 }
00592 
00593 
00594 
00595 
00596 
00597 
00598 
00599 
00600 
00601 OrthoMesh::Cell_It OrthoMesh::end_cell()  const
00602 {
00603   return Cell_It(*_pEndCellH);
00604 }
00605 
00606 
00607 OrthoMesh::Raw_Cell_It OrthoMesh::begin_raw_cell() const 
00608 {
00609   return Raw_Cell_It(*_pBeginCell);
00610 }
00611 
00612 OrthoMesh::Raw_Cell_It OrthoMesh::end_raw_cell() const 
00613 {
00614   return Raw_Cell_It(*_pEndCell);
00615 }
00616 
00617 OrthoMesh::Raw_Vertex_It OrthoMesh::begin_raw_vertice() const 
00618 {
00619   return Raw_Vertex_It(*_pBeginVert);
00620 }
00621 
00622 OrthoMesh::Raw_Vertex_It OrthoMesh::end_raw_vertice() const 
00623 {
00624   return Raw_Vertex_It(*_pEndVert);
00625 }
00626 
00627 
00628 
00629 OrthoMesh::Face_It OrthoMesh::begin_face() const 
00630 {
00631   return Face_It(*_pBeginFaceH);
00632   
00633 }
00638 OrthoMesh::Face_It OrthoMesh::get_face_from_raw_index(unsigned raw_index) const 
00639 {
00640   OrthoFaceAccessorWithHoles acc(*this,raw_index);
00641   return Face_It(acc);
00642     
00643 }
00644 
00649 OrthoMesh::Face_It OrthoMesh::get_face(unsigned index) 
00650 {
00651   OrthoFaceAccessorWithHoles acc(*this,_delFaceLst.index_to_raw_index(index));
00652   return Face_It(acc);
00653 }
00659 OrthoMesh::Cell_It OrthoMesh::get_cell(unsigned index)  const
00660 {
00661   OrthoCellAccessorWithHoles acc(*this,_delCellLst.index_to_raw_index(index));
00662   return Cell_It(acc);
00663 
00664 }
00669 OrthoMesh::Vertex_It OrthoMesh::get_vertice(unsigned index) const 
00670 {
00671   OrthoVerticeAccessorWithHoles acc(*this,_delVertLst.index_to_raw_index(index));
00672   return Vertex_It(acc);
00673 
00674 }
00675 
00676 OrthoMesh::Face_It OrthoMesh::end_face() const 
00677 {
00678   return Face_It(*_pEndFaceH);
00679 }
00680 
00681 OrthoMesh::Raw_Face_It OrthoMesh::begin_raw_face() const 
00682 {
00683   return Raw_Face_It(*_pBeginFace);
00684 }
00685 
00686 OrthoMesh::Inner_Raw_Face_It OrthoMesh::begin_inner_raw_face() 
00687 {
00688   return Inner_Raw_Face_It(*_pBeginInnerFace);
00689 }
00690 
00691 OrthoMesh::Inner_Raw_Face_It OrthoMesh::end_inner_raw_face() 
00692 {
00693   return Inner_Raw_Face_It(*_pEndFace);
00694 }
00695 
00696 
00697 
00698 
00699 
00700 OrthoMesh::Raw_Face_It OrthoMesh::end_raw_face() const 
00701 {
00702   return Raw_Face_It(*_pEndFace);
00703 }
00704 
00705 OrthoMesh::Cash_Face_It OrthoMesh::begin_cash_face()  
00706 {
00707   if (m_vFCash.size() ==0)
00708   {
00709     this->buildCashInfo();
00710   }
00711   return CashOrthoFaceAccessor(*this,m_vFCash.begin());
00712 }
00713 
00714 OrthoMesh::Cash_Face_It OrthoMesh::end_cash_face() const 
00715 {
00716   assert(m_vFCash.size() !=0);
00717   return CashOrthoFaceAccessor(*this,m_vFCash.end());
00718 }
00719 
00720 
00721 OrthoMesh::Vertex_It OrthoMesh::begin_vertice() const 
00722 {
00723   return Vertex_It(*_pBeginVertH);
00724 }
00725 
00726 OrthoMesh::Vertex_It OrthoMesh::end_vertice() const 
00727 {
00728   return Vertex_It(*_pEndVertH);
00729 }
00730 
00731 
00732 
00736 void OrthoMesh::print(std::string fileName) 
00737 {
00738   
00739   FILE *file;
00740   if (fileName.empty())
00741     file=stdout;
00742   else
00743     file=fopen(fileName.c_str(),"w");
00744 
00745   Cell_It cell = begin_cell();
00746   Cell_It endc = end_cell();
00747   fprintf(file,"Domain: %g x %g x %g - %g %g %g\n",getP()[0],getP()[1],getP()[2],getQ()[0],getQ()[1],getQ()[2]); 
00748   fprintf(file,"Endc: %d\n",endc->index());
00749   fprintf(file,"NumRawCells: %d\n",numRawCells());
00750   fprintf(file,"NumCells: %d\n",numCells());
00751   fprintf(file,"Elems: %u x %u x %u\n",numElemX(),numElemY(),numElemZ());
00752   fprintf(file,"Steps: %g x %g x %g\n",getDX()[0],getDX()[1],getDX()[2]);
00753   for(;cell!=endc;cell++)
00754   {
00755     Point3D bc;
00756     cell->barycenter(bc);
00757     fprintf(file,"\n\n=====(%d, %d, %d)==============\n",cell->getI(),cell->getJ(),cell->getK());
00758     fprintf(file,"Cell: %u BC <%g, %g,%g>\n",cell->index(),bc[0],bc[1],bc[2]);
00759     fprintf(file,"VERTICE_000: %u\n",cell->vertex_index(VERTEX_000));
00760     fprintf(file,"VERTICE_100: %u\n",cell->vertex_index(VERTEX_100));
00761     fprintf(file,"VERTICE_010: %u\n",cell->vertex_index(VERTEX_010));
00762     fprintf(file,"VERTICE_110: %u\n",cell->vertex_index(VERTEX_110));
00763     fprintf(file,"VERTICE_001: %u\n",cell->vertex_index(VERTEX_001));
00764     fprintf(file,"VERTICE_101: %u\n",cell->vertex_index(VERTEX_101));
00765     fprintf(file,"VERTICE_011: %u\n",cell->vertex_index(VERTEX_011));
00766     fprintf(file,"VERTICE_111: %u\n",cell->vertex_index(VERTEX_111));
00767     fprintf(file,"LEFT_FACE  : %u\n",cell->face_index(LEFT_FACE));
00768     fprintf(file,"RIGHT_FACE : %u\n",cell->face_index(RIGHT_FACE));
00769     fprintf(file,"UP_FACE    : %u\n",cell->face_index(UP_FACE));
00770     fprintf(file,"BOTTOM_FACE: %u\n",cell->face_index(BOTTOM_FACE));
00771     fprintf(file,"FRONT_FACE : %u\n",cell->face_index(FRONT_FACE));
00772     fprintf(file,"BACK_FACE  : %u\n",cell->face_index(BACK_FACE));
00773   }
00774 
00775   fprintf(file,"===================\n   FACES\n=================\n");
00776   //Now iterate the faces
00777   Face_It face = begin_face();
00778   Face_It endf = end_face();
00779 
00780   for(;face!=endf;face++)
00781   {
00782     unsigned c1,c2;
00783     face->getAdjCellIndices(c1,c2);
00784     Point3D bc;
00785     face->barycenter(bc);
00786     fprintf(file,"\n\n=====(%d, %d, %d)==============\n",face->getI(),face->getJ(),face->getK());
00787     fprintf(file,"Face: %u\nBC <%g, %g,%g>\n",face->index(),bc[0],bc[1],bc[2]);
00788     if (c1 == OrthoMesh::INVALID_INDEX)
00789       fprintf(file,"CELL1: none\n");
00790     else
00791       fprintf(file,"CELL1: %u\n",c1);
00792     if (c2 == OrthoMesh::INVALID_INDEX)
00793       fprintf(file,"CELL2: none\n");
00794     else
00795       fprintf(file,"CELL2: %u\n",c2);
00796     fprintf(file,"Normal: ");
00797     switch (face->getNormalOrientation())
00798     {
00799     case NORMAL_X:
00800       fprintf(file,"X\n");
00801       break;
00802     case NORMAL_Y:
00803       fprintf(file,"Y\n");
00804       break;
00805     case NORMAL_Z:
00806       fprintf(file,"Z\n");
00807       break;
00808     default:
00809       abort();
00810     }
00811     fprintf(file,"At Boundary: ");
00812     if (face->at_boundary())
00813       fprintf(file,"true\n");
00814     else
00815       fprintf(file,"false\n");
00816     
00817   }
00818   if (!fileName.empty())
00819     fclose(file);
00820   
00821 }
00822 
00823   
00824 
00825 
00832 void OrthoMesh::setCentralValuesFromFunction(Function3D &f,VecDouble &cValues,unsigned component) 
00833 {
00834   assert(cValues.size() == numCells());
00835          
00836   Cell_It cell = this->begin_cell(),
00837     endc = this->end_cell();
00838   Point<DIM> BC;     
00839   for (;cell!=endc;cell++)
00840   {
00841     BC = cell->barycenter();
00842     cValues(cell->index())=f(BC,component);
00843   }
00844 
00845 }
00846 void OrthoMesh::setCentralValuesFromFunction(Function3D &f,Matrix &MValues) 
00855 {
00856   assert(MValues.m() == numCells());
00857   assert(MValues.n() == f.n_components());
00858          
00859   Cell_It cell = this->begin_cell(),
00860     endc = this->end_cell();
00861   Point<DIM> BC;     
00862   for (;cell!=endc;cell++)
00863   {
00864     BC = cell->barycenter();
00865     for (unsigned cmp=0;cmp<f.n_components();cmp++)
00866       MValues(cell->index(),cmp)=f(BC,cmp);
00867   }
00868 
00869 
00870 }
00871 
00886 bool OrthoMesh::nonZeroInwardNormalComponenentIsPositive(FaceDirection3D face_no) const 
00887 {
00888   //Remember that LEFT_FACE=0,RIGHT_FACE=1,FRONT_FACE=2,BACK_FACE=3,BOTTOM_FACE=4,UP_FACE=5
00889   static bool a[] = {1,0,1,0,1,0};
00890   assert(face_no < 6);
00891   return a[face_no];
00892 
00893   
00894 }
00895 
00896 
00897 
00905 void OrthoMesh::projectCentralValuesAtVertices(const VecDouble &cValues,VecDouble &vValues) const
00906 {
00907   //Check the size of the vectors;
00908   assert(vValues.size() == numVertices());
00909   assert(cValues.size() == numCells());
00910 
00911   //This vector is a counter of the number of cells contributing to each vertice
00912   //It is initialized as zero.
00913   std::vector<char> nCellsPerVertice(numVertices(),0);
00914 
00915   //Initialize the values of the vector vValues.
00916   vValues=0.0;
00917 
00918   Cell_It endc = end_cell();
00919   for (Cell_It cell = begin_cell();cell != endc;cell++)
00920   {
00921     vValues(cell->vertex_index(VERTEX_000)) += cValues(cell->index());
00922     vValues(cell->vertex_index(VERTEX_001)) += cValues(cell->index());
00923     vValues(cell->vertex_index(VERTEX_010)) += cValues(cell->index());
00924     vValues(cell->vertex_index(VERTEX_011)) += cValues(cell->index());
00925     vValues(cell->vertex_index(VERTEX_100)) += cValues(cell->index());
00926     vValues(cell->vertex_index(VERTEX_101)) += cValues(cell->index());
00927     vValues(cell->vertex_index(VERTEX_110)) += cValues(cell->index());
00928     vValues(cell->vertex_index(VERTEX_111)) += cValues(cell->index());
00929 
00930     nCellsPerVertice[cell->vertex_index(VERTEX_000)]++;
00931     nCellsPerVertice[cell->vertex_index(VERTEX_001)]++;
00932     nCellsPerVertice[cell->vertex_index(VERTEX_010)]++;
00933     nCellsPerVertice[cell->vertex_index(VERTEX_011)]++;
00934     nCellsPerVertice[cell->vertex_index(VERTEX_100)]++;
00935     nCellsPerVertice[cell->vertex_index(VERTEX_101)]++;
00936     nCellsPerVertice[cell->vertex_index(VERTEX_110)]++;
00937     nCellsPerVertice[cell->vertex_index(VERTEX_111)]++;
00938 
00939   }
00940 
00941   //Now in vValues we have the sum of all sunrounding cells
00942   //and in nCellsPerVertice the number of cells sunrounding the vertices.
00943   //Rember that these number of cells sunrounding a vertice can differ
00944   //for vertices at boundary.
00945   for (unsigned i =0;i<vValues.size();i++)
00946       vValues(i)/=nCellsPerVertice[i];
00947 }
00948 
00949 
00950 
00951 
00960 void OrthoMesh::setFacesValuesFromFunction(Function3D &f,VecDouble &fValues,unsigned deg,bool bOnlyPrescribedBoundary) 
00961 {
00962   assert(fValues.size() == numFaces());
00963   OrthoMesh::Face_It face    = begin_face();
00964   OrthoMesh::Face_It endFace = end_face();
00965   Point3D P;
00966   for (;face!=endFace;face++)
00967   {
00968     if (bOnlyPrescribedBoundary && !face->at_boundary()) 
00969       continue;
00970 
00971     face->barycenter(P);
00972     if (!f.isInDomain(P,deg))
00973       continue;
00974     else
00975       fValues(face->index())=f(P,deg);
00976   }
00977 }
00978 
00979 
00986 double OrthoMesh::getIntegralAtCells(const VecDouble &cValues) const
00987 {
00988   assert(cValues.size() == numCells());
00989   double acum=0.0;
00990   for (unsigned i=0;i<cValues.size();i++)
00991   {
00992     acum+=cValues(i);
00993   }
00994   return acum*cellVolume();
00995 }
00996 
00997 
00998 
01006 void OrthoMesh::setFacesValuesFromFunction(Function3D &f,Matrix &MValues,const bool bOnlyAtBoundary) 
01007 {
01008   assert(MValues.m() == numFaces());
01009   assert(MValues.n() == f.n_components());
01010   
01011   Face_It face    = begin_face();
01012   Face_It endFace = end_face();
01013   
01014   for (;face!=endFace;face++)
01015   {
01016     if (bOnlyAtBoundary && !face->at_boundary()) 
01017       continue;
01018 
01019     Point<3> P = face->barycenter();
01020     unsigned faceIndex = face->index();
01021     for (unsigned deg=0;deg<f.n_components();deg++)
01022     {
01023       MValues(faceIndex,deg) = f(P,deg);
01024     }
01025   }
01026 }
01027 
01028 
01029 
01030  
01031 
01032 
01033 void OrthoMesh::putWells(VecWellInfo &wells) 
01034 {
01035   m_wells=wells;
01036   Point3D DX_4 = this->getDX();
01037   Point3D DX = this->getDX();
01038   DX_4[0]/=4.0;
01039   DX_4[1]/=4.0;
01040   DX_4[2]/=4.0;
01041   std::vector<Index> vInnerWellCells;
01042   
01043   DelIndexLst delCellLst,delVertLst,delFaceLst;
01044   for (unsigned i=0;i<wells.size();i++)
01045   {
01046       WellInfo well = wells[i];
01047       well.P+=DX_4;
01048       well.Q-=DX_4;
01049       Raw_Cell_It endc = end_raw_cell();
01050       Point3D p;
01051       for(Raw_Cell_It cell = begin_raw_cell();cell != endc;cell++)
01052       {
01053         cell->barycenter(p);
01054         if (well.isPointInWell(p))
01055           delCellLst.addIndex(cell->index(),wells[i].isInnerCell(p,DX));
01056       }
01057 
01058       Raw_Vertex_It endv = end_raw_vertice();
01059       for(Raw_Vertex_It vert = begin_raw_vertice();vert != endv;vert++)
01060       {
01061         if (well.isPointInWell(vert->getPoint()))
01062             delVertLst.addIndex(vert->index());
01063       }
01064       
01065       Raw_Face_It endf = end_raw_face();
01066       for (Raw_Face_It face = begin_raw_face();face != endf;face++)
01067       {
01068         if (well.isPointInWell(face->barycenter()))
01069             delFaceLst.addIndex(face->index());
01070       }
01071   }
01072 
01073   delCellLst.close();
01074   delVertLst.close();
01075   delFaceLst.close();
01076   
01077   _delCellLst = delCellLst;
01078   _delVertLst = delVertLst;
01079   _delFaceLst = delFaceLst;
01080   
01081   _numFaces   -= _delFaceLst.size();
01082   _numCells   -= _delCellLst.size();
01083   _numVertices-= _delVertLst.size();
01084 
01085   delete _pBeginCellH;
01086   delete _pEndCellH;
01087   _pBeginCellH  = new OrthoCellAccessorWithHoles(*this,0); 
01088    _pEndCellH = new OrthoCellAccessorWithHoles(*this,this->numRawCells()-1);
01089    (*_pEndCellH)++;
01090 
01091    delete _pBeginFaceH;
01092    delete _pEndFaceH;
01093    _pBeginFaceH  = new OrthoFaceAccessorWithHoles(*this,0); 
01094    _pEndFaceH = new OrthoFaceAccessorWithHoles(*this,this->numRawFaces()-1);
01095    (*_pEndFaceH)++;
01096 
01097 
01098    delete _pBeginVertH;
01099    delete _pEndVertH;
01100    _pBeginVertH  = new OrthoVerticeAccessorWithHoles(*this,0); 
01101    _pEndVertH = new OrthoVerticeAccessorWithHoles(*this,this->numVertices()-1);
01102    (*_pEndVertH)++;
01103 
01104   
01105    _HGhostLeftOff  = _numRawCells-delCellLst.getNumInnerGhostCells();
01106    _HGhostRightOff = _HGhostLeftOff   + _nElemYZ;
01107    _HGhostFrontOff = _HGhostRightOff  + _nElemYZ;
01108    _HGhostBackOff  = _HGhostFrontOff  + _nElemXZ;
01109    _HGhostBottomOff= _HGhostBackOff   + _nElemXZ;
01110    _HGhostUpOff    = _HGhostBottomOff + _nElemXY;
01111    _HGhostWellOff=_numCells;
01112 
01113 
01114    _numHGhostCells = 2*(_nElemYZ + _nElemXY + _nElemXZ) + delCellLst.size() - delCellLst.getNumInnerGhostCells();
01115    
01116 }
01117 
01118 
01140 void OrthoMesh::getFacesInYZSlab(double X,VecIndex &fIndices) 
01141 {
01142   OrthoMesh::Cell_It cell = begin_cell();
01143   OrthoMesh::Cell_It endc = end_cell();
01144 
01145   double x_tol = getDX()[0]/4.0;
01146   fIndices.reserve(numElemY()*numElemZ());
01147   
01148   //First find the first face in the slab
01149   OrthoMesh::Face_It face = cell->face(LEFT_FACE);
01150   FaceDirection3D faceDir;
01151   if (NumericMethods::a_equal(face->barycenter()[0],X,x_tol))
01152     faceDir=LEFT_FACE;
01153   else
01154   {
01155     faceDir = RIGHT_FACE;
01156     face=cell->face(RIGHT_FACE);
01157     while (!NumericMethods::a_equal(face->barycenter()[0],X,x_tol))
01158     {
01159       cell = cell->neighbor(RIGHT_CELL);
01160       if (cell != endc)
01161       {
01162         face = cell->face(RIGHT_FACE);
01163       }
01164       else
01165       {
01166         throw new Exception("OrthoMesh::getFacesInYZSlab(): Cant find the slab X= %g in the domain",X);
01167       }
01168     }
01169   }
01170 
01171   //Found the slab
01172   //Now run in Z direction 
01173   do {
01174     //For each Cell in Z, run in Y direction
01175     for(OrthoMesh::Cell_It cellY=cell;cellY != endc;cellY=cellY->neighbor(BACK_CELL))
01176       fIndices.push_back(cellY->face_index(faceDir));
01177   } while(cell->advance(UP_CELL));
01178   
01179   
01180   
01181 }
01182 
01183 
01184 
01185 
01186 
01187 void OrthoMesh::getFacesInSlab(double Const,VecIndex &fIndices,unsigned dim) 
01188 {
01189   OrthoMesh::Cell_It cell = begin_cell();
01190   OrthoMesh::Cell_It endc = end_cell();
01191   static FaceDirection3D OrthoFace[3][2] = { {LEFT_FACE  ,RIGHT_FACE},
01192                                               {FRONT_FACE ,BACK_FACE },
01193                                               {BOTTOM_FACE,UP_FACE   }};
01194   static CellDirection3D CellDir[3] = {RIGHT_CELL,BACK_CELL,UP_CELL};
01195 
01196   static CellDirection3D RunDim[3][2] = { {BACK_CELL  ,UP_CELL},
01197                                           {RIGHT_CELL  ,UP_CELL},
01198                                           {RIGHT_CELL  ,BACK_CELL}};
01199 
01200     
01201   double _tol = getDX()[dim]/4.0;
01202   fIndices.reserve(numElems()[(dim+1)%3]*numElems()[(dim-1)%3]);
01203   
01204   //First find the first face in the slab
01205   FaceDirection3D faceDir=OrthoFace[dim][0];
01206   OrthoMesh::Face_It face = cell->face(faceDir);
01207   if (!NumericMethods::a_equal(face->barycenter()[dim],Const,_tol))
01208   {
01209     faceDir = OrthoFace[dim][1];
01210     face=cell->face(faceDir);
01211     while (!NumericMethods::a_equal(face->barycenter()[dim],Const,_tol))
01212     {
01213       cell = cell->neighbor(CellDir[dim]);
01214       if (cell != endc)
01215       {
01216         face = cell->face(faceDir);
01217       }
01218       else
01219       {
01220         throw new Exception("OrthoMesh::getFacesInYZSlab(): Cant find the slab Dir %d = %g in the domain",dim,Const);
01221       }
01222     }
01223   }
01224 
01225   do {
01226     for(OrthoMesh::Cell_It cell2=cell;cell2 != endc;cell2=cell2->neighbor(RunDim[dim][0]))
01227       fIndices.push_back(cell2->face_index(faceDir));
01228   } while(cell->advance(RunDim[dim][1]));
01229 
01230 }
01231 
01232 
01233 
01234 
01235 
01236 
01237 
01238 
01239 
01240 
01241 
01247 void OrthoMesh::getCellsIndicesInDomain(Function3D &f,VecIndex &vec) 
01248 {
01249   vec.clear();
01250   Cell_It cell = this->begin_cell(),
01251     endc = this->end_cell();
01252   for (;cell!=endc;cell++)
01253   {
01254     Point3D barycenter = cell->barycenter(); 
01255     if (f.isInDomain(barycenter))
01256     {
01257       vec.push_back(cell->index());
01258     }
01259     
01260   }
01261 }
01262 
01267 OrthoMesh::Cell_It OrthoMesh::getCellAt(Point3D &p) 
01268 {
01269   Cell_It cell = this->begin_cell(),
01270     endc = this->end_cell();
01271   for (;cell!=endc;cell++)
01272   {
01273     if (cell->containPoint(p)) 
01274       return cell;
01275   }
01276   return endc;
01277 }
01278 
01279 
01280 
01287 void OrthoMesh::tagFacesInDomain(Function3D &f,VecTag &tags,char value)
01288 {
01289   assert(tags.size() == numFaces());
01290   OrthoMesh::Face_It face    = begin_face();
01291   OrthoMesh::Face_It endFace = end_face();
01292 
01293 
01294   
01295   for (;face!=endFace;face++)
01296   {
01297     Point<3> P = face->barycenter();
01298     if (f.isInDomain(P,0))
01299       tags[face->index()] = value;
01300   }
01301 }
01302 
01303 void OrthoMesh::getFacesIndicesInDomain(Function3D &f, VecIndex &vec) 
01304 {
01305   OrthoMesh::Face_It face    = begin_face();
01306   OrthoMesh::Face_It endFace = end_face();
01307   vec.clear();
01308 
01309   for (;face!=endFace;face++)
01310   {
01311     Point<3> P = face->barycenter();
01312     if (f.isInDomain(P,0))
01313       vec.push_back(face->index());
01314   }
01315 }
01316 
01317 
01318 
01319 void OrthoMesh::printCells(VecIndex &vec,std::ostream &out) 
01320 {
01321   for (unsigned i = 0;i<vec.size();i++)
01322   {
01323     out <<  i << ")";
01324     printCell(get_cell(vec[i]),out);
01325   }
01326 }
01327 
01328 void OrthoMesh::printCell(Cell_It cell,std::ostream &out) 
01329 {
01330   char str[500];
01331   sprintf(str,"==== Cell %d ====\n",cell->index());
01332   out << str;
01333   out << "Barycenter: " << cell->barycenter() << std::endl;
01334   //     cout << "0 = " <<  cell->vertex(0) << endl;
01335   //     cout << "1 = " <<  cell->vertex(1) << endl;
01336   //     cout << "2 = " <<  cell->vertex(2) << endl;
01337   //     cout << "3 = " <<  cell->vertex(3) << endl;
01338   //     cout << "4 = " <<  cell->vertex(4) << endl;
01339   //     cout << "5 = " <<  cell->vertex(5) << endl;
01340   //     cout << "6 = " <<  cell->vertex(6) << endl;
01341   //     cout << "7 = " <<  cell->vertex(7) << endl;
01342 
01343   out << "V000 = "  << cell->vertex_index(VERTEX_000) <<  ") " << cell->vertex(VERTEX_000)  << std::endl;
01344   out << "V100 = " << cell->vertex_index(VERTEX_100) <<  ") " << cell->vertex(VERTEX_100) << std::endl;
01345   out << "V110 = " << cell->vertex_index(VERTEX_110) <<  ") " << cell->vertex(VERTEX_110) << std::endl;
01346   out << "V010 = " << cell->vertex_index(VERTEX_010) <<  ") " << cell->vertex(VERTEX_010) << std::endl;
01347   out << "V001 = " << cell->vertex_index(VERTEX_001) <<  ") " << cell->vertex(VERTEX_001) << std::endl;
01348   out << "V101 = " << cell->vertex_index(VERTEX_101) <<  ") " << cell->vertex(VERTEX_101) << std::endl;
01349   out << "V111 = " << cell->vertex_index(VERTEX_111) <<  ") " << cell->vertex(VERTEX_111) << std::endl;
01350   out << "V011 = " << cell->vertex_index(VERTEX_011) <<  ") " << cell->vertex(VERTEX_011) << std::endl;
01351 
01352   out << "RIGHT_FACE  " << cell->face_index (RIGHT_FACE) << ")" << std::endl;
01353   out << "LEFT_FACE   " << cell->face_index  (LEFT_FACE) << ")"  << std::endl;
01354   out << "FRONT_FACE  " << cell->face_index (FRONT_FACE) << ")" << std::endl;
01355   out << "BACK_FACE   " << cell->face_index  (BACK_FACE) << ")"  << std::endl;
01356   out << "BOTTOM_FACE " << cell->face_index(BOTTOM_FACE) << ")"<< std::endl;
01357   out << "UP_FACE     " << cell->face_index    (UP_FACE) << ")"    << std::endl;
01358 
01359 
01360 }
01361 
01362 
01363 
01364 
01365 void OrthoMesh::buildCashInfo() 
01366 {
01367   
01368   m_vFCash.clear();
01369   m_vFCash.reserve(numFaces());
01370   Face_It endf = end_face();
01371   for (Face_It face = begin_face();face!=endf;face++)
01372   {
01373     struct FaceInfoCash info;
01374     face->getAdjCellIndices(info.c1,info.c2);
01375     info.normal=face->getNormalOrientation();
01376     m_vFCash.push_back(info);
01377   }
01378 }
01379 
01380 
01381 
01382 
01383 
01384 
01385 
01386 
01387 
01390 OrthoMesh::Vertex_It OrthoMesh::getNearestVertice(Point3D p) 
01391 {
01392   assert(NumericMethods::isInCube(p,getP(),getQ()));
01393   p-=getP();
01394   unsigned i = round(p[0]/getDX()[0]);
01395   unsigned j = round(p[1]/getDX()[1]);
01396   unsigned k = round(p[2]/getDX()[2]);
01397   return OrthoVerticeAccessorWithHoles(*this,this->getVerticeIndexFromPosition(i,j,k));
01398 }
01399 
01411 double OrthoMesh::getOutwardNormalOrientation(FaceDirection3D faceDir) 
01412 {
01413   static double tbl[] = {-1,1,-1,1,-1,1}; //THE TABLE OF CONVERSIONS
01414   return tbl[faceDir];
01415 }
01416 
01417 
01418 
01419 
01420 void OrthoMesh::projectVelocitiesInFacesToCells(Matrix &vel,const VecDouble &vNormalFaceVel) 
01421 {
01422   assert(vel.m() == numCells() && vel.n() == 3);
01423   assert(vNormalFaceVel.size() == numFaces());
01424   OrthoMesh::Cell_It endc = end_cell();
01425 
01426   for (OrthoMesh::Cell_It cell = begin_cell();cell!=endc;cell++)
01427   {
01428     int i = cell->index();
01429     vel(i,0) = (vNormalFaceVel(cell->face_index(LEFT_FACE))    + vNormalFaceVel(cell->face_index(RIGHT_FACE)))/2.0;
01430     vel(i,1) = (vNormalFaceVel(cell->face_index(FRONT_FACE))   + vNormalFaceVel(cell->face_index(BACK_FACE)))/2.0;
01431     vel(i,2) = (vNormalFaceVel(cell->face_index(BOTTOM_FACE))  + vNormalFaceVel(cell->face_index(UP_FACE)))/2.0;
01432   }
01433 
01434 }
01435 
01436 
01437 
01438 
01439 
01440 
01441 
01442 
01443 
01444 
01445 void OrthoMesh::projectVelocitiesInFacesToVertices(Matrix &Vel,const VecDouble &vNormalFaceVel) 
01446 {
01447   assert(Vel.m() == numVertices());
01448   assert(Vel.n() == 3);
01449   assert(vNormalFaceVel.size() == numFaces());
01450   char M[numVertices()][3];
01451   memset(M,0,3*numVertices());
01452   OrthoMesh::Face_It endf = end_face();
01453   Vel=0.0;
01454   for (OrthoMesh::Face_It face = begin_face();face!=endf;face++)
01455   {
01456     Vel(face->vertex_index_00(),face->getNormalOrientation())+=vNormalFaceVel(face->index());
01457     Vel(face->vertex_index_01(),face->getNormalOrientation())+=vNormalFaceVel(face->index());
01458     Vel(face->vertex_index_10(),face->getNormalOrientation())+=vNormalFaceVel(face->index());
01459     Vel(face->vertex_index_11(),face->getNormalOrientation())+=vNormalFaceVel(face->index());
01460 
01461     M[face->vertex_index_00()][face->getNormalOrientation()]+=1;
01462     M[face->vertex_index_01()][face->getNormalOrientation()]+=1;
01463     M[face->vertex_index_10()][face->getNormalOrientation()]+=1;
01464     M[face->vertex_index_11()][face->getNormalOrientation()]+=1;
01465   }
01466   for (unsigned i=0;i<numVertices();i++)
01467   {
01468     Vel(i,0)/=M[i][0];
01469     Vel(i,1)/=M[i][1];
01470     Vel(i,2)/=M[i][2];
01471   }
01472 }
01473 
01474 void OrthoMesh::getMaxAbsVelocitiesInFacesByComponents(VecDouble &vel, const VecDouble &vNormalFaceVel) 
01475 {
01476   assert(vel.size() == 3);
01477   assert(vNormalFaceVel.size() == numFaces());
01478   OrthoMesh::Face_It endf = end_face();
01479   vel=0.0;
01480   for (OrthoMesh::Face_It face = begin_face();face!=endf;face++)
01481   {
01482     unsigned normal = face->getNormalOrientation();
01483     double dd = fabs(vNormalFaceVel(face->index()));
01484     if ( vel(normal) < dd)
01485       vel(normal)=dd;
01486   }
01487 }
01488 
01489 
01490 
01491 
01492 
01493 
01494 
01495 
01496 
01497 
01498 
01499 
01500 double OrthoMesh::getIntegralAtVertices(const VecDouble &vValues) const 
01501 {
01502   double acum=0.0;
01503   OrthoMesh::Cell_It cell = begin_cell();
01504   OrthoMesh::Cell_It endc = end_cell();
01505 
01506   for (  ;cell!=endc;cell++)
01507   {
01508   
01509     acum+= (vValues(cell->vertex_index(VERTEX_000)) +
01510             vValues(cell->vertex_index(VERTEX_001)) +
01511             vValues(cell->vertex_index(VERTEX_010)) +
01512             vValues(cell->vertex_index(VERTEX_011)) +
01513             vValues(cell->vertex_index(VERTEX_100)) +
01514             vValues(cell->vertex_index(VERTEX_101)) +
01515             vValues(cell->vertex_index(VERTEX_110)) +
01516             vValues(cell->vertex_index(VERTEX_111)))/8.0; 
01517     
01518   }
01519   return acum*cellVolume();
01520 }