Hello ,
I want to save the data on the meshgrid for the complex geometry (for ML modelling), which is imported from gmsh. But I am getting discrepancies in compute_midpoints between the in-build mesh and mesh generated using Gmsh.
Here is the simple MWE I experimented as an comparison on simple square of 5*5 quad elements.
import numpy as np
import meshio
from mpi4py import MPI
from dolfinx import mesh, io
from dolfinx.mesh import CellType, GhostMode, create_unit_square, compute_midpoints
from dolfinx.io import XDMFFile
from dolfinx.fem import functionspace
import matplotlib.pyplot as plt
comm = MPI.COMM_WORLD
rank = comm.rank
def create_gmsh_mesh():
# Load and process the GMSH mesh
domain = meshio.read("gmshsquare.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.astype(np.int32)]})
return out_mesh
if rank == 0:
triangle_mesh = create_mesh(domain, "quad", prune_z=True)
meshio.write("mesh_gmshsquare.xdmf", triangle_mesh)
MPI.COMM_WORLD.barrier()
with XDMFFile(MPI.COMM_WORLD, "mesh_gmshsquare.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
return domain
def create_dolfinx_mesh():
domain = create_unit_square(MPI.COMM_WORLD, 5, 5,
CellType.quadrilateral, ghost_mode=GhostMode.shared_facet)
with XDMFFile(MPI.COMM_WORLD, "mesh_dolfinxsquare.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
return domain
def process_mesh(domain, mesh_type):
# Example of processing the mesh to get number of elements and nodes
num_elements = domain.topology.index_map(2).size_local
num_nodes = domain.topology.index_map(0).size_local
num_elements_global = MPI.COMM_WORLD.allreduce(num_elements, op=MPI.SUM)
num_nodes_global = MPI.COMM_WORLD.allreduce(num_nodes, op=MPI.SUM)
print(f"Number of local elements: {num_elements}")
print(f"Number of local nodes: {num_nodes}")
print(f"Number of elements: {num_elements_global}")
print(f"Number of nodes: {num_nodes_global}")
# Compute midpoints
Vh0 = functionspace(domain, ("DG", 0))
# midpoints = dolfinx.mesh.compute_midpoints(domain, tdim, np.arange(num_elements, dtype=np.int32))
midpoints = Vh0.tabulate_dof_coordinates()
xcoords = midpoints[:, 0]
ycoords = midpoints[:, 1]
E = np.arange(1, len(xcoords) + 1)
# Ensure E fits the dimensions of the grid
unique_x = np.unique(xcoords)
unique_y = np.unique(ycoords)
grid_x, grid_y = np.meshgrid(unique_x, unique_y)
E_reshaped = np.resize(E, grid_x.shape)
return xcoords, ycoords, grid_x, grid_y, E_reshaped
def plot_comparison(grid_x_gmsh, grid_y_gmsh, E_gmsh, grid_x_dolfinx, grid_y_dolfinx, E_dolfinx):
fig, ax = plt.subplots(1, 2, figsize=(20, 8))
# Plot GMSH Mesh
scatter_gmsh = ax[0].scatter(grid_x_gmsh, grid_y_gmsh, c=E_gmsh.flatten(), cmap='viridis', marker='o')
fig.colorbar(scatter_gmsh, ax=ax[0], label='E Value')
ax[0].set_title('GMSH Mesh Grid with E')
ax[0].set_xlabel('X Coordinate')
ax[0].set_ylabel('Y Coordinate')
ax[0].grid(True)
# Plot Dolfinx Mesh
scatter_dolfinx = ax[1].scatter(grid_x_dolfinx, grid_y_dolfinx, c=E_dolfinx.flatten(), cmap='viridis', marker='o')
fig.colorbar(scatter_dolfinx, ax=ax[1], label='E Value')
ax[1].set_title('Dolfinx Mesh Grid with E')
ax[1].set_xlabel('X Coordinate')
ax[1].set_ylabel('Y Coordinate')
ax[1].grid(True)
# Save the figure
fig.savefig("comparisonplot.png", bbox_inches='tight') # Save the figure with tight bounding box
# Adjust spacing between subplots
fig.subplots_adjust(wspace=1) # Adjust wspace as needed
plt.show()
def print_comparison(xcoords_gmsh, ycoords_gmsh, E_gmsh, xcoords_dolfinx, ycoords_dolfinx, E_dolfinx):
print("GMSH Mesh Coordinates and E Values:")
print("X Coordinate | Y Coordinate | E Value")
print("-------------------------------------")
for x, y, e in zip(xcoords_gmsh, ycoords_gmsh, E_gmsh.flatten()):
print(f"{x:12.1f} | {y:12.1f} | {e:6}")
print("\nDolfinx Mesh Coordinates and E Values:")
print("X Coordinate | Y Coordinate | E Value")
print("-------------------------------------")
for x, y, e in zip(xcoords_dolfinx, ycoords_dolfinx, E_dolfinx.flatten()):
print(f"{x:12.1f} | {y:12.1f} | {e:6}")
if __name__ == '__main__':
# Create and process both meshes
domain_gmsh = create_gmsh_mesh()
xcoords_gmsh, ycoords_gmsh, grid_x_gmsh, grid_y_gmsh, E_gmsh = process_mesh(domain_gmsh, 'gmsh')
domain_dolfinx = create_dolfinx_mesh()
xcoords_dolfinx, ycoords_dolfinx, grid_x_dolfinx, grid_y_dolfinx, E_dolfinx = process_mesh(domain_dolfinx,
'dolfinx')
# Plot side by side
plot_comparison(grid_x_gmsh, grid_y_gmsh, E_gmsh, grid_x_dolfinx, grid_y_dolfinx, E_dolfinx)
# Print comparison
print_comparison(xcoords_gmsh, ycoords_gmsh, E_gmsh, xcoords_dolfinx, ycoords_dolfinx, E_dolfinx)
This is plot of comparison
and the corresponding values on the midpoints
GMSH Mesh Coordinates and E Values:
X Coordinate | Y Coordinate | E Value
-------------------------------------
0.1 | 0.9 | 1
0.3 | 0.9 | 2
0.1 | 0.7 | 3
0.5 | 0.9 | 4
0.3 | 0.7 | 5
0.1 | 0.5 | 6
0.7 | 0.9 | 7
0.5 | 0.7 | 8
0.3 | 0.5 | 9
0.1 | 0.3 | 10
0.9 | 0.9 | 11
0.7 | 0.7 | 12
0.5 | 0.5 | 13
0.3 | 0.3 | 14
0.1 | 0.1 | 15
0.9 | 0.7 | 16
0.7 | 0.5 | 17
0.5 | 0.3 | 18
0.3 | 0.1 | 19
0.9 | 0.5 | 20
0.7 | 0.3 | 21
0.5 | 0.1 | 22
0.9 | 0.3 | 23
0.7 | 0.1 | 24
0.9 | 0.1 | 25
Dolfinx Mesh Coordinates and E Values:
X Coordinate | Y Coordinate | E Value
-------------------------------------
0.1 | 0.1 | 1
0.1 | 0.3 | 2
0.3 | 0.1 | 3
0.1 | 0.5 | 4
0.3 | 0.3 | 5
0.5 | 0.1 | 6
0.1 | 0.7 | 7
0.3 | 0.5 | 8
0.5 | 0.3 | 9
0.7 | 0.1 | 10
0.1 | 0.9 | 11
0.3 | 0.7 | 12
0.5 | 0.5 | 13
0.7 | 0.3 | 14
0.9 | 0.1 | 15
0.3 | 0.9 | 16
0.5 | 0.7 | 17
0.7 | 0.5 | 18
0.9 | 0.3 | 19
0.5 | 0.9 | 20
0.7 | 0.7 | 21
0.9 | 0.5 | 22
0.7 | 0.9 | 23
0.9 | 0.7 | 24
0.9 | 0.9 | 25
Attaching the .msh file
The inbuilt mesh is performing as expected, but not the gmsh mesh. I really appreciate if anyone can help me to lift this is issue sooner.
Thanks!