2D structured mesh for hollow cylinder

1

Hi,

I am trying to generate a structured mesh for a 2-D complete hollow cylinder geometry. It is possible to do it by using mesh mapping. My problem is exactly the same with the one in this closed old topic:

https://fenicsproject.org/qa/9235/merging-nodes/

The boundary lines that closes the loop must be merged but I understand that it was not possible. Not sure if it is still not possible.

Eventhough it has an answer, I don’t actually understand what is meant by "iterating over the rectangle mesh and creating new vertices and cells by using MeshEditor. An example for that or an alternative solution to this specific problem would be highly appreciated.

Note: I cannot use another software for generating mesh for this specific problem at the moment.

2

What Garth means as you are trying to remove vertices, and change the connectivity of each cell, you need to create a new mesh. This can be done with MeshEditor.
How to use MeshEditor has been covered in for instance:

1 Like
3

Hi dokken,

Thank you for your help. I checked the link and repeated the codes you provided. But now I have to ask for more. I understand that since I have double vertices on that line, you suggest me to remove those vertices and re-add vertices exactly to the same spots so that my new mesh would become a continous one.

I change my mesh iteratively in my actual code. So, this procedure is going to be in a loop in my case. But I hope I can handle that once I can do it for a single case.

Here is the working code for the mesh:

from fenics import *
import numpy as np
%matplotlib inline

Theta = 2*pi

nr = 10
nt = 40

a_x = 1
a_y = 0
b_x = 2
b_y = 1

a = Point(a_x,a_y)
b = Point(b_x,b_y)
mesh = RectangleMesh(a, b, nr, nt)

x = mesh.coordinates()[:,0]
y = mesh.coordinates()[:,1]

s = 1

def denser(x, y):
    return [a_x + (b_x-a_x)*((x-a_x)/(b_x-a_x))**s, y]

x_bar, y_bar = denser(x, y)
xy_bar_coor = np.array([x_bar, y_bar]).transpose()
mesh.coordinates()[:] = xy_bar_coor

def cylinder(r, s):
    return [r*np.cos(Theta*s), r*np.sin(Theta*s)]

x_hat, y_hat = cylinder(x_bar, y_bar)
xy_hat_coor = np.array([x_hat, y_hat]).transpose()
mesh.coordinates()[:] = xy_hat_coor

plot(mesh)

Here they are constants but actually nt and nr are variables in my actual code. That is why I need to know those vertices everytime I construct a new mesh. Here are my questions:

1- How can I locate vertices of a mesh? I hope there is a single easy function for it (also, if there is an easy way: how can I locate vertices on some specific interface/line etc.)

2- How can I remove vertices?

I think I can then use your method to add necessary vertices that I already would’ve located and removed (doubles).

4

Following code does what you want:

from IPython import embed
import matplotlib.pyplot as plt
from fenics import *
import numpy as np

Theta = 2*pi

nr = 10
nt = 10

a_x = 1
a_y = 0
b_x = 2
b_y = 1

a = Point(a_x, a_y)
b = Point(b_x, b_y)
mesh = RectangleMesh(a, b, nr, nt)
mesh.init(mesh.topology().dim(), 0)
c_to_vertex = mesh.topology()(mesh.topology().dim(), 0)
x = mesh.coordinates()[:, 0]
y = mesh.coordinates()[:, 1]


mf = MeshFunction("size_t", mesh, 0, 0)


class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], a[1])


class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], b[1])


# Create map between top and bottom on the square before deforming it
top().mark(mf, 1)
bottom().mark(mf, 2)
top_vertices = np.flatnonzero(mf.array() == 1)
bottom_vertices = np.flatnonzero(mf.array() == 2)
top_to_bottom = {}
for top in top_vertices:
    bottom = np.flatnonzero(x[bottom_vertices] == x[top])[0]
    top_to_bottom[top] = bottom_vertices[bottom]


def denser(x, y):
    return [a_x + (b_x-a_x)*((x-a_x)/(b_x-a_x))**s, y]


def cylinder(r, s):
    return [r*np.cos(Theta*s), r*np.sin(Theta*s)]


# Convert the mesh to hollow circle
s = 2.
x_bar, y_bar = denser(x, y)
xy_bar_coor = np.array([x_bar, y_bar]).transpose()
mesh.coordinates()[:] = xy_bar_coor
x_hat, y_hat = cylinder(x_bar, y_bar)
xy_hat_coor = np.array([x_hat, y_hat]).transpose()
mesh.coordinates()[:] = xy_hat_coor


# Create mesh object and open editor
new_mesh = Mesh()
editor = MeshEditor()
topological_dim = 2
geometrical_dim = 2
num_local_vertices = mesh.num_entities(0) - len(top_vertices)
num_global_vertices = num_local_vertices  # True if run in serial
num_local_cells = mesh.num_cells()
num_global_cells = num_local_cells
editor.open(new_mesh, "triangle", topological_dim, geometrical_dim)
editor.init_vertices_global(num_local_vertices, num_global_vertices)
editor.init_cells_global(num_local_cells, num_global_cells)

# Add new vertices (renumber as there are fewer vertices in new mesh)
map_vertex = {}
v_count = 0
for i, coordinate in enumerate(mesh.coordinates()):
    if not (i in top_vertices):
        editor.add_vertex(v_count, coordinate)
        map_vertex[i] = v_count
        v_count += 1
    else:
        map_vertex[i] = top_to_bottom[i]

# Add the cells
for i in range(num_global_cells):
    vertices = c_to_vertex(i)
    new_vertices = []
    for vertex in vertices:
        if vertex in top_vertices:
            new_vertices.append(map_vertex[map_vertex[vertex]])
        else:
            new_vertices.append(map_vertex[vertex])
    editor.add_cell(i, np.asarray(new_vertices, dtype='uint'))

# Close editor
editor.close()

print(new_mesh.num_entities(0), mesh.num_entities(0))
print(new_mesh.num_cells(), mesh.num_cells())

plot(new_mesh, color="r")
plt.savefig("hollow.png")
2 Likes
5

This is amazing! Thank you very much. But there is still a problem with the interface. Looks like that boundary is seen as a wall when I run my code. I guess I need to specify that new interlayer is a continous part of the domain.

Do you know why that is happening?

6

As you have not provided a minimal code example that illustrate what issues you are having, i cannot really help you.

1 Like
7

I will try to figure out an example code and a solution if I can. That won’t be easy for me right away. But thank you so much, you have done an amazing favor! I am grateful.

I know the problem is still with that interface because I have checked my solution with standard Circle(), generate_mesh() etc. in Fenics, and got the correct results. I will work on it and hopefuly write later.

8

Make sure you use the new_mesh, and not mesh the old mesh. You could start by creating a very coarse mesh and tabulating the dof coordinates, and visually inspect them.

1 Like
9

This is very funny but it is better that I made that mistake cause others who read this post may make the same. Since you updated the mesh I thought everything was fine and kept solving my problem with the “mesh”, just plugging-in the code. This was super funny and thank you again for guessing my mistake :smiley:

I changed the mesh to new_mesh you wrote in my code and the solution is perfect now. Thank you! :smiley: