00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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 <vnl/algo/vnl_lsqr.h>
00034 #include <math.h>
00035
00036
00037 namespace itk {
00038 namespace fem {
00039
00040
00063 class SolverCrankNicolson : public Solver
00064 {
00065 public:
00067
00071 void InitializeForSolution();
00072
00077 void AssembleKandM();
00078
00086 void AssembleFforTimeStep(int dim=0);
00087
00091 void Solve();
00092
00097 void AddToDisplacements(Float optimum=1.0);
00098 void AverageLastTwoDisplacements(Float t=0.5);
00099 void ZeroVector(int which=0);
00100 void PrintDisplacements();
00101 void PrintForce();
00103
00105 inline void SetAlpha(Float a = 0.5) { m_alpha=a; }
00106
00108 inline void SetDeltatT(Float T) { m_deltaT=T; }
00109
00111 inline void SetRho(Float rho) { m_rho=rho; }
00112
00116 void RecomputeForceVector(unsigned int index);
00117
00118
00119 void FindBracketingTriplet(Float* a,Float* b,Float* c);
00120
00124 Float GoldenSection(Float tol=0.01,unsigned int MaxIters=25);
00125
00126 Float BrentsMethod(Float tol=0.01,unsigned int MaxIters=25);
00127 Float EvaluateResidual(Float t=1.0);
00128 Float GetDeformationEnergy(Float t=1.0);
00129 inline Float GSSign(Float a,Float b) { return (b > 0.0 ? vcl_fabs(a) : -1.*vcl_fabs(a)); }
00130 inline Float GSMax(Float a,Float b) { return (a > b ? a : b); }
00132
00133 void SetEnergyToMin(Float xmin);
00134 inline LinearSystemWrapper* GetLS(){ return m_ls;}
00135
00136 Float GetCurrentMaxSolution() { return m_CurrentMaxSolution; }
00137
00140 void PrintMinMaxOfSolution();
00141
00146 SolverCrankNicolson()
00147 {
00148 m_deltaT=0.5;
00149 m_rho=1.;
00150 m_alpha=0.5;
00151
00152 ForceTIndex=0;
00153 ForceTMinus1Index=2;
00154 SolutionVectorTMinus1Index=3;
00155 DiffMatrixBySolutionTMinus1Index=4;
00156 ForceTotalIndex=5;
00157 SolutionTIndex=0;
00158 TotalSolutionIndex=1;
00159 SolutionTMinus1Index=2;
00160 SumMatrixIndex=0;
00161 DifferenceMatrixIndex=1;
00162 m_CurrentMaxSolution=1.0;
00163 }
00164
00165
00166 ~SolverCrankNicolson() { }
00167
00168 Float m_deltaT;
00169 Float m_rho;
00170 Float m_alpha;
00171 Float m_CurrentMaxSolution;
00172
00173 unsigned int ForceTIndex;
00174 unsigned int ForceTotalIndex;
00175 unsigned int ForceTMinus1Index;
00176 unsigned int SolutionTIndex;
00177 unsigned int SolutionTMinus1Index;
00178 unsigned int SolutionVectorTMinus1Index;
00179 unsigned int TotalSolutionIndex;
00180 unsigned int DifferenceMatrixIndex;
00181 unsigned int SumMatrixIndex;
00182 unsigned int DiffMatrixBySolutionTMinus1Index;
00183
00184 };
00185
00186 }}
00187
00188 #endif // #ifndef __itkFEMSolverCrankNicolson_h
00189