ITK  4.6.0
Insight Segmentation and Registration Toolkit
itkFEMElementBase.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 
19 #ifndef __itkFEMElementBase_h
20 #define __itkFEMElementBase_h
21 
22 #include "itkFEMLightObject.h"
23 #include "itkFEMPArray.h"
24 #include "itkFEMMaterialBase.h"
25 #include "itkFEMSolution.h"
26 #include "itkVectorContainer.h"
27 #include "vnl/vnl_matrix.h"
28 #include "vnl/vnl_vector.h"
29 
30 #include <set>
31 #include <vector>
32 
33 namespace itk
34 {
35 namespace fem
36 {
72 class Element : public FEMLightObject
73 {
74 public:
76  typedef Element Self;
80 
82  itkTypeMacro(Element, FEMLightObject);
83 
87  typedef double Float;
88  typedef unsigned long ElementIdentifier;
89 
93  // FIXME - Remove FEMPArray Type and replace with VectorContainer version
96 
100  typedef vnl_matrix<Float> MatrixType;
101 
105  typedef vnl_vector<Float> VectorType;
106 
121 
126  typedef unsigned int DegreeOfFreedomIDType;
127 
133  enum { InvalidDegreeOfFreedomID = 0xffffffff };
134 
135 
147  class Node : public FEMLightObject
148  {
149  public:
150 
152  typedef Node Self;
156 
158  // itkNewMacro(Self);
159  static Pointer New(void)
160  {
163 
164  if( smartPtr.IsNull() )
165  {
166  smartPtr = static_cast<Pointer>(new Self);
167  }
168  smartPtr->UnRegister();
169  return smartPtr;
170  }
171 
173  itkTypeMacro(Node, FEMLightObject);
174 
177  virtual::itk::LightObject::Pointer CreateAnother(void) const ITK_OVERRIDE
178  {
180  Pointer copyPtr = Self::New();
182 
183  copyPtr->m_coordinates = this->m_coordinates;
184  copyPtr->m_dof = this->m_dof;
185  copyPtr->m_elements = this->m_elements;
186  copyPtr->SetGlobalNumber( this->GetGlobalNumber() );
187 
188  smartPtr = static_cast<Pointer>(copyPtr);
189 
190  return smartPtr;
191  }
192 
196  typedef double Float;
197 
202 
207  {
208  }
209 
214  {
215  this->ClearDegreesOfFreedom();
216  this->m_elements.clear();
217  }
219 
224  const VectorType & GetCoordinates(void) const
225  {
226  return m_coordinates;
227  }
228 
232  void SetCoordinates(const VectorType & coords)
233  {
234  m_coordinates = coords;
235  }
236 
241  {
242  if( i >= m_dof.size() )
243  {
245  }
246  return m_dof[i];
247  }
249 
253  void SetDegreeOfFreedom(unsigned int i, DegreeOfFreedomIDType dof) const
254  {
255  if( i >= m_dof.size() )
256  {
257  m_dof.resize(i + 1, InvalidDegreeOfFreedomID);
258  }
259  m_dof[i] = dof;
260  }
262 
263  virtual void ClearDegreesOfFreedom(void) const
264  {
265  m_dof.clear();
266  }
267 
268  public:
273  typedef std::set<Element *> SetOfElements;
275  protected:
276  virtual void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE
277  {
278  Superclass::PrintSelf(os, indent);
279  // os << indent << "DOF: " << this->m_dof << std::endl;
280  // os << indent << "Coordinates: " << this->m_coordinates << std::endl;
281  // os << indent << "Elements: " << this->m_elements << std::endl;
282  }
284 
285  private:
290 
295  mutable std::vector<DegreeOfFreedomIDType> m_dof;
296  }; // end class Node
297 
298 // ////////////////////////////////////////////////////////////////////////
299 /*
300  * Methods related to the physics of the problem.
301  */
302 
303  virtual VectorType GetStrainsAtPoint(const VectorType & pt, const Solution & sol, unsigned int index) const;
304 
305  virtual VectorType GetStressesAtPoint(const VectorType & pt, const VectorType & e, const Solution & sol,
306  unsigned int index) const;
307 
328  virtual void GetStiffnessMatrix(MatrixType & Ke) const;
329 
339  virtual Float GetElementDeformationEnergy(MatrixType & LocalSolution) const;
340 
352  virtual void GetMassMatrix(MatrixType & Me) const;
353 
365  virtual void GetLandmarkContributionMatrix(float eta, MatrixType & Le) const;
366 
374  virtual void GetStrainDisplacementMatrix(MatrixType & B, const MatrixType & shapeDgl) const = 0;
375 
381  virtual void GetMaterialMatrix(MatrixType & D) const = 0;
382 
394  virtual VectorType InterpolateSolution(const VectorType & pt,
395  const Solution & sol,
396  unsigned int solutionIndex = 0) const;
397 
411  virtual Float InterpolateSolutionN(const VectorType & pt, const Solution & sol, unsigned int f,
412  unsigned int solutionIndex = 0) const;
413 
420  DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int local_dof) const
421  {
422  if( local_dof > this->GetNumberOfDegreesOfFreedom() )
423  {
425  }
426  return this->GetNode(local_dof /
428  ->GetDegreeOfFreedom(local_dof % this->GetNumberOfDegreesOfFreedomPerNode() );
429  }
431 
449  {
450  return ITK_NULLPTR;
451  }
452 
462  {
463  }
464 
465  // ////////////////////////////////////////////////////////////////////////
492  virtual void GetIntegrationPointAndWeight(unsigned int i,
493  VectorType & pt,
494  Float & w,
495  unsigned int order = 0) const = 0;
496 
507  virtual unsigned int GetNumberOfIntegrationPoints(unsigned int order = 0) const = 0;
508 
517  enum { gaussMaxOrder = 10 };
518 
531 
538 
539 // ////////////////////////////////////////////////////////////////////////
540 /*
541  * Methods related to the geometry of an element
542  */
543 
549 
553  virtual unsigned int GetNumberOfNodes(void) const = 0;
554 
558  virtual NodeIDType GetNode(unsigned int n) const = 0;
559 
563  virtual void SetNode(unsigned int n, NodeIDType node) = 0;
564  virtual void SetNode(unsigned int n, Node::Pointer node)
565  {
566  this->SetNode(n,NodeIDType(node.GetPointer()));
567  }
568 
574  virtual const VectorType & GetNodeCoordinates(unsigned int n) const = 0;
575 
581  virtual VectorType GetGlobalFromLocalCoordinates(const VectorType & pt) const;
582 
589  virtual bool GetLocalFromGlobalCoordinates(const VectorType & globalPt, VectorType & localPt) const = 0;
590 
596  virtual unsigned int GetNumberOfSpatialDimensions() const = 0;
597 
605  virtual VectorType ShapeFunctions(const VectorType & pt) const = 0;
606 
622  virtual void ShapeFunctionDerivatives(const VectorType & pt, MatrixType & shapeD) const = 0;
623 
643  virtual void ShapeFunctionGlobalDerivatives(const VectorType & pt, MatrixType & shapeDgl, const MatrixType *pJ = ITK_NULLPTR,
644  const MatrixType *pshapeD = ITK_NULLPTR) const;
645 
667  virtual void Jacobian(const VectorType & pt, MatrixType & J, const MatrixType *pshapeD = ITK_NULLPTR) const;
668 
678  virtual Float JacobianDeterminant(const VectorType & pt, const MatrixType *pJ = ITK_NULLPTR) const;
679 
691  virtual void JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixType *pJ = ITK_NULLPTR) const;
692 
698  virtual unsigned int GetNumberOfDegreesOfFreedom(void) const
699  {
700  return this->GetNumberOfNodes() * this->GetNumberOfDegreesOfFreedomPerNode();
701  }
702 
706  virtual std::vector<std::vector<int> > GetEdgeIds(void) const
707  {
708  return this->m_EdgeIds;
709  }
710 
718  virtual unsigned int GetNumberOfDegreesOfFreedomPerNode(void) const = 0;
719 
721  virtual void PopulateEdgeIds(void) = 0;
722 
723 protected:
724 
725  // to store edge connectivity data
726  std::vector<std::vector<int> > m_EdgeIds;
727 
728  virtual void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE;
729 
730 };
731 
732 }
733 } // end namespace itk::fem
734 
735 #endif // #ifndef __itkFEMElementBase_h
DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int i) const
DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int local_dof) const
virtual NodeIDType GetNode(unsigned int n) const =0
static Pointer New(void)
virtual void Jacobian(const VectorType &pt, MatrixType &J, const MatrixType *pshapeD=ITK_NULLPTR) const
virtual void ShapeFunctionGlobalDerivatives(const VectorType &pt, MatrixType &shapeDgl, const MatrixType *pJ=ITK_NULLPTR, const MatrixType *pshapeD=ITK_NULLPTR) const
virtual void SetMaterial(Material::ConstPointer)
virtual Float JacobianDeterminant(const VectorType &pt, const MatrixType *pJ=ITK_NULLPTR) const
virtual void GetStrainDisplacementMatrix(MatrixType &B, const MatrixType &shapeDgl) const =0
virtual VectorType GetStressesAtPoint(const VectorType &pt, const VectorType &e, const Solution &sol, unsigned int index) const
LoadType::Pointer LoadPointer
int GetGlobalNumber() const
virtual void GetLandmarkContributionMatrix(float eta, MatrixType &Le) const
FEMLightObject LoadType
virtual unsigned int GetNumberOfDegreesOfFreedomPerNode(void) const =0
virtual void PrintSelf(std::ostream &os, Indent indent) const ITK_OVERRIDE
std::vector< DegreeOfFreedomIDType > m_dof
unsigned long ElementIdentifier
virtual std::vector< std::vector< int > > GetEdgeIds(void) const
static const Float gaussWeight[gaussMaxOrder+1][gaussMaxOrder]
std::vector< std::vector< int > > m_EdgeIds
virtual void PopulateEdgeIds(void)=0
std::set< Element * > SetOfElements
SmartPointer< const Self > ConstPointer
virtual Float GetElementDeformationEnergy(MatrixType &LocalSolution) const
virtual Float InterpolateSolutionN(const VectorType &pt, const Solution &sol, unsigned int f, unsigned int solutionIndex=0) const
static T::Pointer Create()
virtual void GetMassMatrix(MatrixType &Me) const
ObjectType * GetPointer() const
void SetCoordinates(const VectorType &coords)
static const double e
The base of the natural logarithm or Euler&#39;s number
Definition: itkMath.h:45
virtual unsigned int GetNumberOfIntegrationPoints(unsigned int order=0) const =0
Class that stores information required to define a node.
virtual void GetMaterialMatrix(MatrixType &D) const =0
vnl_matrix< Float > MatrixType
virtual void PrintSelf(std::ostream &os, Indent indent) const ITK_OVERRIDE
virtual void GetIntegrationPointAndWeight(unsigned int i, VectorType &pt, Float &w, unsigned int order=0) const =0
virtual void PrintSelf(std::ostream &os, Indent indent) const ITK_OVERRIDE
virtual void JacobianInverse(const VectorType &pt, MatrixType &invJ, const MatrixType *pJ=ITK_NULLPTR) const
unsigned int DegreeOfFreedomIDType
Provides functions to access the values of the solution vector.
virtual bool GetLocalFromGlobalCoordinates(const VectorType &globalPt, VectorType &localPt) const =0
VectorContainer< ElementIdentifier, Element::Pointer > ArrayType1
virtual unsigned int GetNumberOfNodes(void) const =0
Base class for all classes that define the FEM system.
virtual Material::ConstPointer GetMaterial(void) const
vnl_vector< Float > VectorType
SmartPointer< Self > Pointer
virtual unsigned int GetNumberOfSpatialDimensions() const =0
virtual void GetStiffnessMatrix(MatrixType &Ke) const
virtual ::itk::LightObject::Pointer CreateAnother(void) const ITK_OVERRIDE
SmartPointer< Self > Pointer
virtual void SetNode(unsigned int n, Node::Pointer node)
virtual VectorType GetStrainsAtPoint(const VectorType &pt, const Solution &sol, unsigned int index) const
bool IsNull() const
virtual VectorType ShapeFunctions(const VectorType &pt) const =0
const VectorType & GetCoordinates(void) const
Array for FEMP objects.
Definition: itkFEMPArray.h:40
Abstract base element class.
virtual void ClearDegreesOfFreedom(void) const
static const Float gaussPoint[gaussMaxOrder+1][gaussMaxOrder]
Node::ConstPointer NodeIDType
virtual void ShapeFunctionDerivatives(const VectorType &pt, MatrixType &shapeD) const =0
Define a front-end to the STL &quot;vector&quot; container that conforms to the IndexedContainerInterface.
SmartPointer< const Self > ConstPointer
Control indentation during Print() invocation.
Definition: itkIndent.h:49
virtual const VectorType & GetNodeCoordinates(unsigned int n) const =0
FEMPArray< Self > ArrayType
virtual void SetNode(unsigned int n, NodeIDType node)=0
virtual VectorType InterpolateSolution(const VectorType &pt, const Solution &sol, unsigned int solutionIndex=0) const
FEMPArray< Element > ArrayType
void SetDegreeOfFreedom(unsigned int i, DegreeOfFreedomIDType dof) const
FEMLightObject Superclass
virtual VectorType GetGlobalFromLocalCoordinates(const VectorType &pt) const
virtual unsigned int GetNumberOfDegreesOfFreedom(void) const