From 3ce232d9c4738a47be0597736eddfd68255da4a9 Mon Sep 17 00:00:00 2001 From: Nicolas Vuaille <nicolas.vuaille@kitware.com> Date: Tue, 27 Mar 2018 17:20:28 +0200 Subject: [PATCH] Introduce Molecule Append filter * Introduce vtkAppendMolecule filter to append input molecules into one (with atomData and bondData) * Since vtkMolecule is in Common/DataModel, move vtkMoleculeAlgorithm to Common/ExecutionModel --- Common/DataModel/vtkMolecule.cxx | 123 ++++---- Common/DataModel/vtkMolecule.h | 72 ++++- Common/ExecutionModel/CMakeLists.txt | 1 + .../ExecutionModel}/vtkMoleculeAlgorithm.cxx | 0 .../ExecutionModel}/vtkMoleculeAlgorithm.h | 4 +- Domains/Chemistry/CMakeLists.txt | 1 - Filters/Core/CMakeLists.txt | 1 + Filters/Core/Testing/Cxx/CMakeLists.txt | 1 + .../Core/Testing/Cxx/TestAppendMolecule.cxx | 283 ++++++++++++++++++ Filters/Core/vtkMoleculeAppend.cxx | 282 +++++++++++++++++ Filters/Core/vtkMoleculeAppend.h | 82 +++++ 11 files changed, 782 insertions(+), 68 deletions(-) rename {Domains/Chemistry => Common/ExecutionModel}/vtkMoleculeAlgorithm.cxx (100%) rename {Domains/Chemistry => Common/ExecutionModel}/vtkMoleculeAlgorithm.h (96%) create mode 100644 Filters/Core/Testing/Cxx/TestAppendMolecule.cxx create mode 100644 Filters/Core/vtkMoleculeAppend.cxx create mode 100644 Filters/Core/vtkMoleculeAppend.h diff --git a/Common/DataModel/vtkMolecule.cxx b/Common/DataModel/vtkMolecule.cxx index b1697621283..5b46c40df42 100644 --- a/Common/DataModel/vtkMolecule.cxx +++ b/Common/DataModel/vtkMolecule.cxx @@ -18,6 +18,8 @@ PURPOSE. See the above copyright notice for more information. #include "vtkEdgeListIterator.h" #include "vtkGraphInternals.h" #include "vtkIdTypeArray.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" #include "vtkMatrix3x3.h" #include "vtkNew.h" #include "vtkObjectFactory.h" @@ -40,7 +42,9 @@ vtkMolecule::vtkMolecule() Lattice(nullptr), LatticeOrigin(0., 0., 0.), AtomGhostArray(nullptr), - BondGhostArray(nullptr) + BondGhostArray(nullptr), + AtomicNumberArrayName(nullptr), + BondOrdersArrayName(nullptr) { this->Initialize(); } @@ -56,9 +60,10 @@ void vtkMolecule::Initialize() vertData->AllocateArrays(1); // atomic nums // Atomic numbers + this->SetAtomicNumberArrayName("Atomic Numbers"); vtkNew<vtkUnsignedShortArray> atomicNums; atomicNums->SetNumberOfComponents(1); - atomicNums->SetName(vtkMolecule::GetAtomicNumberArrayName()); + atomicNums->SetName(this->GetAtomicNumberArrayName()); vertData->SetScalars(atomicNums); // Nuclear coordinates @@ -70,9 +75,10 @@ void vtkMolecule::Initialize() vtkDataSetAttributes *edgeData = this->GetEdgeData(); edgeData->AllocateArrays(1); // Bond orders + this->SetBondOrdersArrayName("Bond Orders"); vtkNew<vtkUnsignedShortArray> bondOrders; bondOrders->SetNumberOfComponents(1); - bondOrders->SetName("Bond Orders"); + bondOrders->SetName(this->GetBondOrdersArrayName()); edgeData->SetScalars(bondOrders); this->UpdateBondList(); @@ -87,6 +93,8 @@ void vtkMolecule::Initialize() vtkMolecule::~vtkMolecule() { this->SetElectronicData(nullptr); + delete [] this->AtomicNumberArrayName; + delete [] this->BondOrdersArrayName; } //---------------------------------------------------------------------------- @@ -131,6 +139,9 @@ void vtkMolecule::PrintSelf(ostream &os, vtkIndent indent) { os << subIndent << "Not set.\n"; } + + os << indent << "Atomic number array name : " << this->GetAtomicNumberArrayName() << "\n"; + os << indent << "Bond orders array name : " << this->GetBondOrdersArrayName(); } //---------------------------------------------------------------------------- @@ -228,8 +239,8 @@ vtkIdType vtkMolecule::GetNumberOfAtoms() vtkBond vtkMolecule::AppendBond(const vtkIdType atom1, const vtkIdType atom2, const unsigned short order) { - vtkUnsignedShortArray *bondOrders = vtkArrayDownCast<vtkUnsignedShortArray> - (this->GetEdgeData()->GetScalars()); + vtkUnsignedShortArray *bondOrders = this->GetBondOrdersArray(); + assert(bondOrders); vtkEdgeType edgeType; @@ -258,9 +269,7 @@ void vtkMolecule::SetBondOrder(vtkIdType bondId, unsigned short order) { assert(bondId >= 0 && bondId < this->GetNumberOfBonds()); - vtkUnsignedShortArray *bondOrders = vtkArrayDownCast<vtkUnsignedShortArray> - (this->GetEdgeData()->GetScalars()); - + vtkUnsignedShortArray *bondOrders = this->GetBondOrdersArray(); assert(bondOrders); this->Modified(); @@ -272,12 +281,9 @@ unsigned short vtkMolecule::GetBondOrder(vtkIdType bondId) { assert(bondId >= 0 && bondId < this->GetNumberOfBonds()); - vtkUnsignedShortArray *bondOrders = vtkArrayDownCast<vtkUnsignedShortArray> - (this->GetEdgeData()->GetScalars()); - - assert(bondOrders); + vtkUnsignedShortArray *bondOrders = this->GetBondOrdersArray(); - return bondOrders->GetValue(bondId); + return bondOrders ? bondOrders->GetValue(bondId) : 0; } //---------------------------------------------------------------------------- @@ -307,13 +313,20 @@ vtkPoints * vtkMolecule::GetAtomicPositionArray() vtkUnsignedShortArray * vtkMolecule::GetAtomicNumberArray() { vtkUnsignedShortArray *atomicNums = vtkArrayDownCast<vtkUnsignedShortArray> - (this->GetVertexData()->GetScalars(vtkMolecule::GetAtomicNumberArrayName())); + (this->GetVertexData()->GetScalars(this->GetAtomicNumberArrayName())); assert(atomicNums); return atomicNums; } +//---------------------------------------------------------------------------- +vtkUnsignedShortArray * vtkMolecule::GetBondOrdersArray() +{ + return vtkArrayDownCast<vtkUnsignedShortArray> + (this->GetBondData()->GetScalars(this->GetBondOrdersArrayName())); +} + //---------------------------------------------------------------------------- vtkIdType vtkMolecule::GetNumberOfBonds() { @@ -637,6 +650,10 @@ void vtkMolecule::AllocateAtomGhostArray() ghosts->FillComponent(0, 0); this->GetVertexData()->AddArray(ghosts); } + else + { + this->GetAtomGhostArray()->SetNumberOfTuples(this->GetNumberOfAtoms()); + } } //---------------------------------------------------------------------------- @@ -658,25 +675,50 @@ void vtkMolecule::AllocateBondGhostArray() ghosts->FillComponent(0, 0); this->GetEdgeData()->AddArray(ghosts); } + else + { + this->GetBondGhostArray()->SetNumberOfTuples(this->GetNumberOfBonds()); + } } //---------------------------------------------------------------------------- int vtkMolecule::Initialize(vtkPoints* atomPositions, - vtkUnsignedShortArray* atomicNumberArray, + vtkDataArray* atomicNumberArray, vtkDataSetAttributes* atomData) { // Start with default initialization the molecule this->Initialize(); + // if no atomicNumberArray given, look for one in atomData + if (!atomicNumberArray && atomData) + { + atomicNumberArray = atomData->GetArray(this->GetAtomicNumberArrayName()); + } + + // ensure it is a short array + if (atomicNumberArray && !vtkUnsignedShortArray::SafeDownCast(atomicNumberArray)) + { + vtkNew<vtkUnsignedShortArray> newAtomicNumberShortArray; + vtkIdType nbPoints = atomicNumberArray->GetNumberOfTuples(); + newAtomicNumberShortArray->SetNumberOfComponents(1); + newAtomicNumberShortArray->SetNumberOfTuples(nbPoints); + newAtomicNumberShortArray->SetName(atomicNumberArray->GetName()); + for (vtkIdType i = 0; i < nbPoints; i++) + { + newAtomicNumberShortArray->SetTuple1(i, atomicNumberArray->GetTuple1(i)); + } + atomicNumberArray->ShallowCopy(newAtomicNumberShortArray); + } + if (!atomPositions && !atomicNumberArray) { - vtkWarningMacro(<< "Atom position and atomic numbers were not found: skip atomic data."); + vtkDebugMacro(<< "Atom position and atomic numbers were not found: skip atomic data."); return 1; } if (!atomPositions || !atomicNumberArray) { - vtkErrorMacro(<< "Empty atoms or atomic numbers."); + vtkDebugMacro(<< "Empty atoms or atomic numbers."); return 0; } @@ -692,7 +734,7 @@ int vtkMolecule::Initialize(vtkPoints* atomPositions, return 0; } - static const std::string atomicNumberName = vtkMolecule::GetAtomicNumberArrayName(); + static const std::string atomicNumberName = this->GetAtomicNumberArrayName(); // update atoms positions this->ForceOwnership(); @@ -749,51 +791,26 @@ int vtkMolecule::Initialize(vtkPoints* atomPositions, } //---------------------------------------------------------------------------- -int vtkMolecule::Initialize(vtkPoints* atomPositions, - vtkDataArray* atomicNumberArray, - vtkDataSetAttributes* atomData) +int vtkMolecule::Initialize(vtkMolecule* molecule) { - vtkUnsignedShortArray* atomicNumberShortArray = - vtkUnsignedShortArray::SafeDownCast(atomicNumberArray); - vtkNew<vtkUnsignedShortArray> newAtomicNumberShortArray; - if (!atomicNumberShortArray && atomicNumberArray) + if (molecule == nullptr) { - vtkIdType nbPoints = atomicNumberArray->GetNumberOfTuples(); - newAtomicNumberShortArray->SetNumberOfComponents(1); - newAtomicNumberShortArray->SetNumberOfTuples(nbPoints); - newAtomicNumberShortArray->SetName(atomicNumberArray->GetName()); - for (vtkIdType i = 0; i < nbPoints; i++) - { - newAtomicNumberShortArray->SetTuple1(i, atomicNumberArray->GetTuple1(i)); - } - atomicNumberShortArray = newAtomicNumberShortArray; + this->Initialize(); + return 1; } - return this->Initialize(atomPositions, atomicNumberShortArray, atomData); + return this->Initialize( + molecule->GetPoints(), molecule->GetAtomicNumberArray(), molecule->GetVertexData()); } //---------------------------------------------------------------------------- -int vtkMolecule::Initialize(vtkPoints* atomPositions, vtkDataSetAttributes* atomData) +vtkMolecule *vtkMolecule::GetData(vtkInformation *info) { - vtkDataArray* atomicNumberArray = atomData->GetArray(vtkMolecule::GetAtomicNumberArrayName()); - if (!atomicNumberArray) - { - atomicNumberArray = atomData->GetScalars(); - } - - return this->Initialize(atomPositions, atomicNumberArray, atomData); + return info? vtkMolecule::SafeDownCast(info->Get(DATA_OBJECT())) : nullptr; } //---------------------------------------------------------------------------- -int vtkMolecule::Initialize(vtkMolecule* molecule) +vtkMolecule *vtkMolecule::GetData(vtkInformationVector *v, int i) { - if (molecule == nullptr) - { - this->Initialize(); - vtkErrorMacro(<< "No input molecule."); - return 0; - } - - return this->Initialize( - molecule->GetPoints(), molecule->GetAtomicNumberArray(), molecule->GetVertexData()); + return vtkMolecule::GetData(v->GetInformationObject(i)); } diff --git a/Common/DataModel/vtkMolecule.h b/Common/DataModel/vtkMolecule.h index 4e7f798f360..cc573888488 100644 --- a/Common/DataModel/vtkMolecule.h +++ b/Common/DataModel/vtkMolecule.h @@ -80,6 +80,8 @@ class vtkAbstractElectronicData; class vtkDataArray; +class vtkInformation; +class vtkInformationVector; class vtkMatrix3x3; class vtkPlane; class vtkPoints; @@ -210,6 +212,7 @@ public: */ vtkPoints * GetAtomicPositionArray(); vtkUnsignedShortArray * GetAtomicNumberArray(); + vtkUnsignedShortArray * GetBondOrdersArray(); //@} //@{ @@ -370,30 +373,71 @@ public: * Initialize a molecule with an atom per input point. * Parameters atomPositions and atomicNumberArray should have the same size. */ - int Initialize(vtkPoints* atomPositions, - vtkUnsignedShortArray* atomicNumberArray, - vtkDataSetAttributes* atomData = nullptr); - - /** - * Overloads Initialize method to allow vtkDataArray instead of vtkUnsignedShortArray. - */ int Initialize(vtkPoints* atomPositions, vtkDataArray* atomicNumberArray, - vtkDataSetAttributes* atomData = nullptr); + vtkDataSetAttributes* atomData); /** - * Overloads Initialize method. Look for an atomic number array in atomData. - * If none found, take the first array. + * Overloads Initialize method. */ int Initialize(vtkPoints* atomPositions, - vtkDataSetAttributes* atomData); + vtkDataSetAttributes* atomData) + { + return this->Initialize(atomPositions, nullptr, atomData); + } /** * Use input molecule points, atomic number and atomic data to initialize the new molecule. */ int Initialize(vtkMolecule* molecule); - static const char* GetAtomicNumberArrayName() {return "Atomic Numbers";} + //@{ + /** + * Retrieve a molecule from an information vector. + */ + static vtkMolecule* GetData(vtkInformation *info); + static vtkMolecule* GetData(vtkInformationVector *v, int i=0); + //@} + + /** + * Return the VertexData of the underlying graph + */ + vtkDataSetAttributes* GetAtomData() + { + return this->GetVertexData(); + } + + /** + * Return the EdgeData of the underlying graph + */ + vtkDataSetAttributes* GetBondData() + { + return this->GetEdgeData(); + } + + /** + * Return the edge id from the underlying graph. + */ + vtkIdType GetBondId(vtkIdType a, vtkIdType b) + { + return this->GetEdgeId(a, b); + } + + //@{ + /** + * Get/Set the atomic number array name. + */ + vtkSetStringMacro(AtomicNumberArrayName); + vtkGetStringMacro(AtomicNumberArrayName); + //@} + + //@{ + /** + * Get/Set the bond orders array name. + */ + vtkSetStringMacro(BondOrdersArrayName); + vtkGetStringMacro(BondOrdersArrayName); + //@} protected: vtkMolecule(); @@ -431,6 +475,10 @@ public: vtkUnsignedCharArray* AtomGhostArray; vtkUnsignedCharArray* BondGhostArray; + + char* AtomicNumberArrayName; + char* BondOrdersArrayName; + private: vtkMolecule(const vtkMolecule&) = delete; void operator=(const vtkMolecule&) = delete; diff --git a/Common/ExecutionModel/CMakeLists.txt b/Common/ExecutionModel/CMakeLists.txt index 5eb4d7a1e4c..6f6bb71adfd 100644 --- a/Common/ExecutionModel/CMakeLists.txt +++ b/Common/ExecutionModel/CMakeLists.txt @@ -28,6 +28,7 @@ SET(Module_SRCS vtkInformationExecutivePortKey.cxx vtkInformationExecutivePortVectorKey.cxx vtkInformationIntegerRequestKey.cxx + vtkMoleculeAlgorithm.cxx vtkMultiBlockDataSetAlgorithm.cxx vtkMultiTimeStepAlgorithm.cxx vtkPassInputTypeAlgorithm.cxx diff --git a/Domains/Chemistry/vtkMoleculeAlgorithm.cxx b/Common/ExecutionModel/vtkMoleculeAlgorithm.cxx similarity index 100% rename from Domains/Chemistry/vtkMoleculeAlgorithm.cxx rename to Common/ExecutionModel/vtkMoleculeAlgorithm.cxx diff --git a/Domains/Chemistry/vtkMoleculeAlgorithm.h b/Common/ExecutionModel/vtkMoleculeAlgorithm.h similarity index 96% rename from Domains/Chemistry/vtkMoleculeAlgorithm.h rename to Common/ExecutionModel/vtkMoleculeAlgorithm.h index 4766d419c3c..6a513e77d25 100644 --- a/Domains/Chemistry/vtkMoleculeAlgorithm.h +++ b/Common/ExecutionModel/vtkMoleculeAlgorithm.h @@ -33,13 +33,13 @@ #ifndef vtkMoleculeAlgorithm_h #define vtkMoleculeAlgorithm_h -#include "vtkDomainsChemistryModule.h" // For export macro +#include "vtkCommonExecutionModelModule.h" // For export macro #include "vtkAlgorithm.h" class vtkDataSet; class vtkMolecule; -class VTKDOMAINSCHEMISTRY_EXPORT vtkMoleculeAlgorithm : public vtkAlgorithm +class VTKCOMMONEXECUTIONMODEL_EXPORT vtkMoleculeAlgorithm : public vtkAlgorithm { public: static vtkMoleculeAlgorithm *New(); diff --git a/Domains/Chemistry/CMakeLists.txt b/Domains/Chemistry/CMakeLists.txt index 04944134b60..5370210eea0 100644 --- a/Domains/Chemistry/CMakeLists.txt +++ b/Domains/Chemistry/CMakeLists.txt @@ -4,7 +4,6 @@ set(Module_SRCS vtkBlueObeliskDataParser.cxx vtkCMLMoleculeReader.cxx vtkGaussianCubeReader2.cxx - vtkMoleculeAlgorithm.cxx vtkMoleculeMapper.cxx vtkMoleculeToAtomBallFilter.cxx vtkMoleculeToBondStickFilter.cxx diff --git a/Filters/Core/CMakeLists.txt b/Filters/Core/CMakeLists.txt index 06b7e3fd678..dbd60ba6e3c 100644 --- a/Filters/Core/CMakeLists.txt +++ b/Filters/Core/CMakeLists.txt @@ -44,6 +44,7 @@ set(Module_SRCS vtkMergeDataObjectFilter.cxx vtkMergeFields.cxx vtkMergeFilter.cxx + vtkMoleculeAppend.cxx vtkMultiObjectMassProperties.cxx vtkPlaneCutter.cxx vtkPointDataToCellData.cxx diff --git a/Filters/Core/Testing/Cxx/CMakeLists.txt b/Filters/Core/Testing/Cxx/CMakeLists.txt index 9ddef7729d7..09642dd7838 100644 --- a/Filters/Core/Testing/Cxx/CMakeLists.txt +++ b/Filters/Core/Testing/Cxx/CMakeLists.txt @@ -1,6 +1,7 @@ vtk_add_test_cxx(vtkFiltersCoreCxxTests tests TestAppendArcLength.cxx,NO_VALID TestAppendFilter.cxx,NO_VALID + TestAppendMolecule.cxx,NO_VALID TestAppendPolyData.cxx,NO_VALID TestAppendSelection.cxx,NO_VALID TestArrayCalculator.cxx,NO_VALID diff --git a/Filters/Core/Testing/Cxx/TestAppendMolecule.cxx b/Filters/Core/Testing/Cxx/TestAppendMolecule.cxx new file mode 100644 index 00000000000..9de1d199796 --- /dev/null +++ b/Filters/Core/Testing/Cxx/TestAppendMolecule.cxx @@ -0,0 +1,283 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: TestAppendMolecule.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "vtkDataSetAttributes.h" +#include "vtkDoubleArray.h" +#include "vtkMolecule.h" +#include "vtkMoleculeAppend.h" +#include "vtkNew.h" +#include "vtkUnsignedCharArray.h" +#include "vtkUnsignedShortArray.h" + +#define CheckNumbers(name, first, second) \ + if (first != second) \ + { \ + cerr << "Error : wrong number of " << #name << ". Got " << first << " but expects " << second \ + << endl; \ + return EXIT_FAILURE; \ + } + +// Used to creates different atoms and data for each molecule +static int NB_OF_MOL = 0; + +// Add two atoms with a bond between them. +void InitSimpleMolecule(vtkMolecule* molecule) +{ + NB_OF_MOL++; + vtkAtom h1 = molecule->AppendAtom(1, 0.5, 1.5, -NB_OF_MOL); + vtkAtom h2 = molecule->AppendAtom(1, 0.5, 1.5, NB_OF_MOL); + molecule->AppendBond(h1, h2, 1); +} + +void AddAtomData(vtkMolecule* molecule, vtkIdType size) +{ + vtkNew<vtkDoubleArray> data; + data->SetName("Data"); + data->SetNumberOfComponents(1); + for (vtkIdType i = 0; i < size; i++) + { + data->InsertNextValue(NB_OF_MOL * 1.01); + } + molecule->GetAtomData()->AddArray(data); +} + +int CheckMolecule(vtkMolecule* molecule, + int nbAtoms, + int nbBonds, + int nbArrays, + vtkDoubleArray* values, + int nbGhostAtoms, + int nbGhostBonds) +{ + CheckNumbers("atoms", molecule->GetNumberOfAtoms(), nbAtoms); + CheckNumbers("bonds", molecule->GetNumberOfBonds(), nbBonds); + CheckNumbers("atom data arrays", molecule->GetAtomData()->GetNumberOfArrays(), nbArrays); + + vtkDataArray* resultData = molecule->GetAtomData()->GetArray("Data"); + if (!resultData) + { + std::cerr << "Error : atoms data array not found in result" << std::endl; + return EXIT_FAILURE; + } + CheckNumbers("atom data array values", resultData->GetNumberOfTuples(), nbAtoms); + + for (vtkIdType i = 0; i < nbAtoms; i++) + { + CheckNumbers("data value", resultData->GetTuple1(i), values->GetValue(i)); + } + + vtkUnsignedShortArray* bondOrderArray = molecule->GetBondOrdersArray(); + if (!bondOrderArray) + { + std::cerr << "Error : bonds data array not found in result" << std::endl; + return EXIT_FAILURE; + } + CheckNumbers("bond data array values", bondOrderArray->GetNumberOfTuples(), nbBonds); + + vtkUnsignedCharArray* ghostAtoms = molecule->GetAtomGhostArray(); + int nbOfGhosts = 0; + for (vtkIdType id = 0; id < nbAtoms; id++) + { + if (ghostAtoms->GetValue(id) == 1) + { + nbOfGhosts++; + } + } + // ghost atom from molecule2 is still ghost in result + CheckNumbers("ghost atoms", nbOfGhosts, nbGhostAtoms); + + vtkUnsignedCharArray* ghostBonds = molecule->GetBondGhostArray(); + nbOfGhosts = 0; + for (vtkIdType id = 0; id < nbBonds; id++) + { + if (ghostBonds->GetValue(id) == 1) + { + nbOfGhosts++; + } + } + // ghost bond from molecule2 is still ghost in result + CheckNumbers("ghost bonds", nbOfGhosts, nbGhostBonds); + + return EXIT_SUCCESS; +} + +int TestAppendMolecule(int, char* []) +{ + // -------------------------------------------------------------------------- + // Simple test : 2 molecules, no data + vtkNew<vtkMolecule> simpleMolecule1; + InitSimpleMolecule(simpleMolecule1); + + vtkNew<vtkMolecule> simpleMolecule2; + InitSimpleMolecule(simpleMolecule2); + + vtkNew<vtkMoleculeAppend> appender; + appender->AddInputData(simpleMolecule1); + appender->AddInputData(simpleMolecule2); + appender->Update(); + vtkMolecule* resultMolecule = appender->GetOutput(); + + int expectedResult = simpleMolecule1->GetNumberOfAtoms() + simpleMolecule2->GetNumberOfAtoms(); + CheckNumbers("atoms", resultMolecule->GetNumberOfAtoms(), expectedResult); + + expectedResult = simpleMolecule1->GetNumberOfBonds() + simpleMolecule2->GetNumberOfBonds(); + CheckNumbers("bonds", resultMolecule->GetNumberOfBonds(), expectedResult); + + // -------------------------------------------------------------------------- + // Full test : ghosts and data + + /** + * Use 3 molecules: + * - fullMolecule1 : 2 atoms and one bond, no ghost + * - fullMolecule2 : 3 atoms and 2 bonds, one ghost atom and one ghost bond + * - fullMolecule3 : 3 atoms and 2 bonds, one ghost atom and one ghost bond + */ + + // INIT + vtkNew<vtkMolecule> fullMolecule1; + InitSimpleMolecule(fullMolecule1); + AddAtomData(fullMolecule1, 2); + vtkNew<vtkMolecule> fullMolecule2; + InitSimpleMolecule(fullMolecule2); + AddAtomData(fullMolecule2, 3); + vtkNew<vtkMolecule> fullMolecule3; + InitSimpleMolecule(fullMolecule3); + AddAtomData(fullMolecule3, 3); + + // duplicate first atom of molecule 2 to be ghost in molecule 3, and vice versa. + vtkAtom firstAtom2 = fullMolecule2->GetAtom(0); + vtkAtom firstAtom3 = fullMolecule3->GetAtom(0); + + vtkAtom ghostAtom2 = + fullMolecule2->AppendAtom(firstAtom3.GetAtomicNumber(), firstAtom3.GetPosition()); + vtkBond ghostBond2 = fullMolecule2->AppendBond(firstAtom2, ghostAtom2, 1); + + vtkAtom ghostAtom3 = + fullMolecule3->AppendAtom(firstAtom2.GetAtomicNumber(), firstAtom2.GetPosition()); + vtkBond ghostBond3 = fullMolecule3->AppendBond(firstAtom3, ghostAtom3, 1); + + // set ghost flag on relevant atoms and bonds. + fullMolecule1->AllocateAtomGhostArray(); + fullMolecule1->AllocateBondGhostArray(); + + fullMolecule2->AllocateAtomGhostArray(); + fullMolecule2->GetAtomGhostArray()->SetValue(ghostAtom2.GetId(), 1); + fullMolecule2->AllocateBondGhostArray(); + fullMolecule2->GetBondGhostArray()->SetValue(ghostBond2.GetId(), 1); + + fullMolecule3->AllocateAtomGhostArray(); + fullMolecule3->GetAtomGhostArray()->SetValue(ghostAtom3.GetId(), 1); + fullMolecule3->AllocateBondGhostArray(); + fullMolecule3->GetBondGhostArray()->SetValue(ghostBond3.GetId(), 1); + + // -------------------------------------------------------------------------- + // First part: 2 molecules, ghost and data + + // APPEND 1 with 2 + vtkNew<vtkMoleculeAppend> appender2; + appender2->AddInputData(fullMolecule1); + appender2->AddInputData(fullMolecule2); + appender2->Update(); + vtkMolecule* resultFullMolecule = appender2->GetOutput(); + + // CHECK RESULT + int nbOfExpectedAtoms = fullMolecule1->GetNumberOfAtoms() + fullMolecule2->GetNumberOfAtoms(); + int nbOfExpectedBonds = fullMolecule1->GetNumberOfBonds() + fullMolecule2->GetNumberOfBonds(); + int nbOfExpectedArrays = fullMolecule1->GetAtomData()->GetNumberOfArrays(); + vtkNew<vtkDoubleArray> expectedResultValues; + expectedResultValues->InsertNextValue( + fullMolecule1->GetAtomData()->GetArray("Data")->GetTuple1(0)); + expectedResultValues->InsertNextValue( + fullMolecule1->GetAtomData()->GetArray("Data")->GetTuple1(1)); + expectedResultValues->InsertNextValue( + fullMolecule2->GetAtomData()->GetArray("Data")->GetTuple1(0)); + expectedResultValues->InsertNextValue( + fullMolecule2->GetAtomData()->GetArray("Data")->GetTuple1(1)); + expectedResultValues->InsertNextValue( + fullMolecule2->GetAtomData()->GetArray("Data")->GetTuple1(2)); + + int res = CheckMolecule(resultFullMolecule, + nbOfExpectedAtoms, + nbOfExpectedBonds, + nbOfExpectedArrays, + expectedResultValues, + 1, + 1); + if (res == EXIT_FAILURE) + { + return EXIT_FAILURE; + } + + // -------------------------------------------------------------------------- + // Second part: 3 molecules, ghost and data, no merge + // APPEND third molecule + appender2->MergeCoincidentAtomsOff(); + appender2->AddInputData(fullMolecule3); + appender2->Update(); + resultFullMolecule = appender2->GetOutput(); + + // CHECK RESULT + nbOfExpectedAtoms = fullMolecule1->GetNumberOfAtoms() + fullMolecule2->GetNumberOfAtoms() + + fullMolecule3->GetNumberOfAtoms(); + nbOfExpectedBonds = fullMolecule1->GetNumberOfBonds() + fullMolecule2->GetNumberOfBonds() + + fullMolecule3->GetNumberOfBonds(); + nbOfExpectedArrays = fullMolecule1->GetAtomData()->GetNumberOfArrays(); + + // Result contains data of non ghost atom. + expectedResultValues->InsertNextValue( + fullMolecule3->GetAtomData()->GetArray("Data")->GetTuple1(0)); + expectedResultValues->InsertNextValue( + fullMolecule3->GetAtomData()->GetArray("Data")->GetTuple1(1)); + expectedResultValues->InsertNextValue( + fullMolecule3->GetAtomData()->GetArray("Data")->GetTuple1(2)); + + res = CheckMolecule(resultFullMolecule, + nbOfExpectedAtoms, + nbOfExpectedBonds, + nbOfExpectedArrays, + expectedResultValues, + 2, + 2); + if (res == EXIT_FAILURE) + { + return EXIT_FAILURE; + } + + // -------------------------------------------------------------------------- + // Third part: 3 molecules, ghost and data, merge coincident atoms + + appender2->MergeCoincidentAtomsOn(); + appender2->Update(); + resultFullMolecule = appender2->GetOutput(); + // the ghost atoms are not duplicated in output. + nbOfExpectedAtoms = fullMolecule1->GetNumberOfAtoms() + fullMolecule2->GetNumberOfAtoms() + + fullMolecule3->GetNumberOfAtoms() - 2; + // the ghost bond is not duplicated in output. + nbOfExpectedBonds = fullMolecule1->GetNumberOfBonds() + fullMolecule2->GetNumberOfBonds() + + fullMolecule3->GetNumberOfBonds() - 1; + expectedResultValues->Resize(nbOfExpectedAtoms); + expectedResultValues->InsertValue( + 4, fullMolecule3->GetAtomData()->GetArray("Data")->GetTuple1(0)); + expectedResultValues->InsertValue( + 4, fullMolecule3->GetAtomData()->GetArray("Data")->GetTuple1(1)); + + return CheckMolecule(resultFullMolecule, + nbOfExpectedAtoms, + nbOfExpectedBonds, + nbOfExpectedArrays, + expectedResultValues, + 0, + 0); +} diff --git a/Filters/Core/vtkMoleculeAppend.cxx b/Filters/Core/vtkMoleculeAppend.cxx new file mode 100644 index 00000000000..8893c878ee0 --- /dev/null +++ b/Filters/Core/vtkMoleculeAppend.cxx @@ -0,0 +1,282 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkMoleculeAppend.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + + =========================================================================*/ +#include "vtkMoleculeAppend.h" + +#include "vtkAlgorithmOutput.h" +#include "vtkDataArray.h" +#include "vtkDataSetAttributes.h" +#include "vtkExecutive.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" +#include "vtkMergePoints.h" +#include "vtkMolecule.h" +#include "vtkNew.h" +#include "vtkObjectFactory.h" +#include "vtkPoints.h" +#include "vtkUnsignedCharArray.h" + +#include <set> +#include <utility> + +vtkStandardNewMacro(vtkMoleculeAppend); + +//---------------------------------------------------------------------------- +vtkMoleculeAppend::vtkMoleculeAppend() : MergeCoincidentAtoms(true) +{ +} + +//---------------------------------------------------------------------------- +vtkDataObject* vtkMoleculeAppend::GetInput(int idx) +{ + if (this->GetNumberOfInputConnections(0) <= idx) + { + return nullptr; + } + return vtkMolecule::SafeDownCast(this->GetExecutive()->GetInputData(0, idx)); +} + +//---------------------------------------------------------------------------- +int vtkMoleculeAppend::RequestData(vtkInformation*, + vtkInformationVector** inputVector, + vtkInformationVector* outputVector) +{ + vtkMolecule* output = vtkMolecule::GetData(outputVector, 0); + vtkDataSetAttributes* outputAtomData = output->GetAtomData(); + vtkDataSetAttributes* outputBondData = output->GetBondData(); + + // ******************** + // Create output data arrays following first input arrays. + vtkMolecule* mol0 = vtkMolecule::SafeDownCast(this->GetInput()); + outputAtomData->CopyStructure(mol0->GetAtomData()); + outputBondData->CopyStructure(mol0->GetBondData()); + output->SetAtomicNumberArrayName(mol0->GetAtomicNumberArrayName()); + output->SetBondOrdersArrayName(mol0->GetBondOrdersArrayName()); + vtkUnsignedCharArray* outputGhostAtoms = output->GetAtomGhostArray(); + vtkUnsignedCharArray* outputGhostBonds = output->GetBondGhostArray(); + + // ******************** + // Initialize unique atoms/bonds containers + vtkNew<vtkMergePoints> uniquePoints; + vtkNew<vtkPoints> uniquePointsList; + double bounds[6] = { 0., 0., 0., 0., 0., 0. }; + uniquePoints->InitPointInsertion(uniquePointsList, bounds, 0); + std::set<std::pair<vtkIdType, vtkIdType> > uniqueBonds; + + // ******************** + // Process each input + for (int idx = 0; idx < this->GetNumberOfInputConnections(0); ++idx) + { + vtkMolecule* input = vtkMolecule::GetData(inputVector[0], idx); + + // -------------------- + // Sanity check on input + int inputNbAtomArrays = input->GetAtomData()->GetNumberOfArrays(); + if (inputNbAtomArrays != outputAtomData->GetNumberOfArrays()) + { + vtkErrorMacro(<< "Input " << idx << ": Wrong number of atom array. Has " << inputNbAtomArrays + << " instead of " + << outputAtomData->GetNumberOfArrays()); + return 0; + } + + int inputNbBondArrays = input->GetBondData()->GetNumberOfArrays(); + if (input->GetNumberOfBonds() > 0 && inputNbBondArrays != outputBondData->GetNumberOfArrays()) + { + vtkErrorMacro(<< "Input " << idx << ": Wrong number of bond array. Has " << inputNbBondArrays + << " instead of " + << outputBondData->GetNumberOfArrays()); + return 0; + } + + for (vtkIdType ai = 0; ai < inputNbAtomArrays; ai++) + { + vtkDataArray* inArray = input->GetAtomData()->GetArray(ai); + if (!this->CheckArrays(inArray, outputAtomData->GetArray(inArray->GetName()))) + { + vtkErrorMacro(<< "Input " << idx << ": atoms arrays do not match with output"); + return 0; + } + } + + for (vtkIdType ai = 0; ai < inputNbBondArrays; ai++) + { + vtkDataArray* inArray = input->GetBondData()->GetArray(ai); + if (!this->CheckArrays(inArray, outputBondData->GetArray(inArray->GetName()))) + { + vtkErrorMacro(<< "Input " << idx << ": bonds arrays do not match with output"); + return 0; + } + } + + // -------------------- + // add atoms and bonds without duplication + + // map from 'input molecule atom ids' to 'output molecule atom ids' + std::vector<vtkIdType> atomIdMap(input->GetNumberOfAtoms(), -1); + + vtkIdType previousNbOfAtoms = output->GetNumberOfAtoms(); + int nbOfAtoms = 0; + for (vtkIdType i = 0; i < input->GetNumberOfAtoms(); i++) + { + double pt[3]; + input->GetAtomicPositionArray()->GetPoint(i, pt); + bool addAtom = true; + if (this->MergeCoincidentAtoms) + { + addAtom = uniquePoints->InsertUniquePoint(pt, atomIdMap[i]) == 1; + } + else + { + atomIdMap[i] = previousNbOfAtoms + nbOfAtoms; + } + + if (addAtom) + { + nbOfAtoms++; + vtkAtom atom = input->GetAtom(i); + output->AppendAtom(atom.GetAtomicNumber(), atom.GetPosition()).GetId(); + if (outputGhostAtoms) + { + outputGhostAtoms->InsertValue(atomIdMap[i], 255); + } + } + } + vtkIdType previousNbOfBonds = output->GetNumberOfBonds(); + int nbOfBonds = 0; + for (vtkIdType i = 0; i < input->GetNumberOfBonds(); i++) + { + vtkBond bond = input->GetBond(i); + // as bonds are undirected, put min atom number at first to avoid duplication. + vtkIdType atom1 = atomIdMap[bond.GetBeginAtomId()]; + vtkIdType atom2 = atomIdMap[bond.GetEndAtomId()]; + auto result = + uniqueBonds.insert(std::make_pair(std::min(atom1, atom2), std::max(atom1, atom2))); + if (result.second) + { + nbOfBonds++; + output->AppendBond(atom1, atom2, bond.GetOrder()); + } + } + + // -------------------- + // Reset arrays size (and allocation if needed) + for (vtkIdType ai = 0; ai < input->GetAtomData()->GetNumberOfArrays(); ai++) + { + vtkDataArray* inArray = input->GetAtomData()->GetArray(ai); + vtkDataArray* outArray = output->GetAtomData()->GetArray(inArray->GetName()); + outArray->Resize(previousNbOfAtoms + nbOfAtoms); + } + + for (vtkIdType ai = 0; ai < input->GetBondData()->GetNumberOfArrays(); ai++) + { + // skip bond orders array as it is auto-filled by AppendBond method + vtkDataArray* inArray = input->GetBondData()->GetArray(ai); + if (!strcmp(inArray->GetName(), input->GetBondOrdersArrayName())) + { + continue; + } + vtkDataArray* outArray = output->GetBondData()->GetArray(inArray->GetName()); + outArray->Resize(previousNbOfBonds + nbOfBonds); + } + + // -------------------- + // Fill DataArrays + for (vtkIdType i = 0; i < input->GetNumberOfAtoms(); i++) + { + for (vtkIdType ai = 0; ai < input->GetAtomData()->GetNumberOfArrays(); ai++) + { + vtkDataArray* inArray = input->GetAtomData()->GetArray(ai); + vtkDataArray* outArray = output->GetAtomData()->GetArray(inArray->GetName()); + // Use Value of non-ghost atom. + if (outputGhostAtoms && outputGhostAtoms->GetValue(atomIdMap[i]) == 0) + { + continue; + } + outArray->InsertTuple(atomIdMap[i], inArray->GetTuple(i)); + } + } + for (vtkIdType i = 0; i < input->GetNumberOfBonds(); i++) + { + vtkBond bond = input->GetBond(i); + vtkIdType outputBondId = + output->GetBondId(atomIdMap[bond.GetBeginAtomId()], atomIdMap[bond.GetBeginAtomId()]); + + for (vtkIdType ai = 0; ai < input->GetBondData()->GetNumberOfArrays(); ai++) + { + // skip bond orders array as it is auto-filled by AppendBond method + vtkDataArray* inArray = input->GetBondData()->GetArray(ai); + if (!strcmp(inArray->GetName(), input->GetBondOrdersArrayName())) + { + continue; + } + vtkDataArray* outArray = output->GetBondData()->GetArray(inArray->GetName()); + outArray->InsertTuple(outputBondId, inArray->GetTuple(i)); + } + } + } + + if (outputGhostBonds) + { + outputGhostBonds->SetNumberOfTuples(output->GetNumberOfBonds()); + outputGhostBonds->Fill(0); + for (vtkIdType bondId = 0; bondId < output->GetNumberOfBonds(); bondId++) + { + vtkIdType atom1 = output->GetBond(bondId).GetBeginAtomId(); + vtkIdType atom2 = output->GetBond(bondId).GetEndAtomId(); + if (outputGhostAtoms->GetValue(atom1) == 1 || outputGhostAtoms->GetValue(atom2) == 1) + { + outputGhostBonds->SetValue(bondId, 1); + } + } + } + + return 1; +} + +//---------------------------------------------------------------------------- +int vtkMoleculeAppend::FillInputPortInformation(int i, vtkInformation* info) +{ + info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1); + return this->Superclass::FillInputPortInformation(i, info); +} + +//---------------------------------------------------------------------------- +bool vtkMoleculeAppend::CheckArrays(vtkDataArray* array1, vtkDataArray* array2) +{ + if (strcmp(array1->GetName(), array2->GetName())) + { + vtkErrorMacro(<< "Execute: input name (" << array1->GetName() << "), must match output name (" + << array2->GetName() + << ")"); + return false; + } + + if (array1->GetDataType() != array2->GetDataType()) + { + vtkErrorMacro(<< "Execute: input ScalarType (" << array1->GetDataType() + << "), must match output ScalarType (" + << array2->GetDataType() + << ")"); + return false; + } + + if (array1->GetNumberOfComponents() != array2->GetNumberOfComponents()) + { + vtkErrorMacro("Components of the inputs do not match"); + return false; + } + + return true; +} diff --git a/Filters/Core/vtkMoleculeAppend.h b/Filters/Core/vtkMoleculeAppend.h new file mode 100644 index 00000000000..78695afb7ba --- /dev/null +++ b/Filters/Core/vtkMoleculeAppend.h @@ -0,0 +1,82 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkMoleculeAppend.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + + =========================================================================*/ + +/** + * @class vtkMoleculeAppend + * @brief Appends one or more molecules into a single molecule + * + * vtkMoleculeAppend appends molecule into a single molecule. It also appends + * the associated atom data and edge data. + * Note that input data arrays should match (same number of arrays with same names in each input) + * + * Option MergeCoincidentAtoms specifies if coincident atoms should be merged or not. + * This may be usefull in Parallel mode to remove ghost atoms when gather molecule on a rank. + * When merging, use the data of the non ghost atom. If none, use the data of the last coincident atom. + * This option is active by default. + */ + +#ifndef vtkMoleculeAppend_h +#define vtkMoleculeAppend_h + +#include "vtkFiltersCoreModule.h" // For export macro +#include "vtkMoleculeAlgorithm.h" + +class VTKFILTERSCORE_EXPORT vtkMoleculeAppend : public vtkMoleculeAlgorithm +{ +public: + static vtkMoleculeAppend* New(); + vtkTypeMacro(vtkMoleculeAppend, vtkMoleculeAlgorithm); + + //@{ + /** + * Get one input to this filter. This method is only for support of + * old-style pipeline connections. When writing new code you should + * use vtkAlgorithm::GetInputConnection(0, num). + */ + vtkDataObject* GetInput(int num); + vtkDataObject* GetInput() { return this->GetInput(0); } + //@} + + //@{ + /** + * Specify if coincident atoms (atom with excatly the same position) + * should be merged into one. + * True by default. + */ + vtkGetMacro(MergeCoincidentAtoms, bool); + vtkSetMacro(MergeCoincidentAtoms, bool); + vtkBooleanMacro(MergeCoincidentAtoms, bool); + // @} + +protected: + vtkMoleculeAppend(); + ~vtkMoleculeAppend() override = default; + + int RequestData(vtkInformation*, vtkInformationVector**, vtkInformationVector*) override; + + // see vtkAlgorithm for docs. + int FillInputPortInformation(int, vtkInformation*) override; + + // Check arrays information : name, type and number of components. + bool CheckArrays(vtkDataArray* array1, vtkDataArray* array2); + + bool MergeCoincidentAtoms; + +private: + vtkMoleculeAppend(const vtkMoleculeAppend&) = delete; + void operator=(const vtkMoleculeAppend&) = delete; +}; + +#endif -- GitLab