Meshing issue- likely in how I'm constructing my meshes in GMSH

mesh
1

Hello!

I’m having some difficulty getting linear systems to converge (or simply be accepted as valid) in a mixed FE space when I use a mesh I’ve created in GMSH. (When I use a UnitCubeMesh my code works just fine.)

My goal is to create a unit cube mesh with tight discretizations near the six faces, which relax to a coarser discretization as you move further away.

In 2D problems, I have little trouble getting a method like this working by creating Surfaces within Surfaces, and specifying mesh widths as I please. In 3D, I presumed a Volume within a Volume would work the same way.

I’ll present my GMSH code below, but if a user has an idea that would satisfy my requirements without feeling the need to get into the weeds here, I’d be happy to try it.

iSpacing = 1.0/8;
oSpacing = 1.0/20;
th = 3*oSpacing;

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



Point(13) = {(1-th), th,     th,     iSpacing};
Point(14) = {(1-th), (1-th), th,     iSpacing};
Point(15) = {th,     (1-th), th,     iSpacing};
Point(16) = {th,     th,     (1-th), iSpacing};
Point(17) = {(1-th), th,     (1-th), iSpacing};
Point(18) = {(1-th), (1-th), (1-th), iSpacing};
Point(19) = {th,     (1-th), (1-th), iSpacing};
Point(20) = {th,     th,     th,     iSpacing};
Line(13) = {19, 18};
Line(14) = {18, 17};
Line(15) = {17, 13};
Line(16) = {13, 20};
Line(17) = {20, 15};
Line(18) = {15, 19};
Line(19) = {19, 16};
Line(20) = {16, 20};
Line(21) = {16, 17};
Line(22) = {14, 13};
Line(23) = {14, 18};
Line(24) = {14, 15};

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

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

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

Volume(1) = {1,2};
Volume(2) = {2};

When I create this mesh in gmsh, and export to .msh (v2 or v4 ascii) I can get meshio-convert to create an xml file, which dolfin can still read. I know this format isn’t preferred anymore, but I’m actually having difficulty finding a functioning code to replace it. I’ve found many posts by dokken since 2019 emphatically pleading with users to adopt the new way, and I’ve have tried many of them (as they’ve changed here and there with context and package updates since 2019) but none of the ones I’ve found manage to save, and read an xdmf mesh without some variety of errors (in my case.)

My current builds are all, presumably, the latest given that I’m working with a fresh install as of 48hrs ago.

Here is a selection of code that fails

path_to_mesh = r"/home/sean/Dropbox/Research/3D Mesh/Unit Cubes/nested_cubes_8_20.xml"
datapath = r"/home/sean/Dropbox/Research/BLM/"


mesh = Mesh(path_to_mesh)
#mesh = UnitCubeMesh(15,15,15)

K = 3
vspace = 'P'
pspace = 'P' 

VELO = VectorElement(vspace, tetrahedron, K, dim=3)
PRES = FiniteElement(pspace, tetrahedron, K-1)

V    = FunctionSpace(mesh,MIXD) 
Vec  = FunctionSpace(mesh,VELO) 
Pre  = FunctionSpace(mesh,PRES)

There error I currently get is:

ufl.log.UFLException: Non-matching cell of finite element and domain.

After creating other (similarly prescribed) meshes, sometimes I’ll get past this error and instead just get a divergence error from the solver. None of this happens with a uniform mesh, so I’m confident the problem lies in either my mesh itself, or how I’m getting that mesh into dolfin.

As always- thank you all for being so helpful!

2

Hi @Sean_Breckling,

I believe the issue is the physical groups and the default behavior of Gmsh when physical groups are defined:

Since you have defined a physical group of surfaces, but not of volumes, your volume elements are not being exported. Hence the error when you try to define a tetrahedron finite element, since your mesh contains only triangles. To fix this, you can add Physical Volume(1) = {1,2}; to the end of your .geo file.

To convert the mesh, you can use the meshio Python package (note: not the meshio-convert script). Copy the following to a .py file (e.g. convert_mesh.py) and run using python convert_mesh.py:

import meshio

mesh_name = "mesh"

msh = meshio.read(mesh_name+".msh")

tet_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
meshio.write(mesh_name+".xdmf",
    meshio.Mesh(points=msh.points,
        cells={"tetra": msh.cells_dict["tetra"]},
        cell_data={"dom_marker": [tet_data]}
    )
)

tri_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
meshio.write(mesh_name+"_surf.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"triangle": msh.cells_dict["triangle"]},
        cell_data={"bnd_marker": [tri_data]}
    )
)

This will give you two files, one containing the domain markers and one containing the surface markers. To read the mesh, use

from fenics import *

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

P.S. There’s a birthday cake emoji appearing beside your username today. Happy birthday!

1 Like
3

This might be off into the weeds, but you might find Gmsh size fields to be interesting. A roughly similar mesh can be rather elegantly achieved with:

SetFactory("OpenCASCADE");

iSpacing = 1.0/8;
oSpacing = 1.0/20;
th = 3*oSpacing;

Box(1) = {0,0,0,1,1,1};

Field[1] = Distance;
Field[1].SurfacesList = {1,2,3,4,5,6};
Field[1].NumPointsPerCurve = 100;

Field[2] = Threshold;
Field[2].InField = 1;
Field[2].SizeMin = oSpacing;
Field[2].SizeMax = iSpacing;
Field[2].DistMin = 0;
Field[2].DistMax = th;

Background Field = 2;
Mesh.MeshSizeExtendFromBoundary = 0;
Mesh.MeshSizeFromPoints = 0;
Mesh.MeshSizeFromCurvature = 0;
2 Likes
4

Thank you! This worked!

1 Like