Coverage for dragonfly_sim/utils/mesh_io

1"""

2Mesh input: dolfinx meshes from XDMF files or from structured CGNS grids.

4A structured multi-block CGNS grid of a 2D case is often stored as a 3D grid

5one cell thick in span (e.g. ONERA's OAT15A grids). read_cgns_planar reads it

6with PyVista (VTK's CGNS reader), keeps one span plane of every block as

7quadrilaterals, merges the nodes the blocks share, and finds the wall and far

8field edges from the CGNS boundary condition families. read_cgns_volume reads

9a 3D grid (e.g. the Simple Transonic Wing's wing_vol_L*.cgns) into one

10hexahedral mesh the same way. cgns_to_dolfinx_mesh turns either into a

11distributed dolfinx mesh in memory: no intermediate file.

13Boundary tags: load_dolfinx_mesh also reads the boundary facet tags a mesh

14file carries (an XDMF "Facet tags" grid, the physical groups of a gmsh .msh,

15the Wall and Farfield families of a planar CGNS grid) and stores them on the

16returned mesh with attach_facet_tags, where core/mesh_manager.py:Mesh finds

17them. They can define the wall and symmetry boundaries in place of a function

18of the facet midpoints (DG_windtunnel_model's wall_tags and symmetry_tags).

20Block interface nodes are merged by EXACT coordinate match, which holds for

21grids whose 1-to-1 interfaces were written point-matched (Pointwise exports

22are); a grid that needs a tolerance merge fails the boundary check below.

23"""

24import os

25import xml.etree.ElementTree as ET

27import numpy as np

28from mpi4py import MPI

29import basix.ufl

30import dolfinx

33def _plane_points(grid, plane_index):

34 """Points of k-plane `plane_index` of a StructuredGrid, shape (nj, ni, 3)."""

35 ni, nj, nk = grid.dimensions

36 if nk > 2:

37 raise ValueError("block has {} nodes in its third index direction; only grids one cell "

38 "thick (or planar) can be reduced to 2D".format(nk))

39 k = min(plane_index, nk - 1)

40 # VTK structured points run i fastest, then j, then k

41 return np.asarray(grid.points).reshape(nk, nj, ni, 3)[k]

44def _block_quads(nj, ni, offset):

45 """Quads of an nj x ni node block in dolfinx ordering [(i,j), (i+1,j), (i,j+1), (i+1,j+1)]."""

46 idx = offset + np.arange(nj * ni).reshape(nj, ni)

47 return np.stack([idx[:-1, :-1], idx[:-1, 1:], idx[1:, :-1], idx[1:, 1:]], axis=-1).reshape(-1, 4)

50def _signed_areas(points, quads):

51 """Signed areas of quads in dolfinx ordering (positive: counterclockwise)."""

52 p0, p1, p2, p3 = (points[quads[:, i]] for i in range(4))

53 d1, d2 = p3 - p0, p2 - p1

54 return 0.5 * (d1[:, 0] * d2[:, 1] - d1[:, 1] * d2[:, 0])

57def _quad_edges(quads):

58 """The four edges of every quad (dolfinx ordering), as sorted node pairs."""

59 e = np.concatenate([quads[:, [0, 1]], quads[:, [1, 3]], quads[:, [3, 2]], quads[:, [2, 0]]])

60 return np.sort(e, axis=1)

63def _family_patch_lines(reader, family, plane_index):

64 """Points of every boundary patch of `family`, cut at the span plane: a list

65 of (n, 3) arrays, one node line per patch, in order along the line."""

66 reader.disable_all_families()

67 reader.reader.SetFamilyArrayStatus(family, 1)

68 lines = []

69 for zone in reader.read()["Base"]:

70 if zone is None or "Patches" not in zone.keys():

71 continue

72 for patch in zone["Patches"]:

73 pts = _plane_points(patch, plane_index)

74 if min(pts.shape[:2]) != 1:

75 raise ValueError("a {} patch is not a line at the span plane (node grid {} x {}): "

76 "it lies in the plane itself".format(family, *pts.shape[:2]))

77 lines.append(pts.reshape(-1, 3))

78 return lines

81def merge_structured_blocks(blocks):

82 """

83 Merge 2D structured node blocks into one quad mesh.

85 blocks: list of (nj, ni, 2) node arrays. Returns (points (N, 2), quads

86 (M, 4) in dolfinx ordering, counterclockwise). Nodes shared by blocks are

87 merged by exact coordinate match; a block whose (i, j) axes are left-handed

88 has its quads flipped.

89 """

90 quads, offset = [], 0

91 raw = np.concatenate([b.reshape(-1, 2) for b in blocks])

92 for b in blocks:

93 nj, ni = b.shape[:2]

94 q = _block_quads(nj, ni, offset)

95 if _signed_areas(raw, q).sum() < 0:

96 q = q[:, [1, 0, 3, 2]]

97 quads.append(q)

98 offset += nj * ni

99 points, inverse = np.unique(raw, axis=0, return_inverse=True)

100 quads = inverse.reshape(-1)[np.concatenate(quads)]

101 return points, quads

102

103

104def read_cgns_planar(cgns_path, plane_index=0, wall_family="Wall", farfield_family="Farfield",

105 scale_to_chord=True):

106 """

107 Read a structured multi-block CGNS grid that is one cell thick (or planar)

108 and reduce it to a 2D quadrilateral mesh. Serial.

109

110 plane_index: which span plane of every block to keep (0 or 1). The span

111 direction is the coordinate axis that is constant on that plane; the other

112 two, in ascending order, become (x, y) -- an X-Z airfoil grid maps

113 (X, Z) -> (x, y).

114

115 scale_to_chord: shift the points so the wall's leading edge (its smallest

116 x) is at x = 0 and scale them by 1 / chord, the wall's x extent.

117

118 Returns a dict: points (N, 2), quads (M, 4) (dolfinx ordering,

119 counterclockwise), wall_edges and farfield_edges (node index pairs),

120 chord (in the file's units), n_raw_points (before merging).

121

122 Raises ValueError if a cell has non-positive area, an edge has more than

123 two cells, or the boundary edges are not exactly the wall and far field

124 edges.

125 """

126 import pyvista as pv

127

128 if not os.path.isfile(cgns_path):

129 raise FileNotFoundError(cgns_path)

130 reader = pv.get_reader(cgns_path)

131 planes = [_plane_points(zone["Internal"], plane_index) for zone in reader.read()["Base"]]

132

133 span_axes = [ax for ax in range(3) if all(np.ptp(p[..., ax]) == 0 for p in planes)]

134 if len(span_axes) != 1 or len({p[0, 0, span_axes[0]] for p in planes}) != 1:

135 raise ValueError("the span plane is not one coordinate plane shared by all blocks")

136 in_plane = [ax for ax in range(3) if ax != span_axes[0]]

137 points, quads = merge_structured_blocks([p[..., in_plane] for p in planes])

138 n_raw = sum(p.shape[0] * p.shape[1] for p in planes)

139

140 # boundary edges from the CGNS families, as merged node index pairs

141 lookup = {tuple(x): i for i, x in enumerate(points)}

142 reader.load_boundary_patch = True

143 edges = {}

144 for name, family in (("wall", wall_family), ("farfield", farfield_family)):

145 pairs = []

146 for line in _family_patch_lines(reader, family, plane_index):

147 try:

148 ids = np.array([lookup[tuple(x)] for x in line[:, in_plane]])

149 except KeyError:

150 raise ValueError("a {} patch node is not a grid node".format(family)) from None

151 pairs.append(np.column_stack([ids[:-1], ids[1:]]))

152 if not pairs:

153 raise ValueError("no boundary patches of family '{}' in {}".format(family, cgns_path))

154 edges[name] = np.concatenate(pairs)

155

156 # checks

157 areas = _signed_areas(points, quads)

158 if not np.all(areas > 0):

159 raise ValueError("{} cells have non-positive area".format(int(np.sum(areas <= 0))))

160 all_edges, counts = np.unique(_quad_edges(quads), axis=0, return_counts=True)

161 if counts.max() > 2:

162 raise ValueError("{} edges have more than two cells".format(int(np.sum(counts > 2))))

163 bdry = {tuple(e) for e in all_edges[counts == 1]}

164 tagged = [{tuple(e) for e in np.sort(edges[k], axis=1)} for k in ("wall", "farfield")]

165 if tagged[0] & tagged[1] or tagged[0] | tagged[1] != bdry:

166 raise ValueError("the boundary edges are not exactly the {} and {} edges ({} boundary, {} "

167 "{}, {} {})".format(wall_family, farfield_family, len(bdry), len(tagged[0]),

168 wall_family, len(tagged[1]), farfield_family))

169

170 wall_x = points[np.unique(edges["wall"]), 0]

171 chord = float(wall_x.max() - wall_x.min())

172 if scale_to_chord:

173 points = (points - np.array([wall_x.min(), 0.0])) / chord

174 return dict(points=points, quads=quads, wall_edges=edges["wall"],

175 farfield_edges=edges["farfield"], chord=chord, n_raw_points=n_raw)

176

177

178# The six faces of a hexahedron in dolfinx ordering, as local node indices

179_HEX_FACES = np.array([[0, 1, 2, 3], [4, 5, 6, 7], [0, 1, 4, 5],

180 [2, 3, 6, 7], [0, 2, 4, 6], [1, 3, 5, 7]])

181

182

183def _block_hexes(nk, nj, ni, offset):

184 """Hexes of an nk x nj x ni node block in dolfinx ordering (i fastest, then

185 j, then k: [(i,j,k), (i+1,j,k), (i,j+1,k), (i+1,j+1,k), the same at k+1])."""

186 idx = offset + np.arange(nk * nj * ni).reshape(nk, nj, ni)

187 bottom = [idx[:-1, :-1, :-1], idx[:-1, :-1, 1:], idx[:-1, 1:, :-1], idx[:-1, 1:, 1:]]

188 top = [idx[1:, :-1, :-1], idx[1:, :-1, 1:], idx[1:, 1:, :-1], idx[1:, 1:, 1:]]

189 return np.stack(bottom + top, axis=-1).reshape(-1, 8)

190

191

192def _center_jacobians(points, hexes):

193 """Jacobian determinants of the trilinear map at the cell centers, hexes in

194 dolfinx ordering (positive: right-handed; 8x the volume for a parallelepiped)."""

195 v = [points[hexes[:, i]] for i in range(8)]

196 d_i = (v[1] - v[0]) + (v[3] - v[2]) + (v[5] - v[4]) + (v[7] - v[6])

197 d_j = (v[2] - v[0]) + (v[3] - v[1]) + (v[6] - v[4]) + (v[7] - v[5])

198 d_k = (v[4] - v[0]) + (v[5] - v[1]) + (v[6] - v[2]) + (v[7] - v[3])

199 return np.einsum("ij,ij->i", d_i, np.cross(d_j, d_k)) / 64.0

200

201

202def _hex_faces(hexes):

203 """The six faces of every hex, as sorted node quadruples."""

204 return np.sort(hexes[:, _HEX_FACES].reshape(-1, 4), axis=1)

205

206

207def merge_structured_volume_blocks(blocks):

208 """

209 Merge 3D structured node blocks into one hexahedral mesh.

210

211 blocks: list of (nk, nj, ni, 3) node arrays. Returns (points (N, 3), hexes

212 (M, 8) in dolfinx ordering, right-handed). Nodes shared by blocks are

213 merged by exact coordinate match; a block whose (i, j, k) axes are

214 left-handed has its hexes mirrored in i.

215 """

216 hexes, offset = [], 0

217 raw = np.concatenate([b.reshape(-1, 3) for b in blocks])

218 for b in blocks:

219 nk, nj, ni = b.shape[:3]

220 h = _block_hexes(nk, nj, ni, offset)

221 if _center_jacobians(raw, h).sum() < 0:

222 h = h[:, [1, 0, 3, 2, 5, 4, 7, 6]]

223 hexes.append(h)

224 offset += nk * nj * ni

225 points, inverse = np.unique(raw, axis=0, return_inverse=True)

226 hexes = inverse.reshape(-1)[np.concatenate(hexes)]

227 return points, hexes

228

229

230def read_cgns_volume(cgns_path, interface_rel_tol=1e-10):

231 """

232 Read a structured multi-block 3D CGNS grid into one hexahedral mesh. Serial.

233

234 The blocks' nodes are merged by exact coordinate match (see

235 merge_structured_volume_blocks), in the file's units. The CGNS boundary

236 conditions are not read -- VTK skips BC_t types such as BCWallViscous and

237 BCFarfield -- so the solver finds its boundaries geometrically, as for an

238 untagged XDMF mesh.

239

240 interface_rel_tol: two boundary faces whose centroids are closer than this

241 times the grid's bounding-box diagonal are an unmerged block interface.

242

243 Returns a dict: points (N, 3), hexes (M, 8) (dolfinx ordering,

244 right-handed), boundary_faces (node index quadruples), n_raw_points

245 (before merging), n_blocks.

246

247 Raises ValueError if a cell is inverted or degenerate at its center, a

248 face has more than two cells, or block interfaces did not merge.

249 """

250 import pyvista as pv

251 import vtk

252 from scipy.spatial import cKDTree

253

254 if not os.path.isfile(cgns_path):

255 raise FileNotFoundError(cgns_path)

256 reader = pv.get_reader(cgns_path)

257 # VTK warns once per block and BC that it skips the BC_t node

258 warnings_on = vtk.vtkObject.GetGlobalWarningDisplay()

259 vtk.vtkObject.GlobalWarningDisplayOff()

260 try:

261 bases = reader.read()

262 finally:

263 vtk.vtkObject.SetGlobalWarningDisplay(warnings_on)

264 blocks = []

265 for base in bases:

266 for zone in base:

267 if zone is None:

268 continue

269 grid = zone["Internal"]

270 ni, nj, nk = grid.dimensions

271 if min(ni, nj, nk) < 2:

272 raise ValueError("a block has {} x {} x {} nodes: not a volume grid".format(ni, nj, nk))

273 # VTK structured points run i fastest, then j, then k

274 blocks.append(np.asarray(grid.points, dtype=np.float64).reshape(nk, nj, ni, 3))

275 if not blocks:

276 raise ValueError("no structured blocks in {}".format(cgns_path))

277 points, hexes = merge_structured_volume_blocks(blocks)

278 n_raw = sum(b.shape[0] * b.shape[1] * b.shape[2] for b in blocks)

279

280 jac = _center_jacobians(points, hexes)

281 if not np.all(jac > 0):

282 raise ValueError("{} cells are inverted or degenerate".format(int(np.sum(jac <= 0))))

283 faces, counts = np.unique(_hex_faces(hexes), axis=0, return_counts=True)

284 if counts.max() > 2:

285 raise ValueError("{} faces have more than two cells".format(int(np.sum(counts > 2))))

286 boundary_faces = faces[counts == 1]

287 diagonal = np.linalg.norm(np.ptp(points, axis=0))

288 pairs = cKDTree(points[boundary_faces].mean(axis=1)).query_pairs(interface_rel_tol * diagonal)

289 if pairs:

290 raise ValueError("{} pairs of boundary faces coincide: block interfaces whose nodes do not "

291 "match exactly".format(len(pairs)))

292 return dict(points=points, hexes=hexes, boundary_faces=boundary_faces, n_raw_points=n_raw,

293 n_blocks=len(blocks))

294

295

296def attach_facet_tags(mesh, facet_tags, tag_names=None):

297 """

298 Store boundary facet tags on a dolfinx mesh, where core/mesh_manager.py:Mesh

299 looks for them: `mesh.facet_tags`, a dolfinx MeshTags of the facets (or

300 None), and `mesh.facet_tag_names`, a dict {name: tag value} of the tags that

301 have a name (gmsh physical groups, CGNS families). Returns the mesh.

302 """

303 if facet_tags is not None and facet_tags.dim != mesh.topology.dim - 1:

304 raise ValueError("facet tags are defined on entities of dimension {}, not on facets "

305 "(dimension {})".format(facet_tags.dim, mesh.topology.dim - 1))

306 mesh.facet_tags = facet_tags

307 mesh.facet_tag_names = dict(tag_names or {})

308 return mesh

309

310

311def _cgns_facet_tags(mesh, edge_sets, comm):

312 """

313 MeshTags of the facets of a planar CGNS mesh from node index pairs, one

314 array per tag value ({value: (n, 2) array on rank 0, empty elsewhere}); the

315 pairs are in the input (merged-block) node numbering.

316 """

317 import dolfinx.io.gmsh

318 entities = np.concatenate([e for e in edge_sets.values()]).astype(np.int64)

319 values = np.concatenate([np.full(len(e), v, dtype=np.int32) for v, e in edge_sets.items()])

320 local_entities, local_values = dolfinx.io.gmsh.distribute_entity_data(mesh, 1, entities, values)

321 mesh.topology.create_connectivity(1, 2)

322 tags = dolfinx.mesh.meshtags_from_entities(

323 mesh, 1, dolfinx.graph.adjacencylist(local_entities), local_values.astype(np.int32, copy=False))

324 tags.name = "facet_tags"

325 return tags

326

327

328def cgns_to_dolfinx_mesh(cgns_path, comm=MPI.COMM_WORLD, tdim=2, **kwargs):

329 """

330 A distributed dolfinx mesh from a structured CGNS grid: rank 0 reads it

331 and dolfinx partitions it, ghosting cells across shared facets as

332 XDMFFile.read_mesh does (the interior-facet forms need it). Collective.

333

334 tdim=2: a grid one cell thick, reduced to quadrilaterals (read_cgns_planar,

335 kwargs passed on). info holds chord, n_points, n_cells, n_raw_points,

336 n_wall_edges and n_farfield_edges.

337 tdim=3: a volume grid, as hexahedra (read_cgns_volume, kwargs passed on).

338 info holds n_points, n_cells, n_raw_points, n_blocks and n_boundary_faces.

339

340 Returns (mesh, info), info being the same on every rank.

341 """

342 if tdim not in (2, 3):

343 raise ValueError("tdim must be 2 (planar grid) or 3 (volume grid), got {}".format(tdim))

344 if comm.rank == 0:

345 try:

346 data = (read_cgns_planar if tdim == 2 else read_cgns_volume)(cgns_path, **kwargs)

347 error = None

348 except Exception as e: # raise on every rank, not only rank 0

349 data, error = None, e

350 else:

351 data, error = None, None

352 error = comm.bcast(error, root=0)

353 if error is not None:

354 raise error

355 cells_key, nodes_per_cell, cell_name = (("quads", 4, "quadrilateral") if tdim == 2

356 else ("hexes", 8, "hexahedron"))

357 if comm.rank == 0:

358 cells, x = data[cells_key].astype(np.int64), data["points"]

359 if tdim == 2:

360 info = dict(chord=data["chord"], n_points=len(x), n_cells=len(cells),

361 n_raw_points=data["n_raw_points"], n_wall_edges=len(data["wall_edges"]),

362 n_farfield_edges=len(data["farfield_edges"]))

363 else:

364 info = dict(n_points=len(x), n_cells=len(cells), n_raw_points=data["n_raw_points"],

365 n_blocks=data["n_blocks"], n_boundary_faces=len(data["boundary_faces"]))

366 else:

367 cells, x, info = (np.empty((0, nodes_per_cell), dtype=np.int64), np.empty((0, tdim)),

368 None)

369 info = comm.bcast(info, root=0)

370 element = basix.ufl.element("Lagrange", cell_name, 1, shape=(tdim,))

371 partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet, 2)

372 mesh = dolfinx.mesh.create_mesh(comm, cells, element, x, partitioner=partitioner)

373 if tdim == 2:

374 # the families the planar reader found are the boundary tags: wall 1, far field 2

375 empty = np.empty((0, 2), dtype=np.int64)

376 edge_sets = {1: data["wall_edges"] if comm.rank == 0 else empty,

377 2: data["farfield_edges"] if comm.rank == 0 else empty}

378 names = {kwargs.get("wall_family", "Wall"): 1, kwargs.get("farfield_family", "Farfield"): 2}

379 attach_facet_tags(mesh, _cgns_facet_tags(mesh, edge_sets, comm), names)

380 else:

381 attach_facet_tags(mesh, None)

382 return mesh, info

383

384

385def _xdmf_grid_names(path):

386 """The names of the Grid nodes of an XDMF file (read by every rank)."""

387 return [g.get("Name") for g in ET.parse(path).getroot().iter("Grid")]

388

389

390def _read_xdmf(path, comm, mesh_name, facet_tags_name):

391 with dolfinx.io.XDMFFile(comm, path, "r") as xdmf:

392 mesh = xdmf.read_mesh(name=mesh_name)

393 mesh.topology.create_entities(mesh.topology.dim - 1)

394 grids = _xdmf_grid_names(path)

395 if facet_tags_name is not None and facet_tags_name not in grids:

396 raise ValueError("{} has no grid '{}' of facet tags (grids: {})".format(

397 path, facet_tags_name, grids))

398 name = facet_tags_name if facet_tags_name is not None else (

399 "Facet tags" if "Facet tags" in grids else None)

400 tags = xdmf.read_meshtags(mesh, name) if name is not None else None

401 return attach_facet_tags(mesh, tags)

402

403

404def _read_msh(path, comm, tdim):

405 import dolfinx.io.gmsh

406 partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet, 2)

407 data = dolfinx.io.gmsh.read_from_msh(path, comm, rank=0, gdim=tdim, partitioner=partitioner)

408 names = {name: group.tag for name, group in data.physical_groups.items()

409 if group.dim == data.mesh.topology.dim - 1}

410 return attach_facet_tags(data.mesh, data.facet_tags, names)

411

412

413def load_dolfinx_mesh(path, comm=MPI.COMM_WORLD, mesh_name="mesh", tdim=2, facet_tags_name=None,

414 **cgns_kwargs):

415 """

416 A dolfinx mesh from an XDMF file (grid `mesh_name`), a gmsh .msh file or a

417 structured CGNS grid (cgns_to_dolfinx_mesh with `tdim`, cgns_kwargs passed

418 on), by file extension. Returns (mesh, info); info is

419 cgns_to_dolfinx_mesh's dict, or {} for the others. Collective.

420

421 The boundary facet tags the file carries are stored on the mesh

422 (attach_facet_tags): mesh.facet_tags is a MeshTags or None, mesh.facet_tag_names

423 maps tag names to values where the format has names. XDMF: the grid named

424 `facet_tags_name` (default: "Facet tags" if the file has it); .msh: the

425 physical groups of facets, `tdim` being the mesh dimension; CGNS: the Wall

426 (1) and Farfield (2) families of a planar grid, none for a volume grid.

427 """

428 ext = os.path.splitext(path)[1].lower()

429 if ext == ".cgns":

430 return cgns_to_dolfinx_mesh(path, comm, tdim=tdim, **cgns_kwargs)

431 if ext == ".xdmf":

432 return _read_xdmf(path, comm, mesh_name, facet_tags_name), {}

433 if ext == ".msh":

434 return _read_msh(path, comm, tdim), {}

435 raise ValueError("unsupported mesh file '{}': expected .xdmf, .msh or .cgns".format(path))

Coverage for dragonfly_sim/utils/mesh_io_utils.py: 42%

227 statements