ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkFEMElementBase.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 
00019 #ifndef __itkFEMElementBase_h
00020 #define __itkFEMElementBase_h
00021 
00022 #include "itkFEMLightObject.h"
00023 #include "itkFEMPArray.h"
00024 #include "itkFEMMaterialBase.h"
00025 #include "itkFEMSolution.h"
00026 #include "itkVectorContainer.h"
00027 #include "vnl/vnl_matrix.h"
00028 #include "vnl/vnl_vector.h"
00029 
00030 #include <set>
00031 #include <vector>
00032 
00033 namespace itk
00034 {
00035 namespace fem
00036 {
00072 class Element : public FEMLightObject
00073 {
00074 public:
00076   typedef Element                  Self;
00077   typedef FEMLightObject           Superclass;
00078   typedef SmartPointer<Self>       Pointer;
00079   typedef SmartPointer<const Self> ConstPointer;
00080 
00082   itkTypeMacro(Element, FEMLightObject);
00083 
00087   typedef double        Float;
00088   typedef unsigned long ElementIdentifier;
00089 
00093   // FIXME - Remove FEMPArray Type and replace with VectorContainer version
00094   typedef FEMPArray<Element>                                   ArrayType;
00095   typedef VectorContainer<ElementIdentifier, Element::Pointer> ArrayType1;
00096 
00100   typedef vnl_matrix<Float> MatrixType;
00101 
00105   typedef vnl_vector<Float> VectorType;
00106 
00118   typedef FEMLightObject    LoadType;
00119   typedef LoadType::Pointer LoadPointer;
00121 
00126   typedef unsigned int DegreeOfFreedomIDType;
00127 
00133   enum { InvalidDegreeOfFreedomID = 0xffffffff };
00134 
00135 
00147   class Node : public FEMLightObject
00148   {
00149   public:
00150 
00152     typedef Node                     Self;
00153     typedef FEMLightObject           Superclass;
00154     typedef SmartPointer<Self>       Pointer;
00155     typedef SmartPointer<const Self> ConstPointer;
00156 
00158     // itkNewMacro(Self);
00159     static Pointer New(void)
00160       {
00161         Pointer smartPtr = ::itk::ObjectFactory<Self>::Create();
00163 
00164         if( smartPtr.IsNull() )
00165           {
00166           smartPtr = static_cast<Pointer>(new Self);
00167           }
00168         smartPtr->UnRegister();
00169         return smartPtr;
00170       }
00171 
00173     itkTypeMacro(Node, FEMLightObject);
00174 
00177     virtual::itk::LightObject::Pointer CreateAnother(void) const
00178       {
00179         ::itk::LightObject::Pointer smartPtr;
00180         Pointer copyPtr = Self::New();
00182 
00183         copyPtr->m_coordinates = this->m_coordinates;
00184         copyPtr->m_dof = this->m_dof;
00185         copyPtr->m_elements = this->m_elements;
00186         copyPtr->SetGlobalNumber( this->GetGlobalNumber() );
00187 
00188         smartPtr = static_cast<Pointer>(copyPtr);
00189 
00190         return smartPtr;
00191       }
00192 
00196     typedef double Float;
00197 
00201     typedef FEMPArray<Self> ArrayType;
00202 
00206     Node()
00207       {
00208       }
00209 
00213     ~Node()
00214       {
00215         this->ClearDegreesOfFreedom();
00216         this->m_elements.clear();
00217       }
00219 
00224     const VectorType & GetCoordinates(void) const
00225       {
00226         return m_coordinates;
00227       }
00228 
00232     void SetCoordinates(const VectorType & coords)
00233       {
00234         m_coordinates = coords;
00235       }
00236 
00240     DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int i) const
00241       {
00242         if( i >= m_dof.size() )
00243           {
00244           return InvalidDegreeOfFreedomID;
00245           }
00246         return m_dof[i];
00247       }
00249 
00253     void SetDegreeOfFreedom(unsigned int i, DegreeOfFreedomIDType dof) const
00254       {
00255         if( i >= m_dof.size() )
00256           {
00257           m_dof.resize(i + 1, InvalidDegreeOfFreedomID);
00258           }
00259         m_dof[i] = dof;
00260       }
00262 
00263     virtual void ClearDegreesOfFreedom(void) const
00264       {
00265         m_dof.clear();
00266       }
00267 
00268   public:
00273     typedef std::set<Element *> SetOfElements;
00274     mutable SetOfElements m_elements;
00275   protected:
00276     virtual void PrintSelf(std::ostream& os, Indent indent) const
00277       {
00278         Superclass::PrintSelf(os, indent);
00279         // os << indent << "DOF: " << this->m_dof << std::endl;
00280         // os << indent << "Coordinates: " << this->m_coordinates << std::endl;
00281         // os << indent << "Elements: " << this->m_elements << std::endl;
00282       }
00284 
00285   private:
00289     VectorType m_coordinates;
00290 
00295     mutable std::vector<DegreeOfFreedomIDType> m_dof;
00296   };  // end class Node
00297 
00298 // ////////////////////////////////////////////////////////////////////////
00299 /*
00300  * Methods related to the physics of the problem.
00301  */
00302 
00303   virtual VectorType GetStrainsAtPoint(const VectorType & pt, const Solution & sol, unsigned int index) const;
00304 
00305   virtual VectorType GetStressesAtPoint(const VectorType & pt, const VectorType & e, const Solution & sol,
00306                                         unsigned int index) const;
00307 
00328   virtual void GetStiffnessMatrix(MatrixType & Ke) const;
00329 
00339   virtual Float GetElementDeformationEnergy(MatrixType & LocalSolution) const;
00340 
00352   virtual void GetMassMatrix(MatrixType & Me) const;
00353 
00365   virtual void GetLandmarkContributionMatrix(float eta, MatrixType & Le) const;
00366 
00374   virtual void GetStrainDisplacementMatrix(MatrixType & B, const MatrixType & shapeDgl) const = 0;
00375 
00381   virtual void GetMaterialMatrix(MatrixType & D) const = 0;
00382 
00394   virtual VectorType InterpolateSolution(const VectorType & pt,
00395                                          const Solution & sol,
00396                                          unsigned int solutionIndex = 0) const;
00397 
00411   virtual Float InterpolateSolutionN(const VectorType & pt, const Solution & sol, unsigned int f,
00412                                      unsigned int solutionIndex = 0) const;
00413 
00420   DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int local_dof) const
00421     {
00422       if( local_dof > this->GetNumberOfDegreesOfFreedom() )
00423         {
00424         return InvalidDegreeOfFreedomID;
00425         }
00426       return this->GetNode(local_dof /
00427                            this->GetNumberOfDegreesOfFreedomPerNode() )
00428         ->GetDegreeOfFreedom(local_dof % this->GetNumberOfDegreesOfFreedomPerNode() );
00429     }
00431 
00448   virtual Material::ConstPointer GetMaterial(void) const
00449     {
00450       return 0;
00451     }
00452 
00461   virtual void SetMaterial(Material::ConstPointer)
00462     {
00463     }
00464 
00465   // ////////////////////////////////////////////////////////////////////////
00492   virtual void GetIntegrationPointAndWeight(unsigned int i,
00493                                             VectorType & pt,
00494                                             Float & w,
00495                                             unsigned int order = 0) const = 0;
00496 
00507   virtual unsigned int GetNumberOfIntegrationPoints(unsigned int order = 0) const = 0;
00508 
00517   enum { gaussMaxOrder = 10 };
00518 
00530   static const Float gaussPoint[gaussMaxOrder + 1][gaussMaxOrder];
00531 
00537   static const Float gaussWeight[gaussMaxOrder + 1][gaussMaxOrder];
00538 
00539 // ////////////////////////////////////////////////////////////////////////
00540 /*
00541  * Methods related to the geometry of an element
00542  */
00543 
00548   typedef Node::ConstPointer NodeIDType;
00549 
00553   virtual unsigned int GetNumberOfNodes(void) const = 0;
00554 
00558   virtual NodeIDType GetNode(unsigned int n) const = 0;
00559 
00563   virtual void SetNode(unsigned int n, NodeIDType node) = 0;
00564   virtual void SetNode(unsigned int n, Node::Pointer node)
00565     {
00566       this->SetNode(n,NodeIDType(node.GetPointer()));
00567     }
00568 
00574   virtual const VectorType & GetNodeCoordinates(unsigned int n) const = 0;
00575 
00581   virtual VectorType GetGlobalFromLocalCoordinates(const VectorType & pt) const;
00582 
00589   virtual bool GetLocalFromGlobalCoordinates(const VectorType & globalPt, VectorType & localPt) const = 0;
00590 
00596   virtual unsigned int GetNumberOfSpatialDimensions() const = 0;
00597 
00605   virtual VectorType ShapeFunctions(const VectorType & pt) const = 0;
00606 
00622   virtual void ShapeFunctionDerivatives(const VectorType & pt, MatrixType & shapeD) const = 0;
00623 
00643   virtual void ShapeFunctionGlobalDerivatives(const VectorType & pt, MatrixType & shapeDgl, const MatrixType *pJ = 0,
00644                                               const MatrixType *pshapeD = 0) const;
00645 
00667   virtual void Jacobian(const VectorType & pt, MatrixType & J, const MatrixType *pshapeD = 0) const;
00668 
00678   virtual Float JacobianDeterminant(const VectorType & pt, const MatrixType *pJ = 0) const;
00679 
00691   virtual void JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixType *pJ = 0) const;
00692 
00698   virtual unsigned int GetNumberOfDegreesOfFreedom(void) const
00699     {
00700       return this->GetNumberOfNodes() * this->GetNumberOfDegreesOfFreedomPerNode();
00701     }
00702 
00706   virtual std::vector<std::vector<int> > GetEdgeIds(void) const
00707     {
00708       return this->m_EdgeIds;
00709     }
00710 
00718   virtual unsigned int GetNumberOfDegreesOfFreedomPerNode(void) const = 0;
00719 
00721   virtual void PopulateEdgeIds(void) = 0;
00722 
00723 protected:
00724 
00725   // to store edge connectivity data
00726   std::vector<std::vector<int> > m_EdgeIds;
00727 
00728   virtual void PrintSelf(std::ostream& os, Indent indent) const;
00729 
00730 };
00731 
00732 }
00733 }  // end namespace itk::fem
00734 
00735 #endif // #ifndef __itkFEMElementBase_h
00736