mixedhybridcompressiblenewton.cpp

Go to the documentation of this file.

#include "mixedhybridcompressiblenewton.h"


MixedHybridCompressibleNewton::MixedHybridCompressibleNewton(OrthoMesh &mesh,Function3D &fPrescribedVelocities, Function3D &fPrescribedPression, const VecDouble &K,  Function1D &g, Function1D &dg,LinearSolver &solver) 
  :m_mesh(mesh),
   m_g(g),
   m_dg(dg),
   m_vK(K),
   m_grav(grav),
   m_solver(solver),
   _debugLevel(debugLevel)
{
   this->setBoundaryCondition(m_bc,m_mesh,fPrescribedVelocities,fPrescribedPression,true);
  this->buildSparsityPattern(mesh,m_sparsePattern);
  m_A.reinit(m_sparsePattern);
  m_B.reinit(m_mesh.numCells());
  m_Pn.reinit(m_mesh.numCells(),NAN); 
  m_Pnn.reinit(m_mesh.numCells(),NAN);

}

 void MixedHybridCompressibleNewton::iterate(TransportBase &trans) 
{
  //Zero entries 
  A = 0;
  B = 0;
  //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())
    {
      const double P1 = Pn(c1);
      const double P2 = Pn(c2);
      const double K1 = K(c1);
      const double K2 = K(c2);
      const double Kh1 = K1*g(P1);
      const double Kh2 = K2*g(P2);
      const double Mh = 2*Kh1*Kh2/(Kh1+Kh2); //Harmonic average
      const double Keff =Mh*Inv_DX2;
      const double vv = Keff*(P1-P2);
      const double den = (Kh1+Kh2)*(Kh1+Kh2); 
      const double d1 = Kh2*Kh2*K1*dg(P1)/den*(P1-P2)*Inv_DX2;
      const double d2 = Kh1*Kh1*K2*dg(P2)/den*(P1-P2)*Inv_DX2;

 
      
      A.add(c1,c1,  d1+Keff);
      A.add(c1,c2,  d2-Keff);
      A.add(c2,c2, -d2+Keff);
      A.add(c2,c1, -d1-Keff);
      
      B(c1)-=vv;
      B(c2)+=vv;
    }
  }
  //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 vv 
        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;
      }
    }

}

---