Using facet tags to define boundary conditions - GMSH - FEniCSx - dolfinx

1

Suppose we have the 2D rectangular domain defined as follows:

import numpy as np
import gmsh
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI

import pyvista
from dolfinx import plot

from ufl import (FacetNormal, FiniteElement, Identity,TestFunction, TrialFunction, VectorElement,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)
from dolfinx.fem import Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
from petsc4py import PETSc


A = np.array([0,0])
B = np.array([3,0])
C = np.array([3,1])
D = np.array([0,1])

gmsh.initialize()
#gmsh.option.setNumber("General.Terminal",0) #To hide the mesh output values

mesh_size = 0.1
point_1 = gmsh.model.geo.add_point(A[0],A[1], 0, mesh_size)
point_2 = gmsh.model.geo.add_point(B[0],B[1], 0, mesh_size)
point_3 = gmsh.model.geo.add_point(C[0],C[1], 0, mesh_size)
point_4 = gmsh.model.geo.add_point(D[0],D[1], 0, mesh_size)

line_1 = gmsh.model.geo.add_line(point_1, point_2)
line_2 = gmsh.model.geo.add_line(point_2, point_3)
line_3 = gmsh.model.geo.add_line(point_3, point_4)
line_4 = gmsh.model.geo.add_line(point_4, point_1)

curve_loop    = gmsh.model.geo.add_curve_loop([line_1,line_2,line_3,line_4])
plane_surface = gmsh.model.geo.add_plane_surface([curve_loop])

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(2, [plane_surface], name = "fluid")

gmsh.model.addPhysicalGroup(1, [line_1,line_3], name = "walls")
gmsh.model.addPhysicalGroup(1, [line_4], name = "in_flow")
gmsh.model.addPhysicalGroup(1, [line_2], name = "out_flow")


gmsh.model.mesh.generate()

mesh, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0,gdim=2)

gmsh.finalize()

You can draw it via:

topology, cell_types, geometry_for_plotting = plot.create_vtk_mesh(mesh, 2)
#If you have the word geometry in place of geometry_for_plotting, it might conflict with the boundingBoxTree statement in the get coordinate function
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry_for_plotting)

pyvista.set_jupyter_backend("pythreejs")

plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()

if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    pyvista.start_xvfb()
    figure = plotter.screenshot("fundamentals_mesh.png")

Suppose to solve the Navier Stokes Poiseuille flow, we can define the boundary conditions as mentioned in the tutorial

t = 0
T = 10
num_steps = 500
dt = T/num_steps

v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

def inflow(x):
    return np.isclose(x[0], A[0])
inflow_dofs = locate_dofs_geometrical(Q, inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)

def walls(x):
    return np.logical_or(np.isclose(x[1],A[1]), np.isclose(x[1],C[1]))
wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)


def outflow(x):
    return np.isclose(x[0], C[0])
outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]

But, can you please help me to utilize the physical groups I have defined through GMSH directly; For example
gmsh.model.addPhysicalGroup(1, [line_4], name = "in_flow") for the inflow boundary condition

I’m trying as below:

fdim = mesh.topology.dim - 1
bc_inflow = dirichletbc(PETSc.ScalarType(8), locate_dofs_topological(Q, fdim, facet_tags.find(in_flow)))

I tried to use what they have done here-,To%20tag,-the%20different%20surfaces) but it was not successful.

I really appreciate your help

1 Like
2

You cannot use strings in your physical groups, as DOLFINx reads in integer markers, see for instance:
https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code2.html#mesh-generation

3

Thank you for the comment. But I’m a bit confused with the syntax. In the tutorial, they use open cascade kernal. But in my case it was .geo

So, even if I have it as:
gmsh.model.addPhysicalGroup(1, [line_4], in_flow_marker) where in_flow_marker = 100 it still gives an error. (I apologize if this is a very basic question)

For your convenience, I providing the entire code below so that it can be directly run:

A = np.array([0,0])
B = np.array([3,0])
C = np.array([3,1])
D = np.array([0,1])
in_flow_marker = 100
gmsh.initialize()
#gmsh.option.setNumber("General.Terminal",0) #To hide the mesh output values

mesh_size = 0.1
point_1 = gmsh.model.geo.add_point(A[0],A[1], 0, mesh_size)
point_2 = gmsh.model.geo.add_point(B[0],B[1], 0, mesh_size)
point_3 = gmsh.model.geo.add_point(C[0],C[1], 0, mesh_size)
point_4 = gmsh.model.geo.add_point(D[0],D[1], 0, mesh_size)

line_1 = gmsh.model.geo.add_line(point_1, point_2)
line_2 = gmsh.model.geo.add_line(point_2, point_3)
line_3 = gmsh.model.geo.add_line(point_3, point_4)
line_4 = gmsh.model.geo.add_line(point_4, point_1)

curve_loop    = gmsh.model.geo.add_curve_loop([line_1,line_2,line_3,line_4])
plane_surface = gmsh.model.geo.add_plane_surface([curve_loop])

gmsh.model.geo.synchronize()

# gmsh.model.addPhysicalGroup(1, [line_1,line_2,line_3,line_4], name = "added physical curve")

gmsh.model.addPhysicalGroup(2, [plane_surface], name = "fluid") # You need this for dolfinx

# ---- New edit --------

gmsh.model.addPhysicalGroup(1, [line_1,line_3], name = "walls")
gmsh.model.addPhysicalGroup(1, [line_4], in_flow_marker)
gmsh.model.addPhysicalGroup(1, [line_2], name = "out_flow")

# ---- end --------

gmsh.model.mesh.generate()

mesh, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0,gdim=2)

gmsh.finalize()

And for the boundary conditions:

t = 0
T = 10
num_steps = 500
dt = T/num_steps

v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

def walls(x):
    return np.logical_or(np.isclose(x[1],A[1]), np.isclose(x[1],C[1]))
wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)

def inflow(x):
    return np.isclose(x[0], A[0])
inflow_dofs = locate_dofs_geometrical(Q, inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)

fdim = mesh.topology.dim - 1
bc_inflow = dirichletbc(PETSc.ScalarType(8), locate_dofs_topological(Q, fdim, facet_tags.find(in_flow_marker)))



def outflow(x):
    return np.isclose(x[0], C[0])
outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]
4

With the first code snippet, you can inspect facet_tags.values and facet_tags.indices,
or write these to file

 with dolfinx.io.XDMFFile(mesh.comm, "ft.xdmf", "w") as xdmf:
     xdmf.write_mesh(mesh)
     xdmf.write_meshtags(facet_tags)

you can inspect the boundaries

ft
ft1280×720 9.33 KB

However, in the second part of your script (the DOLFINx code), you are using

where you should use locate_dofs_topological, as shown in
https://jorgensd.github.io/dolfinx-tutorial/chapter1/fundamentals_code.html?highlight=locate_dofs_topological#defining-the-boundary-conditions
You can for instance do:

locate_dofs_topological(Q, mesh.topology.dim-1, facet_tags.find(100))

if the inflow is at the dofs marked with 100

5

Thank you so much Dokken!