nschloe / pygmsh

:spider_web: Gmsh for Python
GNU General Public License v3.0
858 stars 162 forks source link

Geom.revolve throws Index of Range Exception for Surface Mesh #588

Open asmithAE22 opened 2 weeks ago

asmithAE22 commented 2 weeks ago

Attempting to revolve a 2D cross section into a surface of revolution with geom.revovle() fails with error in CommonGeometry._revolve :

    top, surf, _ = geom.revolve(
                   ^^^^^^^^^^^^^
  File "<redacted>\anaconda3\envs\<redacted>\Lib\site-packages\pygmsh\geo\geometry.py", line 48, in revolve
    return super()._revolve(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<redacted>\anaconda3\envs\<redacted>\Lib\site-packages\pygmsh\common\geometry.py", line 218, in _revolve
    extruded = Dummy(*out_dim_tags[1])
                      ~~~~~~~~~~~~^^^
IndexError: list index out of range

It seems the extruded and lateral values arent evaluated. The code to replicate is shown below:

R = data[:,0]
Z = data[:,1]
components = data[:,2]

regen = True
if regen:
    with pygmsh.geo.Geometry() as geom:
        axis = [0,0,1]

        # Define the RZ profile points
        points = []
        mesh_size = .001
        for i,(r, z) in enumerate(zip(R, Z)):
            if i != 0:
                mesh_size = float(np.linalg.norm(data[i,0:2]-(data[i-1,0:2])))
            if i < 10:
                mesh_size = np.min([mesh_size, .005])
            else:
                mesh_size = np.min([mesh_size, .01])
            mesh_size = 2* np.min([mesh_size, .01])
            points.append(geom.add_point([r, 0, z],mesh_size=mesh_size))
        points.append(geom.add_point([0, 0, Z[-1]], mesh_size=.001))
        points.append(points[0])

        # Create a line through the points to form a curve
        lines = []
        for i in range(len(points) - 1):
            lines.append(geom.add_line(points[i], points[i + 1]))

        # line_ids = np.array(lines[:-1])
        # for comp in np.unique(components):
        #     mask = components == comp
        #     geom.add_physical(list(line_ids[mask]), str(comp))

        curve_loop = geom.add_curve_loop(lines)
        for k, p in enumerate(lines):
            n = 4
            for _ in range(n):
                top, surf, _ = geom.revolve(
                    p,
                    rotation_axis=[0,0,1],
                    point_on_axis=[0,0,0],
                    angle=2*np.pi/n,
                )
                # print(surf, dir(surf))
                p = top

        mesh = geom.generate_mesh(2)

        # Save the mesh to a file
        mesh.write(str(Path.cwd() / "surface_mesh.vtk"))