Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkFEMElementBase.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkFEMElementBase.h,v $
00005   Language:  C++
00006   Date:      $Date: 2009-01-29 21:28:16 $
00007   Version:   $Revision: 1.43 $
00008 
00009   Copyright (c) Insight Software Consortium. All rights reserved.
00010   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00011 
00012      This software is distributed WITHOUT ANY WARRANTY; without even 
00013      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00014      PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 #ifndef __itkFEMElementBase_h
00019 #define __itkFEMElementBase_h
00020 
00021 #include "itkFEMLightObject.h"
00022 #include "itkFEMPArray.h"
00023 #include "itkFEMMaterialBase.h"
00024 #include "itkFEMSolution.h"
00025 #include "itkVisitorDispatcher.h"
00026 #include "vnl/vnl_matrix.h"
00027 #include "vnl/vnl_vector.h"
00028 #include <set>
00029 #include <vector>
00030 
00031 namespace itk {
00032 namespace fem {
00033 
00034 // FIXME: Write better documentation
00067 #define HANDLE_ELEMENT_LOADS() \
00068 
00069  \
00070   typedef void (*LoadImplementationFunctionPointer)(ConstPointer,Element::LoadPointer, Element::VectorType& ); \
00071   virtual void GetLoadVector( Element::LoadPointer l, Element::VectorType& Fe ) const \
00072   { VisitorDispatcher<Self,Element::LoadType, LoadImplementationFunctionPointer>::Visit(l)(this,l,Fe); }
00074 
00075 class Element : public FEMLightObject
00076 {
00077   FEM_ABSTRACT_CLASS(Element,FEMLightObject)
00078 public:
00079 
00083   typedef double Float;
00084 
00088   typedef FEMPArray<Element> ArrayType;
00089 
00093   typedef vnl_matrix<Float> MatrixType;
00094 
00098   typedef vnl_vector<Float> VectorType;
00099 
00111   typedef FEMLightObject    LoadType;
00112   typedef LoadType::Pointer LoadPointer;
00114 
00118   typedef unsigned int DegreeOfFreedomIDType;
00119 
00125   enum{ InvalidDegreeOfFreedomID = 0xffffffff };
00126 
00135   class Node : public FEMLightObject
00136     {
00137     FEM_CLASS(Node,FEMLightObject)
00138      public:
00139 
00143     typedef double Float;
00144 
00148     typedef FEMPArray<Self> ArrayType;
00149 
00150 
00151     /* Windows visualization */
00152   #ifdef FEM_BUILD_VISUALIZATION
00153 
00154     void Draw(CDC* pDC, Solution::ConstPointer sol) const;
00155 
00157     static double& DC_Scale;
00158   #endif
00159 
00163     Node() {}
00164 
00168     Node(Float x, Float y) : m_coordinates(VectorType(2))
00169     { m_coordinates[0]=x; m_coordinates[1]=y; }
00170 
00174     Node(Float x, Float y, Float z) : m_coordinates(VectorType(3))
00175       { m_coordinates[0]=x; m_coordinates[1]=y; m_coordinates[2]=z;}
00176 
00181     const VectorType& GetCoordinates( void ) const
00182       { return m_coordinates; }
00183 
00187     void SetCoordinates( const VectorType& coords )
00188       { m_coordinates=coords; }
00189 
00193     DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int i) const
00194       {
00195       if( i>=m_dof.size() ) { return InvalidDegreeOfFreedomID; }
00196       return m_dof[i];
00197       }
00199 
00203     void SetDegreeOfFreedom(unsigned int i, DegreeOfFreedomIDType dof) const
00204       {
00205       if( i>=m_dof.size() ) { m_dof.resize(i+1, InvalidDegreeOfFreedomID); }
00206       m_dof[i]=dof;
00207       }
00209 
00210     virtual void ClearDegreesOfFreedom( void ) const
00211       {
00212       m_dof.clear();
00213       }
00214 
00215     virtual void Read(  std::istream& f, void* info );
00216     virtual void Write( std::ostream& f ) const;
00217 
00218   public:
00223     typedef std::set<Element*> SetOfElements;
00224     mutable SetOfElements m_elements;
00225 
00226   private:
00230     VectorType m_coordinates;
00231 
00236     mutable std::vector<DegreeOfFreedomIDType> m_dof;
00237 
00238   }; // end class Node
00239 
00241   /*
00242    * Methods related to the physics of the problem.
00243    */
00244 
00245   virtual VectorType GetStrainsAtPoint(const VectorType& pt, const Solution& sol, unsigned int index) const;
00246 
00247   virtual VectorType GetStressesAtPoint(const VectorType& pt, const VectorType& e, const Solution& sol, unsigned int index) const;
00248 
00269   virtual void GetStiffnessMatrix( MatrixType& Ke ) const;
00270 
00280   virtual Float GetElementDeformationEnergy( MatrixType& LocalSolution ) const;
00281 
00293   virtual void GetMassMatrix( MatrixType& Me ) const;
00294 
00306   virtual void GetLandmarkContributionMatrix(float eta, MatrixType& Le) const;
00307 
00337   virtual void GetLoadVector( LoadPointer l, VectorType& Fe ) const = 0;
00338 
00346   virtual void GetStrainDisplacementMatrix( MatrixType& B, const MatrixType& shapeDgl ) const = 0;
00347 
00353   virtual void GetMaterialMatrix( MatrixType& D ) const = 0;
00354 
00366   virtual VectorType InterpolateSolution( const VectorType& pt, const Solution& sol , unsigned int solutionIndex=0 ) const;
00367 
00381   virtual Float InterpolateSolutionN( const VectorType& pt, const Solution& sol, unsigned int f , unsigned int solutionIndex=0 ) const;
00382 
00389   DegreeOfFreedomIDType GetDegreeOfFreedom( unsigned int local_dof ) const
00390     {
00391     if(local_dof>this->GetNumberOfDegreesOfFreedom()) { return InvalidDegreeOfFreedomID; }
00392     return this->GetNode(local_dof/this->GetNumberOfDegreesOfFreedomPerNode())->GetDegreeOfFreedom(local_dof%this->GetNumberOfDegreesOfFreedomPerNode());
00393     }
00395 
00412   virtual Material::ConstPointer GetMaterial(void) const { return 0; }
00413 
00422   virtual void SetMaterial(Material::ConstPointer) {} // FIXME: maybe we should throw an exception instead
00423 
00425 
00451   virtual void GetIntegrationPointAndWeight( unsigned int i, VectorType& pt, Float& w, unsigned int order=0 ) const = 0;
00452 
00463   virtual unsigned int GetNumberOfIntegrationPoints( unsigned int order=0 ) const = 0;
00464 
00473   enum { gaussMaxOrder=10 };
00474 
00486   static const Float gaussPoint[gaussMaxOrder+1][gaussMaxOrder];
00487 
00493   static const Float gaussWeight[gaussMaxOrder+1][gaussMaxOrder];
00494 
00495 
00497   /*
00498    * Methods related to the geometry of an element
00499    */
00500 
00505   typedef Node::ConstPointer NodeIDType;
00506 
00510   virtual unsigned int GetNumberOfNodes( void ) const = 0;
00511 
00515   virtual NodeIDType GetNode(unsigned int n) const = 0;
00516 
00520   virtual void SetNode(unsigned int n, NodeIDType node) = 0;
00521 
00527   virtual const VectorType& GetNodeCoordinates( unsigned int n ) const = 0;
00528 
00534   virtual VectorType GetGlobalFromLocalCoordinates( const VectorType& pt ) const;
00535 
00542   virtual bool GetLocalFromGlobalCoordinates( const VectorType& globalPt , VectorType& localPt ) const = 0;
00543 
00549   virtual unsigned int GetNumberOfSpatialDimensions() const = 0;
00550 
00558   virtual VectorType ShapeFunctions( const VectorType& pt ) const = 0;
00559 
00575   virtual void ShapeFunctionDerivatives( const VectorType& pt, MatrixType& shapeD ) const = 0;
00576 
00596   virtual void ShapeFunctionGlobalDerivatives( const VectorType& pt, MatrixType& shapeDgl, const MatrixType* pJ=0, const MatrixType* pshapeD=0 ) const;
00597 
00619   virtual void Jacobian( const VectorType& pt, MatrixType& J, const MatrixType* pshapeD = 0 ) const;
00620 
00630   virtual Float JacobianDeterminant( const VectorType& pt, const MatrixType* pJ = 0 ) const;
00631 
00643   virtual void JacobianInverse( const VectorType& pt, MatrixType& invJ, const MatrixType* pJ = 0 ) const;
00644 
00650   virtual unsigned int GetNumberOfDegreesOfFreedom( void ) const
00651     {
00652     return this->GetNumberOfNodes() * this->GetNumberOfDegreesOfFreedomPerNode();
00653     }
00654 
00662   virtual unsigned int GetNumberOfDegreesOfFreedomPerNode( void ) const = 0;
00663 
00665 
00669 #ifdef FEM_BUILD_VISUALIZATION
00670 
00673   virtual void Draw(CDC* pDC, Solution::ConstPointer sol) const {}
00674 
00676   static double DC_Scale;
00677 #endif
00678 
00679 };
00680 
00681 // Make sure that Element::Node class is registered with the object factory.
00682 static INITClass Initializer_ElementNode(Element::Node::CLID());
00683 
00684 // Alias for Element::Node class
00685 typedef Element::Node Node;
00686 
00698 class ReadInfoType
00699 {
00700 public:
00701 
00702   typedef Node::ArrayType::ConstPointer     NodeArrayPointer;
00703   typedef Element::ArrayType::ConstPointer  ElementArrayPointer;
00704   typedef Material::ArrayType::ConstPointer MaterialArrayPointer;
00705 
00707   NodeArrayPointer m_node;
00708 
00710   ElementArrayPointer m_el;
00711 
00713   MaterialArrayPointer m_mat;
00714 
00716   ReadInfoType( NodeArrayPointer node_, ElementArrayPointer el_, MaterialArrayPointer mat_) :
00717     m_node(node_), m_el(el_), m_mat(mat_) {}
00718 };
00720 
00721 }} // end namespace itk::fem
00722 
00723 #endif // #ifndef __itkFEMElementBase_h
00724 

Generated at Fri Apr 16 18:18:59 2010 for ITK by doxygen 1.6.1 written by Dimitri van Heesch, © 1997-2000