diff --git a/examples/xml/LoadDataArrays.scn b/examples/xml/LoadDataArrays.scn new file mode 100644 index 0000000..6a9572f --- /dev/null +++ b/examples/xml/LoadDataArrays.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + diff --git a/src/sofa/vtk/BaseVTKLoader.cpp b/src/sofa/vtk/BaseVTKLoader.cpp index 89dff35..009c5f4 100644 --- a/src/sofa/vtk/BaseVTKLoader.cpp +++ b/src/sofa/vtk/BaseVTKLoader.cpp @@ -1,5 +1,11 @@ #include +#include #include +#include +#include +#include +#include +#include #include #include @@ -7,7 +13,7 @@ namespace { template -vtkSmartPointer getDataSet(std::string fileName) +vtkSmartPointer getDataSet(const std::string& fileName) { vtkNew reader; reader->SetFileName(fileName.c_str()); @@ -15,11 +21,162 @@ vtkSmartPointer getDataSet(std::string fileName) return reader->GetOutput(); } +// Map VTK array value types to SOFA-registered Data types. +// float/double → SReal; small signed integers → int; unsigned integers → sofa::Index; +// long long / unsigned long long are kept as-is (both registered). +// long/unsigned long are platform-dependent and folded into the appropriate fixed-width type. +template struct SofaType { using type = T; }; +template<> struct SofaType { using type = SReal; }; +template<> struct SofaType { using type = SReal; }; +template<> struct SofaType { using type = int; }; +template<> struct SofaType { using type = int; }; +template<> struct SofaType { using type = sofa::Index; }; +template<> struct SofaType { using type = sofa::Index; }; +template<> struct SofaType { using type = sofa::Index; }; +template<> struct SofaType { using type = std::conditional_t; }; +template<> struct SofaType { using type = std::conditional_t; }; + +template +using SofaType_t = typename SofaType::type; + +struct ScalarDataWorker +{ + sofavtk::BaseVTKLoader& loader; + const std::string& arrayName; + std::unique_ptr result; + + template + void operator()(ArrayT* array) + { + using T = SofaType_t>; + + auto dataPtr = std::make_unique>>(); + dataPtr->setName(arrayName); + dataPtr->setHelp("Data array loaded from VTK file"); + + { + auto accessor = sofa::helper::getWriteOnlyAccessor(*dataPtr); + auto& vec = accessor.wref(); + vec.resize(array->GetNumberOfTuples()); + auto values = vtk::DataArrayValueRange<1>(array); + const vtkIdType n = static_cast(array->GetNumberOfTuples()); + for (vtkIdType i = 0; i < n; ++i) + vec[i] = values[i]; + } + + loader.addData(dataPtr.get(), arrayName); + result = std::move(dataPtr); + } +}; + +struct MultiComponentDataWorker +{ + sofavtk::BaseVTKLoader& loader; + const std::string& arrayName; + int numComponents; + std::unique_ptr result; + + template + void operator()(ArrayT* array) + { + using T = SofaType_t>; + dispatchN(array); + } + +private: + // Recursively dispatch fill function specialization based on the number of components until a + // match is found or the max supported is reached. + template + void dispatchN(ArrayT* array) + { + if (numComponents == N) + fill(array); + else if constexpr (N < 9) + dispatchN(array); + } + + template + void fill(ArrayT* array) + { + auto dataPtr = std::make_unique>>>(); + dataPtr->setName(arrayName); + dataPtr->setHelp("Data array loaded from VTK file"); + + { + auto accessor = sofa::helper::getWriteOnlyAccessor(*dataPtr); + auto& vec = accessor.wref(); + vec.resize(array->GetNumberOfTuples()); + auto tuples = vtk::DataArrayTupleRange(array); + const vtkIdType n = static_cast(array->GetNumberOfTuples()); + for (vtkIdType i = 0; i < n; ++i) + { + for (int c = 0; c < N; ++c) + vec[i][c] = tuples[i][c]; + } + } + + loader.addData(dataPtr.get(), arrayName); + result = std::move(dataPtr); + } +}; + } namespace sofavtk { +void BaseVTKLoader::loadDataArrayByName(vtkFieldData* fieldData, const std::string& arrayName, + std::map>& storage) +{ + vtkDataArray* array = fieldData->GetArray(arrayName.c_str()); + if (!array) + { + msg_warning() << fieldData->GetClassName() << " array '" << arrayName << "' not found in VTK file"; + return; + } + + const int numComponents = array->GetNumberOfComponents(); + + // Explicitly supported value types. Each maps to a fixed-size C++ type through vtkArrayDispatch + // long and unsigned long are included and remapped via CanonicalLong_t to int/long long based + // on their size on the current platform. + using SupportedTypes = vtkTypeList::Create< + float, double, + signed char, unsigned char, + int, unsigned int, + long, unsigned long, + long long, unsigned long long>; + using Dispatcher = vtkArrayDispatch::DispatchByValueType; + + auto dispatch = [&](auto& worker) { + if (!Dispatcher::Execute(array, worker)) + worker(array); + if (worker.result) + { + msg_info() << "Loaded " << fieldData->GetClassName() << " '" << arrayName + << "' with " << array->GetNumberOfTuples() << " entries"; + storage[arrayName] = std::move(worker.result); + } + }; + + if (numComponents == 1) + { + ScalarDataWorker worker{*this, arrayName}; + dispatch(worker); + return; + } + + if (numComponents > 9) + { + msg_warning() << fieldData->GetClassName() << " array '" << arrayName + << "' has " << numComponents << " components (max supported: 9)"; + return; + } + + MultiComponentDataWorker worker{*this, arrayName, numComponents}; + dispatch(worker); +} + bool BaseVTKLoader::doLoad() { const auto& fileName = d_filename.getFullPath(); @@ -36,13 +193,33 @@ bool BaseVTKLoader::doLoad() sofavtk::extractCells(*this, dataSet); + auto* cellData = dataSet->GetCellData(); + auto* pointData = dataSet->GetPointData(); + + for (const auto& arrayName : d_cellDataNames.getValue()) + loadDataArrayByName(cellData, arrayName, m_cellData); + + for (const auto& arrayName : d_pointDataNames.getValue()) + loadDataArrayByName(pointData, arrayName, m_pointData); + return true; } return false; } -void BaseVTKLoader::doClearBuffers() {} +void BaseVTKLoader::doClearBuffers() +{ + auto clearMap = [&](auto& map) { + for (const auto& [name, data] : map) + this->removeData(data.get()); + map.clear(); + }; + + clearMap(m_cellData); + clearMap(m_pointData); +} + } // namespace sofavtk diff --git a/src/sofa/vtk/BaseVTKLoader.h b/src/sofa/vtk/BaseVTKLoader.h index 28ad021..c5f4fa2 100644 --- a/src/sofa/vtk/BaseVTKLoader.h +++ b/src/sofa/vtk/BaseVTKLoader.h @@ -2,6 +2,7 @@ #include #include #include +#include #include namespace sofavtk @@ -11,11 +12,25 @@ struct SOFA_VTK_API BaseVTKLoader : sofa::core::loader::MeshLoader { SOFA_ABSTRACT_CLASS(BaseVTKLoader, sofa::core::loader::MeshLoader); + sofa::core::objectmodel::Data> d_cellDataNames{ + initData(&d_cellDataNames, "cellDataNames", + "Names of cell data arrays to load from the VTK file")}; + + sofa::core::objectmodel::Data> d_pointDataNames{ + initData(&d_pointDataNames, "pointDataNames", + "Names of point data arrays to load from the VTK file")}; + private: bool doLoad() final; void doClearBuffers() final; + void loadDataArrayByName(vtkFieldData* fieldData, const std::string& arrayName, + std::map>& storage); + + std::map> m_cellData; + std::map> m_pointData; + protected: virtual vtkSmartPointer getDataSet(const sofa::core::objectmodel::DataFileName& fileName) = 0; diff --git a/src/sofa/vtk/UnstructuredGridVTKLoader.cpp b/src/sofa/vtk/UnstructuredGridVTKLoader.cpp index 3c0c5bc..1707f88 100644 --- a/src/sofa/vtk/UnstructuredGridVTKLoader.cpp +++ b/src/sofa/vtk/UnstructuredGridVTKLoader.cpp @@ -1,9 +1,6 @@ #include #include -#include -#include - #include #include #include diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7696cb6..bc98c0c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,13 +1,20 @@ cmake_minimum_required(VERSION 3.12) project(SOFAVTK_test VERSION 1.0) find_package(Sofa.Testing REQUIRED) +find_package(VTK COMPONENTS CommonCore CommonDataModel REQUIRED) set(SOURCE_FILES test.cpp + TestCellPointData.cpp ) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) -target_link_libraries(${PROJECT_NAME} Sofa.Testing SOFA.VTK) +target_link_libraries(${PROJECT_NAME} Sofa.Testing SOFA.VTK ${VTK_LIBRARIES}) + +vtk_module_autoinit( + TARGETS ${PROJECT_NAME} + MODULES ${VTK_LIBRARIES} +) add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME}) diff --git a/tests/TestCellPointData.cpp b/tests/TestCellPointData.cpp new file mode 100644 index 0000000..e85d56b --- /dev/null +++ b/tests/TestCellPointData.cpp @@ -0,0 +1,264 @@ +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +// Fixture: 22 points, 10 cells — three spatial zones: +// Zone A z=0 2 quads (cells 0-1) + 2 triangles (cells 2-3) +// Zone B z=1→2 2 hexahedra (cells 4-5) +// Zone C z=2→3 4 tetrahedra (cells 6-9) +// +// All data values are functions of coordinates, verifiable by selecting +// an element in ParaView and computing the expected value directly: +// +// Cell data (10 entries): +// pressure float32 centroid z [0.0, 0.0, 0.0, 0.0, 1.5, 1.5, 2.25, 2.25, 2.25, 2.25] +// material_id int32 zone id [1, 1, 1, 1, 2, 2, 3, 3, 3, 3] +// fiber_direction float64 normalised centroid vector +// int8_data int8 cell index [1 .. 10] +// int32_data int32 index * 100 [100 .. 1000] +// int64_data int64 index * 10000 [10000 .. 100000] +// too_many float64 10 components, all 0 (>9-component guard) +// +// Point data (22 entries): +// temperature float64 x + y + z [0.0 .. 5.0] +// velocity float64 (x, y, z) radial from origin +static std::string fixtureVtuPath() +{ + std::filesystem::path here = std::filesystem::path(__FILE__).parent_path(); + return (here / "fixtures" / "test_mesh.vtu").string(); +} + +namespace +{ + +// Returns a pointer to the inner vector of a typed Data field on obj, or +// nullptr if the field is missing or has a different type. +template +const sofa::type::vector* getVec(sofa::core::objectmodel::Base* obj, + const std::string& name) +{ + auto* base = obj->findData(name); + if (!base) + return nullptr; + auto* typed = + dynamic_cast>*>(base); + return typed ? &typed->getValue() : nullptr; +} + +} // namespace + +class CellPointDataTest : public sofa::testing::BaseSimulationTest +{ +protected: + std::unique_ptr m_scene; + + sofavtk::UnstructuredGridVTKLoader* makeLoader( + const std::vector& cellNames = {}, + const std::vector& pointNames = {}) + { + m_scene = std::make_unique(); + auto loader = + sofa::core::objectmodel::New(); + m_scene->root->addObject(loader); + + loader->d_filename.setValue(fixtureVtuPath()); + loader->d_cellDataNames.setValue( + sofa::type::vector(cellNames.begin(), cellNames.end())); + loader->d_pointDataNames.setValue( + sofa::type::vector(pointNames.begin(), pointNames.end())); + + m_scene->initScene(); + return loader.get(); + } +}; + +// Fixture has 10 cells (2 quads + 2 tris + 2 hexas + 4 tets) +// and 22 points. + +// pressure = centroid z: 0.0 for surface cells, 1.5 for hexas, 2.25 for tets +TEST_F(CellPointDataTest, CellData_ScalarFloat) +{ + auto* loader = makeLoader({"pressure"}); + const auto* v = getVec(loader, "pressure"); + ASSERT_NE(v, nullptr); + ASSERT_EQ(v->size(), 10u); + // Zone A: centroid z = 0.0 + EXPECT_NEAR((*v)[0], 0.0, 1e-6); + EXPECT_NEAR((*v)[3], 0.0, 1e-6); + // Zone B: centroid z = 1.5 + EXPECT_NEAR((*v)[4], 1.5, 1e-6); + EXPECT_NEAR((*v)[5], 1.5, 1e-6); + // Zone C: centroid z = 2.25 + EXPECT_NEAR((*v)[6], 2.25, 1e-6); + EXPECT_NEAR((*v)[9], 2.25, 1e-6); +} + +// material_id = zone id: 1=surface, 2=hexa, 3=tet +TEST_F(CellPointDataTest, CellData_ScalarInt) +{ + auto* loader = makeLoader({"material_id"}); + const auto* v = getVec(loader, "material_id"); + ASSERT_NE(v, nullptr); + ASSERT_EQ(v->size(), 10u); + EXPECT_EQ((*v)[0], 1); // quad + EXPECT_EQ((*v)[2], 1); // triangle + EXPECT_EQ((*v)[4], 2); // hexa + EXPECT_EQ((*v)[6], 3); // tet +} + +// fiber_direction = normalised centroid vector +TEST_F(CellPointDataTest, CellData_Vec3Double) +{ + using Vec3r = sofa::type::Vec<3, SReal>; + auto* loader = makeLoader({"fiber_direction"}); + const auto* v = getVec(loader, "fiber_direction"); + ASSERT_NE(v, nullptr); + ASSERT_EQ(v->size(), 10u); + + // Every entry must be a unit vector + for (std::size_t i = 0; i < v->size(); ++i) + { + const SReal mag = (*v)[i].norm(); + EXPECT_NEAR(mag, 1.0, 1e-6) << "cell " << i << " fiber_direction is not unit"; + } + + // Cell 0 centroid = (0.5, 0.5, 0) → normalised (1/√2, 1/√2, 0) + const SReal inv_sqrt2 = SReal(1) / std::sqrt(SReal(2)); + EXPECT_NEAR((*v)[0][0], inv_sqrt2, 1e-6); + EXPECT_NEAR((*v)[0][1], inv_sqrt2, 1e-6); + EXPECT_NEAR((*v)[0][2], 0.0, 1e-6); + + // Zone C tets: centroid z >> x,y so z-component dominates + EXPECT_GT((*v)[6][2], SReal(0.9)); +} + +TEST_F(CellPointDataTest, CellData_Int8) +{ + auto* loader = makeLoader({"int8_data"}); + const auto* v = getVec(loader, "int8_data"); + ASSERT_NE(v, nullptr); + ASSERT_EQ(v->size(), 10u); + // int8_data = cell index + 1 + for (std::size_t i = 0; i < v->size(); ++i) + EXPECT_EQ(static_cast((*v)[i]), static_cast(i + 1)); +} + +TEST_F(CellPointDataTest, CellData_Int32) +{ + auto* loader = makeLoader({"int32_data"}); + const auto* v = getVec(loader, "int32_data"); + ASSERT_NE(v, nullptr); + ASSERT_EQ(v->size(), 10u); + // int32_data = (cell index + 1) * 100 + for (std::size_t i = 0; i < v->size(); ++i) + EXPECT_EQ((*v)[i], static_cast((i + 1) * 100)); +} + +TEST_F(CellPointDataTest, CellData_Int64) +{ + auto* loader = makeLoader({"int64_data"}); + const auto* v = getVec(loader, "int64_data"); + ASSERT_NE(v, nullptr); + ASSERT_EQ(v->size(), 10u); + // int64_data = (cell index + 1) * 10000 + for (std::size_t i = 0; i < v->size(); ++i) + EXPECT_EQ((*v)[i], static_cast((i + 1) * 10000)); +} + +// temperature = x + y + z per point, range 0.0 (P0 at origin) to 5.0 (P19) +TEST_F(CellPointDataTest, PointData_ScalarDouble) +{ + auto* loader = makeLoader({}, {"temperature"}); + const auto* v = getVec(loader, "temperature"); + ASSERT_NE(v, nullptr); + ASSERT_EQ(v->size(), 22u); + EXPECT_NEAR((*v)[0], 0.0, 1e-6); // P0 (0,0,0) + EXPECT_NEAR((*v)[1], 1.0, 1e-6); // P1 (1,0,0) + EXPECT_NEAR((*v)[5], 3.0, 1e-6); // P5 (2,1,0) + EXPECT_NEAR((*v)[19], 5.0, 1e-6); // P19 (2,1,2) + EXPECT_NEAR((*v)[21], 4.0, 1e-6); // P21 (0.5,0.5,3) +} + +// velocity = (x, y, z) per point — direction and magnitude encode position +TEST_F(CellPointDataTest, PointData_Vec3Double) +{ + using Vec3r = sofa::type::Vec<3, SReal>; + auto* loader = makeLoader({}, {"velocity"}); + const auto* v = getVec(loader, "velocity"); + ASSERT_NE(v, nullptr); + ASSERT_EQ(v->size(), 22u); + + // P0 (0,0,0) → zero velocity + EXPECT_NEAR((*v)[0][0], 0.0, 1e-6); + EXPECT_NEAR((*v)[0][1], 0.0, 1e-6); + EXPECT_NEAR((*v)[0][2], 0.0, 1e-6); + + // P5 (2,1,0) → velocity = coordinates + EXPECT_NEAR((*v)[5][0], 2.0, 1e-6); + EXPECT_NEAR((*v)[5][1], 1.0, 1e-6); + EXPECT_NEAR((*v)[5][2], 0.0, 1e-6); + + // P21 (0.5,0.5,3) → velocity = coordinates + EXPECT_NEAR((*v)[21][0], 0.5, 1e-6); + EXPECT_NEAR((*v)[21][1], 0.5, 1e-6); + EXPECT_NEAR((*v)[21][2], 3.0, 1e-6); +} + +TEST_F(CellPointDataTest, MissingArrayName) +{ + auto* loader = makeLoader({"nonexistent"}); + EXPECT_EQ(loader->findData("nonexistent"), nullptr); +} + +TEST_F(CellPointDataTest, TooManyComponents) +{ + auto* loader = makeLoader({"too_many"}); + EXPECT_EQ(loader->findData("too_many"), nullptr); +} + +TEST_F(CellPointDataTest, MultipleArraysAtOnce) +{ + auto* loader = makeLoader({"pressure", "material_id"}, {"temperature"}); + + const auto* pressure = getVec(loader, "pressure"); + ASSERT_NE(pressure, nullptr); + EXPECT_EQ(pressure->size(), 10u); + + const auto* matId = getVec(loader, "material_id"); + ASSERT_NE(matId, nullptr); + EXPECT_EQ(matId->size(), 10u); + + const auto* temp = getVec(loader, "temperature"); + ASSERT_NE(temp, nullptr); + EXPECT_EQ(temp->size(), 22u); +} + +TEST_F(CellPointDataTest, Reload) +{ + auto* loader = makeLoader({"pressure"}); + ASSERT_NE(getVec(loader, "pressure"), nullptr); + EXPECT_EQ(loader->findData("material_id"), nullptr); + + loader->d_cellDataNames.setValue({"material_id"}); + loader->load(); + + EXPECT_EQ(loader->findData("pressure"), nullptr); + const auto* matId = getVec(loader, "material_id"); + ASSERT_NE(matId, nullptr); + EXPECT_EQ(matId->size(), 10u); +} + +// temperature is point data; requesting it as cell data must not load it +TEST_F(CellPointDataTest, CellVsPoint_Separation) +{ + auto* loader = makeLoader({"temperature"}, {}); + EXPECT_EQ(loader->findData("temperature"), nullptr); +} diff --git a/tests/fixtures/test_mesh.vtu b/tests/fixtures/test_mesh.vtu new file mode 100644 index 0000000..701819a --- /dev/null +++ b/tests/fixtures/test_mesh.vtu @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + _AQAAAACAAACwAAAAIQAAAA==eJxjYEAGH+yhDAfsfA4H7Hxi1QkQEEfni2CIAwBGOgiaAQAAAACAAAAQAgAATAAAAA==eJxjYMAHPtjjkHDAqw1DH4yPYZ4DdvEHUP4PNPEfOMRx2Y/LHejuJ+ReXDS6+wmFFy7/Y4QnLvfhss8BuzyMDwvPBwTEORwAa3ElpA==AQAAAACAAAAoAAAAEwAAAA==eJxjYEAHB+whWMABGQMAMIIDPw==AQAAAACAAAAoAAAAEwAAAA==eJxjZGBgYETCTFDMjIYBAYAAFQ==AQAAAACAAADwAAAApgAAAA==eJw7Y12fNm/BM/szUJoBCqbl5gjPjX5nz+/pkSRsewUu7uWl5fXT9KY945q2a4sy38LFnRmvTol3eml/uEj8ZNrSx3BxtQxR22/el+1htKXMguiPH9/Yv+g6/0+O7Zm92U6RlXWxZ+F89+aIXW9nnYbSu+03TF5c6M333n7nQRH3JWcv239fI72lNf2U/TSl+VllW97C+TB5mDhMP8w8mDkAnZpoUw==AQAAAACAAAAKAAAAEgAAAA==eJxjZGJmYWVj5+DkAgAA5gA4AQAAAACAAAAoAAAAJwAAAA==eJxLYWBgOAHEOowMDBOA+AsQRzAxMOwBYgVmBoYWIH4BxAB0vAWNAQAAAACAAABQAAAALgAAAA==eJwTUGcAAwU/CG1QCqEd5kDogMMQOuEVhC4QZATTDRYQekI8hF7QBqEBJmYH1g==AQAAAACAAAAgAwAAEAAAAA==eJxjYBgFo2AU4AIAAyAAAQ==AQAAAACAAAAIAQAAQwAAAA==eJxtjgEKACAIA/ez9rSe5tMqKrkkQZyo5yRGb2is72zV3DP0rHF0QPOejMuvXOblV1/8nz7J88vIG2+fKtoeXtAXpA==AQAAAACAAABwAQAAUQAAAA==eJx1zjcSwCAMRFFwxgEnuP9VKfhqdkY0b5QJob+IE46SH3B2+ixeJLb+FTdMeOKBGW/88JW5HS/ZY3MP/rJH9xeszl2t63+0bvns1BsnxwIuAQAAAACAAABQAAAAIwAAAA==eJxjYYAADijNDaX5oLQYlJaD0kpQWg1Ka0FpPSgNAB5QAPo=AQAAAACAAAAKAAAAEAAAAA==eJzj5GRl5eHhAgIAAegAXQ== + +