BlockMatrixModule Class Reference

#include <blockmatrixmodule.h>

Inheritance diagram for BlockMatrixModule:
Inheritance graph
[legend]
Collaboration diagram for BlockMatrixModule:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 BlockMatrixModule (OrthoMesh &mesh, unsigned n_fracture_cell, VecDouble &vPcInitial, VecDouble &vPor, VecDouble &k, Function1D &dfpc, Function1D &fpcinv, Function1D &fmobt, Function1D &fmobrW)
void construct_bases (VecDouble &bc_pc, double dt)
double get_source (unsigned blockId, double new_bc_pc)
 ~BlockMatrixModule ()
virtual void iterate (TransportBase &trans)
virtual void printOutput ()

Protected Member Functions

void create_bases_functions (double bc_pc, const VecDouble &previous_sw, const VecDouble &previous_si, double por, double k, VecDouble &newSi, VecDouble &newHomSi)
void buildSparsityPattern ()
void assembleSystem (const VecDouble &Sol_old, const VecDouble &sat_old, double _k, double _por, double _pc_bc, double dt)
void Sol_Si ()
void Sol_Sihom ()

Private Attributes

SparsityPattern sp
SparseMatrix< double > M
Matrix m_MSi
Matrix m_MHomSi
OrthoMesh & _mesh
unsigned _m
unsigned n_fracture_cells
VecDouble _BcPcNew
VecDouble _BcPcOld
VecDouble _SwIntegral
VecDouble _vPcInitial
VecDouble rhs_Si
VecDouble rhs_Sihom
VecDouble Sol
VecDouble Sol_hom
VecDouble Sol_old
VecDouble Sol_hom_old
VecDouble sat_old
VecDouble Solution1
VecDouble Solution2
double _dt
double _dx
double _dy
double _dz
double _Bx
VecDouble _vPor
VecDouble _vK
Function1D & _fdpc
Function1D & _fpcinv
Function1D & _fmobt
Function1D & _fmobrW

Detailed Description

Definition at line 11 of file blockmatrixmodule.h.

Constructor & Destructor Documentation

BlockMatrixModule::BlockMatrixModule ( OrthoMesh &  mesh,
unsigned  n_fracture_cell,
VecDouble &  vPcInitial,
VecDouble &  vPor,
VecDouble &  k,
Function1D &  dfpc,
Function1D &  fpcinv,
Function1D &  fmobt,
Function1D &  fmobrW 
)

Definition at line 524 of file blockmatrixmodule.cpp.

00524                                                                                                                                                                                                                 : _mesh(mesh),n_fracture_cells(n_fracture_cell),_BcPcOld(vPcInitial), _vPor(vPor), _vK(vK) ,_fdpc(fdpc), _fpcinv(fpcinv), _fmobt(fmobt), _fmobrW(fmobrW)
00525 {
00526   _dt=1.0;
00527   _dx = _mesh.get_dx();
00528   _dy = _mesh.get_dy();
00529   _dz = _mesh.get_dz();
00530   //  vPcInitial.print();
00531   // cout<<"dx \t"<<_dx<<_dy<<_dz<<endl;
00532   _m =_mesh.numRawVertices();
00533   //Initialize internal structures
00534   m_MSi.reinit(n_fracture_cells,_m);
00535   m_MHomSi.reinit(n_fracture_cells,_m);
00536   //Initialize the vectors to contain the Pc from the fractures
00537   _SwIntegral.reinit(n_fracture_cells);
00538   _BcPcNew=vPcInitial; // Initial Pc
00539 
00540   _Bx = _mesh.meshVolume(); // volume of the matrix block   
00541   m_MHomSi=0;
00542   for (unsigned i=0; i < n_fracture_cells; i++)
00543   {
00544     m_MSi.fill_row(i,_BcPcOld(i));
00545   }
00546   sp.reinit(2*_m, 2*_m,14);
00547   rhs_Si.reinit(2*_m);
00548   rhs_Sihom.reinit(2*_m);
00549   Sol.reinit(_m);
00550   Sol_hom.reinit(_m);
00551   Sol_old.reinit(_m);
00552   Sol_hom_old.reinit(_m);
00553   sat_old.reinit(_m);
00554   Solution1.reinit(2*_m);
00555   Solution2.reinit(2*_m);
00556   buildSparsityPattern();
00557 
00558 }

BlockMatrixModule::~BlockMatrixModule (  )  [inline]

Definition at line 46 of file blockmatrixmodule.h.

00046 {}

Member Function Documentation

void BlockMatrixModule::assembleSystem ( const VecDouble &  Sol_old,
const VecDouble &  sat_old,
double  _k,
double  _por,
double  _pc_bc,
double  dt 
) [protected]

Definition at line 280 of file blockmatrixmodule.cpp.

00281 {
00282   
00283    OrthoMesh::Raw_Vertex_It vt=_mesh.begin_raw_vertice();
00284    OrthoMesh::Raw_Vertex_It vend=_mesh.end_raw_vertice();
00285    // SparseMatrix<double> M;
00286    // cout<<"Assemble System"<<endl;
00287       M.reinit(sp);
00288 
00289       //     rhs_Si.reinit(2*m);
00290     
00291            
00292       while (vt != vend)
00293       {
00294         
00295             OrthoVerticeAccessor ortho(_mesh,vt->index());
00296             if(ortho.at_boundary())
00297             {
00298               // printf("Boundary \n");
00299               M.add(vt->index(),vt->index(),1);
00300               M.add(_m+vt->index(),_m+vt->index(),1);
00301               // cout<<"index \t"<<vt->index()<<endl;
00302               // if (vt->index() <_m)
00303               //  {
00304                 rhs_Si(vt->index())=_bc_pc; //Eq. 46 boundary condition
00305                 rhs_Sihom(vt->index())=1;  //Eq. 47
00306                 // }
00307               // else
00308               // {
00309                  rhs_Si(_m+vt->index())=0;
00310                  rhs_Sihom(_m+vt->index())=0;
00311                  //       }
00312                 
00313             }
00314             else
00315             {
00316              //  printf("Not boundary \n");
00317              //  cout<<"here\t"<<vt->index()<<endl;
00318               //              Sol_old.print();
00319                M.add(vt->index(),vt->index(), _por /( _fdpc(sat_old(vt->index())) *_dt ) ); // si_c eq. 46
00320                //saturation at the midpoint
00321                // cout<<"sat old \t "<<vt->index()<<'\t'<<sat_old(vt->index())<<endl;
00322               
00323                double L = 0.5*( sat_old( vt->index()) + sat_old(vt->index_left()) );
00324                double R = 0.5*(  sat_old( vt->index()) + sat_old(vt->index_right()) );
00325                double F = 0.5*( sat_old( vt->index()) + sat_old(vt->index_front()) );
00326                double B = 0.5*( sat_old( vt->index()) + sat_old(vt->index_back()) );
00327                double U  = 0.5*(  sat_old( vt->index()) + sat_old(vt->index_up()) );
00328                double Bt = 0.5*( sat_old( vt->index()) + sat_old(vt->index_bottom()) );
00329                double _Dx = _dx * _dx;
00330                double _Dy = _dy * _dy;
00331                double _Dz = _dz * _dz;
00332                //              cout<<"L \t"<<L<<'\t'<<_fmobt(L)-_fmobrW(L)<<'\t'<<endl;
00333               
00334                 M.add(vt->index(),_m+vt->index_left(), - _k *  _fmobrW(L)/_Dx );
00335                 M.add(vt->index(),_m+vt->index_right(),-_k *  _fmobrW(R)/_Dx );
00336                 M.add(vt->index(),_m+vt->index_front(),-_k *  _fmobrW(F)/_Dy );
00337                 M.add(vt->index(),_m+vt->index_back(),- _k *  _fmobrW(B)/_Dy );
00338                 M.add(vt->index(),_m+vt->index_bottom(),- _k * _fmobrW(Bt)/_Dz );
00339                 M.add(vt->index(),_m+vt->index_up(), - _k * _fmobrW(U)/_Dz );
00340                 M.add(vt->index(),_m+vt->index(), _k *(  (_fmobrW(L) + _fmobrW(R))/_Dx + (_fmobrW(F) + _fmobrW(B))/_Dy + (_fmobrW(Bt) + _fmobrW(U))/_Dz ) );
00341                 rhs_Si(vt->index())=_por *Sol_old(vt->index()) /( _fdpc( sat_old(vt->index()) )*_dt);
00342                 //cout<<"old solution" << sat_old(vt->index())<<endl;
00343                 rhs_Sihom(vt->index())=0;
00344                 //      ------------Eq. 47----------------------
00345                 //      Si_c
00346                 //      cout<<"ver \t"<<vt->index()<<'\t'<<_m+vt->index()<<'\t'<<_fmobrW(L) + _fmobrW(R) + _fmobrW(F) + _fmobrW(B) + _fmobrW(Bt) + _fmobrW(U)<<endl;
00347                 M.add(_m+vt->index(),vt->index(), _k *( (_fmobt(L)-_fmobrW(L))/_Dx +  (_fmobt(R)-_fmobrW(R))/_Dx +  (_fmobt(F)-_fmobrW(F))/_Dy + (_fmobt(B)-_fmobrW(B))/_Dy +  (_fmobt(Bt)-_fmobrW(Bt))/_Dz + (_fmobt(U)-_fmobrW(U))/_Dz ) );
00348                 M.add(_m+vt->index(),vt->index_left(),- _k * (_fmobt(L)- _fmobrW(L))/_Dx );
00349                 M.add(_m+vt->index(),vt->index_right(),- _k  * (_fmobt(R)- _fmobrW(R))/_Dx );
00350                 M.add(_m+vt->index(),vt->index_front(), -_k * (_fmobt(F)- _fmobrW(F))/_Dy );
00351                 M.add(_m+vt->index(),vt->index_back(), - _k * (_fmobt(B)- _fmobrW(B))/_Dy );
00352                 M.add(_m+vt->index(),vt->index_bottom(),- _k  * (_fmobt(Bt)- _fmobrW(Bt))/_Dz );
00353                 M.add(_m+vt->index(),vt->index_up(), - _k * (_fmobt(U)- _fmobrW(U))/_Dz );
00354 
00355                 //Si_w
00356                 M.add(_m+vt->index(),_m+vt->index(), _k *( (_fmobt(L)+_fmobt(R))/_Dx + (_fmobt(F)+_fmobt(B))/_Dy + (_fmobt(Bt)+_fmobt(U))/_Dz ));
00357                 M.add(_m+vt->index(),_m+vt->index_left(), -_k *_fmobt(L)/_Dx);
00358                 M.add(_m+vt->index(),_m+vt->index_right(), -_k *_fmobt(R)/_Dx );
00359                 M.add(_m+vt->index(),_m+vt->index_front(), -_k *_fmobt(F)/_Dy);
00360                 M.add(_m+vt->index(),_m+vt->index_back(), -_k *_fmobt(B)/_Dy );
00361                 M.add(_m+vt->index(),_m+vt->index_bottom(), -_k *_fmobt(Bt)/_Dz );
00362                 M.add(_m+vt->index(),_m+vt->index_up(), -_k *_fmobt(U)/_Dz );
00363                 
00364                 
00365         
00366             } 
00367             vt++;
00368        }
00369 
00370       // FILE *fdata;
00371       // fdata= fopen("fdata","w");
00372       //    for (unsigned i=0; i< 2*_m; i++)
00373       //    {
00374       //           for (unsigned j=0; j< 2*_m; j++)
00375       //           {
00376       //             //printf("(%u,%u)%g \n",i,j, M.el(i,j));
00377              
00378       //              fprintf(fdata,"%g ",M.el(i,j));
00379       //              //    if (sp.exists(i,j))
00380       //              //        cout<<'\t'<<i<<','<<j<<'\t'<<M.operator()(i,j);
00381       //           }
00382       //           //      cout<<endl;
00383       //         }
00384       //         cout<<"total number of rows"<<_m<<endl;
00385       
00386       Sol_Si();
00387        Sol_Sihom();
00388        // rhs_Si.print();
00389        // rhs_Sihom.print();
00390            //     // cout<<endl;
00391        //   Solution1.print();
00392       //cout<<endl;
00393          //     Solution2.print();
00394       //    cout<<Solution1(40)<<endl;
00395        //   for (unsigned i=0; i< 2*_m ;i++)
00396   // {
00397   //    cout<<"rhs \t"<<i<<'\t'<<rhs_Sihom(i)<<endl;
00398   //       cout<<i<<'\t'<<Solution1(i)<<endl;
00399   // }            
00400 }

void BlockMatrixModule::buildSparsityPattern (  )  [protected]

Definition at line 30 of file blockmatrixmodule.cpp.

00031 {
00032   // SparsityPattern sp;
00033   // unsigned m=_mesh.numRawVertices();
00034   // sp.reinit(2*m,2*m,14);
00035  
00036   // printf(" Vertices: sparsity %u \n", 2*_m);
00037 
00038   OrthoMesh::Raw_Vertex_It end = _mesh.end_raw_vertice();
00039 
00040   unsigned ix, iy, iz, ver_xy, _ver_xy,ver_tol;
00041   _ver_xy=(_mesh.numElemX()+1)*(_mesh.numElemY()+1); // total # of vertices in xy-plane
00042   ver_xy= (_mesh.numElemX()+1)*(_mesh.numElemY()+1)-(_mesh.numElemX()+1); // number used for front and back face boundaries.
00043   // printf("xy %u \n", ver_xy);
00044   ver_tol = (_mesh.numElemX()+1)*(_mesh.numElemY()+1)*(_mesh.numElemY()+1);  
00045 
00046  
00047   /*Front face boundary*/
00048   OrthoMesh::Raw_Vertex_It vf = _mesh.begin_raw_vertice();
00049 
00050   for (iz=0; iz<_mesh.numElemZ()+1; iz++)
00051   {
00052     //  printf("check vertex %u \n ", vf->index());
00053     for (ix=0; ix<_mesh.numElemX()+1; ix++)
00054     {
00055       sp.add(vf->index(),vf->index());
00056       sp.add(_m+vf->index(),_m+vf->index());
00057       //      printf("check vertex %u \n ", vf->index());
00058       vf++;      
00059     }
00060     
00061     for (unsigned k=0; k < ver_xy; k++)
00062     {
00063           vf++;
00064           //  printf("check vertex %u \n ", vf->index());
00065     }
00066   }
00067  
00068   /*Back face boundary*/
00069   OrthoMesh::Raw_Vertex_It vb = _mesh.begin_raw_vertice();
00070   // printf("here %u \n", vb->index());
00071  
00072   for (iz=0; iz<_mesh.numElemZ()+1; iz++)
00073   {
00074      for (unsigned k=0; k < ver_xy; k++)
00075           vb++;
00076      //  printf("check vertex vb  %u \n ", vb->index());
00077     for (ix=0; ix<_mesh.numElemX()+1; ix++)
00078     {
00079       sp.add(vb->index(),vb->index());
00080       sp.add(_m+vb->index(),_m+vb->index());
00081       //  printf("check vertex vb  %u \n ", vb->index());
00082       vb++;     
00083     }
00084   }
00085 
00086   /*Left face boundary*/
00087   OrthoMesh::Raw_Vertex_It vl = _mesh.begin_raw_vertice();
00088    
00089   for (iz=0; iz<_mesh.numElemZ()+1; iz++)
00090   {    
00091     //  for (iy=0; iy<_mesh.numElemY()+1; iy++)
00092     for (iy=0; iy<_mesh.numElemY()+1; iy++)
00093     {  
00094       sp.add(vl->index(),vl->index());
00095       sp.add(_m+vl->index(),_m+vl->index());
00096       //  printf("vl %u \n", vl->index()); 
00097       for (unsigned k=0; k <_mesh.numElemX()+1; k++)
00098          vl++;
00099           
00100     }
00101   }
00102 
00103   /*Right face boundary*/  
00104   OrthoMesh::Raw_Vertex_It vr = _mesh.begin_raw_vertice();
00105   // printf("here vr %u \n", vr->index());
00106   for (unsigned k=0; k <_mesh.numElemX(); k++)
00107          vr++;
00108   //     printf("VR: out %u \n", vr->index());
00109      
00110   for (iz=0; iz<_mesh.numElemZ()+1; iz++)
00111   {
00112     for (iy=0; iy<_mesh.numElemY()+1; iy++)
00113     {    
00114       sp.add(vr->index(),vr->index());
00115       sp.add(_m+vr->index(),_m+vr->index());
00116       //       cout<<"vr \t"<< vr->index()<<endl;
00117       for (unsigned k=0; k <_mesh.numElemX()+1; k++)
00118       {
00119         
00120         if (vr->index() != _m-1)
00121         {
00122         vr++;
00123         //      cout<<"vertices inside\t"<< vr->index()<<endl;
00124         }
00125       }   
00126     }
00127   }
00128 
00129 
00130   /*Bottom face boundary*/ 
00131   OrthoMesh::Raw_Vertex_It vbt = _mesh.begin_raw_vertice();
00132  
00133   iz=0;
00134  
00135   for (iy=0; iy<_mesh.numElemY()+1; iy++)
00136   {    
00137     for (ix=0; ix<_mesh.numElemX()+1; ix++)
00138     {  
00139       sp.add(vbt->index(),vbt->index());
00140       sp.add(_m+vbt->index(),_m+vbt->index());
00141       //        printf("vbt %u \n", vbt->index());        
00142      vbt++;
00143     }
00144   }
00145  
00146   /*Top (up) face boundary*/  
00147   OrthoMesh::Raw_Vertex_It vu = _mesh.begin_raw_vertice();
00148   
00149   for (unsigned k=0 ; k <_m-(_ver_xy); k++)
00150     vu++;
00151   // printf("VU %u \n", vu->index());
00152   
00153   iz=_mesh.numElemZ()+1;
00154   //   if (I == -1)
00155   for (iy=0; iy<_mesh.numElemY()+1; iy++)
00156   {   
00157     for (ix=0; ix<_mesh.numElemX()+1; ix++)
00158     {
00159       sp.add(vu->index(),vu->index());
00160       sp.add(_m+vu->index(),_m+vu->index());
00161       //  printf("vu %u \n", vu->index());        
00162       vu++;
00163     }
00164   }
00165 
00166   /*Center vertices*/
00167   OrthoMesh::Raw_Vertex_It vc = _mesh.begin_raw_vertice();
00168   
00169   for (int i=0; i < (_mesh.numElemX()+1) * (_mesh.numElemY()+1) + (_mesh.numElemX()+2); i++)
00170   {
00171     vc++;
00172   }
00173   // printf("vc %u\n", vc->index());
00174   //--------------------------------------------------------------------
00175    for (iz=1; iz < _mesh.numElemZ(); iz++)
00176    {
00177     for (iy=1; iy < _mesh.numElemY(); iy++)
00178     {
00179       for (ix=1; ix < _mesh.numElemX(); ix++)
00180       {
00181         //      printf("vc %u\n", vc->index());
00182           //si_c 1st eq.
00183           sp.add(vc->index(),vc->index());
00184           //si_w 1st eq.
00185           sp.add(vc->index(),_m+vc->index());
00186           sp.add(vc->index(),_m+vc->index_left());
00187           sp.add(vc->index(), _m+ vc->index_right());
00188           sp.add(vc->index(),_m+vc->index_bottom());
00189           sp.add(vc->index(),_m+vc->index_up());
00190           sp.add(vc->index(),_m+ vc->index_back());
00191           sp.add(vc->index(),_m+ vc->index_front());
00192           //si_c 2nd eq.
00193           sp.add(_m+vc->index(),vc->index());
00194           sp.add(_m+vc->index(),vc->index_left());
00195           sp.add(_m+vc->index(), vc->index_right());
00196           sp.add(_m+vc->index(),vc->index_bottom());
00197           sp.add(_m+vc->index(),vc->index_up());
00198           sp.add(_m+vc->index(), vc->index_back());
00199           sp.add(_m+vc->index(), vc->index_front());
00200           //si_w 2nd eq.
00201           sp.add(_m+vc->index(),_m+vc->index());
00202           sp.add(_m+vc->index(),_m+vc->index_left());
00203           sp.add(_m+vc->index(),_m+ vc->index_right());
00204           sp.add(_m+vc->index(),_m+vc->index_bottom());
00205           sp.add(_m+vc->index(),_m+vc->index_up());
00206           sp.add(_m+vc->index(),_m+ vc->index_back());
00207           sp.add(_m+vc->index(),_m+ vc->index_front());       
00208           vc++;
00209         }
00210 
00211       for (unsigned i=0; i<2; i++)
00212         vc++;
00213     }
00214 
00215     for (unsigned i=0; i< 2*(_mesh.numElemX()+1); i++)
00216         vc++;
00217  
00218     } 
00219  
00220 
00221 
00222 
00223   //-----------------------------------------------------------------------
00224 
00225   /*  
00226    for (iz=0; iz < _mesh.numElemZ()+1; iz++)
00227     for (iy=0; iy < _mesh.numElemY()+1; iy++)
00228       for (ix=0; ix < _mesh.numElemX()+1; ix++)
00229       {
00230         //      cout<<"here\t"<<vc->index()<<endl;
00231         if (ix != 0 && iy != 0 && iz != 0 && ix != _mesh.numElemX() &&  iy != _mesh.numElemY() &&  iz != _mesh.numElemZ())
00232         {
00233           // printf("vc %u %u %u %u %u\n", vc->index(), vc->index_bottom(), ix, iy, iz);
00234           //si_c 1st eq.
00235           sp.add(vc->index(),vc->index());
00236           //si_w 1st eq.
00237           sp.add(vc->index(), _m+vc->index());
00238           sp.add(vc->index(), _m+vc->index_left());
00239           sp.add(vc->index(), _m+vc->index_right());
00240           sp.add(vc->index(), _m+vc->index_bottom());
00241           sp.add(vc->index(), _m+vc->index_up());
00242           sp.add(vc->index(), _m+vc->index_back());
00243           sp.add(vc->index(), _m+vc->index_front());
00244           //si_c 2nd eq.
00245           sp.add(_m+vc->index(),vc->index());
00246           sp.add(_m+vc->index(),vc->index_left());
00247           sp.add(_m+vc->index(), vc->index_right());
00248           sp.add(_m+vc->index(),vc->index_bottom());
00249           sp.add(_m+vc->index(),vc->index_up());
00250           sp.add(_m+vc->index(), vc->index_back());
00251           sp.add(_m+vc->index(), vc->index_front());
00252           //si_w 2nd eq.
00253           sp.add(_m+vc->index(),_m+vc->index());
00254           sp.add(_m+vc->index(),_m+vc->index_left());
00255           sp.add(_m+vc->index(),_m+ vc->index_right());
00256           sp.add(_m+vc->index(),_m+vc->index_bottom());
00257           sp.add(_m+vc->index(),_m+vc->index_up());
00258           sp.add(_m+vc->index(),_m+ vc->index_back());
00259           sp.add(_m+vc->index(),_m+ vc->index_front());
00260           printf("vc %u \n", vc->index());
00261         }
00262         vc++;
00263         
00264         
00265       }
00266   */
00267      sp.compress();
00268   
00269      //      printf("Sparsity end: %u \n", sp.n_nonzero_elements());
00270        
00271      // std::ofstream out ("sp1");
00272      //  sp.print_gnuplot (out);
00273      /* to save the image:
00274         set terminal png
00275         set output 'sp1.png'
00276         plot "sp1" with points */
00277            
00278 }

void BlockMatrixModule::construct_bases ( VecDouble &  bc_pc,
double  dt 
)

Definition at line 428 of file blockmatrixmodule.cpp.

00429 {
00430   //_dt=0.1;
00431   VecDouble SiPrev(_m);
00432   VecDouble SwPrev;
00433  
00434   //For each fracture block 
00435   for (unsigned blockId =0; blockId < n_fracture_cells ; blockId++)
00436   {
00437     //SiPrev = Si + (PcNew-PcOld)*SiHom
00438     SiPrev.copyRow(blockId,m_MSi); //prev capillary pressure solution
00439     //  cout<<"construct basis \n";
00440    
00441     VecDoubleRef SiHom(&(m_MHomSi(blockId,0)),_m);
00442   
00443     SiPrev.sadd(1,_BcPcNew(blockId)-_BcPcOld(blockId),SiHom);
00444   
00445     //Calculate the saturation corresponding to the SiPrev      
00446     SwPrev = SiPrev;
00447     _fpcinv.evaluate(SwPrev); // saturation at previous time level
00448 
00449     //  SwPrev.print();
00450     //Store the integral of the sw field
00451     _SwIntegral(blockId)=_mesh.getIntegralAtVertices(SwPrev);
00452     
00453     VecDoubleRef Si(&(m_MSi(blockId,0)),_m);
00454 
00455     create_bases_functions(bc_pc(blockId),SwPrev,SiPrev,_vPor(blockId),_vK(blockId),Si,SiHom);
00456     _BcPcOld(blockId)=bc_pc(blockId);
00457 
00458     HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
00459     hdf5.writeScalarField(Si,"Psi","MBmesh");
00460 
00461   }
00462     
00463     //  for (unsigned i=0; i<n_fracture_cells; i++)
00464     // {
00465     //   for (unsigned j=0; j<_m; j++)
00466     //     cout<<j<<'\t'<<m_MSi(i,j)<<endl;
00467     //   cout<<endl<<"new row"<<endl;
00468     // }
00469 
00470   // FILE *fout;
00471   //    fout= fopen("source","w");
00472            
00473   //      for (unsigned i=0; i< 200; i++)
00474   //      {
00475   //     fprintf(fout, "%u %g \n",i,get_source(0,i)  );
00476   //      }
00477 }

void BlockMatrixModule::create_bases_functions ( double  bc_pc,
const VecDouble &  previous_sw,
const VecDouble &  previous_si,
double  por,
double  k,
VecDouble &  newSi,
VecDouble &  newHomSi 
) [protected]

Create basis functions

Parameters:
dt Time step
bc Boundary Pc (capital si)
por Porosity
k Permeability
previous Previous solution
newSi To contain the non homogeneous solution
newHomSi To contain the homegenous
previous_sw : previous time level saturation
previous_si : solution at the previous time level

Definition at line 491 of file blockmatrixmodule.cpp.

00492 {
00493 //Initialize the vectors to contain the Pc from the fractures
00494 // _SwIntegral.reinit(n_fracture_cells);
00495   
00496   assembleSystem(previous_sw, previous_si,k, por, bc_pc, _dt);
00497   // store newSi and newHomSi in m_ MSi and m_HomMSi
00498    newSi=Sol;
00499    newHomSi=Sol_hom;
00500 }

double BlockMatrixModule::get_source ( unsigned  blockId,
double  new_bc_pc 
)

Definition at line 502 of file blockmatrixmodule.cpp.

00503 {
00504   VecDoubleRef PSi(&(m_MSi(blockId,0)),_m);
00505   VecDoubleRef PSiHom(&(m_MHomSi(blockId,0)),_m);
00506   
00507   _BcPcNew(blockId)=new_bc_pc;
00508   
00509   VecDouble Sol;
00510 
00511   //Sol = PSi + (PcNew - PcOld) *PSiHom
00512   Sol=PSi;
00513   
00514   Sol.sadd(1,new_bc_pc - _BcPcOld(blockId),PSiHom);
00515 
00516   //Sol =  PcInv(Sol)
00517   _fpcinv.evaluate(Sol);
00518   
00519   //Now Sol vector contains  the saturation
00520   double SwInt=_mesh.getIntegralAtVertices(Sol);
00521   return -(SwInt-_SwIntegral(blockId))*_vPor(blockId)/(_dt*_Bx);
00522 }

void BlockMatrixModule::iterate ( TransportBase &  trans  )  [virtual]

Implements DynamicBase.

Definition at line 13 of file blockmatrixmodule.cpp.

00014 {
00015   //  n_fracture_cells=125 * 2;
00016   // cout<<"frac cell "<<n_fracture_cells<<endl;
00017   VecDouble P;
00018   P.reinit(n_fracture_cells);
00019   for (unsigned i=0; i < n_fracture_cells; i++)
00020     P(i)=0.2;
00021 
00022   
00023   construct_bases (P, _dt);
00024   
00025   // std::ofstream out ("sp1");
00026   // sp.print_gnuplot (out);
00027   // abort();
00028 }

virtual void BlockMatrixModule::printOutput (  )  [inline, virtual]

Implements DynamicBase.

Definition at line 48 of file blockmatrixmodule.h.

00048 {}

void BlockMatrixModule::Sol_Si (  )  [protected]

Definition at line 403 of file blockmatrixmodule.cpp.

00404 {
00405   
00406   UMFPACKSolver umfp;
00407   umfp.solve(M,Solution1,rhs_Si);
00408 
00409   //  Solution
00410   for (unsigned i=0; i< _m ;i++)
00411   {
00412      Sol(i)=Solution1(i);
00413     //     cout<<Sol(i)<<endl;
00414   }
00415 
00416     
00417 }

void BlockMatrixModule::Sol_Sihom (  )  [protected]

Definition at line 419 of file blockmatrixmodule.cpp.

00420 {
00421   UMFPACKSolver umfp;
00422   umfp.solve(M,Solution2,rhs_Sihom);
00423   for (unsigned i=0; i< _m ; i++)
00424     Sol_hom(i)=Solution2(i);
00425 
00426  }

Member Data Documentation

VecDouble BlockMatrixModule::_BcPcNew [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::_BcPcOld [private]

Definition at line 22 of file blockmatrixmodule.h.

double BlockMatrixModule::_Bx [private]

Definition at line 25 of file blockmatrixmodule.h.

double BlockMatrixModule::_dt [private]

Definition at line 24 of file blockmatrixmodule.h.

double BlockMatrixModule::_dx [private]

Definition at line 24 of file blockmatrixmodule.h.

double BlockMatrixModule::_dy [private]

Definition at line 24 of file blockmatrixmodule.h.

double BlockMatrixModule::_dz [private]

Definition at line 24 of file blockmatrixmodule.h.

Function1D& BlockMatrixModule::_fdpc [private]

Definition at line 28 of file blockmatrixmodule.h.

Function1D & BlockMatrixModule::_fmobrW [private]

Definition at line 28 of file blockmatrixmodule.h.

Function1D & BlockMatrixModule::_fmobt [private]

Definition at line 28 of file blockmatrixmodule.h.

Function1D & BlockMatrixModule::_fpcinv [private]

Definition at line 28 of file blockmatrixmodule.h.

unsigned BlockMatrixModule::_m [private]

Definition at line 19 of file blockmatrixmodule.h.

OrthoMesh& BlockMatrixModule::_mesh [private]

Definition at line 18 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::_SwIntegral [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::_vK [private]

Definition at line 27 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::_vPcInitial [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::_vPor [private]

Definition at line 26 of file blockmatrixmodule.h.

SparseMatrix<double> BlockMatrixModule::M [private]

Definition at line 15 of file blockmatrixmodule.h.

Matrix BlockMatrixModule::m_MHomSi [private]

Definition at line 17 of file blockmatrixmodule.h.

Matrix BlockMatrixModule::m_MSi [private]

Definition at line 17 of file blockmatrixmodule.h.

unsigned BlockMatrixModule::n_fracture_cells [private]

Definition at line 20 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::rhs_Si [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::rhs_Sihom [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::sat_old [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::Sol [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::Sol_hom [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::Sol_hom_old [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::Sol_old [private]

Definition at line 22 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::Solution1 [private]

Definition at line 23 of file blockmatrixmodule.h.

VecDouble BlockMatrixModule::Solution2 [private]

Definition at line 23 of file blockmatrixmodule.h.

SparsityPattern BlockMatrixModule::sp [private]

Definition at line 14 of file blockmatrixmodule.h.

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Sun Apr 8 23:12:58 2012 for CO2INJECTION by  doxygen 1.6.3