CO2INJECTION: mixedhybridbase.cpp Source File

00001 #include "mixedhybridbase.h"
00002 #include "numericmethods.h"
00003 using NumericMethods::_div;
00004 
00005 void MixedHybridBase::buildSparsityPattern(OrthoMesh &mesh,SparsityPattern &sparsePattern) 
00006 {
00007   OrthoMesh::Cell_It cell = mesh.begin_cell();
00008   OrthoMesh::Cell_It endc = mesh.end_cell();
00009   sparsePattern.reinit(mesh.numCells(),mesh.numCells(),7);
00010 
00011   
00012   for(;cell!=endc;cell++)
00013   {
00014     int ci = cell->index();
00015     sparsePattern.add(ci,ci);
00016     //For each boundary face of the cell
00017     for (unsigned int face_no = 0; face_no < GeometryInfo<DIM>::faces_per_cell;++face_no)
00018     {
00019       unsigned adj_index  = cell->neighbor_index(face_no);
00020       if (adj_index != OrthoMesh::INVALID_INDEX)
00021         sparsePattern.add(ci,adj_index);
00022     }
00023   }
00024   sparsePattern.compress();
00025 }
00026 
00027 
00028 
00056 void MixedHybridBase::assemblySystem(SparseMatrix<double> &A, VecDouble &B,OrthoMesh &mesh, HybridFaceBC &bc,const VecDouble &K, const VecDouble &vMobT, const VecDouble &vGravTerm) 
00057 {
00058 
00059 
00060   
00061   //Get the face iterator and the number of faces.
00062   OrthoMesh::Face_It face = mesh.begin_face();
00063   unsigned nFaces = mesh.numFaces();
00064     
00065   //For each face
00066   for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++)
00067   {
00068     unsigned c1,c2;
00069     face->getAdjCellIndices(c1,c2);
00070     
00071     //Make sure the face is not at boundary
00072     if (!face->at_boundary())
00073     {
00074       double K1 = K(c1)*vMobT(c1);
00075       double K2 = K(c2)*vMobT(c2);
00076       double Mh_2 = K1*K2/(K1+K2); //Harmonic media divided by two
00077       double Keff =2*Mh_2 * face->area()*face->areaPerCellVol();
00078 
00079  
00080       
00081       A.add(c1,c1, Keff);
00082       A.add(c1,c2,-Keff);
00083       A.add(c2,c2, Keff);
00084       A.add(c2,c1,-Keff);
00085       if (face->getNormalOrientation() == OrthoMesh::NORMAL_Z)
00086       {
00087         double aux=+Mh_2*(vGravTerm(c1)+vGravTerm(c2))*face->area();
00088         B(c1) += -aux;
00089         B(c2) +=  aux;
00090       }
00091     }
00092   }
00093   //Now process the boundary conditions
00094     
00095     HybridFaceBC::iterator bcEnd = bc.end();
00096     for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
00097     {
00098       OrthoMesh::Face_It face = bcIt->face;
00099       assert(face->at_boundary());
00100       
00101       unsigned c1,c2;
00102       face->getAdjCellIndices(c1,c2);
00103 
00104       if (bcIt->bcType == Dirichlet)
00105       {
00106         double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2;
00107         double Keff = 2*K(c)*vMobT(c)*face->area()*face->areaPerCellVol();
00108         A.add(c,c,Keff);
00109         B(c) += Keff*bcIt->value;
00110       }
00111       else //Neumann condition
00112       {
00113         if (c1 != OrthoMesh::invalidIndex())
00114           B(c1) += -bcIt->value;
00115         else
00116           B(c2) += +bcIt->value;
00117       }
00118     }
00119 }
00120 
00121 
00122 
00123 
00124 void MixedHybridBase::assemblySystem(SparseMatrix<double> &A, VecDouble &B,OrthoMesh &mesh, HybridFaceBC &bc,const VecDouble &K) 
00125 {
00126   //Get the face iterator and the number of faces.
00127   OrthoMesh::Face_It face = mesh.begin_face();
00128   unsigned nFaces = mesh.numFaces();
00129     
00130   //For each face
00131   for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++)
00132   {
00133     unsigned c1,c2;
00134     face->getAdjCellIndices(c1,c2);
00135     
00136     //Make sure the face is not at boundary
00137     if (!face->at_boundary())
00138     {
00139       double K1 = K(c1);
00140       double K2 = K(c2);
00141       double Mh_2 = _div(K1*K2,(K1+K2)); //Harmonic media divided by two
00142       double Keff =2*Mh_2 * face->area()*face->areaPerCellVol();
00143 
00144  
00145       
00146       A.add(c1,c1, Keff);
00147       A.add(c1,c2,-Keff);
00148       A.add(c2,c2, Keff);
00149       A.add(c2,c1,-Keff);
00150     }
00151   }
00152   //Now process the boundary conditions
00153     
00154     HybridFaceBC::iterator bcEnd = bc.end();
00155     for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
00156     {
00157       OrthoMesh::Face_It face = bcIt->face;
00158       assert(face->at_boundary());
00159       
00160       unsigned c1,c2;
00161       face->getAdjCellIndices(c1,c2);
00162 
00163       if (bcIt->bcType == Dirichlet)
00164       {
00165         double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2;
00166         double Keff = 2*K(c)*face->area()*face->areaPerCellVol();
00167         A.add(c,c,Keff);
00168         B(c) += Keff*bcIt->value;
00169       }
00170       else //Neumann condition
00171       {
00172         if (c1 != OrthoMesh::invalidIndex())
00173           B(c1) += -bcIt->value;
00174         else
00175           B(c2) += +bcIt->value;
00176       }
00177     }
00178 
00179 }
00180 
00181 
00185 void MixedHybridBase::assemblySystem(SparseMatrix<double> &A, VecDouble &B,OrthoMesh &mesh, ScalarBC &bc,const VecDouble &vCoeff) 
00186 {
00187   //Get the face iterator and the number of faces.
00188   OrthoMesh::Face_It face = mesh.begin_face();
00189   unsigned nFaces = mesh.numFaces();
00190     
00191   //For each face
00192   for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++)
00193   {
00194     unsigned c1,c2;
00195     face->getAdjCellIndices(c1,c2);
00196     
00197     //Make sure the face is not at boundary
00198     if (!face->at_boundary())
00199     {
00200       double K1 = vCoeff(c1);
00201       double K2 = vCoeff(c2);
00202       double Mh_2 = NumericMethods::_div(K1*K2,(K1+K2)); //Harmonic media divided by two
00203       double Keff =2*Mh_2 * face->area()*face->areaPerCellVol();
00204 
00205  
00206       
00207       A.add(c1,c1, Keff);
00208       A.add(c1,c2,-Keff);
00209       A.add(c2,c2, Keff);
00210       A.add(c2,c1,-Keff);
00211     }
00212   }
00213   //Now process the boundary conditions
00214   ScalarBC::iterator bcEnd = bc.end();
00215   for(ScalarBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
00216   {
00217     OrthoMesh::Face_It face = bcIt->face;
00218     assert(face->at_boundary());
00219       
00220     unsigned c1,c2;
00221     face->getAdjCellIndices(c1,c2);
00222 
00223     double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2;
00224     double Keff = 2*vCoeff(c)*face->area()*face->areaPerCellVol();
00225     A.add(c,c,Keff);
00226     B(c) += Keff*bcIt->value;
00227   }
00228 
00229 }
00230 
00231 
00241 void MixedHybridBase::addPcTerm(VecDouble &B,OrthoMesh &mesh, ScalarBC &bc,const VecDouble &KPc, const VecDouble &vPc) 
00242 {
00243   //Get the face iterator and the number of faces.
00244   OrthoMesh::Face_It face = mesh.begin_face();
00245   unsigned nFaces = mesh.numFaces();
00246     
00247   //For each face
00248   for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++)
00249   {
00250     unsigned c1,c2;
00251     face->getAdjCellIndices(c1,c2);
00252     
00253     //Make sure the face is not at boundary
00254     if (!face->at_boundary())
00255     {
00256       double K1 = KPc(c1);
00257       double K2 = KPc(c2);
00258       double Mh_2 = _div(K1*K2,(K1+K2)); //Harmonic media divided by two
00259       double Keff =2*Mh_2 * face->area()*face->areaPerCellVol();
00260       double q = -Keff*(vPc(c2)-vPc(c1));
00261       B(c1) += -q;
00262       B(c2) += +q;
00263     }
00264   }
00265 
00266   ScalarBC::iterator bcEnd = bc.end();
00267   for(ScalarBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
00268     {
00269       
00270       OrthoMesh::Face_It face = bcIt->face;
00271       assert(face->at_boundary());
00272       
00273       unsigned c1,c2;
00274       face->getAdjCellIndices(c1,c2);
00275 
00276       double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2;
00277       double Keff = 2*KPc(c)*face->area()*face->areaPerCellVol();
00278       B(c) += Keff*(bcIt->value-vPc(c));
00279     }
00280 
00281   
00282 }
00283 
00284 
00285 void MixedHybridBase::getNormalVelocitiesFromPressure(VecDouble &Vq,OrthoMesh &mesh,HybridFaceBC &bc,const VecDouble &K,const VecDouble &vP,const VecDouble &vMobT, const VecDouble &vGravTerm) 
00286 {
00287   assert(Vq.size() == mesh.numFaces());
00288   assert(vP.size() == mesh.numCells());
00289   assert(vMobT.size() == mesh.numCells());
00290   assert(vGravTerm.size() == mesh.numCells());
00291   assert(K.size() == mesh.numCells());
00292   
00293   
00294   //Get the face iterator and the number of faces.
00295   OrthoMesh::Face_It face = mesh.begin_face();
00296   OrthoMesh::Face_It endf = mesh.end_face();
00297 
00298     
00299   //For each face
00300   for (;face!=endf;face++)
00301   {
00302     if (!face->at_boundary())
00303     {
00304       unsigned c1,c2;
00305       face->getAdjCellIndices(c1,c2);
00306       double K1 = K(c1)*vMobT(c1);
00307       double K2 = K(c2)*vMobT(c2);
00308       double Keff = 2*K1*K2/(K1+K2)*face->areaPerCellVol();
00309       Vq(face->index()) = - Keff * (vP(c2) - vP(c1));
00310       if (face->getNormalOrientation() == OrthoMesh::NORMAL_Z)
00311       {
00312         Vq(face->index()) += +K1*K2/(K1+K2)*(vGravTerm(c1)+vGravTerm(c2));
00313       }
00314       
00315     }
00316   }
00317 
00318   //Now for boundary conditions
00319     
00320   HybridFaceBC::iterator bcEnd = bc.end();
00321   for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
00322   {
00323     OrthoMesh::Face_It face = bcIt->face;
00324     assert(face->at_boundary());
00325       
00326     unsigned c1,c2;
00327     face->getAdjCellIndices(c1,c2);
00328 
00329     if (bcIt->bcType == Dirichlet)
00330     {
00331       if (c1 != OrthoMesh::invalidIndex())
00332       {
00333         double Keff = 2*K(c1)*vMobT(c1)*face->areaPerCellVol();
00334         Vq(face->index()) = - Keff*(bcIt->value - vP(c1)); 
00335         //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
00336       }
00337       else
00338       {
00339         double Keff = 2*K(c2)*vMobT(c2)*face->areaPerCellVol();
00340         Vq(face->index()) = Keff*(bcIt->value - vP(c2));
00341         //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
00342       }
00343     }
00344     else
00345     {
00346       if (c1 != OrthoMesh::invalidIndex())
00347         Vq(face->index()) = bcIt->value/face->area();
00348       else
00349         Vq(face->index()) = bcIt->value/face->area();
00350     }
00351   }
00352 }
00353 
00354 
00355 
00356 void MixedHybridBase::getNormalVelocitiesFromPressure(VecDouble &Vq,OrthoMesh &mesh,HybridFaceBC &bc,const VecDouble &K,const VecDouble &vP) 
00357 {
00358   assert(Vq.size() == mesh.numFaces());
00359   assert(vP.size() == mesh.numCells());
00360   assert(K.size() == mesh.numCells());
00361   
00362   
00363   //Get the face iterator and the number of faces.
00364   OrthoMesh::Face_It face = mesh.begin_face();
00365   OrthoMesh::Face_It endf = mesh.end_face();
00366 
00367     
00368   //For each face
00369   for (;face!=endf;face++)
00370   {
00371     if (!face->at_boundary())
00372     {
00373       unsigned c1,c2;
00374       face->getAdjCellIndices(c1,c2);
00375       double K1 = K(c1);
00376       double K2 = K(c2);
00377       double Keff = 2*_div(K1*K2,(K1+K2))*face->areaPerCellVol();
00378       Vq(face->index()) = - Keff * (vP(c2) - vP(c1));
00379     }
00380   }
00381 
00382   //Now for boundary conditions
00383     
00384   HybridFaceBC::iterator bcEnd = bc.end();
00385   for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
00386   {
00387     OrthoMesh::Face_It face = bcIt->face;
00388     assert(face->at_boundary());
00389       
00390     unsigned c1,c2;
00391     face->getAdjCellIndices(c1,c2);
00392 
00393     if (bcIt->bcType == Dirichlet)
00394     {
00395       if (c1 != OrthoMesh::invalidIndex())
00396       {
00397         double Keff = 2*K(c1)*face->areaPerCellVol();
00398         Vq(face->index()) = - Keff*(bcIt->value - vP(c1)); 
00399         //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
00400       }
00401       else
00402       {
00403         double Keff = 2*K(c2)*face->areaPerCellVol();
00404         Vq(face->index()) = Keff*(bcIt->value - vP(c2));
00405         //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
00406       }
00407     }
00408     else
00409     {
00410       if (c1 != OrthoMesh::invalidIndex())
00411         Vq(face->index()) = bcIt->value/face->area();
00412       else
00413         Vq(face->index()) = bcIt->value/face->area();
00414     }
00415   }
00416 
00417 }
00418 
00419 
00420 
00421 
00422 
00423 
00424 void MixedHybridBase::getNormalVelocitiesFromPressure(VecDouble &Vq,OrthoMesh &mesh,HybridFaceBC &bc,const VecDouble &K,const VecDouble &vP,const VecDouble &vMobT, const VecDouble &vGravTerm, ScalarBC &vScalarBC,const VecDouble &KPc, const VecDouble &vPc) 
00425 {
00426   assert(Vq.size() >= mesh.numFaces());
00427   assert(vP.size() >= mesh.numCells());
00428   assert(vMobT.size() >= mesh.numCells());
00429   assert(vGravTerm.size() >=  mesh.numCells());
00430   assert(KPc.size() >= mesh.numCells());
00431   assert(vPc.size() >= mesh.numCells());
00432   assert(K.size() >= mesh.numCells());
00433   
00434   
00435   //Get the face iterator and the number of faces.
00436   OrthoMesh::Face_It face = mesh.begin_face();
00437   OrthoMesh::Face_It endf = mesh.end_face();
00438 
00439     
00440   //For each face
00441   for (;face!=endf;face++)
00442   {
00443     if (!face->at_boundary())
00444     {
00445       unsigned c1,c2;
00446       face->getAdjCellIndices(c1,c2);
00447       double K1 = K(c1)*vMobT(c1);
00448       double K2 = K(c2)*vMobT(c2);
00449       double Keff = 2*K1*K2/(K1+K2)*face->areaPerCellVol();
00450       Vq(face->index()) = - Keff * (vP(c2) - vP(c1));
00451 
00452       //Now for cappillary pressure.
00453       Keff=2*_div(KPc(c1)*KPc(c2),(KPc(c1)+KPc(c2)))*face->areaPerCellVol();
00454       Vq(face->index()) += -Keff*(vPc(c2)-vPc(c1));
00455       if (face->getNormalOrientation() == OrthoMesh::NORMAL_Z)
00456       {
00457         Vq(face->index()) += +K1*K2/(K1+K2)*(vGravTerm(c1)+vGravTerm(c2));
00458       }      //printf("%g ",-Keff*(vPc(c2)-vPc(c1)));
00459 
00460       
00461     }
00462   }
00463 
00464   //Now for boundary conditions
00465     
00466   HybridFaceBC::iterator bcEnd = bc.end();
00467   for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
00468   {
00469     OrthoMesh::Face_It face = bcIt->face;
00470     assert(face->at_boundary());
00471       
00472     unsigned c1,c2;
00473     face->getAdjCellIndices(c1,c2);
00474 
00475     if (bcIt->bcType == Dirichlet)
00476     {
00477       if (c1 != OrthoMesh::invalidIndex())
00478       {
00479         double Keff = 2*K(c1)*vMobT(c1)*face->areaPerCellVol();
00480         Vq(face->index()) = - Keff*(bcIt->value - vP(c1)); 
00481         //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
00482       }
00483       else
00484       {
00485         double Keff = 2*K(c2)*vMobT(c2)*face->areaPerCellVol();
00486         Vq(face->index()) = Keff*(bcIt->value - vP(c2));
00487         //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
00488       }
00489     }
00490     else
00491     {
00492       if (c1 != OrthoMesh::invalidIndex())
00493         Vq(face->index()) = bcIt->value/face->area();
00494       else
00495         Vq(face->index()) = bcIt->value/face->area();
00496     }
00497   }
00498 
00499 
00500   ScalarBC::iterator ScalarBCEnd = vScalarBC.end();
00501   for(ScalarBC::iterator bcIt = vScalarBC.begin();bcIt!=ScalarBCEnd;bcIt++)
00502   {
00503     OrthoMesh::Face_It face = bcIt->face;
00504     assert(face->at_boundary());
00505       
00506     unsigned c1,c2;
00507     face->getAdjCellIndices(c1,c2);
00508 
00509       if (c1 != OrthoMesh::invalidIndex())
00510       {
00511         double Keff = 2*KPc(c1)*face->areaPerCellVol();
00512         Vq(face->index()) += - Keff*(bcIt->value - vP(c1)); 
00513         //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
00514       }
00515       else
00516       {
00517         double Keff = 2*KPc(c2)*face->areaPerCellVol();
00518         Vq(face->index()) += - Keff*(vPc(c2) - bcIt->value);
00519         //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
00520       }
00521   }
00522   
00523 }
00524 
00525 
00526 
00527 
00528 void MixedHybridBase::setBoundaryCondition(HybridFaceBC &vBC,OrthoMesh &mesh,Function3D &fInwardVel,Function3D &fP,bool bCheckNoBC) 
00529 {
00530   vBC.reserve(mesh.numFacesAtBoundary());
00531   Point3D p;  
00532 
00533   OrthoMesh::Cell_It cell = mesh.begin_cell();
00534   OrthoMesh::Cell_It endc = mesh.end_cell();
00535 
00536   _HybridFaceBC bc;
00537   
00538    //For eahc cell in the Mesh
00539    for(;cell != endc;cell++)
00540    {
00541      //For each boundary face of the cell
00542      for (unsigned int face_no = 0; face_no < GeometryInfo<3>::faces_per_cell;++face_no)
00543      {
00544        OrthoMesh::Face_It face = cell->face((FaceDirection3D) face_no);
00545        if (face->at_boundary())
00546        {
00547          face->barycenter(p);
00548          if (fInwardVel.isInDomain(p)) //Neumman condition
00549          {
00550            //Set boundary condition entry 
00551            bc.face = face; //store the face index
00552            bc.bcType = Neumann; //Store the type of condition (Neumman)
00553            if (mesh.nonZeroInwardNormalComponenentIsPositive((FaceDirection3D) face_no))
00554              bc.value = fInwardVel(p)*face->area();
00555            else
00556              bc.value =- fInwardVel(p)*face->area();
00557            vBC.push_back(bc);
00558          }
00559          else if (fP.isInDomain(p)) //Dirichlet condition
00560          {
00561            //Set boundary condition entry 
00562            bc.face = face;
00563            bc.bcType = Dirichlet;
00564            bc.value= fP(p);
00565            vBC.push_back(bc);
00566          }
00567          else if (bCheckNoBC)
00568            throw new Exception("Error in HybridPressionModule::setBoundaryCondition(), face %d in point <%g,%g,%g> does not have any boundary condition",face->index(),p[0],p[1],p[2]);
00569          
00570        }
00571      }
00572    }
00573 
00574 }
00575 
00576 
00577 
00578 
00579 
00580 
00581 
00582 
00583 
00584