GeomSimplices.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/GeomSimplex.hpp"
 
// ======================================================
// 0D element
 
template<>
CepsReal GeomPointCell::
getMeasure () const
{
  CEPS_ABORT ("cannot return 0-dim element measure");
  return 0.0;
}
 
template<>
void GeomPointCell::
computeInverseJacobianMatrix ()
{
  checkJacobianDeterminant ();
  //   const JacType  &J    = m_jacobian;
  //   InvJacType     &Jinv = m_invJacobian;
  //   const CepsReal &det  = m_detJacobian;
  m_invJcomputed = true;
  return;
}
 
template<>
void GeomPointCell::
computeJacobianDeterminant ()
{
  m_detJacobian = 1.0;
  return;
}
 
// ======================================================
// 1D element
 
template<>
CepsReal GeomCableCell::
getMeasure () const
{
  return getJacobianDeterminant ();
}
 
template<>
void GeomCableCell::
computeInverseJacobianMatrix ()
{
  checkJacobianDeterminant ();
  const JacobianMatrixType  &J    = m_jacobian;
  InverseJacobianMatrixType &Jinv = m_invJacobian;
  const CepsReal            &det  = m_detJacobian;
 
  Jinv.row (0) = J.col (0) / det / det;
 
  m_invJcomputed = true;
  return;
}
 
template<>
void GeomCableCell::
computeJacobianDeterminant ()
{
  m_detJacobian = getJacobianMatrix ().col (0).norm ();
  return;
}
 
// ======================================================
// 2D element
 
template<>
CepsReal GeomTriangleCell::
getMeasure () const
{
  return getJacobianDeterminant () / 2.0;
}
 
template<>
void GeomTriangleCell::
computeInverseJacobianMatrix ()
{
  checkJacobianDeterminant ();
  const JacobianMatrixType  &J    = m_jacobian;
  InverseJacobianMatrixType &Jinv = m_invJacobian;
 
  CepsReal a       = J (0, 0) * J (0, 0) + J (1, 0) * J (1, 0) + J (2, 0) * J (2, 0);
  CepsReal b       = J (0, 0) * J (0, 1) + J (1, 0) * J (1, 1) + J (2, 0) * J (2, 1);
  CepsReal c       = b;
  CepsReal d       = J (0, 1) * J (0, 1) + J (1, 1) * J (1, 1) + J (2, 1) * J (2, 1);
  CepsReal det_1   = a * d - b * c;
  CepsReal a_m_inv = +d / det_1;
  CepsReal b_m_inv = -b / det_1;
  CepsReal c_m_inv = -c / det_1;
  CepsReal d_m_inv = +a / det_1;
 
  Jinv (0, 0) = a_m_inv * J (0, 0) + b_m_inv * J (0, 1);
  Jinv (1, 0) = c_m_inv * J (0, 0) + d_m_inv * J (0, 1);
  Jinv (0, 1) = a_m_inv * J (1, 0) + b_m_inv * J (1, 1);
  Jinv (1, 1) = c_m_inv * J (1, 0) + d_m_inv * J (1, 1);
  Jinv (0, 2) = a_m_inv * J (2, 0) + b_m_inv * J (2, 1);
  Jinv (1, 2) = c_m_inv * J (2, 0) + d_m_inv * J (2, 1);
 
  m_invJcomputed = true;
  return;
}
 
template<>
void GeomTriangleCell::
computeJacobianDeterminant ()
{
  const JacobianMatrixType &J = getJacobianMatrix ();
  m_detJacobian    = J.col (0).cross (J.col (1)).norm ();
  return;
}
 
// ======================================================
// 3D element
 
template<>
CepsReal GeomTetraCell::
getMeasure () const
{
  return getJacobianDeterminant () / 6.0;
}
 
template<>
void GeomTetraCell::
computeInverseJacobianMatrix ()
{
  checkJacobianDeterminant ();
  const JacobianMatrixType  &J    = m_jacobian;
  InverseJacobianMatrixType &Jinv = m_invJacobian;
  const CepsReal            &det  = m_detJacobian;
 
  Jinv (0, 0) = +(J (1, 1) * J (2, 2) - J (1, 2) * J (2, 1)) / det;
  Jinv (1, 0) = -(J (1, 0) * J (2, 2) - J (1, 2) * J (2, 0)) / det;
  Jinv (2, 0) = +(J (1, 0) * J (2, 1) - J (1, 1) * J (2, 0)) / det;
  Jinv (0, 1) = -(J (0, 1) * J (2, 2) - J (0, 2) * J (2, 1)) / det;
  Jinv (1, 1) = +(J (0, 0) * J (2, 2) - J (0, 2) * J (2, 0)) / det;
  Jinv (2, 1) = -(J (0, 0) * J (2, 1) - J (0, 1) * J (2, 0)) / det;
  Jinv (0, 2) = +(J (0, 1) * J (1, 2) - J (0, 2) * J (1, 1)) / det;
  Jinv (1, 2) = -(J (0, 0) * J (1, 2) - J (0, 2) * J (1, 0)) / det;
  Jinv (2, 2) = +(J (0, 0) * J (1, 1) - J (0, 1) * J (1, 0)) / det;
 
  m_invJcomputed = true;
  return;
}
 
template<>
void GeomTetraCell::
computeJacobianDeterminant ()
{
  const JacobianMatrixType &J = getJacobianMatrix ();
  m_detJacobian               = J.determinant ();
  return;
}