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: 2003-09-10 14:29:41 $
00007   Version:   $Revision: 1.42 $
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 
00035 
00036 // FIXME: Write better documentation
00069 #define HANDLE_ELEMENT_LOADS() \
00070 
00071  \
00072   typedef void (*LoadImplementationFunctionPointer)(ConstPointer,Element::LoadPointer, Element::VectorType& ); \
00073   virtual void GetLoadVector( Element::LoadPointer l, Element::VectorType& Fe ) const \
00074   { VisitorDispatcher<Self,Element::LoadType, LoadImplementationFunctionPointer>::Visit(l)(this,l,Fe); }
00076 
00077 class Element : public FEMLightObject
00078 {
00079 FEM_ABSTRACT_CLASS(Element,FEMLightObject)
00080 public:
00081 
00085   typedef double Float;
00086 
00090   typedef FEMPArray<Element> ArrayType;
00091 
00095   typedef vnl_matrix<Float> MatrixType;
00096 
00100   typedef vnl_vector<Float> VectorType;
00101 
00113   typedef FEMLightObject LoadType;
00114   typedef LoadType::Pointer LoadPointer;
00116 
00120   typedef unsigned int DegreeOfFreedomIDType;
00121 
00127   enum{ InvalidDegreeOfFreedomID = 0xffffffff };
00128 
00137   class Node : public FEMLightObject
00138   {
00139   FEM_CLASS(Node,FEMLightObject)
00140   public:
00141 
00145     typedef double Float;
00146 
00150     typedef FEMPArray<Self> ArrayType;
00151 
00152 
00153     /* Windows visualization */
00154   #ifdef FEM_BUILD_VISUALIZATION
00155 
00156     void Draw(CDC* pDC, Solution::ConstPointer sol) const;
00157 
00159     static double& DC_Scale;
00160   #endif
00161 
00165     Node() {}
00166 
00170     Node(Float x, Float y) : m_coordinates(VectorType(2))
00171     { m_coordinates[0]=x; m_coordinates[1]=y; }
00172 
00176         Node(Float x, Float y, Float z) : m_coordinates(VectorType(3))
00177     { m_coordinates[0]=x; m_coordinates[1]=y; m_coordinates[2]=z;}
00178 
00183     const VectorType& GetCoordinates( void ) const
00184     { return m_coordinates; }
00185 
00189     void SetCoordinates( const VectorType& coords )
00190     { m_coordinates=coords; }
00191 
00195     DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int i) const
00196     {
00197       if( i>=m_dof.size() ) { return InvalidDegreeOfFreedomID; }
00198       return m_dof[i];
00199     }
00201 
00205     void SetDegreeOfFreedom(unsigned int i, DegreeOfFreedomIDType dof) const
00206     {
00207       if( i>=m_dof.size() ) { m_dof.resize(i+1, InvalidDegreeOfFreedomID); }
00208       m_dof[i]=dof;
00209     }
00211 
00212     virtual void ClearDegreesOfFreedom( void ) const
00213     {
00214       m_dof.clear();
00215     }
00216 
00217     virtual void Read(  std::istream& f, void* info );
00218     virtual void Write( std::ostream& f ) const;
00219 
00220   public:
00225     typedef std::set<Element*> SetOfElements;
00226     mutable SetOfElements m_elements;
00227 
00228   private:
00232     VectorType m_coordinates;
00233 
00238     mutable std::vector<DegreeOfFreedomIDType> m_dof;
00239 
00240   }; // end class Node
00241 
00242 
00243 
00244 
00246   /*
00247    * Methods related to the physics of the problem.
00248    */
00249 
00250   virtual VectorType GetStrainsAtPoint(const VectorType& pt, const Solution& sol, unsigned int index) const;
00251 
00252   virtual VectorType GetStressesAtPoint(const VectorType& pt, const VectorType& e, const Solution& sol, unsigned int index) const;
00253 
00274   virtual void GetStiffnessMatrix( MatrixType& Ke ) const;
00275 
00285   virtual Float GetElementDeformationEnergy( MatrixType& LocalSolution ) const;
00286 
00298   virtual void GetMassMatrix( MatrixType& Me ) const;
00299 
00311   virtual void GetLandmarkContributionMatrix(float eta, MatrixType& Le) const;
00312 
00342   virtual void GetLoadVector( LoadPointer l, VectorType& Fe ) const = 0;
00343 
00351   virtual void GetStrainDisplacementMatrix( MatrixType& B, const MatrixType& shapeDgl ) const = 0;
00352 
00358   virtual void GetMaterialMatrix( MatrixType& D ) const = 0;
00359 
00371   virtual VectorType InterpolateSolution( const VectorType& pt, const Solution& sol , unsigned int solutionIndex=0 ) const;
00372 
00386   virtual Float InterpolateSolutionN( const VectorType& pt, const Solution& sol, unsigned int f , unsigned int solutionIndex=0 ) const;
00387 
00394   DegreeOfFreedomIDType GetDegreeOfFreedom( unsigned int local_dof ) const
00395   {
00396     if(local_dof>this->GetNumberOfDegreesOfFreedom()) { return InvalidDegreeOfFreedomID; }
00397     return this->GetNode(local_dof/this->GetNumberOfDegreesOfFreedomPerNode())->GetDegreeOfFreedom(local_dof%this->GetNumberOfDegreesOfFreedomPerNode());
00398   }
00400 
00417   virtual Material::ConstPointer GetMaterial(void) const { return 0; }
00418 
00427   virtual void SetMaterial(Material::ConstPointer) {} // FIXME: maybe we should throw an exception instead
00428 
00429 
00430 
00431 
00433   /*
00434    * Methods related to numeric integration
00435    */
00436 
00459   virtual void GetIntegrationPointAndWeight( unsigned int i, VectorType& pt, Float& w, unsigned int order=0 ) const = 0;
00460 
00471   virtual unsigned int GetNumberOfIntegrationPoints( unsigned int order=0 ) const = 0;
00472 
00481   enum { gaussMaxOrder=10 };
00482 
00494   static const Float gaussPoint[gaussMaxOrder+1][gaussMaxOrder];
00495 
00501   static const Float gaussWeight[gaussMaxOrder+1][gaussMaxOrder];
00502 
00503 
00505   /*
00506    * Methods related to the geometry of an element
00507    */
00508 
00513   typedef Node::ConstPointer NodeIDType;
00514 
00518   virtual unsigned int GetNumberOfNodes( void ) const = 0;
00519 
00523   virtual NodeIDType GetNode(unsigned int n) const = 0;
00524 
00528   virtual void SetNode(unsigned int n, NodeIDType node) = 0;
00529 
00535   virtual const VectorType& GetNodeCoordinates( unsigned int n ) const = 0;
00536 
00542   virtual VectorType GetGlobalFromLocalCoordinates( const VectorType& pt ) const;
00543 
00550   virtual bool GetLocalFromGlobalCoordinates( const VectorType& globalPt , VectorType& localPt ) const = 0;
00551 
00557   virtual unsigned int GetNumberOfSpatialDimensions() const = 0;
00558 
00566   virtual VectorType ShapeFunctions( const VectorType& pt ) const = 0;
00567 
00583   virtual void ShapeFunctionDerivatives( const VectorType& pt, MatrixType& shapeD ) const = 0;
00584 
00604   virtual void ShapeFunctionGlobalDerivatives( const VectorType& pt, MatrixType& shapeDgl, const MatrixType* pJ=0, const MatrixType* pshapeD=0 ) const;
00605 
00627   virtual void Jacobian( const VectorType& pt, MatrixType& J, const MatrixType* pshapeD = 0 ) const;
00628 
00638   virtual Float JacobianDeterminant( const VectorType& pt, const MatrixType* pJ = 0 ) const;
00639 
00651   virtual void JacobianInverse( const VectorType& pt, MatrixType& invJ, const MatrixType* pJ = 0 ) const;
00652 
00658   virtual unsigned int GetNumberOfDegreesOfFreedom( void ) const
00659   {
00660     return this->GetNumberOfNodes() * this->GetNumberOfDegreesOfFreedomPerNode();
00661   }
00662 
00670   virtual unsigned int GetNumberOfDegreesOfFreedomPerNode( void ) const = 0;
00671 
00672 
00673 
00674 
00676   /*
00677    * Methods and classes related to IO and drawing
00678    */
00679 
00680 #ifdef FEM_BUILD_VISUALIZATION
00681 
00684   virtual void Draw(CDC* pDC, Solution::ConstPointer sol) const {}
00685 
00687   static double DC_Scale;
00688 #endif
00689 
00690 };
00691 
00692 // Make sure that Element::Node class is registered with the object factory.
00693 static INITClass Initializer_ElementNode(Element::Node::CLID());
00694 
00695 // Alias for Element::Node class
00696 typedef Element::Node Node;
00697 
00698 
00699 
00700 
00712 class ReadInfoType
00713 {
00714 public:
00715 
00716   typedef Node::ArrayType::ConstPointer NodeArrayPointer;
00717   typedef Element::ArrayType::ConstPointer ElementArrayPointer;
00718   typedef Material::ArrayType::ConstPointer MaterialArrayPointer;
00719 
00721   NodeArrayPointer m_node;
00722 
00724   ElementArrayPointer m_el;
00725 
00727   MaterialArrayPointer m_mat;
00728 
00730   ReadInfoType( NodeArrayPointer node_, ElementArrayPointer el_, MaterialArrayPointer mat_) :
00731     m_node(node_), m_el(el_), m_mat(mat_) {}
00732 };
00734 
00735 
00736 
00737 
00738 }} // end namespace itk::fem
00739 
00740 #endif // #ifndef __itkFEMElementBase_h
00741 

Generated at Mon Apr 14 12:23:02 2008 for ITK by doxygen 1.5.1 written by Dimitri van Heesch, © 1997-2000