How to define a mesh in colab as the one imported?

Hi,

I am new on Fenics and i face a problem when I want to import meshes. I use Fenics on Google Colaboratory. In the following lines, i present what i have done :

step 1. I created the geometry (.geo) of my mesh using the script which follows…

Jt = 0.2 ;
Ht = 0.2 ;

D = 1. ;
B = 3*D ;
L = 5*D ;
S = D ;
H = Ht*D ;
jmax = Sqrt( ((D/2)^2)+H^2) ;
j = Jt*jmax ;

xI = -(((j/2)/jmax)*D)/2 ;
yI = -Sqrt(((j/2)^2)-(xI^2)) ;
xJ = -xI ;
yJ = -yI ;

xF = xI-Sqrt(((jmax-j)^2)-(H-2*yJ)^2) ;
yF = yI-(H-2*yJ) ;
xM = xJ + Sqrt(((jmax-j)^2)-(H-2*yJ)^2) ;
yM = yJ + (H-2*yJ) ;

xL = xM-D ;
yL = yM ;
xG = xF+D ;
yG = yF ;

xE = -(L/2) ;
yE = yF ;
xH = L/2 ;
yH = yF ;
xK = xE ;
yK = yL ;
xN = xH ;
yN = yL ;
 
xC = xE ;
yC = yE-S  ;
xD = xH  ;
yD = yH-S  ;
xO = xK  ;
yO = yK+S  ;
xP = xN  ;
yP = yN+S  ;

xA = xC  ;
yA = yC-(B-S) ;
xB = xD  ;
yB = yD-(B-S) ;
xQ = xO ;
yQ = yO+(B-S) ;
xR = xP ;
yR = yP+(B-S) ;

Point(1) = {xA, yA, 0, 1};
Point(2) = {xB, yB, 0, 1};
Point(3) = {xC, yC, 0, 1};
Point(4) = {xD, yD, 0, 1};
Point(5) = {xE, yE, 0, 1};
Point(6) = {xF, yF, 0, 1};
Point(7) = {xG, yG, 0, 1};
Point(8) = {xH, yH, 0, 1};
Point(9) = {xI, yI, 0, 1};
Point(10) = {xJ, yJ, 0, 1};
Point(11) = {xK, yK, 0, 1};
Point(12) = {xL, yL, 0, 1};
Point(13) = {xM, yM, 0, 1};
Point(14) = {xN, yN, 0, 1};
Point(15) = {xO, yO, 0, 1};
Point(16) = {xP, yP, 0, 1};
Point(17) = {xQ, yQ, 0, 1};
Point(18) = {xR, yR, 0, 1};


Line(1) = {1, 2};
Line(2) = {2, 4};
Line(3) = {8, 7};
Line(4) = {7, 10};
Line(5) = {10, 6};
Line(6) = {6, 5};
Line(7) = {3, 1};
Line(8) = {9, 12};
Line(9) = {12, 11};
Line(10) = {15, 17};
Line(11) = {17, 18};
Line(12) = {18, 16};
Line(13) = {14, 13};
Line(14) = {13, 9};
Line(15) = {15, 16};
Line(16) = {4, 3};
Line(17) = {11, 15};
Line(18) = {14, 16};
Line(19) = {5, 3};
Line(20) = {4, 8};
Line(21) = {12, 13};
Line(22) = {7, 6};
Line(55) = {10, 13};
Line(144) = {9, 6};


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


Plane Surface(29) = {23};
Plane Surface(30) = {24};
Plane Surface(31) = {25};
Plane Surface(32) = {26};
Plane Surface(33) = {27};
Plane Surface(34) = {28};


Transfinite Line {1, 11} = 5 ; 
Transfinite Line {-12, 10, 7, -2} = 2*5 Using Progression 1;
Transfinite Line {15, 16} = 500 ;
Transfinite Line {17, 18, 19, 20, 21, 22, 8, 14, 5, 4} = 100;
Transfinite Line {9, 13, 6, 3} = 200;

Physical Line(35) = {1} ;
Physical Line(36) = {2} ;
Physical Line(37) = {20} ;
Physical Line(38) = {3} ;
Physical Line(39) = {6} ;
Physical Line(40) = {19} ;
Physical Line(41) = {7} ;
Physical Line(42) = {10} ;
Physical Line(43) = {11} ;
Physical Line(44) = {12} ;
Physical Line(45) = {18} ;
Physical Line(46) = {13} ;
Physical Line(47) = {9} ;
Physical Line(48) = {17} ;
Physical Line(49) = {8} ;
Physical Line(50) = {4} ;
Physical Line(51) = {55} ;
Physical Line(52) = {144} ;

Physical Surface(53) = {29};
Physical Surface(54) = {30};
Physical Surface(56) = {31};
Physical Surface(57) = {32};
Physical Surface(58) = {33};
Physical Surface(59) = {34};

step 2. I created the mesh (.msh) from the previous geometry using the “2D mesh” button in Gmsh.

step 3. I would like to use the mesh (.msh) previously created in fenics (colab). I saw in some previous topics that I have to convert my (.msh) file to an (.xdmf) file so that it could be handled by Fenics. So i tried the following in order to install meshio and convert the file into xdmf :

!pip install -U meshio
!apt-get install python-lxml
!meshio-convert perfect_contact.msh perfect_contact.xdmf

step 4. Then I tried this in order to define the mesh for the FE simulation and plot the mesh in colab :

mesh = Mesh(“perfect_contact.xdmf”)
plt.figure(figsize=(10, 10))
plot(mesh)

Question : Please do you know a simple way for me to do this in google colab ? The error message is “Unable to open the file […] Unknown file type” I am a newbie, so please tell me on which step you give the correction… Or specify “everything” if i am completely wrong :frowning: Please, find enclosed a picture of the domain…

You need to use the XDMFFile to read the mesh:

mesh = Mesh()
with XDMFFile("perfect_contact.xdmf") as infile:
     infile.read(mesh)

which reads the mesh data into the mesh variable

Thank you. I tried this like you said :

# import and plot mesh 
!meshio-convert perfect_contact.msh perfect_contact.xdmf
mesh = Mesh()
with XDMFFile("perfect_contact.xdmf") as infile:
     infile.read(mesh)

but I faced this issue:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
     34 mesh = Mesh()
     35 with XDMFFile("perfect_contact.xdmf") as infile:
---> 36      infile.read(mesh)
     37 

RuntimeError: 

*** -------------------------------------------------------------------------
*** 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 recognise cell type.
*** Reason:  Unknown value "mixed".
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Please do you know how to fix this ?

You need to split the mesh data into one file for the cells, and one for the facets, see: Extracting mesh regions for gmsh mesh - #16 by dokken
You can load the facet data as described here:
How to extract and use physical labels from GMSH surfaces on FEniCS (2D mesh) - #12 by dokken

I wrote :

# import and plot mesh 

# part 1
mesh_from_file = !meshio-read("perfect_contact.msh")

def create_mesh(mesh, cell_type, prune_z=False):
    cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
    cell_data = np.hstack([mesh.cell_data_dict["gmsh:physical"][key]
                         for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type])
    # Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
    points= mesh.points
    if prune_z:
        points = points[:,:2]
    mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"perfect_contact":[cell_data]})
    return mesh_new
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
!meshio-write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
!meshio-write("mesh.xdmf", triangle_mesh)


# part 2

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

because I imported numpy as np and the lines meshio.read etc… didn’t worked on colab, I changed them to !meshio-read etc…, they seem to work now.

I face one last problem. I have this error message:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
     45     mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"perfect_contact":[cell_data]})
     46     return mesh_new
---> 47 line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
     48 get_ipython().system('meshio-write("facet_mesh.xdmf", line_mesh)')
     49 

<ipython-input-19-ba6cba785f07> in create_mesh(mesh, cell_type, prune_z)
     36 
     37 def create_mesh(mesh, cell_type, prune_z=False):
---> 38     cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
     39     cell_data = np.hstack([mesh.cell_data_dict["gmsh:physical"][key]
     40                          for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type])

AttributeError: 'SList' object has no attribute 'cells'

it is due to the line just after define_mesh :

cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])

Do you know a way to solve this last problem please ?

Why doesn’t import meshio and meshio.read work? You claim that your other command works, But it clearly doesn’t as it does not read in the mesh.
Try for instance:

print(mesh_from_file) 

What does this yield you? Does it give you a meshio.Mesh object?

with :

print(mesh_from_file)

I got :

[‘/bin/bash: -c: line 0: syntax error near unexpected token "perfect_contact.msh"\'', '/bin/bash: -c: line 0: meshio-read(“perfect_contact.msh”)'’]

I think my other command didn’t work as you said. I tried back what you sent me and I have this issue:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
     47     return mesh_new
     48 line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
---> 49 meshio.write("facet_mesh.xdmf", line_mesh)
     50 
     51 triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)

4 frames
/usr/local/lib/python3.6/dist-packages/meshio/_helpers.py in write(filename, mesh, file_format, **kwargs)
    145 
    146     # Write
--> 147     return writer(filename, mesh, **kwargs)

/usr/local/lib/python3.6/dist-packages/meshio/xdmf/main.py in write(*args, **kwargs)
    529 
    530 def write(*args, **kwargs):
--> 531     XdmfWriter(*args, **kwargs)
    532 
    533 

/usr/local/lib/python3.6/dist-packages/meshio/xdmf/main.py in __init__(self, filename, mesh, data_format, compression, compression_opts)
    342         if data_format == "HDF":
    343             self.h5_filename = os.path.splitext(self.filename)[0] + ".h5"
--> 344             self.h5_file = h5py.File(self.h5_filename, "w")
    345 
    346         xdmf_file = ET.Element("Xdmf", Version="3.0")

/usr/local/lib/python3.6/dist-packages/h5py/_hl/files.py in __init__(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, **kwds)
    403 
    404             with phil:
--> 405                 fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0, **kwds)
    406                 fid = make_fid(name, mode, userblock_size,
    407                                fapl, fcpl=make_fcpl(track_order=track_order),

/usr/local/lib/python3.6/dist-packages/h5py/_hl/files.py in make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0, **kwds)
    109         # we default to earliest
    110         low, high = h5f.LIBVER_EARLIEST, h5f.LIBVER_LATEST
--> 111     plist.set_libver_bounds(low, high)
    112 
    113     cache_settings = list(plist.get_cache())

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/h5p.pyx in h5py.h5p.PropFAID.set_libver_bounds()

ValueError: Invalid high library version bound (invalid high library version bound)

The solution to this error message can be found at: Gmsh 4.4.1 in FEniCS? Meshio - #9 by SantiagoOrtiz

I saw two solutions:

This one

FROM quay.io/fenicsproject/stable:current
USER root
RUN apt-get -qq update && \
    apt-get -y upgrade && \
    apt-get clean && \
    apt-get -y install python3-h5py && \
    pip install --upgrade pip && \
    pip install meshio[all] && \
    rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*
USER root

and this one

HDF5_DISABLE_VERSION_CHECK=1 python filename.py

But those commands didn’t directly worked on colabs, I had an error message. I tried this lines just to see the result:

!wget https://quay.io/fenicsproject/stable:current

!apt-get -qq update && \

!apt-get -y upgrade && \

!apt-get clean && \

!apt-get -y install python3-h5py && \

!pip install --upgrade pip && \

!pip install meshio[all] && \

!rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*

!pip install -U meshio

!apt-get install python-lxml

and I had those two messages after this line and what we have done…

2019.2.0.dev0
--2021-02-13 00:13:18--  https://quay.io/fenicsproject/stable:current
Resolving quay.io (quay.io)... 18.210.212.75, 52.72.182.67, 18.211.164.173, ...
Connecting to quay.io (quay.io)|18.210.212.75|:443... connected.
HTTP request sent, awaiting response... 308 PERMANENT REDIRECT
Location: https://quay.io/fenicsproject/stable:current/ [following]
--2021-02-13 00:13:18--  https://quay.io/fenicsproject/stable:current/
Reusing existing connection to quay.io:443.
HTTP request sent, awaiting response... 302 FOUND
Location: https://quay.io/repository/fenicsproject/stable?tag=current&tab=tags [following]
--2021-02-13 00:13:18--  https://quay.io/repository/fenicsproject/stable?tag=current&tab=tags
Reusing existing connection to quay.io:443.
HTTP request sent, awaiting response... 200 OK
Length: 74395 (73K) [text/html]
Saving to: ‘stable:current’

stable:current      100%[===================>]  72.65K  --.-KB/s    in 0.08s   

2021-02-13 00:13:18 (856 KB/s) - ‘stable:current’ saved [74395/74395]

/bin/bash: !apt-get: command not found

and

Requirement already up-to-date: meshio in /usr/local/lib/python3.6/dist-packages (4.3.8)
Requirement already satisfied, skipping upgrade: importlib-metadata; python_version < "3.8" in /usr/local/lib/python3.6/dist-packages (from meshio) (3.4.0)
Requirement already satisfied, skipping upgrade: numpy in /usr/local/lib/python3.6/dist-packages (from meshio) (1.19.5)
Requirement already satisfied, skipping upgrade: typing-extensions>=3.6.4; python_version < "3.8" in /usr/local/lib/python3.6/dist-packages (from importlib-metadata; python_version < "3.8"->meshio) (3.7.4.3)
Requirement already satisfied, skipping upgrade: zipp>=0.5 in /usr/local/lib/python3.6/dist-packages (from importlib-metadata; python_version < "3.8"->meshio) (3.4.0)
Reading package lists... Done
Building dependency tree       
Reading state information... Done
python-lxml is already the newest version (4.2.1-1ubuntu0.3).
0 upgraded, 0 newly installed, 0 to remove and 19 not upgraded.
Requirement already satisfied: h5py in /usr/local/lib/python3.6/dist-packages (2.10.0)
Requirement already satisfied: numpy>=1.7 in /usr/local/lib/python3.6/dist-packages (from h5py) (1.19.5)
Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from h5py) (1.15.0)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-d31982fdbecb> in <module>()
     39 import meshio
     40 
---> 41 mesh_from_file = meshio.read("perfect_contact.msh")
     42 
     43 def create_mesh(mesh, cell_type, prune_z=False):

4 frames
/usr/local/lib/python3.6/dist-packages/meshio/gmsh/_gmsh41.py in _read_nodes(f, is_ascii, data_size)
    187         # Store the point densely and in the order in which they appear in the file.
    188         # x(double) y(double) z(double) (* numNodes)
--> 189         points[ixx] = fromfile(f, c_double, num_nodes * 3).reshape((num_nodes, 3))
    190 
    191         # Entity tag and entity dimension of the nodes. Stored as point-data.

ValueError: cannot reshape array of size 61508 into shape (57395,3)

So first of all you need to get your installation right.
I have made a minimal recipe on how to use google co-lab with fenics and meshio (using a smaller mesh called mesh.msh.
See: https://colab.research.google.com/drive/1OBhNSF_EcmVaLiLnkcWJ4nNAi0kx6cXa?usp=sharing

Thank you.

I think I didn’t get this point :

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

do i have to modify something there ? because I have an error message, underlining this part :

infile.read(mvc, "name_to_read")

and this is the error message :

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-3-1d3b05d9b64d> in <module>()
     61 mvc = MeshValueCollection("size_t", mesh, 1)
     62 with XDMFFile("facet_mesh.xdmf") as infile:
---> 63     infile.read(mvc, "name_to_read")
     64 mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
     65 print(mesh.num_cells())

RuntimeError: 

*** -------------------------------------------------------------------------
*** 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.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

This is dependent on your original mesh, that you have converted to xdmf format. You should inspect the facet_mesh file in for instance paraview, to see if you spot anything. I would also suggest to look at the generated mesh, mesh.xdmf with paraview to make sure it looks Like what you expect.
See for instance: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #62 by dokken

Thank you, I would be able to fix my program with all your hints !

Hi, finally with your help i imported the mesh. This is the part that i added :

import meshio
mesh_from_file = meshio.read(“wear.msh”)
def create_mesh(mesh, cell_type, prune_z=False):
cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
cell_data = np.hstack([mesh.cell_data_dict[“gmsh:physical”][key]
for key in mesh.cell_data_dict[“gmsh:physical”].keys() if key==cell_type])
# Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
points= mesh.points
if prune_z:
points = points[:,:2]
mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})
return mesh_new
line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True)
meshio.write(“facet_mesh.xdmf”, line_mesh)

triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True)
meshio.write(“mesh.xdmf”, triangle_mesh)

then

mesh = Mesh()
with XDMFFile(“wear.xdmf”) as infile:
infile.read(mesh)
ndim = mesh.topology().dim()
plt.figure(figsize=(10, 10))
plot(mesh)

Please i have one last question. Can you say me how to correctly define my boundary conditions now ? i used to do

bottom = CompiledSubDomain(“near(x[1], %s, DOLFIN_EPS)”%(yF-B))
top = CompiledSubDomain(“near(x[1], %s, DOLFIN_EPS)”%(yL+B))

(the values (yF-B) and (yL+B) just refers to the region at the bottom and the top of the entire geometry. The geometry is presented at the beginning of the post).

I think i have to use “facet_mesh.xdmf” somewhere for the boundary conditions but i don’t know how… and this is my .geo file : (i haven’t define any region for the bc, is it correct like this ?)

Jt = 0.2 ;
Ht = 0.2 ;

D = 1. ;
B = 3D ;
L = 5
D ;
S = D ;
H = HtD ;
jmax = Sqrt( ((D/2)^2)+H^2) ;
j = Jt
jmax ;

xI = -(((j/2)/jmax)*D)/2 ;
yI = -Sqrt(((j/2)^2)-(xI^2)) ;
xJ = -xI ;
yJ = -yI ;

xF = xI-Sqrt(((jmax-j)^2)-(H-2yJ)^2) ;
yF = yI-(H-2
yJ) ;
xM = xJ + Sqrt(((jmax-j)^2)-(H-2yJ)^2) ;
yM = yJ + (H-2
yJ) ;

xL = xM-D ;
yL = yM ;
xG = xF+D ;
yG = yF ;

xE = -(L/2) ;
yE = yF ;
xH = L/2 ;
yH = yF ;
xK = xE ;
yK = yL ;
xN = xH ;
yN = yL ;

xC = xE ;
yC = yE-S ;
xD = xH ;
yD = yH-S ;
xO = xK ;
yO = yK+S ;
xP = xN ;
yP = yN+S ;

xA = xC ;
yA = yC-(B-S) ;
xB = xD ;
yB = yD-(B-S) ;
xQ = xO ;
yQ = yO+(B-S) ;
xR = xP ;
yR = yP+(B-S) ;

Point(1) = {xA, yA, 0, 1};
Point(2) = {xB, yB, 0, 1};
Point(3) = {xC, yC, 0, 1};
Point(4) = {xD, yD, 0, 1};
Point(5) = {xE, yE, 0, 1};
Point(6) = {xF, yF, 0, 1};
Point(7) = {xG, yG, 0, 1};
Point(8) = {xH, yH, 0, 1};
Point(9) = {xI, yI, 0, 1};
Point(10) = {xJ, yJ, 0, 1};
Point(11) = {xK, yK, 0, 1};
Point(12) = {xL, yL, 0, 1};
Point(13) = {xM, yM, 0, 1};
Point(14) = {xN, yN, 0, 1};
Point(15) = {xO, yO, 0, 1};
Point(16) = {xP, yP, 0, 1};
Point(17) = {xQ, yQ, 0, 1};
Point(18) = {xR, yR, 0, 1};

Line(1) = {1, 2};
Line(2) = {2, 4};
Line(3) = {8, 7};
Line(4) = {7, 10};
Line(5) = {10, 6};
Line(6) = {6, 5};
Line(7) = {3, 1};
Line(8) = {9, 12};
Line(9) = {12, 11};
Line(10) = {15, 17};
Line(11) = {17, 18};
Line(12) = {18, 16};
Line(13) = {14, 13};
Line(14) = {13, 9};
Line(15) = {15, 16};
Line(16) = {4, 3};
Line(17) = {11, 15};
Line(18) = {14, 16};
Line(19) = {5, 3};
Line(20) = {4, 8};
Line(21) = {12, 13};
Line(22) = {7, 6};
Line(55) = {10, 13};
Line(144) = {9, 6};

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

Plane Surface(29) = {23};
Plane Surface(30) = {24};
Plane Surface(31) = {25};
Plane Surface(32) = {26};
Plane Surface(33) = {27};
Plane Surface(34) = {28};

Transfinite Line {1, 11} = 5 ;
Transfinite Line {-12, 10, 7, -2} = 2*5 Using Progression 1;
Transfinite Line {15, 16} = 500 ;
Transfinite Line {17, 18, 19, 20, 21, 22, 8, 14, 5, 4} = 100;
Transfinite Line {9, 13, 6, 3} = 200;

Physical Line(35) = {1} ;
Physical Line(36) = {2} ;
Physical Line(37) = {20} ;
Physical Line(38) = {3} ;
Physical Line(39) = {6} ;
Physical Line(40) = {19} ;
Physical Line(41) = {7} ;
Physical Line(42) = {10} ;
Physical Line(43) = {11} ;
Physical Line(44) = {12} ;
Physical Line(45) = {18} ;
Physical Line(46) = {13} ;
Physical Line(47) = {9} ;
Physical Line(48) = {17} ;
Physical Line(49) = {8} ;
Physical Line(50) = {4} ;
Physical Line(51) = {55} ;
Physical Line(52) = {144} ;

Physical Surface(53) = {29};
Physical Surface(54) = {30};
Physical Surface(56) = {31};
Physical Surface(57) = {32};
Physical Surface(58) = {33};
Physical Surface(59) = {34};

You need to read the facet mesh into a mesh value collection, Which is in turn used in a MeshFunction, that can be used to restrict integral measures or apply boundary conditions.
See for instance the link i posted earlier:

I used MeshValueCollection like this :

import meshio
mesh_from_file = meshio.read(“wear.msh”)
def create_mesh(mesh, cell_type, prune_z=False):
cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
cell_data = np.hstack([mesh.cell_data_dict[“gmsh:physical”][key]
for key in mesh.cell_data_dict[“gmsh:physical”].keys() if key==cell_type])
# Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
points= mesh.points
if prune_z:
points = points[:,:2]
mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})
return mesh_new
line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True)
meshio.write(“facet_mesh.xdmf”, line_mesh)

triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True)
meshio.write(“mesh.xdmf”, triangle_mesh)

mesh = Mesh()
with XDMFFile(“mesh.xdmf”) as infile:
infile.read(mesh)

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

But i had this error :


RuntimeError Traceback (most recent call last)

in () 23 mvc = MeshValueCollection(“size_t”, mesh, 2) 24 with XDMFFile(“facet_mesh.xdmf”) as infile: —> 25 infile.read(mvc, “name_to_read”) 26 mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

RuntimeError: *** ------------------------------------------------------------------------- *** 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.2.0.dev0 *** Git changeset: unknown *** -------------------------------------------------------------------------

Even with your exemple on colab (Google Colab) and my .geo i got exactly the same error…

But my facet seems to be correct since i had this in paraview

Please, do you have an idea on what i have to do now ?

The dimension of the MVC should be 1, not 2, as you are reading in an entity of dimension 1 (line segmentd)

Thank you. I corrected that point, but i have the same error. The underlined line is

infile.read(mvc, “name_to_read”)

Please make a minimal failing example. This means reducing the complexity of your mesh file. Start by creating a single rectangle and tag each of the boundaries with an unique physical marker. If you are able to convert and read in that data, start expanding the mesh file until you hit an issue.

i am on it. Please could you just tell me the next step once i fix the mvc part ? In particular, how do i change the lines

bottom = CompiledSubDomain(“near(x[1], %s, DOLFIN_EPS)”%(yF-B))
top = CompiledSubDomain(“near(x[1], %s, DOLFIN_EPS)”%(yL+B))

in order to take into account the facet and mvc we defined earlier ?
The two bc are just the lines up and down in the picture of facet i presented earlier…