itkFEMSolver.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#ifndef itkFEMSolver_h
#define itkFEMSolver_h

#include "itkProcessObject.h"
#include "itkFEMObject.h"

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

#include "itkImage.h"

namespace itk
{
namespace fem
{
template <unsigned int VDimension = 3>
class Solver : public ProcessObject
{
public:

  typedef Solver                   Self;
  typedef ProcessObject            Superclass;
  typedef SmartPointer<Self>       Pointer;
  typedef SmartPointer<const Self> ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(Solver, ProcessObject);

  itkStaticConstMacro(FEMDimension, unsigned int, VDimension);
  itkStaticConstMacro(MaxDimensions, unsigned int, 3);

  typedef typename itk::fem::FEMObject<VDimension> FEMObjectType;
  typedef typename FEMObjectType::Pointer          FEMObjectPointer;
  typedef typename FEMObjectType::ConstPointer     FEMObjectConstPointer;
  typedef typename DataObject::Pointer             DataObjectPointer;

  typedef Element::Float           Float;
  typedef Element::VectorType      VectorType;
  typedef Element::Node::ArrayType NodeArray;
  typedef Element::ArrayType       ElementArray;
  typedef Load::ArrayType          LoadArray;
  typedef Material::ArrayType      MaterialArray;

  typedef typename itk::Image<Element::ConstPointer, VDimension> InterpolationGridType;
  typedef typename InterpolationGridType::Pointer                InterpolationGridPointerType;
  typedef typename InterpolationGridType::SizeType               InterpolationGridSizeType;
  typedef typename InterpolationGridType::RegionType             InterpolationGridRegionType;
  typedef typename InterpolationGridType::PointType              InterpolationGridPointType;
  typedef typename InterpolationGridType::SpacingType            InterpolationGridSpacingType;
  typedef typename InterpolationGridType::IndexType              InterpolationGridIndexType;
  typedef typename InterpolationGridType::DirectionType          InterpolationGridDirectionType;

  /*
   * Get/Set the Interpolation Grid Origin
   */
  itkSetMacro(Origin, InterpolationGridPointType);
  itkGetMacro(Origin, InterpolationGridPointType);

  /*
   * Get/Set the Interpolation Grid Spacing
   */
  itkSetMacro(Spacing, InterpolationGridSpacingType);
  itkGetMacro(Spacing, InterpolationGridSpacingType);

  /*
   * Get/Set the Interpolation Grid Region
   */
  itkSetMacro(Region, InterpolationGridRegionType);
  itkGetMacro(Region, InterpolationGridRegionType);

  /*
   * Get/Set the Interpolation Grid Direction
   */
  itkSetMacro(Direction, InterpolationGridDirectionType);
  itkGetMacro(Direction, InterpolationGridDirectionType);

  virtual Float GetTimeStep() const;

  virtual void SetTimeStep(Float dt);

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

  using Superclass::SetInput;
  virtual void SetInput( FEMObjectType *fem);

  virtual void SetInput( unsigned int, FEMObjectType * fem);

  FEMObjectType * GetInput();

  FEMObjectType * GetInput(unsigned int idx);

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

  Float GetDeformationEnergy(unsigned int SolutionIndex = 0);

  void SetLinearSystemWrapper(LinearSystemWrapper::Pointer ls);

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

  void InitializeInterpolationGrid(const InterpolationGridSizeType & size, const InterpolationGridPointType & bb1,
                                   const InterpolationGridPointType & bb2);

  void InitializeInterpolationGrid(const InterpolationGridSizeType & size)
  {
    InterpolationGridPointType bb1;

    bb1.Fill(0.0);

    InterpolationGridPointType bb2;
    for( unsigned int i = 0; i < FEMDimension; i++ )
    {
      bb2[i] = size[i] - 1.0;
    }
    InitializeInterpolationGrid(size, bb1, bb2);
  }

  void InitializeInterpolationGrid(const InterpolationGridRegionType & region,
                                   const InterpolationGridPointType & origin,
                                   const InterpolationGridSpacingType & spacing,
                                   const InterpolationGridDirectionType & direction);

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


  typedef ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType;
  using Superclass::MakeOutput;
  virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType) ITK_OVERRIDE;

  FEMObjectType * GetOutput();

  FEMObjectType * GetOutput(unsigned int idx);

protected:
  Solver();
  virtual ~Solver();
  virtual void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE;

  virtual void  GenerateData() ITK_OVERRIDE;


  // void GenerateGFN();

  void AssembleK();

  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::ConstPointer e, float);

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

  void AssembleF(int dim = 0);

  void DecomposeK();

  virtual void RunSolver();

  void UpdateDisplacements();

  void FillInterpolationGrid();

  virtual void InitializeLinearSystemWrapper();

  unsigned int m_NGFN;

  unsigned int m_NMFC;

  LinearSystemWrapper::Pointer m_ls;

  LinearSystemWrapperVNL m_lsVNL;

  InterpolationGridPointerType m_InterpolationGrid;

  FEMObjectPointer m_FEMObject;

private:
  Solver(const Self &) ITK_DELETE_FUNCTION;
  void operator=(const Self &) ITK_DELETE_FUNCTION;

  /*
   * Properties of the Interpolation Grid
   */
  InterpolationGridRegionType       m_Region;
  InterpolationGridPointType        m_Origin;
  InterpolationGridSpacingType      m_Spacing;
  InterpolationGridDirectionType    m_Direction;

};
}  // end namespace fem
}  // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFEMSolver.hxx"
#endif

#endif // #ifndef itkFEMSolver_h