How do I tag and use boundary surfaces from gmsh correctly?

TomM · November 5, 2024, 1:14pm

Hi all,
I think there is something wrong with either how I define the boundary surfaces or with how I access them.
I Initialize gmsh with a horizontal cylinder with circular faces on the yz-plane in 3D.

This is how I get the boundary surfaces (left and right circular faces, and the lateral surface):

surfaces = gmsh.model.getBoundary([(3, cylinder_tag)], oriented=False, recursive=False)

# Right face
right_face_tag = surfaces[0][1]
gmsh.model.addPhysicalGroup(gdim - 1, [right_face_tag], 1)
gmsh.model.setPhysicalName(gdim - 1, 1, "Right")

# Skin surface
lateral_surface_tags = [surf[1] for surf in surfaces[2:]]
gmsh.model.addPhysicalGroup(gdim - 1, lateral_surface_tags, 3)
gmsh.model.setPhysicalName(gdim - 1, 3, "Lateral")

domain, _, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=gdim)

dx = Measure("dx", domain=domain)   
ds = Measure("ds", domain=domain, subdomain_data=facets)  

a = (
    inner(sigma(u), epsilon(v)) * dx +
    (1 / dt) * inner(u, v) * dx +
    friction_coefficient * inner(u, v) * ds(3) # friction at the skin of the cylinder
)

# force is applied at the right face
L = inner(f, v) * ds(1) + (1 / dt) * inner(u_n, v) * dx

When trying to apply the force at the right face using ds(1) nothing happens. If I instead apply it in the whole domain using dx the whole body moves (like expected). I also think ds(3) for friction at the skin is not working as I want either.

For this I worked with Thermo-elastic evolution problem (full coupling) — Computational Mechanics Numerical Tours with FEniCSx
and Diffusion of a Gaussian function — FEniCSx tutorial.

Thanks for any help or guidance in the right direction.

dokken · November 5, 2024, 2:22pm

TomM:

surfaces = gmsh.model.getBoundary([(3, cylinder_tag)], oriented=False, recursive=False)

# Right face
right_face_tag = surfaces[0][1]
gmsh.model.addPhysicalGroup(gdim - 1, [right_face_tag], 1)
gmsh.model.setPhysicalName(gdim - 1, 1, "Right")

# Skin surface
lateral_surface_tags = [surf[1] for surf in surfaces[2:]]
gmsh.model.addPhysicalGroup(gdim - 1, lateral_surface_tags, 3)
gmsh.model.setPhysicalName(gdim - 1, 3, "Lateral")

domain, _, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=gdim

have you visualized facets in Paraview, to check that they are marked as expected? Could you provide a figure with the facets, i.e.

with dolfinx.io.XDMFFile(domain.comm, "facets.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facets, domain.geometry)

and open that in Paraview and take a screenshot?

Note that we cannot give you any further help as you have not provided a minimal reproducible example.

TomM · November 5, 2024, 2:54pm

Here is my full code simplified as much as I could:

import numpy as np  # type: ignore
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction  # type: ignore
from mpi4py import MPI  # type: ignore
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from petsc4py import PETSc # type: ignore
import gmsh

# Initialize Gmsh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gdim = 3

# Create a horizontal cylinder with circular faces on the yz-plane
gmsh.model.add("cylinder")
cylinder_tag = gmsh.model.occ.addCylinder(0, 0, 0, 10.0, 0, 0, 1.0)

gmsh.model.occ.synchronize()

# Retrieve boundary surfaces (left and right circular faces, and the lateral surface)
surfaces = gmsh.model.getBoundary([(3, cylinder_tag)], oriented=False, recursive=False)

# Right face
right_face_tag = surfaces[0][1] 
gmsh.model.addPhysicalGroup(gdim - 1, [right_face_tag], 1)
gmsh.model.setPhysicalName(gdim - 1, 1, "Right")

# Left face
left_face_tag = surfaces[1][1] 
gmsh.model.addPhysicalGroup(gdim - 1, [left_face_tag], 2)
gmsh.model.setPhysicalName(gdim - 1, 2, "Left")

# Lateral surface
lateral_surface_tags = [surf[1] for surf in surfaces[2:]]  
gmsh.model.addPhysicalGroup(gdim - 1, lateral_surface_tags, 3)
gmsh.model.setPhysicalName(gdim - 1, 3, "Lateral")

# Generate the 3D mesh
gmsh.model.mesh.generate(gdim)

domain, _, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=gdim)

gmsh.finalize()


shape = (gdim,) 

# Define the function space
V = fem.functionspace(domain, ("P", 2, shape))
u = TrialFunction(V)
v = TestFunction(V)

u_h = fem.Function(V, name="Displacement")

# Lamé parameters (for a given material)
E = 10.0  # Young's modulus
nu = 0.3  # Poisson's ratio

lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return lambda_ * tr(epsilon(u)) * Identity(len(u)) + 2 * mu * epsilon(u)


friction_coefficient = 0.9
f = fem.Constant(domain, np.array([10.0, 0, 0]))  

dx = Measure("dx", domain=domain)       # Volume measure
ds = Measure("ds", domain=domain, subdomain_data=facets)  # Surface measure for boundaries

bcs = []

a = (
    inner(sigma(u), epsilon(v)) * dx +
    friction_coefficient * inner(u, v) * ds(3)
)

L = inner(f, v) * ds(1)

bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form, bcs=bcs)
A.assemble()
b = create_vector(linear_form)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# define the problem
problem = fem.petsc.LinearProblem(
    a, L, u=u_h, bcs=bcs,
    petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)

vtk = io.VTKFile(domain.comm, "results/dynamic/cylinder.pvd", "w")

with io.XDMFFile(domain.comm, "facets.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facets, domain.geometry)

# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
    loc_b.set(0)
assemble_vector(b, linear_form)

# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)

# Solve linear problem
solver.solve(b, u_h.x.petsc_vec)
u_h.x.scatter_forward()

# Write solution to file
vtk.write_function(u_h, 0)

vtk.close()

And this the paraview screenshot

dokken · November 5, 2024, 3:45pm

You need to do slightly more, i.e. choose the relevant visualization of the facet marker, as you are currently just showing the mesh with no color.

Use “Extract blocks”, then choose the relevant tag from the list on the left, and set the coloring to be of that marker.

TomM · November 6, 2024, 8:34am

I am new to paraview so maybe that’s why I cant find it, but I think there are no tags. I would expect them to be under Properties where I have Bock0 selected

TomM · November 6, 2024, 9:15am

I found that the vector b which is generated from L is zero, meaning ds(1) must be 0 as well

L = inner(f, v) * ds(1)
linear_form = fem.form(L)
print("Vector b norm:", b.norm())
# returns: Vector b norm: 0.0

Edit: If I use dx instead of ds(1) it is also 0 so probably not as meaningful as I thought

dokken · November 6, 2024, 10:48am

So let’s start with the mesh generation and inspection of your tags:

from mpi4py import MPI
import dolfinx
import gmsh

# Initialize Gmsh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gdim = 3

# Create a horizontal cylinder with circular faces on the yz-plane
gmsh.model.add("cylinder")
cylinder_tag = gmsh.model.occ.addCylinder(0, 0, 0, 10.0, 0, 0, 1.0)

gmsh.model.occ.synchronize()

# Retrieve boundary surfaces (left and right circular faces, and the lateral surface)
surfaces = gmsh.model.getBoundary(
    [(3, cylinder_tag)], oriented=False, recursive=False)

# Right face
right_face_tag = surfaces[0][1]
gmsh.model.addPhysicalGroup(gdim - 1, [right_face_tag], 1)
gmsh.model.setPhysicalName(gdim - 1, 1, "Right")

# Left face
left_face_tag = surfaces[1][1]
gmsh.model.addPhysicalGroup(gdim - 1, [left_face_tag], 2)
gmsh.model.setPhysicalName(gdim - 1, 2, "Left")

# Lateral surface
lateral_surface_tags = [surf[1] for surf in surfaces[2:]]
gmsh.model.addPhysicalGroup(gdim - 1, lateral_surface_tags, 3)
gmsh.model.setPhysicalName(gdim - 1, 3, "Lateral")

# Generate the 3D mesh
gmsh.model.mesh.generate(gdim)

domain, _, facets = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, MPI.COMM_WORLD, 0, gdim=gdim)

gmsh.finalize()
with dolfinx.io.XDMFFile(domain.comm, "facets.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    domain.topology.create_connectivity(
        domain.topology.dim-1, domain.topology.dim)
    xdmf.write_meshtags(facets, domain.geometry)

which yields

then extracting the meshtags yields:

which indicates that something is up.
Considering the topological dimension of your mesh:
print(domain.topology.dim), which yields 2, indicates that only the surface mesh has been loaded.
This is because you haven’t made a Physical Volume marker in GMSH
Adding:

gmsh.model.addPhysicalGroup(gdim, [cylinder_tag], 0)

will resolve your issue and show

Topic		Replies	Views
How to select complex surface for boundary conditions mesh	11	963	February 2, 2024
How to extract and use physical labels from GMSH surfaces on FEniCS (2D mesh) mesh	29	6502	July 24, 2022
Facets orientation in inner surface with fenicsX	3	820	April 13, 2022
Creating subdomains and boundary conditions on imported gmsh mesh General	13	238	June 17, 2024
How to use gmsh tags to assign material properties? mesh	7	1591	September 13, 2022

How do I tag and use boundary surfaces from gmsh correctly?

Related Topics