From 31ed00991d03470b08be5fafe704fe8467ab759d Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 28 Aug 2026 13:08:32 -0500 Subject: [PATCH 01/14] VTU output includes the TimeValue field --- Code/Source/solver/VtkData.cpp | 42 +++++++++++++++++++++++++++++++--- Code/Source/solver/VtkData.h | 12 ++++++++++ Code/Source/solver/uris.cpp | 2 ++ Code/Source/solver/vtk_xml.cpp | 5 ++++ 4 files changed, 58 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/VtkData.cpp b/Code/Source/solver/VtkData.cpp index ba9fd3511..d941eb52d 100644 --- a/Code/Source/solver/VtkData.cpp +++ b/Code/Source/solver/VtkData.cpp @@ -5,9 +5,12 @@ #include "Array.h" #include "DebugMsg.h" -#include #include "vtkCellArray.h" #include "vtkCellData.h" +#include +#include +#include +#include #include #include #include @@ -19,13 +22,28 @@ #include #include #include -#include -#include ///////////////////////////////////////////////////////////////// // I n t e r n a l I m p l e m e n t a t i o n // ///////////////////////////////////////////////////////////////// +/// @brief Add a 'TimeValue' field data array to a VTK data object. +/// +/// The VTK XML readers use a single-tuple Float64 field data array named +/// 'TimeValue' as the pipeline time of the data object. An array of that name +/// already present in the field data is replaced. +/// +/// @param[in,out] data The data object being written. +/// @param[in] time The time value to associate with the data object. +static void set_time_value_field_data(vtkDataObject *data, const double time) { + auto time_array = vtkSmartPointer::New(); + time_array->SetName("TimeValue"); + time_array->SetNumberOfComponents(1); + time_array->SetNumberOfTuples(1); + time_array->SetValue(0, time); + data->GetFieldData()->AddArray(time_array); +} + ///////////////////////////////////////////////////////////////// // V t k V t p D a t a // ///////////////////////////////////////////////////////////////// @@ -39,6 +57,7 @@ class VtkVtpData::VtkVtpDataImpl { void set_connectivity(const int nsd, const Array& conn, const int pid); void set_point_data(const std::string& data_name, const Vector& data); void set_points(const Array& points); + void set_time_value(const double time); void write(const std::string& file_name); vtkSmartPointer vtk_polydata; @@ -246,6 +265,10 @@ void VtkVtpData::VtkVtpDataImpl::set_points(const Array& points) vtk_polydata->GetPointData()->AddArray(node_ids); } +void VtkVtpData::VtkVtpDataImpl::set_time_value(const double time) { + set_time_value_field_data(vtk_polydata, time); +} + void VtkVtpData::VtkVtpDataImpl::write(const std::string& file_name) { #ifdef debug_vtk_data @@ -276,6 +299,7 @@ class VtkVtuData::VtkVtuDataImpl { void set_point_data(const std::string& data_name, const Vector& data); void set_points(const Array& points); + void set_time_value(const double time); void write(const std::string& file_name); template @@ -526,6 +550,10 @@ void VtkVtuData::VtkVtuDataImpl::set_points(const Array& points) vtk_ugrid->SetPoints(node_coords); } +void VtkVtuData::VtkVtuDataImpl::set_time_value(const double time) { + set_time_value_field_data(vtk_ugrid, time); +} + void VtkVtuData::VtkVtuDataImpl::write(const std::string& file_name) { auto writer = vtkSmartPointer::New(); @@ -937,6 +965,10 @@ void VtkVtpData::set_points(const Array& points) impl->set_points(points); } +void VtkVtpData::set_time_value(const double time) { + impl->set_time_value(time); +} + void VtkVtpData::write() { impl->write(file_name); @@ -1287,6 +1319,10 @@ void VtkVtuData::set_points(const Array& points) impl->set_points(points); } +void VtkVtuData::set_time_value(const double time) { + impl->set_time_value(time); +} + void VtkVtuData::write() { impl->write(file_name); diff --git a/Code/Source/solver/VtkData.h b/Code/Source/solver/VtkData.h index 945a276f0..ca6c8fc89 100644 --- a/Code/Source/solver/VtkData.h +++ b/Code/Source/solver/VtkData.h @@ -33,6 +33,16 @@ class VtkData { virtual void set_points(const Array& points) = 0; virtual void set_connectivity(const int nsd, const Array& conn, const int pid = 0) = 0; + /// @brief Store a time value as field data, using the VTK convention for + /// time meta-data in XML files. + /// + /// The value is written as a single-tuple Float64 field data array named + /// 'TimeValue'. VTK XML readers such as ParaView turn this array into the + /// pipeline time of the data object. + /// + /// @param[in] time The time value to associate with the data. + virtual void set_time_value(const double time) = 0; + virtual bool has_cell_data(const std::string& data_name) = 0; virtual bool has_point_data(const std::string& data_name) = 0; @@ -103,6 +113,7 @@ class VtkVtpData : public VtkData { virtual void set_point_data(const std::string& data_name, const Vector& data) override; virtual void set_points(const Array& points) override; + virtual void set_time_value(const double time) override; virtual void write() override; private: @@ -151,6 +162,7 @@ class VtkVtuData : public VtkData { virtual void set_point_data(const std::string& data_name, const Vector& data) override; virtual void set_points(const Array& points) override; + virtual void set_time_value(const double time) override; virtual void write() override; private: diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 5c9926985..2ad9e6398 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -1107,6 +1107,8 @@ void uris_write_vtus(ComMod& com_mod) { std::string fName = com_mod.saveName + "_uris_" + uris_obj.name + "_" + fName_num + ".vtu"; auto vtk_writer = VtkData::create_writer(fName); + vtk_writer->set_time_value(com_mod.time); + // Writing the position data int iOut = 0; int s = outS(iOut); diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index c9c94f476..a638ee73c 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1467,6 +1467,11 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b fName = com_mod.saveName + "_" + fName + ".vtu"; auto vtk_writer = VtkData::create_writer(fName); + // No time field is assigned for time-averaged output. + if (!lAve) { + vtk_writer->set_time_value(com_mod.time); + } + // Writing the position data // int iOut = 0; From 93083d1de0d8d70b00dce4e9cabe0ecdc83fd417 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 28 Aug 2026 13:26:07 -0500 Subject: [PATCH 02/14] Export initial state to VTU file before starting iterations --- Code/Source/solver/main.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Code/Source/solver/main.cpp b/Code/Source/solver/main.cpp index 4cfb35fb4..e00a2ba49 100644 --- a/Code/Source/solver/main.cpp +++ b/Code/Source/solver/main.cpp @@ -266,6 +266,13 @@ void iterate_solution(Simulation* simulation) //Array::write_enabled = true; //Array3::write_enabled = true; + // Write the initial condition, which the Integrator has copied into the + // current solution. cTS is non-zero when restarting from a file or when + // continuing after remeshing, where the time step has already been written. + if (com_mod.saveVTK && cTS == 0 && com_mod.saveATS == 0) { + vtk_xml::write_vtus(simulation, solutions, /* lAvg = */ false); + } + // Outer loop for marching in time. When entering this loop, all old // variables are completely set and satisfy BCs. // From e9504764cf1a49634e6adbd4106295c098f8af53 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 28 Aug 2026 13:35:27 -0500 Subject: [PATCH 03/14] cep_init correctly assigns initial potential to old velocity solution vector --- Code/Source/solver/cep_ion.cpp | 16 ++++++++++------ Code/Source/solver/cep_ion.h | 14 +++++++++++++- Code/Source/solver/initialize.cpp | 2 +- 3 files changed, 24 insertions(+), 8 deletions(-) diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 479e77610..f9a92080c 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -11,12 +11,7 @@ namespace cep_ion { -/// @brief Modifies: -/// \code {.cpp} -/// cep_mod.Xion -/// \endcode -// -void cep_init(Simulation *simulation) { +void cep_init(Simulation *simulation, SolutionStates &solutions) { using namespace consts; auto &com_mod = simulation->com_mod; @@ -107,6 +102,15 @@ void cep_init(Simulation *simulation) { } } } + + // Copy the action potential into the solution vector, as cep_integ does at + // every time step. Without this the initial solution holds a zero potential + // while the ionic model rests at its own initial value. + auto &Yo = solutions.old.get_velocity(); + + for (int Ac = 0; Ac < tnNo; Ac++) { + Yo(eq.e, Ac) = cep_mod.Xion(0, Ac); + } } } diff --git a/Code/Source/solver/cep_ion.h b/Code/Source/solver/cep_ion.h index 754dec77f..8e15ec4db 100644 --- a/Code/Source/solver/cep_ion.h +++ b/Code/Source/solver/cep_ion.h @@ -16,7 +16,19 @@ namespace cep_ion { -void cep_init(Simulation *simulation); +/// @brief Initialize the ionic model state variables and the action potential +/// of the initial solution. +/// +/// The state variables of each CEP domain's ionic model are stored in +/// cep_mod.Xion. Its first row, the membrane potential, is also copied into the +/// row of the old velocity reserved for the CEP equation, so that the initial +/// solution agrees with the ionic model. +/// +/// @param[in,out] simulation The simulation, whose cep_mod.Xion is filled with +/// the initial ionic model state variables. +/// @param[in,out] solutions The solution states, whose old velocity receives +/// the initial action potential. +void cep_init(Simulation *simulation, SolutionStates &solutions); void cep_integ(Simulation *simulation, const int iEq, const int iDof, SolutionStates &solutions, const Vector &I4f); diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 7baf4f330..c797b912b 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -682,7 +682,7 @@ void initialize(Simulation* simulation, Vector& timeP) // if (cep_mod.cepEq) { cep_mod.Xion.resize(cep_mod.nXion,tnNo); - cep_ion::cep_init(simulation); + cep_ion::cep_init(simulation, initial_solutions); } // Electromechanics. From 061bd8a785fd3b42fbdbabfb9b240c14680769d4 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 28 Aug 2026 13:48:27 -0500 Subject: [PATCH 04/14] Add option to save Domain_ID field in every file, not just first output --- Code/Source/solver/ComMod.cpp | 50 +++++++++++++++---------------- Code/Source/solver/ComMod.h | 4 +++ Code/Source/solver/Parameters.cpp | 2 ++ Code/Source/solver/Parameters.h | 1 + Code/Source/solver/Simulation.cpp | 1 + Code/Source/solver/distribute.cpp | 1 + Code/Source/solver/vtk_xml.cpp | 8 +++-- 7 files changed, 39 insertions(+), 28 deletions(-) diff --git a/Code/Source/solver/ComMod.cpp b/Code/Source/solver/ComMod.cpp index 93d8d8be7..04dfc899a 100644 --- a/Code/Source/solver/ComMod.cpp +++ b/Code/Source/solver/ComMod.cpp @@ -11,31 +11,31 @@ // ComMod::ComMod() { - mvMsh = false; - - stFileFlag = false; - stFileRepl = false; - - bin2VTK = false; - saveAve = false; - sepOutput = false; - saveATS = 1; - saveIncr = 10; - nITs = 0; - startTS = 0; - stFileName = "stFile"; - iniFilePath = ""; - stopTrigName = "STOP_SIM"; - ichckIEN = true; - zeroAve = false; - cmmInit = false; - cmmVarWall = false; - shlEq = false; - pstEq = false; - sstEq = false; - ibFlag = false; - risFlag = false; - + mvMsh = false; + + stFileFlag = false; + stFileRepl = false; + + bin2VTK = false; + saveAve = false; + alwaysSaveDomainID = false; + sepOutput = false; + saveATS = 1; + saveIncr = 10; + nITs = 0; + startTS = 0; + stFileName = "stFile"; + iniFilePath = ""; + stopTrigName = "STOP_SIM"; + ichckIEN = true; + zeroAve = false; + cmmInit = false; + cmmVarWall = false; + shlEq = false; + pstEq = false; + sstEq = false; + ibFlag = false; + risFlag = false; } //--------- diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 2e9f60068..3538995b2 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1594,6 +1594,10 @@ class ComMod { /// @brief Whether to averaged results bool saveAve = false; + /// @brief Whether to save the domain ID to every VTK file rather than to + /// the first one only + bool alwaysSaveDomainID = false; + /// @brief Whether to save to VTK files bool saveVTK = false; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 8f9eed3d6..fdf713f4d 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2828,6 +2828,8 @@ GeneralSimulationParameters::GeneralSimulationParameters() { set_parameter("Save_averaged_results", false, !required, save_averaged_results); + set_parameter("Save_domain_ID_in_every_file", false, !required, + save_domain_id_in_every_file); set_parameter("Save_results_in_folder", "", !required, save_results_in_folder); set_parameter("Save_results_to_VTK_format", false, required, diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index b18bd4661..c77c5af54 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1819,6 +1819,7 @@ class GeneralSimulationParameters : public ParameterLists Parameter debug; Parameter overwrite_restart_file; Parameter save_averaged_results; + Parameter save_domain_id_in_every_file; Parameter save_results_to_vtk_format; Parameter simulation_requires_remeshing; Parameter start_averaging_from_zero; diff --git a/Code/Source/solver/Simulation.cpp b/Code/Source/solver/Simulation.cpp index 7e6d56523..1800c2d69 100644 --- a/Code/Source/solver/Simulation.cpp +++ b/Code/Source/solver/Simulation.cpp @@ -69,6 +69,7 @@ void Simulation::set_module_parameters() com_mod.saveIncr = general.increment_in_saving_vtk_files.value(); com_mod.saveATS = general.start_saving_after_time_step.value(); com_mod.saveAve = general.save_averaged_results.value(); + com_mod.alwaysSaveDomainID = general.save_domain_id_in_every_file.value(); com_mod.zeroAve = general.start_averaging_from_zero.value(); com_mod.stFileRepl = general.overwrite_restart_file.value(); com_mod.stFileName = chnl_mod.appPath + general.restart_file_name.value(); diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 00986613e..cabf6c774 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -320,6 +320,7 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &com_mod.saveATS); cm.bcast(cm_mod, &com_mod.saveAve); + cm.bcast(cm_mod, &com_mod.alwaysSaveDomainID); cm.bcast(cm_mod, &com_mod.saveVTK); cm.bcast(cm_mod, &com_mod.bin2VTK); diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index a638ee73c..790b1742a 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -183,7 +183,9 @@ void int_msh_data(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, // int m = nOute; - if (!com_mod.savedOnce || com_mod.nMsh > 1) { + // This condition selects the same xe columns as the one in write_vtus that + // reads them back. The two must agree. + if (!com_mod.savedOnce || com_mod.nMsh > 1 || com_mod.alwaysSaveDomainID) { if (com_mod.savedOnce) { m = m + 1; } else { @@ -230,7 +232,7 @@ void int_msh_data(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, // If files have not been written or there are multiple meshes. // - if (!com_mod.savedOnce || com_mod.nMsh > 1) { + if (!com_mod.savedOnce || com_mod.nMsh > 1 || com_mod.alwaysSaveDomainID) { #ifdef debug_int_msh_data dmsg << "!com_mod.savedOnce || com_mod.nMsh > 1 "; dmsg << "com_mod.dmnId.size(): " << com_mod.dmnId.size(); @@ -1538,7 +1540,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b // int ne = -1; - if (!com_mod.savedOnce || nMsh > 1) { + if (!com_mod.savedOnce || nMsh > 1 || com_mod.alwaysSaveDomainID) { Array tmpI(1,nEl); // Write the domain ID From 7d88503275a2d932f4d02779c0086cb63087ab9b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 28 Aug 2026 13:57:21 -0500 Subject: [PATCH 05/14] Force Increment_in_saving_VTK_files and Increment_in_saving_restart_files to be greater than 0 --- Code/Source/solver/Parameters.cpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index fdf713f4d..ef7ae02f6 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2807,9 +2807,9 @@ GeneralSimulationParameters::GeneralSimulationParameters() { set_parameter("Debug", false, !required, debug); set_parameter("Include_xml", "", !required, include_xml); - set_parameter("Increment_in_saving_restart_files", 0, !required, + set_parameter("Increment_in_saving_restart_files", 1, !required, increment_in_saving_restart_files); - set_parameter("Increment_in_saving_VTK_files", 0, !required, + set_parameter("Increment_in_saving_VTK_files", 1, !required, increment_in_saving_vtk_files); set_parameter("Name_prefix_of_saved_VTK_files", "", !required, @@ -2912,9 +2912,22 @@ void GeneralSimulationParameters::set_values(tinyxml2::XMLElement *xml_element, item = item->NextSiblingElement(); } - // Check that required parameters have been set. if (!from_external_xml) { + // Check that required parameters have been set. check_required(); + + // The saving increments select the time steps to save by taking the + // remainder of the time step number, so they cannot be zero. + svmp::check( + increment_in_saving_restart_files.value() >= 1, + "The GeneralSimulationParameters element " + "'Increment_in_saving_restart_files' must be greater than or equal " + "to 1."); + + svmp::check( + increment_in_saving_vtk_files.value() >= 1, + "The GeneralSimulationParameters element " + "'Increment_in_saving_VTK_files' must be greater than or equal to 1."); } } From e52d9aeb7960a2ad373d5558e5d1e99adbe6a2e1 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 28 Aug 2026 13:58:39 -0500 Subject: [PATCH 06/14] Remove unused ComMod::sepOutput --- Code/Source/solver/ComMod.cpp | 1 - Code/Source/solver/ComMod.h | 3 --- 2 files changed, 4 deletions(-) diff --git a/Code/Source/solver/ComMod.cpp b/Code/Source/solver/ComMod.cpp index 04dfc899a..d44249a7a 100644 --- a/Code/Source/solver/ComMod.cpp +++ b/Code/Source/solver/ComMod.cpp @@ -19,7 +19,6 @@ ComMod::ComMod() bin2VTK = false; saveAve = false; alwaysSaveDomainID = false; - sepOutput = false; saveATS = 1; saveIncr = 10; nITs = 0; diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 3538995b2..3fa0dd1ba 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1604,9 +1604,6 @@ class ComMod { /// @brief Whether any file being saved bool savedOnce = false; - /// @brief Whether to use separator in output - bool sepOutput = false; - /// @brief Whether start from beginning or from simulations bool stFileFlag = false; From 340cf04c8a8d86206ccd1ceff4632b70bf6cc747 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 1 Sep 2026 11:27:24 -0500 Subject: [PATCH 07/14] Fix typo in XML files --- tests/cases/stokes/manufactured_solution/P1P1/N008/solver.xml | 2 +- tests/cases/stokes/manufactured_solution/P2P1/N016/solver.xml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/cases/stokes/manufactured_solution/P1P1/N008/solver.xml b/tests/cases/stokes/manufactured_solution/P1P1/N008/solver.xml index 1776496e0..6d887a244 100644 --- a/tests/cases/stokes/manufactured_solution/P1P1/N008/solver.xml +++ b/tests/cases/stokes/manufactured_solution/P1P1/N008/solver.xml @@ -47,7 +47,7 @@ 1e-12 false - "Constant" > + 1.0 diff --git a/tests/cases/stokes/manufactured_solution/P2P1/N016/solver.xml b/tests/cases/stokes/manufactured_solution/P2P1/N016/solver.xml index f23b88d20..0ba63c71b 100644 --- a/tests/cases/stokes/manufactured_solution/P2P1/N016/solver.xml +++ b/tests/cases/stokes/manufactured_solution/P2P1/N016/solver.xml @@ -47,7 +47,7 @@ 1e-9 true - "Constant" > + 1.0 From 6788f74fb39467c2c8add7b17193fa0aef44c59e Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 1 Sep 2026 11:28:09 -0500 Subject: [PATCH 08/14] Integration tests also check the TimeValue field --- tests/conftest.py | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/tests/conftest.py b/tests/conftest.py index a0e8ce9be..e155f2ddd 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,10 +1,12 @@ import numpy as np +import math import pytest import os import shutil import subprocess import sys +import xml.etree.ElementTree as ET import meshio this_file_dir = os.path.abspath(os.path.dirname(__file__)) @@ -101,10 +103,35 @@ def _detect_oversubscribe_flag(): "Active_tension_normal": 1.0e-10, } +# Relative tolerance for the TimeValue field data. The solver accumulates the +# time as time += dt while the test computes t_max * dt, so the two differ only +# by floating point round-off. +RTOL_TIME_VALUE = 1.0e-12 + # Number of processors to test PROCS = [1, 3, 4] +def read_time_step_size(name_inp): + """ + Read the time step size from a svMultiPhysics input file + Args: + name_inp: path to the svMultiPhysics input file (.xml) + + Returns: + Time step size + """ + general = ET.parse(name_inp).getroot().find("GeneralSimulationParameters") + if general is None: + raise ValueError("No GeneralSimulationParameters in " + name_inp) + + time_step_size = general.find("Time_step_size") + if time_step_size is None: + raise ValueError("No Time_step_size in " + name_inp) + + return float(time_step_size.text) + + # Fixture to parametrize the number of processors for all tests @pytest.fixture(params=PROCS) def n_proc(request): @@ -178,6 +205,7 @@ def run_with_reference( t_max=1, name_ref=None, name_inp="solver.xml", + check_time_value=True, ): """ Run a test case and compare it to a stored reference solution @@ -188,6 +216,8 @@ def run_with_reference( t_max: time step to compare name_inp: name of svMultiPhysics input file (.xml) name_ref: name of refence file (.vtu) + check_time_value: whether to compare the TimeValue field data against + the time reached at time step t_max """ # default reference name if not name_ref: @@ -203,6 +233,22 @@ def run_with_reference( # check results msg = "" + + # check the time attached to the result as field data. This assumes a + # constant time step size, which does not hold if the case sets + # Number_of_initialization_time_steps. + if check_time_value: + if "TimeValue" not in res.field_data.keys(): + raise ValueError("Field data TimeValue not in simulation result") + + time_value = res.field_data["TimeValue"][0] + time_expected = t_max * read_time_step_size(os.path.join(folder, name_inp)) + + if not math.isclose(time_value, time_expected, rel_tol=RTOL_TIME_VALUE): + msg += "Test failed in field data TimeValue." + msg += " Result is " + str(time_value) + msg += " instead of " + str(time_expected) + ".\n" + for f in fields: # extract field if f not in res.point_data.keys(): From 1765aec5acd6874575853a325199cec9e0b78574 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 1 Sep 2026 17:10:32 -0500 Subject: [PATCH 09/14] Unify and harmonize classes VtkVtpData and VtkVtuData --- .../solver/CoupledBoundaryCondition.cpp | 2 +- Code/Source/solver/VtkData.cpp | 1437 +++++------------ Code/Source/solver/VtkData.h | 614 +++++-- Code/Source/solver/vtk_xml.cpp | 6 +- 4 files changed, 838 insertions(+), 1221 deletions(-) diff --git a/Code/Source/solver/CoupledBoundaryCondition.cpp b/Code/Source/solver/CoupledBoundaryCondition.cpp index 2d2ed353d..13518ae6d 100644 --- a/Code/Source/solver/CoupledBoundaryCondition.cpp +++ b/Code/Source/solver/CoupledBoundaryCondition.cpp @@ -993,7 +993,7 @@ void CappingSurface::load_from_vtp(const std::string& vtp_file_path, const faceT int vtk_cell_type = 0; try { conn = vtp_data.get_connectivity(); - eNoN = vtp_data.np_elem(); + eNoN = vtp_data.num_points_per_elem(); vtk_cell_type = vtp_data.elem_type(); } catch (const std::exception& e) { throw CappingSurfaceVtpException("Failed to get connectivity from VTP file '" + vtp_file_path + "': " + diff --git a/Code/Source/solver/VtkData.cpp b/Code/Source/solver/VtkData.cpp index d941eb52d..97d8fa97a 100644 --- a/Code/Source/solver/VtkData.cpp +++ b/Code/Source/solver/VtkData.cpp @@ -2,842 +2,411 @@ // SPDX-License-Identifier: BSD-3-Clause #include "VtkData.h" -#include "Array.h" -#include "DebugMsg.h" -#include "vtkCellArray.h" -#include "vtkCellData.h" -#include +#include +#include #include -#include + +#include +#include +#include #include +#include #include #include #include -#include -#include -#include -#include +#include #include #include #include #include -///////////////////////////////////////////////////////////////// -// I n t e r n a l I m p l e m e n t a t i o n // -///////////////////////////////////////////////////////////////// - -/// @brief Add a 'TimeValue' field data array to a VTK data object. -/// -/// The VTK XML readers use a single-tuple Float64 field data array named -/// 'TimeValue' as the pipeline time of the data object. An array of that name -/// already present in the field data is replaced. -/// -/// @param[in,out] data The data object being written. -/// @param[in] time The time value to associate with the data object. -static void set_time_value_field_data(vtkDataObject *data, const double time) { - auto time_array = vtkSmartPointer::New(); - time_array->SetName("TimeValue"); - time_array->SetNumberOfComponents(1); - time_array->SetNumberOfTuples(1); - time_array->SetValue(0, time); - data->GetFieldData()->AddArray(time_array); -} +void VtkData::read_file(const std::string &file_name) { + file_name_ = file_name; -///////////////////////////////////////////////////////////////// -// V t k V t p D a t a // -///////////////////////////////////////////////////////////////// - -class VtkVtpData::VtkVtpDataImpl { - public: - VtkVtpDataImpl(); - VtkVtpDataImpl(const VtkVtpDataImpl& other); - VtkVtpDataImpl& operator=(const VtkVtpDataImpl& other); - void read_file(const std::string& file_name); - void set_connectivity(const int nsd, const Array& conn, const int pid); - void set_point_data(const std::string& data_name, const Vector& data); - void set_points(const Array& points); - void set_time_value(const double time); - void write(const std::string& file_name); - - vtkSmartPointer vtk_polydata; - int elem_type; - int num_elems; - int np_elem; - int num_points; -}; - -VtkVtpData::VtkVtpDataImpl::VtkVtpDataImpl() -{ - vtk_polydata = vtkSmartPointer::New(); -} + read_file_internal(file_name); -VtkVtpData::VtkVtpDataImpl::VtkVtpDataImpl(const VtkVtpDataImpl& other) - : elem_type(other.elem_type) - , num_elems(other.num_elems) - , np_elem(other.np_elem) - , num_points(other.num_points) -{ - vtk_polydata = vtkSmartPointer::New(); - vtk_polydata->DeepCopy(other.vtk_polydata); -} + // Extract metadata (number of cells, points, cell type). + { + num_elems_ = vtk_data->GetNumberOfCells(); + num_points_ = vtk_data->GetPoints()->GetNumberOfPoints(); + if (num_points_ == 0) { + throw std::runtime_error("Error reading the VTK file '" + file_name + + "'."); + } -VtkVtpData::VtkVtpDataImpl& VtkVtpData::VtkVtpDataImpl::operator=(const VtkVtpDataImpl& other) -{ - if (this != &other) { - elem_type = other.elem_type; - num_elems = other.num_elems; - np_elem = other.np_elem; - num_points = other.num_points; - vtk_polydata = vtkSmartPointer::New(); - vtk_polydata->DeepCopy(other.vtk_polydata); + // Get the cell type. + auto cell = vtkGenericCell::New(); + vtk_data->GetCell(0, cell); + num_points_per_elem_ = cell->GetNumberOfPoints(); + elem_type_ = cell->GetCellType(); } - return *this; } -void VtkVtpData::VtkVtpDataImpl::read_file(const std::string& file_name) -{ - auto reader = vtkSmartPointer::New(); - reader->SetFileName(file_name.c_str()); - reader->Update(); - vtk_polydata= reader->GetOutput(); - num_elems = vtk_polydata->GetNumberOfCells(); - num_points = vtk_polydata->GetPoints()->GetNumberOfPoints(); - if (num_points == 0) { - throw std::runtime_error("Error reading the VTK VTP file '" + file_name + "'."); - } +Array VtkData::get_connectivity() const { + Array connectivity(num_points_per_elem_, num_elems_); - // Get the cell type. auto cell = vtkGenericCell::New(); - vtk_polydata->GetCell(0, cell); - np_elem = cell->GetNumberOfPoints(); - elem_type = cell->GetCellType(); -} - -void VtkVtpData::VtkVtpDataImpl::set_connectivity(const int nsd, const Array& conn, const int pid) -{ - #define n_debug_vtk_data - #ifdef debug_vtk_data - DebugMsg dmsg(__func__, 0); - dmsg << "[VtkVtpData.set_connectivity] vtk_polydata: " << vtk_polydata; - dmsg << "[VtkVtpData.set_connectivity] nsd: " << nsd; - #endif - - int num_elems = conn.ncols(); - int np_elem = conn.nrows(); - unsigned char vtk_cell_type; - - #ifdef debug_vtk_data - dmsg << "[VtkVtpData.set_connectivity] num_elems: " << num_elems; - dmsg << "[VtkVtpData.set_connectivity] np_elem: " << np_elem; - #endif - - if (nsd == 2) { - if (np_elem == 4) { - vtk_cell_type = VTK_QUAD; - } - else if (np_elem == 3) { - vtk_cell_type = VTK_TRIANGLE; - } - else if (np_elem == 6) { - vtk_cell_type = VTK_QUADRATIC_TRIANGLE; - } - else if (np_elem == 8) { - vtk_cell_type = VTK_QUADRATIC_QUAD; - } - else if (np_elem == 9) { - vtk_cell_type = VTK_BIQUADRATIC_QUAD; - } - - } else if (nsd == 3) { - if (np_elem == 4) { - vtk_cell_type = VTK_QUAD; - } - else if (np_elem == 3) { - vtk_cell_type = VTK_TRIANGLE; - #ifdef debug_vtk_data - dmsg << "[VtkVtpData.set_connectivity] vtk_cell_type = VTK_TRIANGLE"; - #endif - } - else if (np_elem == 6) { - vtk_cell_type = VTK_QUADRATIC_TRIANGLE; - } - else if (np_elem == 8) { - vtk_cell_type = VTK_HEXAHEDRON; - } - else if (np_elem == 10) { - vtk_cell_type = VTK_QUADRATIC_TETRA; - } - else if (np_elem == 20) { - vtk_cell_type = VTK_QUADRATIC_HEXAHEDRON; - } - else if (np_elem == 27) { - vtk_cell_type = VTK_TRIQUADRATIC_HEXAHEDRON; - } - } - - if (np_elem == 2) { - vtk_cell_type = VTK_LINE; - } - - auto elem_nodes = vtkSmartPointer::New(); - elem_nodes->Allocate(np_elem); - elem_nodes->Initialize(); - elem_nodes->SetNumberOfIds(np_elem); - - auto elem_ids = vtkSmartPointer::New(); - elem_ids->SetNumberOfComponents(1); - elem_ids->Allocate(num_elems,1000); - elem_ids->SetNumberOfTuples(num_elems); - elem_ids->SetName("GlobalElementID"); + for (int i = 0; i < num_elems_; i++) { + vtk_data->GetCell(i, cell); - vtkSmartPointer element_cells = vtkSmartPointer::New(); + const int num_cell_pts = cell->GetNumberOfPoints(); - for (int i = 0; i < num_elems; i++) { - #ifdef debug_vtk_data - dmsg << "[VtkVtpData.set_connectivity] ---------- i " << i; - #endif - for (int j = 0; j < np_elem; j++) { - #ifdef debug_vtk_data - dmsg << "[VtkVtpData.set_connectivity] ----- j " << j; - dmsg << "[VtkVtpData.set_connectivity] conn(j,i): " << conn(j,i); - #endif - elem_nodes->SetId(j, conn(j,i)); + for (int j = 0; j < num_cell_pts; ++j) { + connectivity(j, i) = cell->PointIds->GetId(j); } - element_cells->InsertNextCell(elem_nodes); - //vtk_polydata->InsertNextCell(vtk_cell_type, elem_nodes); - elem_ids->SetTuple1(i,i+1); - } - - vtk_polydata->SetPolys(element_cells); - vtk_polydata->GetCellData()->AddArray(elem_ids); -} - -void VtkVtpData::VtkVtpDataImpl::set_point_data(const std::string& data_name, const Vector& data) -{ - int num_vals = data.size(); - auto data_array = vtkSmartPointer::New(); - data_array->SetNumberOfComponents(1); - data_array->Allocate(num_vals); - data_array->SetName(data_name.c_str()); - - for (int i = 0; i < num_vals; i++) { - data_array->InsertNextTuple1(data(i)); } - vtk_polydata->GetPointData()->AddArray(data_array); + return connectivity; } -/// @brief Set the 3D point (coordinate) data for the polydata. -// -void VtkVtpData::VtkVtpDataImpl::set_points(const Array& points) -{ - #ifdef debug_vtk_data - DebugMsg dmsg(__func__, 0); - dmsg << "[VtkVtpData.set_points] vtk_polydata: " << vtk_polydata; - #endif - - int num_coords = points.ncols(); - if (num_coords == 0) { - throw std::runtime_error("Error in VTK VTP set_points: the number of points is zero."); - } - - #ifdef debug_vtk_data - dmsg << "[VtkVtpData.set_points] num_coords: " << num_coords; - #endif - - auto node_coords = vtkSmartPointer::New(); - node_coords->Allocate(num_coords, 1000); - node_coords->SetNumberOfPoints(num_coords); - - auto node_ids = vtkSmartPointer::New(); - node_ids->SetNumberOfComponents(1); - node_ids->Allocate(num_coords,1000); - node_ids->SetNumberOfTuples(num_coords); - node_ids->SetName("GlobalNodeID"); - - for (int i = 0; i < num_coords; i++ ) { - node_coords->SetPoint(i, points(0,i), points(1,i), points(2,i)); - node_ids->SetTuple1(i,i+1); - } - - vtk_polydata->SetPoints(node_coords); - vtk_polydata->GetPointData()->AddArray(node_ids); -} - -void VtkVtpData::VtkVtpDataImpl::set_time_value(const double time) { - set_time_value_field_data(vtk_polydata, time); -} - -void VtkVtpData::VtkVtpDataImpl::write(const std::string& file_name) -{ - #ifdef debug_vtk_data - DebugMsg dmsg(__func__, 0); - dmsg << "[VtkVtpData.write] file_name: " << file_name; - #endif - auto writer = vtkSmartPointer::New(); - writer->SetInputDataObject(vtk_polydata); - writer->SetFileName(file_name.c_str()); - writer->Write(); -} +Array VtkData::get_points() const { + auto vtk_points = vtk_data->GetPoints(); + auto num_points = vtk_points->GetNumberOfPoints(); -///////////////////////////////////////////////////////////////// -// V t k V t u D a t a // -///////////////////////////////////////////////////////////////// - -class VtkVtuData::VtkVtuDataImpl { - public: - void create_grid(); - void read_file(const std::string& file_name); - void set_connectivity(const int nsd, const Array& conn, const int pid); - - void set_element_data(const std::string& data_name, const Array& data); - void set_element_data(const std::string& data_name, const Array& data); - - void set_point_data(const std::string& data_name, const Array& data); - void set_point_data(const std::string& data_name, const Array& data); - void set_point_data(const std::string& data_name, const Vector& data); - - void set_points(const Array& points); - void set_time_value(const double time); - void write(const std::string& file_name); - - template - void set_element_data(const std::string& data_name, const T1& data, T2& data_array) - { - int num_vals = data.ncols(); - int num_comp = data.nrows(); - data_array->SetNumberOfComponents(num_comp); - data_array->Allocate(num_vals,1000); - data_array->SetNumberOfTuples(num_vals); - data_array->SetName(data_name.c_str()); - for (int i = 0; i < num_vals; i++) { - for (int j = 0; j < num_comp; j++) { - data_array->SetComponent(i, j, data(j,i)); - } - } - vtk_ugrid->GetCellData()->AddArray(data_array); - }; - - vtkSmartPointer vtk_ugrid; - int elem_type; - int num_elems; - int np_elem; - int num_points; -}; - -void VtkVtuData::VtkVtuDataImpl::create_grid() -{ - vtk_ugrid = vtkSmartPointer::New(); -} + Array points_array(3, num_points); -void VtkVtuData::VtkVtuDataImpl::read_file(const std::string& file_name) -{ - auto reader = vtkSmartPointer::New(); - reader->SetFileName(file_name.c_str()); - reader->Update(); - vtk_ugrid = reader->GetOutput(); - num_elems = vtk_ugrid->GetNumberOfCells(); - num_points = vtk_ugrid->GetPoints()->GetNumberOfPoints(); - if (num_points == 0) { - throw std::runtime_error("Error reading the VTK VTU file '" + file_name + "'."); + double point[3]; + for (int i = 0; i < num_points; i++) { + vtk_points->GetPoint(i, point); + points_array(0, i) = point[0]; + points_array(1, i) = point[1]; + points_array(2, i) = point[2]; } - // Get the cell type. - auto cell = vtkGenericCell::New(); - vtk_ugrid->GetCell(0, cell); - np_elem = cell->GetNumberOfPoints(); - elem_type = cell->GetCellType(); + return points_array; } -void VtkVtuData::VtkVtuDataImpl::set_connectivity(const int nsd, const Array& conn, const int pid) -{ - int num_elems = conn.ncols(); - int np_elem = conn.nrows(); - int num_coords = vtk_ugrid->GetPoints()->GetNumberOfPoints(); - unsigned char vtk_cell_type; - /* - std::cout << "[VtkVtuData.set_connectivity] " << std::endl; - std::cout << "[VtkVtuData.set_connectivity] nsd: " << nsd << std::endl; - std::cout << "[VtkVtuData.set_connectivity] num_elems: " << num_elems << std::endl; - std::cout << "[VtkVtuData.set_connectivity] np_elem: " << np_elem << std::endl; - std::cout << "[VtkVtuData.set_connectivity] num_coords: " << num_coords << std::endl; - */ - - if (nsd == 2) { - - if (np_elem == 3) { - vtk_cell_type = VTK_TRIANGLE; - - } else if (np_elem == 4) { - vtk_cell_type = VTK_QUAD; - - } else if (np_elem == 6) { - vtk_cell_type = VTK_QUADRATIC_TRIANGLE; - - } else if (np_elem == 8) { - vtk_cell_type = VTK_QUADRATIC_QUAD; - - } else if (np_elem == 9) { - vtk_cell_type = VTK_BIQUADRATIC_QUAD; - } - - } else if (nsd == 3) { - - if (np_elem == 3) { - vtk_cell_type = VTK_TRIANGLE; - - } else if (np_elem == 4) { - vtk_cell_type = VTK_TETRA; +int VtkData::num_elems() const { return num_elems_; } - } else if (np_elem == 6) { - vtk_cell_type = VTK_WEDGE; +int VtkData::elem_type() const { return elem_type_; } - } else if (np_elem == 8) { - vtk_cell_type = VTK_HEXAHEDRON; +int VtkData::num_points_per_elem() const { return num_points_per_elem_; } - } else if (np_elem == 10) { - vtk_cell_type = VTK_QUADRATIC_TETRA; +int VtkData::num_points() const { return num_points_; } - } else if (np_elem == 20) { - vtk_cell_type = VTK_QUADRATIC_HEXAHEDRON; +void VtkData::set_element_data(const std::string &data_name, + const Array &data) { + const int num_vals = data.ncols(); + const int num_components = data.nrows(); - } else if (np_elem == 27) { - vtk_cell_type = VTK_TRIQUADRATIC_HEXAHEDRON; + auto data_array = vtkSmartPointer::New(); + data_array->SetNumberOfComponents(num_components); + data_array->Allocate(num_vals, 1000); + data_array->SetNumberOfTuples(num_vals); + data_array->SetName(data_name.c_str()); + for (int i = 0; i < num_vals; ++i) { + for (int j = 0; j < num_components; ++j) { + data_array->SetComponent(i, j, data(j, i)); } - - } - - if (np_elem == 2) { - vtk_cell_type = VTK_LINE; } - auto elem_nodes = vtkSmartPointer::New(); - elem_nodes->Allocate(np_elem); - elem_nodes->Initialize(); - elem_nodes->SetNumberOfIds(np_elem); - //std::cout << "[VtkVtuData.set_connectivity] Set conn ... " << std::endl; - - for (int i = 0; i < num_elems; i++) { - //std::cout << "[VtkVtuData.set_connectivity] " << i << ": "; - for (int j = 0; j < np_elem; j++) { - //std::cout << conn(j,i) << " "; - if ((conn(j,i) < 0) || (conn(j,i) >= num_coords)) { - throw std::runtime_error("[VtkVtuData.set_connectivity] Element " + std::to_string(i+1) + - " has the non-valid node ID " + std::to_string(conn(j,i)) + "."); - } - elem_nodes->SetId(j, conn(j,i)); - } - //std::cout << std::endl; - vtk_ugrid->InsertNextCell(vtk_cell_type, elem_nodes); - } + vtk_data->GetCellData()->AddArray(data_array); } -//------------------ -// set_element_data -//------------------ -// -/* -void VtkVtuData::VtkVtuDataImpl::set_element_data(const std::string& data_name, const Array& data) -{ - int num_vals = data.num_cols(); - int num_comp = data.num_rows(); +void VtkData::set_element_data(const std::string &data_name, + const Array &data) { + const int num_vals = data.ncols(); + const int num_components = data.nrows(); - auto data_array = vtkSmartPointer::New(); - data_array->SetNumberOfComponents(num_comp); - data_array->Allocate(num_vals,1000); + auto data_array = vtkSmartPointer::New(); + data_array->SetNumberOfComponents(num_components); + data_array->Allocate(num_vals, 1000); data_array->SetNumberOfTuples(num_vals); data_array->SetName(data_name.c_str()); - - for (int i = 0; i < num_vals; i++) { - for (int j = 0; j < num_comp; j++) { - data_array->SetComponent(i, j, data(j,i)); + for (int i = 0; i < num_vals; ++i) { + for (int j = 0; j < num_components; ++j) { + data_array->SetComponent(i, j, data(j, i)); } } - vtk_ugrid->GetCellData()->AddArray(data_array); + vtk_data->GetCellData()->AddArray(data_array); } -//------------------ -// set_element_data -//------------------ -// -void VtkVtuData::VtkVtuDataImpl::set_element_data(const std::string& data_name, const Array& data) -{ - int num_vals = data.num_cols(); - int num_comp = data.num_rows(); +void VtkData::set_point_data(const std::string &data_name, + const Array &data) { + const int num_vals = data.ncols(); + const int num_comp = data.nrows(); - auto data_array = vtkSmartPointer::New(); + auto data_array = vtkSmartPointer::New(); data_array->SetNumberOfComponents(num_comp); - data_array->Allocate(num_vals,1000); + data_array->Allocate(num_vals, 1000); data_array->SetNumberOfTuples(num_vals); data_array->SetName(data_name.c_str()); for (int i = 0; i < num_vals; i++) { for (int j = 0; j < num_comp; j++) { - data_array->SetComponent(i, j, data(j,i)); + data_array->SetComponent(i, j, data(j, i)); } } - vtk_ugrid->GetCellData()->AddArray(data_array); + vtk_data->GetPointData()->AddArray(data_array); } -*/ - -//---------------- -// set_point_data -//---------------- -// -void VtkVtuData::VtkVtuDataImpl::set_point_data(const std::string& data_name, const Array& data) -{ + +void VtkData::set_point_data(const std::string &data_name, + const Array &data) { int num_vals = data.ncols(); int num_comp = data.nrows(); - auto data_array = vtkSmartPointer::New(); + auto data_array = vtkSmartPointer::New(); data_array->SetNumberOfComponents(num_comp); - data_array->Allocate(num_vals,1000); + data_array->Allocate(num_vals, 1000); data_array->SetNumberOfTuples(num_vals); data_array->SetName(data_name.c_str()); for (int i = 0; i < num_vals; i++) { for (int j = 0; j < num_comp; j++) { - data_array->SetComponent(i, j, data(j,i)); + data_array->SetComponent(i, j, data(j, i)); } } - vtk_ugrid->GetPointData()->AddArray(data_array); + vtk_data->GetPointData()->AddArray(data_array); } -void VtkVtuData::VtkVtuDataImpl::set_point_data(const std::string& data_name, const Array& data) -{ - int num_vals = data.ncols(); - int num_comp = data.nrows(); +void VtkData::set_point_data(const std::string &data_name, + const Vector &data) { + const int num_vals = data.size(); auto data_array = vtkSmartPointer::New(); - data_array->SetNumberOfComponents(num_comp); - data_array->Allocate(num_vals,1000); - data_array->SetNumberOfTuples(num_vals); + data_array->SetNumberOfComponents(1); + data_array->Allocate(num_vals); data_array->SetName(data_name.c_str()); for (int i = 0; i < num_vals; i++) { - for (int j = 0; j < num_comp; j++) { - data_array->SetComponent(i, j, data(j,i)); - } + data_array->InsertNextTuple1(data(i)); } - vtk_ugrid->GetPointData()->AddArray(data_array); + vtk_data->GetPointData()->AddArray(data_array); } -void VtkVtuData::VtkVtuDataImpl::set_point_data(const std::string& data_name, const Vector& data) -{ - throw std::runtime_error("[VtkVtuData] set_point_data for Vector not implemented."); -} +void VtkData::set_points(const Array &points) { + const int num_coords = points.ncols(); + if (num_coords == 0) { + throw std::runtime_error( + "Error in vtkData::set_points: the number of points is zero."); + } -/// @brief Set the 3D points (coordinates) data for the unstructured grid. -// -void VtkVtuData::VtkVtuDataImpl::set_points(const Array& points) -{ - int num_coords = points.ncols(); auto node_coords = vtkSmartPointer::New(); - node_coords->Allocate(num_coords ,1000); + node_coords->Allocate(num_coords, 1000); node_coords->SetNumberOfPoints(num_coords); - //std::cout << "[VtkVtuData.set_points] " << std::endl; - //std::cout << "[VtkVtuData.set_points] num_coords: " << num_coords << std::endl; - for (int i = 0; i < num_coords; i++ ) { - node_coords->SetPoint(i, points(0,i), points(1,i), points(2,i)); + for (int i = 0; i < num_coords; i++) { + node_coords->SetPoint(i, points(0, i), points(1, i), points(2, i)); } - vtk_ugrid->SetPoints(node_coords); -} - -void VtkVtuData::VtkVtuDataImpl::set_time_value(const double time) { - set_time_value_field_data(vtk_ugrid, time); -} - -void VtkVtuData::VtkVtuDataImpl::write(const std::string& file_name) -{ - auto writer = vtkSmartPointer::New(); - writer->SetInputDataObject(vtk_ugrid); - writer->SetFileName(file_name.c_str()); - writer->Write(); + vtk_data->SetPoints(node_coords); } -///////////////////////////////////////////////////////////////// -// V t k D a t a I m p l e m e n t a t i o n // -///////////////////////////////////////////////////////////////// - -VtkData::VtkData() -{ -} +void VtkData::set_connectivity(const int nsd, const Array &conn) { + int num_elems = conn.ncols(); + int np_elem = conn.nrows(); -VtkData::~VtkData() -{ -} + auto elem_nodes = vtkSmartPointer::New(); + elem_nodes->Allocate(np_elem); + elem_nodes->Initialize(); + elem_nodes->SetNumberOfIds(np_elem); -VtkData* VtkData::create_reader(const std::string& file_name) -{ - auto file_ext = file_name.substr(file_name.find_last_of(".") + 1); - if (file_ext == "vtp") { - return new VtkVtpData(file_name); - } else if (file_ext == "vtu") { - return new VtkVtuData(file_name); - } -} + for (int i = 0; i < num_elems; i++) { + for (int j = 0; j < np_elem; j++) { + elem_nodes->SetId(j, conn(j, i)); + } -VtkData* VtkData::create_writer(const std::string& file_name) -{ - auto file_ext = file_name.substr(file_name.find_last_of(".") + 1); - bool reader = false; - if (file_ext == "vtp") { - return new VtkVtpData(file_name, reader); - } else if (file_ext == "vtu") { - return new VtkVtuData(file_name, reader); + insert_cell(cell_type(nsd, np_elem), elem_nodes); } } -//void VtkData::write(const std::string& file_name) -//{ -//} +void VtkData::set_time_value(const double time) { + auto time_array = vtkSmartPointer::New(); -///////////////////////////////////////////////////////////////// -// V t k V t p D a t a I m p l e m e n t a t i o n // -///////////////////////////////////////////////////////////////// + time_array->SetName("TimeValue"); + time_array->SetNumberOfComponents(1); + time_array->SetNumberOfTuples(1); + time_array->SetValue(0, time); -VtkVtpData::VtkVtpData() -{ - impl = new VtkVtpDataImpl; + vtk_data->GetFieldData()->AddArray(time_array); } -VtkVtpData::VtkVtpData(const std::string& file_name, bool reader) -{ - this->file_name = file_name; - impl = new VtkVtpDataImpl; - if (reader) { - read_file(file_name); - } -} +bool VtkData::has_cell_data(const std::string &data_name) const { + const int num_arrays = vtk_data->GetCellData()->GetNumberOfArrays(); -VtkVtpData::~VtkVtpData() -{ - delete impl; -} + for (int i = 0; i < num_arrays; i++) { + if (!strcmp(vtk_data->GetCellData()->GetArrayName(i), data_name.c_str())) { + return true; + } + } -VtkVtpData::VtkVtpData(const VtkVtpData& other) -{ - impl = new VtkVtpDataImpl(*other.impl); - file_name = other.file_name; + return false; } -VtkVtpData& VtkVtpData::operator=(const VtkVtpData& other) -{ - if (this != &other) { - delete impl; - impl = new VtkVtpDataImpl(*other.impl); - file_name = other.file_name; +bool VtkData::has_point_data(const std::string &data_name) const { + const int num_arrays = vtk_data->GetPointData()->GetNumberOfArrays(); + + for (int i = 0; i < num_arrays; i++) { + if (!strcmp(vtk_data->GetPointData()->GetArrayName(i), data_name.c_str())) { + return true; + } } - return *this; + + return false; } -Array VtkVtpData::get_connectivity() const -{ - int num_elems = impl->num_elems; - int np_elem = impl->np_elem; - Array conn(np_elem, num_elems); +void VtkData::copy_points(Array &points) const { + auto vtk_points = vtk_data->GetPoints(); + auto num_points = vtk_points->GetNumberOfPoints(); - auto cell = vtkGenericCell::New(); - for (int i = 0; i < num_elems; i++) { - impl->vtk_polydata->GetCell(i, cell); - auto num_cell_pts = cell->GetNumberOfPoints(); - for (int j = 0; j < num_cell_pts; j++) { - auto id = cell->PointIds->GetId(j); - conn(j,i) = id; - } + double point[3]; + for (int i = 0; i < num_points; i++) { + vtk_points->GetPoint(i, point); + points(0, i) = point[0]; + points(1, i) = point[1]; + points(2, i) = point[2]; } - - return conn; } -/// @brief Copy an array of cell data from a polydata mesh into the given Array. -// -void VtkVtpData::copy_cell_data(const std::string& data_name, Array& mesh_data) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_polydata->GetCellData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { +void VtkData::copy_point_data(const std::string &data_name, + Array &mesh_data) const { + const auto vtk_array = vtkDoubleArray::SafeDownCast( + vtk_data->GetPointData()->GetArray(data_name.c_str())); + if (vtk_array == nullptr) { + // @todo[michelebucelli] This should probably be an exception. return; } - int num_data = vtk_data->GetNumberOfTuples(); + const int num_data = vtk_array->GetNumberOfTuples(); if (num_data == 0) { - return; - } + // @todo[michelebucelli] This should probably be an exception. + return; + } - int num_comp = vtk_data->GetNumberOfComponents(); + const int num_comp = vtk_array->GetNumberOfComponents(); // Set the data. for (int i = 0; i < num_data; i++) { - auto tuple = vtk_data->GetTuple(i); + const auto tuple = vtk_array->GetTuple(i); for (int j = 0; j < num_comp; j++) { mesh_data(j, i) = tuple[j]; } } } -void VtkVtpData::copy_cell_data(const std::string& data_name, Vector& mesh_data) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_polydata->GetCellData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { +void VtkData::copy_point_data(const std::string &data_name, + Vector &mesh_data) const { + const auto vtk_array = vtkDoubleArray::SafeDownCast( + vtk_data->GetPointData()->GetArray(data_name.c_str())); + if (vtk_array == nullptr) { return; } - int num_data = vtk_data->GetNumberOfTuples(); + int num_data = vtk_array->GetNumberOfTuples(); if (num_data == 0) { return; } // Set the data. for (int i = 0; i < num_data; i++) { - mesh_data(i) = vtk_data->GetValue(i); + mesh_data(i) = vtk_array->GetValue(i); } } -/// @brief Copy an array of cell data from a polydata mesh into the given Vector. -// -void VtkVtpData::copy_cell_data(const std::string& data_name, Vector& mesh_data) -{ - auto vtk_data = vtkIntArray::SafeDownCast(impl->vtk_polydata->GetCellData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { +void VtkData::copy_point_data(const std::string &data_name, + Vector &mesh_data) const { + const auto vtk_array = vtkIntArray::SafeDownCast( + vtk_data->GetPointData()->GetArray(data_name.c_str())); + if (vtk_array == nullptr) { return; } - int num_data = vtk_data->GetNumberOfTuples(); + int num_data = vtk_array->GetNumberOfTuples(); if (num_data == 0) { return; } - int num_comp = vtk_data->GetNumberOfComponents(); - + // Set the data. for (int i = 0; i < num_data; i++) { - auto tuple = vtk_data->GetTuple(i); - for (int j = 0; j < num_comp; j++) { - mesh_data(i) = tuple[j]; - } + mesh_data(i) = vtk_array->GetValue(i); } } -std::pair VtkVtpData::get_cell_data_dimensions(const std::string& data_name) const -{ - auto vtk_array = impl->vtk_polydata->GetCellData()->GetArray(data_name.c_str()); +void VtkData::copy_cell_data(const std::string &data_name, + Array &mesh_data) const { + const auto vtk_array = vtkDoubleArray::SafeDownCast( + vtk_data->GetCellData()->GetArray(data_name.c_str())); if (vtk_array == nullptr) { - return std::make_pair(0, 0); - } - - return std::make_pair(vtk_array->GetNumberOfComponents(), vtk_array->GetNumberOfTuples()); -} - -/// @brief Copy an array of point data from an polydata mesh into the given Array. -// -void VtkVtpData::copy_point_data(const std::string& data_name, Array& mesh_data) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_polydata->GetPointData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { return; } - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return; + const int num_data = vtk_array->GetNumberOfTuples(); + if (num_data == 0) { + return; } - int num_comp = vtk_data->GetNumberOfComponents(); + const int num_comp = vtk_array->GetNumberOfComponents(); // Set the data. for (int i = 0; i < num_data; i++) { - auto tuple = vtk_data->GetTuple(i); + const auto tuple = vtk_array->GetTuple(i); for (int j = 0; j < num_comp; j++) { mesh_data(j, i) = tuple[j]; } } } -void VtkVtpData::copy_point_data(const std::string& data_name, Vector& mesh_data) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_polydata->GetPointData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { +void VtkData::copy_cell_data(const std::string &data_name, + Vector &mesh_data) const { + const auto vtk_array = vtkDoubleArray::SafeDownCast( + vtk_data->GetCellData()->GetArray(data_name.c_str())); + if (vtk_array == nullptr) { return; } - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return; + const int num_data = vtk_array->GetNumberOfTuples(); + if (num_data == 0) { + return; } - int num_comp = vtk_data->GetNumberOfComponents(); - // Set the data. for (int i = 0; i < num_data; i++) { - mesh_data(i) = vtk_data->GetValue(i); + mesh_data(i) = vtk_array->GetValue(i); } } -void VtkVtpData::copy_point_data(const std::string& data_name, Vector& mesh_data) -{ - auto vtk_data = vtkIntArray::SafeDownCast(impl->vtk_polydata->GetPointData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { +void VtkData::copy_cell_data(const std::string &data_name, + Vector &mesh_data) const { + const auto vtk_array = vtkIntArray::SafeDownCast( + vtk_data->GetCellData()->GetArray(data_name.c_str())); + if (vtk_array == nullptr) { return; } - int num_data = vtk_data->GetNumberOfTuples(); + const int num_data = vtk_array->GetNumberOfTuples(); if (num_data == 0) { return; } - int num_comp = vtk_data->GetNumberOfComponents(); - // Set the data. for (int i = 0; i < num_data; i++) { - mesh_data(i) = vtk_data->GetValue(i); + mesh_data(i) = vtk_array->GetValue(i); } } -/// @brief Copy points into the given array. -// -void VtkVtpData::copy_points(Array& points) -{ - auto vtk_points = impl->vtk_polydata->GetPoints(); - auto num_points = vtk_points->GetNumberOfPoints(); - Array points_array(3, num_points); - - double point[3]; - for (int i = 0; i < num_points; i++) { - vtk_points->GetPoint(i, point); - points(0,i) = point[0]; - points(1,i) = point[1]; - points(2,i) = point[2]; - } - - return; -} - -/// @brief Get an array of point data from an unstructured grid. -// -Array VtkVtpData::get_point_data(const std::string& data_name) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_polydata->GetPointData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { +Array VtkData::get_point_data(const std::string &data_name) const { + auto vtk_array = vtkDoubleArray::SafeDownCast( + vtk_data->GetPointData()->GetArray(data_name.c_str())); + if (vtk_array == nullptr) { return Array(); } - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { + int num_data = vtk_array->GetNumberOfTuples(); + if (num_data == 0) { return Array(); } - int num_comp = vtk_data->GetNumberOfComponents(); + int num_comp = vtk_array->GetNumberOfComponents(); // Set the data. Array data(num_data, num_comp); for (int i = 0; i < num_data; i++) { - auto tuple = vtk_data->GetTuple(i); + auto tuple = vtk_array->GetTuple(i); for (int j = 0; j < num_comp; j++) { data(i, j) = tuple[j]; } @@ -846,485 +415,223 @@ Array VtkVtpData::get_point_data(const std::string& data_name) return data; } -/// @brief Get a list of point data names. -std::vector VtkVtpData::get_point_data_names() -{ - std::vector data_names; - int num_arrays = impl->vtk_polydata->GetPointData()->GetNumberOfArrays(); +std::vector VtkData::get_point_data_names() const { + std::vector data_names; + const int num_arrays = vtk_data->GetPointData()->GetNumberOfArrays(); for (int i = 0; i < num_arrays; i++) { - auto array_name = impl->vtk_polydata->GetPointData()->GetArrayName(i); - data_names.push_back(array_name); + data_names.push_back(vtk_data->GetPointData()->GetArrayName(i)); } - return data_names; + return data_names; } -/// @brief Get an array of point data from an unstructured grid. -// -Array VtkVtpData::get_points() const -{ - auto vtk_points = impl->vtk_polydata->GetPoints(); - auto num_points = vtk_points->GetNumberOfPoints(); - Array points_array(3, num_points); - - double point[3]; - for (int i = 0; i < num_points; i++) { - vtk_points->GetPoint(i, point); - points_array(0,i) = point[0]; - points_array(1,i) = point[1]; - points_array(2,i) = point[2]; +std::pair +VtkData::get_cell_data_dimensions(const std::string &data_name) const { + auto vtk_array = vtk_data->GetCellData()->GetArray(data_name.c_str()); + if (vtk_array == nullptr) { + return std::make_pair(0, 0); } - return points_array; + return std::make_pair(vtk_array->GetNumberOfComponents(), + vtk_array->GetNumberOfTuples()); } -bool VtkVtpData::has_cell_data(const std::string& data_name) -{ - int num_arrays = impl->vtk_polydata->GetCellData()->GetNumberOfArrays(); - - for (int i = 0; i < num_arrays; i++) { - if (!strcmp(impl->vtk_polydata->GetCellData()->GetArrayName(i), data_name.c_str())) { - return true; - } +VtkData *VtkData::create_reader(const std::string &file_name) { + auto file_ext = file_name.substr(file_name.find_last_of(".") + 1); + if (file_ext == "vtp") { + return new VtkVtpData(file_name); + } else if (file_ext == "vtu") { + return new VtkVtuData(file_name); } - return false; + throw std::runtime_error( + "Error in VtkData::create_reader: the file '" + file_name + + "' has the extension '" + file_ext + "', which is not 'vtp' or 'vtu'."); } -bool VtkVtpData::has_point_data(const std::string& data_name) -{ - int num_arrays = impl->vtk_polydata->GetPointData()->GetNumberOfArrays(); - - for (int i = 0; i < num_arrays; i++) { - if (!strcmp(impl->vtk_polydata->GetPointData()->GetArrayName(i), data_name.c_str())) { - return true; - } +VtkData *VtkData::create_writer(const std::string &file_name) { + auto file_ext = file_name.substr(file_name.find_last_of(".") + 1); + bool reader = false; + if (file_ext == "vtp") { + return new VtkVtpData(file_name, reader); + } else if (file_ext == "vtu") { + return new VtkVtuData(file_name, reader); } - return false; -} - -int VtkVtpData::elem_type() const -{ - return impl->elem_type; + throw std::runtime_error( + "Error in VtkData::create_writer: the file '" + file_name + + "' has the extension '" + file_ext + "', which is not 'vtp' or 'vtu'."); } -int VtkVtpData::num_elems() const -{ - return impl->num_elems; -} - -int VtkVtpData::np_elem() const -{ - return impl->np_elem; -} - -int VtkVtpData::num_points() const -{ - return impl->num_points; -} - -void VtkVtpData::read_file(const std::string& file_name) -{ - impl->read_file(file_name); -} - -void VtkVtpData::set_connectivity(const int nsd, const Array& conn, const int pid) -{ - impl->set_connectivity(nsd, conn, pid); -} - -void VtkVtpData::set_element_data(const std::string& data_name, const Array& data) -{ - throw std::runtime_error("[VtkVtpData] set_element_data not implemented."); -} - -void VtkVtpData::set_element_data(const std::string& data_name, const Array& data) -{ - throw std::runtime_error("[VtkVtpData] set_element_data not implemented."); -} - -void VtkVtpData::set_point_data(const std::string& data_name, const Array& data) -{ - throw std::runtime_error("[VtkVtpData] set_point_data for Array not implemented."); -} +VtkVtpData::VtkVtpData() { create_grid(); } -void VtkVtpData::set_point_data(const std::string& data_name, const Array& data) -{ - throw std::runtime_error("[VtkVtpData] set_point_data Array not implemented."); -} - -void VtkVtpData::set_point_data(const std::string& data_name, const Vector& data) -{ - impl->set_point_data(data_name, data); -} - -void VtkVtpData::set_points(const Array& points) -{ - impl->set_points(points); -} - -void VtkVtpData::set_time_value(const double time) { - impl->set_time_value(time); -} +VtkVtpData::VtkVtpData(const std::string &file_name, bool reader) { + file_name_ = file_name; -void VtkVtpData::write() -{ - impl->write(file_name); -} - -///////////////////////////////////////////////////////////////// -// V t k V t u D a t a I m p l e m e n t a t i o n // -///////////////////////////////////////////////////////////////// - - -VtkVtuData::VtkVtuData() -{ - impl = new VtkVtuDataImpl; -} - -VtkVtuData::VtkVtuData(const std::string& file_name, bool reader) -{ - this->file_name = file_name; - impl = new VtkVtuDataImpl; if (reader) { - read_file(file_name); + read_file(file_name); } else { - impl->create_grid(); + create_grid(); } } -VtkVtuData::~VtkVtuData() -{ - delete impl; -} - -Array VtkVtuData::get_connectivity() const -{ - int num_elems = impl->num_elems; - int np_elem = impl->np_elem; - - Array conn(np_elem, num_elems); - - auto cell = vtkGenericCell::New(); - for (int i = 0; i < num_elems; i++) { - impl->vtk_ugrid->GetCell(i, cell); - auto num_cell_pts = cell->GetNumberOfPoints(); - for (int j = 0; j < num_cell_pts; j++) { - auto id = cell->PointIds->GetId(j); - conn(j,i) = id; - } - } - return conn; -} - -/// @brief Get a list of point data names. -std::vector VtkVtuData::get_point_data_names() -{ - std::vector data_names; - int num_arrays = impl->vtk_ugrid->GetPointData()->GetNumberOfArrays(); - - for (int i = 0; i < num_arrays; i++) { - auto array_name = impl->vtk_ugrid->GetPointData()->GetArrayName(i); - data_names.push_back(array_name); - } - - return data_names; -} - -void VtkVtuData::copy_cell_data(const std::string& data_name, Array& mesh_data) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_ugrid->GetCellData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { - return; - } - - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return; - } - - int num_comp = vtk_data->GetNumberOfComponents(); +void VtkVtpData::create_grid() { + vtk_polydata = vtkSmartPointer::New(); + vtk_data = vtk_polydata; - // Set the data. - for (int i = 0; i < num_data; i++) { - auto tuple = vtk_data->GetTuple(i); - for (int j = 0; j < num_comp; j++) { - mesh_data(j, i) = tuple[j]; - } - } + elem_type_ = -1; + num_elems_ = 0; + num_points_per_elem_ = 0; + num_points_ = 0; } -void VtkVtuData::copy_cell_data(const std::string& data_name, Vector& mesh_data) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_ugrid->GetCellData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { - return; - } - - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return; - } - - // Set the data. - for (int i = 0; i < num_data; i++) { - mesh_data[i] = vtk_data->GetValue(i); - } +void VtkVtpData::write() const { + auto writer = vtkSmartPointer::New(); + writer->SetInputDataObject(vtk_polydata); + writer->SetFileName(file_name_.c_str()); + writer->Write(); } -/// @brief Copy an array of cell data from an unstructured mesh into the given Vector. -// -void VtkVtuData::copy_cell_data(const std::string& data_name, Vector& mesh_data) -{ - auto vtk_data = vtkIntArray::SafeDownCast(impl->vtk_ugrid->GetCellData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { - return; - } - - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return; - } - - int num_comp = vtk_data->GetNumberOfComponents(); - - for (int i = 0; i < num_data; i++) { - auto tuple = vtk_data->GetTuple(i); - for (int j = 0; j < num_comp; j++) { - mesh_data(i) = tuple[j]; - } - } +void VtkVtpData::read_file_internal(const std::string &file_name) { + auto reader = vtkSmartPointer::New(); + reader->SetFileName(file_name.c_str()); + reader->Update(); + vtk_polydata = reader->GetOutput(); + vtk_data = vtk_polydata; } -/// @brief Copy an array of point data from an unstructured grid into the given Array. -// -void VtkVtuData::copy_point_data(const std::string& data_name, Array& mesh_data) -{ - - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_ugrid->GetPointData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { - return; - } - - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return; +int VtkVtpData::cell_type(int nsd, int np_elem) const { + if (np_elem == 2) { + return VTK_LINE; } - int num_comp = vtk_data->GetNumberOfComponents(); - - // Set the data. - for (int i = 0; i < num_data; i++) { - auto tuple = vtk_data->GetTuple(i); - for (int j = 0; j < num_comp; j++) { - mesh_data(j, i) = tuple[j]; + if (nsd == 2) { + switch (np_elem) { + case 3: + return VTK_TRIANGLE; + case 4: + return VTK_QUAD; + case 6: + return VTK_QUADRATIC_TRIANGLE; + case 8: + return VTK_QUADRATIC_QUAD; + case 9: + return VTK_BIQUADRATIC_QUAD; } - } -} - -void VtkVtuData::copy_point_data(const std::string& data_name, Vector& mesh_data) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_ugrid->GetPointData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { - return; - } - - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return; - } - - int num_comp = vtk_data->GetNumberOfComponents(); - - // Set the data. - for (int i = 0; i < num_data; i++) { - mesh_data[i] = vtk_data->GetValue(i); - } -} - -void VtkVtuData::copy_point_data(const std::string& data_name, Vector& mesh_data) -{ - auto vtk_data = vtkIntArray::SafeDownCast(impl->vtk_ugrid->GetPointData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { - return; - } - - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return; - } - - int num_comp = vtk_data->GetNumberOfComponents(); - - // Set the data. - for (int i = 0; i < num_data; i++) { - mesh_data[i] = vtk_data->GetValue(i); - } -} - -/// @brief Copy points into the given array. -// -void VtkVtuData::copy_points(Array& points) -{ - auto vtk_points = impl->vtk_ugrid->GetPoints(); - auto num_points = vtk_points->GetNumberOfPoints(); - Array points_array(3, num_points); - - double point[3]; - for (int i = 0; i < num_points; i++) { - vtk_points->GetPoint(i, point); - points(0,i) = point[0]; - points(1,i) = point[1]; - points(2,i) = point[2]; - } - - return; -} - -bool VtkVtuData::has_cell_data(const std::string& data_name) -{ - int num_arrays = impl->vtk_ugrid->GetCellData()->GetNumberOfArrays(); - - for (int i = 0; i < num_arrays; i++) { - if (!strcmp(impl->vtk_ugrid->GetCellData()->GetArrayName(i), data_name.c_str())) { - return true; + } else if (nsd == 3) { + switch (np_elem) { + case 3: + return VTK_TRIANGLE; + case 4: + return VTK_QUAD; + case 6: + return VTK_QUADRATIC_TRIANGLE; + case 8: + return VTK_HEXAHEDRON; + case 10: + return VTK_QUADRATIC_TETRA; + case 20: + return VTK_QUADRATIC_HEXAHEDRON; + case 27: + return VTK_TRIQUADRATIC_HEXAHEDRON; } } - return false; + throw std::runtime_error( + "Error in VtkVtpData::cell_type: no cell type for an element with " + + std::to_string(np_elem) + " points in " + std::to_string(nsd) + + " dimensions."); } -bool VtkVtuData::has_point_data(const std::string& data_name) -{ - int num_arrays = impl->vtk_ugrid->GetPointData()->GetNumberOfArrays(); - - for (int i = 0; i < num_arrays; i++) { - if (!strcmp(impl->vtk_ugrid->GetPointData()->GetArrayName(i), data_name.c_str())) { - return true; - } - } - - return false; +void VtkVtpData::insert_cell(int vtk_cell_type, + vtkSmartPointer elem_nodes) { + // A vtkPolyData derives the type of a polygon from its number of points, so + // vtk_cell_type is not needed here. + vtk_polydata->GetPolys()->InsertNextCell(elem_nodes); } -/// @brief Get an array of point data from an unstructured grid. -// -Array VtkVtuData::get_point_data(const std::string& data_name) -{ - auto vtk_data = vtkDoubleArray::SafeDownCast(impl->vtk_ugrid->GetPointData()->GetArray(data_name.c_str())); - if (vtk_data == nullptr) { - return Array(); - } +VtkVtuData::VtkVtuData() { create_grid(); } - int num_data = vtk_data->GetNumberOfTuples(); - if (num_data == 0) { - return Array(); - } - - int num_comp = vtk_data->GetNumberOfComponents(); - - // Set the data. - Array data(num_data, num_comp); - for (int i = 0; i < num_data; i++) { - auto tuple = vtk_data->GetTuple(i); - for (int j = 0; j < num_comp; j++) { - data(i, j) = tuple[j]; - } - } - - return data; -} - -Array VtkVtuData::get_points() const -{ - auto vtk_points = impl->vtk_ugrid->GetPoints(); - auto num_points = vtk_points->GetNumberOfPoints(); - Array points_array(3, num_points); +VtkVtuData::VtkVtuData(const std::string &file_name, bool reader) { + file_name_ = file_name; - double point[3]; - for (int i = 0; i < num_points; i++) { - vtk_points->GetPoint(i, point); - points_array(0,i) = point[0]; - points_array(1,i) = point[1]; - points_array(2,i) = point[2]; + if (reader) { + read_file(file_name); + } else { + create_grid(); } - - return points_array; -} - -int VtkVtuData::elem_type() const -{ - return impl->elem_type; -} - -int VtkVtuData::num_elems() const -{ - return impl->num_elems; -} - -int VtkVtuData::np_elem() const -{ - return impl->np_elem; } -int VtkVtuData::num_points() const -{ - return impl->num_points; -} - - -void VtkVtuData::read_file(const std::string& file_name) -{ - impl->read_file(file_name); -} - -void VtkVtuData::set_connectivity(const int nsd, const Array& conn, const int pid) -{ - impl->set_connectivity(nsd, conn, pid); -} - -void VtkVtuData::set_element_data(const std::string& data_name, const Array& data) -{ - auto data_array = vtkSmartPointer::New(); - impl->set_element_data(data_name, data, data_array); - //impl->set_element_data(data_name, data); -} - -void VtkVtuData::set_element_data(const std::string& data_name, const Array& data) -{ - auto data_array = vtkSmartPointer::New(); - impl->set_element_data(data_name, data, data_array); - //impl->set_element_data(data_name, data); -} - -void VtkVtuData::set_point_data(const std::string& data_name, const Array& data) -{ - impl->set_point_data(data_name, data); -} +void VtkVtuData::create_grid() { + vtk_ugrid = vtkSmartPointer::New(); + vtk_data = vtk_ugrid; -void VtkVtuData::set_point_data(const std::string& data_name, const Array& data) -{ - impl->set_point_data(data_name, data); + elem_type_ = -1; + num_elems_ = 0; + num_points_per_elem_ = 0; + num_points_ = 0; } -void VtkVtuData::set_point_data(const std::string& data_name, const Vector& data) -{ - impl->set_point_data(data_name, data); +void VtkVtuData::write() const { + auto writer = vtkSmartPointer::New(); + writer->SetInputDataObject(vtk_ugrid); + writer->SetFileName(file_name_.c_str()); + writer->Write(); } -void VtkVtuData::set_points(const Array& points) -{ - impl->set_points(points); +void VtkVtuData::read_file_internal(const std::string &file_name) { + auto reader = vtkSmartPointer::New(); + reader->SetFileName(file_name.c_str()); + reader->Update(); + vtk_ugrid = reader->GetOutput(); + vtk_data = vtk_ugrid; } -void VtkVtuData::set_time_value(const double time) { - impl->set_time_value(time); -} +int VtkVtuData::cell_type(int nsd, int np_elem) const { + if (np_elem == 2) { + return VTK_LINE; + } -void VtkVtuData::write() -{ - impl->write(file_name); + if (nsd == 2) { + switch (np_elem) { + case 3: + return VTK_TRIANGLE; + case 4: + return VTK_QUAD; + case 6: + return VTK_QUADRATIC_TRIANGLE; + case 8: + return VTK_QUADRATIC_QUAD; + case 9: + return VTK_BIQUADRATIC_QUAD; + } + } else if (nsd == 3) { + switch (np_elem) { + case 3: + return VTK_TRIANGLE; + case 4: + return VTK_TETRA; + case 6: + return VTK_WEDGE; + case 8: + return VTK_HEXAHEDRON; + case 10: + return VTK_QUADRATIC_TETRA; + case 20: + return VTK_QUADRATIC_HEXAHEDRON; + case 27: + return VTK_TRIQUADRATIC_HEXAHEDRON; + } + } + + throw std::runtime_error( + "Error in VtkVtuData::cell_type: no cell type for an element with " + + std::to_string(np_elem) + " points in " + std::to_string(nsd) + + " dimensions."); +} + +void VtkVtuData::insert_cell(int vtk_cell_type, + vtkSmartPointer elem_nodes) { + vtk_ugrid->InsertNextCell(vtk_cell_type, elem_nodes); } - diff --git a/Code/Source/solver/VtkData.h b/Code/Source/solver/VtkData.h index ca6c8fc89..2eadbee87 100644 --- a/Code/Source/solver/VtkData.h +++ b/Code/Source/solver/VtkData.h @@ -7,167 +7,477 @@ #include "Array.h" #include "Vector.h" -#include #include - +#include +#include + +#include +#include +#include +#include +#include + +/** + * @brief A mesh stored in one of the VTK XML file formats. + * + * The mesh is held as a vtkPointSet, and consists of the point coordinates, + * the element connectivity, and any number of named data arrays associated + * with the points, with the elements, or with the mesh as a whole (field + * data). The class provides the operations needed to read those from a file, + * to build them up before writing a file, and to copy them to and from the + * Array and Vector types used by the solver. + * + * The file format determines how the mesh is represented, which is what the + * derived classes provide: VtkVtpData holds a surface mesh as a vtkPolyData, + * VtkVtuData a volume mesh as a vtkUnstructuredGrid. Use create_reader() and + * create_writer() to obtain an object of the type matching a given file name. + */ class VtkData { public: - VtkData(); - virtual ~VtkData(); - - virtual Array get_connectivity() const = 0; - virtual Array get_points() const = 0; - virtual int num_elems() const = 0; - virtual int elem_type() const = 0; - virtual int np_elem() const = 0; - virtual int num_points() const = 0; - virtual void read_file(const std::string& file_name) = 0; - - virtual void set_element_data(const std::string& data_name, const Array& data) = 0; - virtual void set_element_data(const std::string& data_name, const Array& data) = 0; - - virtual void set_point_data(const std::string& data_name, const Array& data) = 0; - virtual void set_point_data(const std::string& data_name, const Array& data) = 0; - virtual void set_point_data(const std::string& data_name, const Vector& data) = 0; - - virtual void set_points(const Array& points) = 0; - virtual void set_connectivity(const int nsd, const Array& conn, const int pid = 0) = 0; - - /// @brief Store a time value as field data, using the VTK convention for - /// time meta-data in XML files. - /// - /// The value is written as a single-tuple Float64 field data array named - /// 'TimeValue'. VTK XML readers such as ParaView turn this array into the - /// pipeline time of the data object. - /// - /// @param[in] time The time value to associate with the data. - virtual void set_time_value(const double time) = 0; - - virtual bool has_cell_data(const std::string& data_name) = 0; - virtual bool has_point_data(const std::string& data_name) = 0; - - virtual void copy_points(Array& points) = 0; - - virtual void copy_point_data(const std::string& data_name, Array& mesh_data) = 0; - virtual void copy_point_data(const std::string& data_name, Vector& mesh_data) = 0; - virtual void copy_point_data(const std::string& data_name, Vector& mesh_data) = 0; - - virtual void copy_cell_data(const std::string& data_name, Array& mesh_data) = 0; - virtual void copy_cell_data(const std::string& data_name, Vector& mesh_data) = 0; - virtual void copy_cell_data(const std::string& data_name, Vector& mesh_data) = 0; - - virtual Array get_point_data(const std::string& data_name) = 0; - virtual std::vector get_point_data_names() = 0; - - virtual void write() = 0; - - static VtkData* create_reader(const std::string& file_name); - static VtkData* create_writer(const std::string& file_name); - - std::string file_name; + /** + * @brief Default constructor. + */ + VtkData() = default; + + /** + * @brief Virtual destructor. + */ + virtual ~VtkData() = default; + + /** + * @brief Read the mesh data from a VTK file. + */ + virtual void read_file(const std::string &file_name); + + /** + * @brief Create an empty grid. + * + * Derived classes need to override this function by implementing the + * appropriate logic to create an empty grid. This will include initializing + * the correct vtkPointSet object. + */ + virtual void create_grid() = 0; + + /** + * @brief Write the mesh data to a VTK file. + */ + virtual void write() const = 0; + + /** + * @brief Get the connectivity of the mesh elements. + * + * @return An array of size (num_points_per_elem, num_elems) containing the + * connectivity of the mesh elements. Each column corresponds to an + * element, and each row corresponds to a point index in that element. + */ + Array get_connectivity() const; + + /** + * @brief Get the points of the mesh. + * + * @return An array of size (3, num_points) containing the coordinates of + * the mesh points. Each column corresponds to a point, and each row + * corresponds to a coordinate (x, y, z). + */ + Array get_points() const; + + /** + * @brief Get the number of elements in the mesh. + */ + int num_elems() const; + + /** + * @brief Get the element type. + */ + int elem_type() const; + + /** + * @brief Get the number of points per element. + */ + int num_points_per_elem() const; + + /** + * @brief Get the number of points in the mesh. + */ + int num_points() const; + + /** + * @brief Set a double-valued element data array. + * + * @param[in] data_name The name of the data array to set. + * @param[in] data The data array to set. + */ + void set_element_data(const std::string &data_name, + const Array &data); + + /** + * @brief Set an int-valued element data array. + * + * @param[in] data_name The name of the data array to set. + * @param[in] data The data array to set. + */ + void set_element_data(const std::string &data_name, const Array &data); + + /** + * @brief Set a double-valued point data array. + * + * @param[in] data_name The name of the data array to set. + * @param[in] data The data array to set. + */ + void set_point_data(const std::string &data_name, + const Array &data); + + /** + * @brief Set an int-valued point data array. + * + * @param[in] data_name The name of the data array to set. + * @param[in] data The data array to set. + */ + void set_point_data(const std::string &data_name, const Array &data); + + /** + * @brief Set an int-valued point data array. + * + * @param[in] data_name The name of the data array to set. + * @param[in] data The data array to set. + */ + void set_point_data(const std::string &data_name, const Vector &data); + + /** + * @brief Set the point coordinates of the mesh. + */ + void set_points(const Array &points); + + /** + * @brief Set the mesh connectivity to define the elements. + * + * The elements are appended to those already defined, so a mesh made of + * several parts is built up by calling this once per part. + * + * @param[in] nsd The number of spatial dimensions, which together with the + * number of points per element determines the element type. + * @param[in] conn The connectivity, of size (num_points_per_elem, + * num_elems). Each column holds the point indices of one element. + */ + void set_connectivity(const int nsd, const Array &conn); + + /** + * @brief Store a time value as field data. + * + * The value is written as a single-tuple Float64 field data array named + * 'TimeValue'. VTK XML readers such as ParaView turn this array into the + * pipeline time of the data object. + * + * @param[in] time The time value to associate with the data. + */ + void set_time_value(const double time); + + /** + * @brief Check if a given cell data array exists. + */ + bool has_cell_data(const std::string &data_name) const; + + /** + * @brief Check if a given point data array exists. + */ + bool has_point_data(const std::string &data_name) const; + + /** + * @brief Copy the mesh points to an Array. + * + * @param[out] points The array to copy the mesh points into. It must be + * of size (3, num_points). + */ + void copy_points(Array &points) const; + + /** + * @brief Copy an array of point data from the mesh into the given Array. + * + * @param[in] data_name The name of the point data array to copy. + * @param[out] mesh_data The array to copy the point data into. It must be + * of size (num_components, num_points). + */ + void copy_point_data(const std::string &data_name, + Array &mesh_data) const; + + /** + * @brief Copy an array of point data from the mesh into the given Vector. + * + * @param[in] data_name The name of the point data array to copy. + * @param[out] mesh_data The vector to copy the point data into. It must be + * of size (num_points). + */ + void copy_point_data(const std::string &data_name, + Vector &mesh_data) const; + + /** + * @brief Copy an array of int-valued point data from the mesh into the + * given Vector. + * + * @param[in] data_name The name of the point data array to copy. + * @param[out] mesh_data The vector to copy the point data into. It must be + * of size (num_points). + */ + void copy_point_data(const std::string &data_name, + Vector &mesh_data) const; + + /** + * @brief Copy an array of cell data from the mesh into the given Array. + * + * @param[in] data_name The name of the cell data array to copy. + * @param[out] mesh_data The array to copy the cell data into. It must be + * of size (num_components, num_cells). + */ + void copy_cell_data(const std::string &data_name, + Array &mesh_data) const; + + /** + * @brief Copy an array of cell data from the mesh into the given Vector. + * + * @param[in] data_name The name of the cell data array to copy. + * @param[out] mesh_data The vector to copy the cell data into. It must be + * of size (num_points). + */ + void copy_cell_data(const std::string &data_name, + Vector &mesh_data) const; + + /** + * @brief Copy an array of int-valued cell data from the mesh into the given + * Vector. + * + * @param[in] data_name The name of the cell data array to copy. + * @param[out] mesh_data The vector to copy the cell data into. It must be + * of size (num_points). + */ + void copy_cell_data(const std::string &data_name, + Vector &mesh_data) const; + + /** + * @brief Get an array of point data from the mesh. + * + * @todo[michelebucelli] This should fall back onto copy_point_data. + */ + Array get_point_data(const std::string &data_name) const; + + /** + * @brief Get a list of point data names. + * + * @return A vector of strings containing the names of the point data + * arrays. + */ + std::vector get_point_data_names() const; + + /** + * @brief Get the dimensions of a cell data array. + * + * @param[in] data_name The name of the cell data array to get the + * dimensions of. + * + * @return A pair of integers representing the number of components and the + * number of tuples in the array. + */ + std::pair + get_cell_data_dimensions(const std::string &data_name) const; + + /** + * @brief Create an object to read a mesh from a VTK file. + * + * The concrete type is selected from the file extension, and the mesh is + * read as part of the construction. The file extension must be 'vtp' or + * 'vtu'. + * + * @param[in] file_name The name of the VTK file to read. + * + * @return A pointer to a newly allocated object holding the mesh read from + * the file. The caller owns the object and must delete it. + * + * @throws std::runtime_error if the file extension is not 'vtp' or 'vtu' + */ + static VtkData *create_reader(const std::string &file_name); + + /** + * @brief Create an object to write a mesh to a VTK file. + * + * The concrete type is selected from the file extension, and the object is + * created holding an empty mesh. The file extension must be 'vtp' or 'vtu'. + * The file itself is written by write(), once the mesh has been defined. + * + * @param[in] file_name The name of the VTK file to write. + * + * @return A pointer to a newly allocated object holding an empty mesh. The + * caller owns the object and must delete it. + * + * @throws std::runtime_error if the file extension is not 'vtp' or 'vtu' + */ + static VtkData *create_writer(const std::string &file_name); + + protected: + /** + * @brief Read the mesh data from a file. + * + * Derived classes need to override this function by implementing the + * appropriate reading logic. This will include selecting the correct VTK + * reader, and initializing vtk_data. + */ + virtual void read_file_internal(const std::string &file_name) = 0; + + /** + * @brief Get the VTK cell type for a given number of spatial dimensions and + * number of points per element. + * + * Derived classes must override this to implement the appropriate mapping + * from the number of spatial dimensions and number of points per element to + * the VTK cell type. + */ + virtual int cell_type(int nsd, int np_elem) const = 0; + + /** + * @brief Insert a new cell into the VTK data object. + * + * Derived classes must override this to implement the appropriate insertion + * call. + * + * @param[in] vtk_cell_type The VTK cell type of the element. + * @param[in] elem_nodes The list of point IDs that define the element. + */ + virtual void insert_cell(int vtk_cell_type, + vtkSmartPointer elem_nodes) = 0; + /** + * Pointer to the underlying VTK data object. The concrete type will be + * either vtkPolyData, for VTP files, or vtkUnstructuredGrid, for VTU files. + */ + vtkSmartPointer vtk_data; + + /// Filename that the mesh is read from, or written to. + std::string file_name_; + + /// Type of elements in the mesh. + int elem_type_ = -1; + + /// Number of elements. + int num_elems_ = 0; + + /// Number of points per element. + int num_points_per_elem_ = 0; + + /// Number of points. + int num_points_ = 0; }; +/** + * @brief A mesh stored in the VTK XML polygonal data format ('.vtp' files). + * + * The mesh is held as a vtkPolyData, and its elements are polygons: the + * surface meshes that define the faces of a volume mesh, and the line meshes + * used for one-dimensional domains. + */ class VtkVtpData : public VtkData { - public: - VtkVtpData(); - VtkVtpData(const std::string& file_name, bool reader=true); - ~VtkVtpData(); - - // Copy constructor - VtkVtpData(const VtkVtpData& other); - // Copy assignment operator - VtkVtpData& operator=(const VtkVtpData& other); - - virtual Array get_connectivity() const override; - virtual Array get_points() const override; - virtual int elem_type() const override; - virtual int num_elems() const override; - virtual int np_elem() const override; - virtual int num_points() const override; - virtual void read_file(const std::string& file_name) override; - - virtual void copy_points(Array& points) override; - - virtual void copy_point_data(const std::string& data_name, Array& mesh_data) override; - virtual void copy_point_data(const std::string& data_name, Vector& mesh_data) override; - virtual void copy_point_data(const std::string& data_name, Vector& mesh_data) override; - - virtual void copy_cell_data(const std::string& data_name, Array& mesh_data) override; - virtual void copy_cell_data(const std::string& data_name, Vector& mesh_data) override; - virtual void copy_cell_data(const std::string& data_name, Vector& mesh_data) override; - std::pair get_cell_data_dimensions(const std::string& data_name) const; - - virtual Array get_point_data(const std::string& data_name) override; - virtual std::vector get_point_data_names() override; - - virtual bool has_cell_data(const std::string& data_name) override; - virtual bool has_point_data(const std::string& data_name) override; - - virtual void set_connectivity(const int nsd, const Array& conn, const int pid = 0) override; - - virtual void set_element_data(const std::string& data_name, const Array& data) override; - virtual void set_element_data(const std::string& data_name, const Array& data) override; - - virtual void set_point_data(const std::string& data_name, const Array& data) override; - virtual void set_point_data(const std::string& data_name, const Array& data) override; - virtual void set_point_data(const std::string& data_name, const Vector& data) override; - - virtual void set_points(const Array& points) override; - virtual void set_time_value(const double time) override; - virtual void write() override; - - private: - class VtkVtpDataImpl; - VtkVtpDataImpl* impl; +public: + /** + * @brief Default constructor. + * + * Creates an empty VtkVtpData object with an empty vtkPolyData grid. + */ + VtkVtpData(); + + /** + * @brief Constructor. + * + * @param[in] file_name The name of the VTP file to read from or write to. + * @param[in] reader If true, the constructor reads the mesh data from the + * given file. If false, it creates an empty grid. + */ + VtkVtpData(const std::string &file_name, bool reader = true); + + /** + * @brief Create an empty grid. + */ + virtual void create_grid() override; + + /** + * @brief Write the mesh data to a file. + */ + virtual void write() const override; + +protected: + /** + * @brief Read the mesh data from a file. + */ + virtual void read_file_internal(const std::string &file_name) override; + + /** + * @brief Get the VTK cell type of a surface element. + */ + virtual int cell_type(int nsd, int np_elem) const override; + + /** + * @brief Insert a new cell into the vtkPolyData object. + */ + virtual void insert_cell(int vtk_cell_type, + vtkSmartPointer elem_nodes) override; + + /** + * The mesh, represented as a vtkPolyData object. + */ + vtkSmartPointer vtk_polydata; }; +/** + * @brief A mesh stored in the VTK XML unstructured grid format ('.vtu' files). + * + * The mesh is held as a vtkUnstructuredGrid, whose elements may be of any VTK + * cell type: the volume meshes the equations are solved on, and the results + * written at each saved time step. + */ class VtkVtuData : public VtkData { - public: - VtkVtuData(); - VtkVtuData(const std::string& file_name, bool reader=true); - ~VtkVtuData(); - - virtual Array get_connectivity() const override; - virtual int elem_type() const override; - virtual int num_elems() const override; - virtual int np_elem() const override; - virtual int num_points() const override; - virtual void read_file(const std::string& file_name) override; - - virtual void copy_points(Array& points) override; - - virtual void copy_point_data(const std::string& data_name, Array& mesh_data) override; - virtual void copy_point_data(const std::string& data_name, Vector& mesh_data) override; - virtual void copy_point_data(const std::string& data_name, Vector& mesh_data) override; - - virtual void copy_cell_data(const std::string& data_name, Array& mesh_data) override; - virtual void copy_cell_data(const std::string& data_name, Vector& mesh_data) override; - virtual void copy_cell_data(const std::string& data_name, Vector& mesh_data) override; - - virtual Array get_point_data(const std::string& data_name) override; - virtual std::vector get_point_data_names() override; - - virtual Array get_points() const override; - - virtual bool has_cell_data(const std::string& data_name) override; - virtual bool has_point_data(const std::string& data_name) override; - - virtual void set_connectivity(const int nsd, const Array& conn, const int pid = 0) override; - - virtual void set_element_data(const std::string& data_name, const Array& data) override; - virtual void set_element_data(const std::string& data_name, const Array& data) override; - - virtual void set_point_data(const std::string& data_name, const Array& data) override; - virtual void set_point_data(const std::string& data_name, const Array& data) override; - virtual void set_point_data(const std::string& data_name, const Vector& data) override; - - virtual void set_points(const Array& points) override; - virtual void set_time_value(const double time) override; - virtual void write() override; - - private: - class VtkVtuDataImpl; - VtkVtuDataImpl* impl; +public: + /** + * @brief Default constructor. + * + * Creates an empty VtkVtuData object with an empty vtkUnstructuredGrid grid. + */ + VtkVtuData(); + + /** + * @brief Constructor. + * + * @param[in] file_name The name of the VTP file to read from or write to. + * @param[in] reader If true, the constructor reads the mesh data from the + * given file. If false, it creates an empty grid. + */ + VtkVtuData(const std::string &file_name, bool reader = true); + + /** + * @brief Create an empty grid. + */ + virtual void create_grid() override; + + /** + * @brief Write the mesh data to a file. + */ + virtual void write() const override; + +protected: + /** + * @brief Read the mesh data from a file. + */ + virtual void read_file_internal(const std::string &file_name) override; + + /** + * @brief Get the VTK cell type of a volume element. + */ + virtual int cell_type(int nsd, int np_elem) const override; + + /** + * @brief Insert a new cell into the vtkUnstructuredGrid object. + */ + virtual void insert_cell(int vtk_cell_type, + vtkSmartPointer elem_nodes) override; + + /** + * The mesh, represented as a vtkUnstructuredGrid object. + */ + vtkSmartPointer vtk_ugrid; }; #endif diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 790b1742a..abfa5de4f 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -581,8 +581,8 @@ void read_vtu(const std::string& file_name, mshType& mesh) #define n_read_vtu_use_VtkData #ifdef read_vtu_use_VtkData auto vtk_data = VtkData::create_reader(file_name); - int num_elems = vtk_data->num_elems(); - int np_elem = vtk_data->np_elem(); + int num_elems = vtk_data->num_elems(); + int np_elem = vtk_data->num_points_per_elem(); int elem_type = vtk_data->elem_type(); // Set mesh data. @@ -629,7 +629,7 @@ void read_precomputed_solution_vtu(const std::string& file_name, const std::stri #ifdef read_vtu_use_VtkData auto vtk_data = VtkData::create_reader(file_name); int num_elems = vtk_data->num_elems(); - int np_elem = vtk_data->np_elem(); + int np_elem = vtk_data->num_points_per_elem(); // Set mesh data. mesh.nEl = num_elems; From 373e9cb3c8cc43afca53534b77d1c18aea414b1e Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 1 Sep 2026 18:27:35 -0500 Subject: [PATCH 10/14] Add and document assertions to VtkData, VtkVtpData and VtkVtuData --- Code/Source/solver/VtkData.cpp | 295 ++++++++++++++++++++++----------- Code/Source/solver/VtkData.h | 77 ++++++++- 2 files changed, 266 insertions(+), 106 deletions(-) diff --git a/Code/Source/solver/VtkData.cpp b/Code/Source/solver/VtkData.cpp index 97d8fa97a..47861d8ec 100644 --- a/Code/Source/solver/VtkData.cpp +++ b/Code/Source/solver/VtkData.cpp @@ -3,14 +3,17 @@ #include "VtkData.h" +#include "Core/Exception.h" +#include "FE/Common/FEException.h" + #include -#include #include #include #include #include #include +#include #include #include #include @@ -28,12 +31,15 @@ void VtkData::read_file(const std::string &file_name) { // Extract metadata (number of cells, points, cell type). { + // A failed read leaves the data object empty, without any points. + const auto points = vtk_data->GetPoints(); + num_points_ = (points == nullptr) ? 0 : points->GetNumberOfPoints(); + svmp::check(num_points_ != 0, file_name, + "The file has no points."); + num_elems_ = vtk_data->GetNumberOfCells(); - num_points_ = vtk_data->GetPoints()->GetNumberOfPoints(); - if (num_points_ == 0) { - throw std::runtime_error("Error reading the VTK file '" + file_name + - "'."); - } + svmp::check(num_elems_ != 0, file_name, + "The file has no elements."); // Get the cell type. auto cell = vtkGenericCell::New(); @@ -90,6 +96,12 @@ void VtkData::set_element_data(const std::string &data_name, const int num_vals = data.ncols(); const int num_components = data.nrows(); + svmp::check( + num_vals == vtk_data->GetNumberOfCells(), + "The element data array named '" + data_name + "' holds " + + std::to_string(num_vals) + " values, while the mesh has " + + std::to_string(vtk_data->GetNumberOfCells()) + " elements."); + auto data_array = vtkSmartPointer::New(); data_array->SetNumberOfComponents(num_components); data_array->Allocate(num_vals, 1000); @@ -109,6 +121,12 @@ void VtkData::set_element_data(const std::string &data_name, const int num_vals = data.ncols(); const int num_components = data.nrows(); + svmp::check( + num_vals == vtk_data->GetNumberOfCells(), + "The element data array named '" + data_name + "' holds " + + std::to_string(num_vals) + " values, while the mesh has " + + std::to_string(vtk_data->GetNumberOfCells()) + " elements."); + auto data_array = vtkSmartPointer::New(); data_array->SetNumberOfComponents(num_components); data_array->Allocate(num_vals, 1000); @@ -128,6 +146,12 @@ void VtkData::set_point_data(const std::string &data_name, const int num_vals = data.ncols(); const int num_comp = data.nrows(); + svmp::check( + num_vals == vtk_data->GetNumberOfPoints(), + "The point data array named '" + data_name + "' holds " + + std::to_string(num_vals) + " values, while the mesh has " + + std::to_string(vtk_data->GetNumberOfPoints()) + " points."); + auto data_array = vtkSmartPointer::New(); data_array->SetNumberOfComponents(num_comp); data_array->Allocate(num_vals, 1000); @@ -148,6 +172,12 @@ void VtkData::set_point_data(const std::string &data_name, int num_vals = data.ncols(); int num_comp = data.nrows(); + svmp::check( + num_vals == vtk_data->GetNumberOfPoints(), + "The point data array named '" + data_name + "' holds " + + std::to_string(num_vals) + " values, while the mesh has " + + std::to_string(vtk_data->GetNumberOfPoints()) + " points."); + auto data_array = vtkSmartPointer::New(); data_array->SetNumberOfComponents(num_comp); data_array->Allocate(num_vals, 1000); @@ -167,6 +197,12 @@ void VtkData::set_point_data(const std::string &data_name, const Vector &data) { const int num_vals = data.size(); + svmp::check( + num_vals == vtk_data->GetNumberOfPoints(), + "The point data array named '" + data_name + "' holds " + + std::to_string(num_vals) + " values, while the mesh has " + + std::to_string(vtk_data->GetNumberOfPoints()) + " points."); + auto data_array = vtkSmartPointer::New(); data_array->SetNumberOfComponents(1); data_array->Allocate(num_vals); @@ -181,10 +217,14 @@ void VtkData::set_point_data(const std::string &data_name, void VtkData::set_points(const Array &points) { const int num_coords = points.ncols(); - if (num_coords == 0) { - throw std::runtime_error( - "Error in vtkData::set_points: the number of points is zero."); - } + svmp::check( + num_coords != 0, "The number of points is zero."); + + svmp::check( + points.nrows() >= 3, + "The point coordinates are given as an array of " + + std::to_string(points.nrows()) + + " rows, while three coordinates per point are needed."); auto node_coords = vtkSmartPointer::New(); node_coords->Allocate(num_coords, 1000); @@ -201,6 +241,8 @@ void VtkData::set_connectivity(const int nsd, const Array &conn) { int num_elems = conn.ncols(); int np_elem = conn.nrows(); + const vtkIdType num_points = vtk_data->GetNumberOfPoints(); + auto elem_nodes = vtkSmartPointer::New(); elem_nodes->Allocate(np_elem); elem_nodes->Initialize(); @@ -208,7 +250,18 @@ void VtkData::set_connectivity(const int nsd, const Array &conn) { for (int i = 0; i < num_elems; i++) { for (int j = 0; j < np_elem; j++) { - elem_nodes->SetId(j, conn(j, i)); + const int node_id = conn(j, i); + + // The check is written out to keep the error message from being composed + // for every point of every element. + if (node_id < 0 || node_id >= num_points) { + svmp::raise( + "Element " + std::to_string(i) + " refers to point " + + std::to_string(node_id) + ", which is not among the " + + std::to_string(num_points) + " points of the mesh."); + } + + elem_nodes->SetId(j, node_id); } insert_cell(cell_type(nsd, np_elem), elem_nodes); @@ -230,7 +283,8 @@ bool VtkData::has_cell_data(const std::string &data_name) const { const int num_arrays = vtk_data->GetCellData()->GetNumberOfArrays(); for (int i = 0; i < num_arrays; i++) { - if (!strcmp(vtk_data->GetCellData()->GetArrayName(i), data_name.c_str())) { + const char *array_name = vtk_data->GetCellData()->GetArrayName(i); + if (array_name != nullptr && !strcmp(array_name, data_name.c_str())) { return true; } } @@ -242,7 +296,8 @@ bool VtkData::has_point_data(const std::string &data_name) const { const int num_arrays = vtk_data->GetPointData()->GetNumberOfArrays(); for (int i = 0; i < num_arrays; i++) { - if (!strcmp(vtk_data->GetPointData()->GetArrayName(i), data_name.c_str())) { + const char *array_name = vtk_data->GetPointData()->GetArrayName(i); + if (array_name != nullptr && !strcmp(array_name, data_name.c_str())) { return true; } } @@ -254,6 +309,13 @@ void VtkData::copy_points(Array &points) const { auto vtk_points = vtk_data->GetPoints(); auto num_points = vtk_points->GetNumberOfPoints(); + svmp::check( + points.nrows() >= 3 && points.ncols() >= num_points, + "The " + std::to_string(num_points) + " points of the VTK file '" + + file_name_ + "' do not fit in an array of " + + std::to_string(points.nrows()) + " rows and " + + std::to_string(points.ncols()) + " columns."); + double point[3]; for (int i = 0; i < num_points; i++) { vtk_points->GetPoint(i, point); @@ -267,19 +329,23 @@ void VtkData::copy_point_data(const std::string &data_name, Array &mesh_data) const { const auto vtk_array = vtkDoubleArray::SafeDownCast( vtk_data->GetPointData()->GetArray(data_name.c_str())); - if (vtk_array == nullptr) { - // @todo[michelebucelli] This should probably be an exception. - return; - } + svmp::check( + vtk_array != nullptr, + "There is no double-valued point data array named '" + data_name + + "' in the VTK file '" + file_name_ + "'."); const int num_data = vtk_array->GetNumberOfTuples(); - if (num_data == 0) { - // @todo[michelebucelli] This should probably be an exception. - return; - } - const int num_comp = vtk_array->GetNumberOfComponents(); + svmp::check( + num_data <= mesh_data.ncols() && num_comp <= mesh_data.nrows(), + "The point data array named '" + data_name + "' of the VTK file '" + + file_name_ + "' has " + std::to_string(num_comp) + + " components and " + std::to_string(num_data) + + " tuples, which do not fit in an array of " + + std::to_string(mesh_data.nrows()) + " rows and " + + std::to_string(mesh_data.ncols()) + " columns."); + // Set the data. for (int i = 0; i < num_data; i++) { const auto tuple = vtk_array->GetTuple(i); @@ -293,14 +359,19 @@ void VtkData::copy_point_data(const std::string &data_name, Vector &mesh_data) const { const auto vtk_array = vtkDoubleArray::SafeDownCast( vtk_data->GetPointData()->GetArray(data_name.c_str())); - if (vtk_array == nullptr) { - return; - } + svmp::check( + vtk_array != nullptr, + "There is no double-valued point data array named '" + data_name + + "' in the VTK file '" + file_name_ + "'."); int num_data = vtk_array->GetNumberOfTuples(); - if (num_data == 0) { - return; - } + + svmp::check( + num_data <= mesh_data.size(), + "The point data array named '" + data_name + "' of the VTK file '" + + file_name_ + "' has " + std::to_string(num_data) + + " values, which do not fit in a vector of size " + + std::to_string(mesh_data.size()) + "."); // Set the data. for (int i = 0; i < num_data; i++) { @@ -312,14 +383,19 @@ void VtkData::copy_point_data(const std::string &data_name, Vector &mesh_data) const { const auto vtk_array = vtkIntArray::SafeDownCast( vtk_data->GetPointData()->GetArray(data_name.c_str())); - if (vtk_array == nullptr) { - return; - } + svmp::check( + vtk_array != nullptr, "There is no int-valued point data array named '" + + data_name + "' in the VTK file '" + file_name_ + + "'."); int num_data = vtk_array->GetNumberOfTuples(); - if (num_data == 0) { - return; - } + + svmp::check( + num_data <= mesh_data.size(), + "The point data array named '" + data_name + "' of the VTK file '" + + file_name_ + "' has " + std::to_string(num_data) + + " values, which do not fit in a vector of size " + + std::to_string(mesh_data.size()) + "."); // Set the data. for (int i = 0; i < num_data; i++) { @@ -331,17 +407,23 @@ void VtkData::copy_cell_data(const std::string &data_name, Array &mesh_data) const { const auto vtk_array = vtkDoubleArray::SafeDownCast( vtk_data->GetCellData()->GetArray(data_name.c_str())); - if (vtk_array == nullptr) { - return; - } + svmp::check( + vtk_array != nullptr, + "There is no double-valued element data array named '" + data_name + + "' in the VTK file '" + file_name_ + "'."); const int num_data = vtk_array->GetNumberOfTuples(); - if (num_data == 0) { - return; - } - const int num_comp = vtk_array->GetNumberOfComponents(); + svmp::check( + num_data <= mesh_data.ncols() && num_comp <= mesh_data.nrows(), + "The element data array named '" + data_name + "' of the VTK file '" + + file_name_ + "' has " + std::to_string(num_comp) + + " components and " + std::to_string(num_data) + + " tuples, which do not fit in an array of " + + std::to_string(mesh_data.nrows()) + " rows and " + + std::to_string(mesh_data.ncols()) + " columns."); + // Set the data. for (int i = 0; i < num_data; i++) { const auto tuple = vtk_array->GetTuple(i); @@ -355,14 +437,19 @@ void VtkData::copy_cell_data(const std::string &data_name, Vector &mesh_data) const { const auto vtk_array = vtkDoubleArray::SafeDownCast( vtk_data->GetCellData()->GetArray(data_name.c_str())); - if (vtk_array == nullptr) { - return; - } + svmp::check( + vtk_array != nullptr, + "There is no double-valued element data array named '" + data_name + + "' in the VTK file '" + file_name_ + "'."); const int num_data = vtk_array->GetNumberOfTuples(); - if (num_data == 0) { - return; - } + + svmp::check( + num_data <= mesh_data.size(), + "The element data array named '" + data_name + "' of the VTK file '" + + file_name_ + "' has " + std::to_string(num_data) + + " values, which do not fit in a vector of size " + + std::to_string(mesh_data.size()) + "."); // Set the data. for (int i = 0; i < num_data; i++) { @@ -374,14 +461,19 @@ void VtkData::copy_cell_data(const std::string &data_name, Vector &mesh_data) const { const auto vtk_array = vtkIntArray::SafeDownCast( vtk_data->GetCellData()->GetArray(data_name.c_str())); - if (vtk_array == nullptr) { - return; - } + svmp::check( + vtk_array != nullptr, + "There is no int-valued element data array named '" + data_name + + "' in the VTK file '" + file_name_ + "'."); const int num_data = vtk_array->GetNumberOfTuples(); - if (num_data == 0) { - return; - } + + svmp::check( + num_data <= mesh_data.size(), + "The element data array named '" + data_name + "' of the VTK file '" + + file_name_ + "' has " + std::to_string(num_data) + + " values, which do not fit in a vector of size " + + std::to_string(mesh_data.size()) + "."); // Set the data. for (int i = 0; i < num_data; i++) { @@ -392,15 +484,12 @@ void VtkData::copy_cell_data(const std::string &data_name, Array VtkData::get_point_data(const std::string &data_name) const { auto vtk_array = vtkDoubleArray::SafeDownCast( vtk_data->GetPointData()->GetArray(data_name.c_str())); - if (vtk_array == nullptr) { - return Array(); - } + svmp::check( + vtk_array != nullptr, + "There is no double-valued point data array named '" + data_name + + "' in the VTK file '" + file_name_ + "'."); int num_data = vtk_array->GetNumberOfTuples(); - if (num_data == 0) { - return Array(); - } - int num_comp = vtk_array->GetNumberOfComponents(); // Set the data. @@ -445,7 +534,7 @@ VtkData *VtkData::create_reader(const std::string &file_name) { return new VtkVtuData(file_name); } - throw std::runtime_error( + svmp::raise( "Error in VtkData::create_reader: the file '" + file_name + "' has the extension '" + file_ext + "', which is not 'vtp' or 'vtu'."); } @@ -459,7 +548,7 @@ VtkData *VtkData::create_writer(const std::string &file_name) { return new VtkVtuData(file_name, reader); } - throw std::runtime_error( + svmp::raise( "Error in VtkData::create_writer: the file '" + file_name + "' has the extension '" + file_ext + "', which is not 'vtp' or 'vtu'."); } @@ -490,13 +579,25 @@ void VtkVtpData::write() const { auto writer = vtkSmartPointer::New(); writer->SetInputDataObject(vtk_polydata); writer->SetFileName(file_name_.c_str()); - writer->Write(); + + const int status = writer->Write(); + svmp::check( + status != 0, + "Error writing the VTK file '" + file_name_ + "': " + + vtkErrorCode::GetStringFromErrorCode(writer->GetErrorCode()) + ".", + svmp::StatusCode::IOError); } void VtkVtpData::read_file_internal(const std::string &file_name) { auto reader = vtkSmartPointer::New(); reader->SetFileName(file_name.c_str()); reader->Update(); + + svmp::check( + reader->GetErrorCode() == vtkErrorCode::NoError, file_name, + std::string("Error reading VTK file: ") + + vtkErrorCode::GetStringFromErrorCode(reader->GetErrorCode()) + "."); + vtk_polydata = reader->GetOutput(); vtk_data = vtk_polydata; } @@ -506,40 +607,24 @@ int VtkVtpData::cell_type(int nsd, int np_elem) const { return VTK_LINE; } - if (nsd == 2) { - switch (np_elem) { - case 3: - return VTK_TRIANGLE; - case 4: - return VTK_QUAD; - case 6: - return VTK_QUADRATIC_TRIANGLE; - case 8: - return VTK_QUADRATIC_QUAD; - case 9: - return VTK_BIQUADRATIC_QUAD; - } - } else if (nsd == 3) { - switch (np_elem) { - case 3: - return VTK_TRIANGLE; - case 4: - return VTK_QUAD; - case 6: - return VTK_QUADRATIC_TRIANGLE; - case 8: - return VTK_HEXAHEDRON; - case 10: - return VTK_QUADRATIC_TETRA; - case 20: - return VTK_QUADRATIC_HEXAHEDRON; - case 27: - return VTK_TRIQUADRATIC_HEXAHEDRON; - } - } - - throw std::runtime_error( - "Error in VtkVtpData::cell_type: no cell type for an element with " + + // A vtkPolyData holds surface elements only, whose type is determined by the + // number of points alone. + switch (np_elem) { + case 3: + return VTK_TRIANGLE; + case 4: + return VTK_QUAD; + case 6: + return VTK_QUADRATIC_TRIANGLE; + case 8: + return VTK_QUADRATIC_QUAD; + case 9: + return VTK_BIQUADRATIC_QUAD; + } + + svmp::raise( + "Error in VtkVtpData::cell_type: no surface cell type for an element " + "with " + std::to_string(np_elem) + " points in " + std::to_string(nsd) + " dimensions."); } @@ -577,13 +662,25 @@ void VtkVtuData::write() const { auto writer = vtkSmartPointer::New(); writer->SetInputDataObject(vtk_ugrid); writer->SetFileName(file_name_.c_str()); - writer->Write(); + + const int status = writer->Write(); + svmp::check( + status != 0, + "Error writing the VTK file '" + file_name_ + "': " + + vtkErrorCode::GetStringFromErrorCode(writer->GetErrorCode()) + ".", + svmp::StatusCode::IOError); } void VtkVtuData::read_file_internal(const std::string &file_name) { auto reader = vtkSmartPointer::New(); reader->SetFileName(file_name.c_str()); reader->Update(); + + svmp::check( + reader->GetErrorCode() == vtkErrorCode::NoError, file_name, + std::string("Error reading VTK file: ") + + vtkErrorCode::GetStringFromErrorCode(reader->GetErrorCode()) + "."); + vtk_ugrid = reader->GetOutput(); vtk_data = vtk_ugrid; } @@ -625,7 +722,7 @@ int VtkVtuData::cell_type(int nsd, int np_elem) const { } } - throw std::runtime_error( + svmp::raise( "Error in VtkVtuData::cell_type: no cell type for an element with " + std::to_string(np_elem) + " points in " + std::to_string(nsd) + " dimensions."); diff --git a/Code/Source/solver/VtkData.h b/Code/Source/solver/VtkData.h index 2eadbee87..074544948 100644 --- a/Code/Source/solver/VtkData.h +++ b/Code/Source/solver/VtkData.h @@ -46,6 +46,9 @@ class VtkData { /** * @brief Read the mesh data from a VTK file. + * + * @throws svmp::FileFormatException if the file cannot be read, or if it + * contains no points or no elements */ virtual void read_file(const std::string &file_name); @@ -60,6 +63,8 @@ class VtkData { /** * @brief Write the mesh data to a VTK file. + * + * @throws svmp::CoreException if the file cannot be written */ virtual void write() const = 0; @@ -105,7 +110,10 @@ class VtkData { * @brief Set a double-valued element data array. * * @param[in] data_name The name of the data array to set. - * @param[in] data The data array to set. + * @param[in] data The data array to set, holding one value per element. + * + * @throws svmp::FE::InvalidArgumentException if the number of values + * differs from the number of elements of the mesh */ void set_element_data(const std::string &data_name, const Array &data); @@ -114,7 +122,10 @@ class VtkData { * @brief Set an int-valued element data array. * * @param[in] data_name The name of the data array to set. - * @param[in] data The data array to set. + * @param[in] data The data array to set, holding one value per element. + * + * @throws svmp::FE::InvalidArgumentException if the number of values + * differs from the number of elements of the mesh */ void set_element_data(const std::string &data_name, const Array &data); @@ -122,7 +133,10 @@ class VtkData { * @brief Set a double-valued point data array. * * @param[in] data_name The name of the data array to set. - * @param[in] data The data array to set. + * @param[in] data The data array to set, holding one value per point. + * + * @throws svmp::FE::InvalidArgumentException if the number of values + * differs from the number of points of the mesh */ void set_point_data(const std::string &data_name, const Array &data); @@ -131,7 +145,10 @@ class VtkData { * @brief Set an int-valued point data array. * * @param[in] data_name The name of the data array to set. - * @param[in] data The data array to set. + * @param[in] data The data array to set, holding one value per point. + * + * @throws svmp::FE::InvalidArgumentException if the number of values + * differs from the number of points of the mesh */ void set_point_data(const std::string &data_name, const Array &data); @@ -139,12 +156,20 @@ class VtkData { * @brief Set an int-valued point data array. * * @param[in] data_name The name of the data array to set. - * @param[in] data The data array to set. + * @param[in] data The data array to set, holding one value per point. + * + * @throws svmp::FE::InvalidArgumentException if the number of values + * differs from the number of points of the mesh */ void set_point_data(const std::string &data_name, const Vector &data); /** * @brief Set the point coordinates of the mesh. + * + * @param[in] points The coordinates, of size (3, num_points). + * + * @throws svmp::FE::InvalidArgumentException if there are no points or + * fewer than three coordinates are given for each point */ void set_points(const Array &points); @@ -158,6 +183,9 @@ class VtkData { * number of points per element determines the element type. * @param[in] conn The connectivity, of size (num_points_per_elem, * num_elems). Each column holds the point indices of one element. + * + * @throws svmp::FE::InvalidArgumentException if a point index does not + * refer to one of the points of the mesh */ void set_connectivity(const int nsd, const Array &conn); @@ -187,6 +215,9 @@ class VtkData { * * @param[out] points The array to copy the mesh points into. It must be * of size (3, num_points). + * + * @throws svmp::FE::InvalidArgumentException if the points do not fit in + * the array */ void copy_points(Array &points) const; @@ -196,6 +227,10 @@ class VtkData { * @param[in] data_name The name of the point data array to copy. * @param[out] mesh_data The array to copy the point data into. It must be * of size (num_components, num_points). + * + * @throws svmp::FE::InvalidArgumentException if the mesh has no + * double-valued point data array with the given name, or if its values do + * not fit in mesh_data */ void copy_point_data(const std::string &data_name, Array &mesh_data) const; @@ -206,6 +241,10 @@ class VtkData { * @param[in] data_name The name of the point data array to copy. * @param[out] mesh_data The vector to copy the point data into. It must be * of size (num_points). + * + * @throws svmp::FE::InvalidArgumentException if the mesh has no + * double-valued point data array with the given name, or if its values do + * not fit in mesh_data */ void copy_point_data(const std::string &data_name, Vector &mesh_data) const; @@ -217,6 +256,10 @@ class VtkData { * @param[in] data_name The name of the point data array to copy. * @param[out] mesh_data The vector to copy the point data into. It must be * of size (num_points). + * + * @throws svmp::FE::InvalidArgumentException if the mesh has no int-valued + * point data array with the given name, or if its values do not fit in + * mesh_data */ void copy_point_data(const std::string &data_name, Vector &mesh_data) const; @@ -227,6 +270,10 @@ class VtkData { * @param[in] data_name The name of the cell data array to copy. * @param[out] mesh_data The array to copy the cell data into. It must be * of size (num_components, num_cells). + * + * @throws svmp::FE::InvalidArgumentException if the mesh has no + * double-valued cell data array with the given name, or if its values do + * not fit in mesh_data */ void copy_cell_data(const std::string &data_name, Array &mesh_data) const; @@ -237,6 +284,10 @@ class VtkData { * @param[in] data_name The name of the cell data array to copy. * @param[out] mesh_data The vector to copy the cell data into. It must be * of size (num_points). + * + * @throws svmp::FE::InvalidArgumentException if the mesh has no + * double-valued cell data array with the given name, or if its values do + * not fit in mesh_data */ void copy_cell_data(const std::string &data_name, Vector &mesh_data) const; @@ -248,6 +299,10 @@ class VtkData { * @param[in] data_name The name of the cell data array to copy. * @param[out] mesh_data The vector to copy the cell data into. It must be * of size (num_points). + * + * @throws svmp::FE::InvalidArgumentException if the mesh has no int-valued + * cell data array with the given name, or if its values do not fit in + * mesh_data */ void copy_cell_data(const std::string &data_name, Vector &mesh_data) const; @@ -255,6 +310,9 @@ class VtkData { /** * @brief Get an array of point data from the mesh. * + * @throws svmp::FE::InvalidArgumentException if the mesh has no + * double-valued point data array with the given name + * * @todo[michelebucelli] This should fall back onto copy_point_data. */ Array get_point_data(const std::string &data_name) const; @@ -291,7 +349,8 @@ class VtkData { * @return A pointer to a newly allocated object holding the mesh read from * the file. The caller owns the object and must delete it. * - * @throws std::runtime_error if the file extension is not 'vtp' or 'vtu' + * @throws svmp::FE::InvalidArgumentException if the file extension is not + * 'vtp' or 'vtu' */ static VtkData *create_reader(const std::string &file_name); @@ -307,7 +366,8 @@ class VtkData { * @return A pointer to a newly allocated object holding an empty mesh. The * caller owns the object and must delete it. * - * @throws std::runtime_error if the file extension is not 'vtp' or 'vtu' + * @throws svmp::FE::InvalidArgumentException if the file extension is not + * 'vtp' or 'vtu' */ static VtkData *create_writer(const std::string &file_name); @@ -328,6 +388,9 @@ class VtkData { * Derived classes must override this to implement the appropriate mapping * from the number of spatial dimensions and number of points per element to * the VTK cell type. + * + * @throws svmp::FE::InvalidArgumentException if there is no cell type that + * the file format can hold for the given element */ virtual int cell_type(int nsd, int np_elem) const = 0; From fd23aef30aa78114de194588dfc1f46e35a56f58 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 1 Sep 2026 18:28:38 -0500 Subject: [PATCH 11/14] Remove unnecessary downcast from VtkData methods operating on Array --- Code/Source/solver/VtkData.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/Code/Source/solver/VtkData.cpp b/Code/Source/solver/VtkData.cpp index 47861d8ec..858f3aa8d 100644 --- a/Code/Source/solver/VtkData.cpp +++ b/Code/Source/solver/VtkData.cpp @@ -327,8 +327,7 @@ void VtkData::copy_points(Array &points) const { void VtkData::copy_point_data(const std::string &data_name, Array &mesh_data) const { - const auto vtk_array = vtkDoubleArray::SafeDownCast( - vtk_data->GetPointData()->GetArray(data_name.c_str())); + const auto vtk_array = vtk_data->GetPointData()->GetArray(data_name.c_str()); svmp::check( vtk_array != nullptr, "There is no double-valued point data array named '" + data_name + @@ -405,8 +404,7 @@ void VtkData::copy_point_data(const std::string &data_name, void VtkData::copy_cell_data(const std::string &data_name, Array &mesh_data) const { - const auto vtk_array = vtkDoubleArray::SafeDownCast( - vtk_data->GetCellData()->GetArray(data_name.c_str())); + const auto vtk_array = vtk_data->GetCellData()->GetArray(data_name.c_str()); svmp::check( vtk_array != nullptr, "There is no double-valued element data array named '" + data_name + From 521466f9997efaae3f3ae69c8f0892ae7c361829 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 2 Sep 2026 09:17:23 -0500 Subject: [PATCH 12/14] Fix typo in comment --- tests/conftest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/conftest.py b/tests/conftest.py index e155f2ddd..67baf394e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -215,7 +215,7 @@ def run_with_reference( n_proc: number of processors t_max: time step to compare name_inp: name of svMultiPhysics input file (.xml) - name_ref: name of refence file (.vtu) + name_ref: name of reference file (.vtu) check_time_value: whether to compare the TimeValue field data against the time reached at time step t_max """ From 3d22c4e0213a9533d2f2bd13d6bf4c9ddc7a5af0 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 2 Sep 2026 09:21:27 -0500 Subject: [PATCH 13/14] Use vtkSmartPointer to allocate cell in VtkData::read_file --- Code/Source/solver/VtkData.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Source/solver/VtkData.cpp b/Code/Source/solver/VtkData.cpp index 858f3aa8d..4108e2c63 100644 --- a/Code/Source/solver/VtkData.cpp +++ b/Code/Source/solver/VtkData.cpp @@ -42,7 +42,7 @@ void VtkData::read_file(const std::string &file_name) { "The file has no elements."); // Get the cell type. - auto cell = vtkGenericCell::New(); + auto cell = vtkSmartPointer::New(); vtk_data->GetCell(0, cell); num_points_per_elem_ = cell->GetNumberOfPoints(); elem_type_ = cell->GetCellType(); From 7d67876a45dbe46ebba204483fd1abf57511fc98 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 2 Sep 2026 09:37:45 -0500 Subject: [PATCH 14/14] write_vtp uses VtkData::set_element_data to assign the cell-based array GlobalElementID --- Code/Source/solver/VtkData.cpp | 22 ++++++++++++++++++++++ Code/Source/solver/VtkData.h | 12 ++++++++++++ Code/Source/solver/vtk_xml.cpp | 2 +- 3 files changed, 35 insertions(+), 1 deletion(-) diff --git a/Code/Source/solver/VtkData.cpp b/Code/Source/solver/VtkData.cpp index 4108e2c63..eab486ecf 100644 --- a/Code/Source/solver/VtkData.cpp +++ b/Code/Source/solver/VtkData.cpp @@ -141,6 +141,28 @@ void VtkData::set_element_data(const std::string &data_name, vtk_data->GetCellData()->AddArray(data_array); } +void VtkData::set_element_data(const std::string &data_name, + const Vector &data) { + const int num_vals = data.size(); + + svmp::check( + num_vals == vtk_data->GetNumberOfCells(), + "The element data array named '" + data_name + "' holds " + + std::to_string(num_vals) + " values, while the mesh has " + + std::to_string(vtk_data->GetNumberOfCells()) + " elements."); + + auto data_array = vtkSmartPointer::New(); + data_array->SetNumberOfComponents(1); + data_array->Allocate(num_vals); + data_array->SetName(data_name.c_str()); + + for (int i = 0; i < num_vals; ++i) { + data_array->InsertNextTuple1(data(i)); + } + + vtk_data->GetCellData()->AddArray(data_array); +} + void VtkData::set_point_data(const std::string &data_name, const Array &data) { const int num_vals = data.ncols(); diff --git a/Code/Source/solver/VtkData.h b/Code/Source/solver/VtkData.h index 074544948..8c8c26e80 100644 --- a/Code/Source/solver/VtkData.h +++ b/Code/Source/solver/VtkData.h @@ -129,6 +129,18 @@ class VtkData { */ void set_element_data(const std::string &data_name, const Array &data); + /** + * @brief Set an int-valued element data vector. + * + * @param[in] data_name The name of the data array to set. + * @param[in] data The data vector to set, holding one value per element. + * + * @throws svmp::FE::InvalidArgumentException if the number of values + * differs from the number of elements of the mesh + */ + void set_element_data(const std::string &data_name, + const Vector &data); + /** * @brief Set a double-valued point data array. * diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index abfa5de4f..34b8b22e8 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -841,7 +841,7 @@ void write_vtp(ComMod& com_mod, faceType& lFa, const std::string& fName) } if (lFa.gE.size() != 0) { - vtk_writer->set_point_data("GlobalElementID", lFa.gE); + vtk_writer->set_element_data("GlobalElementID", lFa.gE); } vtk_writer->write();