mappedfefacevalues.cpp

Go to the documentation of this file.

#include "mappedfefacevalues.h"
#include "unitcube.h"
#include "exception.h"
#include "numericmethods.h"


MappedFEFaceValues::MappedFEFaceValues(FiniteElementInterface &fe,Index faceId) 
  :MappedFEValuesBase(fe,fe.face_to_cell_local_map(faceId).size())
{
  _dir = faceId;//Store the face of the face values

  //Store the mapping between the local face dofs
  //and the local cell dofs. This is obtained
  //fron the finite element.
  _face_to_cell_map= fe.face_to_cell_local_map(faceId);

  UnitCube ref;
  ref.face_base_vectors(faceId,&e1,&e2);

  
}

void MappedFEFaceValues::setPoints(const ArrayOfVecDouble &v,const VecDouble &weights) 
{
  
  //The array of points must be 2D points
  assert(v.vecs_size() ==2);

  ArrayOfVecDouble pts;
  pts.reinit(v.size(),3);

  VecDoubleRef pt3D(3);
  VecDouble    pt2D(2);
  UnitCube refDomain;

  for (unsigned i=0;i<v.size();i++)
  {
    pts.getVecValues(i,&pt3D);
    v.copyVecValues(i,&pt2D);
    refDomain.project_2d_points_into_faces(pt2D,&pt3D,static_cast<FaceDirection3D>(_dir));
  }
  pts.print();
  MappedFEValuesBase::setPoints(pts,weights);
}




void MappedFEFaceValues::compute_values() 
{
  unsigned iEnd = _pts.size();
  VecDoubleRef pt(3);
  unsigned n_local_dofs = _face_to_cell_map.size();

  for (unsigned i=0;i<iEnd;i++) //for all points in _pts list
  {
    _pts.getVecValues(i,&pt); //Get the point
    for (unsigned j=0;j<n_local_dofs;j++) //for all dofs
    {
      //evaluat the base function j on the point pt
      //and store in the matrix this->_Values;
      _Values(i,j)=_fe.shape_value(_face_to_cell_map[j],pt);
    }
  }

}


void MappedFEFaceValues::compute_grads(MappingOrthoMesh &map) 
{
  throw new Exception("MappedFEFaceValues::compute_grads()  not implemented yet\n");
}



 void MappedFEFaceValues::compute_J(MappingOrthoMesh &map) 
{
  static Matrix GradT(3,3); //A matrix 3x3
  VecDoubleRef pt(3);
  static VecDouble aux1(3);
  static VecDouble aux2(3);

  for(unsigned i=0;i<_pts.size();i++)
  {
    //get pt i
    _pts.getVecValues(i,&pt); //Get the point

    //get GradT
    map.T_grad(pt,GradT);

    GradT.vmult(aux1,e1);
    GradT.vmult(aux2,e2);
    
    _J(i) = NumericMethods::abs_vector_product(aux1,aux2);
  }

}


void MappedFEFaceValues::compute_JxW(MappingOrthoMesh &map,const  VecDouble &weights) 
{
  assert(weights.size() == _pts.size());
  compute_J(map);
  _JxW.term_mult(_J,weights);
}


void MappedFEFaceValues::setQuadrature(const dealii::Quadrature<2> &quad) 
{
  //printf("Caralho\n");
  MappedFEFaceValues::setPoints(quad.get_points(),quad.get_weights());
}


 void MappedFEFaceValues::compute_JxW(MappingOrthoMesh &map) 
{
  compute_JxW(map,_weights);
}

---