00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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
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
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 };
00239
00241
00242
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) {}
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
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
00682 static INITClass Initializer_ElementNode(Element::Node::CLID());
00683
00684
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 }}
00722
00723 #endif // #ifndef __itkFEMElementBase_h
00724