assemblyorthomeshsystem.h

Go to the documentation of this file.

#ifndef  _MY_AssemblySystem_
#define  _MY_AssemblySystem_
#include "numericmethods.h"
#include "assemblysystembase.h"



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

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

  protected:
 void distribute_local_to_global(Matriz&GM,FullMatrix<double> &lM,const std::vector<unsigned> &local_dof_indices,const bool checkZero,double tolerance=0.0)
  {
    const unsigned dofs_per_cell = local_dof_indices.size();
    if (checkZero)
    {
      //Distribute the global matrix`s nodes. 
      for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
        {
          if (fabs(lM(i,j))<=tolerance)
            continue;
          else
          {
#ifdef DEBUG
            if (GM.get_sparsity_pattern().exists(local_dof_indices[i],local_dof_indices[j]))
            {
              GM.add (local_dof_indices[i],local_dof_indices[j],lM(i,j));
            }
            else
              printf("Trying to put value %g in invalid entry (%d,%d)\n",lM(i,j),local_dof_indices[i],local_dof_indices[j]);
#else
              GM.add (local_dof_indices[i],local_dof_indices[j],lM(i,j));
#endif  
          }
        }
    }
    else
    {
      //Distribute the local matrix`s nodes. 
      for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
            GM.add (local_dof_indices[i],local_dof_indices[j],lM(i,j));
    }

  }


  public:
    AssemblyOrthoMeshSystem(){}
    ~AssemblyOrthoMeshSystem(){}
  
  void buildGlobalMatrix(Matriz& GM,OrthoMesh &mesh, EqVarType& eqvar,const bool checkZero=false,double tolerance =0.0)
  {


    unsigned i,j;
    const unsigned int   dofs_per_cell = eqvar.n_dofs_per_cell();
    const unsigned n_quadrature_points = eqvar.n_quadrature_points();

    //Inicializar matriz global e local.
    GM = 0;
    m_Ml.reinit(dofs_per_cell,dofs_per_cell);
   
    

    //Obter referencia de alguns campos de EqVarT
 

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

    for (; cell!=endc; cell++)
    {
      eqvar.setCell(cell);
      eqvar.changedCellEvent();
      //Build local matrix.
      m_Ml=0.0;
      for (Index qPoint=0; qPoint<n_quadrature_points; qPoint++)
      {
        for (i=0; i<dofs_per_cell; ++i)
        {
          for (j=0; j<dofs_per_cell; ++j)
            m_Ml(i,j) += eqvar.localEqVar(i,j,qPoint);
        }
      }
      //Distribua os nos na matriz global
      this->distribute_local_to_global(GM,m_Ml,eqvar.local_to_global_dofs_map(),checkZero,tolerance);
    }
    return;
  }
    
    
  template<class SparsePattern>
  void buildSparsityPattern(SparsePattern &pattern,OrthoMesh &mesh,EqVarType& eqvar)
  {
    unsigned i,j;

    //Initialize local matrix.
    m_Ml.reinit(eqvar.n_dofs_per_cell(),eqvar.n_dofs_per_cell());
   
    OrthoMesh::Cell_It cell = mesh.begin_cell(),
    endc = mesh.end_cell();

      //Build local matrix.
      m_Ml=0.0;
    unsigned n_dofs_per_cell=eqvar.n_dofs_per_cell();

    for (i=0; i<n_dofs_per_cell; ++i)
    {
      for (j=0; j<n_dofs_per_cell; ++j)
        m_Ml(i,j) += (eqvar.nullProduct(i,j)) ? 0 : 1;
    }


    for (; cell!=endc; cell++)
    {
      eqvar.setCell(cell);
      for (unsigned int i=0; i<n_dofs_per_cell; ++i)
        for (unsigned int j=0; j<n_dofs_per_cell; ++j)
        {
          if (m_Ml(i,j) == 0.0)
          {
            //printf("Not adding position (%d,%d) for value %g\n",local_dof_indices[i],local_dof_indices[j],m_Ml(i,j));
            continue;
          }
          else
            pattern.add (eqvar.local_to_global_dofs_map()[i],eqvar.local_to_global_dofs_map()[j]);
        }
    }
  }
  
    void  buildGlobalRhs(Vetor& vRhs,OrthoMesh &mesh,EqVarType& eqvar)
    {
      unsigned i;

      const unsigned int   n_dofs_per_cell = eqvar.n_dofs_per_cell();
      unsigned n_quadrature_points = eqvar.n_quadrature_points();


      assert(vRhs.size() == eqvar.n_dofs());
      vRhs = 0;



      m_Bl.reinit(n_dofs_per_cell);



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

      for (; cell!=endc; cell++)
      {
        eqvar.setCell(cell);
        m_Bl = 0;
        for (unsigned qPoint=0; qPoint<n_quadrature_points; qPoint++)
        {
          for (i=0; i<n_dofs_per_cell; ++i)
          {
            m_Bl(i) += eqvar.localRHS(i,qPoint);
          }
        }
        //Distribua os nos na matriz global
        for (unsigned int i=0; i<n_dofs_per_cell; ++i)
        {
          vRhs(eqvar.local_to_global_dofs_map()[i]) += m_Bl(i);
        }
      }
    }


  void    addNeummanCondition(const OrthoMesh &mesh,Vetor &vRHS,EqVarType& eqvar,Function3D &fN)
  {
    Point3D pt;
    const unsigned int   dofs_per_cell = eqvar.n_dofs_per_cell();//Number of dofs per cell
    const unsigned int n_face_quadrature_pts = eqvar.n_face_quadrature_points();//Number of QPoints per face
    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.setCell(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_no=0;face_no<OrthoMesh::FACES_PER_CELL;face_no++)
        {
          OrthoMesh::Face_It face = cell->face(static_cast<FaceDirection3D>(face_no));
          if (face->at_boundary()) //found face at boundary  
          {
            if (fN.isInDomain(face->barycenter()))
            {
              eqvar.setFace(face); //Set the face;
              //For all local dofs of the cell
              for (Index local_dof=0;local_dof<dofs_per_cell;local_dof++)
              {
                //get what is the dof number of the face with respect to its local number
                //in terms of the cell
                Index cell_dof = eqvar.local_to_global_dofs_map()[local_dof]; 
                double acum=0;
                for (unsigned fQPoint=0;fQPoint < n_face_quadrature_pts;fQPoint++)
                {
                  vRHS(cell_dof)+=eqvar.localNeummanCondition(face_no,local_dof,fQPoint,fN);
                  acum+=eqvar.localNeummanCondition(face_no,local_dof,fQPoint,fN);
                }
                printf("Acum: %g\n",acum);
              }
            }
          }
        }
      }
    }
  }


  void setBoundaryValue(const OrthoMesh &mesh,EqVarType& eqvar,Function3D &f,MapIntDouble &map_boundary)
  {
    Point3D pt;
      const unsigned int   dofs_per_face = eqvar.n_dofs_per_face();
    std::vector<unsigned> face_dofs(dofs_per_face);
  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.setCell(cell);//If it is, set the cell and find faces at boundary

        for (unsigned face_no=0;face_no<OrthoMesh::FACES_PER_CELL;face_no++)
        {
          
          if (cell->face(static_cast<FaceDirection3D>(face_no))->at_boundary())
          {
            if (f.isInDomain(cell->face(static_cast<FaceDirection3D>(face_no))->barycenter()))
            //For all local dofs at face
            for (Index local_face_dof=0;local_face_dof<dofs_per_face;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 = eqvar.local_face_to_local_cell_map(face_no,local_face_dof); 
              pt=eqvar.get_global_dof_position(cell_dof);
              map_boundary[eqvar.local_to_global_dofs_map()[cell_dof]]=f(pt,eqvar.local_dof_component(cell_dof));
              
            }
          }
        }
      }
    }
  }

  

  };


#endif

---