itkFEMElementBase.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMElementBase.h,v $
  Language:  C++
  Date:      $Date: 2009-01-29 21:28:16 $
  Version:   $Revision: 1.43 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#ifndef __itkFEMElementBase_h
#define __itkFEMElementBase_h

#include "itkFEMLightObject.h"
#include "itkFEMPArray.h"
#include "itkFEMMaterialBase.h"
#include "itkFEMSolution.h"
#include "itkVisitorDispatcher.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector.h"
#include <set>
#include <vector>

namespace itk {
namespace fem {

// FIXME: Write better documentation
#define HANDLE_ELEMENT_LOADS() \

 \
  typedef void (*LoadImplementationFunctionPointer)(ConstPointer,Element::LoadPointer, Element::VectorType& ); \
  virtual void GetLoadVector( Element::LoadPointer l, Element::VectorType& Fe ) const \
  { VisitorDispatcher<Self,Element::LoadType, LoadImplementationFunctionPointer>::Visit(l)(this,l,Fe); }

class Element : public FEMLightObject
{
  FEM_ABSTRACT_CLASS(Element,FEMLightObject)
public:

  typedef double Float;

  typedef FEMPArray<Element> ArrayType;

  typedef vnl_matrix<Float> MatrixType;

  typedef vnl_vector<Float> VectorType;

  typedef FEMLightObject    LoadType;
  typedef LoadType::Pointer LoadPointer;

  typedef unsigned int DegreeOfFreedomIDType;

  enum{ InvalidDegreeOfFreedomID = 0xffffffff };

  class Node : public FEMLightObject
    {
    FEM_CLASS(Node,FEMLightObject)
     public:

    typedef double Float;

    typedef FEMPArray<Self> ArrayType;


    /* Windows visualization */
  #ifdef FEM_BUILD_VISUALIZATION

    void Draw(CDC* pDC, Solution::ConstPointer sol) const;

    static double& DC_Scale;
  #endif

    Node() {}

    Node(Float x, Float y) : m_coordinates(VectorType(2))
    { m_coordinates[0]=x; m_coordinates[1]=y; }

    Node(Float x, Float y, Float z) : m_coordinates(VectorType(3))
      { m_coordinates[0]=x; m_coordinates[1]=y; m_coordinates[2]=z;}

    const VectorType& GetCoordinates( void ) const
      { return m_coordinates; }

    void SetCoordinates( const VectorType& coords )
      { m_coordinates=coords; }

    DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int i) const
      {
      if( i>=m_dof.size() ) { return InvalidDegreeOfFreedomID; }
      return m_dof[i];
      }

    void SetDegreeOfFreedom(unsigned int i, DegreeOfFreedomIDType dof) const
      {
      if( i>=m_dof.size() ) { m_dof.resize(i+1, InvalidDegreeOfFreedomID); }
      m_dof[i]=dof;
      }

    virtual void ClearDegreesOfFreedom( void ) const
      {
      m_dof.clear();
      }

    virtual void Read(  std::istream& f, void* info );
    virtual void Write( std::ostream& f ) const;

  public:
    typedef std::set<Element*> SetOfElements;
    mutable SetOfElements m_elements;

  private:
    VectorType m_coordinates;

    mutable std::vector<DegreeOfFreedomIDType> m_dof;

  }; // end class Node

  /*
   * Methods related to the physics of the problem.
   */

  virtual VectorType GetStrainsAtPoint(const VectorType& pt, const Solution& sol, unsigned int index) const;

  virtual VectorType GetStressesAtPoint(const VectorType& pt, const VectorType& e, const Solution& sol, unsigned int index) const;

  virtual void GetStiffnessMatrix( MatrixType& Ke ) const;

  virtual Float GetElementDeformationEnergy( MatrixType& LocalSolution ) const;

  virtual void GetMassMatrix( MatrixType& Me ) const;

  virtual void GetLandmarkContributionMatrix(float eta, MatrixType& Le) const;

  virtual void GetLoadVector( LoadPointer l, VectorType& Fe ) const = 0;

  virtual void GetStrainDisplacementMatrix( MatrixType& B, const MatrixType& shapeDgl ) const = 0;

  virtual void GetMaterialMatrix( MatrixType& D ) const = 0;

  virtual VectorType InterpolateSolution( const VectorType& pt, const Solution& sol , unsigned int solutionIndex=0 ) const;

  virtual Float InterpolateSolutionN( const VectorType& pt, const Solution& sol, unsigned int f , unsigned int solutionIndex=0 ) const;

  DegreeOfFreedomIDType GetDegreeOfFreedom( unsigned int local_dof ) const
    {
    if(local_dof>this->GetNumberOfDegreesOfFreedom()) { return InvalidDegreeOfFreedomID; }
    return this->GetNode(local_dof/this->GetNumberOfDegreesOfFreedomPerNode())->GetDegreeOfFreedom(local_dof%this->GetNumberOfDegreesOfFreedomPerNode());
    }

  virtual Material::ConstPointer GetMaterial(void) const { return 0; }

  virtual void SetMaterial(Material::ConstPointer) {} // FIXME: maybe we should throw an exception instead


  virtual void GetIntegrationPointAndWeight( unsigned int i, VectorType& pt, Float& w, unsigned int order=0 ) const = 0;

  virtual unsigned int GetNumberOfIntegrationPoints( unsigned int order=0 ) const = 0;

  enum { gaussMaxOrder=10 };

  static const Float gaussPoint[gaussMaxOrder+1][gaussMaxOrder];

  static const Float gaussWeight[gaussMaxOrder+1][gaussMaxOrder];


  /*
   * Methods related to the geometry of an element
   */

  typedef Node::ConstPointer NodeIDType;

  virtual unsigned int GetNumberOfNodes( void ) const = 0;

  virtual NodeIDType GetNode(unsigned int n) const = 0;

  virtual void SetNode(unsigned int n, NodeIDType node) = 0;

  virtual const VectorType& GetNodeCoordinates( unsigned int n ) const = 0;

  virtual VectorType GetGlobalFromLocalCoordinates( const VectorType& pt ) const;

  virtual bool GetLocalFromGlobalCoordinates( const VectorType& globalPt , VectorType& localPt ) const = 0;

  virtual unsigned int GetNumberOfSpatialDimensions() const = 0;

  virtual VectorType ShapeFunctions( const VectorType& pt ) const = 0;

  virtual void ShapeFunctionDerivatives( const VectorType& pt, MatrixType& shapeD ) const = 0;

  virtual void ShapeFunctionGlobalDerivatives( const VectorType& pt, MatrixType& shapeDgl, const MatrixType* pJ=0, const MatrixType* pshapeD=0 ) const;

  virtual void Jacobian( const VectorType& pt, MatrixType& J, const MatrixType* pshapeD = 0 ) const;

  virtual Float JacobianDeterminant( const VectorType& pt, const MatrixType* pJ = 0 ) const;

  virtual void JacobianInverse( const VectorType& pt, MatrixType& invJ, const MatrixType* pJ = 0 ) const;

  virtual unsigned int GetNumberOfDegreesOfFreedom( void ) const
    {
    return this->GetNumberOfNodes() * this->GetNumberOfDegreesOfFreedomPerNode();
    }

  virtual unsigned int GetNumberOfDegreesOfFreedomPerNode( void ) const = 0;


#ifdef FEM_BUILD_VISUALIZATION

  virtual void Draw(CDC* pDC, Solution::ConstPointer sol) const {}

  static double DC_Scale;
#endif

};

// Make sure that Element::Node class is registered with the object factory.
static INITClass Initializer_ElementNode(Element::Node::CLID());

// Alias for Element::Node class
typedef Element::Node Node;

class ReadInfoType
{
public:

  typedef Node::ArrayType::ConstPointer     NodeArrayPointer;
  typedef Element::ArrayType::ConstPointer  ElementArrayPointer;
  typedef Material::ArrayType::ConstPointer MaterialArrayPointer;

  NodeArrayPointer m_node;

  ElementArrayPointer m_el;

  MaterialArrayPointer m_mat;

  ReadInfoType( NodeArrayPointer node_, ElementArrayPointer el_, MaterialArrayPointer mat_) :
    m_node(node_), m_el(el_), m_mat(mat_) {}
};

}} // end namespace itk::fem

#endif // #ifndef __itkFEMElementBase_h