ITK  4.3.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 
101 
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
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
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 0;
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 = 0,
644  const MatrixType *pshapeD = 0) const;
645 
667  virtual void Jacobian(const VectorType & pt, MatrixType & J, const MatrixType *pshapeD = 0) const;
668 
678  virtual Float JacobianDeterminant(const VectorType & pt, const MatrixType *pJ = 0) const;
679 
691  virtual void JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixType *pJ = 0) 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;
729 
730 };
731 
732 }
733 } // end namespace itk::fem
734 
735 #endif // #ifndef __itkFEMElementBase_h
736