FiniteElement.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _FINITE_ELEMENT_H
#define _FINITE_ELEMENT_H
 
#include "BasisFunctionVals.h"
#include "ItgPoint.h"
#include "MatVec.h"
#include "Vec3.h"
#include "Tensor.h"
 
 
class FiniteElement : public ItgPoint
{
public:
  explicit FiniteElement(size_t n = 0, size_t i = 0) : ItgPoint(i), N(n), Te(3)
  { p = q = r = 0; age = h = 0.0; detJxW = 1.0; }
 
  virtual size_t getNoBasis() const { return 1; }
 
  virtual const Vector& basis(char) const { return N; }
  virtual Vector& basis(char) { return N; }
 
  virtual const Matrix& grad(char) const { return dNdX; }
  virtual const Matrix3D& hess(char) const { return d2NdX2; }
  virtual const Matrix4D& hess2(char) const { return d3NdX3; }
 
protected:
  virtual Matrix& grad(char) { return dNdX; }
  virtual Matrix3D& hess(char) { return d2NdX2; }
 
  virtual std::ostream& write(std::ostream& os) const;
 
  friend std::ostream& operator<<(std::ostream& os, const FiniteElement& fe)
  {
    return fe.write(os);
  }
 
public:
  // Gauss point quantities
  double   detJxW; 
  Vector     N;    
  Matrix    dNdX;  
  Matrix3D d2NdX2; 
  Matrix4D d3NdX3; 
  Matrix     G;    
  Matrix     H;    
 
  Matrix P;
 
  Matrix dPdX;
 
  // Element quantities
  short int           p;    
  short int           q;    
  short int           r;    
  double              age;  
  double              h;    
  Vec3Vec             XC;   
  Vector              Navg; 
  Matrix              Xn;   
  Tensor              Te;   
  std::vector<Tensor> Tn;   
  Vec3Vec             En;   
};
 
 
class MxFiniteElement : public FiniteElement
{
public:
  explicit MxFiniteElement(const std::vector<size_t>& n, size_t ip = 0);
 
  size_t getNoBasis() const override { return 1+M.size(); }
 
  const Vector& basis(char b) const override { return b < 2 ? N : M[b-2]; }
  Vector& basis(char b) override { return b < 2 ? N : M[b-2]; }
 
  const Matrix& grad(char b) const override { return b < 2 ? dNdX : dMdX[b-2]; }
  const Matrix3D& hess(char b) const override { return b < 2 ? d2NdX2 : d2MdX2[b-2]; }
  const Matrix4D& hess2(char b) const override { return b < 2 ? d3NdX3 : d3MdX3[b-2]; }
 
  bool Jacobian(Matrix& Jac, const Matrix& Xnod,
                unsigned short int gBasis,
                const BasisValuesPtrs& bfs);
 
  bool Jacobian(Matrix& Jac, Vec3& n, const Matrix& Xnod,
                unsigned short int gBasis,
                const BasisValuesPtrs& bfs,
                size_t t1, size_t t2, size_t nBasis = 0,
                const Matrix* Xnod2 = nullptr);
 
  bool Hessian(Matrix3D& Hess, const Matrix& Jac, const Matrix& Xnod,
               unsigned short int gBasis,
               const BasisValuesPtrs& bfs);
 
  void piolaMapping (const double detJ,
                     const Matrix& Ji,
                     const Matrix& Xnod,
                     const BasisValuesPtrs& bfs);
 
  void piolaBasis(const double detJ, const Matrix& J);
 
protected:
  void piolaGradient(const double detJ, const Matrix& J, const Matrix& Ji,
                     const Matrix& Xnod, const BasisValuesPtrs& bfs);
 
  Matrix& grad(char b) override { return b < 2 ? dNdX : dMdX[b-2]; }
  Matrix3D& hess(char b) override { return b < 2 ? d2NdX2 : d2MdX2[b-2]; }
 
  std::ostream& write(std::ostream& os) const override;
 
private:
  std::vector<Vector>     M;    
  std::vector<Matrix>    dMdX;  
  std::vector<Matrix3D> d2MdX2; 
  std::vector<Matrix4D> d3MdX3; 
};
 
#endif