diff --git a/cpp/dolfinx/io/ADIOS2Writers.h b/cpp/dolfinx/io/ADIOS2Writers.h index ab4030c1885..44d16164f56 100644 --- a/cpp/dolfinx/io/ADIOS2Writers.h +++ b/cpp/dolfinx/io/ADIOS2Writers.h @@ -638,29 +638,33 @@ std::stringstream create_vtk_schema(const std::vector& point_data, /// Extract name of functions and split into real and imaginary component template -std::vector +std::tuple,std::vector> extract_function_names(const typename adios2_writer::U& u) { - std::vector names; + std::vector names, dg0_names; for (auto& v : u) { std::visit( - [&names](auto&& u) + [&names,&dg0_names](auto&& u) { using U = std::decay_t; using X = typename U::element_type; + std::vector* fnames = &names; + if (impl::is_cellwise(*(u->function_space()->element()))) { + fnames = &dg0_names; + } if constexpr (std::is_floating_point_v) - names.push_back(u->name); + fnames->push_back(u->name); else { - names.push_back(u->name + impl_adios2::field_ext[0]); - names.push_back(u->name + impl_adios2::field_ext[1]); + fnames->push_back(u->name + impl_adios2::field_ext[0]); + fnames->push_back(u->name + impl_adios2::field_ext[1]); } }, v); } - return names; + return {names, dg0_names}; } /// Given a Function, write the coefficient to file using ADIOS2. @@ -905,7 +909,7 @@ class VTXWriter : public ADIOS2Writer std::shared_ptr> mesh, std::string engine = "BPFile") : ADIOS2Writer(comm, filename, "VTX mesh writer", engine), _mesh(mesh), - _mesh_reuse_policy(VTXMeshPolicy::update), _is_piecewise_constant(false) + _mesh_reuse_policy(VTXMeshPolicy::update), _has_piecewise_constant(false) { // Define VTK scheme attribute for mesh std::string vtk_scheme = impl_vtx::create_vtk_schema({}, {}).str(); @@ -932,7 +936,7 @@ class VTXWriter : public ADIOS2Writer VTXMeshPolicy mesh_policy = VTXMeshPolicy::update) : ADIOS2Writer(comm, filename, "VTX function writer", engine), _mesh(impl_adios2::extract_common_mesh(u)), - _mesh_reuse_policy(mesh_policy), _u(u), _is_piecewise_constant(false) + _mesh_reuse_policy(mesh_policy), _u(u), _has_piecewise_constant(false) { if (u.empty()) throw std::runtime_error("VTXWriter fem::Function list is empty."); @@ -941,6 +945,23 @@ class VTXWriter : public ADIOS2Writer auto V0 = std::visit([](auto& u) { return u->function_space().get(); }, u.front()); assert(V0); + bool has_V0_changed = false; + for (auto& v : u) + { + std::visit( + [&V0,&has_V0_changed](auto& u) + { + auto V = u->function_space().get(); + assert(V); + if (!impl::is_cellwise(*V->element())) + { + V0 = V; + has_V0_changed = true; + } + }, + v); + if (has_V0_changed) break; + } auto element0 = V0->element().get(); assert(element0); @@ -961,19 +982,17 @@ class VTXWriter : public ADIOS2Writer "supported. Interpolate Functions before output."); } - // Check if function is DG 0 - if (element0->space_dimension() / element0->block_size() == 1) - _is_piecewise_constant = true; - // Check that all functions come from same element type for (auto& v : _u) { std::visit( - [V0](auto& u) + [V0,this](auto& u) { auto element = u->function_space()->element(); assert(element); - if (*element != *V0->element().get()) + bool is_piecewise_constant = impl::is_cellwise(*element); + _has_piecewise_constant = _has_piecewise_constant || is_piecewise_constant; + if (*element != *V0->element().get() and !is_piecewise_constant) { throw std::runtime_error("All functions in VTXWriter must have " "the same element type."); @@ -981,10 +1000,11 @@ class VTXWriter : public ADIOS2Writer #ifndef NDEBUG auto dmap0 = V0->dofmap()->map(); auto dmap = u->function_space()->dofmap()->map(); - if (dmap0.size() != dmap.size() + if ((dmap0.size() != dmap.size() or !std::equal(dmap0.data_handle(), dmap0.data_handle() + dmap0.size(), dmap.data_handle())) + and !is_piecewise_constant) { throw std::runtime_error( "All functions must have the same dofmap for VTXWriter."); @@ -995,12 +1015,9 @@ class VTXWriter : public ADIOS2Writer } // Define VTK scheme attribute for set of functions - std::vector names = impl_vtx::extract_function_names(u); + auto [names, dg0_names] = impl_vtx::extract_function_names(u); std::string vtk_scheme; - if (_is_piecewise_constant) - vtk_scheme = impl_vtx::create_vtk_schema({}, names).str(); - else - vtk_scheme = impl_vtx::create_vtk_schema(names, {}).str(); + vtk_scheme = impl_vtx::create_vtk_schema(names, dg0_names).str(); impl_adios2::define_attribute(*_io, "vtk.xml", vtk_scheme); } @@ -1054,17 +1071,9 @@ class VTXWriter : public ADIOS2Writer _engine->template Put(var_step, t); // If we have no functions or DG functions write the mesh to file - if (_is_piecewise_constant or _u.empty()) + if (_has_piecewise_constant or _u.empty()) { impl_vtx::vtx_write_mesh(*_io, *_engine, *_mesh); - if (_is_piecewise_constant) - { - for (auto& v : _u) - { - std::visit([&](auto& u) - { impl_vtx::vtx_write_data(*_io, *_engine, *u); }, v); - } - } } else { @@ -1093,12 +1102,12 @@ class VTXWriter : public ADIOS2Writer _engine->PerformPuts(); } - // Write function data for each function to file - for (auto& v : _u) - { - std::visit([&](auto& u) - { impl_vtx::vtx_write_data(*_io, *_engine, *u); }, v); - } + } + // Write function data for each function to file + for (auto& v : _u) + { + std::visit([&](auto& u) + { impl_vtx::vtx_write_data(*_io, *_engine, *u); }, v); } _engine->EndStep(); @@ -1115,7 +1124,7 @@ class VTXWriter : public ADIOS2Writer std::vector _x_ghost; // Special handling of piecewise constant functions - bool _is_piecewise_constant; + bool _has_piecewise_constant; }; /// Type deduction diff --git a/cpp/dolfinx/io/VTKFile.cpp b/cpp/dolfinx/io/VTKFile.cpp index 0f9b7d59b14..39d2a401fea 100644 --- a/cpp/dolfinx/io/VTKFile.cpp +++ b/cpp/dolfinx/io/VTKFile.cpp @@ -27,6 +27,7 @@ #include using namespace dolfinx; +namespace impl = dolfinx::io::impl; namespace { @@ -34,15 +35,6 @@ namespace /// field constexpr std::array field_ext = {"_real", "_imag"}; -/// Return true if element is a cell-wise constant, otherwise false -template -bool is_cellwise(const fem::FiniteElement& e) -{ - return e.space_dimension() / e.block_size() == 1; -} - -//---------------------------------------------------------------------------- - /// Get counter string to include in filename std::string get_counter(const pugi::xml_node& node, const std::string& name) { @@ -338,7 +330,7 @@ void write_function( { auto V = v.get().function_space(); assert(V); - if (!is_cellwise(*V->element())) + if (!impl::is_cellwise(*V->element())) { V0 = V; break; @@ -377,7 +369,7 @@ void write_function( } // Check that pointwise elements are the same (up to the block size) - if (!is_cellwise(*e)) + if (!impl::is_cellwise(*e)) { if (*e != *element0) { @@ -412,7 +404,7 @@ void write_function( std::vector x_ghost; std::vector cells; std::array cshape; - if (is_cellwise(*V0->element())) + if (impl::is_cellwise(*V0->element())) { std::vector tmp; std::tie(tmp, cshape) = io::extract_vtk_connectivity( @@ -455,7 +447,7 @@ void write_function( assert(_u.get().function_space()); auto e = _u.get().function_space()->element(); assert(e); - auto data_type = is_cellwise(*e) ? "CellData" : "PointData"; + auto data_type = impl::is_cellwise(*e) ? "CellData" : "PointData"; if (piece_node.child(data_type).empty()) piece_node.append_child(data_type); @@ -487,7 +479,7 @@ void write_function( if (rank > 0) component_vector[0] = num_components; - if (is_cellwise(*e)) + if (impl::is_cellwise(*e)) { // -- Cell-wise data @@ -633,7 +625,7 @@ void write_function( grid_node.append_attribute("GhostLevel") = 1; for (auto _u : u) { - if (auto e = _u.get().function_space()->element(); is_cellwise(*e)) + if (auto e = _u.get().function_space()->element(); impl::is_cellwise(*e)) { if (grid_node.child("PCellData").empty()) grid_node.append_child("PCellData"); @@ -655,7 +647,7 @@ void write_function( assert(V); auto e = V->element(); assert(e); - std::string d_type = is_cellwise(*e) ? "PCellData" : "PPointData"; + std::string d_type = impl::is_cellwise(*e) ? "PCellData" : "PPointData"; pugi::xml_node data_pnode = grid_node.child(d_type.c_str()); // Pad to 3D if vector/tensor is product of dimensions is smaller than diff --git a/cpp/dolfinx/io/vtk_utils.h b/cpp/dolfinx/io/vtk_utils.h index 5bcb7d368da..cd296c5f433 100644 --- a/cpp/dolfinx/io/vtk_utils.h +++ b/cpp/dolfinx/io/vtk_utils.h @@ -169,6 +169,14 @@ tabulate_lagrange_dof_coordinates(const fem::FunctionSpace& V) return {std::move(coords), cshape, std::move(x_id), std::move(id_ghost)}; } +/// Return true if element is a cell-wise constant, otherwise false +/// This could return a constexpr +template +bool is_cellwise(const fem::FiniteElement& e) +{ + return e.space_dimension() / e.block_size() == 1; +} + } // namespace impl /// @brief Given a FunctionSpace, create a topology and geometry based