Help Needed: Correctly Imposing Dirichlet Boundary Conditions with GMSH and Meshio in FEniCSx

General
1

Hello everyone,

I am currently working on solving a Poisson problem using FEniCSx, and I am encountering issues with correctly imposing Dirichlet boundary conditions. Despite following several tutorials, I find that the Dirichlet conditions are being applied to nodes that do not correspond to the boundaries of the hole or the plate in my mesh.

Here is a brief overview of my setup:

  • I generate a mesh using GMSH, which includes a square plate with a circular hole in the center.
  • I then read the mesh using Meshio and attempt to apply Dirichlet boundary conditions using FEniCSx.

Specifically, I need to impose a Dirichlet condition of u = 5 on the outer boundary of the plate and u = -1 on the boundary of the hole. However, these conditions are not being correctly applied to the intended nodes.

I suspect that there might be an issue with how I am reading the mesh or locating the degrees of freedom for the boundary conditions. I can provide the .msh file or the code used to generate the mesh if needed.

Here are the relevant parts of my code:

domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=2)

…

Define Dirichlet boundary conditions

u_bc_outer = fem.Function(V)
u_bc_outer.interpolate(lambda x: np.full_like(x[0], 5))

u_bc_hole = fem.Function(V)
u_bc_hole.interpolate(lambda x: np.full_like(x[0], -1))

Apply Dirichlet conditions on the outer boundary

outer_boundary_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isin(facet_markers.values, [3]))
bc_outer = fem.dirichletbc(u_bc_outer, outer_boundary_dofs)

Apply Dirichlet conditions on the hole boundary

hole_boundary_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isin(facet_markers.values, [2]))
bc_hole = fem.dirichletbc(u_bc_hole, hole_boundary_dofs)

I would greatly appreciate any insights or suggestions on how to correctly impose these boundary conditions. If needed, I can attach the .msh file or provide more details about my mesh generation process.

Thank you very much for your assistance!

Best regards,
Diego.

2

Please read Read before posting: How do I get my question answered?

In particular, your code is not properly formatted (not included in a “```” block), and not reproducible.

The general suggestion is not to use meshio, and rely on dolfinx.io.gmshio.read_from_msh.

3

Thank you Francesco!
This is how I construct the mesh:

import gmsh
import np
from mpi4py import MPI
from dolfinx.io import gmshio

# Initialize GMSH
gmsh.initialize()

# Parameters for the plate and the hole
L = 1.0  # Length of the plate side
H = 1.0  # Height of the plate (equal to the length for a square plate)
r = 0.2  # Radius of the hole
c_x = L / 2  # x-coordinate of the center of the hole
c_y = H / 2  # y-coordinate of the center of the hole
gdim = 2  # Dimension of the problem (2D)

# Create the plate and the hole
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H)
hole = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

# Subtract the hole from the plate
plate = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, hole)])
gmsh.model.occ.synchronize()

# Add the plate volume as a physical group
plate_marker = 1
volumes = gmsh.model.getEntities(dim=gdim)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], plate_marker)
gmsh.model.setPhysicalName(volumes[0][0], plate_marker, "Plate")

# Add the outer boundary and hole boundary as physical groups
hole_marker = 2
border_marker = 3
outer_boundary, hole_boundary = [], []

boundaries = gmsh.model.getBoundary(volumes, oriented=False)

for boundary in boundaries:
    center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
    if np.allclose(center_of_mass, [c_x, c_y, 0], atol=r):
        hole_boundary.append(boundary[1])
    else:
        outer_boundary.append(boundary[1])

if outer_boundary:
    gmsh.model.addPhysicalGroup(1, outer_boundary, border_marker)
    gmsh.model.setPhysicalName(1, border_marker, "OuterBoundary")

if hole_boundary:
    gmsh.model.addPhysicalGroup(1, hole_boundary, hole_marker)
    gmsh.model.setPhysicalName(1, hole_marker, "HoleBoundary")

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

# Save the mesh to a .msh file
gmsh.write("square_with_hole.msh")

# Finalize GMSH
gmsh.finalize()

And then, following your idea:

import numpy as np
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx import mesh, fem
from dolfinx.io import XDMFFile
from ufl import TrialFunction, TestFunction, inner, grad, dx
from dolfinx.fem.petsc import LinearProblem


# MPI communicator

comm = MPI.COMM_WORLD

# Read the mesh from the .msh file

filename = "square_with_hole.msh"

domain, cell_tags, facet_tags = gmshio.read_from_msh(filename, comm, gdim=2)
# Create function space
V = fem.functionspace(domain, ("CG", 1))

# Define Dirichlet boundary conditions
u_bc_outer = fem.Function(V)
u_bc_outer.interpolate(lambda x: np.full_like(x[0], 5))

u_bc_hole = fem.Function(V)
u_bc_hole.interpolate(lambda x: np.full_like(x[0], -1))

# Apply Dirichlet conditions on the outer boundary
outer_boundary_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isin(facet_tags.values, [3]))
bc_outer = fem.dirichletbc(u_bc_outer, outer_boundary_dofs)

# Apply Dirichlet conditions on the hole boundary
hole_boundary_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isin(facet_tags.values, [2]))
bc_hole = fem.dirichletbc(u_bc_hole, hole_boundary_dofs)

# Define the variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = fem.Constant(domain, 1.0)  # Source term

a = inner(grad(u), grad(v)) * dx
L = f * v * dx

# Solve the system
u = fem.Function(V)
problem = LinearProblem(a, L, bcs=[bc_outer, bc_hole])
uh = problem.solve()

# Save the solution to an .xdmf file
with XDMFFile(MPI.COMM_WORLD, "solution.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(uh)

The code runs and I have a solution but it’s incorrect, since the bcs are not applied to the correct dofs (see the picture). Any help is welcome :slight_smile:
Imagen1

4

If you intend to use the tags to set the BCs then you must use fem.locate_dofs_topological, not fem.locate_dofs_geometrical, which is supposed to be used when the boundary has an analytical expression (e.g, you provide the analytical expression of a side of your rectangle). Please see the demos, or search this forum, for plenty of examples on how to use that function.

1 Like
5

Thank you, Francesco! That was the problem. With your guidance I’ve implemented this step (with good results) as:

# Apply Dirichlet conditions on the outer boundary
dofs_out = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_tags.find(3))
bc_out = fem.dirichletbc(value=10., dofs=dofs_out, V=V)

# Apply Dirichlet conditions on the hole boundary
dofs_hole = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_tags.find(2))
bc_hole = fem.dirichletbc(value=3., dofs=dofs_hole, V=V)

Using facet_tags.find() is the best way?

6

I personally use facet_tags.indices[facet_tags.values == 1] rather than facet_tags.find(1), but they ought to be equivalent.

1 Like