jorgensd / jorgensd.github.io

Personal web page
http://jsdokken.com/
0 stars 0 forks source link

trying to create a 3d mesh using meshio #15

Closed Mirialex closed 3 months ago

Mirialex commented 3 months ago

im trying to create a 3d model using @dokken 's doc but for some reason im getting the following error any idea why?

Jørgen S. Dokken Mesh generation and conversion with GMSH and PYGMSH 1 Personal webpage

`

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 build in kernel in GMSH

geometry = pygmsh.geo.Geometry()

Fetch model we would like to add data to

model = geometry.enter()

Add circle

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

-

The next step is to create the channel with the circle as a hole.

+

Add points with finer resolution on left side

points = [model.add_point((0, 0, 0), mesh_size=resolution), model.add_point((L, 0, 0), mesh_size=5resolution), model.add_point((L, H, 0), mesh_size=5resolution), 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])

Call gmsh kernel before add physical entities

model.synchronize()

-

The final step before mesh generation is to mark the different boundaries and the volume mesh. Note that with pygmsh, boundaries with the same tag has to be added simultaneously. In this example this means that we have to add the top and

bottom wall in one function call.

volume_marker = 6 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")

We generate the mesh using the pygmsh function generate_mesh. Generate mesh returns a meshio.Mesh. However, this mesh is tricky to extract physical tags from. Therefore we write the mesh to file using the gmsh.write function.

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

2. How to convert your mesh to XDMF

Now that we have save the mesh to a msh file, we would like to convert it to a format that interfaces with DOLFIN and DOLFINx.

For this I suggest using the XDMF-format as it supports parallel IO.

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

Now that we have loaded the mesh, we need to extract the cells and physical data. We need to create a separate file for the facets (lines), which we will use when we define boundary conditions in DOLFIN/DOLFINx. We do this with the following convenience function. Note that as we would like a 2 dimensional mesh, we need to remove the z-values in the mesh coordinates.

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

With this function at hand, we can save the meshes to XDMF.

+

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)

-

3. How to create a 3D mesh using pygmsh

To create more advanced meshes, such as 3D geometries, using the OpenCASCADE geometry kernel is recommended.

We start by importing this kernel, and creating three objects:

- A box $[0,0,0]\times[1,1,1]$

- A box $[0.5,0.0.5,1]\times[1,1,2]$

- A ball from $[0.5,0.5,0.5]$ with radius $0.25$.

Clear previous model

mesh_size = 0.1 geom = pygmsh.occ.Geometry() model3D = geom.enter() box0 = model3D.add_box([0.0, 0, 0], [1, 1, 1]) box1 = model3D.add_box([0.5, 0.5, 1], [0.5, 0.5, 1]) ball = model3D.add_ball([0.5, 0.5, 0.5], 0.25)

In this demo, we would like to make a mesh that is the union of these three objects.

In addition, we would like the internal boundary of the sphere to be preserved in the final mesh.

We will do this by using boolean operations. First we make a boolean_union of the two boxes (whose internal boundaries will not be preserved). Then, we use boolean fragments to perserve the outer boundary of the sphere.

union = model3D.boolean_union([box0, box1]) union_minus_ball = model3D.boolean_fragments(union, ball) model3D.synchronize()

To create physical markers for the two regions, we use the add_physical function. This function only works nicely if the domain whose boundary should be preserved (the sphere) is fully embedded in the other domain (the union of boxes). For more complex operations, it is recommened to do the tagging of entities in the gmsh GUI, as explained in the GMSH tutorial.

model3D.add_physical(union, "Union") model3D.add_physical(union_minus_ball, "Union minus ball")

We finally generate the 3D mesh, and save both the geo and msh file as in the previous example.

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

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

These XDMF-files can be visualized in Paraview and looks like

#

The 2D mesh and the corresponding facet data visualized in Paraview

#

We use the same strategy for the 3D mesh as the 2D mesh.

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 triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) meshio.write("facet_mesh3D.xdmf", triangle_mesh)

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

` and getting the following error any idea why? Traceback (most recent call last): File "/home/mirialex/PycharmProjects/fibers/from_community.py", line 133, in triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) File "/home/mirialex/PycharmProjects/fibers/from_community.py", line 128, in create_mesh cell_data = mesh.get_cell_data("gmsh:physical", cell_type) File "/usr/lib/python3/dist-packages/meshio/_mesh.py", line 249, in get_cell_data return np.concatenate( File "<__array_function__ internals>", line 5, in concatenate ValueError: need at least one array to concatenate

jorgensd commented 3 months ago

As you can see by the location of the error message, you have successfully converted the 3D mesh.

however, as you have not tagged any physical surfaces of the 3D mesh (ie. 2D triangles), you cannot read them in, as there aren’t any present in your msh file. Please use model3D.add_physical on relevant surface entities that you would like to export

Mirialex commented 3 months ago

hi, do you have a little doc on how to use those add_physical? if you have also some doc on how to use any other command with meshio i would appreciate it as well

‫בתאריך יום ג׳, 28 במאי 2024 ב-19:34 מאת ‪Jørgen Schartum Dokken‬‏ <‪ @.***‬‏>:‬

As you can see by the location of the error message, you have successfully converted the 3D mesh.

however, as you have not tagged any physical surfaces of the 3D mesh (ie. 2D triangles), you cannot read them in, as there aren’t any present in your msh file. Please use model3D.add_physical on relevant surface entities that you would like to export

— Reply to this email directly, view it on GitHub https://github.com/jorgensd/jorgensd.github.io/issues/15#issuecomment-2135681241, or unsubscribe https://github.com/notifications/unsubscribe-auth/AT5PBW46JQKDPZUBMHQAM2TZESWZHAVCNFSM6AAAAABINJAGPSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMZVGY4DCMRUGE . You are receiving this because you authored the thread.Message ID: @.***>

jorgensd commented 3 months ago

This is a pygmsh question. I am not a maintainer or developer of pygmsh. You should consider: https://pygmsh.readthedocs.io/en/latest/

These days i use the Gmsh Python API directly, which there are many examples of in my various tutorials.

Similarly, you would have to read the meshio API https://github.com/nschloe/meshio/tree/v5.3.5