Skip to content

Commit

Permalink
Merge pull request #12967 from bangerth/comsol
Browse files Browse the repository at this point in the history 
More work for COMSOL meshes
  • Loading branch information
@drwells
drwells committed Nov 20, 2021
2 parents 079f8e9 + 9aee080 commit ebd1e52
Show file tree
Hide file tree
Showing 10 changed files with 4,020 additions and 73 deletions.
207 changes: 137 additions & 70 deletions source/grid/grid_in.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1615,61 +1615,93 @@ GridIn<dim, spacedim>::read_comsol_mphtxt(std::istream &in)
whole_file >> n_types;
for (unsigned int type = 0; type < n_types; ++type)
{
std::string object_type;
// The object type is prefixed by the number of characters the
// object type string takes up (e.g., 3 for 'tri' and 5 for
// 'prism'), but we really don't need that.
{
unsigned int dummy;
whole_file >> dummy;
}
whole_file >> object_type;

// Read the object type. Also do a number of safety checks.
std::string object_name;
whole_file >> object_name;

static const std::map<std::string, ReferenceCell> name_to_type = {
{"vtx", ReferenceCells::Vertex},
{"edg", ReferenceCells::Line},
{"tri", ReferenceCells::Triangle},
{"quad", ReferenceCells::Quadrilateral},
{"tet", ReferenceCells::Tetrahedron},
{"prism", ReferenceCells::Wedge}
// TODO: Add hexahedra and pyramids once we have a sample input file
// that contains these
};
AssertThrow(name_to_type.find(object_name) != name_to_type.end(),
ExcMessage("The input file contains a cell type <" +
object_name +
"> that the reader does not "
"current support"));
const ReferenceCell object_type = name_to_type.at(object_name);

unsigned int n_vertices_per_element;
whole_file >> n_vertices_per_element;

unsigned int n_elements;
whole_file >> n_elements;

if (object_type == "vtx")

if (object_type == ReferenceCells::Vertex)
{
AssertThrow(n_vertices_per_element == 1, ExcInternalError());
}
else if (object_type == "edg")
else if (object_type == ReferenceCells::Line)
{
AssertThrow(n_vertices_per_element == 2, ExcInternalError());
if (dim == 1)
cells.resize(n_elements);
else
subcelldata.boundary_lines.resize(n_elements);
}
else if (object_type == "tri")
else if (object_type == ReferenceCells::Triangle)
{
AssertThrow(dim >= 2,
ExcMessage("Triangles should not appear in input files "
"for 1d meshes."));
AssertThrow(n_vertices_per_element == 3, ExcInternalError());
if (dim == 2)
cells.resize(n_elements);
else
subcelldata.boundary_quads.resize(n_elements);
}
else if (object_type == "tet")
else if (object_type == ReferenceCells::Quadrilateral)
{
AssertThrow(dim >= 2,
ExcMessage(
"Quadrilaterals should not appear in input files "
"for 1d meshes."));
AssertThrow(n_vertices_per_element == 4, ExcInternalError());
}
else if (object_type == ReferenceCells::Tetrahedron)
{
AssertThrow(dim >= 3,
ExcMessage("Tetrahedra should not appear in input files "
"for 1d or 2d meshes."));
AssertThrow(n_vertices_per_element == 4, ExcInternalError());
if (dim == 3)
cells.resize(n_elements);
else
Assert(false, ExcInternalError());
}
else if (object_type == ReferenceCells::Wedge)
{
AssertThrow(dim >= 3,
ExcMessage(
"Prisms (wedges) should not appear in input files "
"for 1d or 2d meshes."));
AssertThrow(n_vertices_per_element == 6, ExcInternalError());
}
else
AssertThrow(false, ExcNotImplemented());


// Next, for each element read the vertex numbers. Then we have to decide
// what to do with it. If it is a vertex, we ignore the information.
// If it is a cell, we have to put it into the appropriate object, and the
// same if it is an edge or face.
// Next, for each element read the vertex numbers. Then we have
// to decide what to do with it. If it is a vertex, we ignore
// the information. If it is a cell, we have to put it into the
// appropriate object, and the same if it is an edge or
// face. Since multiple object type blocks can refer to cells or
// faces (e.g., for mixed meshes, or for prisms where there are
// boundary triangles and boundary quads), the element index 'e'
// below does not correspond to the index in the 'cells' or
// 'subcelldata.boundary_*' objects; we just keep pushing
// elements onto the back.
//
// In any case, we adjust vertex indices right after reading them based on
// the starting index read above
Expand All @@ -1684,75 +1716,110 @@ GridIn<dim, spacedim>::read_comsol_mphtxt(std::istream &in)
vertices_for_this_element[v] -= starting_vertex_index;
}

if (object_type == "vtx")
if (object_type == ReferenceCells::Vertex)
; // do nothing
else if (object_type == "edg")
else if (object_type == ReferenceCells::Line)
{
if (spacedim == 1)
cells[e].vertices = vertices_for_this_element;
{
cells.emplace_back();
cells.back().vertices = vertices_for_this_element;
}
else
subcelldata.boundary_lines[e].vertices =
vertices_for_this_element;
{
subcelldata.boundary_lines.emplace_back();
subcelldata.boundary_lines.back().vertices =
vertices_for_this_element;
}
}
else if (object_type == "tri")
else if ((object_type == ReferenceCells::Triangle) ||
(object_type == ReferenceCells::Quadrilateral))
{
if (spacedim == 2)
cells[e].vertices = vertices_for_this_element;
{
cells.emplace_back();
cells.back().vertices = vertices_for_this_element;
}
else
subcelldata.boundary_quads[e].vertices =
vertices_for_this_element;
{
subcelldata.boundary_quads.emplace_back();
subcelldata.boundary_quads.back().vertices =
vertices_for_this_element;
}
}
else if (object_type == "tet")
else if ((object_type == ReferenceCells::Tetrahedron) ||
(object_type == ReferenceCells::Wedge))
{
if (spacedim == 3)
cells[e].vertices = vertices_for_this_element;
{
cells.emplace_back();
cells.back().vertices = vertices_for_this_element;
}
else
Assert(false, ExcInternalError());
}
else
Assert(false, ExcNotImplemented());
}

// Then also read the "geometric entity indices". There need to be as
// many as there were elements to begin with
{
unsigned int dummy;
whole_file >> dummy;
AssertThrow(dummy == n_elements, ExcInternalError());
}

for (unsigned int e = 0; e < n_elements; ++e)
// Then also read the "geometric entity indices". There need to
// be as many as there were elements to begin with, or
// alternatively zero if no geometric entity indices will be set
// at all.
unsigned int n_geom_entity_indices;
whole_file >> n_geom_entity_indices;
AssertThrow((n_geom_entity_indices == 0) ||
(n_geom_entity_indices == n_elements),
ExcInternalError());

// Loop over these objects. Since we pushed them onto the back
// of various arrays before, we need to recalculate which index
// in these array element 'e' corresponds to when setting
// boundary and manifold indicators.
if (n_geom_entity_indices != 0)
{
AssertThrow(whole_file, ExcIO());
unsigned int geometric_entity_index;
whole_file >> geometric_entity_index;
if (object_type == "vtx")
; // do nothing
else if (object_type == "edg")
{
if (spacedim == 1)
cells[e].boundary_id = geometric_entity_index;
else
subcelldata.boundary_lines[e].boundary_id =
geometric_entity_index;
}
else if (object_type == "tri")
{
if (spacedim == 2)
cells[e].boundary_id = geometric_entity_index;
else
subcelldata.boundary_quads[e].boundary_id =
geometric_entity_index;
}
else if (object_type == "tet")
for (unsigned int e = 0; e < n_geom_entity_indices; ++e)
{
if (spacedim == 3)
cells[e].boundary_id = geometric_entity_index;
AssertThrow(whole_file, ExcIO());
unsigned int geometric_entity_index;
whole_file >> geometric_entity_index;
if (object_type == ReferenceCells::Vertex)
; // do nothing
else if (object_type == ReferenceCells::Line)
{
if (spacedim == 1)
cells[cells.size() - n_elements + e].material_id =
geometric_entity_index;
else
subcelldata
.boundary_lines[subcelldata.boundary_lines.size() -
n_elements + e]
.boundary_id = geometric_entity_index;
}
else if ((object_type == ReferenceCells::Triangle) ||
(object_type == ReferenceCells::Quadrilateral))
{
if (spacedim == 2)
cells[cells.size() - n_elements + e].material_id =
geometric_entity_index;
else
subcelldata
.boundary_quads[subcelldata.boundary_quads.size() -
n_elements + e]
.boundary_id = geometric_entity_index;
}
else if ((object_type == ReferenceCells::Tetrahedron) ||
(object_type == ReferenceCells::Wedge))
{
if (spacedim == 3)
cells[cells.size() - n_elements + e].material_id =
geometric_entity_index;
else
Assert(false, ExcInternalError());
}
else
Assert(false, ExcInternalError());
Assert(false, ExcNotImplemented());
}
else
Assert(false, ExcNotImplemented());
}
}
AssertThrow(whole_file, ExcIO());
Expand Down
6 changes: 3 additions & 3 deletions tests/grid/grid_in_comsol_mphtxt_01.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@



// check whether we can read in with the gmsh format
// Check whether we can read in a mesh in COMSOL's .mphtxt format

#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
Expand All @@ -29,7 +29,7 @@

template <int dim>
void
gmsh_grid(const char *name)
comsol_grid(const char *name)
{
Triangulation<dim> tria;
GridIn<dim> grid_in;
Expand Down Expand Up @@ -117,7 +117,7 @@ main()

try
{
gmsh_grid<3>(SOURCE_DIR "/grids/comsol/mesh_1.mphtxt");
comsol_grid<3>(SOURCE_DIR "/grids/comsol/mesh_1.mphtxt");
}
catch (std::exception &exc)
{
Expand Down

0 comments on commit ebd1e52

Please sign in to comment.