itkFEMSolverCrankNicolson.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 itkFEMSolverCrankNicolson_h
#define itkFEMSolverCrankNicolson_h

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

#include "vnl/vnl_sparse_matrix.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector.h"
#include "vnl/algo/vnl_svd.h"
#include "vnl/algo/vnl_cholesky.h"
#include <vnl/vnl_sparse_matrix_linear_system.h>
#include <cmath>

namespace itk
{
namespace fem
{

template <unsigned int TDimension = 3>
class SolverCrankNicolson : public Solver<TDimension>
{
public:
  typedef SolverCrankNicolson      Self;
  typedef Solver<TDimension>       Superclass;
  typedef SmartPointer<Self>       Pointer;
  typedef SmartPointer<const Self> ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(SolverCrankNicolson, Solver<TDimension> );

  typedef Element::Float Float;

  itkSetMacro(UseMassMatrix, bool);
  itkGetMacro(UseMassMatrix, bool);

  itkGetConstMacro(Iterations, unsigned int);

  void ResetIterations(void)
  {
    m_Iterations = 0;
  }

  void AddToDisplacements(Float optimum = 1.0);

  void AverageLastTwoDisplacements(Float t = 0.5);

  void ZeroVector(int which = 0);

  void PrintDisplacements();

  void PrintForce();

  itkGetMacro(TotalSolutionIndex, unsigned int);

  itkGetMacro(SolutionTMinus1Index, unsigned int);

  itkSetMacro(Alpha, Float);
  itkGetMacro(Alpha, Float);

  itkSetMacro(Rho, Float);
  itkGetMacro(Rho, Float);

  virtual Float GetTimeStep(void) const ITK_OVERRIDE
  {
    return m_TimeStep;
  }

  virtual void SetTimeStep(Float dt) ITK_OVERRIDE
  {
    m_TimeStep = dt;
  }

  void RecomputeForceVector(unsigned int index);


  /* Finds a triplet that brackets the energy minimum.  From Numerical
    Recipes.*/
  void FindBracketingTriplet(Float *a, Float *b, Float *c);

  Float GoldenSection(Float tol = 0.01, unsigned int MaxIters = 25);

  /* Brents method from Numerical Recipes. */
  Float BrentsMethod(Float tol = 0.01, unsigned int MaxIters = 25);

  Float EvaluateResidual(Float t = 1.0);

  Float GetDeformationEnergy(Float t = 1.0);

  inline Float GSSign(Float a, Float b)
  {
    return b > 0.0 ? std::fabs(a) : -1. * std::fabs(a);
  }
  inline Float GSMax(Float a, Float b)
  {
    return a > b ? a : b;
  }

  void SetEnergyToMin(Float xmin);

  inline LinearSystemWrapper * GetLS()
  {
    return this->m_ls;
  }

  Float GetCurrentMaxSolution()
  {
    return m_CurrentMaxSolution;
  }

  void PrintMinMaxOfSolution();

protected:

  SolverCrankNicolson();
  ~SolverCrankNicolson() { }

  void  GenerateData() ITK_OVERRIDE;

  virtual void RunSolver(void) ITK_OVERRIDE;

  void InitializeForSolution();

  void AssembleKandM();

  void AssembleFforTimeStep(int dim = 0);

  Float m_TimeStep;
  Float m_Rho;
  Float m_Alpha;
  Float m_CurrentMaxSolution;

  bool         m_UseMassMatrix;
  unsigned int m_Iterations;

  unsigned int m_ForceTIndex;
  unsigned int m_ForceTotalIndex;
  unsigned int m_ForceTMinus1Index;
  unsigned int m_SolutionTIndex;
  unsigned int m_SolutionTMinus1Index;
  unsigned int m_SolutionVectorTMinus1Index;
  unsigned int m_TotalSolutionIndex;
  unsigned int m_DifferenceMatrixIndex;
  unsigned int m_SumMatrixIndex;
  unsigned int m_DiffMatrixBySolutionTMinus1Index;

private:
  SolverCrankNicolson(const Self &); // purposely not implemented
  void operator=(const Self &);      // purposely not implemented

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

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

#endif // #ifndef itkFEMSolverCrankNicolson_h