MixedHybridBase Class Reference

#include <mixedhybridbase.h>

Inheritance diagram for MixedHybridBase:

[legend]

List of all members.

Classes
struct	_HybridFaceBC
struct	_ScalarBC
Protected Types
enum	BCType { Dirichlet, Neumann }
typedef std::vector< _ScalarBC >	ScalarBC
typedef std::vector < _HybridFaceBC >	HybridFaceBC
Protected Member Functions
void	setBoundaryCondition (HybridFaceBC &vBC, OrthoMesh &mesh, Function3D &fInwardVel, Function3D &fP, bool bCheckNoBC)
void	addPcTerm (VecDouble &B, OrthoMesh &mesh, ScalarBC &bc, const VecDouble &KPc, const VecDouble &vPc)
void	assemblySystem (SparseMatrix< double > &A, VecDouble &B, OrthoMesh &mesh, HybridFaceBC &bc, const VecDouble &K, const VecDouble &vMobT, const VecDouble &vGravTerm)
void	assemblySystem (SparseMatrix< double > &A, VecDouble &B, OrthoMesh &mesh, HybridFaceBC &bc, const VecDouble &K)
void	assemblySystem (SparseMatrix< double > &A, VecDouble &B, OrthoMesh &mesh, ScalarBC &bc, const VecDouble &vCoeff)
void	getNormalVelocitiesFromPressure (VecDouble &Vq, OrthoMesh &mesh, HybridFaceBC &bc, const VecDouble &K, const VecDouble &vP, const VecDouble &vMobT, const VecDouble &vGravTerm)
void	getNormalVelocitiesFromPressure (VecDouble &Vq, OrthoMesh &mesh, HybridFaceBC &bc, const VecDouble &K, const VecDouble &vP)
void	getNormalVelocitiesFromPressure (VecDouble &Vq, OrthoMesh &mesh, HybridFaceBC &bc, const VecDouble &K, const VecDouble &vP, const VecDouble &vMobT, const VecDouble &vGravTerm, ScalarBC &vPcBC, const VecDouble &KPc, const VecDouble &vPc)
void	buildSparsityPattern (OrthoMesh &mesh, SparsityPattern &sparsePattern)
	MixedHybridBase ()
	~MixedHybridBase ()

---

Detailed Description

This class implements the hybrid aproach for the pressure equation of the biphasic problem

div(Vdt) = 0 Vdt = -K mobT(Sw) Grad(P) + K mobT(Sw) mobGrav(Sw) *Grad(Z)

Where Vdt is total percolation velocity, P is the global pression K is the permeability, Sw is the saturation of water, mobT() is the total mobility function and mobGrav(Sw) is the mobility function belonging to the gravity effect and in the biphasic flow it is written in termos of the mobilities of the two fluids (mobW,mobO) by the expression

mobGrav(Sw) = (po mobO(Sw) + pw mobW(Sw))*g

where pw,po is the density of the two fluids and g is the acceleration of gravity.

The variational equation guives

[Vdt / (K mobT(Sw) ), v] - [p,div(v)] = -<p,v> + [ mobGrav(Sw) Grad(Z), v]

[div(Vdt), q] = 0

Where q,v are the test functions of pression and velocity Vdt, <p,v> means the surface integral of the vector field p*v around the boundary (Neumman Condition) and [u,v] means integral of u*v function over all the problem's domain

For pratical reasons we always assume that the flux due to gravity terms in the boundary is always zero.

Definition at line 53 of file mixedhybridbase.h.

---

Member Typedef Documentation

typedef std::vector<_HybridFaceBC> MixedHybridBase::HybridFaceBC [protected]

Definition at line 107 of file mixedhybridbase.h.

typedef std::vector<_ScalarBC> MixedHybridBase::ScalarBC [protected]

Definition at line 65 of file mixedhybridbase.h.

---

Member Enumeration Documentation

enum MixedHybridBase::BCType [protected]

Enumerator:

Dirichlet
Neumann

Definition at line 68 of file mixedhybridbase.h.

{Dirichlet,Neumann};

---

Constructor & Destructor Documentation

MixedHybridBase::MixedHybridBase

(

)

[inline, protected]

Definition at line 131 of file mixedhybridbase.h.

{}

MixedHybridBase::~MixedHybridBase

(

)

[inline, protected]

Definition at line 132 of file mixedhybridbase.h.

{}

---

Member Function Documentation

void MixedHybridBase::addPcTerm	(	VecDouble &	B,
		OrthoMesh &	mesh,
		ScalarBC &	bc,
		const VecDouble &	KPc,
		const VecDouble &	vPc
	)			[protected]

Add the term - KPc* Grad Pc to the system equation. Since we know Pc it goes to the RHS as + KPc Grad Pc

Parameters:

	B	RHS
	mesh
	bc	Boundary Conditions for Cappilar pressure
	KeffPc	Permeability KPc for each cell.
	vPc	Cappilar pressure for each cell.

Returns:

Definition at line 241 of file mixedhybridbase.cpp.

{
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = mesh.begin_face();
  unsigned nFaces = mesh.numFaces();
    
  //For each face
  for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++)
  {
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);
    
    //Make sure the face is not at boundary
    if (!face->at_boundary())
    {
      double K1 = KPc(c1);
      double K2 = KPc(c2);
      double Mh_2 = _div(K1*K2,(K1+K2)); //Harmonic media divided by two
      double Keff =2*Mh_2 * face->area()*face->areaPerCellVol();
      double q = -Keff*(vPc(c2)-vPc(c1));
      B(c1) += -q;
      B(c2) += +q;
    }
  }

  ScalarBC::iterator bcEnd = bc.end();
  for(ScalarBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
    {
      
      OrthoMesh::Face_It face = bcIt->face;
      assert(face->at_boundary());
      
      unsigned c1,c2;
      face->getAdjCellIndices(c1,c2);

      double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2;
      double Keff = 2*KPc(c)*face->area()*face->areaPerCellVol();
      B(c) += Keff*(bcIt->value-vPc(c));
    }

  
}

void MixedHybridBase::assemblySystem	(	SparseMatrix< double > &	A,
		VecDouble &	B,
		OrthoMesh &	mesh,
		ScalarBC &	bc,
		const VecDouble &	vCoeff
	)			[protected]

void MixedHybridBase::assemblySystem	(	SparseMatrix< double > &	A,
		VecDouble &	B,
		OrthoMesh &	mesh,
		HybridFaceBC &	bc,
		const VecDouble &	vCoeff
	)			[protected]

Assembly the mixed finite element system given by div(-Coeff Grad P) = 0

Definition at line 124 of file mixedhybridbase.cpp.

{
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = mesh.begin_face();
  unsigned nFaces = mesh.numFaces();
    
  //For each face
  for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++)
  {
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);
    
    //Make sure the face is not at boundary
    if (!face->at_boundary())
    {
      double K1 = K(c1);
      double K2 = K(c2);
      double Mh_2 = _div(K1*K2,(K1+K2)); //Harmonic media divided by two
      double Keff =2*Mh_2 * face->area()*face->areaPerCellVol();

 
      
      A.add(c1,c1, Keff);
      A.add(c1,c2,-Keff);
      A.add(c2,c2, Keff);
      A.add(c2,c1,-Keff);
    }
  }
  //Now process the boundary conditions
    
    HybridFaceBC::iterator bcEnd = bc.end();
    for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
    {
      OrthoMesh::Face_It face = bcIt->face;
      assert(face->at_boundary());
      
      unsigned c1,c2;
      face->getAdjCellIndices(c1,c2);

      if (bcIt->bcType == Dirichlet)
      {
        double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2;
        double Keff = 2*K(c)*face->area()*face->areaPerCellVol();
        A.add(c,c,Keff);
        B(c) += Keff*bcIt->value;
      }
      else //Neumann condition
      {
        if (c1 != OrthoMesh::invalidIndex())
          B(c1) += -bcIt->value;
        else
          B(c2) += +bcIt->value;
      }
    }

}

void MixedHybridBase::assemblySystem	(	SparseMatrix< double > &	A,
		VecDouble &	B,
		OrthoMesh &	mesh,
		HybridFaceBC &	bc,
		const VecDouble &	K,
		const VecDouble &	vMobT,
		const VecDouble &	vGravTerm
	)			[protected]

for each cell we calculate the sum q0 + q1 + q2 + q3 + q4 + q5 where qi means the out flux of the cell for each face i. The flux is obtained from the velocity field by qi = Vi*Nout*dS where Vi is the velocity, Nout the outward normal of the face and dS is the face's area. Vi is obtained from the pressure by the equation Vi = -Keff (p2'-p1)/h - Keff (GravTerm(c1) + GravTerm(c2))/2 where p1 is the pressure of the positive cell (nonzero component of the outward normal is positite), p2 is the pressure of the negative cell and c1, c2 are the indices of the positive and negative cells. Instead of iterate using the cells, we asseble the matrix iterating with the faces and putting the contribution of each face in both cells at once

---*---* | | | | p1| p2| ---*---* NOTE THAT WE DONT ZERO THE A AND B entries. WE JUST ADD TERMS TO IT.

Parameters:

	A	Matrix to be mounted
	B	RHS
	mesh	The mesh
	bc	Structure containing the boundary conditions
	K	The permeability for each cell
	vMobT	The vector containing the total mobility for each cell.
	vGravTerm	The vector containing the gravity term for each cell

Returns:

Definition at line 56 of file mixedhybridbase.cpp.

{


  
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = mesh.begin_face();
  unsigned nFaces = mesh.numFaces();
    
  //For each face
  for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++)
  {
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);
    
    //Make sure the face is not at boundary
    if (!face->at_boundary())
    {
      double K1 = K(c1)*vMobT(c1);
      double K2 = K(c2)*vMobT(c2);
      double Mh_2 = K1*K2/(K1+K2); //Harmonic media divided by two
      double Keff =2*Mh_2 * face->area()*face->areaPerCellVol();

 
      
      A.add(c1,c1, Keff);
      A.add(c1,c2,-Keff);
      A.add(c2,c2, Keff);
      A.add(c2,c1,-Keff);
      if (face->getNormalOrientation() == OrthoMesh::NORMAL_Z)
      {
        double aux=+Mh_2*(vGravTerm(c1)+vGravTerm(c2))*face->area();
        B(c1) += -aux;
        B(c2) +=  aux;
      }
    }
  }
  //Now process the boundary conditions
    
    HybridFaceBC::iterator bcEnd = bc.end();
    for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
    {
      OrthoMesh::Face_It face = bcIt->face;
      assert(face->at_boundary());
      
      unsigned c1,c2;
      face->getAdjCellIndices(c1,c2);

      if (bcIt->bcType == Dirichlet)
      {
        double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2;
        double Keff = 2*K(c)*vMobT(c)*face->area()*face->areaPerCellVol();
        A.add(c,c,Keff);
        B(c) += Keff*bcIt->value;
      }
      else //Neumann condition
      {
        if (c1 != OrthoMesh::invalidIndex())
          B(c1) += -bcIt->value;
        else
          B(c2) += +bcIt->value;
      }
    }
}

void MixedHybridBase::buildSparsityPattern	(	OrthoMesh &	mesh,
		SparsityPattern &	sparsePattern
	)			[protected]

Definition at line 5 of file mixedhybridbase.cpp.

{
  OrthoMesh::Cell_It cell = mesh.begin_cell();
  OrthoMesh::Cell_It endc = mesh.end_cell();
  sparsePattern.reinit(mesh.numCells(),mesh.numCells(),7);

  
  for(;cell!=endc;cell++)
  {
    int ci = cell->index();
    sparsePattern.add(ci,ci);
    //For each boundary face of the cell
    for (unsigned int face_no = 0; face_no < GeometryInfo<DIM>::faces_per_cell;++face_no)
    {
      unsigned adj_index  = cell->neighbor_index(face_no);
      if (adj_index != OrthoMesh::INVALID_INDEX)
        sparsePattern.add(ci,adj_index);
    }
  }
  sparsePattern.compress();
}

void MixedHybridBase::getNormalVelocitiesFromPressure	(	VecDouble &	Vq,
		OrthoMesh &	mesh,
		HybridFaceBC &	bc,
		const VecDouble &	K,
		const VecDouble &	vP,
		const VecDouble &	vMobT,
		const VecDouble &	vGravTerm,
		ScalarBC &	vPcBC,
		const VecDouble &	KPc,
		const VecDouble &	vPc
	)			[protected]

Definition at line 424 of file mixedhybridbase.cpp.

{
  assert(Vq.size() >= mesh.numFaces());
  assert(vP.size() >= mesh.numCells());
  assert(vMobT.size() >= mesh.numCells());
  assert(vGravTerm.size() >=  mesh.numCells());
  assert(KPc.size() >= mesh.numCells());
  assert(vPc.size() >= mesh.numCells());
  assert(K.size() >= mesh.numCells());
  
  
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = mesh.begin_face();
  OrthoMesh::Face_It endf = mesh.end_face();

    
  //For each face
  for (;face!=endf;face++)
  {
    if (!face->at_boundary())
    {
      unsigned c1,c2;
      face->getAdjCellIndices(c1,c2);
      double K1 = K(c1)*vMobT(c1);
      double K2 = K(c2)*vMobT(c2);
      double Keff = 2*K1*K2/(K1+K2)*face->areaPerCellVol();
      Vq(face->index()) = - Keff * (vP(c2) - vP(c1));

      //Now for cappillary pressure.
      Keff=2*_div(KPc(c1)*KPc(c2),(KPc(c1)+KPc(c2)))*face->areaPerCellVol();
      Vq(face->index()) += -Keff*(vPc(c2)-vPc(c1));
      if (face->getNormalOrientation() == OrthoMesh::NORMAL_Z)
      {
        Vq(face->index()) += +K1*K2/(K1+K2)*(vGravTerm(c1)+vGravTerm(c2));
      }      //printf("%g ",-Keff*(vPc(c2)-vPc(c1)));

      
    }
  }

  //Now for boundary conditions
    
  HybridFaceBC::iterator bcEnd = bc.end();
  for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
  {
    OrthoMesh::Face_It face = bcIt->face;
    assert(face->at_boundary());
      
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);

    if (bcIt->bcType == Dirichlet)
    {
      if (c1 != OrthoMesh::invalidIndex())
      {
        double Keff = 2*K(c1)*vMobT(c1)*face->areaPerCellVol();
        Vq(face->index()) = - Keff*(bcIt->value - vP(c1)); 
        //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
      }
      else
      {
        double Keff = 2*K(c2)*vMobT(c2)*face->areaPerCellVol();
        Vq(face->index()) = Keff*(bcIt->value - vP(c2));
        //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
      }
    }
    else
    {
      if (c1 != OrthoMesh::invalidIndex())
        Vq(face->index()) = bcIt->value/face->area();
      else
        Vq(face->index()) = bcIt->value/face->area();
    }
  }


  ScalarBC::iterator ScalarBCEnd = vScalarBC.end();
  for(ScalarBC::iterator bcIt = vScalarBC.begin();bcIt!=ScalarBCEnd;bcIt++)
  {
    OrthoMesh::Face_It face = bcIt->face;
    assert(face->at_boundary());
      
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);

      if (c1 != OrthoMesh::invalidIndex())
      {
        double Keff = 2*KPc(c1)*face->areaPerCellVol();
        Vq(face->index()) += - Keff*(bcIt->value - vP(c1)); 
        //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
      }
      else
      {
        double Keff = 2*KPc(c2)*face->areaPerCellVol();
        Vq(face->index()) += - Keff*(vPc(c2) - bcIt->value);
        //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
      }
  }
  
}

void MixedHybridBase::getNormalVelocitiesFromPressure	(	VecDouble &	Vq,
		OrthoMesh &	mesh,
		HybridFaceBC &	bc,
		const VecDouble &	K,
		const VecDouble &	vP
	)			[protected]

Definition at line 356 of file mixedhybridbase.cpp.

{
  assert(Vq.size() == mesh.numFaces());
  assert(vP.size() == mesh.numCells());
  assert(K.size() == mesh.numCells());
  
  
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = mesh.begin_face();
  OrthoMesh::Face_It endf = mesh.end_face();

    
  //For each face
  for (;face!=endf;face++)
  {
    if (!face->at_boundary())
    {
      unsigned c1,c2;
      face->getAdjCellIndices(c1,c2);
      double K1 = K(c1);
      double K2 = K(c2);
      double Keff = 2*_div(K1*K2,(K1+K2))*face->areaPerCellVol();
      Vq(face->index()) = - Keff * (vP(c2) - vP(c1));
    }
  }

  //Now for boundary conditions
    
  HybridFaceBC::iterator bcEnd = bc.end();
  for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
  {
    OrthoMesh::Face_It face = bcIt->face;
    assert(face->at_boundary());
      
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);

    if (bcIt->bcType == Dirichlet)
    {
      if (c1 != OrthoMesh::invalidIndex())
      {
        double Keff = 2*K(c1)*face->areaPerCellVol();
        Vq(face->index()) = - Keff*(bcIt->value - vP(c1)); 
        //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
      }
      else
      {
        double Keff = 2*K(c2)*face->areaPerCellVol();
        Vq(face->index()) = Keff*(bcIt->value - vP(c2));
        //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
      }
    }
    else
    {
      if (c1 != OrthoMesh::invalidIndex())
        Vq(face->index()) = bcIt->value/face->area();
      else
        Vq(face->index()) = bcIt->value/face->area();
    }
  }

}

void MixedHybridBase::getNormalVelocitiesFromPressure	(	VecDouble &	Vq,
		OrthoMesh &	mesh,
		HybridFaceBC &	bc,
		const VecDouble &	K,
		const VecDouble &	vP,
		const VecDouble &	vMobT,
		const VecDouble &	vGravTerm
	)			[protected]

Definition at line 285 of file mixedhybridbase.cpp.

{
  assert(Vq.size() == mesh.numFaces());
  assert(vP.size() == mesh.numCells());
  assert(vMobT.size() == mesh.numCells());
  assert(vGravTerm.size() == mesh.numCells());
  assert(K.size() == mesh.numCells());
  
  
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = mesh.begin_face();
  OrthoMesh::Face_It endf = mesh.end_face();

    
  //For each face
  for (;face!=endf;face++)
  {
    if (!face->at_boundary())
    {
      unsigned c1,c2;
      face->getAdjCellIndices(c1,c2);
      double K1 = K(c1)*vMobT(c1);
      double K2 = K(c2)*vMobT(c2);
      double Keff = 2*K1*K2/(K1+K2)*face->areaPerCellVol();
      Vq(face->index()) = - Keff * (vP(c2) - vP(c1));
      if (face->getNormalOrientation() == OrthoMesh::NORMAL_Z)
      {
        Vq(face->index()) += +K1*K2/(K1+K2)*(vGravTerm(c1)+vGravTerm(c2));
      }
      
    }
  }

  //Now for boundary conditions
    
  HybridFaceBC::iterator bcEnd = bc.end();
  for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
  {
    OrthoMesh::Face_It face = bcIt->face;
    assert(face->at_boundary());
      
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);

    if (bcIt->bcType == Dirichlet)
    {
      if (c1 != OrthoMesh::invalidIndex())
      {
        double Keff = 2*K(c1)*vMobT(c1)*face->areaPerCellVol();
        Vq(face->index()) = - Keff*(bcIt->value - vP(c1)); 
        //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
      }
      else
      {
        double Keff = 2*K(c2)*vMobT(c2)*face->areaPerCellVol();
        Vq(face->index()) = Keff*(bcIt->value - vP(c2));
        //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
      }
    }
    else
    {
      if (c1 != OrthoMesh::invalidIndex())
        Vq(face->index()) = bcIt->value/face->area();
      else
        Vq(face->index()) = bcIt->value/face->area();
    }
  }
}

void MixedHybridBase::setBoundaryCondition	(	HybridFaceBC &	vBC,
		OrthoMesh &	mesh,
		Function3D &	fInwardVel,
		Function3D &	fP,
		bool	bCheckNoBC
	)			[protected]

Definition at line 528 of file mixedhybridbase.cpp.

{
  vBC.reserve(mesh.numFacesAtBoundary());
  Point3D p;  

  OrthoMesh::Cell_It cell = mesh.begin_cell();
  OrthoMesh::Cell_It endc = mesh.end_cell();

  _HybridFaceBC bc;
  
   //For eahc cell in the Mesh
   for(;cell != endc;cell++)
   {
     //For each boundary face of the cell
     for (unsigned int face_no = 0; face_no < GeometryInfo<3>::faces_per_cell;++face_no)
     {
       OrthoMesh::Face_It face = cell->face((FaceDirection3D) face_no);
       if (face->at_boundary())
       {
         face->barycenter(p);
         if (fInwardVel.isInDomain(p)) //Neumman condition
         {
           //Set boundary condition entry 
           bc.face = face; //store the face index
           bc.bcType = Neumann; //Store the type of condition (Neumman)
           if (mesh.nonZeroInwardNormalComponenentIsPositive((FaceDirection3D) face_no))
             bc.value = fInwardVel(p)*face->area();
           else
             bc.value =- fInwardVel(p)*face->area();
           vBC.push_back(bc);
         }
         else if (fP.isInDomain(p)) //Dirichlet condition
         {
           //Set boundary condition entry 
           bc.face = face;
           bc.bcType = Dirichlet;
           bc.value= fP(p);
           vBC.push_back(bc);
         }
         else if (bCheckNoBC)
           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]);
         
       }
     }
   }

}

---

The documentation for this class was generated from the following files:

• mixedhybridbase.h
• mixedhybridbase.cpp