Define Subdomains created from .geo into XDMF file into dolfin legacy

General
1

Hey everyone. I’ve been trying to import a .geo file generated mesh into legacy FeniCs, in order to solve a problem that concerns subdomains. However, I haven´t had any luck with defining the subdomains. I took an answer from a while ago and attempted to modify it, put I havent been able to do so.

The code looks like this:

  import meshio
  import numpy as np
  from dolfin import *
  
  msh = meshio.read("idea3.msh")
  
  line_cells = []
  for cell in msh.cells:
      if cell.type == "triangle":
          triangle_cells = cell.data
      elif  cell.type == "line":
          if len(line_cells) == 0:
              line_cells = cell.data
          else:
              line_cells = np.vstack([line_cells, cell.data])
  
  line_data = []
  for key in msh.cell_data_dict["gmsh:physical"].keys():
      if key == "line":
          if len(line_data) == 0:
              line_data = msh.cell_data_dict["gmsh:physical"][key]
          else:
              line_data = np.vstack([line_data, msh.cell_data_dict["gmsh:physical"][key]])
      elif key == "triangle":
          triangle_data = msh.cell_data_dict["gmsh:physical"][key]
  
  triangle_mesh = meshio.Mesh(points=msh.points, cells={"triangle": triangle_cells})
  line_mesh =meshio.Mesh(points=msh.points,
                             cells=[("line", line_cells)],
                             cell_data={"name_to_read":[line_data]})
  meshio.write("mesh.xdmf", triangle_mesh)
  meshio.xdmf.write("mf.xdmf", line_mesh)
  
  mesh = Mesh()
  with XDMFFile("mesh.xdmf") as infile:
      infile.read(mesh)
  mvc = MeshValueCollection("size_t", mesh, 2)
  with XDMFFile("mf.xdmf") as infile:
      infile.read(mvc, "name_to_read")
  mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

and the error it reads is as follows:

  *** Remember to include the error message listed below and, if possible,
  *** include a *minimal* running example to reproduce the error.
  ***
  *** -------------------------------------------------------------------------
  *** Error:   Unable to find entity in map.
  *** Reason:  Error reading MeshValueCollection.
  *** Where:   This error was encountered inside HDF5File.cpp.
  *** Process: 0
  *** 
  *** DOLFIN version: 2019.2.0.dev0
  *** Git changeset:  fd662efc1c0ddefa341c1dac753f717f6b6292a8

The .geo file is a simple file:

  SetFactory("OpenCASCADE");
  //+
  Point(1) = {0, 0, 0, 1.0};
  //+
  Point(2) = {1, 0, 0, 1.0};
  Point(3) = {1, 1, 0, 1.0};
  Point(4) = {0, 1, 0, 1.0};
  //+
  Point(5) = {0.5, 1, 0, 1.0};
  //+
  Point(6) = {0, 0.5, 0, 1.0};
  //+
  Point(7) = {1, 0.5, 0, 1.0};
  //+
  Point(8) = {0.5, -0, 0, 1.0};
  
  //+
  Point(9) = {0.5, 0.5, 0, 1.0};
  //+
  
  //+
  Line(1) = {4, 6};
  //+
  Line(2) = {6, 1};
  //+
  Line(3) = {1, 8};
  //+
  Line(4) = {8, 2};
  //+
  Line(5) = {2, 7};
  //+
  Line(6) = {7, 3};
  //+
  Line(7) = {3, 5};
  //+
  Line(8) = {5, 4};
  //+
  Line(9) = {5, 9};
  //+
  Line(10) = {9, 6};
  //+
  Line(11) = {9, 8};
  //+
  Line(12) = {9, 7};
  //+
  Curve Loop(1) = {10, 2, 3, -11};
  //+
  Plane Surface(1) = {1};
  //+
  Curve Loop(2) = {12, -5, -4, -11};
  //+
  Plane Surface(2) = {2};
  //+
  Curve Loop(3) = {7, 9, 12, 6};
  //+
  Plane Surface(3) = {3};
  //+
  Curve Loop(4) = {8, 1, -10, -9};
  //+
  Plane Surface(4) = {4};
  //+
  Physical Curve("inlet", 10000) = {2, 1};
  //+
  Physical Curve("outlet", 10001) = {6, 5};
  //+
  Physical Curve("wall", 10002) = {7, 8, 3, 4};
  //+
  Physical Surface("surface1", 1) = {1};
  //+
  Physical Surface("surface2", 2) = {2};
  //+
  Physical Surface("surface3", 3) = {4};
  //+

And Im working with the following versions:

dolfin: 2019.2.0.dev0
meshio: 5.3.4

Any help into how to read my .geo as a mesh, and access to its defined mesh will be highly appreciated

2

Your mesh generation is inconsistent.
If we consider the following code:

import numpy as np
from IPython import embed
import dolfin
import numpy
import meshio


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           "name_to_read": [cell_data]})
    return out_mesh


mesh_from_file = meshio.read("test_mesh.msh")

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("mf.xdmf", line_mesh)
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

assert np.allclose(line_mesh.points, triangle_mesh.points)
vertices_in_triangles = np.unique(
    triangle_mesh.cells_dict["triangle"].reshape(-1))
print(vertices_in_triangles)

vertices_in_lines = np.unique(line_mesh.cells_dict["line"].reshape(-1))
print(vertices_in_lines)

you see that you get
[ 0 1 3 4 5 6 7 8 9 10 11]
as the indices of the vertices in any triangle, while you have
[0 1 2 3 4 5 6 7]
in your facets (line segments).
As vertex 2 is not part of the cells, you cannot find the facet with [2,4] in the mesh:

In [1]: line_mesh.cells_dict["line"]
Out[1]: 
array([[3, 5],
       [5, 0],
       [0, 7],
       [7, 1],
       [1, 6],
       [6, 2],
       [2, 4],
       [4, 3]])

In general the number of points in the mesh should be equal to the number of unique nodes in the mesh topology, i.e.
len(np.unique(triangle_mesh.cells_dict["triangle"].reshape(-1)))
should be equal to
triangle_mesh.points.shape[0]

This is probably because you haven’t added all the physical surfaces (as you add 4 plane surfaces above, but only 3 physical surfaces

3

Dear Dokken, my apologies, I changed it to this .geo and I’m still getting the same error in FeniCs

  SetFactory("OpenCASCADE");
  //+
  Point(1) = {0, 0, 0, 1.0};
  //+
  Point(2) = {1, 0, 0, 1.0};
  Point(3) = {1, 1, 0, 1.0};
  Point(4) = {0, 1, 0, 1.0};
  //+
  Point(5) = {0.5, 1, 0, 1.0};
  //+
  Point(6) = {0, 0.5, 0, 1.0};
  //+
  Point(7) = {1, 0.5, 0, 1.0};
  //+
  Point(8) = {0.5, -0, 0, 1.0};
  
  //+
  Point(9) = {0.5, 0.5, 0, 1.0};
  //+
  
  //+
  Line(1) = {4, 6};
  //+
  Line(2) = {6, 1};
  //+
  Line(3) = {1, 8};
  //+
  Line(4) = {8, 2};
  //+
  Line(5) = {2, 7};
  //+
  Line(6) = {7, 3};
  //+
  Line(7) = {3, 5};
  //+
  Line(8) = {5, 4};
  //+
  Line(9) = {5, 9};
  //+
  Line(10) = {9, 6};
  //+
  Line(11) = {9, 8};
  //+
  Line(12) = {9, 7};
  //+
  Curve Loop(1) = {10, 2, 3, -11};
  //+
  Plane Surface(1) = {1};
  //+
  Curve Loop(2) = {12, -5, -4, -11};
  //+
  Plane Surface(2) = {2};
  //+
  Curve Loop(3) = {7, 9, 12, 6};
  //+
  Plane Surface(3) = {3};
  //+
  Curve Loop(4) = {8, 1, -10, -9};
  //+
  Plane Surface(4) = {4};
  //+
  Physical Curve("inlet", 10000) = {2, 1};
  //+
  Physical Curve("outlet", 10001) = {6, 5};
  //+
  Physical Curve("wall", 10002) = {7, 8, 3, 4};
  //+
  Physical Surface("surface1", 1) = {1};
  //+
  Physical Surface("surface2", 2) = {2};
  //+
  Physical Surface("surface3", 3) = {3};
  //+
  Physical Surface("surface4", 4) = {4};
4

This not current, as the last argument should be 1.
Your new mesh runs fine on my system with:

import numpy as np
from IPython import embed
import dolfin
import numpy
import meshio


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           "name_to_read": [cell_data]})
    return out_mesh


mesh_from_file = meshio.read("test_mesh.msh")

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("mf.xdmf", line_mesh)
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

assert np.allclose(line_mesh.points, triangle_mesh.points)
vertices_in_triangles = np.unique(
    triangle_mesh.cells_dict["triangle"].reshape(-1))

mesh = dolfin.Mesh()
with dolfin.XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = dolfin.MeshValueCollection("size_t", mesh, 1)
with dolfin.XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
5

Worked correctly. Thanks, Dokken.