assemblyorthomeshsystemII.h

Go to the documentation of this file.

#ifndef  _MY_AssemblySystemII_
#define  _MY_AssemblySystemII_
#include "numericmethods.h"
#include "assemblysystembase.h"
#include "matrix.h"
#include "mappedfevalues.h"
#include "unitcube.h"


/*----------------- CLASS AssemblyOrthoMeshSystem --------------
  ----------------------------------------*/

template <class EqVarType,class Matriz = SparseMatrix<double>, class Vetor = VecDouble>
  class AssemblyOrthoMeshSystemII : public AssemblySystemBase
  {
  private:
    Matrix m_Ml;
    VecDouble  m_Bl;
  

  protected:


  public:
    AssemblyOrthoMeshSystemII(){}
    ~AssemblyOrthoMeshSystemII(){}
  
  void buildGlobalMatrix(Matriz& GM, EqVarType& eqvar)
  {
    GM=0;
    assert(eqvar.getMesh());
    OrthoMesh &mesh=*(eqvar.getMesh());


    //Get the sparsity pattern.
    const Matrix &MSparsity = eqvar.get_local_sparsity();

    assert(MSparsity.m() == MSparsity.n());
    unsigned n_rank = MSparsity.m();
    //Initialize  global and local matrix.
    //The size of the local system is equal to the MSparsity
    //rank
    m_Ml.reinit(n_rank,n_rank);
    //For all cells
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();
    for (; cell!=endc; cell++)
    {
      eqvar.reinit(cell);
      const VecIndex &local_global_map = eqvar.getFE().get_global_map(cell);

      eqvar.assembly_local_system(cell,m_Ml);


      //Distribute the local matrix into the global matrix
      for (unsigned i=0;i<n_rank;i++)
      {
        for (unsigned j=0;j<n_rank;j++)
        {
          if (MSparsity(i,j) != 0.0) 
            GM.add(local_global_map[i],local_global_map[j],m_Ml(i,j));
        }
      }
    }
    m_Ml.print(std::cout);
  }

      void  buildGlobalRHS(Vetor& vRHS,EqVarType& eqvar)
    {

    assert(eqvar.getMesh());
    OrthoMesh &mesh=*(eqvar.getMesh());
    vRHS=0;

    //Get the rank of the local system
    unsigned n_rank = eqvar.get_local_sparsity().m();



    //Initialize  local RHS
    m_Bl.reinit(n_rank);

    //For all cells
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();
    for (; cell!=endc; cell++)
    {
      eqvar.reinit(cell);
      const VecIndex &local_global_map = eqvar.getFE().get_global_map(cell);

      //Build local rhs
      eqvar.assembly_local_rhs(cell,m_Bl);


      //Distribute the local matrix into the global matrix

      for (unsigned i=0;i<n_rank;i++)
      {
        vRHS(local_global_map[i])+=m_Bl(i);
      }

    }

    }


  Index global_system_rank(EqVarType &eqvar)
  {
    
    Index max=0;
    assert(eqvar.getMesh() != NULL);
    OrthoMesh &mesh = *(eqvar.getMesh());
    //For all cells
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();
    for (; cell!=endc; cell++)
    {
      eqvar.reinit(cell);
      const VecIndex &local_global_map = eqvar.getFE().get_global_map(cell);
      for (unsigned i=0;i<local_global_map.size();i++)
      {
        if (max < local_global_map[i])
          max = local_global_map[i];
      } 
    }
    return max+1;
  }


 
  Index max_nz_per_row(EqVarType& eqvar)
  {
    Point3D P1(0,0,0),P2(1,1,1);
    
    //Create a small mesh
    OrthoMesh mesh(P1,P2,3,3,3);

    //Save the mesh already assigned to the eqvar
    OrthoMesh *backMesh = eqvar.getMesh();

    //Now initialize the finite element with the new mesh.
    eqvar.reinit(mesh);

    //Get the global system ranks
    //for the small system of a mesh 3x3x3
    Index global_rank = global_system_rank(eqvar);

    //Initialize the sparse pattern as a full sparse pattern
    //(i.e The number of non zero columns per row is the
    //number of columns)
    SparsityPattern pattern;
    pattern.reinit(global_rank,global_rank,global_rank);
    
    //Get the local sparsity
    const Matrix &localSparsity = eqvar.get_local_sparsity();
    assert(localSparsity.m() == localSparsity.n());
    unsigned n_rank=localSparsity.m();
    //Now a loop to build the sparsity pattern
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();
    for (; cell!=endc; cell++)
    {
      eqvar.reinit(cell);
      const VecIndex &local_global_map = eqvar.getFE().get_global_map(cell);
      for (unsigned i=0;i<n_rank;i++)
        for (unsigned j=0;j<n_rank;j++)
        {
          if (localSparsity(i,j) != 0.0)
            pattern.add(local_global_map[i],local_global_map[j]);
        }
    }
    //Ok, now revert the finite element with the previous mesh
    eqvar.reinit(*backMesh);


    //Finally compress the pattern
    pattern.compress();

    //After compressed the pattern contains
    //the exact number of  non zero columns per row
    //return this number
    return pattern.max_entries_per_row();
  }
  
  template<class SparsePattern>
  void buildSparsityPattern(SparsePattern &pattern,EqVarType& eqvar)
  {

    assert(eqvar.getMesh());
    OrthoMesh &mesh=*(eqvar.getMesh());

    //Get the global system ranks
    Index global_rank = global_system_rank(eqvar);

    //Initialize the sparse pattern calculating the
    //max non zero per rows with the function max_nz_per_row()
    pattern.reinit(global_rank,
                   global_rank,
                   this->max_nz_per_row(eqvar));
    
    //Get the local sparsity
    const Matrix &localSparsity = eqvar.get_local_sparsity();
    assert(localSparsity.m() == localSparsity.n());
    unsigned n_rank=localSparsity.m();

    
    //Now a loop to build the sparsity pattern
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();
    for (; cell!=endc; cell++)
    {
      eqvar.reinit(cell);
      const VecIndex &local_global_map = eqvar.getFE().get_global_map(cell);
      for (unsigned i=0;i<n_rank;i++)
        for (unsigned j=0;j<n_rank;j++)
        {
          if (localSparsity(i,j) != 0.0)
            pattern.add(local_global_map[i],local_global_map[j]);
        }
    }


    //Finally compress the pattern
    pattern.compress();

    //Done
  }
  


  void    addNeumannCondition(Vetor &vRHS,EqVarType& eqvar,Function3D &fN)
  {
    assert(eqvar.getMesh());
    OrthoMesh &mesh=*(eqvar.getMesh());
    
    Point3D pt;
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();

    VecDouble Bl(eqvar.getFE().n_dofs_per_face());
    
    for (; cell!=endc; cell++)//For all cells
    {
      if (cell->at_boundary()) //Verify if cell is at boundary
      {
        eqvar.reinit(cell);//If it is, set the cell and find faces at boundary

        //Iterate in all faces of the cell finding the cells at boundary
        for (unsigned face_dir=0;face_dir<OrthoMesh::FACES_PER_CELL;face_dir++)
        {
          OrthoMesh::Face_It face = cell->face(static_cast<FaceDirection3D>(face_dir));
          if (face->at_boundary()) //found face at boundary  
          {
            if (fN.isInDomain(face->barycenter()))
            {
              const VecIndex& global_map = eqvar.getFE().get_global_map(cell);
              const VecIndex& face_local_map = eqvar.getFE().face_to_cell_local_map(face_dir);

              //get the neuwmann condition
              eqvar.neumann_condition(Bl,cell,face,face_dir,fN);


              //For all local dofs of the face
              for (unsigned face_dof=0;face_dof<face_local_map.size();face_dof++)
              {
                //get what is the dof number of the face with respect to its local number
                //in terms of the cell
                unsigned local_cell_dof=face_local_map[face_dof];
                vRHS(global_map[local_cell_dof])+=Bl(face_dof);

              }
            }
          }
        }
      }
    }
  }


  void setBoundaryValue(const OrthoMesh &mesh,EqVarType& eqvar,Function3D &f,MapIntDouble &map_boundary)
  {
    assert(&mesh == eqvar.getMesh());
    VecDouble pt(3);
    VecDouble ptRef(3);
    MappingOrthoMesh map;
    map.reinit(mesh);
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();
    for (; cell!=endc; cell++)//For all cells
    {
      if (cell->at_boundary()) //Verify if cell is at boundary
      {
        eqvar.reinit(cell);//If it is, set the cell 

        //For all faces of this cell
        for (unsigned face_no=0;face_no<OrthoMesh::FACES_PER_CELL;face_no++)
        {
          //if face is at boundary
          if (cell->face(static_cast<FaceDirection3D>(face_no))->at_boundary())
          {
            //See if the barycenter of the face belongs to the domain
            //of the boundary condition function
            if (f.isInDomain(cell->face(static_cast<FaceDirection3D>(face_no))->barycenter()))
            {
              map.reinit(cell);
              const VecIndex &global_map=eqvar.getFE().get_global_map(cell);
              const VecIndex &face_dofs= eqvar.getFE().face_to_cell_local_map(face_no);
              //For all local dofs at face
              for (Index local_face_dof=0;local_face_dof<face_dofs.size();local_face_dof++)
              {
                //get what is the dof number of the face with respect to its local number
                //in terms of the cell
                Index cell_dof = face_dofs[local_face_dof]; 

                //Get the support point for this dof in ptRef
                ptRef=eqvar.getFE().get_support_points()[cell_dof];

                //Map to the point in the  original domain
                map.T(ptRef,pt);
                
                //Set the value of the dof with the value of the
                //boundary condition function evaluated in the
                //support point of the deformed domain
                map_boundary[global_map[cell_dof]]=f(pt,eqvar.getFE().local_dof_field(cell_dof));
              }
            }
          }
        }
      }
    }
  }

    
  void getDOFSolAtVertices(VecDouble &vSol,const VecDouble &dofSol,EqVarType &eqvar,unsigned component=0)
  {
    //Set a list of points in the vertices
    
    MappedFEValues feValues(eqvar.getFE());;
    UnitCube refDomain;
    VecDouble weights(refDomain.n_vertices());
    weights=1;
    feValues.setPoints(refDomain.vertices(),weights);
    feValues.compute_values();
    

    //Initially put the sol values with INFINITY value;
    vSol=INFINITY;

    //For all cells of the eqvar's mesh
    assert(eqvar.getMesh());
    OrthoMesh &mesh = *(eqvar.getMesh());
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();
    for (; cell!=endc; cell++)//For all cells
    {
      //get the global map;
      const VecIndex &global_map=  eqvar.getFE().get_global_map(cell);

      //for all vertices
      for (unsigned v=0;v<refDomain.n_vertices();v++)
      {
        //Get vertice on the mesh
        unsigned vIndex=cell->vertex_index(static_cast<VertexDirection3D>(v));
        //Check if the solution was alread set
        if (vSol(vIndex) == INFINITY)
        {
          //Ok the position for this vertice was not set
          //do it.
          vSol(vIndex)=0;
          for (unsigned fId=0;fId<global_map.size();fId++)
          {
            vSol(vIndex)+=feValues.function_value(v,fId)*dofSol(global_map[fId]);
          }
        }
      }

    }
  }
    
  



  };


#endif

---