Refactoring itk::FEM framework - V4

This page outlines the proposed changes to the itk::FEM framework. A number of these changes will break backwards compatibility.

Submission Of Code

http://review.source.kitware.com/#change,1993

Authors

The proposed changes for the itk::FEM framework are being proposed by the University of Iowa. This group is made up of investigators from the Departments of Radiology, Biomedical Engineering, and Psychiatry. Questions and comments regarding these changes should be sent to

• Vincent Magnotta (vincent-magnotta -at- uiowa.edu) - Radiology

Other team members include

• Nicole Grosland - Biomedical Engineering, MIMX
• Kiran Shivanna - MIMX
• Hans Johnson - Psychiatry
• Kent Williams - Psychiatry

Changes to the FEM Framework Classes

itk::Solver

1. Move the following class variables from public to protected:
   • node
   • el
   • load
   • mat
2. Add the following member functions:
   • bool AddNextElement(ElementType *element) – Add the next element to the mesh data structure. This will be used to replace the solver->el.push_back()
   • bool InsertElement(ElementType *element, unsigned long index) – Insert the element to the position within the mesh data structure.
   • bool AddNextNode(NodeType *node) – Add the next node to the mesh data structure. This will be used to replace the solver->node.push_back()
   • bool InsertNode(NodeType *node, unsigned long index) – Insert the node into the desired location within the mesh data structure.
   • bool AddNextMaterial(itk::fem::FEMP<itk::fem::Material>) – Add the next material type to the solver.
   • bool InsertMaterial(itk::fem::FEMP<itk::fem::Material>, int) – Insert the material property to the specified index in the solver material property array.
   • bool AddNextLoad(itk::fem::FEMP<itk::fem::Load>) – Add the next load to the solver.
   • bool InsertLoad(itk::fem::FEMP<itk::fem::Load>, int) – Insert the load to the specified index in the solver load array.
   • bool GetElement(unsigned long index, ElementType *element) – Get the specified element from the mesh data structure
   • bool GetNode(unsigned long index, NodeType *node) – Get the specified node from the mesh data structure
   • bool GetMaterial(itk::fem::FEMP<itk::fem::Material>, int) – Get the specified material from the solver.
   • bool GetLoad(itk::fem::FEMP<itk::fem::Load>, int) – Get the specified load from the solver.
   • unsigned int GetNumberOfNodes() – Returns the number of nodes stored in the solver.
   • unsigned int GetNumberOfElements()– Returns the number of element stored in the solver.
   • unsigned int GetNumberOfMaterials()– Returns the number of materials stored in the solver.
   • unsigned int GetNumberOfLoads() – Returns the number of loads stored in the solver
   • bool ClearMaterialArray() – Insert the load to the specified index in the solver load array.
   • bool ClearloadArray() – Insert the load to the specified index in the solver load array.
   • bool RemoveMaterial(int) – Remove the material property at the specified index from the solver Material array.
   • bool RemoveLoad(int) – Remove the load at the specified index from the solver load array.
   • void SetMesh(itk::Mesh) – Specify the entire mesh to be used by the solver.
   • itk::Mesh GetMesh() – Returns the mesh used by the solver.
   • MeshType::Pointer GetDeformedMesh() – Returns the resulting mesh after applying the loading and boundary conditions on the mesh.
   • itk::Array<double>::Pointer GetMeshStresses() – Returns the stresses for the mesh with each element of the array corresponding to an element in the mesh.
   • itk::Array<double>::Pointer GetMeshStrains()– Returns the strains for the mesh with each element of the array corresponding to an element in the mesh.

itk::fem::Load

itk::fem::LightObject

1. Move the following class variables from public to protected:
   • GN
2. Add public methods to get/set the class variables
   • void SetGlobalNumber(int) – Set the global number for the FEM object. The default value is -1.
   • int GetGlobalNumber() – Get the global number for the FEM object.

itk::fem::LoadBC

1. Move the following class variables from public to protected:
   • m_dof
   • m_element
   • m_value
   • CLID
2. The following class functions will be added, allowing the user to get and set the class variables:
   • void SetDegreesOfFreedom(unsigned int dof) – Set the degrees of freedom for the local element for which the boundary condition is being applied.
   • void SetElement(Element::ConstPointer element) – Set the element for which the boundary condition is being applied.
   • void SetValue(vnl_vector< Element::Float > value) – Set the boundary condition using a vector representation. This can be used to restrict motion of an element in a particular direction.
   • unsigned int GetDegreesOfFreedom() – Returns the local degrees of freedom for the element on which the boundary condition is being applied.
   • Element::ConstPointer GetElement() – Returns the element on which the boundary condition is being applied.
   • vnl_vector< Element::Float > GetValue() – Returns the assigned boundary condition.

itk::fem::LoadBCMFC

1. Move the following member variables from public to protected:
   • Index
   • lhs
   • rhs
2. Add class methods to get and set the class variables.
   • void SetIndex(int) – Set the index variable for the multi freedom displacement constraint. This is used internally by the itk::Fem::Solver.
   • void AddLeftHandSideTerm(LoadBCMFC::MFCTerm) – Add terms to the right hand side of the multi freedom displacement constraint. Used to replace loadBCMFC.lhs.push_back().
   • void AddRightHandSideTerm(vnl_vector< Element::Float >) – Set the right hand side of the linear equation that defines the constraints. Replaces loadBCMFC.rhs = a;
   • int GetIndex – Get the index variable for the multi freedom displacement constraint.
   • int GetNumberOfLeftHandSideTerms() – Returns the number of terms used to define the left hand side of the multi freedom displacement constraint.
   • itk::LoadBFMC GetLeftHandSideTerm(int index) – Returns the specified left hand side term.
   • int GetNumberOfRightHandSideTerms()– Returns the number of terms used to define the left hand side of the multi freedom displacement constraint.
   • Element::Float GetRightHandSideTerm(int index) – Returns the specified right hand side term.

itk::fem::LoadEdge

1. Move the following member variables from public to protected signatures:
   • m_Edge
   • m_Force
2. Provide public interface to get/set class variables
   • void SetEdge(int) – Set the edge to apply the desired force.
   • void SetForce(vnl_matrix< Float >) – Set the force to be applied to an edge.
   • int GetEdge() – Get the edge for the applied force.
   • vnl_matrix< Float > GetForce() – Get the force applied.

itk::fem::LoadGravConst

1. Move the following class variables from public to protected signatures.
   • Fg_value
2. Provide class methods to get/set the class variables
   • void SetForce(vnl_vector< Float >) – Set the constant force vector that exists for every point in space.
   • vnl_vector< Float > GetForce() – Return the constant force vector that exists for every point in space.

itk::fem::LoadLandmark

1. Move the following class variables from public to protected signatures.
   • eta
   • m_force
   • m_pt
   • m_Solution
   • m_Source
   • m_Target
2. Add the following public interface
   • void SetEta(double) – Set the square root of the variance.
   • double GetEta(double) – Get the square root of the variance.

itk::fem::LoadNode

1. Move the following class variables from public to protected signatures:
   • F
   • m_element
   • m_pt
2. Add the following class member functions
   • void SetForce(vnl_vector< Float >) – Set the applied force to the node.
   • void SetElement(Element::ConstPointer) – Set the element in the system that contains the degrees of freedom on which the force is applied.
   • void SetNode(unsigned int) – Set the point on which the force is being applied.
   • vnl_vector< Float > &GetForce() – Get the applied force.
   • Element::ConstPointer GetElement() – Get the element in the system that contains the degrees of freedom on which the force is applied.
   • Unsigned int GetNode() – Get the point on which the force is being applied.

itk::fem::LoadPoint

1. Move the following class variables from public to protected signatures:
   • Fp
   • point
2. Add public class methods to access these variables:
   • void SetForce(vnl_vector<Float>) – Set the force to be applied to the specified point location.
   • void SetPoint(vnl_vector<Float>) – Set the point where the force is applied in global coordinates.
   • vnl_vector<Float> & GetForce(vnl_vector<Float>) – Get the applied force.
   • vnl_vector<Float> & GetPoint (vnl_vector<Float>) – Get the point where the force is applied.

itk::fem::MaterialLinearElasticity

1. Move the following class variables from public to protected signatures:
   • A
   • E
   • h
   • I
   • nu
   • RhoC
2. Add the following public class member functions to get/set the class variables.
   • void SetCrossSectionalArea(double) – Set the cross sectional area (used to set variable A)
   • void SetYoungsModulus(double) – Set the Youngs Modulus (used to set the variable E)
   • void SetThickness(double) – Set the cross sectional thickness (used to set the variable h)
   • void SetMomentOfInertia(double) – Set the Moment of Inertia (used to set the variable I)
   • void SetPoissonsRatio(double) – Set the Poisson’s ratio (used to set the variable nu)
   • void SetDensityHeatProduct(double) – Set the Density - Heat Capacity product (used to set the variable RhoC)
   • double GetCrossSectionalArea() – Get the cross sectional area
   • double GetYoungsModulus() – Get the Youngs Modulus
   • double GetThickness() – Get the cross sectional thickness
   • double GetMomentOfInertia() – Get the Moment of Inertia
   • double GetPoissonsRatio() – Get the Poisson’s ratio
   • double GetDensityHeatProduct() – Get the Density - Heat Capacity product

Add New Object itk::FEMObject

Currently the itk::fem framework offers a complicated interface for specifying mesh to be used for analysis. We will create a new ITK object called itk::FEMObject which will define the entire FE model. The itk::FEMObject will inherit from itk::DataObject. This object will allow the user to completely specify the Nodes, Elements, Loads, and Materials for the problem. This decision was based on the decision to make a parallel structure to itk::Mesh without a complete refactoring of the mesh data structure. This also removes the specification of the FE simulation from the solver.

In order to support input/output, we have decided to use SpatialObjects. To support the new FEMObject, a new FEMSpatialObject was created. This also required that additions be added to the metaIO library to support this new spatial object type. The minor disadvantage of this implementation is that new element, load, and material types will require additions in multiple locations. We have decided to use the same FEM model specification that was used previously, but it is now encoded in a spatial object.

Migration Guide

The itk::fem::FEMObject was created to provide support for a first class object in ITK that will define the finite element problem being used for analysis. This used to be defined only in the itk::fem::Solver. All support for file I/O was removed from this class. The user would now define a FEMObject object and use this an input into the solver. Here is an outline on how to define a new problem.

 itk::fem::FEMObject::Pointer fem = itk::fem::FEMObject::New();
 .
 .
 .
 itk::fem::Solver::Pointer solver = itk::fem::Solver::New();
 solver->SetInput( fem );
 solver->Update( );
 
 itk::fem::FEMObject::Pointer result = solver->GetOutput();

The FEMObject has interfaces to define the nodes, elements, materials, and loads used to define the finite element problem.

• Nodes:

• Elements:

• Loads:
  1. AddNextLoad() - Add another load to the FEM problem. This will be added as the last load in the load container
  2. InsertLoad() - Add another load at the specified position in the load container
  3. GetNumberOfLoads() - Get number of loads in the FEM problem
  4. [1] - Get the entire load container
  5. GetLoadWithGlobalNumber() - Get the load with the specified global number. This is independent of the order in the load container, and is defined by the user when defining the FEM problem.

• Materials:

Update the itk::FEMRegistrationFilter

To make the FEMRegistrationFilter similar to the other registration filters the following changes will be made

1. Remove these public member functions. The application programmer will be responsible for loading the images and setting the parameters for the filter. The configuration file will be removed.
   • GetConfigFileName ()
   • GetFixedFile ()
   • GetMovingFile ()
   • GetResultsFileName ()
   • GetWriteDisplacements ()
   • ReadConfigFile (const char *)
   • SetConfigFileName (const char *f)
   • SetDisplacementsFile (const char *r)
   • SetFixedFile (const char *t)
   • SetLandmarkFile (const char *l)
   • SetMovingFile (const char *r)
   • SetResultsFile (const char *r)
   • SetResultsFileName (const char *f)
   • WriteDisplacementField (unsigned int index)
   • WriteDisplacementFieldMultiComponent ()
   • WriteWarpedImage (const char *fn)
2. Add these new public member functions
   • bool AddNextMovingLandmark(PointType) – Add another point for the moving image to the registration.
   • bool AddNextFixedLandmark(PointType) – Add another point from the fixed image to the registration.
   • void ClearMovingLandmarks() – Remove all of the existing landmarks for the moving image.
   • void ClearFixedLandmarks() – Remove all of the existing landmarks for the fixed image.
   • bool InsertMovingLandmark(PointType p, int index) – Insert the moving landmark at the specified index.
   • bool InsertFixedLandmark(PointType p, int index) – Insert the fixed landmark at the specified index.
   • bool RemoveMovingLandmark(int index)
   • bool RemoveFixedLandmark(int index)
   • void SetFEMesh() – Set the mesh to be used for the registration.
   • void GetFEMesh() – Get the mesh used in the registration.
   • void SetUserDefinedMesh(bool) – Specifies if the user will provide a custom mesh or if it will be generated based on the image.

File Formats

We are proposing to support I/O of the finite element models using Spatial Objects. There currently exists the itk::MeshSpatialObject, but this data structure does not completely contain all of the information required for a FEM. Therefore we are proposing to create a new object (itk::FEMObject) in ITK to hold the FE model information and use spatial objects to support the I/O. This is similar to what currently exists for itk::Mesh.

Other Formats

• Netgen supports gmsh format
• tetgen supports VTK
• VTK supports several additional formats

Data Structure

• Spatial Objects -
• ITK::Mesh -
• VTK Unstructured Grids -

Solvers

• SuperLU
• PETSc