Creating subdomains and boundary conditions on imported gmsh mesh

General
1

hi all,

i was wondering if there is any documentation on how to implement subdomains boundary conditions on imported gmsh files (i tried to search by myself but didnt found anything that is suitable for this )

i have shared my gmsh file before, will share again
i need to create few different subdomains in order to set different material properties for each subdoamin (cylinder and box that wraps the cylinder) and of course some boundary conditions on the facets

im using fenics and not fencisx

will be happy to know if anyone has found any doc for something similar to that

// Gmsh project created on Mon Jun 03 18:05:16 2024
//+
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1, 0.2, 0.2};
//+
Cylinder(2) = {-0.1, 0.1, 0.1, 1.5, 0, 0, 0.05, 2*Pi};
//+
BooleanDifference{ Volume{1}; Delete; }{ Volume{2}; Delete; }
//+
Cylinder(2) = {0, 0.1, 0.1, 1, 0, 0, 0.049, 2*Pi};
//+
Physical Surface("left", 19) = {1};
//+
Physical Surface("right", 20) = {7};
//+
Physical Surface("front", 21) = {3};
//+
Physical Surface("back", 22) = {4};
//+
Physical Surface("top", 23) = {5};
//+
Physical Surface("down", 24) = {2};
//+
Physical Volume("box", 25) = {1};
//+
Physical Volume("cylinder", 55) = {2};
2

I might be misunderstanding the question, but Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio (that you probably know already) should have everything you need.

The subdomain and boundary markers are those that you set in Physical Surface and Physical Volume.

3

hi thanks for your reply

would you say that this will be appropriate way to set my boundary conditions and subdomains according to the geo i attached above?
i want the cylinder to be one subdoamin and the box to be another

class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return (pow((x[1]-0.1),2)+pow((x[2]-0.1),2)) <= pow(0.05,2)
class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return (pow((x[1]-0.1),2)+pow((x[2]-0.1),2)) > pow(0.05,2)

subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(0)
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)

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

left = AutoSubDomain(lambda x: near(x[0], 0))
right = AutoSubDomain(lambda x: near(x[0], 1))
up = AutoSubDomain(lambda x: near(x[2], 0.2))
down = AutoSubDomain(lambda x: near(x[2], 0))
front = AutoSubDomain(lambda x: near(x[1], 0))
back = AutoSubDomain(lambda x: near(x[1], 0.2))

left.mark(boundary_parts, 1)
right.mark(boundary_parts, 2)
up.mark(boundary_parts, 5)
down.mark(boundary_parts, 6)
front.mark(boundary_parts, 3)
back.mark(boundary_parts, 4)
4

You shouldn’t need any of that: the link I posted gives you a way to save the subdomain and boundary markers from gmsh to file, if you read them in within a MeshFunction you’ll have them ready without any need of recomputing them.

1 Like
5

is this enough in order to load the cell functions ?


mesh = meshio.read("my_perfect_model.msh")

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

tetra_mesh = create_mesh(mesh, "tetra")
triangle_mesh = create_mesh(mesh, "triangle")
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)

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



#creating subtomains and boundaty conditions



materials = MeshFunction('size_t', mesh, mesh.topology().dim())


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

bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 1)
bc3 = DirichletBC(V.sub(1), Constant(0), boundary_parts, 3)
bc5 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 5)
bc6 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 2)
tDirBC = Expression(('1.0*time_'), time_=0.0, degree=0)
tDirBC4 = Expression(('((1/(time_+1))-1)*0.2'), time_=0.0, degree=0)
bc2 = DirichletBC(V.sub(0), tDirBC, boundary_parts, 7)
bc4 = DirichletBC(V.sub(1), tDirBC4, boundary_parts, 4)

bcs = [bc1, bc2, bc3, bc4, bc5, bc6]
subdomain_ids = [7, 4]
ds = Measure("ds", subdomain_data=boundary_parts, subdomain_id=(7,4))
ds = ds(degree=4)

or do am i missing something? cause when im running my code i get the following notification , i dont have any reference here for the volumes i have created, i dont know how to apply them here, basically the subdomains (volume) should be in the materials meshfunction

  *** Warning: Found no facets matching domain for boundary condition.
6

Hopefully the following will clarify what you need to do. The information was already there in the post I linked, but I recognize that the topic is long and hence it may have been difficult for you to find it. However, in future please try to make an effort in reading previous posts and, if that fails, please put others in the best possible position to help you (for instance: your code didn’t have any import).

Start with the following file. Let’s call it convert_mesh.py.

import os

import dolfin
import meshio
import numpy as np

mesh_name = "my_perfect_model"

# Then, read it in through meshio
meshio_mesh = meshio.read(f"{mesh_name}.msh")

# Extract mesh corresponding to each cell type, and write it out with meshio
def extract_mesh_by_cell_type(meshio_mesh, cell_type):
    cells = meshio_mesh.get_cells_type(cell_type)
    cell_data = meshio_mesh.get_cell_data("gmsh:physical", cell_type)
    return meshio.Mesh(points=meshio_mesh.points, cells={cell_type: cells}, cell_data={"markers": [cell_data]})

tetra_mesh = extract_mesh_by_cell_type(meshio_mesh, "tetra")
meshio.write(f"{mesh_name}_tetra.xdmf", tetra_mesh)
del tetra_mesh
triangle_mesh = extract_mesh_by_cell_type(meshio_mesh, "triangle")
meshio.write(f"{mesh_name}_triangle.xdmf", triangle_mesh)
del triangle_mesh

del meshio_mesh

# Read back in the meshio mesh from XDMF in order to convert to a dolfin mesh
dolfin_mesh = dolfin.Mesh()
dolfin.XDMFFile(f"{mesh_name}_tetra.xdmf").read(dolfin_mesh)

# Read back in the markers in order to attach them to the dolfin mesh
def read_back_in_markers(meshio_xdmf_file, dolfin_mesh):
    dolfin_markers = dolfin.MeshValueCollection("size_t", dolfin_mesh)
    dolfin.XDMFFile(meshio_xdmf_file).read(dolfin_markers, "markers")
    dolfin_mesh_function = dolfin.cpp.mesh.MeshFunctionSizet(dolfin_mesh, dolfin_markers)
    # Note that the values of the MeshFunction are undefined if the entity (cell or facet)
    # is not listed in the MeshValueCollection. In any implementation I've seen this result
    # in sys.maxsize being stored for all unlisted entities: a simple workaround is to replace
    # undefined values with zero
    dolfin_mesh_function_array = dolfin_mesh_function.array()
    defined_markers = np.array(list(set(dolfin_markers.values().values())), dtype=dolfin_mesh_function_array.dtype)
    defined_entries = np.isin(dolfin_mesh_function_array, defined_markers)
    dolfin_mesh_function_array[~defined_entries] = 0
    dolfin_mesh_function.set_values(dolfin_mesh_function_array)
    return dolfin_mesh_function

dolfin_subdomains = read_back_in_markers(f"{mesh_name}_tetra.xdmf", dolfin_mesh)
dolfin_boundaries = read_back_in_markers(f"{mesh_name}_triangle.xdmf", dolfin_mesh)

# Now export them to file
dolfin.XDMFFile(f"{mesh_name}.xdmf").write(dolfin_mesh)
dolfin.XDMFFile(f"{mesh_name}_subdomains.xdmf").write(dolfin_subdomains)
dolfin.XDMFFile(f"{mesh_name}_boundaries.xdmf").write(dolfin_boundaries)

# Please now open this three files in paraview and check that the labels are as expected

# You may remove *_tetra.* and *_triangle.*, they are not needed anymore
os.remove(f"{mesh_name}_tetra.xdmf")
os.remove(f"{mesh_name}_tetra.h5")
os.remove(f"{mesh_name}_triangle.xdmf")
os.remove(f"{mesh_name}_triangle.h5")

Once you have converted the mesh, there is no need to copy and paste the content of the file above in your application! Simply read the mesh back in. A sample file that does that is

import dolfin
import numpy as np

mesh_name = "my_perfect_model"

mesh = dolfin.Mesh()
dolfin.XDMFFile(f"{mesh_name}.xdmf").read(mesh)

subdomains = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim())
dolfin.XDMFFile(f"{mesh_name}_subdomains.xdmf").read(subdomains)
print(f"Subdomains: {np.min(subdomains.array())}, {np.max(subdomains.array())}")

boundaries = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
dolfin.XDMFFile(f"{mesh_name}_boundaries.xdmf").read(boundaries)
print(f"Boundaries: {np.min(boundaries.array())}, {np.min(boundaries.array()[boundaries.array() > 0])}, {np.max(boundaries.array())}")
7

thank you so much!
but i was wondering why do i get only 2 boundaries if i have more surfaces (6) , only 19 and 24
i can see in the input only 2 , and also i can see 0 but i do not have surface marked as zero
the printed output

Subdomains: 25, 55
Boundaries: 0, 19, 24
8

why do i get only 2 boundaries

the final line of the second code should be self-explanatory on why this happens

also i can see 0 but i do not have surface marked as zero

I put a lot of effort in commenting the codes in the first snippet. Please read those comments and answer your own question.

9

hi , regarding this

i saw that you extracted the min and max boundary value and this is why i see only two values, but even when im trying to print the whole boundaries i can see only 19 and 24 and the others are zero. now i have read your comment that every undefined boundary will be set as zero, but my question is why the other physical surfaces are not defined? why dont i see 20,21 etc
also, when i tried to look for the size of the boundaries i saw that it has over 1600 elements, can you pls share why does it count over 1600 elements? maybe im missing something in the understanding.

appreciate your time for helping me

10

Can’t reproduce what you are saying.

Slightly modifying the second snippet

import dolfin
import numpy as np

mesh_name = "my_perfect_model"

mesh = dolfin.Mesh()
dolfin.XDMFFile(f"{mesh_name}.xdmf").read(mesh)

subdomains = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim())
dolfin.XDMFFile(f"{mesh_name}_subdomains.xdmf").read(subdomains)
print(f"Subdomains: {np.unique(subdomains.array())}")

boundaries = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
dolfin.XDMFFile(f"{mesh_name}_boundaries.xdmf").read(boundaries)
print(f"Boundaries: {np.unique(boundaries.array())}")

I get

Subdomains: [25 55]
Boundaries: [ 0 19 20 21 22 23 24]

when i tried to look for the size of the boundaries i saw that it has over 1600 elements, can you pls share why does it count over 1600 elements

1600 doesn’t mean anything to me. I had to generate my own mesh based on your geo file, so the number of elements in my mesh is definitely going to be different from yours. Anyways, the number of entries you see in the boundary markers is the number of faces (not cells! triangles!), so you shouldn’t compare it to the number of tetras.

11

Hi,

I have published the Gmsh and FEniCS codes to do the same. You can find the codes on the following Github link.

Manish

1 Like
13

Hi,

The codes and mesh to done the requested task are available on the following link.

https://drive.google.com/drive/folders/1IWHsS-g1GrKfIncvjasdtOk2xWitD-bB?usp=drive_link

Manish

1 Like
14

hi @manish_07 thank you so much for your time and help

additional question, if i want to create boundary conditions
would this way be an appropriate way to do it using the xml file you attached and the geo?


import os
import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
from ffc.fiatinterface import create_quadrature as cquad
import numpy as np
import meshio
import dolfin 

boundary_parts = MeshFunction('size_t', mesh, "1_facet_region.xml")

bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 19)
bc3 = DirichletBC(V.sub(1), Constant(0), boundary_parts,21)
bc5 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 23)
bc6 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 24)
tDirBC = Expression(('1.0*time_'), time_=0.0, degree=0)
tDirBC4 = Expression(('((1/(time_+1))-1)*0.2'), time_=0.0, degree=0)
bc2 = DirichletBC(V.sub(0), tDirBC, boundary_parts, 20)
bc4 = DirichletBC(V.sub(1), tDirBC4, boundary_parts, 22)

bcs = [bc1, bc2, bc3, bc4, bc5, bc6]
subdomain_ids = [20,22]
ds = Measure("ds", subdomain_data=boundary_parts, subdomain_id=(20,22))
15

hi @Mirialex yes the method you used to apply the boundary conditions is the correct.

Manish