itkFEMSolver.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMSolver.h,v $
  Language:  C++
  Date:      $Date: 2009-01-30 21:53:03 $
  Version:   $Revision: 1.28 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#ifndef __itkFEMSolver_h
#define __itkFEMSolver_h

#include "itkFEMElementBase.h"
#include "itkFEMMaterialBase.h"
#include "itkFEMLoadBase.h"

#include "itkFEMLinearSystemWrapper.h"
#include "itkFEMLinearSystemWrapperVNL.h"

#include "itkImage.h"

namespace itk {
namespace fem {

class Solver 
{
public:

  typedef Element::Float Float;

  typedef Element::ArrayType ElementArray;
  ElementArray el;

  typedef Node::ArrayType NodeArray;
  NodeArray node;

  typedef Load::ArrayType LoadArray;
  LoadArray load;

  typedef Material::ArrayType MaterialArray;
  MaterialArray mat;

  typedef Element::VectorType VectorType;

  itkStaticConstMacro(MaxGridDimensions, unsigned int, 3);

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

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

  void InitializeInterpolationGrid(const VectorType& size)
    {
    InitializeInterpolationGrid(size, VectorType(size.size(),0.0), size-1.0);
    }

  const InterpolationGridType * GetInterpolationGrid(void) const
  { return m_InterpolationGrid.GetPointer(); }

  const Element * GetElementAtPoint(const VectorType& pt) const;

  void Read( std::istream& f );
  
  void Write( std::ostream& f );

  virtual void Clear( void );


  void GenerateGFN( void );
            
  void AssembleK( void );

  virtual void InitializeMatrixForAssembly(unsigned int N);

  virtual void FinalizeMatrixAfterAssembly( void )
    {
    // Apply the boundary conditions to the K matrix
    this->ApplyBC();
    }

  virtual void AssembleElementMatrix(Element::Pointer e);

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

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

  void DecomposeK( void );

  virtual void Solve( void );

  void UpdateDisplacements( void );

  Float GetSolution(unsigned int i,unsigned int which=0)
    {
    return m_ls->GetSolutionValue(i,which);
    }

  unsigned int GetNumberOfDegreesOfFreedom( void )
    {
    return NGFN;
    }

  Float GetDeformationEnergy(unsigned int SolutionIndex=0);

public:
  Solver();

  virtual ~Solver() {}

  void SetLinearSystemWrapper(LinearSystemWrapper::Pointer ls);

  LinearSystemWrapper::Pointer GetLinearSystemWrapper() { return m_ls; }

  virtual void InitializeLinearSystemWrapper(void);

  virtual Float GetTimeStep( void ) const { return 0.0; }

  virtual void SetTimeStep(Float) {}

protected:

  unsigned int NGFN;

  unsigned int NMFC;

  LinearSystemWrapper::Pointer m_ls;

private:

  LinearSystemWrapperVNL m_lsVNL;

  InterpolationGridType::Pointer m_InterpolationGrid;

};

}} // end namespace itk::fem

#endif // #ifndef __itkFEMSolver_h