Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support list of DG-0 + other space in VTXWriter #3375

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Support list of DG-0 + other space in VTXWriter #3375

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
83 changes: 46 additions & 37 deletions cpp/dolfinx/io/ADIOS2Writers.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -638,29 +638,33 @@ std::stringstream create_vtk_schema(const std::vector<std::string>& point_data,

/// Extract name of functions and split into real and imaginary component
template <std::floating_point T>
std::vector<std::string>
std::tuple<std::vector<std::string>,std::vector<std::string>>
extract_function_names(const typename adios2_writer::U<T>& u)
{
std::vector<std::string> names;
std::vector<std::string> names, dg0_names;
for (auto& v : u)
{
std::visit(
[&names](auto&& u)
[&names,&dg0_names](auto&& u)
{
using U = std::decay_t<decltype(u)>;
using X = typename U::element_type;
std::vector<std::string>* fnames = &names;
if (impl::is_cellwise(*(u->function_space()->element()))) {
fnames = &dg0_names;
}
if constexpr (std::is_floating_point_v<typename X::value_type>)
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.
Expand Down Expand Up @@ -905,7 +909,7 @@ class VTXWriter : public ADIOS2Writer
std::shared_ptr<const mesh::Mesh<T>> 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();
Expand All @@ -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<T>(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.");
Expand All @@ -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);

Expand All @@ -961,30 +982,29 @@ 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.");
}
#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.");
Expand All @@ -995,12 +1015,9 @@ class VTXWriter : public ADIOS2Writer
}

// Define VTK scheme attribute for set of functions
std::vector<std::string> names = impl_vtx::extract_function_names<T>(u);
auto [names, dg0_names] = impl_vtx::extract_function_names<T>(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<std::string>(*_io, "vtk.xml", vtk_scheme);
}
Expand Down Expand Up @@ -1054,17 +1071,9 @@ class VTXWriter : public ADIOS2Writer
_engine->template Put<double>(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
{
Expand Down Expand Up @@ -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();
Expand All @@ -1115,7 +1124,7 @@ class VTXWriter : public ADIOS2Writer
std::vector<std::uint8_t> _x_ghost;

// Special handling of piecewise constant functions
bool _is_piecewise_constant;
bool _has_piecewise_constant;
};

/// Type deduction
Expand Down
24 changes: 8 additions & 16 deletions cpp/dolfinx/io/VTKFile.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,22 +27,14 @@
#include <string>

using namespace dolfinx;
namespace impl = dolfinx::io::impl;

namespace
{
/// String suffix for real and complex components of a vector-valued
/// field
constexpr std::array field_ext = {"_real", "_imag"};

/// Return true if element is a cell-wise constant, otherwise false
template <std::floating_point T>
bool is_cellwise(const fem::FiniteElement<T>& 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)
{
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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)
{
Expand Down Expand Up @@ -412,7 +404,7 @@ void write_function(
std::vector<std::uint8_t> x_ghost;
std::vector<std::int64_t> cells;
std::array<std::size_t, 2> cshape;
if (is_cellwise(*V0->element()))
if (impl::is_cellwise(*V0->element()))
{
std::vector<std::int64_t> tmp;
std::tie(tmp, cshape) = io::extract_vtk_connectivity(
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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");
Expand All @@ -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
Expand Down
8 changes: 8 additions & 0 deletions cpp/dolfinx/io/vtk_utils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,6 +169,14 @@ tabulate_lagrange_dof_coordinates(const fem::FunctionSpace<T>& 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 <std::floating_point T>
bool is_cellwise(const fem::FiniteElement<T>& e)
{
return e.space_dimension() / e.block_size() == 1;
}

} // namespace impl

/// @brief Given a FunctionSpace, create a topology and geometry based
Expand Down