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: 2006/03/19 04:37:20 $
00007   Version:   $Revision: 1.18 $
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:
00067 
00068 /*
00069  * helper initialization function before assembly but after generate GFN.
00070  */
00071   void InitializeForSolution(); 
00076   void AssembleKandM();            
00077 
00085   void AssembleFforTimeStep(int dim=0);
00086 
00090   void Solve();
00091 
00096   void AddToDisplacements(Float optimum=1.0);
00097   void AverageLastTwoDisplacements(Float t=0.5);
00098   void ZeroVector(int which=0);
00099   void PrintDisplacements(); 
00100   void PrintForce();
00102 
00104   inline void SetAlpha(Float a = 0.5) { m_alpha=a; }
00105 
00107   inline void SetDeltatT(Float T) { m_deltaT=T; }
00108 
00110   inline void SetRho(Float rho) { m_rho=rho;  }
00111 
00115   void RecomputeForceVector(unsigned int index);
00116 
00117   /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
00118   void FindBracketingTriplet(Float* a,Float* b,Float* c);
00122   Float GoldenSection(Float tol=0.01,unsigned int MaxIters=25);
00123   /* Brents method from Numerical Recipes. */
00124   Float BrentsMethod(Float tol=0.01,unsigned int MaxIters=25);
00125   Float EvaluateResidual(Float t=1.0);
00126   Float GetDeformationEnergy(Float t=1.0);
00127   inline Float GSSign(Float a,Float b) { return (b > 0.0 ? vcl_fabs(a) : -1.*vcl_fabs(a)); }
00128   inline Float GSMax(Float a,Float b) { return (a > b ? a : b); }
00130 
00131   void SetEnergyToMin(Float xmin);
00132   inline LinearSystemWrapper* GetLS(){ return m_ls;}
00133 
00134   Float GetCurrentMaxSolution() { return m_CurrentMaxSolution; }
00135 
00137   void PrintMinMaxOfSolution();
00138 
00143   SolverCrankNicolson() 
00144   { 
00145     m_deltaT=0.5; 
00146     m_rho=1.; 
00147     m_alpha=0.5;
00148     // BUG FIXME NOT SURE IF SOLVER IS USING VECTOR INDEX 1 FOR BCs
00149     ForceTIndex=0;                        // vector
00150     ForceTMinus1Index=2;                  // vector
00151     SolutionVectorTMinus1Index=3;         // vector
00152     DiffMatrixBySolutionTMinus1Index=4;   // vector
00153     ForceTotalIndex=5;                    // vector
00154     SolutionTIndex=0;                   // solution
00155     TotalSolutionIndex=1;               // solution
00156     SolutionTMinus1Index=2;       // solution
00157     SumMatrixIndex=0;                   // matrix
00158     DifferenceMatrixIndex=1;            // matrix
00159     m_CurrentMaxSolution=1.0;
00160   }
00161 
00162  
00163   ~SolverCrankNicolson() { }
00164 
00165   Float m_deltaT;
00166   Float m_rho;
00167   Float m_alpha;
00168   Float m_CurrentMaxSolution;
00169 
00170   unsigned int ForceTIndex;
00171   unsigned int ForceTotalIndex;
00172   unsigned int ForceTMinus1Index;
00173   unsigned int SolutionTIndex;
00174   unsigned int SolutionTMinus1Index;
00175   unsigned int SolutionVectorTMinus1Index;
00176   unsigned int TotalSolutionIndex;
00177   unsigned int DifferenceMatrixIndex;
00178   unsigned int SumMatrixIndex;
00179   unsigned int DiffMatrixBySolutionTMinus1Index;
00180   
00181 };
00182 
00183 
00184 
00185 
00186 }} // end namespace itk::fem
00187 
00188 #endif // #ifndef __itkFEMSolverCrankNicolson_h
00189 

Generated at Sun Sep 23 12:32:59 2007 for ITK by doxygen 1.5.1 written by Dimitri van Heesch, © 1997-2000