vtkdiff.cpp

/**
 * \copyright
 * Copyright (c) 2015-2025, OpenGeoSys Community (http://www.opengeosys.org)
 * Distributed under the BSD 3-Clause "New" or "Revised" license.
 * See accompanying file LICENSE.txt or
 * http://www.opengeosys.org/project/license
 */

#include <tclap/CmdLine.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCommand.h>
#include <vtkDataArray.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkSmartPointer.h>
#include <vtkUnstructuredGrid.h>
#include <vtkVersion.h>
#include <vtkXMLUnstructuredGridReader.h>

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <ios>
#include <iterator>
#include <optional>
#include <sstream>
#include <tuple>
#include <type_traits>

template <typename T>
auto float_to_string(T const& v) -> std::string
{
 static_assert(std::is_floating_point<T>::value,
 "float_to_string requires a floating point input type.");

 std::stringstream double_eps_sstream;
 double_eps_sstream << std::scientific << std::setprecision(16) << v;
 return double_eps_sstream.str();
}

bool stringEndsWith(std::string const& str, std::string const& ending)
{
 if (str.length() < ending.length())
 return false;

 // now the difference is non-negative, no underflow possible.
 auto const string_end_length = str.length() - ending.length();
 return str.compare(string_end_length, ending.length(), ending) == 0;
}

template <typename T>
std::ostream& operator<<(std::ostream& os, std::vector<T> const& vector)
{
 if (vector.empty())
 {
 return os << "[]";
 }

 // print first n-1 elements
 os << "[";
 std::size_t const size = vector.size();
 for (std::size_t i = 0; i < size - 1; ++i)
 {
 os << vector[i] << ", ";
 }
 return os << vector.back() << "]";
}

struct Args
{
 bool const quiet;
 bool const verbose;
 bool const meshcheck;
 double const abs_err_thr;
 double const rel_err_thr;
 std::string const vtk_input_a;
 std::string const vtk_input_b;
 std::string const data_array_a;
 std::string const data_array_b;
};

auto parseCommandLine(int argc, char* argv[]) -> Args
{
 TCLAP::CmdLine cmd(
 "VtkDiff software.\n"
 "Copyright (c) 2015-2025, OpenGeoSys Community "
 "(http://www.opengeosys.org) "
 "Distributed under the BSD 3-Clause 'New' or 'Revised' License. "
 "See accompanying file LICENSE.txt or "
 "http://www.opengeosys.org/project/license",
 ' ',
 "1.0.0");

 TCLAP::UnlabeledValueArg<std::string> vtk_input_a_arg(
 "input-file-a",
 "Path to the VTK unstructured grid input file.",
 true,
 "",
 "VTK FILE");
 cmd.add(vtk_input_a_arg);

 TCLAP::UnlabeledValueArg<std::string> vtk_input_b_arg(
 "input-file-b",
 "Path to the second VTK unstructured grid input file.",
 false,
 "",
 "VTK FILE");
 cmd.add(vtk_input_b_arg);

 TCLAP::ValueArg<std::string> data_array_a_arg(
 "a", "first_data_array", "First data array name for comparison", false,
 "", "NAME");

 TCLAP::ValueArg<std::string> data_array_b_arg(
 "b",
 "second_data_array",
 "Second data array name for comparison",
 false,
 "",
 "NAME");
 cmd.add(data_array_b_arg);

 TCLAP::SwitchArg meshcheck_arg(
 "m", "mesh_check", "Compare mesh geometries using absolute tolerance.");
 cmd.xorAdd(data_array_a_arg, meshcheck_arg);

 TCLAP::SwitchArg quiet_arg("q", "quiet", "Suppress all but error output.");
 cmd.add(quiet_arg);

 TCLAP::SwitchArg verbose_arg("v", "verbose",
 "Also print which values differ.");
 cmd.add(verbose_arg);

 auto const double_eps_string =
 float_to_string(std::numeric_limits<double>::epsilon());

 TCLAP::ValueArg<double> abs_err_thr_arg(
 "",
 "abs",
 "Tolerance for the absolute error in the maximum norm (" +
 double_eps_string + ")",
 false,
 std::numeric_limits<double>::epsilon(),
 "FLOAT");
 cmd.add(abs_err_thr_arg);

 TCLAP::ValueArg<double> rel_err_thr_arg(
 "",
 "rel",
 "Tolerance for the componentwise relative error (" + double_eps_string +
 ")",
 false,
 std::numeric_limits<double>::epsilon(),
 "FLOAT");
 cmd.add(rel_err_thr_arg);

 cmd.parse(argc, argv);

 return Args{quiet_arg.getValue(), verbose_arg.getValue(),
 meshcheck_arg.getValue(), abs_err_thr_arg.getValue(),
 rel_err_thr_arg.getValue(), vtk_input_a_arg.getValue(),
 vtk_input_b_arg.getValue(), data_array_a_arg.getValue(),
 data_array_b_arg.getValue()};
}

template <typename T>
class ErrorCallback : public vtkCommand
{
public:
 vtkTypeMacro(ErrorCallback, vtkCommand);

 static ErrorCallback<T>* New() { return new ErrorCallback<T>; }

 void Execute(vtkObject* caller, unsigned long vtkNotUsed(eventId),
 void* callData) override
 {
 auto* reader = static_cast<T*>(caller);
 std::cerr << "Error reading file `" << reader->GetFileName() << "'\n"
 << static_cast<char*>(callData) << "\nAborting." << std::endl;
 std::exit(2);
 }
};

enum class Domain
{
 Point,
 Cell,
 Field
};

vtkSmartPointer<vtkUnstructuredGrid> readMesh(std::string const& filename)
{
 if (filename.empty())
 {
 return nullptr;
 }

 if (!stringEndsWith(filename, ".vtu"))
 {
 std::cerr << "Error: Expected a file with .vtu extension."
 << "File '" << filename << "' not read.";
 return nullptr;
 }

 vtkSmartPointer<ErrorCallback<vtkXMLUnstructuredGridReader>> errorCallback =
 vtkSmartPointer<ErrorCallback<vtkXMLUnstructuredGridReader>>::New();

 vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
 vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
 reader->AddObserver(vtkCommand::ErrorEvent, errorCallback);
 reader->SetFileName(filename.c_str());
 reader->Update();
 return reader->GetOutput();
}

std::tuple<vtkSmartPointer<vtkUnstructuredGrid>,
 vtkSmartPointer<vtkUnstructuredGrid>>
readMeshes(std::string const& file_a_name, std::string const& file_b_name)
{
 return {readMesh(file_a_name), readMesh(file_b_name)};
}

std::optional<Domain> determineDomain(vtkUnstructuredGrid& mesh,
 std::string const& data_array_name)
{
 if (mesh.GetPointData()->HasArray(data_array_name.c_str()))
 {
 return Domain::Point;
 }
 else if (mesh.GetCellData()->HasArray(data_array_name.c_str()))
 {
 return Domain::Cell;
 }
 else if (mesh.GetFieldData()->HasArray(data_array_name.c_str()))
 {
 return Domain::Field;
 }
 else
 {
 std::cerr << "Error: Scalars data array "
 << "\'" << data_array_name.c_str() << "\'"
 << " neither found in point data nor in cell data nor in "
 "field data.\n";
 return std::nullopt;
 }
}

vtkSmartPointer<vtkDataArray> getDataArray(vtkUnstructuredGrid& mesh,
 std::string const& data_array_name,
 Domain domain)
{
 vtkSmartPointer<vtkDataArray> data_array;

 switch (domain)
 {
 case Domain::Point:
 data_array =
 mesh.GetPointData()->GetScalars(data_array_name.c_str());
 break;
 case Domain::Cell:
 data_array =
 mesh.GetCellData()->GetScalars(data_array_name.c_str());
 break;
 case Domain::Field:
 // GetArray() is not recommended for use, but according to the VTK
 // API docs it is OK here, because we want to get a vtkDataArray.
 data_array = mesh.GetFieldData()->GetArray(data_array_name.c_str());
 break;
 }

 // Check arrays' validity
 if (!data_array)
 {
 std::cerr << "Error: Scalars data array "
 << "\'" << data_array_name.c_str() << "\'"
 << " could not be read.\n";
 }

 return data_array;
}

std::tuple<bool, vtkSmartPointer<vtkDataArray>, vtkSmartPointer<vtkDataArray>>
readDataArraysFromMeshes(
 std::tuple<vtkSmartPointer<vtkUnstructuredGrid>,
 vtkSmartPointer<vtkUnstructuredGrid>> const& meshes,
 std::string const& data_array_a_name,
 std::string const& data_array_b_name)
{
 auto const& [mesh_a_ptr, mesh_b_ptr] = meshes;

 if (mesh_a_ptr == nullptr)
 {
 std::cerr << "First mesh was not read correctly and is a nullptr.\n";
 return {false, nullptr, nullptr};
 }

 auto const opt_domain = determineDomain(*mesh_a_ptr, data_array_a_name);

 if (!opt_domain)
 {
 return {false, nullptr, nullptr};
 }

 auto const domain = *opt_domain;

 // Get arrays
 vtkSmartPointer<vtkDataArray> a =
 getDataArray(*mesh_a_ptr, data_array_a_name, domain);

 if (!a)
 {
 return {false, nullptr, nullptr};
 }

 vtkSmartPointer<vtkDataArray> b;
 if (mesh_b_ptr == nullptr)
 {
 if (data_array_a_name == data_array_b_name)
 {
 std::cerr << "Error: You are trying to compare data array `"
 << data_array_a_name
 << "' from first file to itself. Aborting.\n";
 std::exit(3);
 }

 b = getDataArray(*mesh_a_ptr, data_array_b_name, domain);
 }
 else
 {
 b = getDataArray(*mesh_b_ptr, data_array_b_name, domain);
 }

 if (!b)
 {
 return {false, nullptr, nullptr};
 }

 return {true, a, b};
}

bool compareCellTopology(vtkCellArray* const cells_a,
 vtkCellArray* const cells_b)
{
 vtkIdType const n_cells_a{cells_a->GetNumberOfCells()};
 vtkIdType const n_cells_b{cells_b->GetNumberOfCells()};

 if (n_cells_a != n_cells_b)
 {
 std::cerr << "Number of cells in the first mesh is " << n_cells_a
 << " and differs from the number of cells in the second "
 "mesh, which is "
 << n_cells_b << "\n";
 return false;
 }

 vtkIdType n_cell_points_a, n_cell_points_b;
#if (VTK_MAJOR_VERSION > 8 || VTK_MINOR_VERSION == 90)
 const vtkIdType *cell_points_a, *cell_points_b;
#else
 vtkIdType *cell_points_a, *cell_points_b;
#endif
 cells_a->InitTraversal();
 cells_b->InitTraversal();
 int get_next_cell_a = cells_a->GetNextCell(n_cell_points_a, cell_points_a);
 int get_next_cell_b = cells_b->GetNextCell(n_cell_points_b, cell_points_b);
 int cell_number = 0;
 while (get_next_cell_a == 1 && get_next_cell_b == 1)
 {
 if (n_cell_points_a != n_cell_points_b)
 {
 std::cerr << "Cell " << cell_number << " in first input has "
 << n_cell_points_a << " points but in the second input "
 << n_cell_points_b << " points.\n";
 }

 for (vtkIdType i = 0; i < n_cell_points_a; ++i)
 {
 if (cell_points_a[i] != cell_points_b[i])
 {
 std::cerr << "Point " << i << " of cell " << cell_number
 << " has id " << cell_points_a[i]
 << " in the first input but id " << cell_points_b[i]
 << " in the second input.\n";
 return false;
 }
 }

 get_next_cell_a = cells_a->GetNextCell(n_cell_points_a, cell_points_a);
 get_next_cell_b = cells_b->GetNextCell(n_cell_points_b, cell_points_b);
 cell_number++;
 }

 if (get_next_cell_a != 0)
 {
 std::cerr << "Unexpected return value (" << get_next_cell_a
 << ") for cells_a->GetNextCell() call. Expected 0 for "
 "end-of-list or 1 for no error.\n";
 return false;
 }
 if (get_next_cell_b != 0)
 {
 std::cerr << "Unexpected return value (" << get_next_cell_b
 << ") for cells_b->GetNextCell() call. Expected 0 for "
 "end-of-list or 1 for no error.\n";
 return false;
 }

 return true;
}

bool comparePoints(vtkPoints* const points_a, vtkPoints* const points_b,
 double const eps_squared)
{
 vtkIdType const n_points_a{points_a->GetNumberOfPoints()};
 vtkIdType const n_points_b{points_b->GetNumberOfPoints()};

 if (n_points_a != n_points_b)
 {
 std::cerr << "Number of points in the first mesh is " << n_points_a
 << " and differs from the number of points in the second "
 "mesh, which is "
 << n_points_b << "\n";
 return false;
 }

 for (vtkIdType p = 0; p < n_points_a; ++p)
 {
 auto const a = points_a->GetPoint(p);
 auto const b = points_b->GetPoint(p);
 double const distance2 = vtkMath::Distance2BetweenPoints(a, b);
 if (distance2 >= eps_squared)
 {
 std::cerr << "Point " << p << " with coordinates (" << a[0] << ", "
 << a[1] << ", " << a[2]
 << ") from the first mesh is significantly different "
 "from the same point in the second mesh, which "
 "has coordinates ("
 << b[0] << ", " << b[1] << ", " << b[2]
 << ") with distance between them " << std::sqrt(distance2)
 << "\n";
 return false;
 }
 }
 return true;
}

int main(int argc, char* argv[])
{
 auto const digits10 = std::numeric_limits<double>::digits10;
 auto const args = parseCommandLine(argc, argv);

 // Setup the standard output and error stream numerical formats.
 std::cout << std::scientific << std::setprecision(digits10);
 std::cerr << std::scientific << std::setprecision(digits10);

 auto meshes = readMeshes(args.vtk_input_a, args.vtk_input_b);

 if (args.meshcheck)
 {
 if (args.vtk_input_b.empty())
 {
 std::cerr
 << "Error: Second input file must be provided for mesh check."\
 << "\n";
 return EXIT_FAILURE;
 }

 if (!std::get<0>(meshes) || !std::get<1>(meshes))
 {
 std::cerr << "Error reading mesh file(s)." << std::endl;
 return EXIT_FAILURE;
 }

 if (args.vtk_input_a == args.vtk_input_b)
 {
 std::cout << "Will not compare meshes from same input file.\n";
 return EXIT_SUCCESS;
 }

 if (!comparePoints(std::get<0>(meshes)->GetPoints(),
 std::get<1>(meshes)->GetPoints(),
 args.abs_err_thr * args.abs_err_thr))
 {
 std::cerr << "Error in mesh points' comparison occurred.\n";
 return EXIT_FAILURE;
 }

 if (!compareCellTopology(std::get<0>(meshes)->GetCells(),
 std::get<1>(meshes)->GetCells()))
 {
 std::cerr << "Error in cells' topology comparison occurred.\n";
 return EXIT_FAILURE;
 }
 return EXIT_SUCCESS;
 }

 // Read arrays from input file.
 bool read_successful;
 vtkSmartPointer<vtkDataArray> a;
 vtkSmartPointer<vtkDataArray> b;

 std::tie(read_successful, a, b) =
 readDataArraysFromMeshes(meshes, args.data_array_a, args.data_array_b);

 if (!read_successful)
 return EXIT_FAILURE;

 if (!args.quiet)
 std::cout << "Comparing data array `" << args.data_array_a
 << "' from file `" << args.vtk_input_a << "' to data array `"
 << args.data_array_b << "' from file `" << args.vtk_input_b
 << "'.\n";

 // Check similarity of the data arrays.

 // Is numeric
 if (!a->IsNumeric())
 {
 std::cerr << "Data in data array a is not numeric:\n"
 << "data type is " << a->GetDataTypeAsString() << "\n";

 return EXIT_FAILURE;
 }
 if (!b->IsNumeric())
 {
 std::cerr << "Data in data array b is not numeric.\n"
 << "data type is " << b->GetDataTypeAsString() << "\n";
 return EXIT_FAILURE;
 }

 auto const num_tuples = a->GetNumberOfTuples();
 // Number of components
 if (num_tuples != b->GetNumberOfTuples())
 {
 std::cerr << "Number of tuples differ:\n"
 << num_tuples << " in data array a and "
 << b->GetNumberOfTuples() << " in data array b\n";
 return EXIT_FAILURE;
 }

 auto const num_components = a->GetNumberOfComponents();
 // Number of components
 if (num_components != b->GetNumberOfComponents())
 {
 std::cerr << "Number of components differ:\n"
 << num_components << " in data array a and "
 << b->GetNumberOfComponents() << " in data array b\n";
 return EXIT_FAILURE;
 }

 // Calculate difference of the data arrays.

 // Absolute error and norms.
 std::vector<double> abs_err_norm_l1(num_components);
 std::vector<double> abs_err_norm_2_2(num_components);
 std::vector<double> abs_err_norm_max(num_components);

 // Relative error and norms.
 std::vector<double> rel_err_norm_l1(num_components);
 std::vector<double> rel_err_norm_2_2(num_components);
 std::vector<double> rel_err_norm_max(num_components);

 for (auto tuple_idx = 0; tuple_idx < num_tuples; ++tuple_idx)
 {
 for (auto component_idx = 0; component_idx < num_components;
 ++component_idx)
 {
 auto const a_comp = a->GetComponent(tuple_idx, component_idx);
 auto const b_comp = b->GetComponent(tuple_idx, component_idx);
 auto const abs_err = std::abs(a_comp - b_comp);

 abs_err_norm_l1[component_idx] += abs_err;
 abs_err_norm_2_2[component_idx] += abs_err * abs_err;
 abs_err_norm_max[component_idx] =
 std::max(abs_err_norm_max[component_idx], abs_err);

 // relative error and its norms:
 double rel_err;

 if (abs_err == 0.0)
 {
 rel_err = 0.0;
 }
 else if (a_comp == 0.0 || b_comp == 0.0)
 {
 rel_err = std::numeric_limits<double>::infinity();
 }
 else
 {
 rel_err =
 abs_err / std::min(std::abs(a_comp), std::abs(b_comp));
 }

 rel_err_norm_l1[component_idx] += rel_err;
 rel_err_norm_2_2[component_idx] += rel_err * rel_err;
 rel_err_norm_max[component_idx] =
 std::max(rel_err_norm_max[component_idx], rel_err);

 if (abs_err > args.abs_err_thr && rel_err > args.rel_err_thr &&
 args.verbose)
 {
 std::cout << "tuple: " << std::setw(4) << tuple_idx
 << "component: " << std::setw(2) << component_idx
 << ": abs err = " << std::setw(digits10 + 7)
 << abs_err
 << ", rel err = " << std::setw(digits10 + 7)
 << rel_err << "\n";
 }
 }
 }

 // Error information
 if (!args.quiet)
 {
 std::cout << "Computed difference between data arrays:\n";
 std::cout << "abs l1 norm = " << abs_err_norm_l1 << "\n";
 std::cout << "abs l2-norm^2 = " << abs_err_norm_2_2 << "\n";

 // temporary squared norm vector for output.
 std::vector<double> abs_err_norm_2;
 std::transform(std::begin(abs_err_norm_2_2), std::end(abs_err_norm_2_2),
 std::back_inserter(abs_err_norm_2),
 [](double x) { return std::sqrt(x); });
 std::cout << "abs l2-norm = " << abs_err_norm_2 << "\n";

 std::cout << "abs maximum norm = " << abs_err_norm_max << "\n";
 std::cout << "\n";

 std::cout << "rel l1 norm = " << rel_err_norm_l1 << "\n";
 std::cout << "rel l2-norm^2 = " << rel_err_norm_2_2 << "\n";

 // temporary squared norm vector for output.
 std::vector<double> rel_err_norm_2;
 std::transform(std::begin(rel_err_norm_2_2), std::end(rel_err_norm_2_2),
 std::back_inserter(rel_err_norm_2),
 [](double x) { return std::sqrt(x); });
 std::cout << "rel l2-norm = " << rel_err_norm_2 << "\n";

 std::cout << "rel maximum norm = " << rel_err_norm_max << "\n";
 }

 if (*std::max_element(abs_err_norm_max.begin(), abs_err_norm_max.end()) >
 args.abs_err_thr &&
 *std::max_element(rel_err_norm_max.begin(), rel_err_norm_max.end()) >
 args.rel_err_thr)
 {
 if (!args.quiet)
 std::cout << "Absolute and relative error (maximum norm) are larger"
 " than the corresponding thresholds "
 << args.abs_err_thr << " and " << args.rel_err_thr
 << ".\n";
 return EXIT_FAILURE;
 }

 return EXIT_SUCCESS;
}