GeomCell.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   CEPS is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "geometry/base/GeomCell.hpp"
 
GeomCell::GeomCell() :
    m_cellType     (CepsCellType::Simplex),
    m_jacobian     (JacobianMatrixType (3, 0)),
    m_invJacobian  (InverseJacobianMatrixType (0, 3)),
    m_detJacobian  (0.),
    m_invJcomputed (false),
    m_normal       (CepsMathVertex(0.,0.,0.)),
    m_diameter     (-1.),
    m_nodeSharingCells({})
{}
 
GeomCell::GeomCell(CepsVector<GeomNode*> nodes,const CepsCellGlobalIndex& gID):
  GeomCell()
{
  m_globalIndex = gID;
  m_nodes       = nodes;
}
 
GeomNode*
GeomCell::getNodeWithMinimalIndex() const
{
  CEPS_ABORT_IF(m_nodes.size()==0,"trying to get node with minimal index of cell with no nodes");
 
  const GeomNode*      node = getNodeAt(0);
  CepsNodeGlobalIndex  id   = std::numeric_limits<CepsNodeGlobalIndex>::max();
 
  for (auto n : m_nodes)
    if (ceps::isValidPtr(n) and n->getGlobalIndex()<id)
    {
      node = n;
      id   = node->getGlobalIndex();
    }
 
  CEPS_ABORT_IF(ceps::isNullPtr(node),
    "nullptr node, so we can not deduce the node with minimal index"
  );
  return const_cast<GeomNode*>(node);
}
 
const CepsCellType&
GeomCell::getCellType() const
{
  return m_cellType;
}
 
void
GeomCell::reset()
{
  ceps::HoldsGlobalIndex::reset();
  ceps::HoldsDimension::reset();
  ceps::HoldsBoundary::reset();
  ceps::HoldsAttributes::reset();
  ceps::HoldsProcIds::reset();
  ceps::HoldsNodes<GeomNode>::reset();
  ceps::HoldsCells<GeomCell>::reset();
  return;
}
 
// --------------------------------------------------------------------------------
// Jacobian
// --------------------------------------------------------------------------------
 
CepsReal
GeomCell::getJacobianDeterminant() const
{
  return m_detJacobian;
}
 
const GeomCell::JacobianMatrixType&
GeomCell::getJacobianMatrix()
{
  CEPS_ABORT_IF (this->getDimension()==0,
    "Trying to get jacobian matrix of element with global index #"
    << getGlobalIndex () << "\nwhich is a 0D element"
  );
  return m_jacobian;
}
 
const GeomCell::InverseJacobianMatrixType&
GeomCell::getInverseJacobianMatrix()
{
  
  CEPS_ABORT_IF (this->getDimension()==0,
    "Trying to get jacobian matrix of element with global index #"
    << getGlobalIndex() << std::endl << "which is a 0D element"
  );
  if (not this->m_invJcomputed or this->m_invJacobian.rows() == 0)
    computeInverseJacobianMatrix();
  return m_invJacobian;
}
 
void
GeomCell::refresh()
{
  initialize (getDimension (), getCellType (), getNumberOfNodes ());
 
  computeJacobianMatrix ();
  computeJacobianDeterminant ();
  
  CEPS_ABORT_IF(not checkJacobianDeterminant (),
    "Determinant of geometrical element #" << getGlobalIndex () << "is non-positive."
  );
  return;
}
 
void 
GeomCell::setNormal(CepsMathVertex n)
{
  m_normal = n;
}
 
CepsMathVertex
GeomCell::getNormal() const
{
  return m_normal;
}
 
void
GeomCell::setNodeSharingCellsIndices(const CepsSet<CepsCellGlobalIndex>& nCellIds)
{
  m_nodeSharingCells = nCellIds;
}
 
const CepsSet<CepsCellGlobalIndex>&
GeomCell::getNodeSharingCellsIndices() const
{
  return m_nodeSharingCells;
}
 
CepsReal
GeomCell::getDiameter()
{
  if (m_diameter<0)
    computeDiameter();
  return m_diameter;
}
 
void
GeomCell::initialize(CepsUInt dim, CepsCellType cellType, CepsUInt nNodes)
{
  m_cellType = cellType;
  setDimension (dim);
  m_jacobian.setZero (3, nNodes-1);
  m_invJacobian.setZero (nNodes-1, 3);
  m_invJcomputed = false;
  return;
}
 
void
GeomCell::initializeNodes(GeomNode *const *nodes, CepsUInt nNodes)
{
  m_nodes.resize(nNodes);
  std::copy(nodes, nodes + nNodes, m_nodes.begin ());
  refresh();
  return;
}
 
void
GeomCell::computeDiameter()
{
  CepsUInt size = m_nodes.size();
  for (CepsUInt i=0  ;i<size;i++)
  for (CepsUInt j=i+1;j<size;j++)
  {
    CepsReal3D n1 = getNodeAt(i)->getCoordinates();
    CepsReal3D n2 = getNodeAt(j)->getCoordinates();
    m_diameter = std::max(m_diameter,CepsVertex(n1-n2).norm2());
  }
}
 
CepsStandardArgs
getStandardArgsFrom(GeomCell* cell)
{
  CEPS_ABORT_IF(ceps::isNullPtr(cell),
    "Passed cell is nullptr"
  );
  CepsStandardArgs args;
  args.cellId = cell->getGlobalIndex();
  args.attr   = cell->getAttributes();
  return args;
}