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: 2010-02-26 05:28:25 $
00007   Version:   $Revision: 1.21 $
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 <math.h>
00034 
00035 
00036 namespace itk {
00037 namespace fem {
00038 
00039 
00064 class SolverCrankNicolson : public Solver
00065 {
00066 public:
00068 
00072   void InitializeForSolution(); 
00073 
00078   void AssembleKandM();
00079 
00087   void AssembleFforTimeStep(int dim=0);
00088 
00092   void Solve();
00093 
00098   void AddToDisplacements(Float optimum=1.0);
00099   void AverageLastTwoDisplacements(Float t=0.5);
00100   void ZeroVector(int which=0);
00101   void PrintDisplacements(); 
00102   void PrintForce();
00104 
00106   inline void SetAlpha(Float a = 0.5) { m_alpha=a; }
00107 
00109   inline void SetDeltatT(Float T) { m_deltaT=T; }
00110 
00112   inline void SetRho(Float rho) { m_rho=rho;  }
00113 
00117   void RecomputeForceVector(unsigned int index);
00118 
00119   /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
00120   void FindBracketingTriplet(Float* a,Float* b,Float* c);
00121 
00125   Float GoldenSection(Float tol=0.01,unsigned int MaxIters=25);
00126   /* Brents method from Numerical Recipes. */
00127   Float BrentsMethod(Float tol=0.01,unsigned int MaxIters=25);
00128   Float EvaluateResidual(Float t=1.0);
00129   Float GetDeformationEnergy(Float t=1.0);
00130   inline Float GSSign(Float a,Float b) { return (b > 0.0 ? vcl_fabs(a) : -1.*vcl_fabs(a)); }
00131   inline Float GSMax(Float a,Float b) { return (a > b ? a : b); }
00133 
00134   void SetEnergyToMin(Float xmin);
00135   inline LinearSystemWrapper* GetLS(){ return m_ls;}
00136 
00137   Float GetCurrentMaxSolution() { return m_CurrentMaxSolution; }
00138 
00141   void PrintMinMaxOfSolution();
00142 
00147   SolverCrankNicolson() 
00148     { 
00149     m_deltaT=0.5; 
00150     m_rho=1.; 
00151     m_alpha=0.5;
00152     // BUG FIXME NOT SURE IF SOLVER IS USING VECTOR INDEX 1 FOR BCs
00153     ForceTIndex=0;                        // vector
00154     ForceTMinus1Index=2;                  // vector
00155     SolutionVectorTMinus1Index=3;         // vector
00156     DiffMatrixBySolutionTMinus1Index=4;   // vector
00157     ForceTotalIndex=5;                    // vector
00158     SolutionTIndex=0;                   // solution
00159     TotalSolutionIndex=1;               // solution
00160     SolutionTMinus1Index=2;       // solution
00161     SumMatrixIndex=0;                   // matrix
00162     DifferenceMatrixIndex=1;            // matrix
00163     m_CurrentMaxSolution=1.0;
00164     }
00165 
00166  
00167   ~SolverCrankNicolson() { }
00168 
00169   Float m_deltaT;
00170   Float m_rho;
00171   Float m_alpha;
00172   Float m_CurrentMaxSolution;
00173 
00174   unsigned int ForceTIndex;
00175   unsigned int ForceTotalIndex;
00176   unsigned int ForceTMinus1Index;
00177   unsigned int SolutionTIndex;
00178   unsigned int SolutionTMinus1Index;
00179   unsigned int SolutionVectorTMinus1Index;
00180   unsigned int TotalSolutionIndex;
00181   unsigned int DifferenceMatrixIndex;
00182   unsigned int SumMatrixIndex;
00183   unsigned int DiffMatrixBySolutionTMinus1Index;
00184   
00185 };
00186 
00187 }} // end namespace itk::fem
00188 
00189 #endif // #ifndef __itkFEMSolverCrankNicolson_h
00190 

Generated at Fri Apr 16 18:19:22 2010 for ITK by doxygen 1.6.1 written by Dimitri van Heesch, © 1997-2000