dealbasedofs.cpp

Go to the documentation of this file.

#include "dealbasedofs.h"
#include "numericmethods.h"
#include <vector>
#include <fe/fe_base.h>
#include <grid/grid_generator.h>
DoFCell DealBaseDoFs::getDoFCellAtPoint(const DoFHandler<DIM>& dof,const Point<DIM> &p) 
{
  DoFCell cell = dof.begin_active();
  DoFCell endc = dof.end();
  
  for (;cell!= endc;cell++)
  {
    if (NumericMethods::isInCube(p,cell->vertex(VERTEX_000),cell->vertex(VERTEX_111)))
    {
      return cell;
    }
  }
  return DoFCell();

}

/***************************************************/
DoFCell DealBaseDoFs::getAdjDoFCell(DoFCell cell,CellDirection3D dir) 
{
  if (cell->face(dir)->at_boundary())
    return DoFCell();
  else
    return cell->neighbor(dir);
}



/***************************************************/
double DealBaseDoFs::getDoFVerticeValue(DoFCell cell,VertexDirection3D dir,unsigned deg,const VecDouble &vSol) 
{
  assert(cell->vertex_dof_index(dir,deg) < vSol.size());
  return vSol(cell->vertex_dof_index(dir,deg));
}



void DealBaseDoFs::printCellDofsInformation(const DoFHandler<2>& dofs,std::ostream &out) 
 {
  const  DoFHandler<2> &dof = dofs;
  DoFHandler<2>::active_cell_iterator cell = dof.begin_active(),
                                      endc = dof.end();
  char strOut[500];
  
  for (;cell!=endc; ++cell)
    {
      sprintf(strOut,"\n==================\nCell %d:\nDoFs: %d\n\nVertices:\n",cell->index(),dof.get_fe().dofs_per_cell);
      out << strOut;
      for (unsigned i=0;i<GeometryInfo<DIM>::vertices_per_cell;i++)
      {
        sprintf(strOut,"\n%d) Pt: <%6g,%6g>, Index: %d, Dofs Index:",i,cell->vertex(i)(0),cell->vertex(i)(1),cell->vertex_index(i));
        out << strOut;
        for (unsigned j=0;j<dof.get_fe().dofs_per_vertex;j++)
        {
          out << cell->vertex_dof_index(i,j) << " ";
        }  
      }
    }
}
void DealBaseDoFs::printCellDofsInformation(const DoFHandler<3>& dofs) 
 {
  const  DoFHandler<DIM> &dof = dofs;
  DoFHandler<DIM>::active_cell_iterator cell = dof.begin_active(),
                                      endc = dof.end();
  int dofs_per_face = dof.get_fe().n_dofs_per_face();
  vector<unsigned> dofs_face(dofs_per_face);
  
  for (;cell!=endc; ++cell)
    {
      printf("\n==================\nCell %d at <%3g,%3g,%3g>\n",cell->index(),
             cell->barycenter()[0],
             cell->barycenter()[1],
             cell->barycenter()[2]);
      if (dofs_per_face > 0)
      {
        printf("Dof at Faces:\n");
        
        for (unsigned i=0;i<GeometryInfo<DIM>::faces_per_cell;i++)
        {
          
          printf("%d)  Index: %d Dofs:",i,cell->face(i)->index());
          cell->face(i)->get_dof_indices(dofs_face);
          for (int j=0;j<dofs_per_face;j++)
          {
            printf("%d ",dofs_face[j]);
          }
          printf("\n");
        }
      }
    }
}



void DealBaseDoFs::getVerticeValuesFromDofs(DoFHandler<DIM> &dofs,const VecDouble &vSol,VecDouble &vVerticesValues,unsigned deg) 
{


  assert(dofs.get_fe().is_primitive());
  assert(vSol.size() == dofs.n_dofs());
  assert(vVerticesValues.size() == dofs.get_tria().n_vertices());


  const unsigned int   dofs_per_cell = dofs.get_fe().dofs_per_cell;
  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
  DoFHandler<3>::active_cell_iterator cell = dofs.begin_active(),
                                      endc = dofs.end();
  for (; cell!=endc; ++cell)
  {
    cell->get_dof_indices (local_dof_indices);
    for (unsigned i=0;i<GeometryInfo<2>::vertices_per_cell;i++)
    {
      vVerticesValues(cell->vertex_index(i))= vSol(cell->vertex_dof_index(i,deg));
    }
  }            
  return;
}

void DealBaseDoFs::printFE(FiniteElement<3> &fe) 
{
  cout << "Finite Element Name: " << fe.get_name() << std::endl;
  cout << "Is primitive: " << fe.is_primitive() << std::endl;
  cout << "Num Base Elements: " << fe.n_base_elements() << std::endl;
  cout << "Num components: " << fe.n_components() << std::endl;
  //  cout << "Num DOFS per object: " << fe.n_dofs_per_object() << std::endl;
  cout << "DOFS per vertex: " << fe.n_dofs_per_vertex() << std::endl;
  cout << "DOFS per line: " << fe.n_dofs_per_line() << std::endl;
  cout << "DOFS per quad: " << fe.n_dofs_per_quad() << std::endl;
  cout << "DOFS per hex: " << fe.n_dofs_per_hex() << std::endl;
  cout << "DOFS per face: " << fe.n_dofs_per_face() << std::endl;
  cout << "DOFS per cell: " << fe.n_dofs_per_cell() << std::endl;
  
  cout << "Unit Support Points: " <<  std::endl;

  const std::vector<Point<3> > &vPts = fe.get_unit_support_points();
  for (unsigned i =0; i < vPts.size(); i++)
    printf("%d) %g, %g, %g\n", i, vPts[i](0), vPts[i](1) , vPts[i](2));
  if (vPts.size() == 0)
    printf("None\n");




  //Now to see the values of the shape functions we need a triangulation
  //with one element in it.
  Triangulation<3> tria;
  Point<3> P1(0,0,0);
  Point<3> P2(1,1,1);
  std::vector<unsigned> v(3);
  v[0]=2;
  v[1]=2;
  v[2]=2;
  GridGenerator::subdivided_hyper_rectangle(tria,v,P1,P2);


  //We need a dofHandler too
  DoFHandler<3> dof(tria);
  dof.distribute_dofs(fe);
  DoFHandler<3>::active_cell_iterator cell = dof.begin_active();

  //Define the points where we will extract the values of the shape function.
  //These points are stored in the vector v and passed as argumento to the Quadrature class
  //who is responsible for the gaussian integration.
  std::vector<Point<3> > vQuad;
  vQuad.push_back(Point<3>(0,0,0));
  vQuad.push_back(Point<3>(1,0,0));
  vQuad.push_back(Point<3>(0,1,0));
  vQuad.push_back(Point<3>(1,1,0));
  vQuad.push_back(Point<3>(0,0,1));
  vQuad.push_back(Point<3>(1,0,1));
  vQuad.push_back(Point<3>(0,1,1));
  vQuad.push_back(Point<3>(1,1,1));

  std::vector<double> vWeights(8,1);
  Quadrature<3> quad(vQuad,vWeights);
  FEValues<3> feValues(fe,quad,UpdateFlags(update_values | update_JxW_values));
  feValues.reinit(cell);


  cout << "Shape Values At Vertex: " <<std::endl;
  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
  {
    cout << "Shape Value: " << i <<"\n";
    for (unsigned int q_point=0; q_point < quad.size(); ++q_point)
    {
      cout << vQuad[q_point] << ":\t";
      for (unsigned int component=0; component < fe.n_components(); component++)
      {
        //cout << fe.shape_value_component(i, quad.point(q_point), component) << "\t";
        cout << feValues.shape_value_component(i,q_point,component);
        //cout << feValues.shape_value(i,q_point);

      }
      cout << std::endl;
    }
  }
  
  //  for (unsigned i=0;i<


}


void DealBaseDoFs::getCentralValuesFromDofs(DoFHandler<3> &dof,const VecDouble &sol, VecDouble &cValues,unsigned component) 
{
  assert(sol.size() == dof.n_dofs());
  assert(cValues.size() == dof.get_tria().n_active_cells());
  assert(component < dof.get_fe().n_components());

  QMidpoint<DIM> qMidPoint;
  FEValues<DIM> fe_values(dof.get_fe(),qMidPoint,update_values);

  const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
 
 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
 DoFHandler<3>::active_cell_iterator cell = dof.begin_active(),
   endc = dof.end();

 double acum;
 for (; cell!=endc; ++cell)
 {
   fe_values.reinit (cell);
   acum =0.0;
   cell->get_dof_indices(local_dof_indices);
   for (unsigned int i=0; i<dofs_per_cell; ++i)
   {
     acum+= sol(local_dof_indices[i])*fe_values.shape_value_component(i,0,component);
   }
   cValues(cell->index())=acum;
 }
  
}



















       

---