Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkFEMSolverCrankNicolson.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkFEMSolverCrankNicolson.h,v $
00005   Language:  C++
00006   Date:      $Date: 2003/09/10 14:29:44 $
00007   Version:   $Revision: 1.17 $
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 <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:
00066 
00067 /*
00068  * helper initialization function before assembly but after generate GFN.
00069  */
00070   void InitializeForSolution(); 
00075   void AssembleKandM();            
00076 
00084   void AssembleFforTimeStep(int dim=0);
00085 
00089   void Solve();
00090 
00095   void AddToDisplacements(Float optimum=1.0);
00096   void AverageLastTwoDisplacements(Float t=0.5);
00097   void ZeroVector(int which=0);
00098   void PrintDisplacements(); 
00099   void PrintForce();
00100   
00102   inline void SetAlpha(Float a = 0.5) { m_alpha=a; }
00103 
00105   inline void SetDeltatT(Float T) { m_deltaT=T; }
00106 
00108   inline void SetRho(Float rho) { m_rho=rho;  }
00109 
00113   void RecomputeForceVector(unsigned int index);
00114 
00115   /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
00116   void FindBracketingTriplet(Float* a,Float* b,Float* c);
00120   Float GoldenSection(Float tol=0.01,unsigned int MaxIters=25);
00121   /* Brents method from Numerical Recipes. */
00122   Float BrentsMethod(Float tol=0.01,unsigned int MaxIters=25);
00123   Float EvaluateResidual(Float t=1.0);
00124   Float GetDeformationEnergy(Float t=1.0);
00125   inline Float GSSign(Float a,Float b) { return (b > 0.0 ? fabs(a) : -1.*fabs(a)); }
00126   inline Float GSMax(Float a,Float b) { return (a > b ? a : b); }
00127 
00128   void SetEnergyToMin(Float xmin);
00129   inline LinearSystemWrapper* GetLS(){ return m_ls;}
00130 
00131   Float GetCurrentMaxSolution() { return m_CurrentMaxSolution; }
00132 
00134   void PrintMinMaxOfSolution();
00139   SolverCrankNicolson() 
00140   { 
00141     m_deltaT=0.5; 
00142     m_rho=1.; 
00143     m_alpha=0.5;
00144     // BUG FIXME NOT SURE IF SOLVER IS USING VECTOR INDEX 1 FOR BCs
00145     ForceTIndex=0;                        // vector
00146     ForceTMinus1Index=2;                  // vector
00147     SolutionVectorTMinus1Index=3;         // vector
00148     DiffMatrixBySolutionTMinus1Index=4;   // vector
00149     ForceTotalIndex=5;                    // vector
00150     SolutionTIndex=0;                   // solution
00151     TotalSolutionIndex=1;               // solution
00152     SolutionTMinus1Index=2;       // solution
00153     SumMatrixIndex=0;                   // matrix
00154     DifferenceMatrixIndex=1;            // matrix
00155     m_CurrentMaxSolution=1.0;
00156   }
00157 
00158  
00159   ~SolverCrankNicolson() { }
00160 
00161   Float m_deltaT;
00162   Float m_rho;
00163   Float m_alpha;
00164   Float m_CurrentMaxSolution;
00165 
00166   unsigned int ForceTIndex;
00167   unsigned int ForceTotalIndex;
00168   unsigned int ForceTMinus1Index;
00169   unsigned int SolutionTIndex;
00170   unsigned int SolutionTMinus1Index;
00171   unsigned int SolutionVectorTMinus1Index;
00172   unsigned int TotalSolutionIndex;
00173   unsigned int DifferenceMatrixIndex;
00174   unsigned int SumMatrixIndex;
00175   unsigned int DiffMatrixBySolutionTMinus1Index;
00176   
00177 };
00178 
00179 
00180 
00181 
00182 }} // end namespace itk::fem
00183 
00184 #endif // #ifndef __itkFEMSolverCrankNicolson_h

Generated at Wed May 24 23:07:57 2006 for ITK by doxygen 1.3.5 written by Dimitri van Heesch, © 1997-2000