Problem in plotting the subdomains from file .xdmf

Loan_Bui · November 24, 2022, 7:46am

Dear,
I created a mesh of periodic cell from Gmsh with the following geometry:

// Gmsh project created on Mon Sep 14 09:33:57 2015
a = 1.;
R = 0.2;
t = Pi/2;
b = Sin(t)*a;
bR = Sin(t)*R;
c = a*Cos(t);
cR = R*Cos(t);

d = 1.;

nx = 80;
nR = R*nx;
ns = nx-2*nR;

Point(1) = {0, 0, 0, d};
Point(2) = {a, 0, 0, d};
Point(3) = {a+c, b, 0, d};
Point(4) = {c, b, 0, d};
Point(5) = {R,0,0,d};
Point(6) = {a-R,0,0,d};
Point(7) = {a+cR,bR,0,d};
Point(8) = {a+c-cR,b-bR,0,d};
Point(9) = {a+c-R,b,0,d};
Point(10) = {c+R,b,0,d};
Point(11) = {c-cR,b-bR,0,d};
Point(12) = {cR,bR,0,d};

Line(1) = {1, 5};
Line(2) = {5, 6};
Line(3) = {6, 2};
Line(4) = {2, 7};
Line(5) = {7, 8};
Line(6) = {8, 3};
Line(7) = {3, 9};
Line(8) = {9, 10};
Line(9) = {10, 4};
Line(10) = {4, 11};
Line(11) = {11, 12};
Line(12) = {12, 1};
Circle(13) = {5, 1, 12};
Circle(14) = {7, 2, 6};
Circle(15) = {9, 3, 8};
Circle(16) = {11, 4, 10};
Line Loop(17) = {13, -11, 16, -8, 15, -5, 14, -2};
Plane Surface(18) = {17};
Line Loop(19) = {1, 13, 12};
Plane Surface(20) = {19};
Line Loop(21) = {3, 4, 14};
Plane Surface(22) = {21};
Line Loop(23) = {7, 15, 6};
Plane Surface(24) = {23};
Line Loop(25) = {10, 16, 9};
Plane Surface(26) = {25};

Transfinite Line{1,3,4,6,7,9,10,12} = nR;
Transfinite Line{2,5,8,11} = ns;
Transfinite Line{14,16} = 2*nR;
Transfinite Line{13,15} = nR;

Periodic Line{1} = {9};
Periodic Line{2} = {8};
Periodic Line{3} = {7};
Periodic Line{4} = {12};
Periodic Line{5} = {11};
Periodic Line{6} = {10};

Physical Line(1) = {1,2,3};
Physical Line(2) = {4,5,6};
Physical Line(3) = {7,8,9};
Physical Line(4) = {10,11,12};

Physical Surface(0) = {18};
Physical Surface(1) = {20,22,24,26};

and save into the file: hexag_incl.msh

Then, this file was convert to the file: hexag_incl.xdmf

Then, I try to read and plot in fenics via the following code:

from __future__ import print_function
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

a = 1.         # unit cell width
b = 1. # unit cell height
c = 0        # horizontal offset of top boundary
R = 0.2        # inclusion radius
vol = a*b      # unit cell volume
# we define the unit cell vertices coordinates for later use
vertices = np.array([[0, 0.],
                     [a, 0.],
                     [a+c, b],
                     [c, b]])
mesh = Mesh()
with XDMFFile("hexag_incl1.xdmf") as in_file:
    in_file.read(mesh)

subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
plt.figure()
plot(subdomains)
plt.show()

but the output as follow

in this image, there is only one physical surface so it is not as true as it should be (in Gmsh, there are 2 physical surfaces)

Could you help me to find out the problem
best regards!

Loan_Bui · November 26, 2022, 3:31am

Please help me to solve the problem.
Actually, the created mesh in Gmsh is shown bellow:

Leo · November 26, 2022, 7:41am

Hi,

You need to change the last two lines of the gmsh file as following:

Physical Surface(1) = {18};
Physical Surface(2) = {20,22,24,26};

The above surfaces will be marked as surface#1 and surface#2 (They should not start from #0).
Here is what you get:

fig

Loan_Bui · November 26, 2022, 10:21am

I’ve already edit last two lines of gmsh as your proposal
however, nothing change in results

// Gmsh project created on Mon Sep 14 09:33:57 2015
a = 1.;
R = 0.2;
t = Pi/2;
b = Sin(t)*a;
bR = Sin(t)*R;
c = a*Cos(t);
cR = R*Cos(t);

d = 1.;

nx = 80;
nR = R*nx;
ns = nx-2*nR;

Point(1) = {0, 0, 0, d};
Point(2) = {a, 0, 0, d};
Point(3) = {a+c, b, 0, d};
Point(4) = {c, b, 0, d};
Point(5) = {R,0,0,d};
Point(6) = {a-R,0,0,d};
Point(7) = {a+cR,bR,0,d};
Point(8) = {a+c-cR,b-bR,0,d};
Point(9) = {a+c-R,b,0,d};
Point(10) = {c+R,b,0,d};
Point(11) = {c-cR,b-bR,0,d};
Point(12) = {cR,bR,0,d};

Line(1) = {1, 5};
Line(2) = {5, 6};
Line(3) = {6, 2};
Line(4) = {2, 7};
Line(5) = {7, 8};
Line(6) = {8, 3};
Line(7) = {3, 9};
Line(8) = {9, 10};
Line(9) = {10, 4};
Line(10) = {4, 11};
Line(11) = {11, 12};
Line(12) = {12, 1};
Circle(13) = {5, 1, 12};
Circle(14) = {7, 2, 6};
Circle(15) = {9, 3, 8};
Circle(16) = {11, 4, 10};
Line Loop(17) = {13, -11, 16, -8, 15, -5, 14, -2};
Plane Surface(18) = {17};
Line Loop(19) = {1, 13, 12};
Plane Surface(20) = {19};
Line Loop(21) = {3, 4, 14};
Plane Surface(22) = {21};
Line Loop(23) = {7, 15, 6};
Plane Surface(24) = {23};
Line Loop(25) = {10, 16, 9};
Plane Surface(26) = {25};

Transfinite Line{1,3,4,6,7,9,10,12} = nR;
Transfinite Line{2,5,8,11} = ns;
Transfinite Line{14,16} = 2*nR;
Transfinite Line{13,15} = nR;

Periodic Line{1} = {9};
Periodic Line{2} = {8};
Periodic Line{3} = {7};
Periodic Line{4} = {12};
Periodic Line{5} = {11};
Periodic Line{6} = {10};

Physical Line(1) = {1,2,3};
Physical Line(2) = {4,5,6};
Physical Line(3) = {7,8,9};
Physical Line(4) = {10,11,12};

Physical Surface(1) = {18};
Physical Surface(2) = {20,22,24,26};

the results as bellow:

Did you change any thing in the following fenics code?

```python
from __future__ import print_function
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

a = 1.         # unit cell width
b = 1. # unit cell height
c = 0        # horizontal offset of top boundary
R = 0.2        # inclusion radius
vol = a*b      # unit cell volume
# we define the unit cell vertices coordinates for later use
vertices = np.array([[0, 0.],
                     [a, 0.],
                     [a+c, b],
                     [c, b]])
mesh = Mesh()
with XDMFFile("hexag_incl2.xdmf") as in_file:
    in_file.read(mesh)

subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
plt.figure()
plot(subdomains)
plt.show()

Leo · November 26, 2022, 6:58pm

Alright. I tried two different ways to visualize the marked physical surfaces based on the gmsh file you provided. I go over each of them quickly. After saving your mesh created in gmsh as a file like: mesh.msh

1. Conversion into XDMF and visualization in Paraview:

Here is how you can convert your mesh.msh intor xdmf and visualize the physical surfaces in Paraview:

    from dolfin import *
    import numpy as np
    import meshio
    
    msh = meshio.read("mesh.msh")
    points = msh.points
    cells = np.vstack([cells.data for cells in msh.cells if cells.type == "triangle"])
    cell_data = np.vstack([msh.cell_data_dict["gmsh:physical"][key]for key in msh.cell_data_dict["gmsh:physical"].keys() if key =="triangle"])
    mesh = meshio.Mesh(points=points,cells=[("triangle", cells)],cell_data={"name_to_read": cell_data})
    meshio.write("mesh.xdmf", mesh)
    
    mesh = Mesh()
    with XDMFFile("mesh.xdmf") as infile:
        infile.read(mesh)
        mvcv = MeshValueCollection("size_t", mesh, mesh.topology().dim())
        infile.read(mvcv, "name_to_read")
    subdomains = cpp.mesh.MeshFunctionSizet(mesh, mvcv)
    
    File("physical surface.pvd") << subdomains

After opening physical surface.pvd in Paraview:

2. Conversion into the classic XML format and visualization by Matplotlib:

Here is another way that I tried and it also worked perfectly. You can convert your mesh.msh into the classic XML format by running:

  dolfin-convert mesh.msh MESH.xml

Then you can visualize the marked psychical surfaces as following:

  from dolfin import *
  import matplotlib.pyplot as plt
  
  mesh = Mesh("MESH.xml")
  boundaries = MeshFunction("size_t", mesh, "MESH_facet_region.xml")
  subdomains = MeshFunction("size_t", mesh, "MESH_physical_region.xml")
  
  plt.figure()
  plot(subdomains)
  plt.show()

Here is the output:

Loan_Bui · November 27, 2022, 1:46pm

Did you run this command in the terminal of docker desktop?
,
dolfin-convert mesh.msh MESH.xml
,

Loan_Bui · November 27, 2022, 2:24pm

when I call “dolfin-convert hexag_incl.msh hexag_incl.xml” in the terminal docker desktop, an error happened as following:

dokken · November 27, 2022, 4:16pm

As @Leo says, he assumed that the mesh file generated by Gmsh is called mesh.msh

It Seems like you do not have the msh file
in the directory you are running from

Loan_Bui · November 28, 2022, 1:28am

Yes, I understand, for me, the mesh file generated by Gmsh is save as hexag_incl.msh
This file was saved in the fenics directory as following:

then, I call the command bellow in the terminal docker desktop:
,
dolfin-convert hexag_incl.msh Hexag_incl.xml
,
but it is mentioned that : “No such file or directory: 'hexag_incl.msh”
I dont understand why
Could you help me, please!

dokken · November 28, 2022, 5:40am

Did you share the directory (mount the volume in docker lingo) with the docker container? What is the output of ls in the docker terminal?

Loan_Bui · November 28, 2022, 7:27am

I don’t know, could you tell me how to check this information, please?

dokken · November 28, 2022, 7:46am

Call the command ls in the terminal where you previously called dolfin-convert.
Regarding sharing repositories with containers, see

Loan_Bui · November 28, 2022, 8:33am

Here the output when I call ls in the terminal:

,

ls

Hexag_incl.xml TrackFenics_old WELCOME demo fenics.env.conf local shared tutorial-env untitled.txt

dokken · November 28, 2022, 8:51am

As you can see, the hexag_incl.msh file is not in this directory.

Loan_Bui · November 28, 2022, 9:12am

I’ve loading the hexag_incl.msh in this directory
as you can see, now when I call ls to the terminal, it appears the notices:
,

ls

Hexag_incl.xml TrackFenics_old WELCOME demo fenics.env.conf hexag_incl.msh local shared tutorial-env untitled.txt

However, when I recall the command: “dolfin-convert hexag_incl.msh Hexag_incl.xml”
another errors appears:

conpierce8 · November 28, 2022, 5:22pm

It might be helpful to back up for a minute. The original issue in this thread had nothing to do with mesh conversion; rather, the issue was that the subdomains were not being displayed correctly by the plot command. It looks like you were able to successfully convert the .msh file to .xdmf.

The subdomains were not being displayed correctly because you had not loaded them from the XDMF file. The following code, copied from your first post, loads only the mesh:

mesh = Mesh()
with XDMFFile("hexag_incl1.xdmf") as in_file:
    in_file.read(mesh)

If you change this to the following code, you should see the subdomains plotted correctly. The MeshValueCollection is the key piece that is needed to load the domain markers from the file.

mesh = Mesh()
with XDMFFile("hexag_incl1.xdmf") as in_file:
    in_file.read(mesh)
    mvc_cells = MeshValueCollection("size_t", mesh, mesh.topology().dim())
    in_file.read(mvc_cells, "name_to_read")
subdomains = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)
plot(subdomains)
plt.show()

This avoids the use of dolfin-convert, which is no longer recommended.

Loan_Bui · November 29, 2022, 3:08am

Thank you
the problem is solved by converting to .xdmf file

Topic		Replies	Views
Find coordinates of cells in a subdomain General	3	204	January 11, 2024
Error in extracting surface and boundary markers for a compiled subdomain from its parent mesh	2	419	September 22, 2021
Initialize in subdomain from gmsh	5	1191	May 24, 2021
How I can differ between two subdomains for gmsh file? mesh	2	777	May 7, 2020
Creating subdomains and boundary conditions on imported gmsh mesh General	13	300	June 17, 2024

Problem in plotting the subdomains from file .xdmf

ls

ls

Related topics