ITK  4.8.0
Insight Segmentation and Registration Toolkit
itkFEMSolverCrankNicolson.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 
19 #ifndef itkFEMSolverCrankNicolson_h
20 #define itkFEMSolverCrankNicolson_h
21 
22 #include "itkFEMSolver.h"
23 #include "itkFEMElementBase.h"
24 #include "itkFEMMaterialBase.h"
25 #include "itkFEMLoadBase.h"
27 
28 #include "vnl/vnl_sparse_matrix.h"
29 #include "vnl/vnl_matrix.h"
30 #include "vnl/vnl_vector.h"
31 #include "vnl/algo/vnl_svd.h"
32 #include "vnl/algo/vnl_cholesky.h"
33 #include <vnl/vnl_sparse_matrix_linear_system.h>
34 #include <cmath>
35 
36 namespace itk
37 {
38 namespace fem
39 {
70 
71 template <unsigned int TDimension = 3>
72 class SolverCrankNicolson : public Solver<TDimension>
73 {
74 public:
79 
82 
85 
87 
91  itkSetMacro(UseMassMatrix, bool);
92  itkGetMacro(UseMassMatrix, bool);
94 
98  itkGetConstMacro(Iterations, unsigned int);
99 
106  void ResetIterations(void)
107  {
108  m_Iterations = 0;
109  }
110 
115  void AddToDisplacements(Float optimum = 1.0);
116 
117  void AverageLastTwoDisplacements(Float t = 0.5);
118 
119  void ZeroVector(int which = 0);
120 
121  void PrintDisplacements();
122 
123  void PrintForce();
124 
126  itkGetMacro(TotalSolutionIndex, unsigned int);
127 
129  itkGetMacro(SolutionTMinus1Index, unsigned int);
130 
132  itkSetMacro(Alpha, Float);
133  itkGetMacro(Alpha, Float);
135 
137  itkSetMacro(Rho, Float);
138  itkGetMacro(Rho, Float);
140 
142  virtual Float GetTimeStep(void) const ITK_OVERRIDE
143  {
144  return m_TimeStep;
145  }
146 
152  virtual void SetTimeStep(Float dt) ITK_OVERRIDE
153  {
154  m_TimeStep = dt;
155  }
156 
160  void RecomputeForceVector(unsigned int index);
161 
162 
163  /* Finds a triplet that brackets the energy minimum. From Numerical
164  Recipes.*/
165  void FindBracketingTriplet(Float *a, Float *b, Float *c);
166 
170  Float GoldenSection(Float tol = 0.01, unsigned int MaxIters = 25);
171 
172  /* Brents method from Numerical Recipes. */
173  Float BrentsMethod(Float tol = 0.01, unsigned int MaxIters = 25);
174 
175  Float EvaluateResidual(Float t = 1.0);
176 
178 
179  inline Float GSSign(Float a, Float b)
180  {
181  return b > 0.0 ? std::fabs(a) : -1. * std::fabs(a);
182  }
183  inline Float GSMax(Float a, Float b)
184  {
185  return a > b ? a : b;
186  }
187 
188  void SetEnergyToMin(Float xmin);
189 
191  {
192  return this->m_ls;
193  }
194 
196  {
197  return m_CurrentMaxSolution;
198  }
199 
202  void PrintMinMaxOfSolution();
203 
204 protected:
205 
208 
211  void GenerateData() ITK_OVERRIDE;
212 
216  virtual void RunSolver(void) ITK_OVERRIDE;
217 
221  void InitializeForSolution();
222 
227  void AssembleKandM();
228 
236  void AssembleFforTimeStep(int dim = 0);
237 
242 
244  unsigned int m_Iterations;
245 
246  unsigned int m_ForceTIndex;
247  unsigned int m_ForceTotalIndex;
248  unsigned int m_ForceTMinus1Index;
249  unsigned int m_SolutionTIndex;
252  unsigned int m_TotalSolutionIndex;
254  unsigned int m_SumMatrixIndex;
256 
257 private:
258  SolverCrankNicolson(const Self &); // purposely not implemented
259  void operator=(const Self &); // purposely not implemented
260 
261 };
262 }
263 } // end namespace itk::fem
264 
265 #ifndef ITK_MANUAL_INSTANTIATION
266 #include "itkFEMSolverCrankNicolson.hxx"
267 #endif
268 
269 #endif // #ifndef itkFEMSolverCrankNicolson_h
virtual void SetTimeStep(Float dt) ITK_OVERRIDE
Float GoldenSection(Float tol=0.01, unsigned int MaxIters=25)
Light weight base class for most itk classes.
void SetEnergyToMin(Float xmin)
FEM Solver for time dependent problems; uses Crank-Nicolson implicit discretization scheme...
void RecomputeForceVector(unsigned int index)
Float GetDeformationEnergy(Float t=1.0)
Float EvaluateResidual(Float t=1.0)
itkGetConstMacro(Iterations, unsigned int)
itkTypeMacro(SolverCrankNicolson, Solver< TDimension >)
SmartPointer< const Self > ConstPointer
void GenerateData() ITK_OVERRIDE
void AddToDisplacements(Float optimum=1.0)
LinearSystemWrapper::Pointer m_ls
Definition: itkFEMSolver.h:405
Defines all functions required by Solver class to allocate, assemble and solve a linear system of equ...
void AssembleFforTimeStep(int dim=0)
void AverageLastTwoDisplacements(Float t=0.5)
itkGetMacro(UseMassMatrix, bool)
virtual Float GetTimeStep(void) const ITK_OVERRIDE
itkSetMacro(UseMassMatrix, bool)
FEM solver used to generate a solution for a FE formulation.
Definition: itkFEMSolver.h:69
Float BrentsMethod(Float tol=0.01, unsigned int MaxIters=25)
void FindBracketingTriplet(Float *a, Float *b, Float *c)
virtual void RunSolver(void) ITK_OVERRIDE
void ZeroVector(int which=0)