assemblysystem.h

Go to the documentation of this file.

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



/*----------------- CLASS AssemblySystem --------------
  ----------------------------------------*/
  enum AssemblyOptimization {NO_OPTIMIZATION,EQUAL_LOCAL_SYSTEMS,EQUAL_TOPOLOGY};

template <class EqVarType,class Matriz = SparseMatrix<double>, class Vetor = VecDouble>
  class AssemblySystem : public AssemblySystemBase, EqVarTAccessor
  {



 public:


  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:
    AssemblySystem(){}
    ~AssemblySystem(){}

  void buildGlobalMatrix(Matriz& GM,DoFHandler<DIM> &dof,EqVarType& eqvar,unsigned nQPointsPerDim,AssemblyOptimization  opt=NO_OPTIMIZATION,const bool checkZero=false,double tolerance =0.0)
  {

      setValuesAtElements(eqvar,true);
      setValuesAtFaces(eqvar,false);

      unsigned i,j;
      QGauss<DIM>  quadrature_formula(nQPointsPerDim);
      FEValues<DIM> fe_values (dof.get_fe(), quadrature_formula, 
                               eqvar.getUpdateFlags());
      const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
      const unsigned n_quadrature_points = quadrature_formula.n_quadrature_points;

      //Inicializar matriz global e local.
      GM = 0;
      m_Ml.reinit(dof.get_fe().dofs_per_cell,dof.get_fe().dofs_per_cell);
   
    
      //Configurar fe values da equacao variacional
      eqvar.setFeValue(&fe_values);

      //Obter referencia de alguns campos de EqVarT
      unsigned& qPoint                             = getRefQPoint(eqvar); 
      std::vector<unsigned int>& local_dof_indices = getRefDofIndices(eqvar);

      //Alocar vetor com dofs.
      local_dof_indices.resize(dof.get_fe().dofs_per_cell);

      DoFHandler<DIM>::active_cell_iterator cell = dof.begin_active(),
        endc = dof.end();

      //Temos 3 opcoes de montagem das matrizes.
      //EQUAL_LOCAL_SYSTEMS
      if (opt == EQUAL_LOCAL_SYSTEMS)
      {
        fe_values.reinit(cell);
        cell->get_dof_indices(local_dof_indices);
        setCell(eqvar,cell);
        eqvar.changedCellEvent();

        //Construir Matriz Local.
        m_Ml=0.0;
        for (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);
          }
        }
        for (; cell!=endc; ++cell)
        {
          cell->get_dof_indices(local_dof_indices);
          //Distribua os nos na matriz global
          this->distribute_local_to_global(GM,m_Ml,local_dof_indices,checkZero,tolerance);
        }
        return;
      }
      if (opt == EQUAL_TOPOLOGY)
      {
        fe_values.reinit(cell);
        for (; cell!=endc; ++cell)
        {
          cell->get_dof_indices(local_dof_indices);
          setCell(eqvar,cell);
          eqvar.changedCellEvent();
          //Construir Matriz Local.
          m_Ml=0.0;
          for (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);
            }
          }
          //Distribua os nos na matriz global
          this->distribute_local_to_global(GM,m_Ml,local_dof_indices,checkZero,tolerance);
        }
        return;
      }
      //Ok we dont have any processing savings
      //Do the usual assembly
    
      for (; cell!=endc; ++cell)
      {
        fe_values.reinit(cell);
        cell->get_dof_indices(local_dof_indices);
        setCell(eqvar,cell);
        eqvar.changedCellEvent();
        //Construir Matriz Local.
        m_Ml=0.0;
        for (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);
          }
        }

        //Distribua os nos na matriz global
        this->distribute_local_to_global(GM,m_Ml,local_dof_indices,checkZero,tolerance);
      }
    }

  template<class SparsePattern>
  void buildSparsityPattern(SparsePattern &pattern,DoFHandler<DIM> &dof,EqVarType& eqvar,double tolerance=0.0)
  {
    setValuesAtElements(eqvar,true);
    setValuesAtFaces(eqvar,false);

    unsigned i,j;

    QGauss<DIM>  quadrature_formula(1);
    FEValues<DIM> fe_values (dof.get_fe(), quadrature_formula, 
                             eqvar.getUpdateFlags());
    const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
    const unsigned n_quadrature_points = quadrature_formula.n_quadrature_points;

    //Initialize local matrix.
    m_Ml.reinit(dof.get_fe().dofs_per_cell,dof.get_fe().dofs_per_cell);
   
    
    //Set fe values of the variational equation
    eqvar.setFeValue(&fe_values);

    //Get the reference of some fields of EqVarT
    unsigned& qPoint                             = getRefQPoint(eqvar); 
    std::vector<unsigned int>& local_dof_indices = getRefDofIndices(eqvar);

    //Allocate the vector of dofs indices.
    local_dof_indices.resize(dof.get_fe().dofs_per_cell);

    DoFHandler<DIM>::active_cell_iterator cell = dof.begin_active(),
    endc = dof.end();

    for (; cell!=endc; ++cell)
    {
      fe_values.reinit(cell);
      cell->get_dof_indices(local_dof_indices);
      setCell(eqvar,cell);
      eqvar.changedCellEvent();

      //Construir Matriz Local.
      m_Ml=0.0;
      for (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);
        }
      }
      for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
        {
          if (fabs(m_Ml(i,j)) < tolerance)
          {
            //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 (local_dof_indices[i],local_dof_indices[j]);
        }
    }
  }

    void  buildGlobalRhs(Vetor& vRhs,DoFHandler<DIM> &dof,EqVarType& eqvar, unsigned nQPointsPerDim)
    {
      unsigned i;

      setValuesAtElements(eqvar,true);
      setValuesAtFaces(eqvar,false);

      QGauss<DIM>  quadrature_formula(nQPointsPerDim);
      FEValues<DIM> fe_values (dof.get_fe(), quadrature_formula, 
                               eqvar.getUpdateFlags());
      const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
      unsigned n_quadrature_points = quadrature_formula.n_quadrature_points;


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

      //Obter referencia dos campos
      std::vector<unsigned int>& local_dof_indices = getRefDofIndices(eqvar);
      unsigned &qPoint = getRefQPoint(eqvar);


      eqvar.setFeValue(&fe_values);
      m_Bl.reinit(dofs_per_cell);
      local_dof_indices.resize(dofs_per_cell);



      DoFHandler<DIM>::active_cell_iterator cell = dof.begin_active(),
        endc = dof.end();

      for (; cell!=endc; ++cell)
      {
        fe_values.reinit(cell);
        cell->get_dof_indices(local_dof_indices);
        setCell(eqvar,cell);
        eqvar.changedCellRHSEvent();
        m_Bl = 0;
        for (qPoint=0; qPoint<n_quadrature_points; ++qPoint)
        {
          eqvar.changedQPointRHSEvent();
          for (i=0; i<dofs_per_cell; ++i)
          {
            m_Bl(i) += eqvar.localRhs(i);
          }
        }

        //Distribua os nos na matriz global
        for (unsigned int i=0; i<dofs_per_cell; ++i)
        {
          vRhs(local_dof_indices[i]) += m_Bl(i);
        }
      }
    }


  
    
    void    addNeummanCondition(const DoFHandler<DIM> &dof,Vetor& vRhs,EqVarType& eqvar,unsigned nQPointsPerDim=2,unsigned boundary_value=0)
    {
      setValuesAtElements(eqvar,false);
      setValuesAtFaces(eqvar,true);

      QGauss<DIM-1> face_quadrature(nQPointsPerDim);
      FEFaceValues<DIM> fe_face_values(dof.get_fe(),face_quadrature,
                                       eqvar.getUpdateFlags()|update_normal_vectors|update_q_points ); 
      const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
      unsigned n_quadrature_points = face_quadrature.n_quadrature_points;


      //Obter referencia de campos de eqvar.
      unsigned &qPointFace = getRefQPointFace(eqvar);
      std::vector<unsigned int>& local_dof_indices = getRefDofIndices(eqvar);

    
      eqvar.setFaceValues(&fe_face_values);
      eqvar.setFeValue(NULL);
      local_dof_indices.resize(dofs_per_cell);
    

      DoFHandler<DIM>::active_cell_iterator cell = dof.begin_active(),
        endc = dof.end();
      for (; cell!=endc; ++cell)
      {
        if (cell->at_boundary()) //Celula na fronteira ? 
        {
          cell->get_dof_indices(local_dof_indices);
          setCell(eqvar,cell);
          for (unsigned face=0; face < GeometryInfo<DIM>::faces_per_cell; ++face)
            if (cell->face(face)->boundary_indicator() == boundary_value)
            {
              fe_face_values.reinit(cell,face);
              for (qPointFace=0;qPointFace<n_quadrature_points;qPointFace++)
              {
                for (unsigned i=0;i<dofs_per_cell;i++)
                {
                  vRhs(local_dof_indices[i])+=eqvar.localNeumannBoundary(i);
                }
              }
            }
        }
      }
    }




  void setBoundaryValue(DoFHandler<DIM> &dof,Function3D &f,MapIntDouble &map_boundary)
  {
    assert(dof.get_fe().n_components()== f.n_components());
    const FiniteElement<DIM> &fe           = dof.get_fe();
    std::vector<unsigned int> face_dofs (fe.dofs_per_face,
                                         DoFHandler<DIM>::invalid_dof_index);
    std::vector<Point<DIM> >  dof_locations (face_dofs.size(), Point<DIM>());
    std::vector<double>          dof_values_scalar (fe.dofs_per_face);
    std::vector<Point<DIM-1> > unit_support_points = fe.get_unit_face_support_points();
    assert(unit_support_points.size() !=0);   //Tem que ser os elementos que estou habituado
    assert(fe.is_primitive());
    Quadrature<DIM-1> aux_quad (unit_support_points);
    FEFaceValues<DIM> fe_values (fe, aux_quad, update_q_points);

    DoFHandler<DIM>::active_cell_iterator cell = dof.begin_active(),
    endc = dof.end();

    if (fe.n_components() == 1)
    {
      for (; cell!=endc; ++cell)
        for (unsigned int face_no = 0; face_no < GeometryInfo<DIM>::faces_per_cell;++face_no)
        {
          DoFHandler<DIM>::face_iterator face = cell->face(face_no);
          if (face->at_boundary())
          {
            face->get_dof_indices (face_dofs);
            fe_values.reinit(cell, face_no);
            const std::vector<Point<DIM> > &dof_locations = fe_values.get_quadrature_points ();
            for (unsigned i=0;i<dof_locations.size();i++)
            {
              if (f.isInDomain(dof_locations[i]))
              {
                map_boundary[face_dofs[i]] = f(dof_locations[i]);
              }
            }
          }
        }
    }
    else //Mais de uma componente, temos mais de um grau de liberdade
    {
      for (; cell!=endc; ++cell)
        for (unsigned int face_no = 0; face_no < GeometryInfo<DIM>::faces_per_cell;++face_no)
        {
          DoFHandler<DIM>::face_iterator face = cell->face(face_no);
          if (face->at_boundary())
          {
            face->get_dof_indices (face_dofs);
            fe_values.reinit(cell, face_no);
            const std::vector<Point<DIM> > &dof_locations = fe_values.get_quadrature_points ();
            for (unsigned i=0;i<dof_locations.size();i++)
            {
              unsigned component = fe.face_system_to_component_index(i).first;
              if (f.isInDomain(dof_locations[i],component))
              {
                map_boundary[face_dofs[i]] = f(dof_locations[i],component);
              }
            }
          }
        }
    }
  }


  

  };


#endif

---