Improving Initial Value Problem Solver

This commit improve different features in initial value problem solver
 * vtkRungeKutta4 now set xnext, even when going out of domain, the same way RK2 and RK45 do.
 * vtkInitialValueProblemSolvers now supports mulitple initialization, allowing to change
 the function set of a solver after initializing it.

Common/Math/vtkInitialValueProblemSolver.cxx

@@ -66,12 +66,14 @@ void vtkInitialValueProblemSolver::PrintSelf(ostream& os, vtkIndent indent)
 
 void vtkInitialValueProblemSolver::Initialize()
 {
-  if (!this->FunctionSet || this->Initialized)
+  if (!this->FunctionSet)
     {
     return;
     }
+  delete[] this->Vals;
   this->Vals =
     new double[this->FunctionSet->GetNumberOfIndependentVariables()];
+  delete[] this->Derivs;
   this->Derivs =
     new double[this->FunctionSet->GetNumberOfFunctions()];
   this->Initialized = 1;

Common/Math/vtkRungeKutta4.cxx

@@ -40,13 +40,14 @@ vtkRungeKutta4::~vtkRungeKutta4()
 void vtkRungeKutta4::Initialize()
 {
   this->vtkInitialValueProblemSolver::Initialize();
-  if (!this->Initialized)
+  if (!this->FunctionSet || !this->Initialized)
     {
     return;
     }
   // Allocate memory for temporary derivatives array
   for(int i=0; i<3; i++)
     {
+    delete[] this->NextDerivs[i];
     this->NextDerivs[i] =
       new double[this->FunctionSet->GetNumberOfFunctions()];
     }
@@ -98,6 +99,7 @@ int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext
     }
   else if ( !this->FunctionSet->FunctionValues(this->Vals, this->Derivs) )
     {
+    memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
     return OUT_OF_DOMAIN;
     }
 
@@ -110,6 +112,8 @@ int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext
   // 2
   if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[0]))
     {
+    memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
+    delTActual = delT/2.0; // we've been able to take half a step
     return OUT_OF_DOMAIN;
     }
 
@@ -122,6 +126,8 @@ int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext
   // 3
   if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[1]))
     {
+    memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
+    delTActual = delT/2.0; // we've been able to take half a step
     return OUT_OF_DOMAIN;
     }
 
@@ -134,6 +140,8 @@ int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext
   // 4
   if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[2]))
     {
+    memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
+    delTActual = delT; // we've been able to take a full step but couldn't finish the algorithm
     return OUT_OF_DOMAIN;
     }
 

Common/Math/vtkRungeKutta45.cxx

@@ -73,13 +73,14 @@ vtkRungeKutta45::~vtkRungeKutta45()
 void vtkRungeKutta45::Initialize()
 {
   this->vtkInitialValueProblemSolver::Initialize();
-  if (!this->Initialized)
+  if (!this->FunctionSet || !this->Initialized)
     {
     return;
     }
   // Allocate memory for temporary derivatives array
   for(int i=0; i<6; i++)
     {
+    delete[] this->NextDerivs[i];
     this->NextDerivs[i] =
       new double[this->FunctionSet->GetNumberOfFunctions()];
     }