Generating two dimensional surface mesh pygmsh

46 Views Asked by At

I am new to the stack overflow. I hope you guys can help me.

I am trying to define a two-dimensional surface mesh from given boundary points using the pygmsh package. Ultimately I want to solve the Helmholtz equation (eigenvalue problem) on this mesh using the finite element method from skfem.

I got the program for the mesh generation to work for flat shapes e.g., a circle. However, I ran into problems with non-flat shapes e.g., a sphere. The program never ends. However, it also does not produce an error. Even for less than a dozen discretization steps the program runs for literal hours without any response or error. For some really small number of discretization steps, I do get an error message.

Why does nothing happen for non flat shapes? Is it realistic to have hours of computation time for such a problem? Did I miss something? Could it be my hardware?

FYI my set up is: Macbook Pro Apple M1 2020 Python 3.10.11 pygmsh 7.1.17

The working code is below. For the circle, it works fine. But the moment I introduce a non-trivial z direction it somehow breaks. I tried removing duplicates but that did not help either.

import pygmsh
import pyvista
import numpy as np
from math import pi

'''
### Circle
r=1
nn = 50
GridPoint = []

for ii in range(nn):
    phi=2*pi*ii/nn
    x = r * np.cos(phi)
    y = r * np.sin(phi)
    GridPoint.append([x,y,0])
'''

#### Sphere
r = 1;              # define radius
nn = 50;            # discretization steps
GridPoint = []      # empty list to save grid points

# Define sphere through spherical coordinates
for jj in range(nn):
    for ii in range(nn):
        x = r * np.sin(pi*ii/nn) * np.cos(2*pi*jj/nn)
        y = r * np.sin(pi*ii/nn) * np.sin(2*pi*jj/nn)
        z = r * np.cos(pi*ii/nn)
        GridPoint.append([x,y,z])

with pygmsh.geo.Geometry() as geom:
    geom.add_polygon(GridPoint, mesh_size = 0.2)
    mesh = geom.generate_mesh(dim=2)


p = pyvista.Plotter(window_size=(800, 800))
p.add_mesh(
    mesh=pyvista.from_meshio(mesh),
    show_edges=True,
)
p.show()

However, if set the dimension of the mesh to 1, i.e., mesh = geom.generate_mesh(dim=1) the program generates an output. However, works I need a mesh to solve the PDE. For really a small value of grid points nn<9 I sometimes get the error

>Traceback (most recent call last):
> File "../file.py", line 34, in <module> mesh = geom.generate_mesh(dim=2)
> File "lib/python3.10/pygmsh/common/geometry.py" , line 374, in generate_mesh gmsh.model.mesh.generate(dim)
> File "/lib/python3.10/site-packages/gmsh.py", line 2046, in generate raise Exception(logger.getLastError())
> Exception: Unable to recover the edge 75 (2/3) on curve 11 (on surface 1)

FYI my set up is: Macbook Pro Apple M1 2020 Python 3.10.11 pygmsh 7.1.17

0

There are 0 best solutions below