How to extract and use physical labels from GMSH surfaces on FEniCS (2D mesh)

mesh
1

Hi,

I would like to know how to extract physical surface label for a 2D mesh define on Gmsh.
For example, a square domain composed of 4 sub domains (with physical tag 1,2,3 and 4) : I would like to be able to integrate on each of these subdomain (dx(1), dx(2), dx(3) and dx(4)).

I naively tried (without success) to initialize my subdomain as I will do with the borders and the measure ds

mesh = Mesh()

with XDMFFile("/home/amoreno3/meshTESTplusieursDomaine.xdmf") as infile:
    infile.read(mesh)
File("Dolfin_circle_mesh.pvd").write(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)


with XDMFFile("/home/amoreno3/mfTESTplusieursDomaine.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("Dolfin_circle_facets.pvd").write(mf)


domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure("dx")(subdomain_data=domains)
dx1 = Measure("dx", domain=mesh, subdomain_data=mvc,   subdomain_id=1)
dx2 = Measure("dx", domain=mesh, subdomain_data=mvc, subdomain_id=2)
dx3 = Measure("dx", domain=mesh, subdomain_data=mvc, subdomain_id=3)
dx4 = Measure("dx", domain=mesh, subdomain_data=mvc, subdomain_id=4)

Thank you

1 Like
2

Dear @AMM,

the subdomain data should be the mesh function mf and not the mesh value collection mvc.

Could you please supply the output of

print(assemble(Constant(1)*dx)))
print(assemble(Constant(1)*dx1)))
print(assemble(Constant(1)*dx2)))
print(assemble(Constant(1)*dx3)))
print(assemble(Constant(1)*dx4)))

Additionally, you could the following instead of creating new measures for each data:

dx = Measure("dx",subdomain_data=mf)
print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(1)))
print(assemble(Constant(1)*dx(2)))
print(assemble(Constant(1)*dx(3)))
print(assemble(Constant(1)*dx(4)))
3

Hi Dokken,

Thank you for your answer.
I tried what you asked me to do but the two things don’t work.
For the first print, it tells me

> UFLException: This integral is missing an integration domain.

For the second print, using an integral (for example dx(4)) of the form assemble(Constant(1)*dx(4, domain=mesh, subdomain_data=mf)) to avoid the above error message, I get 0.0 as result.

One could think of error on the definition of my geometry on gmsh but I checked and it does not seem that it is the case. (the geo file is pretty simple as you can see below)

// Gmsh project created on Wed Jan 22 11:30:14 2020
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};
//+
Physical Surface("surface4", 4) = {3};
4

If you change dx to dx(domain=mesh) the first problem is solved.

5

@AMM
Secondly, writing your mesh to xdmf as described in the next snippets works for me:

import meshio
msh = meshio.read("amm.msh")

meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]}))
meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]},
                                    cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))

from dolfin import *
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)

domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure("dx",domain=mesh, subdomain_data=mf)
print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(1)))
print(assemble(Constant(1)*dx(2)))
print(assemble(Constant(1)*dx(3)))
print(assemble(Constant(1)*dx(4)))

as it gives the output:

1.0
0.25
0.25
0.25
0.25
1 Like
6

Ok, I just saw the problem : I write “line” instead of “triangle” for mf in the bash script that help me to create the xdmf file (and to prune them with meshio). it works well now.

By the way, running other code with this rectification (line to triangle for writing mf.xdmf) give me an error on the Dirichlet Boundary condition :
*** Error: Unable to create Dirichlet boundary condition.
*** Reason: User MeshFunction is not a facet MeshFunction (dimension is wrong).
*** Where: This error was encountered inside DirichletBC.cpp.

The code previously worked with “line” for mf.xdmf.

Can you explain to me why?
Do I have to write two mf files like mf1 (which gives the physical labels of surfaces with “triangle”) and mf2 (which gives the physical labels of borders with “line”)

Thank you in advance for your answer

7

You can write both the facet data and cell data to one xdmf file by doing the same thing for lines as for triangle. However, you Need to read it into a separate meshvaluecollection and meshfunction, Which takes in dimension 1 instead of 2.

8

Thank you very much. That’s now working well.

Maybe I can ask you one last question here (tell me if i should open another thread) : do you know how to define a discontinuous coefficient k on the whole domain (ie k = k_i\chi_i where \chi_i is the characteristic function of a subdomain and k_i are constant).

I would like something like for the 4 subdomains :

k_values = [1.5, 50, 2, 3] 
help = np.asarray(subdomains.array(), dtype=np.int32)
k.vector()[:] = np.choose(help, k_values)

But with that, I need to rename all my physical surface from gmsh. Do you have an idea on how to proceed?

Thank you for your help

9

You can do something like:

V = FunctionSpace(mesh, "DG", 0)
f = Function(V)
f.vector()[:] = 1.5*(1==subdomains.array()) + 50*(2==subdomains.array()) + ...
1 Like
10

Dokken,

I really appreciate your help and advice.

Thank you very much and have a good day.

11

Hi,

I’m sorry to bother you again with this question but it I cannot make it work again.

What did you take for “subdomain” in your lines of code below?

V = FunctionSpace(mesh, "DG", 0)
f = Function(V)
f.vector()[:] = 1.5*(1==subdomains.array()) + 50*(2==subdomains.array()) + ...

When I choose “subdomain = MeshFunction(“size_t”, mesh, mesh.topology().dim())”, it doesn’t work.

The thing is that I don’t know how and where my physical tag (for surfaces from gmsh) are stored.

Thank you in advance for your answer

AMM

12

Hi,
You need to load your data as described in the post above.

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

So subdomains would be equal to mf

13

Thank you very much (I wrote your mf file under two different files “mf” and “cf” and I tested only mf…)
That’s working.

Sorry for this trivial question.

14

I think this does not work anymore?

I am trying to repeat this example but stuck with this part of the code. I will write down what I have done so far.

Write xdmf file for surfaces and boundaries with the most recent working script that I can find:

import meshio
import numpy as np
msh = meshio.read("mesh.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[:,:2], 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)

This creates mesh.h5, mesh.xdmf and mf.xdmf files. Now tried to define mf and mvc and calculate the integrals:

from dolfin import * 
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)

But this gives an error:
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** fenics-support@googlegroups.com

*** 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.1.0

15

See the links in: Marking the boundaries of a 2D mesh for a mixed function-space problem - #4 by dokken

16

Actually that is also me :slight_smile: But this time I am trying to define multiple Physical Surfaces inside the domain as in this problem. Physical Lines do not cause a problem but dividing the domain and marking different regions is a problem.

17

My point is that:

  1. You are using an old and very complicated way of converting meshes (that is known to be error prone).
  2. The links in the previous posts shows you how to tag both surfaces and volumes.
  3. Without a reproducible example, there is not much we can do to help you
1 Like
18

I will be using FEniCSx someday for sure. Unfortunately I have to stick with Docker Toolbox and some older scripts for at least a couple of months.

Here is the .geo for a rectangle divided to two in gmsh and exported as Version 4 ASCI II:

Point(1) = {0, 0, 0, 0.1};
Point(2) = {1, 0, 0, 0.1};
Point(3) = {2, 0, 0, 0.1};
Point(4) = {2, 2, 0, 0.1};
Point(5) = {1, 2, 0, 0.1};
Point(6) = {0, 2, 0, 0.1};
//+
Line(7) = {1, 2};
Line(8) = {2, 5};
Line(9) = {5, 6};
Line(10) = {6, 1};
Line(11) = {2, 3};
Line(12) = {3, 4};
Line(13) = {4, 5};
//+
Curve Loop(14) = {7, 8, 9, 10};
Plane Surface(15) = {14};
Curve Loop(16) = {8, -13, -12, -11};
Plane Surface(17) = {16};
//+
Physical Surface(111) = {15};
Physical Surface(222) = {17};

The code above I used in this and the other topic this time gives:
AttributeError: ‘list’ object has no attribute ‘shape’

Method 1: I have achieved a good accuracy for my problem with mshr but it is of no use to me since I need to create my mesh in gmsh and mark my boundaries and areas:
(and I know this is an ancient method now)

from mshr import *
from fenics import *
import numpy as np

# This method is actually fine but I couldn't implement it to gmsh;
# I tried to find a way to implement this to gmsh but so far I couldn't...

domain_rectangle = Rectangle(Point(0, 0), Point(2, 2))
domain_L = Rectangle(Point(0, 0), Point(1, 2))

domain_setter = domain_rectangle

domain_rectangle.set_subdomain(1, domain_setter)
domain_rectangle.set_subdomain(2, domain_L)

mesh = generate_mesh(domain_rectangle, 10)

mf = MeshFunction('size_t',mesh,mesh.topology().dim(), mesh.domains())
dx_sub = Measure('dx', domain=mesh, subdomain_data=mf)

# Check if the area is marked correctly
File("Marked_Area.pvd").write(mf)

# Check if the area is calculated correctly
Method_1_A = assemble(1*dx_sub(1))
Method_1_S = assemble(1*dx_sub(2))
print("Total Area     = ", Method_1_A)
print("Subdomain Area = ", Method_1_S)

Method 2: If I create mesh with a single Physical Surface and use AutoSubdomain as I show below:

region_L = AutoSubDomain(lambda x: x[0] <= 1)
region_L.mark(mf_1, 1)

my integration gives unmeaningful results.

So the Method 1 is so far my best solution. But it is of no use to me.

I would really appreciate if you showed me a very simple example for this because what I am trying to achieve is to mark my boundaries and subdomains before I get my mesh from gmsh so that I can integrate for different areas within my 2D domain. Also, I actually thought this was one of the most recent methods :smiling_face_with_tear:

19

Consider Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken

1 Like
20

Thank you! I will try my best.