dealfewrapper.cpp

Go to the documentation of this file.

#include "dealfewrapper.h"
#include <grid/grid_generator.h>

/*
void DealFEWrapper::updateFE_Values(OrthoMesh &mesh,FEValues<3> &fe_values) 
{
  Triangulation<3> *tria=new Triangulation<3>();
  Point3D P1(0,0,0);
  GridGenerator::hyper_rectangle(*tria,P1,mesh.getDX());
  DoFHandler<3> *dof = new DoFHandler<3>(*tria);
  dof->distribute_dofs(fe_values.get_fe());
  DoFHandler<DIM>::active_cell_iterator *cell = new DoFHandler<DIM>::active_cell_iterator(dof->begin_active());
  fe_values.reinit(*cell);
  m_pCell=NULL;

}

*/
/*
void DealFEWrapper::updateFE_FaceValues(OrthoMesh &mesh,FaceValuesVec &vFaceValues) 
{
  assert(vFaceValues.size() == 6);
  Triangulation<3> tria;
  Point3D P1(0,0,0);
  GridGenerator::hyper_rectangle(tria,P1,mesh.getDX());
  DoFHandler<3> dof(tria);
  dof.distribute_dofs(vFaceValues[0]->get_fe());
  DoFHandler<DIM>::active_cell_iterator cell = dof.begin_active(),
    endc = dof.end();

  vFaceValues[LEFT_FACE]->reinit(cell,LEFT_FACE);
  vFaceValues[RIGHT_FACE]->reinit(cell,LEFT_FACE);
  vFaceValues[UP_FACE]->reinit(cell,LEFT_FACE);
  vFaceValues[BOTTOM_FACE]->reinit(cell,LEFT_FACE);
  vFaceValues[FRONT_FACE]->reinit(cell,LEFT_FACE);
  vFaceValues[BACK_FACE]->reinit(cell,LEFT_FACE);
}
*/

 DealFEWrapper::DealFEWrapper() 
{
  m_fev=NULL;
  m_pTria=NULL;
  m_pDof=NULL;
  m_fe=NULL;
  m_pMesh=NULL;
}


void  DealFEWrapper::initFEValues(OrthoMesh &mesh,const FiniteElement<3> &fe, const Quadrature<3> &quadrature, const Quadrature<2> &face_quadrature,const UpdateFlags update_flags,const UpdateFlags face_update_flags) 
{
  m_pMesh=&mesh;
  m_fe=&fe;
  m_pTria= new Triangulation<3>();
  Point3D P1(0,0,0);
  GridGenerator::hyper_rectangle(*m_pTria,P1,mesh.getDX());
  m_pDof = new DoFHandler<3>(*m_pTria);
  m_pDof->distribute_dofs(fe);
  m_pDoFCell = new DoFHandler<DIM>::active_cell_iterator(m_pDof->begin_active());
  m_fev = new FEValues<3>(fe,quadrature,update_flags);
  m_fev->reinit(*m_pDoFCell);

  
  //INIT VECTOR OF FACE VALUES, One for each face
  m_facesValues.resize(6,NULL);
  m_facesValues[UP_FACE]    =new FEFaceValues<3>(fe,face_quadrature,face_update_flags);
  m_facesValues[BOTTOM_FACE]=new FEFaceValues<3>(fe,face_quadrature,face_update_flags);
  m_facesValues[LEFT_FACE]  =new FEFaceValues<3>(fe,face_quadrature,face_update_flags);
  m_facesValues[RIGHT_FACE] =new FEFaceValues<3>(fe,face_quadrature,face_update_flags);
  m_facesValues[BACK_FACE]  =new FEFaceValues<3>(fe,face_quadrature,face_update_flags);
  m_facesValues[FRONT_FACE] =new FEFaceValues<3>(fe,face_quadrature,face_update_flags);

  //Set the FEFaceValues
  m_facesValues[LEFT_FACE  ]->reinit(*m_pDoFCell,LEFT_FACE  );
  m_facesValues[RIGHT_FACE ]->reinit(*m_pDoFCell,RIGHT_FACE );
  m_facesValues[BOTTOM_FACE]->reinit(*m_pDoFCell,BOTTOM_FACE);
  m_facesValues[UP_FACE    ]->reinit(*m_pDoFCell,UP_FACE    );
  m_facesValues[FRONT_FACE ]->reinit(*m_pDoFCell,FRONT_FACE );
  m_facesValues[BACK_FACE  ]->reinit(*m_pDoFCell,BACK_FACE) ;
}


 void DealFEWrapper::clear() 
{
  assert(m_pDof);
  m_pDof->clear();
  delete m_pDoFCell;
  delete m_pDof;
  delete m_pTria;
  delete m_fev;
  for (unsigned i=0;i<6;i++)
    delete m_facesValues[i];

  m_pDoFCell=NULL;
  m_pDof=NULL;
  m_pTria=NULL;
  m_fev=NULL;
}

DealFEWrapper::~DealFEWrapper() 
{
  if (m_pDof)
  {
    this->clear();
  }
}
void DealFEWrapper::setCell(OrthoMesh::Cell_It &cell) 
{
  m_pCell=&cell;
}




Point3D DealFEWrapper::get_qpoint(unsigned qPoint) 
{
  assert(m_fev);
  Point3D pt = m_fev->get_quadrature_points()[qPoint];
  pt+=getCell()->vertex(VERTEX_000);
  return pt;
}

Point3D DealFEWrapper::get_face_qpoint(unsigned qPoint) 
{
  assert(m_facesValues[m_faceDir]);
  Point3D pt = m_facesValues[m_faceDir]->get_quadrature_points()[qPoint];
  pt+=(*m_pCell)->vertex(VERTEX_000);
  return pt;
}



void DealFEWrapper::setFace(OrthoMesh::Face_It &face) 
{
  m_faceDir = face->getFaceDir();//thids method only makes sense for faces at boundary
  
}




void DealFEWrapper::printDealFEInfo() 
{
  assert(m_fe);
  printf("====Printing Finete Element====\n");
  printf("DOFS PER CELL %d\n",m_fe->n_dofs_per_cell());
  printf("DOFS PER FACE %d\n",m_fe->n_dofs_per_face());
  printf("DOFS PER HEX %d\n",m_fe->n_dofs_per_hex());
  printf("COMPONENTS %d\n",m_fe->n_components());

  printf("LOCAL FACE LOCAL CELL MAPPING\n\n");
  for (unsigned face=0;face < GeometryInfo<3>::faces_per_cell;face++)
  {
    printf("FACE %d\n",face);
    for (unsigned i=0;i<m_fe->n_dofs_per_face();i++)
    {
      printf("\tDOF %d = %d\n",i,local_face_to_local_cell_map(face,i));
    }
  }
    printf("DOF TO COMPONENT MAPPING\n\n");
  for (unsigned dof=0;dof<n_dofs_per_cell();dof++)
  {
    printf("\tDof %d is Component %d\n",dof,local_dof_component(dof));
  }
    
}

---