dealorthomesh.cpp

Go to the documentation of this file.

#include "dealorthomesh.h"
#include "numericmethods.h"


 DealOrthoMesh::DealOrthoMesh(Triangulation<3> &tria)
   :DealBase(tria)
{
  //Since the mesh is uniform, all cells have the same dimensions, so store these
  //values for easy retieval.
  Cell cell = getTriangulation().begin_active();
  m_dX = DealBase::getDX(cell);
  m_dY = DealBase::getDY(cell);
  m_dZ = DealBase::getDZ(cell);
  m_Vol = m_dX * m_dY * m_dZ;
  m_FacesArea[NORMAL_X] = m_dY*m_dZ;
  m_FacesArea[NORMAL_Y] = m_dX*m_dZ;
  m_FacesArea[NORMAL_Z] = m_dX*m_dY;
  P = cell->vertex(VERTEX_000);

  
  m_FacesAreaFromCellDirection[LEFT_FACE]=m_FacesArea[NORMAL_X];
  m_FacesAreaFromCellDirection[RIGHT_FACE]=m_FacesArea[NORMAL_X];
  m_FacesAreaFromCellDirection[FRONT_FACE]=m_FacesArea[NORMAL_Y];
  m_FacesAreaFromCellDirection[BACK_FACE]=m_FacesArea[NORMAL_Y];
  m_FacesAreaFromCellDirection[BOTTOM_FACE]=m_FacesArea[NORMAL_Z];
  m_FacesAreaFromCellDirection[UP_FACE]=m_FacesArea[NORMAL_Z];
  
  m_FacesAreaPerVolume[NORMAL_X] = 1.0/m_dX;
  m_FacesAreaPerVolume[NORMAL_Y] = 1.0/m_dY;
  m_FacesAreaPerVolume[NORMAL_Z] = 1.0/m_dZ;
  //build the information about the faces.
  buildFaceInformation(tria);

  //calculate number of faces at boundary;
  Face_It end1 = end();
  n_boundary_faces=0;
  for(Face_It face = begin();face != end1;face++)
  {
    if (face->at_boundary())
      n_boundary_faces++;
  }

  int nX=0,nY=0,nZ=0;
  for (Cell cell = beginCell();cell.state()==IteratorState::valid;cell=cell->neighbor(RIGHT_CELL))
    nX++;
  for (Cell cell = beginCell();cell.state()==IteratorState::valid;cell=cell->neighbor(BACK_CELL))
    nY++;
  for (Cell cell = beginCell();cell.state()==IteratorState::valid;cell=cell->neighbor(UP_CELL))
    nZ++;
  nElems[0]=nX;
  nElems[1]=nY;
  nElems[2]=nZ;
  Q[0] = P[0] + nElems[0]*getDX();
  Q[1] = P[1] + nElems[1]*getDY();
  Q[2] = P[2] + nElems[2]*getDZ();
  
  
  

}

 DealOrthoMesh::~DealOrthoMesh() 
{

}



void DealOrthoMesh::buildFaceInformation(Triangulation<3> &tria)
{
  m_faces.clear();
  m_faces.resize(tria.n_active_faces());
  //For each cell.
  for (Triangulation<3>::active_cell_iterator cell = tria.begin_active(); cell != tria.end(); cell++)
  {
    m_faces[cell->face_index(LEFT_FACE  )].setNegNormalCell(cell->index());
    m_faces[cell->face_index(LEFT_FACE  )].setNormalDirection(NORMAL_X);

    m_faces[cell->face_index(RIGHT_FACE )].setPosNormalCell(cell->index());
    m_faces[cell->face_index(RIGHT_FACE  )].setNormalDirection(NORMAL_X);

    m_faces[cell->face_index(FRONT_FACE )].setNegNormalCell(cell->index());
    m_faces[cell->face_index(FRONT_FACE )].setNormalDirection(NORMAL_Y);

    m_faces[cell->face_index(BACK_FACE  )].setPosNormalCell(cell->index());
    m_faces[cell->face_index(BACK_FACE  )].setNormalDirection(NORMAL_Y);

    m_faces[cell->face_index(BOTTOM_FACE)].setNegNormalCell(cell->index());
    m_faces[cell->face_index(BOTTOM_FACE)].setNormalDirection(NORMAL_Z);

    m_faces[cell->face_index(UP_FACE    )].setPosNormalCell(cell->index());
    m_faces[cell->face_index(UP_FACE    )].setNormalDirection(NORMAL_Z);

  }
  
}


double FaceInf::normal_multiply(double *v) const
{
        return v[this->m_normal];

}


double DealOrthoMesh::faceArea(const Face_It & face)
{
        return m_FacesArea[face->getNormalOrientation()];
}


double DealOrthoMesh::faceAreaPerCellVol(const Face_It &face) 
{
  return m_FacesAreaPerVolume[face->getNormalOrientation()];
}






double DealOrthoMesh::faceValue( const Face_It & face,const VecDouble & fValues)
{
  assert( (unsigned) (face-begin()) < fValues.size());
  return fValues(face-begin());
}








void DealOrthoMesh::projectCentralValuesAtVertices(const VecDouble &cValues,VecDouble &vValues) 
{
  //Check the size of the vectors;
  assert(vValues.size() == numVertices());
  assert(cValues.size() == numCells());

  //This vector is a counter of the number of cells contributing to each vertice
  //It is initialized as zero.
  std::vector<unsigned> nCellsPerVertice(numVertices(),0);

  //Initialize the values of the vector vValues.
  vValues=0.0;

  Cell endc = endCell();
  for (Cell cell = beginCell();cell != endc;cell++)
  {
    vValues(cell->vertex_index(VERTEX_000)) += cValues(cell->index());
    vValues(cell->vertex_index(VERTEX_001)) += cValues(cell->index());
    vValues(cell->vertex_index(VERTEX_010)) += cValues(cell->index());
    vValues(cell->vertex_index(VERTEX_011)) += cValues(cell->index());
    vValues(cell->vertex_index(VERTEX_100)) += cValues(cell->index());
    vValues(cell->vertex_index(VERTEX_101)) += cValues(cell->index());
    vValues(cell->vertex_index(VERTEX_110)) += cValues(cell->index());
    vValues(cell->vertex_index(VERTEX_111)) += cValues(cell->index());

    nCellsPerVertice[cell->vertex_index(VERTEX_000)]++;
    nCellsPerVertice[cell->vertex_index(VERTEX_001)]++;
    nCellsPerVertice[cell->vertex_index(VERTEX_010)]++;
    nCellsPerVertice[cell->vertex_index(VERTEX_011)]++;
    nCellsPerVertice[cell->vertex_index(VERTEX_100)]++;
    nCellsPerVertice[cell->vertex_index(VERTEX_101)]++;
    nCellsPerVertice[cell->vertex_index(VERTEX_110)]++;
    nCellsPerVertice[cell->vertex_index(VERTEX_111)]++;

  }

  //Now in vValues we have the sum of all sunrounding cells
  //and in nCellsPerVertice the number of cells sunrounding the vertices.
  //Rember that these number of cells sunrounding a vertice can differ
  //for vertices at boundary.
  for (unsigned i =0;i<vValues.size();i++)
      vValues(i)/=nCellsPerVertice[i];
}








    
    

    

double DealOrthoMesh::getL2NormAtCells(const VecDouble &cValues) const
{
  assert(cValues.size() == numCells());
  double cellArea = getCellMeasure();
  Cell endc = endCell();
  double acum=0;
  for (Cell cell=beginCell();cell!=endc;cell++)
  {
    acum+=fabs(getCentralValue(cell,cValues));
  }
  return acum*cellArea;
}






double DealOrthoMesh::getIntegralAtCells(const VecDouble &cValues) const
{
  assert(cValues.size() == numCells());
  double acum=0.0;
  for (unsigned i=0;i<cValues.size();i++)
  {
    acum+=cValues(i);
  }
  return acum*getCellMeasure();
}

void FaceInf::getAdjCells(int &posCell,int &negCell) 
{
  posCell = getPosCellIndex();
  negCell = getNegCellIndex();

  
  if (posCell == -1)
    posCell=negCell;
  else if (negCell == -1)
    negCell=posCell;
}

bool DealOrthoMesh::nonZeroInwardNormalComponenentIsPositive(FaceDirection3D face_no) 
{
  //Remember that LEFT_FACE=0,RIGHT_FACE=1,FRONT_FACE=2,BACK_FACE=3,BOTTOM_FACE=4,UP_FACE=5
  static bool a[] = {1,0,1,0,1,0};
  assert(face_no < 6);
  return a[face_no];

  
}

double DealOrthoMesh::faceArea(FaceDirection3D face_no) 
{
  return m_FacesAreaFromCellDirection[face_no];
  
}


void DealOrthoMesh::getFacesInYZSlab(double X,VecIndex &fIndices) 
{
  Cell cell = beginCell();

  double x_tol = getDX()/4.0;
  fIndices.reserve(numElems(1)*numElems(2));
  
  //First find the first face in the slab
  Face face = cell->face(LEFT_FACE);
  FaceDirection3D faceDir;
  if (NumericMethods::a_equal(getBarycenter(face)[0],X,x_tol))
    faceDir=LEFT_FACE;
  else
  {
    faceDir = RIGHT_FACE;
    face=cell->face(RIGHT_FACE);
    while (!NumericMethods::a_equal(getBarycenter(face)[0],X,x_tol))
    {
      cell = cell->neighbor(RIGHT_CELL);
      if (cell.state() != IteratorState::valid)
        return; //Oooops did not find the slab in the domain
      face = cell->face(RIGHT_FACE);
    }
  }

  //Found the slab
  //Now run in Z direction 
  for (;cell.state()==IteratorState::valid;cell=cell->neighbor(UP_CELL))
  {
    //For each Cell in Z, run in Y direction
    for(Cell cellY=cell;cellY.state()==IteratorState::valid;cellY=cellY->neighbor(BACK_CELL))
      fIndices.push_back(cellY->face_index(faceDir));
  }
  
  
  
}


Point3D DealOrthoMesh::getBarycenter(Face &face) 
{
  Point3D V = face->vertex(BL_VERTEX);
  V+=face->vertex(UR_VERTEX);
  V/=2.0;
  return V;
}

---