Commit 1b84a44d authored by Yohann Bearzi's avatar Yohann Bearzi Committed by Kitware Robot

Merge topic 'cell-collision'

6a81f6dd New overlapping cell detector filter
73c7971b New cell intersection API in `vtkCell`
e999e0ac Adding API to `vtkBoundingBox`
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: Utkarsh Ayachit's avatarUtkarsh Ayachit <utkarsh.ayachit@kitware.com>
Acked-by: T.J. Corona's avatarT.J. Corona <tj.corona@kitware.com>
Merge-request: !6789
parents 0bd205ba 6a81f6dd
......@@ -628,6 +628,38 @@ bool vtkBoundingBox::IntersectPlane(double origin[3], double normal[3])
}
//------------------------------------------------------------------------------
bool vtkBoundingBox::IntersectsSphere(double center[3], double radius) const
{
return center[0] >= this->MinPnt[0] - radius && center[0] <= this->MaxPnt[0] + radius &&
center[1] >= this->MinPnt[1] - radius && center[1] <= this->MaxPnt[1] + radius &&
center[2] >= this->MinPnt[2] - radius && center[2] <= this->MaxPnt[2] + radius;
}
// ---------------------------------------------------------------------------
int vtkBoundingBox::ComputeInnerDimension() const
{
double thickness = this->MaxPnt[0] - this->MinPnt[0];
int dim = 3;
if (std::abs(thickness) <=
std::max(std::fabs(this->MaxPnt[0]), std::fabs(this->MinPnt[0])) * VTK_DBL_EPSILON)
{
--dim;
}
thickness = this->MaxPnt[1] - this->MinPnt[1];
if (std::abs(thickness) <=
std::max(std::fabs(this->MaxPnt[1]), std::fabs(this->MinPnt[1])) * VTK_DBL_EPSILON)
{
--dim;
}
thickness = this->MaxPnt[2] - this->MinPnt[2];
if (std::fabs(thickness) <=
std::max(std::fabs(this->MaxPnt[2]), std::fabs(this->MinPnt[2])) * VTK_DBL_EPSILON)
{
--dim;
}
return dim;
}
// Support ComputeBounds()
namespace
{
......
......@@ -143,6 +143,11 @@ public:
*/
void AddBounds(const double bounds[]);
/**
* Returns true if this instance is entirely contained by bbox.
*/
bool IsSubsetOf(const vtkBoundingBox& bbox) const;
/**
* Intersect this box with bbox. The method returns 1 if both boxes are
* valid and they do have overlap else it will return 0. If 0 is returned
......@@ -162,6 +167,17 @@ public:
*/
bool IntersectPlane(double origin[3], double normal[3]);
/**
* Intersect this box with a sphere.
* Parameters involve the center of the sphere and the squared radius.
*/
bool IntersectsSphere(double center[3], double squaredRadius) const;
/**
* Returns the inner dimension of the bounding box.
*/
int ComputeInnerDimension() const;
/**
* Returns 1 if the min and max points of bbox are contained
* within the bounds of the specified box, else returns 0.
......@@ -376,6 +392,15 @@ inline void vtkBoundingBox::GetCenter(double center[3]) const
center[2] = 0.5 * (this->MaxPnt[2] + this->MinPnt[2]);
}
inline bool vtkBoundingBox::IsSubsetOf(const vtkBoundingBox& bbox) const
{
const double* bboxMaxPnt = bbox.GetMaxPoint();
const double* bboxMinPnt = bbox.GetMinPoint();
return this->MaxPnt[0] < bboxMaxPnt[0] && this->MinPnt[0] > bboxMinPnt[0] &&
this->MaxPnt[1] < bboxMaxPnt[1] && this->MinPnt[1] > bboxMinPnt[1] &&
this->MaxPnt[2] < bboxMaxPnt[2] && this->MinPnt[2] > bboxMinPnt[2];
}
inline void vtkBoundingBox::SetBounds(const double bounds[6])
{
this->SetBounds(bounds[0], bounds[1], bounds[2], bounds[3], bounds[4], bounds[5]);
......
......@@ -14,8 +14,80 @@
=========================================================================*/
#include "vtkCell.h"
#include "vtkDataArrayRange.h"
#include "vtkMath.h"
#include "vtkNew.h"
#include "vtkPoints.h"
#include "vtkPolygon.h"
#include "vtkSmartPointer.h"
#include <vector>
namespace detail
{
//----------------------------------------------------------------------------
// Strategy:
//
// We throw all edges from one cell to an other and look if they intersect.
// In the case of a cell of one point, we just check if it lies inside
// the other cell.
int IntersectWithCell(vtkCell* self, vtkCell* other, double tol)
{
if (!other->GetNumberOfPoints() || !self->GetNumberOfPoints())
{
return 0;
}
double x[3], pcoords[3];
if (other->GetNumberOfPoints() == 1)
{
double closestPoint[3];
double* point = other->GetPoints()->GetPoint(0);
int subId;
double dist2, *weights = new double[self->GetNumberOfPoints()];
self->EvaluatePosition(point, closestPoint, subId, pcoords, dist2, weights);
delete[] weights;
return dist2 <= tol * tol;
}
if (self->GetNumberOfPoints() == 1)
{
double closestPoint[3];
double* point = self->GetPoints()->GetPoint(0);
int subId;
double dist2, *weights = new double[other->GetNumberOfPoints()];
other->EvaluatePosition(point, closestPoint, subId, pcoords, dist2, weights);
delete[] weights;
return dist2 <= tol * tol;
}
double p1[3], p2[3];
for (vtkIdType edgeId = 0; edgeId < self->GetNumberOfEdges(); ++edgeId)
{
double t;
int subId;
vtkCell* edge = self->GetEdge(edgeId);
vtkPoints* ends = edge->GetPoints();
ends->GetPoint(0, p1);
ends->GetPoint(1, p2);
if (other->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
return 1;
}
}
for (vtkIdType edgeId = 0; edgeId < other->GetNumberOfEdges(); ++edgeId)
{
double t;
int subId;
vtkCell* edge = other->GetEdge(edgeId);
vtkPoints* ends = edge->GetPoints();
ends->GetPoint(0, p1);
ends->GetPoint(1, p2);
if (self->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
return 1;
}
}
return 0;
}
}
//------------------------------------------------------------------------------
// Construct cell.
......@@ -89,6 +161,104 @@ void vtkCell::DeepCopy(vtkCell* c)
}
//------------------------------------------------------------------------------
void vtkCell::Inflate(double dist)
{
if (this->GetNumberOfFaces() != 0)
{
vtkWarningMacro(<< "Base version of vtkCell::Inflate only implements cell inflation"
<< " for linear non 3D cells. Class " << this->GetClassName()
<< " needs to overload this method. Ignoring this cell.");
return;
}
std::vector<double> buf(3 * this->Points->GetNumberOfPoints(), 0.0);
vtkIdType pointId = 0;
auto pointRange = vtk::DataArrayTupleRange<3>(this->Points->GetData());
using TupleRef = typename decltype(pointRange)::TupleReferenceType;
auto stlPoints = vtk::DataArrayTupleRange<3>(this->Points->GetData());
auto postPointIt = stlPoints.begin();
++postPointIt;
auto prePointIt = stlPoints.end();
--prePointIt;
double v[3];
for (auto pointIt = stlPoints.begin(); pointIt != stlPoints.end();
++pointIt, ++postPointIt, ++prePointIt, ++pointId)
{
if (postPointIt == stlPoints.end())
{
postPointIt = stlPoints.begin();
}
if (prePointIt == stlPoints.end())
{
prePointIt = stlPoints.begin();
}
double* vIt = v;
for (auto compIt = pointIt->begin(), postCompIt = postPointIt->begin();
compIt != pointIt->end(); ++compIt, ++postCompIt, ++vIt)
{
*vIt = *compIt - *postCompIt;
}
vtkMath::Normalize(v);
buf[pointId * 3] += v[0] * dist;
buf[pointId * 3 + 1] += v[1] * dist;
buf[pointId * 3 + 2] += v[2] * dist;
vIt = v;
for (auto compIt = pointIt->begin(), preCompIt = prePointIt->begin(); compIt != pointIt->end();
++compIt, ++preCompIt, ++vIt)
{
*vIt = *compIt - *preCompIt;
}
vtkMath::Normalize(v);
buf[pointId * 3] += v[0] * dist;
buf[pointId * 3 + 1] += v[1] * dist;
buf[pointId * 3 + 2] += v[2] * dist;
}
auto it = buf.begin();
using TupleRef = typename decltype(stlPoints)::TupleReferenceType;
using CompRef = typename decltype(stlPoints)::ComponentReferenceType;
for (TupleRef point : stlPoints)
{
for (CompRef comp : point)
{
comp += *(it++);
}
}
}
//----------------------------------------------------------------------------
int vtkCell::IntersectWithCell(vtkCell* other, const vtkBoundingBox& boundingBox,
const vtkBoundingBox& otherBoundingBox, double tol)
{
if (!boundingBox.Intersects(otherBoundingBox))
{
return 0;
}
/**
* Given the strategy of detail::IntersectWithCell,
* the intersection detection is likely to be speeded up
* if exchanging other given this condition.
* The implementation first throws edges from first cell
* to look if it intersects with second cell, then it checks
* the other way.
* Since when one intersection is found, algorithm stops,
* we'd rather check embedded bounding box's cell's edges first.
*/
if (otherBoundingBox.IsSubsetOf(boundingBox))
{
return detail::IntersectWithCell(other, this, tol);
}
return detail::IntersectWithCell(this, other, tol);
}
//----------------------------------------------------------------------------
int vtkCell::IntersectWithCell(vtkCell* other, double tol)
{
return this->IntersectWithCell(
other, vtkBoundingBox(this->GetBounds()), vtkBoundingBox(other->GetBounds()), tol);
}
//----------------------------------------------------------------------------
// Compute cell bounding box (xmin,xmax,ymin,ymax,zmin,zmax). Return pointer
// to array of six double values.
double* vtkCell::GetBounds()
......
......@@ -43,8 +43,9 @@
#include "vtkCommonDataModelModule.h" // For export macro
#include "vtkObject.h"
#include "vtkCellType.h" // Needed to define cell types
#include "vtkIdList.h" // Needed for inline methods
#include "vtkBoundingBox.h" // Needed for IntersectWithCell
#include "vtkCellType.h" // Needed to define cell types
#include "vtkIdList.h" // Needed for inline methods
class vtkCellArray;
class vtkCellData;
......@@ -239,6 +240,13 @@ public:
vtkCellArray* connectivity, vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd,
vtkIdType cellId, vtkCellData* outCd, int insideOut) = 0;
/**
* Inflates the cell by dist. Each point is displaced by dist in the direction of every halfedge.
* Input can be negative. If so, the cell will shrink. Artifacts will appear if one shrinks
* with -dist larger than any edge length, probably causing the cell to self-intersect.
*/
virtual void Inflate(double dist);
/**
* Intersect with a ray. Return parametric coordinates (both line and cell)
* and global intersection coordinates, given ray definition p1[3], p2[3] and tolerance tol.
......@@ -250,6 +258,19 @@ public:
virtual int IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
double x[3], double pcoords[3], int& subId) = 0;
//@{
/**
* Intersects with an other cell. Returns 1 if cells intersect, 0 otherwise.
* If an exact intersection detection with regards to floating point precision is wanted,
* tol should be disregarded.
* `vtkBoundingBox` are optional parameters which slightly improve the speed of the computation
* if bounding boxes are already available to user.
*/
virtual int IntersectWithCell(vtkCell* other, double tol = 0.0);
virtual int IntersectWithCell(vtkCell* other, const vtkBoundingBox& boudingBox,
const vtkBoundingBox& otherBoundingBox, double tol = 0.0);
//@}
/**
* Generate simplices of proper dimension. If cell is 3D, tetrahedron are
* generated; if 2D triangles; if 1D lines; if 0D points. The form of the
......
......@@ -16,6 +16,7 @@
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkDoubleArray.h"
#include "vtkMath.h"
#include "vtkOrderedTriangulator.h"
......@@ -25,6 +26,8 @@
#include "vtkPolygon.h"
#include "vtkTetra.h"
#include <vector>
vtkCell3D::vtkCell3D()
{
this->Triangulator = nullptr;
......@@ -74,6 +77,48 @@ bool vtkCell3D::IsInsideOut()
return signedDistanceToCentroid > 0.0;
}
//----------------------------------------------------------------------------
void vtkCell3D::Inflate(double dist)
{
// The strategy is the following:
// - For each point, get points from its one-ring
// - Displace each of theses points following its halfedge direction,
// by a distance of dist.
std::vector<double> buf(3 * this->Points->GetNumberOfPoints(), 0.0);
double sign = this->IsInsideOut() ? -1.0 : 1.0;
double p[3];
vtkIdType pointId = 0;
auto pointRange = vtk::DataArrayTupleRange<3>(this->Points->GetData());
using ConstTupleRef = typename decltype(pointRange)::ConstTupleReferenceType;
for (ConstTupleRef point : pointRange)
{
const vtkIdType* ringIds;
vtkIdType ringSize = this->GetPointToOneRingPoints(pointId, ringIds);
for (vtkIdType id = 0; id < ringSize; ++id)
{
auto it = point->begin();
this->Points->GetPoint(ringIds[id], p);
double v[3] = { p[0] - *(it++), p[1] - *(it++), p[2] - *(it++) };
vtkMath::Normalize(v);
buf[3 * pointId] -= sign * v[0] * dist;
buf[3 * pointId + 1] -= sign * v[1] * dist;
buf[3 * pointId + 2] -= sign * v[2] * dist;
}
++pointId;
}
auto it = buf.begin();
using TupleRef = typename decltype(pointRange)::TupleReferenceType;
using CompRef = typename decltype(pointRange)::ComponentReferenceType;
for (TupleRef point : pointRange)
{
for (CompRef comp : point)
{
comp += *(it++);
}
}
}
//----------------------------------------------------------------------------
void vtkCell3D::Contour(double value, vtkDataArray* cellScalars,
vtkIncrementalPointLocator* locator, vtkCellArray* verts, vtkCellArray* lines,
vtkCellArray* polys, vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId,
......
......@@ -177,6 +177,13 @@ public:
*/
int GetCellDimension() override { return 3; }
/**
* Inflates the cell by dist. Each point is displaced by dist in the direction of every halfedge.
* Input can be negative. If so, the cell will shrink. Artifacts will appear if one shrinks
* with -dist larger than any edge length, probably causing the cell to self-intersect.
*/
void Inflate(double dist) override;
//@{
/**
* Set the tolerance for merging clip intersection points that are near
......
......@@ -16,6 +16,7 @@
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkLine.h"
#include "vtkMarchingSquaresLineCases.h"
......@@ -27,6 +28,9 @@
#include "vtkQuad.h"
#include "vtkTriangle.h"
#include <algorithm>
#include <array>
vtkStandardNewMacro(vtkPixel);
//------------------------------------------------------------------------------
......@@ -158,6 +162,82 @@ void vtkPixel::EvaluateLocation(int& subId, const double pcoords[3], double x[3]
}
//------------------------------------------------------------------------------
int vtkPixel::ComputeNormal(double n[3])
{
double p0[3], p1[3], p2[3];
this->Points->GetPoint(0, p0);
this->Points->GetPoint(1, p1);
this->Points->GetPoint(2, p2);
p1[0] -= p0[0];
p1[1] -= p0[1];
p1[2] -= p0[2];
p2[0] -= p0[0];
p2[1] -= p0[1];
p2[2] -= p0[2];
vtkMath::Cross(p1, p2, n);
if (std::abs(n[0]) < VTK_DBL_EPSILON && std::abs(n[1]) < VTK_DBL_EPSILON &&
std::abs(n[2]) < VTK_DBL_EPSILON)
{
return -1;
}
vtkMath::Normalize(n);
return (std::abs(n[1]) > 0.5) + (std::abs(n[2]) > 0.5) * 2;
}
//----------------------------------------------------------------------------
void vtkPixel::Inflate(double dist)
{
auto range = vtk::DataArrayTupleRange<3>(this->Points->GetData());
using TupleRef = typename decltype(range)::TupleReferenceType;
double n[3];
int normalDirection = this->ComputeNormal(n);
int degeneratePixelDirection = -1;
if (normalDirection == -1)
{
auto it = range.cbegin(), jt = range.cbegin();
auto compIt = it->cbegin(), compJt = jt->cbegin();
++jt;
std::array<double, 3> diff{ *(compJt++) - *(compIt++), *(compJt++) - *(compIt++),
*(compJt++) - *(compIt++) };
auto maxIt = std::max_element(diff.cbegin(), diff.cend(),
[](double a, double b) -> bool { return std::abs(a) < std::abs(b); });
degeneratePixelDirection = std::distance(diff.cbegin(), maxIt);
}
int index = 0;
for (TupleRef point : range)
{
auto it = point.begin();
switch (normalDirection)
{
case 0:
*(++it++) += dist * (index % 2 ? 1.0 : -1.0);
*(it++) += dist * (index / 2 ? 1.0 : -1.0);
break;
case 1:
*(it++) += dist * (index % 2 ? 1.0 : -1.0);
*(++it++) += dist * (index / 2 ? 1.0 : -1.0);
break;
case 2:
*(it++) += dist * (index % 2 ? 1.0 : -1.0);
*(it++) += dist * (index / 2 ? 1.0 : -1.0);
break;
case -1:
if (degeneratePixelDirection > 0)
{
++it;
if (degeneratePixelDirection > 1)
{
++it;
}
}
*it += dist * (index % 2 ? 1.0 : -1.0);
break;
}
++index;
}
}
//----------------------------------------------------------------------------
int vtkPixel::CellBoundary(int vtkNotUsed(subId), const double pcoords[3], vtkIdList* pts)
{
double t1 = pcoords[0] - pcoords[1];
......
......@@ -58,6 +58,7 @@ public:
int EvaluatePosition(const double x[3], double closestPoint[3], int& subId, double pcoords[3],
double& dist2, double weights[]) override;
void EvaluateLocation(int& subId, const double pcoords[3], double x[3], double* weights) override;
void Inflate(double dist) override;
//@}
/**
......@@ -89,6 +90,13 @@ public:
}
//@}
/**
* vtkPixel's normal cannot be computed using vtkPolygon::ComputeNormal because
* its points are not sorted such that circulating on them forms the pixel.
* This is a convenient method so one can compute normals on a pixel.
*/
int ComputeNormal(double n[3]);
protected:
vtkPixel();
~vtkPixel() override;
......
......@@ -17,6 +17,7 @@
#include "vtkBox.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkLine.h"
#include "vtkMarchingCubesTriangleCases.h"
......@@ -188,6 +189,22 @@ void vtkVoxel::EvaluateLocation(
}
//------------------------------------------------------------------------------
void vtkVoxel::Inflate(double dist)
{
int index = 0;
auto range = vtk::DataArrayTupleRange<3>(this->Points->GetData());
using TupleRef = typename decltype(range)::TupleReferenceType;
for (TupleRef point : range)
{
auto it = point.begin();
*(it++) += dist * (index % 2 ? 1.0 : -1.0);
*(it++) += dist * ((index / 2) % 2 ? 1.0 : -1.0);
*(it++) += dist * (index / 4 ? 1.0 : -1.0);
++index;
}
}
//----------------------------------------------------------------------------
//
// Compute Interpolation functions
//
......@@ -420,7 +437,7 @@ static constexpr vtkIdType pointToOneRingPoints[vtkVoxel::NumberOfPoints]
};
}
//------------------------------------------------------------------------------
//----------------------------------------------------------------------------
//
// Marching cubes case table
//
......
......@@ -114,6 +114,7 @@ public:
int Triangulate(int index, vtkIdList* ptIds, vtkPoints* pts) override;
void Derivatives(
int subId, const double pcoords[3], const double* values, int dim, double* derivs) override;
void Inflate(double dist) override;
//@}
static void InterpolationDerivs(const double pcoords[3], double derivs[24]);
......
......@@ -7,6 +7,7 @@ vtk_object_factory_declare(
set(classes
vtkAdaptiveResampleToImage
vtkOverlappingCellsDetector
vtkExtractSubsetWithSeed
vtkGenerateGlobalIds
vtkPResampleToImage
......
......@@ -2,6 +2,7 @@ vtk_module_test_data(
Data/multicomb_0.vts
Data/multicomb_1.vts
Data/multicomb_2.vts
Data/disk_out_ref.ex2)
Data/disk_out_ref.ex2
Data/overlapping_tetras.vtu)
add_subdirectory(Cxx)
......@@ -2,6 +2,7 @@ if (TARGET VTK::mpi)
set (vtkFiltersParallelDIY2CxxTests-MPI_NUMPROCS 2)
vtk_add_test_mpi(vtkFiltersParallelDIY2CxxTests-MPI tests
TESTING_DATA
TestOverlappingCellsDetector.cxx,NO_VALID
TestPResampleToImageCompositeDataSet.cxx
TestPResampleToImage.cxx
TestPResampleWithDataSet2.cxx
......@@ -28,8 +29,9 @@ endif()
# non-mpi tests
vtk_add_test_cxx(vtkFiltersParallelDIY2CxxTests non_mpi_tests
TestExtractSubsetWithSeed.cxx
TestAdaptiveResampleToImage.cxx,NO_VALID
TestExtractSubsetWithSeed.cxx
TestOverlappingCellsDetector.cxx,NO_VALID
TestGenerateGlobalIds.cxx,NO_VALID
TestRedistributeDataSetFilter.cxx,NO_VALID)
vtk_test_cxx_executable(vtkFiltersParallelDIY2CxxTests non_mpi_tests)
/*=========================================================================
Program: Visualization Toolkit
Module: TestOverlappingCellsDetector.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.
=========================================================================*/
#if VTK_MODULE_ENABLE_VTK_ParallelMPI
#include "vtkMPIController.h"
#else
#include "vtkDummyController.h"
#endif
#include "vtkCellData.h"
#include "vtkDataArray.h"
#include "vtkDataArrayRange.h"
#include "vtkDataSet.h"
#include "vtkGenerateGlobalIds.h"
#include "vtkMultiProcessController.h"
#include "vtkNew.h"
#include "vtkOverlappingCellsDetector.h"
#include "vtkRedistributeDataSetFilter.h"
#include "vtkTestUtilities.h"
#include "vtkUnstructuredGrid.h"
#include "vtkXMLUnstructuredGridReader.h"
namespace
{
static constexpr vtkIdType Collisions[72] = { 6, 0, 6, 0, 4, 4, 6, 0, 10, 7, // 0
4, 0, 0, 7, 9, 0, 0, 5, 5, 0, // 10
0, 0, 5, 9, 0, 6, 0, 6, 4, 4, // 20
0, 6, 1, 0, 4, 8, 7, 7, 1, 7, // 30
5, 0, 0, 5, 7, 5, 0, 2, 0, 0, // 40
0, 0, 0, 6, 1, 4, 0, 1, 0, 0, // 50
4, 0, 0, 0, 0, 0, 2, 6, 0, 0, // 60
0, 0 };
}
int TestOverlappingCellsDetector(int argc, char* argv[])
{
#if VTK_MODULE_ENABLE_VTK_ParallelMPI
vtkNew<vtkMPIController> contr;
#else
vtkNew<vtkDummyController> contr;
#endif
contr->Initialize(&argc, &argv);
vtkMultiProcessController::SetGlobalController(contr);
int retVal = EXIT_SUCCESS;
int myrank = contr->GetLocalProcessId();
const char* name =
vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/overlapping_tetras.vtu");
vtkNew<vtkGenerateGlobalIds> globalIds;
if (myrank == 0)
{
vtkNew<vtkXMLUnstructuredGridReader> reader;
reader->SetFileName(name);
globalIds->SetInputConnection(reader->GetOutputPort());
globalIds->Update();
}
else
{
vtkNew<vtkUnstructuredGrid> ug;
ug->Initialize();
globalIds->SetInputDataObject(ug);
}
vtkNew<vtkRedistributeDataSetFilter> redistribute;
redistribute->SetInputConnection(globalIds->GetOutputPort());
vtkNew<vtkOverlappingCellsDetector> detector;
detector->SetInputConnection(redistribute->GetOutputPort());
detector->Update();
vtkDataSet* output = vtkDataSet::SafeDownCast(detector->GetOutput(0));
vtkDataArray* data =
output->GetCellData()->GetArray(detector->GetNumberOfCollisionsPerCellArrayName());
vtkDataArray* ids = output->GetCellData()->GetArray("GlobalCellIds");