SimplexJac.h

Go to the documentation of this file.

// -*- C++ -*-

#if !defined(__geom_SimplexJac_h__)
#define __geom_SimplexJac_h__

#if defined(DEBUG_geom) && !defined(DEBUG_SimplexJac)
#define DEBUG_SimplexJac
#endif

#include "geometry.h"

#include "../../../ads/tensor/SquareMatrix.h"

BEGIN_NAMESPACE_GEOM


template<int N, typename T = double>
class SimplexJac {
public:

  //
  // public typedefs
  //

  typedef T Number;

  typedef ads::FixedArray<N,Number> Vertex;

  typedef geom::Simplex<N,Vertex,Number> Simplex;

  typedef ads::SquareMatrix<N,Number> Matrix;

private:

  //
  // Member data.
  //

  // The Jacobian matrix that maps the identity simplex to this 
  // simplex (after the first vertex has been translated to the origin.)
  Matrix _matrix;

  // The determinant of _matrix.
  Number _determinant;

  // The gradient of the determinant of the Jacobian matrix.
  Vertex _gradientDeterminant;

  //
  // Static member data.
  //

  // The Jacobian matrix that maps the identity simplex to the reference 
  // simplex.  In 2-D, it maps the identity triangle:
  // (0,0), (1,0), (1/2,sqrt(3)/2)
  // to the reference triangle:
  // (0,0), (1,0), (0,1).
  // In 3-D it maps the identity tetrahedron:
  // (0,0,0), (1,0,0), (1/2,sqrt(3)/2,0), (1/2,sqrt(3)/6,sqrt(2)/sqrt(3))
  // to the reference tetrahedron:
  // (0,0,0), (1,0,0), (0,1,0), (0,0,1).
  static Matrix _identityToReference;

  // The gradient of the Jacobian matrix.
  static ads::FixedArray<N,Matrix> _gradientMatrix;

  // The determinant of _identityToReference;
  static Number _determinantIdentityToReference;

  // N!
  static Number _nFactorial;

public:

  //--------------------------------------------------------------------------

  SimplexJac() :
    _matrix(),
    _determinant(),
    _gradientDeterminant()
  {}

  SimplexJac(const SimplexJac& other) :
    _matrix(other._matrix),
    _determinant(other._determinant),
    _gradientDeterminant(other._gradientDeterminant)
  {}


  SimplexJac(const Simplex& s) {
    set(s);
  }

  SimplexJac& 
  operator=(const SimplexJac& other) {
    if (&other != this) {
      _matrix = other._matrix;
      _determinant = other._determinant;
      _gradientDeterminant = other._gradientDeterminant;
    }
    return *this;
  }
  
  ~SimplexJac()
  {}

  //--------------------------------------------------------------------------

  const Matrix&
  getMatrix() const {
    return _matrix;
  }

  const ads::FixedArray<N,Matrix>&
  getGradientMatrix() const {
    return _gradientMatrix;
  }

  Number
  getDeterminant() const {
    return _determinant;
  }

  const Vertex&
  getGradientDeterminant() const {
    return _gradientDeterminant;
  }

  Number
  computeContent() const {
    return _determinant / _nFactorial / _determinantIdentityToReference;
  }

  void
  computeGradientContent(Vertex* grad) const {
    *grad = _gradientDeterminant;
    *grad /= (_nFactorial * _determinantIdentityToReference);
  }

  Vertex
  computeGradientContent() const {
    Vertex grad;
    computeGradientContent(&grad);
    return grad;
  }

  //--------------------------------------------------------------------------

  void
  setFunction(const Simplex& s);

  void
  set(const Simplex& s);

  // CONTINUE: These kind of functions do not know the orientation of the 
  // simplex.  A 3-2 mesh could become tangled and still have high quality.

  void
  setFunction(const geom::Simplex<N,ads::FixedArray<N+1,Number>,Number>& s);


  void
  set(const geom::Simplex<N,ads::FixedArray<N+1,Number>,Number>& s);


};

END_NAMESPACE_GEOM

#define __geom_SimplexJac_ipp__
#include "SimplexJac.ipp"
#undef __geom_SimplexJac_ipp__

#endif

---