00001 /*========================================================================= 00002 00003 Program: Insight Segmentation & Registration Toolkit 00004 Module: $RCSfile: itkFEMSolverCrankNicolson.h,v $ 00005 Language: C++ 00006 Date: $Date: 2010-02-26 05:28:25 $ 00007 Version: $Revision: 1.21 $ 00008 00009 Copyright (c) Insight Software Consortium. All rights reserved. 00010 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. 00011 00012 This software is distributed WITHOUT ANY WARRANTY; without even 00013 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00014 PURPOSE. See the above copyright notices for more information. 00015 00016 =========================================================================*/ 00017 00018 #ifndef __itkFEMSolverCrankNicolson_h 00019 #define __itkFEMSolverCrankNicolson_h 00020 00021 #include "itkFEMSolver.h" 00022 #include "itkFEMElementBase.h" 00023 #include "itkFEMMaterialBase.h" 00024 #include "itkFEMLoadBase.h" 00025 #include "itkFEMLinearSystemWrapperVNL.h" 00026 00027 #include "vnl/vnl_sparse_matrix.h" 00028 #include "vnl/vnl_matrix.h" 00029 #include "vnl/vnl_vector.h" 00030 #include "vnl/algo/vnl_svd.h" 00031 #include "vnl/algo/vnl_cholesky.h" 00032 #include <vnl/vnl_sparse_matrix_linear_system.h> 00033 #include <math.h> 00034 00035 00036 namespace itk { 00037 namespace fem { 00038 00039 00064 class SolverCrankNicolson : public Solver 00065 { 00066 public: 00068 00072 void InitializeForSolution(); 00073 00078 void AssembleKandM(); 00079 00087 void AssembleFforTimeStep(int dim=0); 00088 00092 void Solve(); 00093 00098 void AddToDisplacements(Float optimum=1.0); 00099 void AverageLastTwoDisplacements(Float t=0.5); 00100 void ZeroVector(int which=0); 00101 void PrintDisplacements(); 00102 void PrintForce(); 00104 00106 inline void SetAlpha(Float a = 0.5) { m_alpha=a; } 00107 00109 inline void SetDeltatT(Float T) { m_deltaT=T; } 00110 00112 inline void SetRho(Float rho) { m_rho=rho; } 00113 00117 void RecomputeForceVector(unsigned int index); 00118 00119 /* Finds a triplet that brackets the energy minimum. From Numerical Recipes.*/ 00120 void FindBracketingTriplet(Float* a,Float* b,Float* c); 00121 00125 Float GoldenSection(Float tol=0.01,unsigned int MaxIters=25); 00126 /* Brents method from Numerical Recipes. */ 00127 Float BrentsMethod(Float tol=0.01,unsigned int MaxIters=25); 00128 Float EvaluateResidual(Float t=1.0); 00129 Float GetDeformationEnergy(Float t=1.0); 00130 inline Float GSSign(Float a,Float b) { return (b > 0.0 ? vcl_fabs(a) : -1.*vcl_fabs(a)); } 00131 inline Float GSMax(Float a,Float b) { return (a > b ? a : b); } 00133 00134 void SetEnergyToMin(Float xmin); 00135 inline LinearSystemWrapper* GetLS(){ return m_ls;} 00136 00137 Float GetCurrentMaxSolution() { return m_CurrentMaxSolution; } 00138 00141 void PrintMinMaxOfSolution(); 00142 00147 SolverCrankNicolson() 00148 { 00149 m_deltaT=0.5; 00150 m_rho=1.; 00151 m_alpha=0.5; 00152 // BUG FIXME NOT SURE IF SOLVER IS USING VECTOR INDEX 1 FOR BCs 00153 ForceTIndex=0; // vector 00154 ForceTMinus1Index=2; // vector 00155 SolutionVectorTMinus1Index=3; // vector 00156 DiffMatrixBySolutionTMinus1Index=4; // vector 00157 ForceTotalIndex=5; // vector 00158 SolutionTIndex=0; // solution 00159 TotalSolutionIndex=1; // solution 00160 SolutionTMinus1Index=2; // solution 00161 SumMatrixIndex=0; // matrix 00162 DifferenceMatrixIndex=1; // matrix 00163 m_CurrentMaxSolution=1.0; 00164 } 00165 00166 00167 ~SolverCrankNicolson() { } 00168 00169 Float m_deltaT; 00170 Float m_rho; 00171 Float m_alpha; 00172 Float m_CurrentMaxSolution; 00173 00174 unsigned int ForceTIndex; 00175 unsigned int ForceTotalIndex; 00176 unsigned int ForceTMinus1Index; 00177 unsigned int SolutionTIndex; 00178 unsigned int SolutionTMinus1Index; 00179 unsigned int SolutionVectorTMinus1Index; 00180 unsigned int TotalSolutionIndex; 00181 unsigned int DifferenceMatrixIndex; 00182 unsigned int SumMatrixIndex; 00183 unsigned int DiffMatrixBySolutionTMinus1Index; 00184 00185 }; 00186 00187 }} // end namespace itk::fem 00188 00189 #endif // #ifndef __itkFEMSolverCrankNicolson_h 00190