itkfemelementbase.h

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMElementBase.h,v $
  Language:  C++
  Date:      $Date: 2003-09-10 14:29:41 $
  Version:   $Revision: 1.42 $

  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
/**
 * \class Element
 * \brief Abstract base element class.
 *
 * Derive this class to create new finite element classes.
 * All derived classes must define:
 *
 *    - Ke():    Function to calculate the element stiffness matrix in global coordinate system.
 *    - Fe():    Function to calculate the element force vector in global coordinate system.
 *    - uDOF():  Provide a pointer to storage of i-th DOF displacement in the element.
 *    - Clone(): Function that creates a duplicate of current element and returns a pointer to it.
 *
 * and optionally (if required):
 *    - Read():  Reads element data from the stream f. assume that the stream position is
 *               already where the element data starts. Take care of the error checking.
 *    - Write(): Writes element data to the stream.
 *    - Draw():  Draws the element on the device context (Windows only).
 *
 * The storage of element parameters (geometry...) can't be implemented here, since we don't know yet, 
 * how much memory each element needs. Instead each derived class should take care of the memory 
 * management (declare appropriate data members) for the element parameters and provide access
 * to these parameters (like nodes, materials...).
 */


/**
 * \def HANDLE_ELEMENT_LOADS()
 * \brief Macro that simplifies the the GetLoadVector function definitions.
 *
 * NOTE: This macro must be called in declaration of ALL
 *       derived Element classes.
 */
#define HANDLE_ELEMENT_LOADS() \
  /** Pointer type that specifies functions that can handle loads on this element */ \
  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:

  /**
   * Floating point type used in all Element classes.
   */
  typedef double Float;

  /**
   * Array class that holds special pointers to the Element objects
   */
  typedef FEMPArray<Element> ArrayType;

  /**
   * Class used to store the element stiffness matrix
   */
  typedef vnl_matrix<Float> MatrixType;

  /**
   * Class to store the element load vector
   */
  typedef vnl_vector<Float> VectorType;

  /**
   * Easy and consistent access to LoadElement and LoadElement::Pointer type.
   * This is a pointer to FEMLightObject to avoid cyclic references between
   * LoadElement and Element classes.
   * As a consequence whenever you need to use a pointer to LoadElement class
   * within the element's declaration or definition, ALWAYS use this typedef
   * instead.
   * When calling the GetLoadVector(...) function from outside, you should
   * ALWAYS first convert the argument to Element::LoadPointer. See
   * code of function Solver::AssembleF(...) for more info.
   */
  typedef FEMLightObject LoadType;
  typedef LoadType::Pointer LoadPointer;

  /**
   * Type that stores global ID's of degrees of freedom.
   */
  typedef unsigned int DegreeOfFreedomIDType;

  /**
   * Constant that represents an invalid DegreeOfFreedomID object.
   * If a degree of freedom is assigned this value, this means that
   * that no specific value was (yet) assigned to this DOF.
   */ 
  enum{ InvalidDegreeOfFreedomID = 0xffffffff };

  /**
   * \class Node
   * \brief Class that stores information required to define a node.
   *
   * A node can define a point in space and can hold an arbitrary number 
   * of coordinates and the DOFs. Since the only classes that use nodes
   * are the elements, the node class is defined within an element base class.
   */
  class Node : public FEMLightObject
  {
  FEM_CLASS(Node,FEMLightObject)
  public:

    /**
     * Floating point precision type.
     */
    typedef double Float;

    /**
     * Array class that holds special pointers to the nodes.
     */
    typedef FEMPArray<Self> ArrayType;


    /* Windows visualization */
  #ifdef FEM_BUILD_VISUALIZATION
    /** Draws the node on the DC */
    void Draw(CDC* pDC, Solution::ConstPointer sol) const;
    /** Global scale for drawing on the DC */
    static double& DC_Scale;
  #endif

    /**
     * Default constructor
     */
    Node() {}

    /**
     * Create 2D node.
     */
    Node(Float x, Float y) : m_coordinates(VectorType(2))
    { m_coordinates[0]=x; m_coordinates[1]=y; }

    /**
     * Create 3D node.
     */
        Node(Float x, Float y, Float z) : m_coordinates(VectorType(3))
    { m_coordinates[0]=x; m_coordinates[1]=y; m_coordinates[2]=z;}
    /**
     * Return a reference to a vector that contains coordinates
     * of this node.
     */
    const VectorType& GetCoordinates( void ) const
    { return m_coordinates; }

    /**
     * Set coordinates of a node.
     */
    void SetCoordinates( const VectorType& coords )
    { m_coordinates=coords; }

    /**
     * Get DOF IDs associated with this node.
     */
    DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int i) const
    {
      if( i>=m_dof.size() ) { return InvalidDegreeOfFreedomID; }
      return m_dof[i];
    }

    /**
     * Set DOF IDs associated with this node.
     */
    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:
    /**
     * List of pointers to elements that use this node. External code is
     * responsible for maintaining the list.
     */
    typedef std::set<Element*> SetOfElements;
    mutable SetOfElements m_elements;

  private:
    /**
     * Vector object that holds node coordinates.
     */
    VectorType m_coordinates;

    /**
     * Array that holds IDs of degrees of freedom that are
     * defined at this node.
     */
    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;

  /**
   * Compute and return element stiffnes matrix (Ke) in global coordinate
   * system.
   * The base class provides a general implementation which only computes
   *
   *     b   T
   * int    B(x) D B(x) dx
   *     a
   *
   * using the Gaussian numeric integration method. The function calls
   * GetIntegrationPointAndWeight() / GetNumberOfIntegrationPoints() to obtain
   * the integration points. It also calls the GetStrainDisplacementMatrix()
   * and GetMaterialMatrix() member functions.
   *
   * \param Ke Reference to the resulting stiffnes matrix.
   *
   * \note This is a very generic implementation of the stiffness matrix
   *       that is suitable for any problem/element definition. A specifc
   *       element may override this implementation with its own simple one.
   */
  virtual void GetStiffnessMatrix( MatrixType& Ke ) const;

  /**
   * Compute the physical energy, U, of the deformation (e.g. stress / strain ).
   *
   *      T
   * U = u  Ke u
   *         
   * The matrix LocalSolution contains the solution to use in the energy
   * computation.  Usually, this is the solution at the nodes.
   */
  virtual Float GetElementDeformationEnergy( MatrixType& LocalSolution ) const;

  /**
   * Compute and return element mass matrix (Me) in global coordinate system.
   *
   *     b   T
   * int    N(x) (rho c) N(x) dx
   *     a
   *
   * where (rho c) is constant (element density), which is here assumed to be
   * equal to one. If this is not the case, this function must be overriden in
   * a derived class. Implementation is similar to GetStiffnessMatrix.
   */
  virtual void GetMassMatrix( MatrixType& Me ) const;

  /**
   * Compute and return landmark contribution to element stiffness
   * matrix (Le) in global coordinate system.
   *
   *     b             T
   * int   (1/eta)^2  N(x) N(x) dx
   *     a
   *
   * where (eta ) is the landmark weight.  Implementation is similar
   * to GetMassMatrix.
   */
  virtual void GetLandmarkContributionMatrix(float eta, MatrixType& Le) const;

  /**
   * Compute and return the element load vector for a given external load.
   * The class of load object determines the type of load acting on the
   * elemnent. Basically this is the contribution of this element on the right
   * side of the master matrix equation, due to the specified load. 
   * Returned vector includes only nodal forces that correspond to the given
   * Load object.
   *
   * Visitor design pattern is used in the loads implementation. This function
   * only selects and calls the proper function based on the given class of
   * load object. The code that performs the actual conversion to the
   * corresponding nodal loads is defined elswhere.
   *
   * \note Each derived class must implement its own version of this function.
   *       This is automated by calling the LOAD_FUNCTION() macro within the
   *       class declaration (in the public: block).
   *
   * For example on how to define specific element load, see funtion
   * LoadImplementationPoint_Bar2D.
   *
   * \note: Before a load can be applied to an element, the function that
   *        implements a load must be registered with the VisitorDispactcher
   *        class.
   *
   * \param l Pointer to a load object.
   * \param Fe Reference to vector object that will store nodal forces.
   *
   * \sa VisitorDispatcher
   */
  virtual void GetLoadVector( LoadPointer l, VectorType& Fe ) const = 0;

  /**
   * Compute the strain displacement matrix at local point.
   *
   * \param B Reference to a matrix object that will contain the result
   * \param shapeDgl Matrix that contains derivatives of shape functions
   *                 w.r.t. global coordinates.
   */
  virtual void GetStrainDisplacementMatrix( MatrixType& B, const MatrixType& shapeDgl ) const = 0;

  /**
   * Compute the element material matrix.
   *
   * \param D Reference to a matrix object
   */
  virtual void GetMaterialMatrix( MatrixType& D ) const = 0;

  /**
   * Return interpolated value of all unknown functions at
   * given local point.
   *
   * \param pt Point in local element coordinates.
   * \param sol Reference to the master solution object. This object
   *            is created by the Solver object when the whole FEM problem
   *            is solved and contains the values of unknown functions
   *            at nodes (degrees of freedom).
   * \param solutionIndex We allow more than one solution vector to be stored - this selects which to use in interpolation.
   */
  virtual VectorType InterpolateSolution( const VectorType& pt, const Solution& sol , unsigned int solutionIndex=0 ) const;

  /**