| 761 | //---------------------------------------------------------------------------- |
| 762 | template <std::floating_point U> |
| 763 | void io::VTKFile::write(const mesh::Mesh<U>& mesh, double time) |
| 764 | { |
| 765 | if (!_pvd_xml) |
| 766 | throw std::runtime_error("VTKFile has already been closed"); |
| 767 | |
| 768 | // Get the PVD "Collection" node |
| 769 | pugi::xml_node xml_collections |
| 770 | = _pvd_xml->child("VTKFile").child("Collection"); |
| 771 | assert(xml_collections); |
| 772 | |
| 773 | // Compute counter string |
| 774 | const std::string counter_str = get_counter(xml_collections, "DataSet"); |
| 775 | |
| 776 | // Get mesh data for this rank |
| 777 | auto topology = mesh.topology(); |
| 778 | assert(topology); |
| 779 | const mesh::Geometry<U>& geometry = mesh.geometry(); |
| 780 | auto xmap = geometry.index_map(); |
| 781 | assert(xmap); |
| 782 | const int tdim = topology->dim(); |
| 783 | const std::int32_t num_points = xmap->size_local() + xmap->num_ghosts(); |
| 784 | const std::int32_t num_cells = topology->index_map(tdim)->size_local() |
| 785 | + topology->index_map(tdim)->num_ghosts(); |
| 786 | |
| 787 | // Create a VTU XML object |
| 788 | pugi::xml_document xml_vtu; |
| 789 | pugi::xml_node vtk_node_vtu = xml_vtu.append_child("VTKFile"); |
| 790 | vtk_node_vtu.append_attribute("type") = "UnstructuredGrid"; |
| 791 | vtk_node_vtu.append_attribute("version") = "2.2"; |
| 792 | pugi::xml_node grid_node_vtu = vtk_node_vtu.append_child("UnstructuredGrid"); |
| 793 | |
| 794 | // Add "Piece" node and required metadata |
| 795 | pugi::xml_node piece_node = grid_node_vtu.append_child("Piece"); |
| 796 | piece_node.append_attribute("NumberOfPoints") = num_points; |
| 797 | piece_node.append_attribute("NumberOfCells") = num_cells; |
| 798 | |
| 799 | mesh::CellType cell_type = topology->cell_type(); |
| 800 | |
| 801 | // Add mesh data to "Piece" node |
| 802 | const auto [cells, cshape] |
| 803 | = extract_vtk_connectivity(mesh.geometry().dofmaps().front(), cell_type); |
| 804 | std::array<std::size_t, 2> xshape = {geometry.x().size() / 3, 3}; |
| 805 | std::vector<std::uint8_t> x_ghost(xshape[0], 0); |
| 806 | std::fill(std::next(x_ghost.begin(), xmap->size_local()), x_ghost.end(), 1); |
| 807 | add_mesh(geometry.x(), xshape, geometry.input_global_indices(), x_ghost, |
| 808 | cells, cshape, *topology->index_map(tdim), cell_type, |
| 809 | topology->dim(), piece_node); |
| 810 | |
| 811 | // Create filepath for a .vtu file |
| 812 | auto create_vtu_path = [file_root = _filename.parent_path(), |
| 813 | file_name = _filename.stem(), counter_str](int rank) |
| 814 | { |
| 815 | std::filesystem::path vtu = file_root / file_name; |
| 816 | vtu += +"_p" + std::to_string(rank) + "_" + counter_str; |
| 817 | vtu.replace_extension("vtu"); |
| 818 | return vtu; |
| 819 | }; |
| 820 | |