Commit bbb0449e authored by Mathieu Westphal's avatar Mathieu Westphal Committed by Kitware Robot
Browse files

Merge topic 'LagrangianParticleTrackerSMP'

58e8bc6d Renable Divergence Theorem for vtkVoxel cells.
fce6c638 LagrangianParticleTracker SMP Implementation
566eeb20 Add a UserData parameters to vtkFunctionSet::FunctionValues
86955c82

 Adding a UserData parameter to  vtkInitialValueProblemSolver::ComputeNextStep
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: default avatarJoachim Pouderoux <joachim.pouderoux@kitware.com>
Merge-request: !5335
parents 9b2d712c 58e8bc6d
Pipeline #136712 failed with stage
in 0 seconds
......@@ -63,6 +63,7 @@ public:
*/
static vtkGenericInterpolatedVelocityField *New();
using Superclass::FunctionValues;
/**
* Evaluate the velocity field, f, at (x, y, z, t).
* For now, t is ignored.
......
......@@ -43,8 +43,11 @@ public:
* x and f have to point to valid double arrays of appropriate
* sizes obtained with GetNumberOfFunctions() and
* GetNumberOfIndependentVariables.
* If you inherit this class, make sure to reimplement at least one of the two
* FunctionValues signatures.
*/
virtual int FunctionValues(double* x, double* f) = 0;
virtual int FunctionValues(double* x, double* f) { return this->FunctionValues(x, f, nullptr); }
virtual int FunctionValues(double* x, double* f, void* vtkNotUsed(userData)) { return this->FunctionValues(x, f); }
/**
* Return the number of functions. Note that this is constant for
......
......@@ -68,35 +68,67 @@ public:
virtual int ComputeNextStep(double* xprev, double* xnext, double t,
double& delT, double maxError,
double& error)
{
return this->ComputeNextStep(xprev, xnext, t, delT, maxError, error, nullptr);
}
virtual int ComputeNextStep(double* xprev, double* xnext, double t,
double& delT, double maxError,
double& error, void* userData)
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, nullptr, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
virtual int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double maxError,
double& error)
{
return this->ComputeNextStep(xprev, dxprev, xnext, t, delT, maxError, error, nullptr);
}
virtual int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double maxError,
double& error, void* userData)
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, dxprev, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
virtual int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error)
{
return this->ComputeNextStep(xprev, xnext, t, delT, delTActual, minStep, maxStep, maxError, error, nullptr);
}
virtual int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error, void* userData)
{
return this->ComputeNextStep(xprev, nullptr, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
virtual int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error) = 0;
double maxError, double& error)
{
return this->ComputeNextStep(xprev, dxprev, xnext, t, delT, delTActual, minStep, maxStep, maxError, error, nullptr);
}
virtual int ComputeNextStep(double* vtkNotUsed(xprev), double* vtkNotUsed(dxprev), double* vtkNotUsed(xnext),
double vtkNotUsed(t), double& vtkNotUsed(delT), double& vtkNotUsed(delTActual),
double vtkNotUsed(minStep), double vtkNotUsed(maxStep),
double vtkNotUsed(maxError), double& vtkNotUsed(error), void* vtkNotUsed(userData)) { return 0; }
//@}
//@{
......@@ -139,7 +171,3 @@ private:
};
#endif
......@@ -26,7 +26,7 @@ vtkRungeKutta2::~vtkRungeKutta2() = default;
// Calculate next time step
int vtkRungeKutta2::ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double& delTActual,
double, double, double, double& error)
double, double, double, double& error, void* userData)
{
int i, numDerivs, numVals;
......@@ -61,7 +61,7 @@ int vtkRungeKutta2::ComputeNextStep(double* xprev, double* dxprev, double* xnext
this->Derivs[i] = dxprev[i];
}
}
else if ( !this->FunctionSet->FunctionValues(this->Vals, this->Derivs) )
else if ( !this->FunctionSet->FunctionValues(this->Vals, this->Derivs, userData) )
{
memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
return OUT_OF_DOMAIN;
......@@ -75,7 +75,7 @@ int vtkRungeKutta2::ComputeNextStep(double* xprev, double* dxprev, double* xnext
this->Vals[numVals-1] = t + delT/2.0;
// Obtain the derivatives at x_i + dt/2 * dx_i
if (!this->FunctionSet->FunctionValues(this->Vals, this->Derivs))
if (!this->FunctionSet->FunctionValues(this->Vals, this->Derivs, userData))
{
memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
delTActual = delT/2.0; // we've only taken half of a time step
......
......@@ -42,6 +42,7 @@ public:
*/
static vtkRungeKutta2 *New();
using Superclass::ComputeNextStep;
//@{
/**
* Given initial values, xprev , initial time, t and a requested time
......@@ -58,36 +59,36 @@ public:
*/
int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, nullptr, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, dxprev, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
return this->ComputeNextStep(xprev, nullptr, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error) override;
double maxError, double& error, void* userData) override;
//@}
protected:
......@@ -100,11 +101,4 @@ private:
#endif
// VTK-HeaderTest-Exclude: vtkRungeKutta2.h
......@@ -60,7 +60,7 @@ void vtkRungeKutta4::Initialize()
// (Addison Wesley)
int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double& delTActual,
double, double, double, double& error)
double, double, double, double& error, void* userData)
{
int i, numDerivs, numVals;
......@@ -97,7 +97,7 @@ int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext
this->Derivs[i] = dxprev[i];
}
}
else if ( !this->FunctionSet->FunctionValues(this->Vals, this->Derivs) )
else if ( !this->FunctionSet->FunctionValues(this->Vals, this->Derivs, userData) )
{
memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
return OUT_OF_DOMAIN;
......@@ -110,7 +110,7 @@ int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext
this->Vals[numVals-1] = t + delT/2.0;
// 2
if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[0]))
if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[0], userData))
{
memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
delTActual = delT/2.0; // we've been able to take half a step
......@@ -124,7 +124,7 @@ int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext
this->Vals[numVals-1] = t + delT/2.0;
// 3
if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[1]))
if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[1], userData))
{
memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
delTActual = delT/2.0; // we've been able to take half a step
......@@ -138,7 +138,7 @@ int vtkRungeKutta4::ComputeNextStep(double* xprev, double* dxprev, double* xnext
this->Vals[numVals-1] = t + delT;
// 4
if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[2]))
if (!this->FunctionSet->FunctionValues(this->Vals, this->NextDerivs[2], userData))
{
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
......
......@@ -43,7 +43,7 @@ public:
*/
static vtkRungeKutta4 *New();
using Superclass::ComputeNextStep;
//@{
/**
* Given initial values, xprev , initial time, t and a requested time
......@@ -60,36 +60,36 @@ public:
*/
int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, nullptr, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, dxprev, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
return this->ComputeNextStep(xprev, nullptr, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error) override;
double maxError, double& error, void* userData) override;
//@}
protected:
......@@ -105,11 +105,3 @@ private:
};
#endif
......@@ -91,7 +91,7 @@ int vtkRungeKutta45::ComputeNextStep(double* xprev, double* dxprev,
double* xnext, double t, double& delT,
double& delTActual,
double minStep, double maxStep,
double maxError, double& estErr )
double maxError, double& estErr, void* userData)
{
estErr = VTK_DOUBLE_MAX;
......@@ -112,7 +112,7 @@ int vtkRungeKutta45::ComputeNextStep(double* xprev, double* dxprev,
if ( ((minStep == absDT) && (maxStep == absDT)) ||
(maxError <= 0.0) )
{
int retVal = this->ComputeAStep(xprev, dxprev, xnext, t, delT, delTActual, estErr);
int retVal = this->ComputeAStep(xprev, dxprev, xnext, t, delT, delTActual, estErr, userData);
return retVal;
}
else if ( minStep > maxStep )
......@@ -127,7 +127,7 @@ int vtkRungeKutta45::ComputeNextStep(double* xprev, double* dxprev,
while ( estErr > maxError )
{
if ((retVal =
this->ComputeAStep(xprev, dxprev, xnext, t, delT, delTActual, estErr)))
this->ComputeAStep(xprev, dxprev, xnext, t, delT, delTActual, estErr, userData)))
{
return retVal;
}
......@@ -189,7 +189,7 @@ int vtkRungeKutta45::ComputeNextStep(double* xprev, double* dxprev,
if (shouldBreak)
{
if ( (retVal =
this->ComputeAStep(xprev, dxprev, xnext, t, delT, delTActual, estErr)) )
this->ComputeAStep(xprev, dxprev, xnext, t, delT, delTActual, estErr, userData)) )
{
return retVal;
}
......@@ -204,7 +204,7 @@ int vtkRungeKutta45::ComputeNextStep(double* xprev, double* dxprev,
// Calculate next time step
int vtkRungeKutta45::ComputeAStep(
double* xprev, double* dxprev, double* xnext, double t, double& delT,
double& actualDelT, double& error)
double& actualDelT, double& error, void* userData)
{
int i, j, k, numDerivs, numVals;
......@@ -240,7 +240,7 @@ int vtkRungeKutta45::ComputeAStep(
}
}
else if ( !this->FunctionSet->FunctionValues(this->Vals,
this->NextDerivs[0]) )
this->NextDerivs[0], userData))
{
for(i=0; i<numVals-1; i++)
{
......@@ -266,7 +266,7 @@ int vtkRungeKutta45::ComputeAStep(
this->Vals[numVals-1] = t + delT*A[i-1];
if ( !this->FunctionSet->FunctionValues(this->Vals,
this->NextDerivs[i]) )
this->NextDerivs[i], userData))
{
for(int l = 0; l < numVals - 1; l++)
{
......
......@@ -49,6 +49,7 @@ public:
*/
static vtkRungeKutta45 *New();
using Superclass::ComputeNextStep;
//@{
/**
* Given initial values, xprev , initial time, t and a requested time
......@@ -74,36 +75,36 @@ public:
*/
int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, nullptr, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, dxprev, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error) override
double maxError, double& error, void* userData) override
{
return this->ComputeNextStep(xprev, nullptr, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
minStep, maxStep, maxError, error, userData);
}
int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error) override;
double maxError, double& error, void* userData) override;
//@}
protected:
......@@ -121,7 +122,7 @@ protected:
double* NextDerivs[6];
int ComputeAStep(double* xprev, double* dxprev, double* xnext, double t,
double& delT, double& delTActual, double& error);
double& delT, double& delTActual, double& error, void* userData);
private:
vtkRungeKutta45(const vtkRungeKutta45&) = delete;
......
......@@ -128,9 +128,6 @@ int TestLagrangianIntegrationModel(int, char*[])
int nvar = odeWavelet->GetNumberOfIndependentVariables();
int seedIdx = 13;
vtkLagrangianParticle* part =
new vtkLagrangianParticle(nvar, seedIdx, seedIdx, 0, 0, pd);
odeWavelet->SetInputArrayToProcess(2, 0, 0,
vtkDataObject::FIELD_ASSOCIATION_CELLS, "");
odeWavelet->SetInputArrayToProcess(3, 0, 0,
......@@ -143,7 +140,6 @@ int TestLagrangianIntegrationModel(int, char*[])
vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDiameter");
odeWavelet->SetInputArrayToProcess(7, 1, 0,
vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDensity");
odeWavelet->SetCurrentParticle(part);
vtkNew<vtkCellLocator> locator;
odeWavelet->SetLocator(locator);
odeWavelet->AddDataSet(wavelet->GetOutput());
......@@ -153,7 +149,6 @@ int TestLagrangianIntegrationModel(int, char*[])
if (odeWavelet->GetLocator() != locator)
{
std::cerr << "Problem with locator" << std::endl;
delete part;
return EXIT_FAILURE;
}
......@@ -161,32 +156,30 @@ int TestLagrangianIntegrationModel(int, char*[])
if (!odeWavelet->GetUseInitialIntegrationTime())
{
std::cerr << "Problems with UseInitialIntegrationTime" << std::endl;
delete part;
return EXIT_FAILURE;
}
odeWavelet->UseInitialIntegrationTimeOff();
if (odeWavelet->GetUseInitialIntegrationTime())
{
std::cerr << "Problems with UseInitialIntegrationTime" << std::endl;
delete part;
return EXIT_FAILURE;
}
odeWavelet->UseInitialIntegrationTimeOn();
if (!odeWavelet->GetUseInitialIntegrationTime())
{
std::cerr << "Problems with UseInitialIntegrationTime" << std::endl;
delete part;
return EXIT_FAILURE;
}
odeWavelet->SetUseInitialIntegrationTime(false);
vtkLagrangianParticle part(nvar, seedIdx, seedIdx, 0, 0, pd, odeWavelet->GetWeightsSize());
odeWavelet->InitializeVariablesParticleData(pd);
odeWavelet->InsertVariablesParticleData(part, pd, 0);
odeWavelet->InitializeParticle(part);
if (odeWavelet->CheckFreeFlightTermination(part))
odeWavelet->InsertVariablesParticleData(&part, pd, 0);
odeWavelet->InitializeParticle(&part);
if (odeWavelet->CheckFreeFlightTermination(&part))
{
std::cerr << "CheckFreeFlightTermination should not return true with a matida model" << std::endl;
delete part;
return EXIT_FAILURE;
}
......@@ -194,68 +187,58 @@ int TestLagrangianIntegrationModel(int, char*[])
if (!odeWavelet->GetNonPlanarQuadSupport())
{
std::cerr << "Something went wrong with NonPlanarQuadSupport" << std::endl;
delete part;
return EXIT_FAILURE;
}
if (odeWavelet->GetWeightsSize() != 8)
{
std::cerr << "Incorrect Weights Size" << std::endl;
delete part;
return EXIT_FAILURE;
}
odeWavelet->ParallelManualShift(part);
odeWavelet->ParallelManualShift(&part);
vtkPolyData* tmpPd = nullptr;
vtkDataObject* tmpDo = nullptr;
if (!odeWavelet->FinalizeOutputs(tmpPd, tmpDo))
{
std::cerr << "FinalizeOutputs should be doing nothing and return true with matida model" << std::endl;
delete part;
return EXIT_FAILURE;
}
if (odeWavelet->GetSeedArrayNames()->GetNumberOfValues() != 4)
{
std::cerr << "Unexpected number of seed array names" << std::endl;
delete part;
return EXIT_FAILURE;
}
if (odeWavelet->GetSeedArrayComps()->GetNumberOfValues() != 4)
{
std::cerr << "Unexpected number of seed array comps" << std::endl;
delete part;
return EXIT_FAILURE;
}
if (odeWavelet->GetSeedArrayTypes()->GetNumberOfValues() != 4)
{
std::cerr << "Unexpected number of seed array type" << std::endl;
delete part;
return EXIT_FAILURE;
}
if (odeWavelet->GetSurfaceArrayNames()->GetNumberOfValues() != 1)
{
std::cerr << "Unexpected number of surface array names" << std::endl;
delete part;
return EXIT_FAILURE;
}
if (odeWavelet->GetSurfaceArrayComps()->GetNumberOfValues() != 1)
{
std::cerr << "Unexpected number of surface array comps" << std::endl;
delete part;
return EXIT_FAILURE;
}
if (odeWavelet->GetSurfaceArrayEnumValues()->GetNumberOfValues() != 11)
{
std::cerr << "Unexpected number of surface array enum values" << std::endl;
delete part;
return EXIT_FAILURE;
}
if (odeWavelet->GetSurfaceArrayTypes()->GetNumberOfValues() != 1)
{
std::cerr << "Unexpected number of surface array types" << std::endl;
delete part;
return EXIT_FAILURE;
}
......@@ -282,7 +265,6 @@ int TestLagrangianIntegrationModel(int, char*[])
vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDiameter");
odeTriangle->SetInputArrayToProcess(7, 1, 0,
vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDensity");
odeTriangle->SetCurrentParticle(part);
odeTriangle->SetLocator(locator);
odeTriangle->AddDataSet(triangle->GetOutput());
......@@ -301,7 +283,6 @@ int TestLagrangianIntegrationModel(int, char*[])
vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDiameter");
odeTransform->SetInputArrayToProcess(7, 1, 0,
vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDensity");
odeTransform->SetCurrentParticle(part);
odeTransform->SetLocator(locator);
odeTransform->AddDataSet(transform->GetOutput());
odeTransform->AddDataSet(wavelet->GetOutput());
......@@ -320,33 +301,29 @@ int TestLagrangianIntegrationModel(int, char*[])
{
x[0] = x0;
y[0] = x0 + xTranslation;
bool imageTest = (odeWavelet->FunctionValues(x, f) == 1);
bool locatorsTest = odeWavelet->FindInLocators(x);
bool unstrucTest = (odeTriangle->FunctionValues(x, f) == 1);
bool multiTest = (odeTransform->FunctionValues(y, f) == 1);
bool imageTest = (odeWavelet->FunctionValues(x, f, &part) == 1);
bool locatorsTest = odeWavelet->FindInLocators(x, &part);
bool unstrucTest = (odeTriangle->FunctionValues(x, f, &part) == 1);