Copyright (C) 2020 Garth N. Wells and Jørgen S. Dokken
This file is part of DOLFINX (https://www.fenicsproject.org)
SPDX-License-Identifier: LGPL-3.0-or-later
Mesh generation using Gmsh and pygmsh¶
import numpy as np
import pygmsh
from mpi4py import MPI
from dolfinx import cpp
from dolfinx.cpp.io import perm_gmsh, extract_local_entities
from dolfinx.io import XDMFFile, ufl_mesh_from_gmsh
from dolfinx.mesh import create_mesh, create_meshtags
Generating a mesh on each process rank¶
Generate a mesh on each rank with pygmsh, and create a DOLFIN-X mesh on each rank
geom = pygmsh.opencascade.Geometry()
geom.add_ball([0.0, 0.0, 0.0], 1.0, char_length=0.2)
pygmsh_mesh = pygmsh.generate_mesh(geom)
cells, x = pygmsh_mesh.cells[-1].data, pygmsh_mesh.points
pygmsh_cell = pygmsh_mesh.cells[-1].type
mesh = create_mesh(MPI.COMM_SELF, cells, x,
ufl_mesh_from_gmsh(pygmsh_cell, x.shape[1]))
with XDMFFile(MPI.COMM_SELF, "mesh_rank_{}.xdmf".format(MPI.COMM_WORLD.rank), "w") as file:
file.write_mesh(mesh)
Create a distributed (parallel) mesh with affine geometry¶
Generate mesh on rank 0, then build a distributed mesh
if MPI.COMM_WORLD.rank == 0:
# Generate a mesh
geom = pygmsh.opencascade.Geometry()
ball = geom.add_ball([0.0, 0.0, 0.0], 1.0, char_length=0.2)
box = geom.add_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0])
cut = geom.boolean_difference([ball], [box])
geom.add_raw_code("Physical Surface(1) = {1};")
geom.add_physical(cut, 2)
pygmsh_mesh = pygmsh.generate_mesh(geom)
# Extract the topology and geometry data
cells, x = pygmsh_mesh.cells[-1].data, pygmsh_mesh.points
# Extract marked facets
marked_entities = pygmsh_mesh.cells[-2].data
values = pygmsh_mesh.cell_data["gmsh:physical"][-2]
# Broadcast cell type data and geometric dimension
num_nodes = MPI.COMM_WORLD.bcast(cells.shape[1], root=0)
else:
num_nodes = MPI.COMM_WORLD.bcast(None, root=0)
cells, x = np.empty([0, num_nodes]), np.empty([0, 3])
marked_entities, values = np.empty((0, 3)), np.empty((0,))
mesh = create_mesh(MPI.COMM_WORLD, cells, x, ufl_mesh_from_gmsh("tetra", 3))
mesh.name = "ball_d1"
local_entities, local_values = extract_local_entities(mesh, 2, marked_entities, values)
mesh.topology.create_connectivity(2, 0)
mt = create_meshtags(mesh, 2, cpp.graph.AdjacencyList_int32(local_entities), np.int32(local_values))
mt.name = "ball_d1_surface"
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as file:
file.write_mesh(mesh)
mesh.topology.create_connectivity(2, 3)
file.write_meshtags(mt, geometry_xpath="/Xdmf/Domain/Grid[@Name='ball_d1']/Geometry")
Create a distributed (parallel) mesh with quadratic geometry¶
Generate mesh on rank 0, then build a distributed mesh
if MPI.COMM_WORLD.rank == 0:
geom = pygmsh.opencascade.Geometry()
ball = geom.add_ball([0.0, 0.0, 0.0], 1.0, char_length=0.2)
box = geom.add_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0])
cut = geom.boolean_difference([ball], [box])
geom.add_raw_code("Physical Surface(1) = {1};")
geom.add_physical(cut, 2)
pygmsh_mesh = pygmsh.generate_mesh(geom, extra_gmsh_arguments=["-order", "2"])
# Extract the topology and geometry data
cells, x = pygmsh_mesh.cells[-1].data, pygmsh_mesh.points
pygmsh_cell = pygmsh_mesh.cells[-1].type
# Extract marked facets
marked_entities = pygmsh_mesh.cells[-2].data
values = pygmsh_mesh.cell_data["gmsh:physical"][-2]
# Broadcast cell type data and geometric dimension
num_nodes = MPI.COMM_WORLD.bcast(cells.shape[1], root=0)
else:
num_nodes = MPI.COMM_WORLD.bcast(None, root=0)
cells, x = np.empty([0, num_nodes]), np.empty([0, 3])
marked_entities, values = np.empty((0, 6)), np.empty((0,))
# Permute the topology from GMSH to DOLFIN-X ordering
domain = ufl_mesh_from_gmsh("tetra10", 3)
gmsh_tetra10 = perm_gmsh(cpp.mesh.CellType.tetrahedron, 10)
cells = cells[:, gmsh_tetra10]
mesh = create_mesh(MPI.COMM_WORLD, cells, x, domain)
mesh.name = "ball_d2"
# Permute also entities which are tagged
gmsh_triangle6 = perm_gmsh(cpp.mesh.CellType.triangle, 6)
marked_entities = marked_entities[:, gmsh_triangle6]
local_entities, local_values = extract_local_entities(mesh, 2, marked_entities, values)
mesh.topology.create_connectivity(2, 0)
mt = create_meshtags(mesh, 2, cpp.graph.AdjacencyList_int32(local_entities), np.int32(local_values))
mt.name = "ball_d2_surface"
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "a") as file:
file.write_mesh(mesh)
mesh.topology.create_connectivity(2, 3)
file.write_meshtags(mt, geometry_xpath="/Xdmf/Domain/Grid[@Name='ball_d2']/Geometry")
if MPI.COMM_WORLD.rank == 0:
# Generate a mesh with 2nd-order hexahedral cells using pygmsh
geom = pygmsh.opencascade.Geometry()
geom.add_raw_code("Mesh.RecombineAll = 1;")
geom.add_raw_code("Mesh.CharacteristicLengthFactor = 1.0;")
geom.add_raw_code("Mesh.RecombinationAlgorithm = 2;")
circle = geom.add_disk([0.0, 0.0, 0.0], 1.0)
circle_inner = geom.add_disk([0.0, 0.0, 0.0], 0.5)
cut = geom.boolean_difference([circle], [circle_inner])
_, box, _ = geom.extrude(cut, translation_axis=[0.0, 0.0, 5], num_layers=5, recombine=True)
geom.add_physical(cut, label=1)
geom.add_physical(box, label=2)
pygmsh_mesh = pygmsh.generate_mesh(geom, extra_gmsh_arguments=["-order", "2"])
# Extract the topology and geometry data
cells, x = pygmsh_mesh.cells[-1].data, pygmsh_mesh.points
pygmsh_cell = pygmsh_mesh.cells[-1].type
# Extract marked facets
marked_entities = pygmsh_mesh.cells[-2].data
values = pygmsh_mesh.cell_data["gmsh:physical"][-2]
# Broadcast cell type data and geometric dimension
num_nodes = MPI.COMM_WORLD.bcast(cells.shape[1], root=0)
else:
# Receive cell type data and geometric dimension
num_nodes = MPI.COMM_WORLD.bcast(None, root=0)
cells, x = np.empty((0, num_nodes)), np.empty((0, 3))
marked_entities, values = np.empty((0, 9)), np.empty((0,))
# Permute the mesh topology from GMSH ordering to DOLFIN-X ordering
domain = ufl_mesh_from_gmsh("hexahedron27", 3)
gmsh_hex27 = perm_gmsh(cpp.mesh.CellType.hexahedron, 27)
cells = cells[:, gmsh_hex27]
mesh = create_mesh(MPI.COMM_WORLD, cells, x, domain)
mesh.name = "hex_d2"
# Permute also entities which are tagged
gmsh_quad9 = perm_gmsh(cpp.mesh.CellType.quadrilateral, 9)
marked_entities = marked_entities[:, gmsh_quad9]
local_entities, local_values = extract_local_entities(mesh, 2, marked_entities, values)
mesh.topology.create_connectivity(2, 0)
mt = create_meshtags(mesh, 2, cpp.graph.AdjacencyList_int32(local_entities), np.int32(local_values))
mt.name = "hex_d2_surface"
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "a") as file:
file.write_mesh(mesh)
mesh.topology.create_connectivity(2, 3)
file.write_meshtags(mt, geometry_xpath="/Xdmf/Domain/Grid[@Name='hex_d2']/Geometry")