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: 2002/10/11 00:33:34 $
00007   Version:   $Revision: 1.14 $
00008 
00009   Copyright (c) 2002 Insight 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 <vxl/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 
00132   void PrintMinMaxOfSolution();
00137   SolverCrankNicolson() 
00138   { 
00139     m_deltaT=0.5; 
00140     m_rho=1.; 
00141     m_alpha=0.5;
00142     // BUG FIXME NOT SURE IF SOLVER IS USING VECTOR INDEX 1 FOR BCs
00143     ForceTIndex=0;                        // vector
00144     ForceTMinus1Index=2;                  // vector
00145     SolutionVectorTMinus1Index=3;         // vector
00146     DiffMatrixBySolutionTMinus1Index=4;   // vector
00147     ForceTotalIndex=5;                    // vector
00148     SolutionTIndex=0;                   // solution
00149     TotalSolutionIndex=1;               // solution
00150     SolutionTMinus1Index=2;       // solution
00151     SumMatrixIndex=0;                   // matrix
00152     DifferenceMatrixIndex=1;            // matrix
00153   }
00154 
00155  
00156   ~SolverCrankNicolson() { }
00157 
00158   Float m_deltaT;
00159   Float m_rho;
00160   Float m_alpha;
00161 
00162   unsigned int ForceTIndex;
00163   unsigned int ForceTotalIndex;
00164   unsigned int ForceTMinus1Index;
00165   unsigned int SolutionTIndex;
00166   unsigned int SolutionTMinus1Index;
00167   unsigned int SolutionVectorTMinus1Index;
00168   unsigned int TotalSolutionIndex;
00169   unsigned int DifferenceMatrixIndex;
00170   unsigned int SumMatrixIndex;
00171   unsigned int DiffMatrixBySolutionTMinus1Index;
00172   
00173 };
00174 
00175 
00176 
00177 
00178 }} // end namespace itk::fem
00179 
00180 #endif // #ifndef __itkFEMSolverCrankNicolson_h

Generated at Wed Mar 12 01:12:56 2003 for ITK by doxygen 1.2.15 written by Dimitri van Heesch, © 1997-2000