//============================================================================ // Copyright (c) Kitware, Inc. // All rights reserved. // See LICENSE.txt 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 #include #include #include #include #include namespace { static vtkm::cont::ArrayHandle StructuredGhostCellArray(vtkm::Id nx, vtkm::Id ny, vtkm::Id nz, int numLayers, bool addMidGhost = false) { vtkm::Id numCells = nx * ny; if (nz > 0) numCells *= nz; constexpr vtkm::UInt8 normalCell = vtkm::CellClassification::NORMAL; constexpr vtkm::UInt8 duplicateCell = vtkm::CellClassification::GHOST; vtkm::cont::ArrayHandle ghosts; ghosts.Allocate(numCells); auto portal = ghosts.WritePortal(); for (vtkm::Id i = 0; i < numCells; i++) { if (numLayers == 0) portal.Set(i, normalCell); else portal.Set(i, duplicateCell); } if (numLayers > 0) { //2D case if (nz == 0) { for (vtkm::Id i = numLayers; i < nx - numLayers; i++) for (vtkm::Id j = numLayers; j < ny - numLayers; j++) portal.Set(j * nx + i, normalCell); } else { for (vtkm::Id i = numLayers; i < nx - numLayers; i++) for (vtkm::Id j = numLayers; j < ny - numLayers; j++) for (vtkm::Id k = numLayers; k < nz - numLayers; k++) portal.Set(k * nx * ny + j * nx + i, normalCell); } } if (addMidGhost) { if (nz == 0) { vtkm::Id mi = numLayers + (nx - numLayers) / 2; vtkm::Id mj = numLayers + (ny - numLayers) / 2; portal.Set(mj * nx + mi, duplicateCell); } else { vtkm::Id mi = numLayers + (nx - numLayers) / 2; vtkm::Id mj = numLayers + (ny - numLayers) / 2; vtkm::Id mk = numLayers + (nz - numLayers) / 2; portal.Set(mk * nx * ny + mj * nx + mi, duplicateCell); } } return ghosts; } static vtkm::cont::DataSet MakeUniform(vtkm::Id numI, vtkm::Id numJ, vtkm::Id numK, int numLayers, bool addMidGhost = false) { vtkm::cont::DataSetBuilderUniform dsb; vtkm::cont::DataSet ds; if (numK == 0) ds = dsb.Create(vtkm::Id2(numI + 1, numJ + 1)); else ds = dsb.Create(vtkm::Id3(numI + 1, numJ + 1, numK + 1)); auto ghosts = StructuredGhostCellArray(numI, numJ, numK, numLayers, addMidGhost); ds.AddCellField("vtkmGhostCells", ghosts); return ds; } static vtkm::cont::DataSet MakeRectilinear(vtkm::Id numI, vtkm::Id numJ, vtkm::Id numK, int numLayers, bool addMidGhost = false) { vtkm::cont::DataSetBuilderRectilinear dsb; vtkm::cont::DataSet ds; std::size_t nx(static_cast(numI + 1)); std::size_t ny(static_cast(numJ + 1)); std::vector x(nx), y(ny); for (std::size_t i = 0; i < nx; i++) x[i] = static_cast(i); for (std::size_t i = 0; i < ny; i++) y[i] = static_cast(i); if (numK == 0) ds = dsb.Create(x, y); else { std::size_t nz(static_cast(numK + 1)); std::vector z(nz); for (std::size_t i = 0; i < nz; i++) z[i] = static_cast(i); ds = dsb.Create(x, y, z); } auto ghosts = StructuredGhostCellArray(numI, numJ, numK, numLayers, addMidGhost); ds.AddCellField("vtkmGhostCells", ghosts); return ds; } template static void MakeExplicitCells(const CellSetType& cellSet, vtkm::Vec& dims, vtkm::cont::ArrayHandle& numIndices, vtkm::cont::ArrayHandle& shapes, vtkm::cont::ArrayHandle& conn) { using Connectivity = vtkm::internal::ConnectivityStructuredInternals; vtkm::Id nCells = cellSet.GetNumberOfCells(); vtkm::Id connLen = (NDIM == 2 ? nCells * 4 : nCells * 8); conn.Allocate(connLen); shapes.Allocate(nCells); numIndices.Allocate(nCells); Connectivity structured; structured.SetPointDimensions(dims); vtkm::Id idx = 0; for (vtkm::Id i = 0; i < nCells; i++) { auto ptIds = structured.GetPointsOfCell(i); for (vtkm::IdComponent j = 0; j < NDIM; j++, idx++) conn.WritePortal().Set(idx, ptIds[j]); shapes.WritePortal().Set(i, (NDIM == 4 ? vtkm::CELL_SHAPE_QUAD : vtkm::CELL_SHAPE_HEXAHEDRON)); numIndices.WritePortal().Set(i, NDIM); } } static vtkm::cont::DataSet MakeExplicit(vtkm::Id numI, vtkm::Id numJ, vtkm::Id numK, int numLayers) { using CoordType = vtkm::Vec3f_32; vtkm::cont::DataSet dsUniform = MakeUniform(numI, numJ, numK, numLayers); auto coordData = dsUniform.GetCoordinateSystem(0).GetDataAsMultiplexer(); vtkm::Id numPts = coordData.GetNumberOfValues(); vtkm::cont::ArrayHandle explCoords; explCoords.Allocate(numPts); auto explPortal = explCoords.WritePortal(); auto cp = coordData.ReadPortal(); for (vtkm::Id i = 0; i < numPts; i++) explPortal.Set(i, cp.Get(i)); vtkm::cont::DynamicCellSet cellSet = dsUniform.GetCellSet(); vtkm::cont::ArrayHandle conn; vtkm::cont::ArrayHandle numIndices; vtkm::cont::ArrayHandle shapes; vtkm::cont::DataSet ds; vtkm::cont::DataSetBuilderExplicit dsb; if (cellSet.IsType>()) { vtkm::Id2 dims(numI, numJ); MakeExplicitCells( cellSet.Cast>(), dims, numIndices, shapes, conn); ds = dsb.Create(explCoords, vtkm::CellShapeTagQuad(), 4, conn, "coordinates"); } else if (cellSet.IsType>()) { vtkm::Id3 dims(numI, numJ, numK); MakeExplicitCells( cellSet.Cast>(), dims, numIndices, shapes, conn); ds = dsb.Create(explCoords, vtkm::CellShapeTagHexahedron(), 8, conn, "coordinates"); } auto ghosts = StructuredGhostCellArray(numI, numJ, numK, numLayers); ds.AddCellField("vtkmGhostCells", ghosts); return ds; } void TestGhostCellRemove() { // specify some 2d tests: {numI, numJ, numK, numGhostLayers}. std::vector> tests2D = { { 4, 4, 0, 2 }, { 5, 5, 0, 2 }, { 10, 10, 0, 3 }, { 10, 5, 0, 2 }, { 5, 10, 0, 2 }, { 20, 10, 0, 3 }, { 10, 20, 0, 3 } }; std::vector> tests3D = { { 4, 4, 4, 2 }, { 5, 5, 5, 2 }, { 10, 10, 10, 3 }, { 10, 5, 10, 2 }, { 5, 10, 10, 2 }, { 20, 10, 10, 3 }, { 10, 20, 10, 3 } }; std::vector> tests; tests.insert(tests.end(), tests2D.begin(), tests2D.end()); tests.insert(tests.end(), tests3D.begin(), tests3D.end()); for (auto& t : tests) { vtkm::Id nx = t[0], ny = t[1], nz = t[2]; int nghost = static_cast(t[3]); for (int layer = 0; layer < nghost; layer++) { std::vector dsTypes = { "uniform", "rectilinear", "explicit" }; for (auto& dsType : dsTypes) { vtkm::cont::DataSet ds; if (dsType == "uniform") ds = MakeUniform(nx, ny, nz, layer); else if (dsType == "rectilinear") ds = MakeRectilinear(nx, ny, nz, layer); else if (dsType == "explicit") ds = MakeExplicit(nx, ny, nz, layer); std::vector removeType = { "all", "byType" }; for (auto& rt : removeType) { vtkm::filter::GhostCellRemove ghostCellRemoval; ghostCellRemoval.RemoveGhostField(); if (rt == "all") ghostCellRemoval.RemoveAllGhost(); else if (rt == "byType") ghostCellRemoval.RemoveByType(vtkm::CellClassification::GHOST); auto output = ghostCellRemoval.Execute(ds); vtkm::Id numCells = output.GetNumberOfCells(); //Validate the output. vtkm::Id numCellsReq = (nx - 2 * layer) * (ny - 2 * layer); if (nz != 0) numCellsReq *= (nz - 2 * layer); VTKM_TEST_ASSERT(numCellsReq == numCells, "Wrong number of cells in output"); if (dsType == "uniform" || dsType == "rectilinear") { if (nz == 0) { VTKM_TEST_ASSERT(output.GetCellSet().IsSameType(vtkm::cont::CellSetStructured<2>()), "Wrong cell type for explicit conversion"); } else if (nz > 0) { VTKM_TEST_ASSERT(output.GetCellSet().IsSameType(vtkm::cont::CellSetStructured<3>()), "Wrong cell type for explicit conversion"); } } else { VTKM_TEST_ASSERT(output.GetCellSet().IsType>(), "Wrong cell type for explicit conversion"); } } // For structured, test the case where we have a ghost in the 'middle' of the cells. // This will produce an explicit cellset. if (dsType == "uniform" || dsType == "rectilinear") { if (dsType == "uniform") ds = MakeUniform(nx, ny, nz, layer, true); else if (dsType == "rectilinear") ds = MakeRectilinear(nx, ny, nz, layer, true); vtkm::filter::GhostCellRemove ghostCellRemoval; ghostCellRemoval.RemoveGhostField(); auto output = ghostCellRemoval.Execute(ds); VTKM_TEST_ASSERT(output.GetCellSet().IsType>(), "Wrong cell type for explicit conversion"); } } } } } } int UnitTestGhostCellRemove(int argc, char* argv[]) { return vtkm::cont::testing::Testing::Run(TestGhostCellRemove, argc, argv); }