From: Luca Heltai Date: Thu, 28 Jan 2021 11:32:49 +0000 (+0100) Subject: Added GridOut::write_msh() X-Git-Tag: v9.3.0-rc1~540^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fadac6a42db523cafef38b61a77c14f246c8aae5;p=dealii.git Added GridOut::write_msh() --- diff --git a/include/deal.II/grid/grid_out.h b/include/deal.II/grid/grid_out.h index 6ae4408c9d..09c9478178 100644 --- a/include/deal.II/grid/grid_out.h +++ b/include/deal.II/grid/grid_out.h @@ -1098,6 +1098,53 @@ public: void write_msh(const Triangulation &tria, std::ostream &out) const; +#ifdef DEAL_II_GMSH_WITH_API + /** + * Write the triangulation in any format supported by gmsh API. + * + * Gmsh API allows writing its output in several formats through their C++ + * API. This function translates a Triangulation object into a gmsh collection + * of entities and calls the gmsh::write() method with the filename passed as + * argument. This method generates a different entity for each unique pair of + * non default manifold id and boundary id, and writes a gmsh physical group + * for each unique combination, allowing you to read back the triangulation + * using the GridIn::read_msh() method that takes a string as argument. + * + * In particular, all cell objects with non default boundary id or non + * default manifold id are grouped in a unique physical tag, whose name + * contains the boundary and manifold indicators. The names are constructed + * using Patterns::Tools::to_value() applied to a `std::map` + * where the keys are either `MaterialID`, `BoundaryID`, or `ManifoldID`, + * i.e., a cell with material id 1, and manifold id 3 would be grouped in a + * physical tag (whose numbering is not specified), named + * `MaterialID:1, ManifoldID:3`. + * + * For example, calling the method with a hyper ball grid refined once, + * would results in the following physical tags defined in the output file: + * @code + * MeshFormat + * 4.1 0 8 + * $EndMeshFormat + * $PhysicalNames + * 3 + * 1 1 "ManifoldID:0" + * 1 2 "BoundaryID:-1, ManifoldID:1" + * 2 3 "ManifoldID:1" + * $EndPhysicalNames + * $Entities + * ... + * @endcode + * + * The special boundary id `-1` is used to indicate internal boundaries. The + * internal boundaries must be specified whenever it is necessary to specify + * a non-flat manifold id. + */ + template + void + write_msh(const Triangulation &tria, + const std::string & filename) const; +#endif + /** * Write the triangulation in the ucd format. * diff --git a/source/grid/grid_out.cc b/source/grid/grid_out.cc index 3d39cc3952..74d95b4812 100644 --- a/source/grid/grid_out.cc +++ b/source/grid/grid_out.cc @@ -32,6 +32,10 @@ #include #include +#ifdef DEAL_II_GMSH_WITH_API +# include +#endif + #include #include #include @@ -1470,6 +1474,211 @@ GridOut::write_xfig(const Triangulation<2> &tria, +#ifdef DEAL_II_GMSH_WITH_API +template +void +GridOut::write_msh(const Triangulation &tria, + const std::string & filename) const +{ + // mesh Type renumbering + const std::array dealii_to_gmsh_type = {{15, 1, 2, 3, 4, 7, 6, 5}}; + + // Vertex renumbering, by dealii type + const std::vector> dealii_to_gmsh = { + {{{0}}, + {{0, 1}}, + {{0, 1, 2}}, + {{0, 1, 3, 2}}, + {{0, 1, 2, 3}}, + {{0, 1, 3, 2, 4}}, + {{0, 1, 2, 3, 4, 5}}, + {{0, 1, 3, 2, 4, 5, 7, 6}}}}; + + // Extract all vertices (nodes in gmsh terminology), and store their three + // dimensional coordinates (regardless of dim). + const auto & vertices = tria.get_vertices(); + std::vector coords(3 * vertices.size()); + std::vector nodes(vertices.size()); + + // Each node has a strictly positive tag. We assign simply its index+1. + std::size_t i = 0; + for (const auto &p : vertices) + { + for (unsigned int d = 0; d < spacedim; ++d) + coords[i * 3 + d] = p[d]; + nodes[i] = i + 1; + ++i; + } + + // Construct one entity tag per boundary and manifold id pair. + // We need to be smart here, in order to save some disk space. All cells need + // to be written, but only faces and lines that have non default boundary ids + // and/or manifold ids. We collect them into pairs, and for each unique pair, + // we create a gmsh entity where we store the elements. Pre-count all the + // entities, and make sure we know which pair refers to what entity and + // vice-versa. + using IdPair = std::pair; + std::map id_pair_to_entity_tag; + std::vector all_pairs; + { + std::set set_of_pairs; + for (const auto &cell : tria.active_cell_iterators()) + { + set_of_pairs.insert({cell->material_id(), cell->manifold_id()}); + for (const auto &f : cell->face_iterators()) + if (f->manifold_id() != numbers::flat_manifold_id || + (f->boundary_id() != 0 && + f->boundary_id() != numbers::internal_face_boundary_id)) + set_of_pairs.insert({f->boundary_id(), f->manifold_id()}); + if (dim > 2) + for (const auto &l : cell->line_indices()) + { + const auto &f = cell->line(l); + if (f->manifold_id() != numbers::flat_manifold_id || + (f->boundary_id() != 0 && + f->boundary_id() != numbers::internal_face_boundary_id)) + set_of_pairs.insert({f->boundary_id(), f->manifold_id()}); + } + } + all_pairs = {set_of_pairs.begin(), set_of_pairs.end()}; + + int entity = 1; + for (const auto &p : set_of_pairs) + id_pair_to_entity_tag[p] = entity++; + } + + const auto n_entity_tags = id_pair_to_entity_tag.size(); + + // All elements in the mesh, by entity tag, and by dealii type. + std::vector>> element_ids( + n_entity_tags, std::vector>(8)); + std::vector>> element_nodes( + n_entity_tags, std::vector>(8)); + + // One elment id counter for all dimensions. + std::size_t element_id = 1; + + const auto add_element = [&](const auto &element, const int &entity_tag) { + const auto type = element->reference_cell_type(); + + Assert(entity_tag > 0, ExcInternalError()); + // Add all vertex ids. Make sure we renumber to gmsh, and we add 1 to the + // global index. + for (const auto &v : element->vertex_indices()) + element_nodes[entity_tag - 1][type].emplace_back( + element->vertex_index(dealii_to_gmsh[type][v]) + 1); + + // Save the element id. + element_ids[entity_tag - 1][type].emplace_back(element_id); + ++element_id; + }; + + // Will create a separate gmsh entity, only if it's a cell, or if the + // boundary and/or the manifold ids are not the default ones. + // In the meanwhile, also store each pair of dimension and entity tag that was + // requested. + std::set> dim_entity_tag; + + auto maybe_add_element = + [&](const auto & element, + const types::boundary_id &boundary_or_material_id) { + const auto struct_dim = element->structure_dimension; + const auto manifold_id = element->manifold_id(); + + // Exclude default boundary/manifold id or invalid/flag + const bool non_default_boundary_or_material_id = + (boundary_or_material_id != 0 && + boundary_or_material_id != numbers::internal_face_boundary_id); + const bool non_default_manifold = + manifold_id != numbers::flat_manifold_id; + if (struct_dim == dim || non_default_boundary_or_material_id || + non_default_manifold) + { + const auto entity_tag = + id_pair_to_entity_tag[{boundary_or_material_id, manifold_id}]; + add_element(element, entity_tag); + dim_entity_tag.insert({struct_dim, entity_tag}); + } + }; + + // Loop recursively over all cells, faces, and possibly lines. + for (const auto &cell : tria.active_cell_iterators()) + { + maybe_add_element(cell, cell->material_id()); + for (const auto &face : cell->face_iterators()) + maybe_add_element(face, face->boundary_id()); + if (dim > 2) + for (const auto &l : cell->line_indices()) + maybe_add_element(cell->line(l), cell->line(l)->boundary_id()); + } + + // Now that we collected everything, plug them into gmsh + gmsh::initialize(); + gmsh::option::setNumber("General.Verbosity", 0); + gmsh::model::add("Grid generated in deal.II"); + for (const auto p : dim_entity_tag) + { + gmsh::model::addDiscreteEntity(p.first, p.second); + gmsh::model::mesh::addNodes(p.first, p.second, nodes, coords); + } + + for (unsigned int entity_tag = 0; entity_tag < n_entity_tags; ++entity_tag) + for (unsigned int t = 1; t < 8; ++t) + { + const auto all_element_ids = element_ids[entity_tag][t]; + const auto all_element_nodes = element_nodes[entity_tag][t]; + const auto gmsh_t = dealii_to_gmsh_type[t]; + if (all_element_ids.size() > 0) + gmsh::model::mesh::addElementsByType(entity_tag + 1, + gmsh_t, + all_element_ids, + all_element_nodes); + } + + // Make sure nodes belong to the right entities. + gmsh::model::mesh::reclassifyNodes(); + gmsh::model::mesh::removeDuplicateNodes(); + + // Now for each individual pair of dim and entry, add a physical group, if + // necessary + for (const auto &it : dim_entity_tag) + { + const auto &d = it.first; + const auto &entity_tag = it.second; + const auto &boundary_id = all_pairs[entity_tag - 1].first; + const auto &manifold_id = all_pairs[entity_tag - 1].second; + + std::string physical_name; + if (d == dim && boundary_id != 0) + physical_name += "MaterialID:" + Utilities::int_to_string( + static_cast(boundary_id)); + else if (d < dim && boundary_id != 0) + physical_name += + "BoundaryID:" + + (boundary_id == numbers::internal_face_boundary_id ? + "-1" : + Utilities::int_to_string(static_cast(boundary_id))); + + std::string sep = physical_name != "" ? ", " : ""; + if (manifold_id != numbers::flat_manifold_id) + physical_name += + sep + "ManifoldID:" + + Utilities::int_to_string(static_cast(manifold_id)); + const auto physical_tag = + gmsh::model::addPhysicalGroup(d, {entity_tag}, -1); + if (physical_name != "") + gmsh::model::setPhysicalName(d, physical_tag, physical_name); + } + + + gmsh::write(filename); + gmsh::clear(); + gmsh::finalize(); +} +#endif + + + namespace { /** diff --git a/source/grid/grid_out.inst.in b/source/grid/grid_out.inst.in index dc59139e57..68baf3875a 100644 --- a/source/grid/grid_out.inst.in +++ b/source/grid/grid_out.inst.in @@ -26,6 +26,10 @@ for (deal_II_dimension : DIMENSIONS) template void GridOut::write_msh(const Triangulation &, std::ostream &) const; +#ifdef DEAL_II_GMSH_WITH_API + template void GridOut::write_msh(const Triangulation &, + const std::string &) const; +#endif #if deal_II_dimension != 2 template void GridOut::write_xfig( @@ -83,6 +87,11 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) template void GridOut::write_msh( const Triangulation &, std::ostream &) const; +# ifdef DEAL_II_GMSH_WITH_API + template void GridOut::write_msh( + const Triangulation &, + const std::string &) const; +# endif template void GridOut::write_ucd( const Triangulation &, std::ostream &) const; diff --git a/tests/gmsh/gmsh_api_01.cc b/tests/gmsh/gmsh_api_01.cc new file mode 100644 index 0000000000..e092c79965 --- /dev/null +++ b/tests/gmsh/gmsh_api_01.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Create a hyper ball, and use gmsh to output it. + +#include +#include +#include + +#include "../tests.h" + +int +main() +{ + initlog(); + + const unsigned int dim = 2; + const unsigned int spacedim = 2; + + Triangulation tria; + GridGenerator::hyper_ball(tria); + tria.refine_global(1); + + GridOut go; + go.write_msh(tria, "output.msh"); + + cat_file("output.msh"); +} diff --git a/tests/gmsh/gmsh_api_01.with_gmsh_with_api=on.output b/tests/gmsh/gmsh_api_01.with_gmsh_with_api=on.output new file mode 100644 index 0000000000..3d0a3e1906 --- /dev/null +++ b/tests/gmsh/gmsh_api_01.with_gmsh_with_api=on.output @@ -0,0 +1,142 @@ + +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +3 +1 1 "ManifoldID:0" +1 2 "BoundaryID:-1, ManifoldID:1" +2 3 "ManifoldID:1" +$EndPhysicalNames +$Entities +0 2 2 0 +1 -0.9999999999999999 -0.9999999999999999 0 0.9999999999999999 0.9999999999999999 0 1 1 0 +4 -0.9999999999999999 -0.9999999999999999 0 0.9999999999999999 0.9999999999999999 0 1 2 0 +2 -0.9999999999999999 -0.9999999999999999 0 0.9999999999999999 0.9999999999999999 0 1 3 0 +3 -0.2928932188134525 -0.2928932188134525 0 0.2928932188134525 0.2928932188134525 0 1 4 0 +$EndEntities +$Nodes +4 25 1 25 +1 1 0 8 +1 +2 +7 +8 +9 +11 +13 +19 +-0.7071067811865475 -0.7071067811865475 0 +0.7071067811865475 -0.7071067811865475 0 +-0.7071067811865475 0.7071067811865475 0 +0.7071067811865475 0.7071067811865475 0 +-1.110223024625156e-16 -0.9999999999999999 0 +-0.9999999999999999 -1.110223024625156e-16 0 +0.9999999999999999 -1.110223024625156e-16 0 +-1.110223024625156e-16 0.9999999999999999 0 +1 4 0 12 +10 +12 +14 +15 +16 +17 +18 +20 +21 +22 +24 +25 +-0.5 -0.5 0 +0.5 -0.5 0 +0 -0.2928932188134525 0 +-0.2928932188134525 0 0 +0.2928932188134525 0 0 +0 0.2928932188134525 0 +-0.5 0.5 0 +0.5 0.5 0 +-5.551115123125782e-17 -0.6464466094067262 0 +-0.6464466094067262 -5.551115123125783e-17 0 +0.6464466094067262 -5.551115123125782e-17 0 +-5.551115123125783e-17 0.6464466094067262 0 +2 2 0 4 +3 +4 +5 +6 +-0.2928932188134525 -0.2928932188134525 0 +0.2928932188134525 -0.2928932188134525 0 +-0.2928932188134525 0.2928932188134525 0 +0.2928932188134525 0.2928932188134525 0 +2 3 0 1 +23 +0 0 0 +$EndNodes +$Elements +4 60 1 60 +1 1 1 8 +3 1 9 +7 9 2 +16 1 11 +23 11 7 +35 2 13 +39 13 8 +48 7 19 +55 19 8 +1 4 1 32 +2 9 21 +4 10 21 +6 9 21 +8 21 12 +10 21 14 +11 10 21 +13 21 14 +14 21 12 +17 10 22 +18 11 22 +20 10 22 +21 22 15 +24 22 18 +25 11 22 +27 22 18 +28 22 15 +34 13 24 +36 12 24 +38 13 24 +40 24 20 +42 24 16 +43 12 24 +45 24 16 +46 24 20 +49 18 25 +50 19 25 +52 18 25 +53 25 17 +56 25 20 +57 19 25 +59 25 20 +60 25 17 +2 2 3 16 +1 1 9 21 10 +5 9 2 12 21 +9 10 21 14 3 +12 21 12 4 14 +15 1 10 22 11 +19 10 3 15 22 +22 11 22 18 7 +26 22 15 5 18 +33 2 13 24 12 +37 13 8 20 24 +41 12 24 16 4 +44 24 20 6 16 +47 7 18 25 19 +51 18 5 17 25 +54 19 25 20 8 +58 25 17 6 20 +2 3 3 4 +29 3 14 23 15 +30 14 4 16 23 +31 15 23 17 5 +32 23 16 6 17 +$EndElements + diff --git a/tests/gmsh/gmsh_api_02.cc b/tests/gmsh/gmsh_api_02.cc new file mode 100644 index 0000000000..7ab73247a4 --- /dev/null +++ b/tests/gmsh/gmsh_api_02.cc @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Create all reference cells, and output them with gmsh. + +#include +#include + +#include "../tests.h" + +template +void +test(const std::uint8_t kind, const std::string out = "") +{ + deallog << "Testing kind(" << (int)kind << ") in dimensions " + << "<" << dim << "," << spacedim << ">" << std::endl; + + Triangulation tria; + ReferenceCell::make_triangulation( + ReferenceCell::internal::make_reference_cell_from_int(kind), tria); + GridOut go; + go.write_msh(tria, "output.msh"); + cat_file("output.msh"); + // Uncomment to inspect each file separately + // if (out != "") + // go.write_msh(tria, out); + (void)out; +} + +int +main() +{ + initlog(); + + // Generate and print all reference cells + // Lines + test<1, 1>(1, "line_11.msh"); + test<1, 2>(1, "line_12.msh"); + test<1, 3>(1, "line_13.msh"); + + // Triangles + test<2, 2>(2, "tri_22.msh"); + test<2, 3>(2, "tri_23.msh"); + + // Quads + test<2, 2>(3, "quad_22.msh"); + test<2, 3>(3, "quad_23.msh"); + + // Tet + test<3, 3>(4, "tet_3.msh"); + + // Pyramid + test<3, 3>(5, "pyramid_3.msh"); + + // Wedge + test<3, 3>(6, "wedge_3.msh"); + + // Hex + test<3, 3>(7, "hex_3.msh"); +} diff --git a/tests/gmsh/gmsh_api_02.with_gmsh_with_api=on.output b/tests/gmsh/gmsh_api_02.with_gmsh_with_api=on.output new file mode 100644 index 0000000000..a2d77950a9 --- /dev/null +++ b/tests/gmsh/gmsh_api_02.with_gmsh_with_api=on.output @@ -0,0 +1,303 @@ + +DEAL::Testing kind(1) in dimensions <1,1> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +1 +0 1 "BoundaryID:1" +$EndPhysicalNames +$Entities +1 1 0 0 +2 0 0 0 1 1 +1 0 0 0 1 0 0 1 2 0 +$EndEntities +$Nodes +2 2 1 2 +0 2 0 0 +1 1 0 2 +1 +2 +0 0 0 +1 0 0 +$EndNodes +$Elements +1 1 1 1 +1 1 1 1 +1 1 2 +$EndElements + +DEAL::Testing kind(1) in dimensions <1,2> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +1 +0 1 "BoundaryID:1" +$EndPhysicalNames +$Entities +1 1 0 0 +2 0 0 0 1 1 +1 0 0 0 1 0 0 1 2 0 +$EndEntities +$Nodes +2 2 1 2 +0 2 0 0 +1 1 0 2 +1 +2 +0 0 0 +1 0 0 +$EndNodes +$Elements +1 1 1 1 +1 1 1 1 +1 1 2 +$EndElements + +DEAL::Testing kind(1) in dimensions <1,3> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +1 +0 1 "BoundaryID:1" +$EndPhysicalNames +$Entities +1 1 0 0 +2 0 0 0 1 1 +1 0 0 0 1 0 0 1 2 0 +$EndEntities +$Nodes +2 2 1 2 +0 2 0 0 +1 1 0 2 +1 +2 +0 0 0 +1 0 0 +$EndNodes +$Elements +1 1 1 1 +1 1 1 1 +1 1 2 +$EndElements + +DEAL::Testing kind(2) in dimensions <2,2> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 1 0 +1 0 0 0 1 1 0 1 1 0 +$EndEntities +$Nodes +1 3 1 3 +2 1 0 3 +1 +2 +3 +0 0 0 +1 0 0 +0 1 0 +$EndNodes +$Elements +1 1 1 1 +2 1 2 1 +1 1 2 3 +$EndElements + +DEAL::Testing kind(2) in dimensions <2,3> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 1 0 +1 0 0 0 1 1 0 1 1 0 +$EndEntities +$Nodes +1 3 1 3 +2 1 0 3 +1 +2 +3 +0 0 0 +1 0 0 +0 1 0 +$EndNodes +$Elements +1 1 1 1 +2 1 2 1 +1 1 2 3 +$EndElements + +DEAL::Testing kind(3) in dimensions <2,2> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 1 0 +1 0 0 0 1 1 0 1 1 0 +$EndEntities +$Nodes +1 4 1 4 +2 1 0 4 +1 +2 +3 +4 +0 0 0 +1 0 0 +0 1 0 +1 1 0 +$EndNodes +$Elements +1 1 1 1 +2 1 3 1 +1 1 2 4 3 +$EndElements + +DEAL::Testing kind(3) in dimensions <2,3> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 1 0 +1 0 0 0 1 1 0 1 1 0 +$EndEntities +$Nodes +1 4 1 4 +2 1 0 4 +1 +2 +3 +4 +0 0 0 +1 0 0 +0 1 0 +1 1 0 +$EndNodes +$Elements +1 1 1 1 +2 1 3 1 +1 1 2 4 3 +$EndElements + +DEAL::Testing kind(4) in dimensions <3,3> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 0 1 +1 0 0 0 1 1 1 1 1 0 +$EndEntities +$Nodes +1 4 1 4 +3 1 0 4 +1 +2 +3 +4 +0 0 0 +1 0 0 +0 1 0 +0 0 1 +$EndNodes +$Elements +1 1 1 1 +3 1 4 1 +1 1 2 3 4 +$EndElements + +DEAL::Testing kind(5) in dimensions <3,3> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 0 1 +1 -1 -1 0 1 1 1 1 1 0 +$EndEntities +$Nodes +1 5 1 5 +3 1 0 5 +1 +2 +3 +4 +5 +-1 -1 0 +1 -1 0 +-1 1 0 +1 1 0 +0 0 1 +$EndNodes +$Elements +1 1 1 1 +3 1 7 1 +1 1 2 4 3 5 +$EndElements + +DEAL::Testing kind(6) in dimensions <3,3> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 0 1 +1 0 0 0 1 1 1 1 1 0 +$EndEntities +$Nodes +1 6 1 6 +3 1 0 6 +1 +2 +3 +4 +5 +6 +1 0 0 +0 1 0 +0 0 0 +1 0 1 +0 1 1 +0 0 1 +$EndNodes +$Elements +1 1 1 1 +3 1 6 1 +1 1 2 3 4 5 6 +$EndElements + +DEAL::Testing kind(7) in dimensions <3,3> +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 0 1 +1 0 0 0 1 1 1 1 1 0 +$EndEntities +$Nodes +1 8 1 8 +3 1 0 8 +1 +2 +3 +4 +5 +6 +7 +8 +0 0 0 +1 0 0 +0 1 0 +1 1 0 +0 0 1 +1 0 1 +0 1 1 +1 1 1 +$EndNodes +$Elements +1 1 1 1 +3 1 5 1 +1 1 2 4 3 5 6 8 7 +$EndElements +