Regarding equivalent representation of dolfin.Mesh module expressions in dolfinx.gmshio

19ec94 · December 19, 2023, 4:21pm

Hi,
While porting my dolfin code to dolfinx, I have stumbled upon the following issue. I have a 3D geometry file (input3D.geo, attached below) which I used in Gmsh to generate 3D mesh. This file contains tags for all the physical subdomains (cell_tags) and tags for all the physical boundaries (facet_tags).

To read XDMF mesh file in dolfin, I followed the (outdated) three-step approach described below:

Use gmsh to convert input3D.geo to input3D.msh.
Use dolfin-convert to convert input3D.msh to input3D.xml , input3D_physical_region.xml , input3D_facet_region.xml.

import os
os.system("dolfin-convert {0} {1}".format("input3D.msh", "input3D.xml"))
import dolfin as df
mesh = df.Mesh("input3D.xml")
subdomains = df.MeshFunction("size_t",mesh,"input3D_physical_region.xml")
boundaries = df.MeshFunction("size_t",mesh,"input3D_facet_region.xml")

Use dolfin.HDF5File to convert ________.xml to input3D.h5 .

file = df.HDF5File(mesh.mpi_comm(), "input3D.h5", "w")
file.write(mesh, "/mesh")
file.write(subdomains, "/subdomains")
file.write(boundaries, "/boundaries")
file.close()

Following the above steps, I was successfully able to read the mesh file back into dolfin by

import dolfin as df
mesh = df.Mesh()
hdf = df.HDF5File(self.mesh.mpi_comm(), "input3D.h5", "r")
hdf.read(mesh, "/mesh", False)
dim = mesh.topology().dim()
subdomains = df.MeshFunction("size_t", mesh, dim)
hdf.read(subdomains, "/subdomains")
boundaries = df.MeshFunction("size_t", mesh, dim - 1)
hdf.read(boundaries, "/boundaries")

I went through all the effort to define custom ufl measures, which proved to be very helpful when imposing boundary conditions weakly.

import ufl
ufl.dx = ufl.Measure("dx", domain=mesh, subdomain_data=subdomains)
ufl.ds = ufl.Measure("ds", domain=mesh, subdomain_data=boundaries)
ufl.dS = ufl.Measure("dS", domain=mesh, subdomain_data=boundaries)

I investigated the data by printing the entries in subdomains and boundaries variables. As expected, there is an associated tag-id for each cell. Hence the lengths of subdomain array and boundaries array are equal to the total number of cells in the mesh. The boundaries array contains either zero (for inner cells) or a tag-id (for cells on the boundary) as defined in the gmsh.

print(self.mesh.num_cells())
print(self.subdomains.array().size)
print(self.subdomains.array().size)

I want to achieve the same end result in dolfinx. To this end, I converted the input3D.geo file to input3D.msh using Gmsh software. I then used the following snippet to obtain the numerical tags.

import dolfinx as dfx
from mpi4py import MPI 

mesh, cell_tags, facet_tags = dfx.io.gmshio.read_from_msh(\
        "input3D.msh", MPI.COMM_WORLD, gdim=3)j

Upon investigation the length of facet_tags is not equal to the number of total cells in the mesh. It only returns the non-zero id, which on the first look makes sense.

print(len(cell_tags.indices))
print(len(facet_tags.indices))
print(mesh_dfx.topology.index_map(mesh_dfx.topology.dim).size_local)

In short, I want to use facet_tags just as I have used boundaries in old dolfin. Could anyone confirm if the boundaries in dolfin and facet_tags in dolfinx mentioned above in my explanation are equivalent?
Ultimately I’d like to define some custom ufl measures using these tags like below.

import ufl
ufl.dx = ufl.Measure("dx", domain=mesh_dfx, subdomain_data=cell_tags)
ufl.ds = ufl.Measure("ds", domain=mesh_dfx, subdomain_data=facet_tags)
ufl.dS = ufl.Measure("dS", domain=mesh_dfx, subdomain_data=facet_tags)

Here is the (long and complicated) input3D.geo file for testing in case you’re interested.

//+
SetFactory("OpenCASCADE");
//Mesh.CharacteristicLengthFromPoints = 1;
//Mesh.CharacteristicLengthFromCurvature = 1;
Mesh.CharacteristicLengthMax = 0.5;
lc = 0.1;
//Mesh.RandomFactor = 5.e-10;
// Choose 2D Algorithm
// 1 : MeshAdapt
// 2 : Automatic
// 5 : Delaunay
// 6 : Frontal-Delaunay
Mesh.Algorithm = 5;
// Choose 3D Algorithm
// 1: Delaunay 
Mesh.Algorithm3D = 1; 
//
vane_diagonal = 1.414235;
radius = 2 * vane_diagonal;
xmid = 0.75 * vane_diagonal;
ymid = 0;
zmid = 0.75 * vane_diagonal;
vane_start = xmid - (vane_diagonal/2);

vane_thickness = 0.005; 
lx = vane_start + vane_diagonal; 
ly = vane_diagonal; 
lz = vane_start + vane_diagonal;

// Right Vane
// (Hot) POINTS

p01 = newp; Point(p01) = {+vane_start, +ymid, +vane_thickness, lc};
p02 = newp; Point(p02) = {+xmid, +ly/2, +vane_thickness, lc};
p03 = newp; Point(p03) = {+lx, +ymid, +vane_thickness, lc};
p04 = newp; Point(p04) = {+xmid, -ly/2, +vane_thickness, lc};
// (Cold) Points
p05 = newp; Point(p05) = {+vane_start, +ymid, -vane_thickness, lc};
p06 = newp; Point(p06) = {+xmid, +ly/2, -vane_thickness, lc};
p07 = newp; Point(p07) = {+lx, +ymid, -vane_thickness, lc};
p08 = newp; Point(p08) = {+xmid, -ly/2, -vane_thickness, lc};
// Bottom Vane
// (Hot) POINTS
p09 = newp; Point(p09) = { +vane_thickness,+ymid, -vane_start, lc};
p10 = newp; Point(p10) = { +vane_thickness,+ly/2, -zmid, lc};
p11 = newp; Point(p11) = { +vane_thickness,+ymid, -lz, lc};
p12 = newp; Point(p12) = { +vane_thickness,-ly/2,-zmid, lc};
// (Cold) Points
p13 = newp; Point(p13) = { -vane_thickness, +ymid, -vane_start, lc};
p14 = newp; Point(p14) = { -vane_thickness, +ly/2, -zmid,lc};
p15 = newp; Point(p15) = { -vane_thickness, +ymid, -lz, lc};
p16 = newp; Point(p16) = { -vane_thickness, -ly/2, -zmid, lc};
// Left Vane
//(HOT) Points
p17 = newp; Point(p17) = {-vane_start, +ymid, -vane_thickness, lc};
p18 = newp; Point(p18) = {-xmid, +ly/2, -vane_thickness, lc};
p19 = newp; Point(p19) = {-lx, +ymid, -vane_thickness, lc};
p20 = newp; Point(p20) = {-xmid, -ly/2, -vane_thickness, lc};
//(COLD) POINTS
p21 = newp; Point(p21) = {-vane_start, +ymid, +vane_thickness, lc};
p22 = newp; Point(p22) = {-xmid, +ly/2, +vane_thickness, lc};
p23 = newp; Point(p23) = {-lx, +ymid, +vane_thickness, lc};
p24 = newp; Point(p24) = {-xmid, -ly/2, +vane_thickness, lc};
// Top Vane
// (HOT) Points
p25 = newp; Point(p25) = { -vane_thickness, +ymid, vane_start, lc};
p26 = newp; Point(p26) = { -vane_thickness, +ly/2, zmid,lc};
p27 = newp; Point(p27) = { -vane_thickness, +ymid, lz, lc};
p28 = newp; Point(p28) = { -vane_thickness, -ly/2, zmid, lc};
// (COLD) POINTS
p29 = newp; Point(p29) = { +vane_thickness,+ymid, vane_start, lc};
p30 = newp; Point(p30) = { +vane_thickness,+ly/2, zmid, lc};
p31 = newp; Point(p31) = { +vane_thickness,+ymid, lz, lc};
p32 = newp; Point(p32) = { +vane_thickness,-ly/2, zmid, lc};

//Connecting Lines
// Right vane
// HOT
l01 = newl; Line(l01) = {p01, p02};
l02 = newl; Line(l02) = {p02, p03};
l03 = newl; Line(l03) = {p03, p04};
l04 = newl; Line(l04) = {p04, p01};
// COLD
l05 = newl; Line(l05) = {p05, p06};
l06 = newl; Line(l06) = {p06, p07};
l07 = newl; Line(l07) = {p07, p08};
l08 = newl; Line(l08) = {p08, p05};
// DOWN VANE
// HOT
l09 = newl; Line(l09) = {p09, p10};
l10 = newl; Line(l10) = {p10, p11};
l11 = newl; Line(l11) = {p11, p12};
l12 = newl; Line(l12) = {p12, p09};
// COLD
l13 = newl; Line(l13) = {p13, p14};
l14 = newl; Line(l14) = {p14, p15};
l15 = newl; Line(l15) = {p15, p16};
l16 = newl; Line(l16) = {p16, p13};
// LEFT vane
// HOT
l17 = newl; Line(l17) = {p17, p18};
l18 = newl; Line(l18) = {p18, p19};
l19 = newl; Line(l19) = {p19, p20};
l20 = newl; Line(l20) = {p20, p17};
// COLD
l21 = newl; Line(l21) = {p21, p22};
l22 = newl; Line(l22) = {p22, p23};
l23 = newl; Line(l23) = {p23, p24};
l24 = newl; Line(l24) = {p24, p21};
// UP
// HOT
l25 = newl; Line(l25) = {p25, p26};
l26 = newl; Line(l26) = {p26, p27};
l27 = newl; Line(l27) = {p27, p28};
l28 = newl; Line(l28) = {p28, p25};
// COLD
l29 = newl; Line(l29) = {p29, p30};
l30 = newl; Line(l30) = {p30, p31};
l31 = newl; Line(l31) = {p31, p32};
l32 = newl; Line(l32) = {p32, p29};
// RIGHT HOT -> COLD
l33 = newl; Line(l33) = {p01, p05};
l34 = newl; Line(l34) = {p02, p06};
l35 = newl; Line(l35) = {p03, p07};
l36 = newl; Line(l36) = {p04, p08};
// DOWN HOT -> COLD
l37 = newl; Line(l37) = {p09, p13};
l38 = newl; Line(l38) = {p10, p14};
l39 = newl; Line(l39) = {p11, p15};
l40 = newl; Line(l40) = {p12, p16};
// LEFT HOT -> COLD
l41 = newl; Line(l41) = {p17, p21};
l42 = newl; Line(l42) = {p18, p22};
l43 = newl; Line(l43) = {p19, p23};
l44 = newl; Line(l44) = {p20, p24};
// UP HOT -> COLD
l45 = newl; Line(l45) = {p25, p29};
l46 = newl; Line(l46) = {p26, p30};
l47 = newl; Line(l47) = {p27, p31};
l48 = newl; Line(l48) = {p28, p32};

//CURVES
// RIGHT
rh_cl = newcl; Curve Loop(rh_cl) = {l01, l02, l03, l04};
rc_cl = newcl; Curve Loop(rc_cl) = {l05, l06, l07, l08};
rt1_cl = newcl; Curve Loop(rt1_cl) = {+l01,+l34,-l05,-l33};
rt2_cl = newcl; Curve Loop(rt2_cl) = {+l02,+l35,-l06,-l34};
rb1_cl = newcl; Curve Loop(rb1_cl) = {-l04,+l36,+l08,-l33};
rb2_cl = newcl; Curve Loop(rb2_cl) = {-l03,+l35,+l07,-l36};
// DOWN
dh_cl = newcl; Curve Loop(dh_cl) = {l09, l10, l11, l12};
dc_cl = newcl; Curve Loop(dc_cl) = {l13, l14, l15, l16};
dt1_cl = newcl; Curve Loop(dt1_cl) = {+l09,+l38,-l13,-l37};
dt2_cl = newcl; Curve Loop(dt2_cl) = {+l10,+l39,-l14,-l38};
db1_cl = newcl; Curve Loop(db1_cl) = {-l12,+l40,+l16,-l37};
db2_cl = newcl; Curve Loop(db2_cl) = {-l11,+l39,+l15,-l40};
// LEFT
lh_cl = newcl; Curve Loop(lh_cl) = {l17, l18, l19, l20};
lc_cl = newcl; Curve Loop(lc_cl) = {l21, l22, l23, l24};
lt1_cl = newcl; Curve Loop(lt1_cl) = {+l17,+l42,-l21,-l41};
lt2_cl = newcl; Curve Loop(lt2_cl) = {+l18,+l43,-l22,-l42};
lb1_cl = newcl; Curve Loop(lb1_cl) = {-l20,+l44,+l24,-l41};
lb2_cl = newcl; Curve Loop(lb2_cl) = {-l19,+l43,+l23,-l44};
// UP
uh_cl = newcl; Curve Loop(uh_cl) = {l25, l26, l27, l28};
uc_cl = newcl; Curve Loop(uc_cl) = {l29, l30, l31, l32};
ut_cl = newcl; Curve Loop(ut_cl) = {+l25, +l26, +l47, -l30, -l29, -l45};
ub_cl = newcl; Curve Loop(ub_cl) = {-l28, -l27, +l47, +l31, +l32, -l45};
ut1_cl = newcl; Curve Loop(ut1_cl) = {+l25,+l46,-l29,-l45};
ut2_cl = newcl; Curve Loop(ut2_cl) = {+l26,+l47,-l30,-l46};
ub1_cl = newcl; Curve Loop(ub1_cl) = {-l28,+l48,+l32,-l45};
ub2_cl = newcl; Curve Loop(ub2_cl) = {-l27,+l47,+l31,-l48};
// Surfaces
// RIGHT VANE
rh_s = news; Surface(rh_s) = {rh_cl};
rc_s = news; Surface(rc_s) = {rc_cl};
rt1_s = news; Surface(rt1_s) = {rt1_cl};
rb1_s = news; Surface(rb1_s) = {rb1_cl};
rt2_s = news; Surface(rt2_s) = {rt2_cl};
rb2_s = news; Surface(rb2_s) = {rb2_cl};
// DOWN VANE
dh_s = news; Surface(dh_s) = {dh_cl};
dc_s = news; Surface(dc_s) = {dc_cl};
dt1_s = news; Surface(dt1_s) = {dt1_cl};
db1_s = news; Surface(db1_s) = {db1_cl};
dt2_s = news; Surface(dt2_s) = {dt2_cl};
db2_s = news; Surface(db2_s) = {db2_cl};
// LEFT VANE
lh_s = news; Surface(lh_s) = {lh_cl};
lc_s = news; Surface(lc_s) = {lc_cl};
lt1_s = news; Surface(lt1_s) = {lt1_cl};
lb1_s = news; Surface(lb1_s) = {lb1_cl};
lt2_s = news; Surface(lt2_s) = {lt2_cl};
lb2_s = news; Surface(lb2_s) = {lb2_cl};
// UP VANE
uh_s = news; Surface(uh_s) = {uh_cl};
uc_s = news; Surface(uc_s) = {uc_cl};
ut1_s = news; Surface(ut1_s) = {ut1_cl};
ub1_s = news; Surface(ub1_s) = {ub1_cl};
ut2_s = news; Surface(ut2_s) = {ut2_cl};
ub2_s = news; Surface(ub2_s) = {ub2_cl};
//
// Physical Surfaces
// RIGHT VANE
Physical Surface("RightHot",188) = {rh_s};
Physical Surface("RightCold",183) = {rc_s};
Physical Surface("RightTop",1820) = {rt1_s, rt2_s};
Physical Surface("RightBottom",182) = {rb1_s, rb2_s};


// DOWN VANE
Physical Surface("DownHot", 48) = {dh_s};
Physical Surface("DownCold",43) = {dc_s};
Physical Surface("DownTop",420) = {dt1_s, dt2_s};
Physical Surface("DownBottom",42) = {db1_s, db2_s};
// LEFT VANE
Physical Surface("LeftHot",128) = {lh_s};
Physical Surface("LeftCold",123) = {lc_s};
Physical Surface("LeftTop",1220) = {lt1_s, lt2_s};
Physical Surface("LeftBottom",122) = {lb1_s, lb2_s};
// UP VANE
Physical Surface("UpHot", 218) = {uh_s};
Physical Surface("UpCold",213) = {uc_s};
Physical Surface("UpTop",2120) = {ut1_s, ut2_s};
Physical Surface("UpBottom",212) = {ub1_s, ub2_s};
//+
r_sl = newsl; Surface Loop(r_sl) = {rh_s, rc_s, rt1_s, rt2_s, rb1_s, rb2_s};
d_sl = newsl; Surface Loop(d_sl) = {dh_s, dc_s, dt1_s, dt2_s, db1_s, db2_s};
l_sl = newsl; Surface Loop(l_sl) = {lh_s, lc_s, lt1_s, lt2_s, lb1_s, lb2_s};
u_sl = newsl; Surface Loop(u_sl) = {uh_s, uc_s, ut1_s, ut2_s, ub1_s, ub2_s};
Volume(1000) = {r_sl};
Volume(1001) = {d_sl};
Volume(1002) = {l_sl};
Volume(1003) = {u_sl};
//+
Sphere(900) = {0, 0, 0, radius};
Physical Surface("Wall", 23) = {99};
//
 Physical Volume ("VOLUME",  8888) = { BooleanDifference{ Volume{900}; Delete; }{ Volume{1000}; Volume{1001}; Volume{1002}; Volume{1003}; }
 };
 
//Characteristic Length {188, 183, lh_s, lc_s} = 0.01;

dokken · December 19, 2023, 4:25pm

19ec94:

Upon investigation the length of facet_tags is not equal to the number of total cells in the mesh. It only returns the non-zero id, which on the first look makes sense.
print(len(cell_tags.indices))
print(len(facet_tags.indices))
print(mesh_dfx.topology.index_map(mesh_dfx.topology.dim).size_local)
In short, I want to use facet_tags just as I have used boundaries in old dolfin. Could anyone confirm if the boundaries in dolfin and facet_tags in dolfinx mentioned above in my explanation are equivalent?

As GMSH only defines tags for the facets that are at the exterior interfaces, there are the only ones that are read into DOLFINx. Then, instead of saving a value for every facet, the dolfinx.mesh.MeshTags-object only stores those of interest.

This means that the meshtag-object will work in the same way as the MeshFunction in legacy Dolfin, just storing less data.

19ec94 · December 19, 2023, 4:58pm

I see! Your response is very reassuring. Indeed, after creating this post, it occurred to me. So, I counted the “non-zero” entries from the old dolfin boundaries array. The length of it is the same as facet_tags array, which is a dolfinx.mesh.MeshTags object. @dokken Thanks for the explanation.

Topic		Replies	Views
Create mesh with gmsh, convert it and import it in dolfin, problems mesh	5	10801	April 1, 2019
Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio mesh	273	42726	March 4, 2024
Dolfin-convert for convert msh to xml file mesh	1	1140	May 3, 2020
Read meshtag from xdmf created by <mesh generation with gmsh> I/O mesh , dolfinx , gmsh	2	45	July 20, 2024
Problem converting file.msh (from gmsh) to file.xml I/O	12	4382	February 5, 2019

Regarding equivalent representation of dolfin.Mesh module expressions in dolfinx.gmshio

Related topics