ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkFEMSolverCrankNicolson.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 
00019 #ifndef __itkFEMSolverCrankNicolson_h
00020 #define __itkFEMSolverCrankNicolson_h
00021 
00022 #include "itkFEMSolver.h"
00023 #include "itkFEMElementBase.h"
00024 #include "itkFEMMaterialBase.h"
00025 #include "itkFEMLoadBase.h"
00026 #include "itkFEMLinearSystemWrapperVNL.h"
00027 
00028 #include "vnl/vnl_sparse_matrix.h"
00029 #include "vnl/vnl_matrix.h"
00030 #include "vnl/vnl_vector.h"
00031 #include "vnl/algo/vnl_svd.h"
00032 #include "vnl/algo/vnl_cholesky.h"
00033 #include <vnl/vnl_sparse_matrix_linear_system.h>
00034 #include <math.h>
00035 
00036 namespace itk
00037 {
00038 namespace fem
00039 {
00070 
00071 template <unsigned int TDimension = 3>
00072 class SolverCrankNicolson : public Solver<TDimension>
00073 {
00074 public:
00075   typedef SolverCrankNicolson      Self;
00076   typedef Solver<TDimension>       Superclass;
00077   typedef SmartPointer<Self>       Pointer;
00078   typedef SmartPointer<const Self> ConstPointer;
00079 
00081   itkNewMacro(Self);
00082 
00084   itkTypeMacro(SolverCrankNicolson, Solver<TDimension> );
00085 
00086   typedef Element::Float Float;
00087 
00091   itkSetMacro(UseMassMatrix, bool);
00092   itkGetMacro(UseMassMatrix, bool);
00094 
00098   itkGetConstMacro(Iterations, unsigned int);
00099 
00106   void ResetIterations(void)
00107   {
00108     m_Iterations = 0;
00109   }
00110 
00115   void AddToDisplacements(Float optimum = 1.0);
00116 
00117   void AverageLastTwoDisplacements(Float t = 0.5);
00118 
00119   void ZeroVector(int which = 0);
00120 
00121   void PrintDisplacements();
00122 
00123   void PrintForce();
00124 
00126   itkGetMacro(TotalSolutionIndex, unsigned int);
00127 
00129   itkGetMacro(SolutionTMinus1Index, unsigned int);
00130 
00132   itkSetMacro(Alpha, Float);
00133   itkGetMacro(Alpha, Float);
00135 
00137   itkSetMacro(Rho, Float);
00138   itkGetMacro(Rho, Float);
00140 
00142   virtual Float GetTimeStep(void) const
00143   {
00144     return m_TimeStep;
00145   }
00146 
00152   virtual void SetTimeStep(Float dt)
00153   {
00154     m_TimeStep = dt;
00155   }
00156 
00160   void RecomputeForceVector(unsigned int index);
00161 
00162 
00163   /* Finds a triplet that brackets the energy minimum.  From Numerical
00164     Recipes.*/
00165   void FindBracketingTriplet(Float *a, Float *b, Float *c);
00166 
00170   Float GoldenSection(Float tol = 0.01, unsigned int MaxIters = 25);
00171 
00172   /* Brents method from Numerical Recipes. */
00173   Float BrentsMethod(Float tol = 0.01, unsigned int MaxIters = 25);
00174 
00175   Float EvaluateResidual(Float t = 1.0);
00176 
00177   Float GetDeformationEnergy(Float t = 1.0);
00178 
00179   inline Float GSSign(Float a, Float b)
00180   {
00181     return b > 0.0 ? vcl_fabs(a) : -1. * vcl_fabs(a);
00182   }
00183   inline Float GSMax(Float a, Float b)
00184   {
00185     return a > b ? a : b;
00186   }
00187 
00188   void SetEnergyToMin(Float xmin);
00189 
00190   inline LinearSystemWrapper * GetLS()
00191   {
00192     return this->m_ls;
00193   }
00194 
00195   Float GetCurrentMaxSolution()
00196   {
00197     return m_CurrentMaxSolution;
00198   }
00199 
00202   void PrintMinMaxOfSolution();
00203 
00204 protected:
00205 
00206   SolverCrankNicolson();
00207   ~SolverCrankNicolson() { }
00208 
00211   void  GenerateData();
00212 
00216   virtual void RunSolver(void);
00217 
00221   void InitializeForSolution();
00222 
00227   void AssembleKandM();
00228 
00236   void AssembleFforTimeStep(int dim = 0);
00237 
00238   Float m_TimeStep;
00239   Float m_Rho;
00240   Float m_Alpha;
00241   Float m_CurrentMaxSolution;
00242 
00243   bool         m_UseMassMatrix;
00244   unsigned int m_Iterations;
00245 
00246   unsigned int m_ForceTIndex;
00247   unsigned int m_ForceTotalIndex;
00248   unsigned int m_ForceTMinus1Index;
00249   unsigned int m_SolutionTIndex;
00250   unsigned int m_SolutionTMinus1Index;
00251   unsigned int m_SolutionVectorTMinus1Index;
00252   unsigned int m_TotalSolutionIndex;
00253   unsigned int m_DifferenceMatrixIndex;
00254   unsigned int m_SumMatrixIndex;
00255   unsigned int m_DiffMatrixBySolutionTMinus1Index;
00256 private:
00257   SolverCrankNicolson(const Self &); // purposely not implemented
00258   void operator=(const Self &);      // purposely not implemented
00259 
00260 };
00261 }
00262 }  // end namespace itk::fem
00263 
00264 #ifndef ITK_MANUAL_INSTANTIATION
00265 #include "itkFEMSolverCrankNicolson.hxx"
00266 #endif
00267 
00268 #endif // #ifndef __itkFEMSolverCrankNicolson_h
00269