diff --git a/vtkm/cont/testing/TestingCellLocatorRectilinearGrid.h b/vtkm/cont/testing/TestingCellLocatorRectilinearGrid.h index 985225010..e82501595 100644 --- a/vtkm/cont/testing/TestingCellLocatorRectilinearGrid.h +++ b/vtkm/cont/testing/TestingCellLocatorRectilinearGrid.h @@ -59,26 +59,35 @@ public: // Linear search in the coordinates. vtkm::Id index; /*Get floor X location*/ - for (index = 0; index < this->Dims[0] - 1; index++) - if (xAxis.Get(index) <= point[0] && point[0] < xAxis.Get(index + 1)) - { - logical[0] = index; - break; - } + if (point[0] == xAxis.Get(this->Dims[0] - 1)) + logical[0] = this->Dims[0] - 1; + else + for (index = 0; index < this->Dims[0] - 1; index++) + if (xAxis.Get(index) <= point[0] && point[0] < xAxis.Get(index + 1)) + { + logical[0] = index; + break; + } /*Get floor Y location*/ - for (index = 0; index < this->Dims[1] - 1; index++) - if (yAxis.Get(index) <= point[1] && point[1] < yAxis.Get(index + 1)) - { - logical[1] = index; - break; - } + if (point[1] == yAxis.Get(this->Dims[1] - 1)) + logical[1] = this->Dims[1] - 1; + else + for (index = 0; index < this->Dims[1] - 1; index++) + if (yAxis.Get(index) <= point[1] && point[1] < yAxis.Get(index + 1)) + { + logical[1] = index; + break; + } /*Get floor Z location*/ - for (index = 0; index < this->Dims[2] - 1; index++) - if (zAxis.Get(index) <= point[2] && point[2] < zAxis.Get(index + 1)) - { - logical[2] = index; - break; - } + if (point[2] == zAxis.Get(this->Dims[2] - 1)) + logical[2] = this->Dims[2] - 1; + else + for (index = 0; index < this->Dims[2] - 1; index++) + if (zAxis.Get(index) <= point[2] && point[2] < zAxis.Get(index + 1)) + { + logical[2] = index; + break; + } if (logical[0] == -1 || logical[1] == -1 || logical[2] == -1) return -1; return logical[2] * (Dims[0] - 1) * (Dims[1] - 1) + logical[1] * (Dims[0] - 1) + logical[0]; @@ -136,12 +145,8 @@ public: vtkm::cont::CoordinateSystem coords = dataset.GetCoordinateSystem(); vtkm::cont::DynamicCellSet cellSet = dataset.GetCellSet(); vtkm::Bounds bounds = coords.GetBounds(); - std::cout << "X bounds : " << bounds.X.Min << " to " << bounds.X.Max << std::endl; - std::cout << "Y bounds : " << bounds.Y.Min << " to " << bounds.Y.Max << std::endl; - std::cout << "Z bounds : " << bounds.Z.Min << " to " << bounds.Z.Max << std::endl; vtkm::Id3 dims = cellSet.Cast().GetSchedulingRange(vtkm::TopologyElementTagPoint()); - std::cout << "Dimensions of dataset : " << dims << std::endl; // Generate some sample points. using PointType = vtkm::Vec3f; @@ -180,7 +185,6 @@ public: { VTKM_TEST_ASSERT(matchPortal.Get(index), "Points do not match"); } - std::cout << "Test finished successfully." << std::endl; } void operator()() const diff --git a/vtkm/exec/CellLocatorRectilinearGrid.h b/vtkm/exec/CellLocatorRectilinearGrid.h index 1fabdfd8b..304169a7e 100644 --- a/vtkm/exec/CellLocatorRectilinearGrid.h +++ b/vtkm/exec/CellLocatorRectilinearGrid.h @@ -57,22 +57,23 @@ public: , Coords(coords.PrepareForInput(DeviceAdapter())) , PointDimensions(cellSet.GetPointDimensions()) { - auto contPortal = coords.GetPortalConstControl(); - for (vtkm::IdComponent dim = 0; dim < dimensions; dim++) + this->AxisPortals[0] = Coords.GetFirstPortal(); + this->MinPoint[0] = coords.GetPortalConstControl().GetFirstPortal().Get(0); + this->MaxPoint[0] = coords.GetPortalConstControl().GetFirstPortal().Get(PointDimensions[0] - 1); + + this->AxisPortals[1] = Coords.GetSecondPortal(); + this->MinPoint[1] = coords.GetPortalConstControl().GetSecondPortal().Get(0); + this->MaxPoint[1] = + coords.GetPortalConstControl().GetSecondPortal().Get(PointDimensions[1] - 1); + if (dimensions == 3) { - auto componentContPortal = (dim == 0) ? contPortal.GetFirstPortal() : (dim == 1) - ? contPortal.GetSecondPortal() - : contPortal.GetThirdPortal(); - auto componentExecPortal = (dim == 0) ? Coords.GetFirstPortal() : (dim == 1) - ? Coords.GetSecondPortal() - : Coords.GetThirdPortal(); - this->AxisPortals[dim] = componentExecPortal; - this->MinPoint[dim] = componentContPortal.Get(0); - this->MaxPoint[dim] = componentContPortal.Get(PointDimensions[0] - 1); + this->AxisPortals[2] = Coords.GetThirdPortal(); + this->MinPoint[2] = coords.GetPortalConstControl().GetThirdPortal().Get(0); + this->MaxPoint[2] = + coords.GetPortalConstControl().GetThirdPortal().Get(PointDimensions[2] - 1); } } - VTKM_EXEC_CONT virtual ~CellLocatorRectilinearGrid() noexcept { // This must not be defaulted, since defaulted virtual destructors are @@ -100,6 +101,7 @@ public: vtkm::Vec3f& parametric, const vtkm::exec::FunctorBase& worklet) const override { + (void)worklet; //suppress unused warning if (!IsInside(point)) { cellId = -1; @@ -121,47 +123,35 @@ public: continue; } - bool found = false; - vtkm::FloatDefault minVal = this->AxisPortals[dim].Get(logicalCell[dim]); - const vtkm::Id searchDir = (point[dim] - minVal >= 0.f) ? 1 : -1; - vtkm::FloatDefault maxVal = this->AxisPortals[dim].Get(logicalCell[dim] + 1); - - while (!found) + vtkm::Id minIndex = 0; + vtkm::Id maxIndex = PointDimensions[dim] - 1; + vtkm::FloatDefault minVal; + vtkm::FloatDefault maxVal; + minVal = this->AxisPortals[dim].Get(minIndex); + maxVal = this->AxisPortals[dim].Get(maxIndex); + while (maxIndex > minIndex + 1) { - if (point[dim] >= minVal && point[dim] < maxVal) + vtkm::Id midIndex = (minIndex + maxIndex) / 2; + vtkm::FloatDefault midVal = this->AxisPortals[dim].Get(midIndex); + if (point[dim] <= midVal) { - found = true; - continue; - } - - logicalCell[dim] += searchDir; - vtkm::Id nextCellId = searchDir == 1 ? logicalCell[dim] + 1 : logicalCell[dim]; - vtkm::FloatDefault next = this->AxisPortals[dim].Get(nextCellId); - if (searchDir == 1) - { - minVal = maxVal; - maxVal = next; + maxIndex = midIndex; + maxVal = midVal; } else { - maxVal = minVal; - minVal = next; + minIndex = midIndex; + minVal = midVal; } } + logicalCell[dim] = minIndex; + //printf("Min Index [%d] : %lld\n", dim, minIndex); + //printf("Max Index [%d] : %lld\n", dim, maxIndex); + //printf("Logical [%d] : %lld\n", dim, logicalCell[dim]); + parametric[dim] = (point[dim] - minVal) / (maxVal - minVal); } - // Get the actual cellId, from the logical cell index of the cell cellId = logicalCell[2] * this->PlaneSize + logicalCell[1] * this->RowSize + logicalCell[0]; - - bool success = false; - using IndicesType = typename CellSetPortal::IndicesType; - IndicesType cellPointIndices = this->CellSet.GetIndices(cellId); - vtkm::VecFromPortalPermute cellPoints(&cellPointIndices, - Coords); - auto cellShape = this->CellSet.GetCellShape(cellId); - // Get Parametric Coordinates from the cell, for the point. - parametric = vtkm::exec::WorldCoordinatesToParametricCoordinates( - cellPoints, point, cellShape, success, worklet); } private: