Create a mesh with boundary markers in Python, Fenics

Question

My aim is to create a mesh in Python and setting markers on some 1-dimensional subset.

Up until now, I have always creates a set, for example a rectangle in gmsh, then put for example a circle in it. Then gmsh puts a mesh on my structure and I can mark the boundary of the rectangle and of the circle as facets (as xdmf files). I can then let Python read my mesh and boundary facets and using it, f.e. to solve PDEs. What I want to do now is the following: I still want to have my rectangle with a mesh, but instead of defining the facet markes in gmsh, I want to define them using the image of a function.

More precicely: Instead of creating a circle in gmsh, I want to consider, for example, the function

Then I want to use as my facet and mark it in my mesh.

Is there a way to do this? I feel kind of lost here.

Answer 1

I am not sure if I understood everything in your question correctly, but here is some code that might help you. One thing I do not know is how to add element edges to physical groups in gmsh. But maybe you can figure that out.

So here is the code:

import gmsh
import numpy as np
    
def mapping(point):
    x = point[0]
    y = point[1]
    z = point[2]
    result = [2*x,3*y,z]
    return result

def inverseMapping(point):
    x = point[0]
    y = point[1]
    z = point[2]
    result = [(1/2)*x,(1/3)*y,z]
    return result

def getNodes():
    nodeTags, nodeCoord, _ = gmsh.model.mesh.getNodes()
    nodeCoord = np.reshape(nodeCoord,(len(nodeTags),3))
    return nodeTags, nodeCoord

def getEdgeNodeCoordinates():
    edgeTags, edgeNodes = gmsh.model.mesh.getAllEdges()
    edgeNodes = np.reshape(edgeNodes,(len(edgeTags),2))
    
    nodeTags, nodeCoord = getNodes()
    coord = []
    for i in range(0,len(edgeTags)):
        tagNode1 = edgeNodes[i][0]
        tagNode2 = edgeNodes[i][1]
        
        nodeIndex1 = list(nodeTags).index(tagNode1)
        nodeIndex2 = list(nodeTags).index(tagNode2)
        
        nodeCoord1 = nodeCoord[nodeIndex1]
        nodeCoord2 = nodeCoord[nodeIndex2]
        
        coord.append([nodeCoord1,nodeCoord2])
        
    return edgeTags, edgeNodes, nodeTags, coord
             
def getInverseNodeCoordinates(edgeNodeCoordinates):
    
    coord = []
    for edgeNodes in edgeNodeCoordinates:
        nodeCoord1 = edgeNodes[0]
        nodeCoord2 = edgeNodes[1]
        
        newCoord1 = inverseMapping(nodeCoord1)
        newCoord2 = inverseMapping(nodeCoord2)
        coord.append([newCoord1, newCoord2])
    return coord
        
def checkIntersection(edgeTags, edgeNodeCoordinates, inverseCoordinates):
    
    intersectingEdgeTags = []
    intersectingEdgeNodeCoord = []
    # 1 = inside, 0 = outside
    for i in range(0,len(inverseCoordinates)):
        pair = inverseCoordinates[i]
        coord1 = pair[0]
        coord2 = pair[1]
        
        e1 = 1 if np.linalg.norm(coord1) <= 1 else 0
        e2 = 1 if np.linalg.norm(coord2) <= 1 else 0
        
        s = e1 + e2 # s = 0 --> both nodes outside of manifold
                    # s = 1 --> one node inside and one node outside of manifold
                    # s = 2 --> both nodes inside of manifold
        
        if s == 1:
            intersectingEdgeTags.append(edgeTags[i])
            intersectingEdgeNodeCoord.append(edgeNodeCoordinates[i])
    return intersectingEdgeTags, intersectingEdgeNodeCoord    
        
def visualizeEdges(intersectingEdgeNodeCoord):
    
    for pair in intersectingEdgeNodeCoord:
        p1 = pair[0]
        p2 = pair[1]
        
        t1 = gmsh.model.occ.addPoint(p1[0],p1[1],p1[2])
        t2 = gmsh.model.occ.addPoint(p2[0],p2[1],p2[2])
        
        line = gmsh.model.occ.addLine(t1, t2)
        gmsh.model.occ.synchronize()
        gmsh.model.setColor([(1,line)], 255, 0, 0)
    gmsh.model.occ.synchronize()



gmsh.initialize()

# Create a rectangle which will be meshed later.
tag_vol_1 = gmsh.model.occ.addRectangle(-3, -4, 0, 6, 8)

# Sample the S1 manifold with n_points
S1_sampling_points = []
n_points = 100
maxAngle = 2*np.pi
angle = maxAngle/n_points
z = 0
for i in range(0,n_points):
    x = np.cos(i*angle)
    y = np.sin(i*angle)
    S1_sampling_points.append([x,y,z])

# Map the sampled S1 points to 2*x, 3*y, z.
# This is only for "visualization" via a spline.
mappedPoints = []
mappedPointTags = []
for point in S1_sampling_points:
    mappedPoint = mapping(point)
    tag = gmsh.model.occ.addPoint(mappedPoint[0], mappedPoint[1], mappedPoint[2])
    mappedPointTags.append(tag)

# Here the spline fitting is performed
# You will see it visualized when gmsh opens.
tagMappedS1 = gmsh.model.occ.addSpline(mappedPointTags + [mappedPointTags[0]]) # make the spline periodic by setting the last point equal to the first one
gmsh.model.occ.synchronize()

# Mesh the rectangle and tell gmsh to create edges which we can access.
gmsh.model.mesh.generate(2)
gmsh.model.mesh.createEdges() # You need to call this before using gmsh.model.mesh.getAllEdges()

# Get all these self-explanatory things
edgeTags, edgeNodes, nodeTags, edgeNodeCoordinates = getEdgeNodeCoordinates()

# Calculate the inverse-mapped coordinates of all nodes.
# With this we can just check if the transformed nodes are inside a circle of radius 1 or outside.
# If for every egde one node is inside, and the other node is outside the circle, then the edge is
# intersected by the mapped manifold f(S1) --> (2*x, 3*y, z). We then save the tag of such an edge.
inverseCoordinates = getInverseNodeCoordinates(edgeNodeCoordinates)
intersectingEdgeTags, intersectingEdgeNodeCoord = checkIntersection(edgeTags, edgeNodeCoordinates, inverseCoordinates)

# Here all intersecting edges are created within gmsh so you can see it.
# This is for visualization only. A lot of nodes with the same coordinates are created here twice.
visualizeEdges(intersectingEdgeNodeCoord)



gmsh.fltk.run()
gmsh.finalize()

And the result in gmsh looks like this: