Main FEM solver. More...

#include <itkFEMSolver.h>

Inheritance diagram for itk::fem::Solver:

Collaboration diagram for itk::fem::Solver:

Public Types
typedef Element::ArrayType	ElementArray

typedef Element::Float	Float

typedef itk::Image < Element::ConstPointer, MaxGridDimensions >	InterpolationGridType

typedef Load::ArrayType	LoadArray

typedef Material::ArrayType	MaterialArray

typedef Node::ArrayType	NodeArray

typedef Element::VectorType	VectorType

Public Member Functions
	Solver ()

virtual	~Solver ()

void	ApplyBC (int dim=0, unsigned int matrix=0)

virtual void	AssembleElementMatrix (Element::Pointer e)

void	AssembleF (int dim=0)

void	AssembleK (void)

virtual void	AssembleLandmarkContribution (Element::Pointer e, float)

virtual void	Clear (void)

void	DecomposeK (void)

virtual void	FinalizeMatrixAfterAssembly (void)

void	GenerateGFN (void)

Float	GetDeformationEnergy (unsigned int SolutionIndex=0)

const Element *	GetElementAtPoint (const VectorType &pt) const

const InterpolationGridType *	GetInterpolationGrid (void) const

LinearSystemWrapper::Pointer	GetLinearSystemWrapper ()

unsigned int	GetNumberOfDegreesOfFreedom (void)

Float	GetSolution (unsigned int i, unsigned int which=0)

virtual Float	GetTimeStep (void) const

void	InitializeInterpolationGrid (const VectorType &size, const VectorType &bb1, const VectorType &bb2)

void	InitializeInterpolationGrid (const VectorType &size)

virtual void	InitializeLinearSystemWrapper (void)

virtual void	InitializeMatrixForAssembly (unsigned int N)

void	Read (std::istream &f)

void	SetLinearSystemWrapper (LinearSystemWrapper::Pointer ls)

virtual void	SetTimeStep (Float)

virtual void	Solve (void)

void	UpdateDisplacements (void)

void	Write (std::ostream &f)

Public Attributes
ElementArray	el

LoadArray	load

MaterialArray	mat

NodeArray	node

Static Public Attributes
static const unsigned int	MaxGridDimensions = 3

Protected Attributes
LinearSystemWrapper::Pointer	m_ls

unsigned int	NGFN

unsigned int	NMFC

Private Attributes
InterpolationGridType::Pointer	m_InterpolationGrid

LinearSystemWrapperVNL	m_lsVNL

Detailed Description

Main FEM solver.

This is the main class used for solving the FEM problems. It also stores all the objects that define the specific FEM problem. Normally there is one Solver object for each FEM problem.

Definition at line 41 of file itkFEMSolver.h.

Member Typedef Documentation

typedef Element::ArrayType itk::fem::Solver::ElementArray

Array that holds pointers to all elements. since we want to be able to manipulate the array we have to use special pointers

Definition at line 54 of file itkFEMSolver.h.

typedef Element::Float itk::fem::Solver::Float

Local float type

Definition at line 48 of file itkFEMSolver.h.

typedef itk::Image<Element::ConstPointer,MaxGridDimensions> itk::fem::Solver::InterpolationGridType

Type used to store interpolation grid

Definition at line 95 of file itkFEMSolver.h.

typedef Load::ArrayType itk::fem::Solver::LoadArray

Array that holds special pointers to all external loads

Definition at line 66 of file itkFEMSolver.h.

typedef Material::ArrayType itk::fem::Solver::MaterialArray

Array that holds pointers to the materials

Definition at line 72 of file itkFEMSolver.h.

typedef Node::ArrayType itk::fem::Solver::NodeArray

Array that holds special pointers to the nodes

Definition at line 60 of file itkFEMSolver.h.

typedef Element::VectorType itk::fem::Solver::VectorType

VectorType from the Element base class

Definition at line 78 of file itkFEMSolver.h.

Constructor & Destructor Documentation

itk::fem::Solver::Solver ( )

Default constructor sets Solver to use VNL linear system .

See Also: Solver::SetLinearSystemWrapper

Definition at line 42 of file itkFEMSolver.cxx.

References m_lsVNL, and SetLinearSystemWrapper().

virtual itk::fem::Solver::~Solver ( )

inlinevirtual

Virtual destructor

Definition at line 274 of file itkFEMSolver.h.

Member Function Documentation

void itk::fem::Solver::ApplyBC	(	int	dim = `0`,
		unsigned int	matrix = `0`
	)

Apply the boundary conditions to the system.

Note: This function must be called after AssembleK().

Parameters

matrix	Index of a matrix, to which the BCs should be applied (master stiffness matrix). Normally this is zero, but in derived classes many matrices may be used and this index must be specified.
dim	This is a parameter that can be passed to the function and is normally used with isotropic elements to specify the dimension in which the DOF is fixed.

Apply the boundary conditions to the system.

Store a temporary pointer to load object for later, so that we don't have to access it via the iterator

Apply boundary conditions in form of MFC loads.

We add the multi freedom constraints contribution to the master stiffness matrix using the lagrange multipliers. Basically we only change the last couple of rows and columns in K.

Apply essential boundary conditions

Definition at line 643 of file itkFEMSolver.cxx.

References itk::fem::LinearSystemWrapper::AddVectorValue(), itk::fem::LinearSystemWrapper::DestroyVector(), itk::fem::LinearSystemWrapper::GetColumnsOfNonZeroMatrixElementsInRow(), itk::fem::LinearSystemWrapper::GetMatrixValue(), itk::fem::LinearSystemWrapper::InitializeVector(), itk::fem::LinearSystemWrapper::IsVectorInitialized(), load, m_ls, NGFN, and itk::fem::LinearSystemWrapper::SetMatrixValue().

Referenced by itk::fem::SolverCrankNicolson::AssembleKandM(), and FinalizeMatrixAfterAssembly().

void itk::fem::Solver::AssembleElementMatrix ( Element::Pointer e )

virtual

Copy the element stiffness matrix into the correct position in the master stiffess matrix. Since more complex Solver classes may need to assemble many matrices and may also do some funky stuff to them, this function is virtual and can be overriden in a derived solver class.

Here we finaly update the corresponding element in the master stiffness matrix. We first check if element in Ke is zero, to prevent zeros from being allocated in sparse matrix.

Reimplemented in itk::fem::SolverHyperbolic.

Definition at line 347 of file itkFEMSolver.cxx.

References itk::fem::LinearSystemWrapper::AddMatrixValue(), itk::fem::Element::GetDegreeOfFreedom(), itk::fem::Element::GetNumberOfDegreesOfFreedom(), itk::fem::Element::GetStiffnessMatrix(), m_ls, and NGFN.

Referenced by AssembleK().

void itk::fem::Solver::AssembleF ( int dim = 0 )

Assemble the master force vector.

Parameters

dim	This is a parameter that can be passed to the function and is normally used with isotropic elements to specify the dimension for which the master force vector should be assembled.

Assemble the master force vector

Convert the external loads to the nodal loads and add them to the master force vector F.

Store a temporary pointer to load object for later, so that we don't have to access it via the iterator

Pass the vector to the solution to the Load object.

Here we only handle Nodal loads

If using the extra dim parameter, we can apply the force to different isotropic dimension.

FIXME: We assume that the implementation of force vector inside the LoadNode class is correct for given number of dimensions

Element loads...

If array of element pointers is not empty, we apply the load to all elements in that array

If the list of element pointers in load object is empty, we apply the load to all elements in a system.

Handle boundary conditions in form of MFC loads are handled next.

Handle essential boundary conditions.

If we got here, we were unable to handle that class of Load object. We do nothing...

Adjust the master force vector for essential boundary conditions as required.

Definition at line 389 of file itkFEMSolver.cxx.

References itk::fem::LinearSystemWrapper::AddVectorValue(), el, itk::fem::Element::GetDegreeOfFreedom(), itk::fem::Element::GetLoadVector(), itk::fem::Element::GetNumberOfDegreesOfFreedom(), itk::fem::LinearSystemWrapper::GetVectorValue(), itk::fem::LinearSystemWrapper::InitializeVector(), itk::fem::LinearSystemWrapper::IsVectorInitialized(), load, m_ls, NGFN, NMFC, and itk::fem::LinearSystemWrapper::SetVectorValue().

Referenced by itk::fem::SolverCrankNicolson::AssembleFforTimeStep(), and itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::IterativeSolve().

void itk::fem::Solver::AssembleK ( void )

Assemble the master stiffness matrix (also apply the MFCs to K)

Before we can start the assembly procedure, we need to know, how many boundary conditions if form of MFCs are there in a system.

Now we can assemble the master stiffness matrix from element stiffness matrices.

Since we're using the Lagrange multiplier method to apply the MFC, each constraint adds a new global DOF.

Step over all elements

Step over all the loads again to add the landmark contributions to the appropriate place in the stiffness matrix

Definition at line 241 of file itkFEMSolver.cxx.

References AssembleElementMatrix(), AssembleLandmarkContribution(), el, FinalizeMatrixAfterAssembly(), InitializeMatrixForAssembly(), load, NGFN, and NMFC.

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::MultiResSolve(), and itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::RunRegistration().

void itk::fem::Solver::AssembleLandmarkContribution	(	Element::Pointer	e,
		float	eta
	)

virtual

Add the contribution of the landmark-containing elements to the correct position in the master stiffess matrix. Since more complex Solver classes may need to assemble many matrices and may also do some funky stuff to them, this function is virtual and can be overriden in a derived solver class.

Here we finaly update the corresponding element in the master stiffness matrix. We first check if element in Le is zero, to prevent zeros from being allocated in sparse matrix.

Definition at line 311 of file itkFEMSolver.cxx.

References itk::fem::LinearSystemWrapper::AddMatrixValue(), itk::fem::Element::GetDegreeOfFreedom(), itk::fem::Element::GetLandmarkContributionMatrix(), itk::fem::Element::GetNumberOfDegreesOfFreedom(), m_ls, and NGFN.

Referenced by AssembleK().

void itk::fem::Solver::Clear ( void )

virtual

Cleans all data members, and initializes the solver to initial state.

Definition at line 47 of file itkFEMSolver.cxx.

References el, load, m_lsVNL, mat, NGFN, NMFC, node, and SetLinearSystemWrapper().

void itk::fem::Solver::DecomposeK ( void )

Decompose matrix using svd, qr, whatever ...

Decompose matrix using svd, qr, whatever ... if needed

Definition at line 586 of file itkFEMSolver.cxx.

virtual void itk::fem::Solver::FinalizeMatrixAfterAssembly ( void )

inlinevirtual

This function is called after the assebly has been completed. In this class it is only used to apply the BCs. You may however use it to perform other stuff in derived solver classes.

Reimplemented in itk::fem::SolverHyperbolic.

Definition at line 188 of file itkFEMSolver.h.

References ApplyBC().

Referenced by AssembleK().

void itk::fem::Solver::GenerateGFN ( void )

System solver functions. Call all six functions below (in listed order) to solve system. Assign a global freedom numbers to each DOF in a system. This must be done before any other solve function can be called.

Assign a global freedom number to each DOF in a system.

Assign new ID to every DOF in a system

Definition at line 186 of file itkFEMSolver.cxx.

References el, itk::fem::Element::InvalidDegreeOfFreedomID, NGFN, and node.

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::CreateMesh().

Solver::Float itk::fem::Solver::GetDeformationEnergy ( unsigned int SolutionIndex = 0 )

Get the total deformation energy using the chosen solution

Definition at line 620 of file itkFEMSolver.cxx.

References el, itk::fem::Solution::GetSolutionValue(), and m_ls.

const Element * itk::fem::Solver::GetElementAtPoint ( const VectorType & pt ) const

Returns the pointer to the element which contains global point pt.

Parameters

pt	Point in global coordinate system.

Note: Interpolation grid must be initializes before you can call this function.

Definition at line 883 of file itkFEMSolver.cxx.

References m_InterpolationGrid, and MaxGridDimensions.

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::InterpolateVectorField().

const InterpolationGridType* itk::fem::Solver::GetInterpolationGrid ( void ) const

inline

Returns pointer to interpolation grid, which is an itk::Image of pointers to Element objects. Normally you would use physical coordinates to get specific points (pointers to elements) from the image. You can then use the Elemenet::InterpolateSolution member function on the returned element to obtain the solution at this point.

Note: Physical coordinates in an image correspond to the global coordinate system in which the mesh (nodes) are.

Definition at line 131 of file itkFEMSolver.h.

References itk::SmartPointer< TObjectType >::GetPointer(), and m_InterpolationGrid.

LinearSystemWrapper::Pointer itk::fem::Solver::GetLinearSystemWrapper ( )

inline

Gets the LinearSystemWrapper object.

See Also: SetLinearSystemWrapper

Definition at line 297 of file itkFEMSolver.h.

References m_ls.

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::EnforceDiffeomorphism(), and itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::SampleVectorFieldAtNodes().

unsigned int itk::fem::Solver::GetNumberOfDegreesOfFreedom ( void )

inline

Definition at line 256 of file itkFEMSolver.h.

References NGFN.

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::MultiResSolve(), and itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::RunRegistration().

Float itk::fem::Solver::GetSolution	(	unsigned int	i,
		unsigned int	which = `0`
	)

inline

Definition at line 251 of file itkFEMSolver.h.

References itk::fem::Solution::GetSolutionValue(), and m_ls.

virtual Float itk::fem::Solver::GetTimeStep ( void ) const

inlinevirtual

Returns the time step used for dynamic problems.

Reimplemented in itk::fem::SolverHyperbolic.

Definition at line 308 of file itkFEMSolver.h.

void itk::fem::Solver::InitializeInterpolationGrid	(	const VectorType &	size,
		const VectorType &	bb1,
		const VectorType &	bb2
	)

Initialize the interpolation grid. The interpolation grid is used to find elements that containg specific points in a mesh. The interpolation grid stores pointers to elements for each point on a grid thereby providing a fast way (lookup table) to perform interpolation of results.

Note: Interpolation grid must be reinitialized each time a mesh changes.

Parameters

size	Vector that represents number of points on a grid in each dimension.
bb1	Lower limit of a bounding box of a grid.
bb2	Upper limit of a bounding box of a grid.

See Also: GetInterpolationGrid

Initialize the interpolation grid

Definition at line 758 of file itkFEMSolver.cxx.

References el, itk::ImageConstIterator< TImage >::GetIndex(), itk::ImageRegionConstIterator< TImage >::GoToBegin(), itk::ImageConstIterator< TImage >::IsAtEnd(), m_InterpolationGrid, MaxGridDimensions, itk::Image< TPixel, VImageDimension >::New(), and itk::ImageRegionIterator< TImage >::Set().

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::CreateMesh(), and InitializeInterpolationGrid().

void itk::fem::Solver::InitializeInterpolationGrid ( const VectorType & size )

inline

Same as InitializeInterpolationGrid(size, {0,0...}, size);

Definition at line 116 of file itkFEMSolver.h.

References InitializeInterpolationGrid().

void itk::fem::Solver::InitializeLinearSystemWrapper ( void )

virtual

Performs any initialization needed for LinearSystemWrapper object i.e. sets the maximum number of matrices and vectors.

Reimplemented in itk::fem::SolverHyperbolic.

Definition at line 72 of file itkFEMSolver.cxx.

References m_ls, itk::fem::LinearSystemWrapper::SetNumberOfMatrices(), itk::fem::LinearSystemWrapper::SetNumberOfSolutions(), and itk::fem::LinearSystemWrapper::SetNumberOfVectors().

Referenced by SetLinearSystemWrapper().

void itk::fem::Solver::InitializeMatrixForAssembly ( unsigned int N )

virtual

This function is called before assembling the matrices. You can override it in a derived class to account for special needs.

Parameters

N	Size of the matrix.

Reimplemented in itk::fem::SolverHyperbolic.

Definition at line 303 of file itkFEMSolver.cxx.

References itk::fem::LinearSystemWrapper::InitializeMatrix(), m_ls, and itk::fem::LinearSystemWrapper::SetSystemOrder().

Referenced by AssembleK().

void itk::fem::Solver::Read ( std::istream & f )

Reads the whole system (nodes, materials and elements) from input stream

If CreateFromStream returned 0, we're ok. That was the signal for the end of stream. Just continue reading... and consequently exit the do loop.

Find out what kind of object did we read from stream and store it in the appropriate array

If we got here, something strange was in the file...

Definition at line 84 of file itkFEMSolver.cxx.

References itk::fem::FEMLightObject::CreateFromStream(), el, load, mat, and node.

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::ApplyLoads(), and itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::CreateMesh().

void itk::fem::Solver::SetLinearSystemWrapper ( LinearSystemWrapper::Pointer ls )

Sets the LinearSystemWrapper object that will be used when solving the master equation. If this function is not called, a default VNL linear system representation will be used (class LinearSystemWrapperVNL).

Parameters

ls	Pointer to an object of class which is derived from LinearSystemWrapper.

Note: Once the LinearSystemWrapper object is changed, it is used until the member function SetLinearSystemWrapper is called again. Since LinearSystemWrapper object was created outside the Solver class, it should also be destroyed outside. Solver class will not destroy it when the Solver object is destroyed.

Change the LinearSystemWrapper object used to solve system of equations.

Definition at line 64 of file itkFEMSolver.cxx.

References InitializeLinearSystemWrapper(), and m_ls.

Referenced by Clear(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::MultiResSolve(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::RunRegistration(), and Solver().

virtual void itk::fem::Solver::SetTimeStep ( Float )

inlinevirtual

Sets the time step used for dynamic problems.

Parameters

dt	New time step.

Reimplemented in itk::fem::SolverHyperbolic.

Definition at line 315 of file itkFEMSolver.h.

void itk::fem::Solver::Solve ( void )

virtual

Solve for the displacement vector u. May be overriden in derived classes.

Solve for the displacement vector u

Reimplemented in itk::fem::SolverCrankNicolson, and itk::fem::SolverHyperbolic.

Definition at line 594 of file itkFEMSolver.cxx.

References itk::fem::LinearSystemWrapper::InitializeSolution(), itk::fem::LinearSystemWrapper::IsMatrixInitialized(), itk::fem::LinearSystemWrapper::IsVectorInitialized(), m_ls, and itk::fem::LinearSystemWrapper::Solve().

void itk::fem::Solver::UpdateDisplacements ( void )

Copy solution vector u to the corresponding nodal values, which are stored in node objects). This is standard post processing of the solution

Copy solution vector u to the corresponding nodal values, which are stored in node objects). This is standard post processing of the solution.

Definition at line 616 of file itkFEMSolver.cxx.

void itk::fem::Solver::Write ( std::ostream & f )

Writes everything (nodes, materials and elements) to output stream

Definition at line 155 of file itkFEMSolver.cxx.

References el, load, mat, and node.

Member Data Documentation

ElementArray itk::fem::Solver::el

Definition at line 55 of file itkFEMSolver.h.

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::ApplyLoads(), AssembleF(), AssembleK(), itk::fem::SolverCrankNicolson::AssembleKandM(), Clear(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::EvaluateResidual(), itk::fem::Generate2DRectilinearMesh(), itk::fem::Generate3DRectilinearMesh(), GenerateGFN(), GetDeformationEnergy(), InitializeInterpolationGrid(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::InterpolateVectorField(), Read(), and Write().

LoadArray itk::fem::Solver::load

Definition at line 67 of file itkFEMSolver.h.

Referenced by ApplyBC(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::ApplyImageLoads(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::ApplyLoads(), AssembleF(), itk::fem::SolverCrankNicolson::AssembleFforTimeStep(), AssembleK(), itk::fem::SolverCrankNicolson::AssembleKandM(), Clear(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::EnforceDiffeomorphism(), itk::fem::Generate2DRectilinearMesh(), Read(), and Write().

InterpolationGridType::Pointer itk::fem::Solver::m_InterpolationGrid

private

An Image of pointers to Element objects that represents a grid used for interpolation of solution. Each Pixel in an image is a pointer to an Element object in which that pixel is located.

Definition at line 345 of file itkFEMSolver.h.

Referenced by GetElementAtPoint(), GetInterpolationGrid(), and InitializeInterpolationGrid().

LinearSystemWrapper::Pointer itk::fem::Solver::m_ls

protected

Pointer to LinearSystemWrapper object.

Definition at line 331 of file itkFEMSolver.h.

LinearSystemWrapperVNL itk::fem::Solver::m_lsVNL

private

LinearSystemWrapperVNL object that is used by default in Solver class.

Definition at line 338 of file itkFEMSolver.h.

Referenced by Clear(), and Solver().

MaterialArray itk::fem::Solver::mat

Definition at line 73 of file itkFEMSolver.h.

Referenced by Clear(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::CreateMesh(), Read(), and Write().

const unsigned int itk::fem::Solver::MaxGridDimensions = 3

static

Since the itk::Image is templated over the number of dimensions, we have to know this at compile time. Solver class, however, can handle elements in any number of dimensions. In order to be able to use the Image, we choose the maximum number of space dimension that this function will be able to handle. Any unused dimensions are filled with zero.

For example: If a 2D node coordinates are {1.0,3.0} then the corresponding phisycal point in an image is {1.0,3.0,0.0};

Definition at line 90 of file itkFEMSolver.h.

Referenced by GetElementAtPoint(), and InitializeInterpolationGrid().

unsigned int itk::fem::Solver::NGFN

protected

Number of global degrees of freedom in a system

Definition at line 322 of file itkFEMSolver.h.

unsigned int itk::fem::Solver::NMFC

protected

Number of multi freedom constraints in a system. This member is set in a AssembleK function.

Definition at line 328 of file itkFEMSolver.h.

Referenced by AssembleF(), AssembleK(), itk::fem::SolverCrankNicolson::AssembleKandM(), Clear(), and itk::fem::SolverCrankNicolson::InitializeForSolution().

NodeArray itk::fem::Solver::node

Definition at line 61 of file itkFEMSolver.h.

Referenced by itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::ApplyLoads(), Clear(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::CreateMesh(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::EnforceDiffeomorphism(), itk::fem::Generate2DRectilinearMesh(), itk::fem::Generate3DRectilinearMesh(), GenerateGFN(), Read(), itk::fem::FEMRegistrationFilter< TMovingImage, TFixedImage >::SampleVectorFieldAtNodes(), and Write().

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Public Attributes

Static Public Attributes

Protected Attributes

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation