mappedfevalues.cpp

Go to the documentation of this file.

#include "mappedfevalues.h"
#include "numericmethods.h"


void MappedFEValues::setPoints(const dealii::Quadrature<3> &quad) 
{
  MappedFEValuesBase::setPoints(quad.get_points(),quad.get_weights());
}


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

  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
    {
      //evaluates the base function j on the point pt
      //and store in the matrix this->_Values;
      _Values(i,j)=_fe.shape_value(j,pt);
    }
  }
}


void MappedFEValues::compute_grads(MappingOrthoMesh &map)  
{
  unsigned iEnd = _pts.size();
  unsigned n_local_dofs = _fe.n_local_dofs();
  VecDoubleRef pt(3);
  static VecDouble vGradRef(3);
  Matrix invTGrad(3,3);
  for (unsigned iPt=0;iPt<iEnd;iPt++) //for all points in _pts list
  {
    _pts.getVecValues(iPt,&pt); //Get the point
    map.Inv_T_Grad(pt,invTGrad);

    for (unsigned j=0;j<n_local_dofs;j++) //for all dofs
    {
      _fe.shape_grad(j,pt,vGradRef);
      invTGrad.Tvmult(_Grad[iPt*m_n_dofs + j],vGradRef);
    }
  }
}




void MappedFEValues::compute_J(MappingOrthoMesh &map) 
{
  static Matrix GradT(3,3); //A matrix 3x3
  VecDoubleRef pt(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);
    _J(i)=NumericMethods::det3x3(GradT);
  }
}



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

void MappedFEValues::compute_JxW(MappingOrthoMesh &map){
  assert(_weights.size() == _pts.size());
  compute_J(map);
  _JxW.term_mult(_J,_weights);
}



MappedFEValues::MappedFEValues(FiniteElementInterface &fe) 
  :MappedFEValuesBase(fe,fe.n_local_dofs())
 {

 }


void MappedFEValues::compute(MappingOrthoMesh &map,unsigned compute_flags) 
{
  m_compute_flags=compute_flags;
  if (compute_flags | COMPUTE_VALUES)
    compute_values();
  if  (compute_flags | COMPUTE_GRADS)
    compute_grads(map);
  if  (compute_flags | COMPUTE_PTS)
    compute_points(map);
  if  (compute_flags | COMPUTE_JxW)
    compute_JxW(map,_weights);
  else if (compute_flags | COMPUTE_J)
    compute_J(map);

}


void MappedFEValues::setPoints(const ArrayOfVecDouble &pts,const VecDouble &weights) 
{
  MappedFEValuesBase::setPoints(pts,weights);
}

---