How to specify Boundary conditions when the mesh is imported from GMSH?

[Pratik]

Hi all,
I have created a 2d mesh with triangular elements in GMSH and have imported it in FEniCS.
There are 4 boundaries in my mesh i.e. free, fixed, top and bottom. These names are given in the .geo file.
I want to specify Dirichlet BC for fixed, Robin BC for free, and Neuman BC for top and bottom boundary.
I don’t want to specify BCs as it given in tutorials eg.

# Defining Dirichlet boundary condition
tol = 1E-14 # tolerance for coordinate comparisons

def Dirichlet_boundary0(x, on_boundary):
    return on_boundary and abs(x[0]) < tol

bc = DirichletBC(V, Constant(0), Dirichlet_boundary0)

Because this way will be limited to simple problems, it won’t work if a boundary is a circle or any other shape.
I want to specify the BC by using the physical groups created in the GMSH i.e. fixed, free, top and bottom. For this, I will have to create subdomains for each boundary.
I am attaching the .geo file and the code here.

// Gmsh project created on Wed Jun 16 11:22:49 2021
SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 3.0};
//+
Point(2) = {2, 0, 0, 3.0};
//+
Point(3) = {2, 1, 0, 3.0};
//+
Point(4) = {0, 1, 0, 3.0};
//+
Line(1) = {4, 1};
//+
Line(2) = {1, 2};
//+
Line(3) = {2, 3};
//+
Line(4) = {3, 4};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Physical Curve("fixed", 5) = {1};
//+
Physical Curve("free", 6) = {3};
//+
Physical Curve("top", 7) = {4};
//+
Physical Curve("bottom", 8) = {2};
//+
Physical Surface("fluid", 9) = {1};

This is the code.

from fenics import *
import matplotlib.pyplot as plt

pip install --user meshio
pip install --user h5py
import h5py
import meshio
mesh_from_file = meshio.read("loaded_bar.msh")

import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 2)

Thanks

[dokken]

As you have written the boundary data to xdmf, you can read it in as described in:

Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

This can be done using meshio and writing the explicit data to xdmf file. import meshio msh = meshio.read("test.msh") meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]})) meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]}, cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}})) from dolfin import * mesh = Mesh() with XDMFFile("mesh.xdmf") …

and use the mesh function in the boundary condition as described in

Boundary markers for moving mesh

Here is an example of what you would like to do: from dolfin import * mesh = UnitSquareMesh(10,10) V = VectorFunctionSpace(mesh, "CG", 1) x = SpatialCoordinate(mesh) u = project(as_vector((0.05*sin(pi*x[1]), 0)), V) class LeftOnBoundary(SubDomain): def inside(self, x, on_boundary): return x[0] < DOLFIN_EPS and on_boundary mf = MeshFunction("size_t", mesh, 1, 0) LeftOnBoundary().mark(mf, 1) bc = DirichletBC(V, u, mf, 1) u_boundary = Function(V) bc.apply(u_boundary.vector()) prin…

[Pratik]

I have used markers by creating subdomains before.
eg. for left boundary

# Mark left boundary as subdomain 2
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14   # tolerance for coordinate comparisons
        return on_boundary and abs(x[0]) < tol

Gamma_0 = LeftBoundary()
Gamma_0.mark(boundary_parts, 2)

Here we are creating a subdomain and using markers for defining a boundary. But the problem here is we are specifying that my left boundary is at x=0, which I don’t want to do. Because if my left boundary is not a straight line, instead it is a circle, then I won’t be able to write x=0 for that.

I want something like wherever I have specified curves in the geometry as fixed, free, top, and bottom, I should be able to apply boundary conditions on that without specifying where is my fixed, free, top and bottom. This was my doubt.

It means roughly I want to create 4 subdomains each one for fixed, free, top, and bottom without specifying their locations. How can we do this in FEniCS?

[dokken]

Please read through my post carefully.
What I stated is:

1. If you have marked your boundary in gmsh, you can read in this data to a mesh function as described in: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #3 by dokken
2. You can now use your mesh function (mf In the link) in boundary conditions such as integral measures ds_C = Measure(“ds”, domain=mesh, subdomain_marker=mf) or in Dirichlet conditions bc = DirichletBC(V, u, mf, 1)

[Pratik]

It worked for me now. Thank you very much.
I just have some small doubts.

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

In this line, what is the meaning of ‘size_t’ and ‘mesh.topology().dim-1’
and

infile.read(mvc, "name_to_read")

In this line, what is the meaning of ‘name_to_read’

Thanks

[dokken]

@Pratik

In this line, what is the meaning of ‘size_t’ and ‘mesh.topology().dim-1’

As the core code in dolfin is C++ “size_t” specifies what kind of input type the values in your mesh function have.

mesh.topology().dim()-1 specifies that the input meshfunction is for facets. If this value is set to 0, it would be reading in vertex data, or if set to mesh.topology().dim() it would assume that the data is cell data.

@Pratik

In this line, what is the meaning of ‘name_to_read’

When you converted your mesh and facet mesh to xdmf

@Pratik

out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})

you wrote the data in a dictionary with key “name_to_read”, therefore you need to specify the same key when reading in the data

[lx1998]

What is the corresponding meshio version and dolfin version？

[Pratik]

I am using FEniCS version 2019.1.0 and meshio version 4.4.6.
My issue got solved.

[Notting_X]

hello dokken, i have some problems
My mesh is imported from gmsh into fenics and converted to xdmf format, which is a 2D mesh.
For the line mvc = MeshValueCollection(“size_t”, mesh, 2), for 2D meshes, the parameter dim of meshvaluecollection is 1 or 2

[dokken]

Without a minimal reproducible example I cannot give you any guidance.

I also do not understand if the following is a question or a statement

@Notting_X

For the line mvc = MeshValueCollection(“size_t”, mesh, 2), for 2D meshes, the parameter dim of meshvaluecollection is 1 or 2

[Notting_X]

that’s a question.
I made a 2D unstructured mesh in gmsh and utilized xdmf to read the mesh into fenics,

with fe.XDMFFile(‘./cylinder_mesh.xdmf’) as infile:
infile.read(mesh)

mvc=fe.MeshValueCollection(‘size_t’,mesh,2)
with fe.XDMFFile(‘./cylinder_mf.xdmf’) as infile:
infile.read(mvc,‘name_to_read’)

In the line mvc=fe.MeshValueCollection(‘size_t’,mesh,2), is 2 the parameter dim? If so, should the parameter dim be 1 or 2?

[dokken]

If you are reading in cell markers, it should be 2, which is the topological dimension of the cell.
You can use mesh.topology().dim().

If you are reading in facet markers, you should use 1, or mesh.topology().dim()-1

[Notting_X]

This is a 2D mesh. I want to read edge markers on the top, bottom, left and right sides to define boundary conditions on all four sides

[dokken]

Edges have the topological dimension 1, and corresponds to facets when the mesh is 2D.

[Notting_X]

I also have the other question.
I have tried to modify the code in the “Navier-Stokes for flow around a cylinder” grid section of the fenics tutorial to use the grids generated by gmsh to replace the grids generated by fenics. When T is greater than 5 and num_steps is greater than 5000, the error “Output the size limit. Open the full output data in a text editor” is displayed in the 4984th step.

[dokken]

The error seems to be related to printing/outputting.
But without the full code that can reproduce this it is really hard to give any suggestions as to what to change.

[Notting_X]

Thanks for your answers, I will check the printing and output errors