CO2INJECTION: orthofaceaccessor.cpp Source File

00001 #include "orthomesh.h"
00002 
00003  OrthoFaceAccessor::OrthoFaceAccessor(const OrthoMesh &mesh,unsigned index) 
00004    :m_mesh(&mesh)
00005 {
00006    m_index = index;
00007    m_mesh->getFacePositionFromIndex(index,m_i,m_j,m_k,m_normal);
00008 }
00009 
00010 OrthoFaceAccessor::OrthoFaceAccessor(const OrthoCellAccessor &cell,FaceDirection3D dirEnum) 
00011    :m_mesh(&(cell.getMesh()))
00012 {
00013   m_mesh->getFaceInfoFromCell(dirEnum,cell.getI(),cell.getJ(),cell.getK(),&m_i,&m_j,&m_k,&m_index,&m_normal);
00014    
00015 }
00016 
00017 
00018 
00019 bool OrthoFaceAccessor::at_boundary() 
00020 {
00021   assert(this->isValid());
00022   if (m_normal == OrthoMesh::NORMAL_X)
00023   {
00024     return m_i == 0 || m_i == m_mesh->numElemX();
00025   }
00026   else if (m_normal == OrthoMesh::NORMAL_Y)
00027   {
00028     return m_j == 0 || m_j == m_mesh->numElemY();
00029   }
00030   else //m_normal == OrthoMesh::NORMAL_Z
00031   {
00032     return m_k == 0 || m_k == m_mesh->numElemZ();
00033   }
00034   
00035 }
00036 
00037 
00038 
00039 void OrthoFaceAccessor::barycenter(Point3D &p) 
00040 {
00041   assert(this->isValid());
00042   switch (m_normal)
00043   {
00044   case OrthoMesh::NORMAL_X:
00045     p[0]=m_mesh->_fxP0[0] + m_i*m_mesh->_DX[0];
00046     p[1]=m_mesh->_fxP0[1] + m_j*m_mesh->_DX[1];
00047     p[2]=m_mesh->_fxP0[2] + m_k*m_mesh->_DX[2];
00048     break;
00049   case OrthoMesh::NORMAL_Y:
00050     p[0]=m_mesh->_fyP0[0] + m_i*m_mesh->_DX[0];
00051     p[1]=m_mesh->_fyP0[1] + m_j*m_mesh->_DX[1];
00052     p[2]=m_mesh->_fyP0[2] + m_k*m_mesh->_DX[2];
00053     break;
00054   case OrthoMesh::NORMAL_Z:
00055     p[0]=m_mesh->_fzP0[0] + m_i*m_mesh->_DX[0];
00056     p[1]=m_mesh->_fzP0[1] + m_j*m_mesh->_DX[1];
00057     p[2]=m_mesh->_fzP0[2] + m_k*m_mesh->_DX[2];
00058     break;
00059   }
00060 }
00061 
00062 Point3D OrthoFaceAccessor::barycenter()
00063 {
00064   Point3D p;
00065   barycenter(p);
00066   return p;
00067 }
00068   
00069 
00070 void OrthoFaceAccessor::getAdjCellIndices(unsigned &index1,unsigned &index2)  const
00071 {
00072   assert(this->isValid());
00073 
00074   if (m_normal == OrthoMesh::NORMAL_X)
00075   {
00076     if (m_i == m_mesh->_nElemX)
00077     {
00078       index2 = OrthoMesh::INVALID_INDEX;
00079       index1 = m_index - m_j - m_k*m_mesh->_nElemY -1;
00080     }
00081     else if (m_i == 0)
00082     {
00083       index1 = OrthoMesh::INVALID_INDEX;
00084       index2 = m_index - m_j - m_k*m_mesh->_nElemY;
00085     }
00086     else
00087     {
00088       index2 = m_index - m_j - m_k*m_mesh->_nElemY;
00089       index1 = index2 -1;
00090     }
00091   }
00092   if (m_normal == OrthoMesh::NORMAL_Y)
00093   {
00094     unsigned setIndex = m_index - m_mesh->_fyOff;
00095     
00096     if (m_j == 0)
00097     {
00098       index2 = setIndex - m_k*m_mesh->_nElemX;
00099       index1 = OrthoMesh::INVALID_INDEX;
00100     }
00101     else if (m_j == m_mesh->_nElemY)
00102     {
00103       index1 = setIndex - (m_k+1)*m_mesh->_nElemX;
00104       index2 = OrthoMesh::INVALID_INDEX;
00105     }
00106     else
00107     {
00108       index2 = setIndex - m_k*m_mesh->_nElemX;
00109       index1 = index2 - m_mesh->_nElemX;
00110     }
00111   }
00112   if (m_normal == OrthoMesh::NORMAL_Z)
00113   {
00114     unsigned setIndex = m_index - m_mesh->_fzOff;
00115     
00116     if (m_k == 0)
00117     {
00118       index2 = setIndex;
00119       index1 = OrthoMesh::INVALID_INDEX;
00120     }
00121     else if (m_k == m_mesh->_nElemZ)
00122     {
00123       index2 = OrthoMesh::INVALID_INDEX;
00124       index1 = setIndex - m_mesh->_ZStride;
00125     }
00126     else
00127     {
00128       index2 = setIndex;
00129       index1 = index2 - m_mesh->_ZStride;
00130     }
00131   }
00132 }
00133 
00134 
00135 void OrthoFaceAccessor::ghostAdjCellIndices(unsigned &index1,unsigned &index2) 
00136 {
00137   assert(this->isValid());
00138 
00139   if (m_normal == OrthoMesh::NORMAL_X)
00140   {
00141     if (m_i == m_mesh->_nElemX)
00142     {
00143       index2 = m_mesh->_GhostRightOff +  m_j + m_k*m_mesh->_nElemY;
00144       index1 = m_index - m_j - m_k*m_mesh->_nElemY -1;
00145     }
00146     else if (m_i == 0)
00147     {
00148       index1 = m_mesh->_GhostLeftOff + m_j + m_k*m_mesh->_nElemY;
00149       index2 = m_index - m_j - m_k*m_mesh->_nElemY;
00150     }
00151     else
00152     {
00153       index2 = m_index - m_j - m_k*m_mesh->_nElemY;
00154       index1 = index2 -1;
00155     }
00156   }
00157   if (m_normal == OrthoMesh::NORMAL_Y)
00158   {
00159     unsigned setIndex = m_index - m_mesh->_fyOff;
00160     
00161     if (m_j == 0)
00162     {
00163       index2 = setIndex - m_k*m_mesh->_nElemX;
00164       index1 = m_mesh->_GhostFrontOff + m_i + m_k*m_mesh->_nElemX;
00165     }
00166     else if (m_j == m_mesh->_nElemY)
00167     {
00168       index1 = setIndex - (m_k+1)*m_mesh->_nElemX;
00169       index2 = m_mesh->_GhostBackOff + m_i + m_k*m_mesh->_nElemX;
00170     }
00171     else
00172     {
00173       index2 = setIndex - m_k*m_mesh->_nElemX;
00174       index1 = index2 - m_mesh->_nElemX;
00175     }
00176   }
00177   if (m_normal == OrthoMesh::NORMAL_Z)
00178   {
00179     unsigned setIndex = m_index - m_mesh->_fzOff;
00180     
00181     if (m_k == 0)
00182     {
00183       index2 = setIndex;
00184       index1 = m_mesh->_GhostBottomOff + m_i + m_j*m_mesh->_nElemX;
00185     }
00186     else if (m_k == m_mesh->_nElemZ)
00187     {
00188       index2 = m_mesh->_GhostUpOff + m_i + m_j*m_mesh->_nElemX;
00189       index1 = setIndex - m_mesh->_ZStride;
00190     }
00191     else
00192     {
00193       index2 = setIndex;
00194       index1 = index2 - m_mesh->_ZStride;
00195     }
00196   }
00197   
00198 }
00199 
00200 
00201 bool OrthoFaceAccessor::isValid()  const
00202 {
00203   if (m_mesh == NULL)
00204   {
00205     assert(0);
00206     return false;
00207 
00208   }
00209   if (m_index >= m_mesh->numRawFaces())
00210   {
00211     assert(0);
00212 
00213     return false;
00214   }
00215   if (m_normal == OrthoMesh::NORMAL_X)
00216   {
00217     if (m_i <= m_mesh->_nElemX && m_j < m_mesh->_nElemY && m_k < m_mesh->_nElemZ) 
00218       return true;
00219     else
00220     {
00221       assert(0);
00222       return false;
00223     }
00224   }
00225   else if (m_normal == OrthoMesh::NORMAL_Y)
00226   {
00227     if (m_i < m_mesh->_nElemX && m_j <= m_mesh->_nElemY && m_k < m_mesh->_nElemZ) 
00228       return true;
00229     else
00230     {
00231       assert(0);
00232       return false;
00233     }
00234   }
00235   else //m_normal == NORMAL_Z
00236   {
00237     if (m_i < m_mesh->_nElemX && m_j < m_mesh->_nElemY && m_k <= m_mesh->_nElemZ) 
00238       return true;
00239     else
00240     {
00241       assert(0);
00242       return false;
00243     }
00244   }
00245   
00246   
00247 }
00248 
00249 
00250 double OrthoFaceAccessor::normal_multiply(double *v) const 
00251 {
00252   assert(this->isValid());
00253   return v[m_normal];
00254 }
00255 
00256 
00257 
00258 
00259 void OrthoFaceAccessor::operator++ (int) 
00260 {
00261   assert(this->isValid());
00262   m_index++;
00263   if (m_normal == OrthoMesh::NORMAL_X)
00264   {
00265     if (m_index != m_mesh->_fyOff)
00266       m_mesh->advanceIJK(m_i,m_j,m_k,
00267                         0,m_mesh->_nElemX,
00268                         0,m_mesh->_lastJ);
00269     else
00270     {
00271       m_i = m_j = m_k =0;
00272       m_normal = OrthoMesh::NORMAL_Y;
00273     }
00274   }
00275   else if (m_normal == OrthoMesh::NORMAL_Y)
00276   {
00277     if (m_index != m_mesh->_fzOff)
00278     {
00279       m_mesh->advanceIJK(m_i,m_j,m_k,
00280                         0,m_mesh->_lastI,
00281                         0,m_mesh->_nElemY);
00282     }
00283     else
00284     {
00285       m_i = m_j = m_k = 0;
00286       m_normal = OrthoMesh::NORMAL_Z;
00287     }
00288   }
00289   else //m_normal == NORMAL_Z
00290   {
00291     if (m_index != m_mesh->numRawFaces())
00292     {
00293       m_mesh->advanceIJK(m_i,m_j,m_k,
00294                         0,m_mesh->_lastI,
00295                         0,m_mesh->_lastJ);
00296     }
00297     //else we reach the end, the iterator is in a state equal to the sentinel. 
00298     //hope the user compares the iterator with the sentinel because if the user increments 
00299     //again the iterator, the program will abort. 
00300   }
00301 
00302 }
00303 /*
00304 void OrthoFaceAccessor::advance_inner() 
00305 {
00306   assert(isValid());
00307   assert(!at_boundary());
00308   
00309   
00310   if (m_i < m_mesh->_LastI) 
00311   {
00312     m_i++;
00313     m_index++;
00314   }
00315     
00316   else
00317   {
00318     if (m_j < m_mesh->m_innerFacesLastIJK[m_normal][1])
00319     {
00320       m_i=m_mesh->m_innerReset[m_normal][0];
00321       m_j++;
00322       m_index+=3;
00323     }
00324     else
00325     {
00326       if (m_k < m_mesh->m_innerFacesLastIJK[m_normal][2])
00327       {
00328         m_i=m_mesh->m_innerReset[m_normal][0];
00329         m_j=m_mesh->m_innerReset[m_normal][1];
00330         m_k++;
00331         m_index+=m_mesh->m_innerFacesChangeKOffset[m_normal];  
00332       }
00333       else
00334       {
00335         m_index=m_mesh->m_innerFacesChangeNormalIndex[m_normal];
00336         inc_normal();
00337         m_i=m_mesh->m_innerReset[m_normal][0];
00338         m_j=m_mesh->m_innerReset[m_normal][1];
00339         m_k=m_mesh->m_innerReset[m_normal][2];
00340       }
00341     }
00342   }
00343     }    
00344    
00345 
00346 
00347 
00348 */
00349   
00350     
00351 
00352 void OrthoFaceAccessor::advance_inner() 
00353 {
00354   assert(isValid());
00355   assert(!at_boundary());
00356   
00357   if (m_normal == OrthoMesh::NORMAL_X)
00358   {
00359     if (m_i < m_mesh->_LastI) 
00360       {
00361         m_i++;
00362         m_index++;
00363       }
00364       else
00365       {
00366         if (m_j < m_mesh->_LastJ)
00367         {
00368           m_i=1;
00369           m_j++;
00370           m_index+=3;
00371         }
00372         else
00373         {
00374           if (m_k < m_mesh->_LastK)
00375           {
00376             m_i=1;
00377             m_j=0;
00378             m_k++;
00379             m_index+=3;  
00380           }
00381           else
00382           {
00383             m_index=m_mesh->m_innerFacesChangeNormalIndex[OrthoMesh::NORMAL_X];
00384             m_normal=OrthoMesh::NORMAL_Y;
00385             m_i=0;
00386             m_j=1;
00387             m_k=0;
00388           }  
00389         }
00390       }
00391   }
00392   else if (m_normal == OrthoMesh::NORMAL_Y)
00393   {
00394     if (m_i < m_mesh->_LastI) 
00395       {
00396         m_i++;
00397         m_index++;
00398       }
00399       else
00400       {
00401         if (m_j < m_mesh->_LastJ)
00402         {
00403           m_i=0;
00404           m_j++;
00405           m_index++;
00406         }
00407         else
00408         {
00409           if (m_k < m_mesh->_LastK)
00410           {
00411             m_i=0;
00412             m_j=1;
00413             m_k++;
00414             m_index+=m_mesh->m_innerFacesChangeKOffset[OrthoMesh::NORMAL_Y];  
00415           }
00416           else
00417           {
00418             m_index=m_mesh->m_innerFacesChangeNormalIndex[OrthoMesh::NORMAL_Y];
00419             m_normal=OrthoMesh::NORMAL_Z;
00420             m_i=0;
00421             m_j=0;
00422             m_k=1;
00423           }  
00424         }
00425       }
00426   }
00427   else
00428   {
00429     if (m_i < m_mesh->_LastI) 
00430       {
00431         m_i++;
00432         m_index++;
00433       }
00434       else
00435       {
00436         if (m_j < m_mesh->_LastJ)
00437         {
00438           m_i=0;
00439           m_j++;
00440           m_index++;
00441         }
00442         else
00443         {
00444           if (m_k < m_mesh->_LastK)
00445           {
00446             m_i=0;
00447             m_j=0;
00448             m_k++;
00449             m_index++;  
00450           }
00451           else
00452           {
00453             *this=*(m_mesh->_pEndFace);
00454           }  
00455         }
00456       }
00457   }
00458 }
00459 
00460 
00461 
00462 
00463 
00464 
00465 
00466 
00467 
00468 
00469   
00470 
00471 
00472 
00473   
00474 
00475 
00476  OrthoFaceAccessor::~OrthoFaceAccessor() 
00477 {
00478 
00479 }
00480 
00481 
00482 bool OrthoFaceAccessor::operator!= (const OrthoFaceAccessor& face) const 
00483 {
00484   assert(m_index <= m_mesh->numRawFaces());
00485   assert(face.m_index <= m_mesh->numRawFaces());
00486   return m_index != face.m_index;
00487 }
00488 
00489 
00490 
00491 
00492 double OrthoFaceAccessor::area() 
00493 {
00494   assert(isValid());
00495   return m_mesh->_faceArea[m_normal];
00496 }
00497 
00498 
00499 double OrthoFaceAccessor::areaPerCellVol() 
00500 {
00501   assert(m_index < m_mesh->numRawFaces());
00502   return m_mesh->_faceAreaPerCellVol[m_normal];
00503 }
00504 
00505 
00506 
00507  OrthoFaceAccessor::OrthoFaceAccessor() 
00508 {
00509   m_mesh = NULL;
00510   m_i=m_j=m_k=m_index=OrthoMesh::INVALID_INDEX;
00511   
00512 }
00513 
00514 
00515 bool OrthoFaceAccessor::hasNegCell() const
00516 {
00517   switch  (m_normal)
00518   {
00519   case OrthoMesh::NORMAL_X:
00520     return m_i != 0;
00521   case OrthoMesh::NORMAL_Y:
00522     return m_j !=0;
00523   case OrthoMesh::NORMAL_Z:
00524     return m_k !=0;
00525   default:
00526     assert(0);
00527     return false;
00528   }
00529 }
00530 
00531 
00532 bool OrthoFaceAccessor::hasPosCell() const
00533 {
00534   switch  (m_normal)
00535   {
00536   case OrthoMesh::NORMAL_X:
00537     return m_i != m_mesh->_nElemX;
00538   case OrthoMesh::NORMAL_Y:
00539     return m_j != m_mesh->_nElemY;
00540   case OrthoMesh::NORMAL_Z:
00541     return m_k != m_mesh->_nElemZ;
00542   default:
00543     assert(0);
00544     return false;
00545   }
00546 }
00547 
00548 
00549 
00550 
00551 
00557 FaceDirection3D OrthoFaceAccessor::getFaceDir() 
00558 {
00559   assert(at_boundary());
00560   switch  (m_normal)
00561   {
00562   case OrthoMesh::NORMAL_X:
00563     return (m_i == 0) ? LEFT_FACE : RIGHT_FACE; 
00564   case OrthoMesh::NORMAL_Y:
00565     return (m_j ==0) ? FRONT_FACE : BACK_FACE;
00566   default:
00567     return (m_k == 0 ) ? BOTTOM_FACE : UP_FACE;
00568   }
00569 }
00570 
00571 
00572 
00573 
00574 Index OrthoFaceAccessor::vertex_index_00() 
00575 {
00576   //NORMAL_X
00577   return m_mesh->getVerticeIndexFromPosition(m_i,m_j,m_k);
00578 }
00579 
00580 
00581 Index OrthoFaceAccessor::vertex_index_10() 
00582 {
00583   static unsigned Tbl[3][3] = {{0,1,0},
00584                                {1,0,0},
00585                                {1,0,0}};
00586    
00587    return m_mesh->getVerticeIndexFromPosition(m_i + Tbl[m_normal][0],m_j+Tbl[m_normal][1],m_k+Tbl[m_normal][2]);
00588 }
00589 
00590 
00591 Index OrthoFaceAccessor::vertex_index_01() 
00592 {
00593   static unsigned Tbl[3][3] = {{0,0,1},
00594                                {0,0,1},
00595                                {0,1,0}};
00596    
00597    return m_mesh->getVerticeIndexFromPosition(m_i + Tbl[m_normal][0],m_j+Tbl[m_normal][1],m_k+Tbl[m_normal][2]);
00598 }
00599 
00600 
00601 Index OrthoFaceAccessor::vertex_index_11() 
00602 {
00603   static unsigned Tbl[3][3] = {{0,1,1},
00604                                {1,0,1},
00605                                {1,1,0}};
00606   return m_mesh->getVerticeIndexFromPosition(m_i + Tbl[m_normal][0],m_j+Tbl[m_normal][1],m_k+Tbl[m_normal][2]);
00607 }
00608   
00609 
00610 
00611 
00612 
00616 bool OrthoFaceAccessor::isInner() 
00617 {
00618   unsigned index1,index2;
00619   getAdjCellIndices(index1,index2);
00620   if (index1 == OrthoMesh::INVALID_INDEX || index2 == OrthoMesh::INVALID_INDEX)
00621   {
00622     return false;
00623   }
00624   else
00625   {
00626     OrthoCellAccessor c1(*m_mesh,index1),c2(*m_mesh,index2);
00627     return (!(c1.at_boundary() && c2.at_boundary())); 
00628   }
00629 
00630 }
00631 
00632 
00633 void OrthoFaceAccessor::inc_normal() 
00634 {
00635   static   OrthoMesh::NORMAL_AXIS next_normal[3]={OrthoMesh::NORMAL_Y,OrthoMesh::NORMAL_Z,OrthoMesh::NORMAL_Z};
00636   m_normal=next_normal[m_normal];
00637 }
00638 
00639