diff --git a/Filters/FlowPaths/Testing/Cxx/TestLagrangianIntegrationModel.cxx b/Filters/FlowPaths/Testing/Cxx/TestLagrangianIntegrationModel.cxx index 181439225399d8717346fb89336ad5487a4fb698..0657d68ef333dc8c85fc2294364bd303489f8295 100644 --- a/Filters/FlowPaths/Testing/Cxx/TestLagrangianIntegrationModel.cxx +++ b/Filters/FlowPaths/Testing/Cxx/TestLagrangianIntegrationModel.cxx @@ -172,10 +172,10 @@ int TestLagrangianIntegrationModel(int, char*[]) } odeWavelet->SetUseInitialIntegrationTime(false); - vtkLagrangianParticle part(nvar, seedIdx, seedIdx, 0, 0, pd, odeWavelet->GetWeightsSize()); + vtkLagrangianParticle part(nvar, seedIdx, seedIdx, 0, 0, pd, odeWavelet->GetWeightsSize(), 3); - odeWavelet->InitializeVariablesParticleData(pd); - odeWavelet->InsertVariablesParticleData(&part, pd, 0); + odeWavelet->InitializeParticleData(pd); + odeWavelet->InsertParticleData(&part, pd, 0); odeWavelet->InitializeParticle(&part); if (odeWavelet->CheckFreeFlightTermination(&part)) { @@ -245,7 +245,7 @@ int TestLagrangianIntegrationModel(int, char*[]) double* tmpPt = nullptr; double tmpVald = 0; int tmpVali = 0; - if (odeWavelet->ManualIntegration(tmpPt, tmpPt, 0, tmpVald, tmpVald, 0, 0, 0, tmpVald, tmpVali)) + if (odeWavelet->ManualIntegration(tmpPt, tmpPt, 0, tmpVald, tmpVald, 0, 0, 0, 1, tmpVald, tmpVali)) { std::cerr << "ManualIntegration should do nothing and return false with matida model" << std::endl; } diff --git a/Filters/FlowPaths/Testing/Cxx/TestLagrangianParticle.cxx b/Filters/FlowPaths/Testing/Cxx/TestLagrangianParticle.cxx index 9c8f0e497e2ee02fa08c9dccbcaef60ac1988fdc..fd10c6cfa7d3bf84681f5720bbe96c4825682793 100644 --- a/Filters/FlowPaths/Testing/Cxx/TestLagrangianParticle.cxx +++ b/Filters/FlowPaths/Testing/Cxx/TestLagrangianParticle.cxx @@ -37,7 +37,7 @@ int TestLagrangianParticle(int, char*[]) vtkIdType particleCounter = 0; vtkLagrangianParticle* part = new vtkLagrangianParticle(nvar, seedId, - particleCounter, seedId, 0, pd, 8); + particleCounter, seedId, 0, pd, 8, 3); particleCounter++; if (nvar != part->GetNumberOfVariables()) { @@ -380,10 +380,10 @@ int TestLagrangianParticle(int, char*[]) particleCounter = 0; vtkLagrangianParticle* part4 = new vtkLagrangianParticle(nvar, seedId, - particleCounter, seedId, 0, pd, 8); + particleCounter, seedId, 0, pd, 8, 17); particleCounter++; vtkLagrangianParticle* part5 = vtkLagrangianParticle::NewInstance(nvar, seedId, - particleCounter, seedId, 0.17, pd, 8, 17, 0.13); + particleCounter, seedId, 0.17, pd, 8, 7, 17, 0.13); particleCounter++; if (part4->GetId() != 0) { diff --git a/Filters/FlowPaths/Testing/Cxx/TestLagrangianParticleTracker.cxx b/Filters/FlowPaths/Testing/Cxx/TestLagrangianParticleTracker.cxx index c6eb3daa87c4845990e149f33332ac2a59a78bb5..ffe7ce58fc01311bbaf6e81d75c680db5d5ef41c 100644 --- a/Filters/FlowPaths/Testing/Cxx/TestLagrangianParticleTracker.cxx +++ b/Filters/FlowPaths/Testing/Cxx/TestLagrangianParticleTracker.cxx @@ -191,6 +191,7 @@ int TestLagrangianParticleTracker(int, char*[]) vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDiameter"); integrationModel->SetInputArrayToProcess(7, 1, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDensity"); + integrationModel->SetNumberOfTrackedUserData(13); // Put in tracker vtkNew<vtkLagrangianParticleTracker> tracker; @@ -226,9 +227,9 @@ int TestLagrangianParticleTracker(int, char*[]) tracker->SetCellLengthComputationMode( vtkLagrangianParticleTracker::STEP_CUR_CELL_VEL_DIR); tracker->AdaptiveStepReintegrationOn(); - tracker->UseParticlePathsRenderingThresholdOn(); - tracker->SetParticlePathsRenderingPointsThreshold(100); + tracker->GenerateParticlePathsOutputOff(); tracker->Update(); + tracker->GenerateParticlePathsOutputOn(); tracker->SetInputConnection(ugFlow->GetOutputPort()); tracker->SetMaximumNumberOfSteps(30); tracker->SetCellLengthComputationMode( @@ -245,7 +246,6 @@ int TestLagrangianParticleTracker(int, char*[]) tracker->SetCellLengthComputationMode( vtkLagrangianParticleTracker::STEP_LAST_CELL_VEL_DIR); tracker->AdaptiveStepReintegrationOff(); - tracker->UseParticlePathsRenderingThresholdOff(); tracker->Update(); if (tracker->GetStepFactor() != 0.1) { @@ -282,14 +282,9 @@ int TestLagrangianParticleTracker(int, char*[]) std::cerr << "Incorrect AdaptiveStepReintegration" << std::endl; return EXIT_FAILURE; } - if (tracker->GetUseParticlePathsRenderingThreshold()) + if (!tracker->GetGenerateParticlePathsOutput()) { - std::cerr << "Incorrect UseParticlePathsRenderingThreshold" << std::endl; - return EXIT_FAILURE; - } - if (tracker->GetParticlePathsRenderingPointsThreshold() != 100) - { - std::cerr << "Incorrect ParticlePathsRenderingThreshold" << std::endl; + std::cerr << "Incorrect GenerateParticlePathsOutput" << std::endl; return EXIT_FAILURE; } tracker->Print(cout); diff --git a/Filters/FlowPaths/vtkLagrangianBasicIntegrationModel.cxx b/Filters/FlowPaths/vtkLagrangianBasicIntegrationModel.cxx index 4d061454db9d8e0e22e509815032a3ba8ad29b76..ce870872f280fbaaa069cd99661e228fbecb5019 100644 --- a/Filters/FlowPaths/vtkLagrangianBasicIntegrationModel.cxx +++ b/Filters/FlowPaths/vtkLagrangianBasicIntegrationModel.cxx @@ -25,6 +25,7 @@ #include "vtkIntArray.h" #include "vtkLagrangianParticle.h" #include "vtkLagrangianParticleTracker.h" +#include "vtkLongLongArray.h" #include "vtkMath.h" #include "vtkNew.h" #include "vtkPointData.h" @@ -1584,7 +1585,8 @@ vtkIntArray* vtkLagrangianBasicIntegrationModel::GetSurfaceArrayTypes() bool vtkLagrangianBasicIntegrationModel::ManualIntegration(double* vtkNotUsed(xcur), double* vtkNotUsed(xnext), double vtkNotUsed(t), double& vtkNotUsed(delT), double& vtkNotUsed(delTActual), double vtkNotUsed(minStep), double vtkNotUsed(maxStep), - double vtkNotUsed(maxError), double& vtkNotUsed(error), int& vtkNotUsed(integrationResult)) + double vtkNotUsed(maxError), double vtkNotUsed(cellLength), + double& vtkNotUsed(error), int& vtkNotUsed(integrationResult)) { return false; } @@ -1597,3 +1599,133 @@ void vtkLagrangianBasicIntegrationModel::ComputeSurfaceDefaultValues( (strcmp(arrayName, "SurfaceType") == 0) ? static_cast<double>(SURFACE_TYPE_TERM) : 0.0; std::fill(defaultValues, defaultValues + nComponents, defVal); } + +//--------------------------------------------------------------------------- +void vtkLagrangianBasicIntegrationModel::InitializeParticleData(vtkFieldData* particleData, int maxTuple) +{ + vtkNew<vtkIntArray> particleStepNumArray; + particleStepNumArray->SetName("StepNumber"); + particleStepNumArray->SetNumberOfComponents(1); + particleStepNumArray->Allocate(maxTuple); + particleData->AddArray(particleStepNumArray); + + vtkNew<vtkDoubleArray> particleVelArray; + particleVelArray->SetName("ParticleVelocity"); + particleVelArray->SetNumberOfComponents(3); + particleVelArray->Allocate(maxTuple * 3); + particleData->AddArray(particleVelArray); + + vtkNew<vtkDoubleArray> particleIntegrationTimeArray; + particleIntegrationTimeArray->SetName("IntegrationTime"); + particleIntegrationTimeArray->SetNumberOfComponents(1); + particleIntegrationTimeArray->Allocate(maxTuple); + particleData->AddArray(particleIntegrationTimeArray); +} + +//--------------------------------------------------------------------------- +void vtkLagrangianBasicIntegrationModel::InitializePathData(vtkFieldData* data) +{ + vtkNew<vtkLongLongArray> particleIdArray; + particleIdArray->SetName("Id"); + particleIdArray->SetNumberOfComponents(1); + data->AddArray(particleIdArray); + + vtkNew<vtkLongLongArray> particleParentIdArray; + particleParentIdArray->SetName("ParentId"); + particleParentIdArray->SetNumberOfComponents(1); + data->AddArray(particleParentIdArray); + + vtkNew<vtkLongLongArray> particleSeedIdArray; + particleSeedIdArray->SetName("SeedId"); + particleSeedIdArray->SetNumberOfComponents(1); + data->AddArray(particleSeedIdArray); + + vtkNew<vtkIntArray> particleTerminationArray; + particleTerminationArray->SetName("Termination"); + particleTerminationArray->SetNumberOfComponents(1); + data->AddArray(particleTerminationArray); +} + +//--------------------------------------------------------------------------- +void vtkLagrangianBasicIntegrationModel::InitializeInteractionData(vtkFieldData* data) +{ + vtkNew<vtkIntArray> interactionArray; + interactionArray->SetName("Interaction"); + interactionArray->SetNumberOfComponents(1); + data->AddArray(interactionArray); +} + +//--------------------------------------------------------------------------- +void vtkLagrangianBasicIntegrationModel::InsertSeedData( + vtkLagrangianParticle* particle, vtkFieldData* data) +{ + // Check for max number of tuples in arrays + vtkIdType maxTuples = 0; + for (int i = 0; i < data->GetNumberOfArrays(); i++) + { + maxTuples = std::max(data->GetArray(i)->GetNumberOfTuples(), maxTuples); + } + + // Copy seed data in not yet written array only + // ie not yet at maxTuple + vtkPointData* seedData = particle->GetSeedData(); + for (int i = 0; i < seedData->GetNumberOfArrays(); i++) + { + const char* name = seedData->GetArrayName(i); + vtkDataArray* arr = data->GetArray(name); + if (arr->GetNumberOfTuples() < maxTuples) + { + arr->InsertNextTuple(seedData->GetArray(i)->GetTuple(particle->GetSeedArrayTupleIndex())); + } + } +} + +//--------------------------------------------------------------------------- +void vtkLagrangianBasicIntegrationModel::InsertPathData( + vtkLagrangianParticle* particle, vtkFieldData* data) +{ + vtkLongLongArray::SafeDownCast(data->GetArray("Id"))->InsertNextValue(particle->GetId()); + vtkLongLongArray::SafeDownCast(data->GetArray("ParentId")) + ->InsertNextValue(particle->GetParentId()); + vtkLongLongArray::SafeDownCast(data->GetArray("SeedId"))->InsertNextValue(particle->GetSeedId()); + vtkIntArray::SafeDownCast(data->GetArray("Termination")) + ->InsertNextValue(particle->GetTermination()); +} + +//--------------------------------------------------------------------------- +void vtkLagrangianBasicIntegrationModel::InsertInteractionData( + vtkLagrangianParticle* particle, vtkFieldData* data) +{ + vtkIntArray::SafeDownCast(data->GetArray("Interaction")) + ->InsertNextValue(particle->GetInteraction()); +} + +//--------------------------------------------------------------------------- +void vtkLagrangianBasicIntegrationModel::InsertParticleData( + vtkLagrangianParticle* particle, vtkFieldData* data, int stepEnum) +{ + switch (stepEnum) + { + case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_PREV: + vtkIntArray::SafeDownCast(data->GetArray("StepNumber")) + ->InsertNextValue(particle->GetNumberOfSteps() - 1); + data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetPrevVelocity()); + data->GetArray("IntegrationTime")->InsertNextTuple1(particle->GetPrevIntegrationTime()); + break; + case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_CURRENT: + vtkIntArray::SafeDownCast(data->GetArray("StepNumber")) + ->InsertNextValue(particle->GetNumberOfSteps()); + data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetVelocity()); + data->GetArray("IntegrationTime")->InsertNextTuple1(particle->GetIntegrationTime()); + break; + case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_NEXT: + vtkIntArray::SafeDownCast(data->GetArray("StepNumber")) + ->InsertNextValue(particle->GetNumberOfSteps() + 1); + data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetNextVelocity()); + data->GetArray("IntegrationTime") + ->InsertNextTuple1(particle->GetIntegrationTime() + particle->GetStepTimeRef()); + break; + default: + break; + } +} diff --git a/Filters/FlowPaths/vtkLagrangianBasicIntegrationModel.h b/Filters/FlowPaths/vtkLagrangianBasicIntegrationModel.h index f820f163d46b99449eb2a8daf31335e1a8cd9954..bba953a6ca7f20cf627c2cbdf2fefb00026f1a04 100644 --- a/Filters/FlowPaths/vtkLagrangianBasicIntegrationModel.h +++ b/Filters/FlowPaths/vtkLagrangianBasicIntegrationModel.h @@ -40,9 +40,12 @@ * Inherited classes could reimplement IntersectWithLine to use a specific algorithm * to intersect particles and surface cells. * - * Inherited class could reimplement CheckFreeFlightTermination to set + * Inherited classes could reimplement CheckFreeFlightTermination to set * the way particles terminate in free flight. * + * Inherited classes could reimplement Initialize*Data and Insert*Data in order + * to customize the output of the tracker + * * @sa * vtkLagrangianParticleTracker vtkLagrangianParticle * vtkLagrangianMatidaIntegrationModel @@ -216,39 +219,6 @@ public: virtual bool FindInLocators(double* x, vtkLagrangianParticle* particle = nullptr); //@} - /** - * Empty method to be reimplemented if necessary - * in inherited classes. - * Allows a intehrited class to create - * Specific array in the output point data - * for storing variables. - */ - virtual void InitializeVariablesParticleData( - vtkPointData* vtkNotUsed(particleData), int vtkNotUsed(maxTuples) = 0) {} - - /** - * Empty method to be reimplemented if necessary in inherited classes. - * Allows a inherited class to create specific array in the output point data - * for filling point data with variables. - */ - virtual void InsertVariablesParticleData(vtkLagrangianParticle* vtkNotUsed(particle), - vtkPointData* vtkNotUsed(particleData), int vtkNotUsed(stepEnum)) {} - - /** - * Empty method to be reimplemented if necessary in inherited classes. - * Allows an inherited class to create specific array in the outputs - * field data associated with each particle path. - */ - virtual void InitializeModelPathData(vtkFieldData* vtkNotUsed(data)) {} - - /** - * Empty method to be reimplemented if necessary in inherited classes. - * Allows a inherited class to insert data in initialized - * array in the outputs field data associated with each particle path. - */ - virtual void InsertModelPathData( - vtkLagrangianParticle* vtkNotUsed(particle), vtkFieldData* vtkNotUsed(data)) {} - /** * Initialize a particle by setting user variables and perform any user * model specific operation. empty in basic implementation. @@ -370,8 +340,8 @@ public: * This method is thread-safe, its reimplementation should still be thread-safe. */ virtual bool ManualIntegration(double* xcur, double* xnext, double t, double& delT, - double& delTActual, double minStep, double maxStep, double maxError, double& error, - int& integrationResult); + double& delTActual, double minStep, double maxStep, double maxError, double cellLength, + double& error, int& integrationResult); /** * Method called by parallel algorithm @@ -402,6 +372,66 @@ public: */ virtual vtkAbstractArray* GetSeedArray(int idx, vtkPointData* pointData); + /** + * Set/Get the number of tracked user data attached to the particles. + * Tracked user data are data that are related to each particle position + * but are not integrated like the user variables. + * They are not saved in the particle path. + * Default is 0. + */ + vtkSetMacro(NumberOfTrackedUserData, int); + vtkGetMacro(NumberOfTrackedUserData, int); + + /** + * Method used by the LPT to initialize data insertion in the provided + * vtkFieldData. It initializes Id, ParentID, SeedID and Termination. + * Reimplement as needed in acccordance with InsertPathData. + */ + virtual void InitializePathData(vtkFieldData* data); + + /** + * Method used by the LPT to initialize data insertion in the provided + * vtkFieldData. It initializes Interaction. + * Reimplement as needed in acccordance with InsertInteractionData. + */ + virtual void InitializeInteractionData(vtkFieldData* data); + + /** + * Method used by the LPT to initialize data insertion in the provided + * vtkFieldData. It initializes StepNumber, ParticleVelocity, IntegrationTime. + * Reimplement as needed in acccordance with InsertParticleData. + */ + virtual void InitializeParticleData(vtkFieldData* particleData, int maxTuples = 0); + + /** + * Method used by the LPT to insert data from the partice into + * the provided vtkFieldData. It inserts Id, ParentID, SeedID and Termination. + * Reimplement as needed in acccordance with InitializePathData. + */ + virtual void InsertPathData(vtkLagrangianParticle* particle, vtkFieldData* data); + + /** + * Method used by the LPT to insert data from the partice into + * the provided vtkFieldData. It inserts Interaction. + * Reimplement as needed in acccordance with InitializeInteractionData. + */ + virtual void InsertInteractionData(vtkLagrangianParticle* particle, vtkFieldData* data); + + /** + * Method used by the LPT to insert data from the partice into + * the provided vtkFieldData. It inserts StepNumber, ParticleVelocity, IntegrationTime. + * stepEnum enables to select which data to insert, Prev, Current or Next. + * Reimplement as needed in acccordance with InitializeParticleData. + */ + virtual void InsertParticleData(vtkLagrangianParticle* particle, vtkFieldData* data, int stepEnum); + + /** + * Method used by the LPT to insert data from the partice into + * the provided vtkFieldData. It insert alls array available in the SeedData. + * Reimplement as needed. + */ + virtual void InsertSeedData(vtkLagrangianParticle* particle, vtkFieldData* data); + protected: vtkLagrangianBasicIntegrationModel(); ~vtkLagrangianBasicIntegrationModel() override; @@ -520,7 +550,7 @@ protected: virtual int GetFlowOrSurfaceDataFieldAssociation(int idx); /** - * Methods used by ParaView surface helper to get default + * Method used by ParaView surface helper to get default * values for each leaf of each dataset of surface * nComponents could be retrieved with arrayName but is * given for simplication purposes. @@ -557,6 +587,7 @@ protected: double Tolerance; bool NonPlanarQuadSupport; bool UseInitialIntegrationTime; + int NumberOfTrackedUserData = 0; vtkNew<vtkStringArray> SeedArrayNames; vtkNew<vtkIntArray> SeedArrayComps; diff --git a/Filters/FlowPaths/vtkLagrangianParticle.cxx b/Filters/FlowPaths/vtkLagrangianParticle.cxx index ee18b00589736fbdd1bade6f4d5aac419f82bb72..53c6946140c3827ac154d6313dc07d32fabcd10c 100644 --- a/Filters/FlowPaths/vtkLagrangianParticle.cxx +++ b/Filters/FlowPaths/vtkLagrangianParticle.cxx @@ -26,7 +26,7 @@ std::mutex seedDataMutex; //--------------------------------------------------------------------------- vtkLagrangianParticle::vtkLagrangianParticle(int numberOfVariables, vtkIdType seedId, vtkIdType particleId, vtkIdType seedArrayTupleIndex, double integrationTime, - vtkPointData* seedData, int weightsSize) + vtkPointData* seedData, int weightsSize, int numberOfTrackedUserData) : Id(particleId) , ParentId(-1) , SeedId(seedId) @@ -66,24 +66,30 @@ vtkLagrangianParticle::vtkLagrangianParticle(int numberOfVariables, vtkIdType se // Initialize surface cell cache this->LastSurfaceCellId = -1; this->LastSurfaceDataSet = nullptr; + + // Initialize tracked user data + this->PrevTrackedUserData.resize(numberOfTrackedUserData, 0); + this->TrackedUserData.resize(numberOfTrackedUserData, 0); + this->NextTrackedUserData.resize(numberOfTrackedUserData, 0); } //--------------------------------------------------------------------------- vtkLagrangianParticle* vtkLagrangianParticle::NewInstance(int numberOfVariables, vtkIdType seedId, vtkIdType particleId, vtkIdType seedArrayTupleIndex, double integrationTime, - vtkPointData* seedData, int weightsSize) + vtkPointData* seedData, int weightsSize, int numberOfTrackedUserData) { return new vtkLagrangianParticle(numberOfVariables, seedId, particleId, seedArrayTupleIndex, - integrationTime, seedData, weightsSize); + integrationTime, seedData, weightsSize, numberOfTrackedUserData); } //--------------------------------------------------------------------------- vtkLagrangianParticle* vtkLagrangianParticle::NewInstance(int numberOfVariables, vtkIdType seedId, vtkIdType particleId, vtkIdType seedArrayTupleIndex, double integrationTime, - vtkPointData* seedData, int weightsSize, vtkIdType numberOfSteps, double previousIntegrationTime) + vtkPointData* seedData, int weightsSize, int numberOfTrackedUserData, + vtkIdType numberOfSteps, double previousIntegrationTime) { vtkLagrangianParticle* particle = new vtkLagrangianParticle(numberOfVariables, seedId, particleId, - seedArrayTupleIndex, integrationTime, seedData, weightsSize); + seedArrayTupleIndex, integrationTime, seedData, weightsSize, numberOfTrackedUserData); particle->NumberOfSteps = numberOfSteps; particle->PrevIntegrationTime = previousIntegrationTime; return particle; @@ -108,7 +114,8 @@ vtkLagrangianParticle* vtkLagrangianParticle::NewParticle(vtkIdType particleId) // Create particle and copy members vtkLagrangianParticle* particle = this->NewInstance(this->GetNumberOfVariables(), this->GetSeedId(), particleId, - seedArrayTupleIndex, this->IntegrationTime + this->StepTime, seedData, this->WeightsSize); + seedArrayTupleIndex, this->IntegrationTime + this->StepTime, seedData, this->WeightsSize, + static_cast<int>(this->TrackedUserData.size())); particle->ParentId = this->GetId(); particle->NumberOfSteps = this->GetNumberOfSteps() + 1; @@ -116,6 +123,12 @@ vtkLagrangianParticle* vtkLagrangianParticle::NewParticle(vtkIdType particleId) std::copy(this->EquationVariables.begin(), this->EquationVariables.end(), particle->PrevEquationVariables.begin()); std::copy(this->NextEquationVariables.begin(), this->NextEquationVariables.end(), particle->EquationVariables.begin()); std::fill(particle->NextEquationVariables.begin(), particle->NextEquationVariables.end(), 0); + + // Copy UserData + std::copy(this->TrackedUserData.begin(), this->TrackedUserData.end(), particle->PrevTrackedUserData.begin()); + std::copy(this->NextTrackedUserData.begin(), this->NextTrackedUserData.end(), particle->TrackedUserData.begin()); + std::fill(particle->NextTrackedUserData.begin(), particle->NextTrackedUserData.end(), 0); + return particle; } @@ -124,7 +137,7 @@ vtkLagrangianParticle* vtkLagrangianParticle::CloneParticle() { vtkLagrangianParticle* clone = this->NewInstance(this->GetNumberOfVariables(), this->GetSeedId(), this->GetId(), this->GetSeedArrayTupleIndex(), this->IntegrationTime, this->GetSeedData(), - this->WeightsSize); + this->WeightsSize, static_cast<int>(this->TrackedUserData.size())); clone->Id = this->Id; clone->ParentId = this->ParentId; clone->NumberOfSteps = this->NumberOfSteps; @@ -132,6 +145,9 @@ vtkLagrangianParticle* vtkLagrangianParticle::CloneParticle() std::copy(this->PrevEquationVariables.begin(), this->PrevEquationVariables.end(), clone->PrevEquationVariables.begin()); std::copy(this->EquationVariables.begin(), this->EquationVariables.end(), clone->EquationVariables.begin()); std::copy(this->NextEquationVariables.begin(), this->NextEquationVariables.end(), clone->NextEquationVariables.begin()); + std::copy(this->PrevTrackedUserData.begin(), this->PrevTrackedUserData.end(), clone->PrevTrackedUserData.begin()); + std::copy(this->TrackedUserData.begin(), this->TrackedUserData.end(), clone->TrackedUserData.begin()); + std::copy(this->NextTrackedUserData.begin(), this->NextTrackedUserData.end(), clone->NextTrackedUserData.begin()); clone->StepTime = this->StepTime; return clone; } @@ -342,6 +358,9 @@ void vtkLagrangianParticle::MoveToNextPosition() std::copy(this->EquationVariables.begin(), this->EquationVariables.end(), this->PrevEquationVariables.begin()); std::copy(this->NextEquationVariables.begin(), this->NextEquationVariables.end(), this->EquationVariables.begin()); std::fill(this->NextEquationVariables.begin(), this->NextEquationVariables.end(), 0); + std::copy(this->TrackedUserData.begin(), this->TrackedUserData.end(), this->PrevTrackedUserData.begin()); + std::copy(this->NextTrackedUserData.begin(), this->NextTrackedUserData.end(), this->TrackedUserData.begin()); + std::fill(this->NextTrackedUserData.begin(), this->NextTrackedUserData.end(), 0); this->NumberOfSteps++; this->PrevIntegrationTime = this->IntegrationTime; @@ -368,23 +387,44 @@ void vtkLagrangianParticle::PrintSelf(ostream& os, vtkIndent indent) os << indent << "Interaction: " << this->Interaction << std::endl; os << indent << "PrevEquationVariables:"; - for (int i = 0; i < this->NumberOfVariables; i++) + for (auto var : this->PrevEquationVariables) { - os << indent << " " << this->PrevEquationVariables[i]; + os << indent << " " << var; } os << std::endl; os << indent << "EquationVariables:"; - for (int i = 0; i < this->NumberOfVariables; i++) + for (auto var : this->EquationVariables) { - os << indent << " " << this->EquationVariables[i]; + os << indent << " " << var; } os << std::endl; os << indent << "NextEquationVariables:"; - for (int i = 0; i < this->NumberOfVariables; i++) + for (auto var : this->NextEquationVariables) + { + os << indent << " " << var; + } + os << std::endl; + + os << indent << "PrevTrackedUserData:"; + for (auto var : this->PrevTrackedUserData) + { + os << indent << " " << var; + } + os << std::endl; + + os << indent << "TrackedUserData:"; + for (auto var : this->TrackedUserData) + { + os << indent << " " << var; + } + os << std::endl; + + os << indent << "NextTrackedUserData:"; + for (auto var : this->NextTrackedUserData) { - os << indent << " " << this->NextEquationVariables[i]; + os << indent << " " << var; } os << std::endl; } diff --git a/Filters/FlowPaths/vtkLagrangianParticle.h b/Filters/FlowPaths/vtkLagrangianParticle.h index 99031b7c157f7a8263e9cf5cb454b6e326a8710d..3c3046420b101c4df8dffb673748801372c41064 100644 --- a/Filters/FlowPaths/vtkLagrangianParticle.h +++ b/Filters/FlowPaths/vtkLagrangianParticle.h @@ -99,7 +99,8 @@ public: * particle data is a pointer to the pointData associated to all particles. */ vtkLagrangianParticle(int numberOfVariables, vtkIdType seedId, vtkIdType particleId, - vtkIdType seedArrayTupleIndex, double integrationTime, vtkPointData* seedData, int weightsSize); + vtkIdType seedArrayTupleIndex, double integrationTime, vtkPointData* seedData, + int weightsSize, int numberOfTrackedUserData); /** * Constructor wrapper to create a partially integrated particle in the domain. @@ -107,7 +108,7 @@ public: */ static vtkLagrangianParticle* NewInstance(int numberOfVariables, vtkIdType seedId, vtkIdType particleId, vtkIdType seedArrayTupleIndex, double integrationTime, - vtkPointData* seedData, int weightsSize, vtkIdType numberOfSteps, + vtkPointData* seedData, int weightsSize, int numberOfTrackedUserData, vtkIdType numberOfSteps, double previousIntegrationTime); /** @@ -242,6 +243,47 @@ public: inline double* GetNextUserVariables() { return this->NextUserVariables; } //@} + //@{ + /** + * Get a reference to PrevTrackedUserData + * See GetTrackedUserData for an explanation on how to use it. + */ + inline std::vector<double>& GetPrevTrackedUserData() { return this->PrevTrackedUserData; } + //@} + + //@{ + /** + * Get a reference to TrackedUserData. + * The tracked user data is a vector of double associated with each position of the particle, + * but it is not integrated contrary to the UserVariables and EquationVariables. + * It is, however, automatically tracked from one position to the next, copied when creating + * new particles with NewInstance and CloneParticle and transferred from one node to the next + * when particles move from one domain to the another in parallel. + * If you are using these, you are supposed to compute and set the next tracked user data + * your implementation of FunctionValues in your model. + */ + inline std::vector<double>& GetTrackedUserData() { return this->TrackedUserData; } + //@} + + //@{ + /** + * Get a reference to NextTrackedUserData + * See GetTrackedUserData for an explanation on how to use it. + */ + inline std::vector<double>& GetNextTrackedUserData() { return this->NextTrackedUserData; } + //@} + + //@{ + /** + * Get/Set a pointer to TemporaryUserData + * This data is not tracked and not transferred nor copied + * This can be used to store any kind of data, structure, class instance that you may need. + * Be cautious if the pointer is shared as particle integration can be multithreaded. + */ + inline void* GetTemporaryUserData() { return this->TemporaryUserData; } + inline void SetTemporaryUserData(void * tempUserData) { this->TemporaryUserData = tempUserData; } + //@} + /** * Move the particle to its next position by putting next equation * variable to equation variable and clearing next equation variable. @@ -426,7 +468,8 @@ protected: * Constructor wrapper for internal convenience */ vtkLagrangianParticle* NewInstance(int numberOfVariables, vtkIdType seedId, vtkIdType particleId, - vtkIdType seedArrayTupleIndex, double integrationTime, vtkPointData* seedData, int weightsSize); + vtkIdType seedArrayTupleIndex, double integrationTime, vtkPointData* seedData, int weightsSize, + int numberOfTrackedUserData); vtkLagrangianParticle(const vtkLagrangianParticle&) = delete; vtkLagrangianParticle() = delete; @@ -444,6 +487,12 @@ protected: double* NextVelocity; double* NextUserVariables; + std::vector<double> PrevTrackedUserData; + std::vector<double> TrackedUserData; + std::vector<double> NextTrackedUserData; + + void* TemporaryUserData = nullptr; + vtkIdType Id; vtkIdType ParentId; vtkIdType SeedId; diff --git a/Filters/FlowPaths/vtkLagrangianParticleTracker.cxx b/Filters/FlowPaths/vtkLagrangianParticleTracker.cxx index e9b7a041633b5c5ede4181d6236809b50659a41c..14b5855edb0ffbff5a28c01059e666aa2d7d721a 100644 --- a/Filters/FlowPaths/vtkLagrangianParticleTracker.cxx +++ b/Filters/FlowPaths/vtkLagrangianParticleTracker.cxx @@ -28,7 +28,6 @@ #include "vtkInformationVector.h" #include "vtkLagrangianMatidaIntegrationModel.h" #include "vtkLagrangianParticle.h" -#include "vtkLongLongArray.h" #include "vtkNew.h" #include "vtkObjectFactory.h" #include "vtkPointData.h" @@ -121,9 +120,7 @@ vtkLagrangianParticleTracker::vtkLagrangianParticleTracker() , MaximumNumberOfSteps(100) , MaximumIntegrationTime(-1.0) , AdaptiveStepReintegration(false) - , UseParticlePathsRenderingThreshold(false) , GeneratePolyVertexInteractionOutput(false) - , ParticlePathsRenderingPointsThreshold(100) , ParticleCounter(0) , IntegratedParticleCounter(0) , IntegratedParticleCounterIncrement(1) @@ -174,11 +171,7 @@ void vtkLagrangianParticleTracker::PrintSelf(ostream& os, vtkIndent indent) os << indent << "MaximumNumberOfSteps: " << this->MaximumNumberOfSteps << endl; os << indent << "MaximumIntegrationTime: " << this->MaximumIntegrationTime << endl; os << indent << "AdaptiveStepReintegration: " << this->AdaptiveStepReintegration << endl; - os << indent << "UseParticlePathsRenderingThreshold: " << this->UseParticlePathsRenderingThreshold - << endl; - os << indent - << "ParticlePathsRenderingPointsThreshold: " << this->ParticlePathsRenderingPointsThreshold - << endl; + os << indent << "GenerateParticlePathsOutput: " << this->GenerateParticlePathsOutput << endl; os << indent << "MinimumVelocityMagnitude: " << this->MinimumVelocityMagnitude << endl; os << indent << "MinimumReductionFactor: " << this->MinimumReductionFactor << endl; os << indent << "ParticleCounter: " << this->ParticleCounter << endl; @@ -308,7 +301,7 @@ int vtkLagrangianParticleTracker::RequestData(vtkInformation* vtkNotUsed(request } // Initialize outputs - vtkPolyData* particlePathsOutput; + vtkPolyData* particlePathsOutput = nullptr; vtkDataObject* interactionOutput; if (!this->InitializeOutputs(outputVector, seedData, static_cast<vtkIdType>(particlesQueue.size()), surfaces, particlePathsOutput, @@ -369,14 +362,6 @@ int vtkLagrangianParticleTracker::RequestData(vtkInformation* vtkNotUsed(request return 1; } -//--------------------------------------------------------------------------- -bool vtkLagrangianParticleTracker::CheckParticlePathsRenderingThreshold( - vtkPolyData* particlePathsOutput) -{ - return this->UseParticlePathsRenderingThreshold && - particlePathsOutput->GetNumberOfPoints() > this->ParticlePathsRenderingPointsThreshold; -} - //--------------------------------------------------------------------------- vtkMTimeType vtkLagrangianParticleTracker::GetMTime() { @@ -458,43 +443,44 @@ bool vtkLagrangianParticleTracker::InitializePathsOutput(vtkInformationVector* o { // Prepare path output vtkInformation* particleOutInfo = outputVector->GetInformationObject(0); - particlePathsOutput = vtkPolyData::SafeDownCast(particleOutInfo->Get(vtkPolyData::DATA_OBJECT())); - if (!particlePathsOutput) - { - vtkErrorMacro(<< "Cannot find a vtkPolyData particle paths output. aborting"); - return false; - } - // Set information keys - particlePathsOutput->GetInformation()->Set( - vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(), - particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES())); - particlePathsOutput->GetInformation()->Set( - vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(), - particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER())); - particlePathsOutput->GetInformation()->Set( - vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), - particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS())); - - vtkNew<vtkPoints> particlePathsPoints; - vtkNew<vtkCellArray> particlePaths; - vtkNew<vtkCellArray> particleVerts; - particlePathsOutput->SetPoints(particlePathsPoints); - particlePathsOutput->SetLines(particlePaths); - particlePathsOutput->SetVerts(particleVerts); - - // Prepare particle paths output point data - vtkCellData* particlePathsCellData = particlePathsOutput->GetCellData(); - particlePathsCellData->CopyStructure(seedData); - this->InitializePathData(particlePathsCellData); - this->IntegrationModel->InitializeModelPathData(particlePathsCellData); + if (this->GenerateParticlePathsOutput) + { + particlePathsOutput = vtkPolyData::SafeDownCast(particleOutInfo->Get(vtkPolyData::DATA_OBJECT())); + if (!particlePathsOutput) + { + vtkErrorMacro(<< "Cannot find a vtkPolyData particle paths output. aborting"); + return false; + } - // Initialize Particle Paths Point Data - vtkPointData* particlePathsPointData = particlePathsOutput->GetPointData(); - this->InitializeParticleData(particlePathsPointData, numberOfSeeds); + // Set information keys + particlePathsOutput->GetInformation()->Set( + vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(), + particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES())); + particlePathsOutput->GetInformation()->Set( + vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(), + particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER())); + particlePathsOutput->GetInformation()->Set( + vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), + particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS())); + + vtkNew<vtkPoints> particlePathsPoints; + vtkNew<vtkCellArray> particlePaths; + vtkNew<vtkCellArray> particleVerts; + particlePathsOutput->SetPoints(particlePathsPoints); + particlePathsOutput->SetLines(particlePaths); + particlePathsOutput->SetVerts(particleVerts); + + // Prepare particle paths output point data + vtkCellData* particlePathsCellData = particlePathsOutput->GetCellData(); + particlePathsCellData->CopyStructure(seedData); + this->IntegrationModel->InitializePathData(particlePathsCellData); + + // Initialize Particle Paths Point Data + vtkPointData* particlePathsPointData = particlePathsOutput->GetPointData(); + this->IntegrationModel->InitializeParticleData(particlePathsPointData, numberOfSeeds); + } - // Initialize particle data from integration model, if any - this->IntegrationModel->InitializeVariablesParticleData(particlePathsPointData, numberOfSeeds); return true; } @@ -544,11 +530,9 @@ bool vtkLagrangianParticleTracker::InitializeInteractionOutput(vtkInformationVec vtkNew<vtkPoints> points; pd->SetPoints(points); pd->GetPointData()->CopyStructure(seedData); - this->InitializePathData(pd->GetPointData()); - this->IntegrationModel->InitializeModelPathData(pd->GetPointData()); - this->InitializeInteractionData(pd->GetPointData()); - this->InitializeParticleData(pd->GetPointData()); - this->IntegrationModel->InitializeVariablesParticleData(pd->GetPointData()); + this->IntegrationModel->InitializePathData(pd->GetPointData()); + this->IntegrationModel->InitializeInteractionData(pd->GetPointData()); + this->IntegrationModel->InitializeParticleData(pd->GetPointData()); hdOutput->SetDataSet(iter, pd); } } @@ -565,84 +549,30 @@ bool vtkLagrangianParticleTracker::InitializeInteractionOutput(vtkInformationVec vtkNew<vtkCellArray> cells; pd->SetPoints(points); pd->GetPointData()->CopyStructure(seedData); - this->InitializePathData(pd->GetPointData()); - this->IntegrationModel->InitializeModelPathData(pd->GetPointData()); - this->InitializeInteractionData(pd->GetPointData()); - this->InitializeParticleData(pd->GetPointData()); - this->IntegrationModel->InitializeVariablesParticleData(pd->GetPointData()); + this->IntegrationModel->InitializePathData(pd->GetPointData()); + this->IntegrationModel->InitializeInteractionData(pd->GetPointData()); + this->IntegrationModel->InitializeParticleData(pd->GetPointData()); } return true; } -//--------------------------------------------------------------------------- -void vtkLagrangianParticleTracker::InitializeParticleData(vtkFieldData* particleData, int maxTuple) -{ - vtkNew<vtkIntArray> particleStepNumArray; - particleStepNumArray->SetName("StepNumber"); - particleStepNumArray->SetNumberOfComponents(1); - particleStepNumArray->Allocate(maxTuple); - particleData->AddArray(particleStepNumArray); - - vtkNew<vtkDoubleArray> particleVelArray; - particleVelArray->SetName("ParticleVelocity"); - particleVelArray->SetNumberOfComponents(3); - particleVelArray->Allocate(maxTuple * 3); - particleData->AddArray(particleVelArray); - - vtkNew<vtkDoubleArray> particleIntegrationTimeArray; - particleIntegrationTimeArray->SetName("IntegrationTime"); - particleIntegrationTimeArray->SetNumberOfComponents(1); - particleIntegrationTimeArray->Allocate(maxTuple); - particleData->AddArray(particleIntegrationTimeArray); -} - -//--------------------------------------------------------------------------- -void vtkLagrangianParticleTracker::InitializePathData(vtkFieldData* data) -{ - vtkNew<vtkLongLongArray> particleIdArray; - particleIdArray->SetName("Id"); - particleIdArray->SetNumberOfComponents(1); - data->AddArray(particleIdArray); - - vtkNew<vtkLongLongArray> particleParentIdArray; - particleParentIdArray->SetName("ParentId"); - particleParentIdArray->SetNumberOfComponents(1); - data->AddArray(particleParentIdArray); - - vtkNew<vtkLongLongArray> particleSeedIdArray; - particleSeedIdArray->SetName("SeedId"); - particleSeedIdArray->SetNumberOfComponents(1); - data->AddArray(particleSeedIdArray); - - vtkNew<vtkIntArray> particleTerminationArray; - particleTerminationArray->SetName("Termination"); - particleTerminationArray->SetNumberOfComponents(1); - data->AddArray(particleTerminationArray); -} - -//--------------------------------------------------------------------------- -void vtkLagrangianParticleTracker::InitializeInteractionData(vtkFieldData* data) -{ - vtkNew<vtkIntArray> interactionArray; - interactionArray->SetName("Interaction"); - interactionArray->SetNumberOfComponents(1); - data->AddArray(interactionArray); -} - //--------------------------------------------------------------------------- bool vtkLagrangianParticleTracker::FinalizeOutputs( vtkPolyData* particlePathsOutput, vtkDataObject* interactionOutput) { - // Recover structures - vtkPointData* particlePathsPointData = particlePathsOutput->GetPointData(); - vtkPoints* particlePathsPoints = particlePathsOutput->GetPoints(); - - // Squeeze and resize point data - for (int i = 0; i < particlePathsPointData->GetNumberOfArrays(); i++) + if (particlePathsOutput) { - vtkDataArray* array = particlePathsPointData->GetArray(i); - array->Resize(particlePathsPoints->GetNumberOfPoints()); - array->Squeeze(); + // Recover structures + vtkPointData* particlePathsPointData = particlePathsOutput->GetPointData(); + vtkPoints* particlePathsPoints = particlePathsOutput->GetPoints(); + + // Squeeze and resize point data + for (int i = 0; i < particlePathsPointData->GetNumberOfArrays(); i++) + { + vtkDataArray* array = particlePathsPointData->GetArray(i); + array->Resize(particlePathsPoints->GetNumberOfPoints()); + array->Squeeze(); + } } // Insert interaction poly-vertex cell @@ -684,12 +614,6 @@ bool vtkLagrangianParticleTracker::FinalizeOutputs( // Enable model post processing this->IntegrationModel->FinalizeOutputs(particlePathsOutput, interactionOutput); - - // Optional paths rendering - if (this->CheckParticlePathsRenderingThreshold(particlePathsOutput)) - { - particlePathsOutput->Initialize(); - } return true; } @@ -964,7 +888,7 @@ void vtkLagrangianParticleTracker::GenerateParticles(const vtkBoundingBox* vtkNo { // Create and set a dummy particle so FindInLocators can use caching. vtkLagrangianParticle dummyParticle( - 0, 0, 0, 0, 0, nullptr, this->IntegrationModel->GetWeightsSize()); + 0, 0, 0, 0, 0, nullptr, this->IntegrationModel->GetWeightsSize(), 0); this->ParticleCounter = 0; this->IntegratedParticleCounter = 0; @@ -976,7 +900,7 @@ void vtkLagrangianParticleTracker::GenerateParticles(const vtkBoundingBox* vtkNo initialIntegrationTimes ? initialIntegrationTimes->GetTuple1(i) : 0; vtkIdType particleId = this->GetNewParticleId(); vtkLagrangianParticle* particle = new vtkLagrangianParticle(nVar, particleId, particleId, i, - initialIntegrationTime, seedData, this->IntegrationModel->GetWeightsSize()); + initialIntegrationTime, seedData, this->IntegrationModel->GetWeightsSize(), this->IntegrationModel->GetNumberOfTrackedUserData()); memcpy(particle->GetPosition(), position, 3 * sizeof(double)); initialVelocities->GetTuple(i, particle->GetVelocity()); this->IntegrationModel->InitializeParticle(particle); @@ -1009,8 +933,6 @@ int vtkLagrangianParticleTracker::Integrate(vtkInitialValueProblemSolver* integr return -1; } - vtkIdList* particlePathPointId = particlePath->GetPointIds(); - // Integrate until MaximumNumberOfSteps or MaximumIntegrationTime is reached or special case stops int integrationRes = 0; double stepFactor = this->StepFactor; @@ -1033,7 +955,7 @@ int vtkLagrangianParticleTracker::Integrate(vtkInitialValueProblemSolver* integr // Integrate one step if (!this->ComputeNextStep(integrator, particle->GetEquationVariables(), particle->GetNextEquationVariables(), particle->GetIntegrationTime(), stepTime, - stepTimeActual, stepTimeMin, stepTimeMax, integrationRes, particle)) + stepTimeActual, stepTimeMin, stepTimeMax, cellLength, integrationRes, particle)) { vtkErrorMacro(<< "Integration Error"); break; @@ -1112,10 +1034,12 @@ int vtkLagrangianParticleTracker::Integrate(vtkInitialValueProblemSolver* integr // Particle has been correctly integrated and interacted, record it // Insert Current particle as an output point + + if (particlePathsOutput) { // Mutex Locked Area std::lock_guard<std::mutex> guard(this->ParticlePathsOutputMutex); - this->InsertPathOutputPoint(particle, particlePathsOutput, particlePathPointId); + this->InsertPathOutputPoint(particle, particlePathsOutput, particlePath->GetPointIds()); } // Particle has been terminated by surface @@ -1124,9 +1048,12 @@ int vtkLagrangianParticleTracker::Integrate(vtkInitialValueProblemSolver* integr // Insert last particle path point on surface particle->MoveToNextPosition(); - // Mutex Locked Area - std::lock_guard<std::mutex> guard(this->ParticlePathsOutputMutex); - this->InsertPathOutputPoint(particle, particlePathsOutput, particlePathPointId); + if (particlePathsOutput) + { + // Mutex Locked Area + std::lock_guard<std::mutex> guard(this->ParticlePathsOutputMutex); + this->InsertPathOutputPoint(particle, particlePathsOutput, particlePath->GetPointIds()); + } // stop integration break; @@ -1161,24 +1088,26 @@ int vtkLagrangianParticleTracker::Integrate(vtkInitialValueProblemSolver* integr } } - if (particlePath->GetPointIds()->GetNumberOfIds() == 1) + if (particlePathsOutput) { - particlePath->GetPointIds()->InsertNextId(particlePath->GetPointId(0)); - } + if (particlePath->GetPointIds()->GetNumberOfIds() == 1) + { + particlePath->GetPointIds()->InsertNextId(particlePath->GetPointId(0)); + } - // Duplicate single point particle paths, to avoid degenerated lines. - if (particlePath->GetPointIds()->GetNumberOfIds() > 0) - { - // Mutex Locked Area - std::lock_guard<std::mutex> guard(this->ParticlePathsOutputMutex); + // Duplicate single point particle paths, to avoid degenerated lines. + if (particlePath->GetPointIds()->GetNumberOfIds() > 0) + { + // Mutex Locked Area + std::lock_guard<std::mutex> guard(this->ParticlePathsOutputMutex); - // Add particle path or vertex to cell array - particlePathsOutput->GetLines()->InsertNextCell(particlePath); - this->InsertPathData(particle, particlePathsOutput->GetCellData()); - this->IntegrationModel->InsertModelPathData(particle, particlePathsOutput->GetCellData()); + // Add particle path or vertex to cell array + particlePathsOutput->GetLines()->InsertNextCell(particlePath); + this->IntegrationModel->InsertPathData(particle, particlePathsOutput->GetCellData()); - // Insert data from seed data only in not yet written arrays - this->InsertSeedData(particle, particlePathsOutput->GetCellData()); + // Insert data from seed data only in not yet written arrays + this->IntegrationModel->InsertSeedData(particle, particlePathsOutput->GetCellData()); + } } return integrationRes; @@ -1199,12 +1128,7 @@ void vtkLagrangianParticleTracker::InsertPathOutputPoint(vtkLagrangianParticle* particlePathPointId->InsertNextId(pointId); // Insert particle data - this->InsertParticleData(particle, particlePathsPointData, - prev ? vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_PREV - : vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_CURRENT); - - // Add Variables data - this->IntegrationModel->InsertVariablesParticleData(particle, particlePathsPointData, + this->IntegrationModel->InsertParticleData(particle, particlePathsPointData, prev ? vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_PREV : vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_CURRENT); } @@ -1248,94 +1172,13 @@ void vtkLagrangianParticleTracker::InsertInteractionOutputPoint(vtkLagrangianPar // Fill up interaction point data vtkPointData* pointData = interactionPd->GetPointData(); - this->InsertPathData(particle, pointData); - this->IntegrationModel->InsertModelPathData(particle, pointData); - this->InsertInteractionData(particle, pointData); - this->InsertParticleData( - particle, pointData, vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_NEXT); - - // Add Variables data - this->IntegrationModel->InsertVariablesParticleData( + this->IntegrationModel->InsertPathData(particle, pointData); + this->IntegrationModel->InsertInteractionData(particle, pointData); + this->IntegrationModel->InsertParticleData( particle, pointData, vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_NEXT); // Finally, Insert data from seed data only on not yet written arrays - this->InsertSeedData(particle, pointData); -} - -//--------------------------------------------------------------------------- -void vtkLagrangianParticleTracker::InsertSeedData( - vtkLagrangianParticle* particle, vtkFieldData* data) -{ - // Check for max number of tuples in arrays - vtkIdType maxTuples = 0; - for (int i = 0; i < data->GetNumberOfArrays(); i++) - { - maxTuples = std::max(data->GetArray(i)->GetNumberOfTuples(), maxTuples); - } - - // Copy seed data in not yet written array only - // ie not yet at maxTuple - vtkPointData* seedData = particle->GetSeedData(); - for (int i = 0; i < seedData->GetNumberOfArrays(); i++) - { - const char* name = seedData->GetArrayName(i); - vtkDataArray* arr = data->GetArray(name); - if (arr->GetNumberOfTuples() < maxTuples) - { - arr->InsertNextTuple(seedData->GetArray(i)->GetTuple(particle->GetSeedArrayTupleIndex())); - } - } - // here all arrays from data should have the exact same size -} - -//--------------------------------------------------------------------------- -void vtkLagrangianParticleTracker::InsertPathData( - vtkLagrangianParticle* particle, vtkFieldData* data) -{ - vtkLongLongArray::SafeDownCast(data->GetArray("Id"))->InsertNextValue(particle->GetId()); - vtkLongLongArray::SafeDownCast(data->GetArray("ParentId")) - ->InsertNextValue(particle->GetParentId()); - vtkLongLongArray::SafeDownCast(data->GetArray("SeedId"))->InsertNextValue(particle->GetSeedId()); - vtkIntArray::SafeDownCast(data->GetArray("Termination")) - ->InsertNextValue(particle->GetTermination()); -} - -//--------------------------------------------------------------------------- -void vtkLagrangianParticleTracker::InsertInteractionData( - vtkLagrangianParticle* particle, vtkFieldData* data) -{ - vtkIntArray::SafeDownCast(data->GetArray("Interaction")) - ->InsertNextValue(particle->GetInteraction()); -} - -//--------------------------------------------------------------------------- -void vtkLagrangianParticleTracker::InsertParticleData( - vtkLagrangianParticle* particle, vtkFieldData* data, int stepEnum) -{ - switch (stepEnum) - { - case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_PREV: - vtkIntArray::SafeDownCast(data->GetArray("StepNumber")) - ->InsertNextValue(particle->GetNumberOfSteps() - 1); - data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetPrevVelocity()); - data->GetArray("IntegrationTime")->InsertNextTuple1(particle->GetPrevIntegrationTime()); - break; - case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_CURRENT: - vtkIntArray::SafeDownCast(data->GetArray("StepNumber")) - ->InsertNextValue(particle->GetNumberOfSteps()); - data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetVelocity()); - data->GetArray("IntegrationTime")->InsertNextTuple1(particle->GetIntegrationTime()); - break; - case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_NEXT: - vtkIntArray::SafeDownCast(data->GetArray("StepNumber")) - ->InsertNextValue(particle->GetNumberOfSteps() + 1); - data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetNextVelocity()); - data->GetArray("IntegrationTime") - ->InsertNextTuple1(particle->GetIntegrationTime() + particle->GetStepTimeRef()); - break; - default: - break; - } + this->IntegrationModel->InsertSeedData(particle, pointData); } //--------------------------------------------------------------------------- @@ -1446,12 +1289,12 @@ double vtkLagrangianParticleTracker::ComputeCellLength(vtkLagrangianParticle* pa //--------------------------------------------------------------------------- bool vtkLagrangianParticleTracker::ComputeNextStep(vtkInitialValueProblemSolver* integrator, double* xprev, double* xnext, double t, double& delT, double& delTActual, double minStep, - double maxStep, int& integrationRes, vtkLagrangianParticle* particle) + double maxStep, double cellLength, int& integrationRes, vtkLagrangianParticle* particle) { // Check for potential manual integration double error; if (!this->IntegrationModel->ManualIntegration(xprev, xnext, t, delT, delTActual, minStep, - maxStep, this->IntegrationModel->GetTolerance(), error, integrationRes)) + maxStep, this->IntegrationModel->GetTolerance(), cellLength, error, integrationRes)) { // integrate one step integrationRes = integrator->ComputeNextStep(xprev, xnext, t, delT, delTActual, minStep, diff --git a/Filters/FlowPaths/vtkLagrangianParticleTracker.h b/Filters/FlowPaths/vtkLagrangianParticleTracker.h index 5c012342216b627c3d27d957385768e6cc5e3a0d..2b91d2b42ca445000a57870282887bd10f78d07c 100644 --- a/Filters/FlowPaths/vtkLagrangianParticleTracker.h +++ b/Filters/FlowPaths/vtkLagrangianParticleTracker.h @@ -238,25 +238,12 @@ public: //@{ /** - * Set/Get the Optional Paths Rendering feature, - * it allows to not show the particle paths - * if there is too many points. - * Default is false. + * Set/Get the generation of the particle path output, + * Default is true. */ - vtkSetMacro(UseParticlePathsRenderingThreshold, bool); - vtkGetMacro(UseParticlePathsRenderingThreshold, bool); - vtkBooleanMacro(UseParticlePathsRenderingThreshold, bool); - //@} - - //@{ - /** - * Set/Get the Optional Paths Rendering threshold, - * ie the maximum number of points to show in the particle - * path if the option is activated. - * Default is 100 - */ - vtkSetMacro(ParticlePathsRenderingPointsThreshold, int); - vtkGetMacro(ParticlePathsRenderingPointsThreshold, int); + vtkSetMacro(GenerateParticlePathsOutput, bool); + vtkGetMacro(GenerateParticlePathsOutput, bool); + vtkBooleanMacro(GenerateParticlePathsOutput, bool); //@} //@{ @@ -349,10 +336,6 @@ protected: virtual bool InitializeInteractionOutput(vtkInformationVector* outputVector, vtkPointData* seedData, vtkDataObject* surfaces, vtkDataObject*& interractionOutput); - virtual void InitializeParticleData(vtkFieldData* particleData, int maxTuples = 0); - virtual void InitializePathData(vtkFieldData* data); - virtual void InitializeInteractionData(vtkFieldData* data); - virtual bool FinalizeOutputs(vtkPolyData* particlePathsOutput, vtkDataObject* interractionOutput); static void InsertPolyVertexCell(vtkPolyData* polydata); @@ -373,21 +356,14 @@ protected: void InsertInteractionOutputPoint(vtkLagrangianParticle* particle, unsigned int interactedSurfaceFlatIndex, vtkDataObject* interactionOutput); - void InsertSeedData(vtkLagrangianParticle* particle, vtkFieldData* data); - void InsertPathData(vtkLagrangianParticle* particle, vtkFieldData* data); - void InsertInteractionData(vtkLagrangianParticle* particle, vtkFieldData* data); - void InsertParticleData(vtkLagrangianParticle* particle, vtkFieldData* data, int stepEnum); - double ComputeCellLength(vtkLagrangianParticle* particle); /** * This method is thread safe */ bool ComputeNextStep(vtkInitialValueProblemSolver* integrator, double* xprev, double* xnext, - double t, double& delT, double& delTActual, double minStep, double maxStep, int& integrationRes, - vtkLagrangianParticle* particle); - - virtual bool CheckParticlePathsRenderingThreshold(vtkPolyData* particlePathsOutput); + double t, double& delT, double& delTActual, double minStep, double maxStep, double cellLength, + int& integrationRes, vtkLagrangianParticle* particle); vtkLagrangianBasicIntegrationModel* IntegrationModel; vtkInitialValueProblemSolver* Integrator; @@ -399,9 +375,8 @@ protected: int MaximumNumberOfSteps; double MaximumIntegrationTime; bool AdaptiveStepReintegration; - bool UseParticlePathsRenderingThreshold; + bool GenerateParticlePathsOutput = true; bool GeneratePolyVertexInteractionOutput; - int ParticlePathsRenderingPointsThreshold; std::atomic<vtkIdType> ParticleCounter; std::atomic<vtkIdType> IntegratedParticleCounter; vtkIdType IntegratedParticleCounterIncrement; diff --git a/Filters/ParallelFlowPaths/Testing/Cxx/TestPLagrangianParticleTracker.cxx b/Filters/ParallelFlowPaths/Testing/Cxx/TestPLagrangianParticleTracker.cxx index 452392554234c69347be401b41b9d0877a53de42..f2acd6fe4d2cab2469f62c4aaaf2d2399f377225 100644 --- a/Filters/ParallelFlowPaths/Testing/Cxx/TestPLagrangianParticleTracker.cxx +++ b/Filters/ParallelFlowPaths/Testing/Cxx/TestPLagrangianParticleTracker.cxx @@ -165,6 +165,7 @@ void MainPLagrangianParticleTracker(vtkMultiProcessController *controller, void vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDiameter"); integrationModel->SetInputArrayToProcess(7, 1, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "ParticleDensity"); + integrationModel->SetNumberOfTrackedUserData(17); // Put in tracker vtkNew<vtkLagrangianParticleTracker> tracker; diff --git a/Filters/ParallelFlowPaths/vtkPLagrangianParticleTracker.cxx b/Filters/ParallelFlowPaths/vtkPLagrangianParticleTracker.cxx index 1db4e48a57fb0ba9e246e7210c26ab8e8322ae2d..aa74a72701739a199da22156552569128d208139 100644 --- a/Filters/ParallelFlowPaths/vtkPLagrangianParticleTracker.cxx +++ b/Filters/ParallelFlowPaths/vtkPLagrangianParticleTracker.cxx @@ -125,8 +125,8 @@ public: // Compute StreamSize for one particle // This is strongly linked to Send and Receive code - this->StreamSize = sizeof(int) * 2 + sizeof(double) * 2 + 4 * sizeof(vtkIdType) + - 2 * sizeof(double) + 3 * sizeof(double) * model->GetNumberOfIndependentVariables(); + this->StreamSize = sizeof(int) * 2 + sizeof(double) * 2 + 4 * sizeof(vtkIdType) + sizeof(int) + + 2 * sizeof(double) + 3 * sizeof(double) * (model->GetNumberOfIndependentVariables() + model->GetNumberOfTrackedUserData()); for (int i = 0; i < seedData->GetNumberOfArrays(); i++) { vtkDataArray* array = seedData->GetArray(i); @@ -162,6 +162,7 @@ public: *this->SendStream << particle->GetId(); *this->SendStream << particle->GetParentId(); *this->SendStream << particle->GetNumberOfVariables(); + *this->SendStream << static_cast<int>(particle->GetTrackedUserData().size()); *this->SendStream << particle->GetNumberOfSteps(); *this->SendStream << particle->GetIntegrationTime(); *this->SendStream << particle->GetPrevIntegrationTime(); @@ -179,6 +180,19 @@ public: *this->SendStream << next[i]; } + for (auto data : particle->GetPrevTrackedUserData()) + { + *this->SendStream << data; + } + for (auto data : particle->GetTrackedUserData()) + { + *this->SendStream << data; + } + for (auto data : particle->GetNextTrackedUserData()) + { + *this->SendStream << data; + } + for (int i = 0; i < particle->GetSeedData()->GetNumberOfArrays(); i++) { vtkDataArray* array = particle->GetSeedData()->GetArray(i); @@ -219,7 +233,7 @@ public: vtkMultiProcessController::ANY_SOURCE, LAGRANGIAN_PARTICLE_TAG); // Deserialize particle // This is strongly linked to Constructor and Send method - int nVar, userFlag; + int nVar, userFlag, nTrackedUserData; vtkIdType seedId, particleId, parentId, nSteps; double iTime, prevITime; bool pInsertPrevious, pManualShift; @@ -227,6 +241,7 @@ public: *this->ReceiveStream >> particleId; *this->ReceiveStream >> parentId; *this->ReceiveStream >> nVar; + *this->ReceiveStream >> nTrackedUserData; *this->ReceiveStream >> nSteps; *this->ReceiveStream >> iTime; *this->ReceiveStream >> prevITime; @@ -241,7 +256,7 @@ public: seedTupleIndex = this->SeedData->GetArray(0)->GetNumberOfTuples(); } particle = vtkLagrangianParticle::NewInstance(nVar, seedId, particleId, seedTupleIndex, iTime, - this->SeedData, this->WeightsSize, nSteps, prevITime); + this->SeedData, this->WeightsSize, nTrackedUserData, nSteps, prevITime); particle->SetParentId(parentId); particle->SetUserFlag(userFlag); particle->SetPInsertPreviousPosition(pInsertPrevious); @@ -256,6 +271,22 @@ public: *this->ReceiveStream >> next[i]; } + std::vector<double>& prevTracked = particle->GetPrevTrackedUserData(); + for (auto &var : prevTracked) + { + *this->ReceiveStream >> var; + } + std::vector<double>& tracked = particle->GetTrackedUserData(); + for (auto &var : tracked) + { + *this->ReceiveStream >> var; + } + std::vector<double>& nextTracked = particle->GetNextTrackedUserData(); + for (auto &var : nextTracked) + { + *this->ReceiveStream >> var; + } + for (int i = 0; i < particle->GetSeedData()->GetNumberOfArrays(); i++) { vtkDataArray* array = particle->GetSeedData()->GetArray(i); @@ -747,7 +778,7 @@ void vtkPLagrangianParticleTracker::GenerateParticles(const vtkBoundingBox* boun // Create and set a dummy particle so FindInLocators can use caching. vtkLagrangianParticle dummyParticle( - 0, 0, 0, 0, 0, nullptr, this->IntegrationModel->GetWeightsSize()); + 0, 0, 0, 0, 0, nullptr, this->IntegrationModel->GetWeightsSize(), 0); // Generate particle and distribute the ones not in domain to other nodes for (vtkIdType i = 0; i < seeds->GetNumberOfPoints(); i++) @@ -758,7 +789,7 @@ void vtkPLagrangianParticleTracker::GenerateParticles(const vtkBoundingBox* boun initialIntegrationTimes ? initialIntegrationTimes->GetTuple1(i) : 0; vtkIdType particleId = this->GetNewParticleId(); vtkLagrangianParticle* particle = new vtkLagrangianParticle(nVar, particleId, particleId, i, - initialIntegrationTime, seedData, this->IntegrationModel->GetWeightsSize()); + initialIntegrationTime, seedData, this->IntegrationModel->GetWeightsSize(), this->IntegrationModel->GetNumberOfTrackedUserData()); memcpy(particle->GetPosition(), position, 3 * sizeof(double)); initialVelocities->GetTuple(i, particle->GetVelocity()); this->IntegrationModel->InitializeParticle(particle); @@ -888,7 +919,7 @@ int vtkPLagrangianParticleTracker::Integrate(vtkInitialValueProblemSolver* integ { if (this->Controller && this->Controller->GetNumberOfProcesses() > 1) { - if (particle->GetPInsertPreviousPosition()) + if (particlePathsOutput && particle->GetPInsertPreviousPosition()) { // Mutex Locked Area std::lock_guard<std::mutex> guard(this->ParticlePathsOutputMutex); @@ -948,7 +979,7 @@ void vtkPLagrangianParticleTracker::ReceiveParticles( bool vtkPLagrangianParticleTracker::FinalizeOutputs( vtkPolyData* particlePathsOutput, vtkDataObject* interactionOutput) { - if (this->Controller && this->Controller->GetNumberOfProcesses() > 1) + if (particlePathsOutput && this->Controller && this->Controller->GetNumberOfProcesses() > 1) { // Construct array with all non outofdomains ids and terminations @@ -988,31 +1019,6 @@ bool vtkPLagrangianParticleTracker::FinalizeOutputs( return this->Superclass::FinalizeOutputs(particlePathsOutput, interactionOutput); } -//--------------------------------------------------------------------------- -bool vtkPLagrangianParticleTracker::CheckParticlePathsRenderingThreshold( - vtkPolyData* particlePathsOutput) -{ - if (this->Controller && this->Controller->GetNumberOfProcesses() > 1) - { - if (this->UseParticlePathsRenderingThreshold) - { - // Reduce the totalNumberOfPoints to check if we need to display the particle paths. - vtkIdType localNPoints = particlePathsOutput->GetNumberOfPoints(); - vtkIdType totalNPoints; - this->Controller->AllReduce(&localNPoints, &totalNPoints, 1, vtkCommunicator::SUM_OP); - return totalNPoints > this->ParticlePathsRenderingPointsThreshold; - } - else - { - return false; - } - } - else - { - return this->Superclass::CheckParticlePathsRenderingThreshold(particlePathsOutput); - } -} - //--------------------------------------------------------------------------- bool vtkPLagrangianParticleTracker::UpdateSurfaceCacheIfNeeded(vtkDataObject*& surfaces) { diff --git a/Filters/ParallelFlowPaths/vtkPLagrangianParticleTracker.h b/Filters/ParallelFlowPaths/vtkPLagrangianParticleTracker.h index efa76057d47f84ebc0877107608feac329e17764..9b8ce1331c394b6d34d7e3bfd8401c7a7cdec17d 100644 --- a/Filters/ParallelFlowPaths/vtkPLagrangianParticleTracker.h +++ b/Filters/ParallelFlowPaths/vtkPLagrangianParticleTracker.h @@ -92,8 +92,6 @@ protected: bool FinalizeOutputs( vtkPolyData* particlePathsOutput, vtkDataObject* interractionOutput) override; - bool CheckParticlePathsRenderingThreshold(vtkPolyData* particlePathsOutput) override; - bool UpdateSurfaceCacheIfNeeded(vtkDataObject*& surfaces) override; /**