nschloe / meshio

:spider_web: input/output for many mesh formats
MIT License
1.93k stars 397 forks source link

creating a 3d mesh with meshio #1471

Open Mirialex opened 4 months ago

Mirialex commented 4 months ago

hi all,

need help in creating 3d mesh with meshio, basically my goal is to create a box with a cylinder inside of it, starting with creating a box this is what i tried to do ` import numpy import meshio import gmsh import pygmsh

resolution = 0.01

Channel parameters

L = 2.2 H = 0.41 c = [0.2, 0.2, 0] r = 0.05

Initialize empty geometry using the built-in kernel in GMSH

geometry = pygmsh.geo.Geometry() model = geometry.enter()

Add circle

circle = model.add_circle(c, r, mesh_size=resolution)

Add points with finer resolution on the left side

points = [ model.add_point((0, 0, 0), mesh_size=resolution), model.add_point((L, 0, 0), mesh_size=5 resolution), model.add_point((L, H, 0), mesh_size=5 resolution), model.add_point((0, H, 0), mesh_size=resolution) ]

Add lines between all points creating the rectangle

channel_lines = [model.add_line(points[i], points[i + 1]) for i in range(-1, len(points) - 1)]

Create a line loop and plane surface for meshing

channel_loop = model.add_curve_loop(channel_lines) plane_surface = model.add_plane_surface(channel_loop, holes=[circle.curve_loop])

Synchronize the model

model.synchronize()

Mark different boundaries and the volume mesh

model.add_physical([plane_surface], "Volume") model.add_physical([channel_lines[0]], "Inflow") model.add_physical([channel_lines[2]], "Outflow") model.add_physical([channel_lines[1], channel_lines[3]], "Walls") model.add_physical(circle.curve_loop.curves, "Obstacle")

Generate the mesh and write to file

geometry.generate_mesh(dim=2) gmsh.write("mesh.msh") gmsh.clear() geometry.exit()

Convert the mesh to XDMF format

mesh_from_file = meshio.read("mesh.msh")

def create_mesh(mesh, cell_type, prune_z=False): cells = mesh.get_cells_type(cell_type) cell_data = mesh.get_cell_data("gmsh:physical", cell_type) points = mesh.points[:, :2] if prune_z else mesh.points out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]}) return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True) meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) meshio.write("mesh.xdmf", triangle_mesh)

3D Mesh generation using OpenCASCADE geometry kernel

L = 1.0 H = 1.0 W = 1.0 resolution = 0.1

geom = pygmsh.occ.Geometry() model3D = geom.enter()

points = [ model3D.add_point([0.0, 0.0, 0.0], mesh_size=resolution), model3D.add_point([L, 0.0, 0.0], mesh_size=5 resolution), model3D.add_point([L, H, 0.0], mesh_size=5 resolution), model3D.add_point([0.0, H, 0.0], mesh_size=resolution), model3D.add_point([0.0, 0.0, W], mesh_size=resolution), model3D.add_point([L, 0.0, W], mesh_size=5 resolution), model3D.add_point([L, H, W], mesh_size=5 resolution), model3D.add_point([0.0, H, W], mesh_size=resolution) ]

lines = [model3D.add_line(points[i], points[i + 1]) for i in range(-1, len(points) - 1) ]

Define surfaces (make sure lines are in the correct order to form closed loops)

surfaces = [] surfaces.append(model3D.add_plane_surface(model3D.add_curve_loop([lines[0], lines[1], lines[2], lines[3]]))) # Bottom surfaces.append(model3D.add_plane_surface(model3D.add_curve_loop([lines[4], lines[5], lines[6], lines[7]]))) # Top surfaces.append(model3D.add_plane_surface(model3D.add_curve_loop([lines[0], lines[9], lines[4], lines[8]]))) # Front surfaces.append(model3D.add_plane_surface(model3D.add_curve_loop([lines[1], lines[10], lines[5], lines[9]]))) # Right surfaces.append(model3D.add_plane_surface(model3D.add_curve_loop([lines[2], lines[11], lines[6], lines[10]]))) # Back surfaces.append(model3D.add_plane_surface(model3D.add_curve_loop([lines[3], lines[8], lines[7], lines[11]]))) # Left

Define the volume

surface_loop = model3D.add_surface_loop(surfaces) volume = model3D.add_volume(surface_loop)

Add Physical Groups

model3D.add_physical(volume, label="MyVolume") for i, surface in enumerate(surfaces): model3D.addphysical(surface, label=f"MySurface{i}")

geom.generate_mesh(dim=3) gmsh.write("mesh3D.msh") model3D.exit()

mesh3D_from_file = meshio.read("mesh3D.msh")

triangle_mesh = create_mesh(mesh3D_from_file, "triangle", prune_z=True) meshio.write("facet_mesh3D.xdmf", triangle_mesh)

tetra_mesh = create_mesh(mesh3D_from_file, "tetra", prune_z=True) meshio.write("mesh3D.xdmf", tetra_mesh)

`

but im not sure how to use this add_curve_loop and how to change it in order for it to work if you have any documention about it, would appreciate

error i got /usr/bin/python3.10 /home/mirialex/PycharmProjects/fibers/from_community.py Traceback (most recent call last): File "/home/mirialex/PycharmProjects/fibers/from_community.py", line 94, in surfaces.append(model3D.add_plane_surface(model3D.add_curve_loop([lines[0], lines[1], lines[2], lines[3]]))) # Bottom File "/usr/lib/python3/dist-packages/pygmsh/common/geometry.py", line 80, in add_curve_loop return CurveLoop(self.env, *args, **kwargs) File "/usr/lib/python3/dist-packages/pygmsh/common/curve_loop.py", line 26, in init assert curves[-1].points[-1] == curves[0].points[0] AssertionError

Process finished with exit code 1