How to select complex surface for boundary conditions

Rene_Thierry_Djoumes · March 30, 2023, 9:42am

Dear all,

I am working on a simulation where my geometry is a CT-scan bowel. The geometry is in .stl and it is really difficult to select the inner or outer surface for the boundary conditions.
I have imported the geometry into Gmsh to try and create physical groups but the geometry comes in one surface and I can’t select the surfaces.
I tried to import the geometry into Freecad to get the file in .stp but the geometry arrives in Gmsh with a lot of lines and points impossible to select one by one.

Please do you have an idea on how to select surfaces in FEniCs with stl geometry?

Thanks in advance.

Ekrem_Ekici · March 31, 2023, 11:44am

Hi,

All you need to do is assigning the physical tags automatically. For example, we have a cube.

We need to import that cube into gmsh model as a step file.

import sys
import numpy as np
import gmsh

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

gmsh.model.add("cubes")
gmsh.option.setString('Geometry.OCCTargetUnit', 'M')

gmsh.model.occ.importShapes('cube.step')
gmsh.model.occ.synchronize()

gmsh.model.occ.removeAllDuplicates()
gmsh.model.occ.synchronize()

Now, we have our model in gmsh environment. We need to start assigning physical tags for surfaces and volumes. As we have only one volume (one cube), we can start with volume;

gmsh.model.addPhysicalGroup(3, [1], tag=1)

Then, we can assign tags to each surface randomly;

surfaces = gmsh.model.occ.getEntities(dim=2)

for surface in surfaces:
    gmsh.model.addPhysicalGroup(surface[0], [surface[1]])

gmsh.model.occ.synchronize()

Now, we can mesh the geometry;

lc = 0.005

gmsh.option.setNumber("Mesh.MeshSizeMin", 0.005)
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.option.setNumber("Mesh.Algorithm3D", 10)#10
gmsh.option.setNumber("Mesh.Optimize", 1)
gmsh.option.setNumber("Mesh.OptimizeNetgen", 0)
gmsh.model.mesh.generate(3)

gmsh.option.setNumber("Mesh.SaveAll", 0)
gmsh.write("{}.msh".format("cube_mesh"))

if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

gmsh.finalize()

This will output this view;

Using the same procedure, you can assign the surface tags for your geometry, and you can manually track the tag you are interested in.

You can find the full script here.

I hope that helps.

Rene_Thierry_Djoumes · March 31, 2023, 12:05pm

Hi Ekrem,

Thank you very much for your answer.
Indeed, I had already tried this strategy. It works for a simple geometry.
But for a complex geometry, it considers all faces of each different terahedron and assigns them different tags. So for 20000 faces, we have 20000 tags.
look at the geometry in picture
Thanks

Ekrem_Ekici · March 31, 2023, 2:20pm

Where do you want to tag exactly in this picture?

Rene_Thierry_Djoumes · March 31, 2023, 4:10pm

on the interior surface, exterior surface and on the surfaces at the ends

Ekrem_Ekici · March 31, 2023, 4:21pm

You need to smoothen the surfaces but I don’t know how. Can you generate this geometry as a step file.

Also regarding the end surfaces you can track the surfaces according to their coordinates. For example,

inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
    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, [0, H/2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H/2, 0]):
            outflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L/2, H, 0]) or np.allclose(center_of_mass, [L/2, 0, 0]):
            walls.append(boundary[1])
        else:
            obstacle.append(boundary[1])
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

in this tutorial.

Rene_Thierry_Djoumes · March 31, 2023, 5:03pm

Thank you for your answer. this geometry is already in step file
In fact the original geometry was in .stl, I imported in FreeCAD to have the step file.

you are right we need to find a way to smooth the surface but I don’t know how to do it.

Thanks for your help

Kei_Yamamoto · April 1, 2023, 8:24pm

@Rene_Thierry_Djoumes

VaMPy has functionality for pre-processing vascular geomtry. You can use .stl as an input.

VaMPy uses VTMK for meshing and you can also directly use VMTK for surface smoothing, mesh generation with IDs for inlet/outlet.

http://www.vmtk.org/tutorials/SurfaceForMeshing.html

http://www.vmtk.org/tutorials/MeshGeneration.html

Rene_Thierry_Djoumes · April 3, 2023, 6:47am

Hi Kei,
Thank for your suggestion
I’ll tried and let you know.
Thank you

deshik · May 31, 2023, 7:35pm

Hope the following code helps. The code generates face groups based on angle between neighbouring cell faces.

Requires mesh file as input

import os
import sys
from dolfin import *
import numpy as np

def angle_between(v1, v2):
    v1 = v1/np.linalg.norm(v1)
    v2 = v2/np.linalg.norm(v2)
    angle = np.arccos(np.dot(v1, v2))
    if np.isnan(angle):
        if (v1 == v2).all():
            return 0.0
        else:
            return np.pi
    return angle

def _mark_loop(bmesh, normals, threshold_angles, domains, queue, visited, marker, threshold=30):
    queue_indices = np.where(queue.array())[0]

    num_marked = 0
    emap = bmesh.entity_map(2)
    while len(queue_indices) > 0:
        idx = queue_indices[0]
        facet = Cell(bmesh, idx)
        #FACET = Facet(mesh, emap[idx])
        visited[facet] = True
        queue[facet] = False
        domains[facet] = marker
        N = normals[idx]
        threshold_radians = np.pi*threshold_angles[facet]/180.

        neighbours = []
        for edge in edges(facet):
            for neighbour in cells(edge):
                if not visited[neighbour]:
                    neighbours.append(neighbour)

        for neighbour in neighbours:
            n = normals[neighbour.index()]

            angle = angle_between(N,n)
            if angle < threshold_radians:
                queue[neighbour] = not visited[neighbour]
        queue_indices = np.where(queue.array())[0]
        num_marked += 1

def mark_domains(mesh, threshold_angle=40):
    """Mark domains automatically by feature edges in mesh. Will connect
    facets to the same domain if angle between them is less than threshold_angle.
    Will also detect and mark disconnected parts of boundary.
    """
    assert MPI.size(mesh.mpi_comm()) == 1, "Currently only works on non-partitioned meshes."

    bmesh = BoundaryMesh(mesh, "exterior")
    visited_facets = MeshFunction("bool", bmesh, bmesh.topology().dim())
    visited_facets.set_all(False)
    bmesh.init(1,2)

    queue = MeshFunction("bool", bmesh, bmesh.topology().dim())
    queue.set_all(False)

    domains = MeshFunction("size_t", bmesh,  bmesh.topology().dim())
    domains.set_all(0)

    normals = np.zeros((bmesh.num_cells(), 3), dtype=np.float_)
    for facet in cells(bmesh):
        N = Facet(mesh, bmesh.entity_map(2)[facet.index()]).normal()
        normals[facet.index()][0] = N.x()
        normals[facet.index()][1] = N.y()
        normals[facet.index()][2] = N.z()


    threshold_angles = MeshFunction("size_t", bmesh, bmesh.topology().dim())
    threshold_angles.set_all(threshold_angle)

    unvisited = np.where(visited_facets.array() == False)[0]
    I = 0
    while len(unvisited) > 0:
        I += 1
        queue[unvisited[0]] = True
        _mark_loop(bmesh, normals, threshold_angles, domains, queue, visited_facets, I)
        unvisited = np.where(visited_facets.array() == False)[0]

    domains_full = MeshFunction("size_t", mesh, 2)
    domains_full.set_all(0)
    for facet in cells(bmesh):
        domains_full[bmesh.entity_map(2)[facet.index()]] = domains[facet]
    return domains_full

mesh = Mesh()
h5_file = "input-h5 file"
hdf = HDF5File(mesh.mpi_comm(), h5_file, "r")
hdf.read(mesh, "/mesh", False)
domains = mark_domains(mesh, 40)   # 40 is threshold angle that suits for my mesh
File("boundaries.pvd") << boundaries

credits: fenics forum

Rene_Thierry_Djoumes · January 29, 2024, 9:05am

Dear Deshik @deshik

I’m so sorry for the late reply.
In fact your code allowed me to write a similar one for Salome and another one for Gmsh which works very well as I can’t manage to write a good h5 file. but I managed to separate my stl (associated image). my problem is that FEniCS seems to really recognise my surface. but doesn’t assign the necessary values for a condition to the dirichlet type boundary.
I don’t understand exactly what the problem is.

deshik · February 2, 2024, 5:33pm

Can i see a minimal snippet of the way you are importing the mesh and applying the boundary conditions?
Just to remind, dolfin-convert works only with gmsh version-2 meshes.

Topic		Replies	Views
How to specify Boundary conditions when the mesh is imported from GMSH? mesh	16	4417	October 28, 2022
Create dirichlet boundary conditions from gmsh model mesh	2	533	June 3, 2023
How do I tag and use boundary surfaces from gmsh correctly? General	6	73	November 6, 2024
Trouble with marking a face in 3D STL files mesh	4	1695	March 30, 2023
Transfering complex multidomain geometry from 3DS Max to gmsh mesh	3	658	March 31, 2020

How to select complex surface for boundary conditions

Related topics