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: 2009-01-30 21:53:03 $
00007   Version:   $Revision: 1.19 $
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 
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   /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
00119   void FindBracketingTriplet(Float* a,Float* b,Float* c);
00120 
00124   Float GoldenSection(Float tol=0.01,unsigned int MaxIters=25);
00125   /* Brents method from Numerical Recipes. */
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     // BUG FIXME NOT SURE IF SOLVER IS USING VECTOR INDEX 1 FOR BCs
00152     ForceTIndex=0;                        // vector
00153     ForceTMinus1Index=2;                  // vector
00154     SolutionVectorTMinus1Index=3;         // vector
00155     DiffMatrixBySolutionTMinus1Index=4;   // vector
00156     ForceTotalIndex=5;                    // vector
00157     SolutionTIndex=0;                   // solution
00158     TotalSolutionIndex=1;               // solution
00159     SolutionTMinus1Index=2;       // solution
00160     SumMatrixIndex=0;                   // matrix
00161     DifferenceMatrixIndex=1;            // matrix
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 }} // end namespace itk::fem
00187 
00188 #endif // #ifndef __itkFEMSolverCrankNicolson_h
00189 

Generated at Tue Sep 15 02:55:35 2009 for ITK by doxygen 1.5.8 written by Dimitri van Heesch, © 1997-2000