How to select complex surface for boundary conditions

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.

2 Likes

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;

image

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.

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

Where do you want to tag exactly in this picture?

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

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.

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

@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

1 Like

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

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

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.

image1.png
image

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.