Coverage for dragonfly_sim/core/mesh_manager.py: 92%

1from dataclasses import dataclass 

2import warnings 

3 

4import numpy as np 

5from mpi4py import MPI 

6from petsc4py import PETSc 

7import ufl 

8import ufl.classes 

9from scipy.spatial import KDTree 

10import basix.ufl 

11import dolfinx 

12import dolfinx.fem.petsc 

13 

14from dragonfly_sim.utils.mpi_utils import allgather_numpy, assign_global_dof_offset 

15from dragonfly_sim.utils.meshwarping_utils import find_local_surface_nodes 

16from dragonfly_sim.core.form_manager import FormManager 

17from dragonfly_sim.utils.filewriter import FileWriter 

18 

19 

20@dataclass(frozen=True) 

21class BoundaryNodeSet: 

22 """Owned nodes on tagged boundary facets, numbered globally in rank-block order. 

23 

24 Built by ``Mesh.boundary_node_set``. ``local_idx`` are this rank's OWNED 

25 geometry-node indices; they carry the global ids ``global_idx`` = 

26 ``offset + arange(n_local)``, so the ranks' blocks tile ``[0, n_global)`` 

27 in rank order. ``coords`` and ``geom_gids`` (dolfinx global geometry ids) 

28 are replicated on every rank, row for row in that same order. 

29 shape_parameterization.WallFFD relies on this numbering. 

30 """ 

31 tags: tuple 

32 local_idx: np.ndarray 

33 offset: int 

34 n_global: int 

35 coords: np.ndarray 

36 geom_gids: np.ndarray 

37 

38 @property 

39 def n_local(self): 

40 return int(self.local_idx.size) 

41 

42 @property 

43 def global_idx(self): 

44 return np.arange(self.offset, self.offset + self.n_local, dtype=np.int64) 

45 

46 

47class Mesh(): 

48 # DOLFINx returns quadrilateral entity nodes in tensor ("Z") order, which 

49 # traverses the quad as a bowtie; this permutation makes it cyclic. 

50 _QUAD_CYCLIC_PERM = np.array([0, 1, 3, 2], dtype=np.int64) 

51 

52 def __init__(self, inp_mesh, inner_bdry_function=None, node_match_tol=1e-10, 

53 wall_tags=None, symmetry_tags=None, facet_tags=None): 

54 """ 

55 The wall (the body) is found in one of two ways: 

56 

57 inner_bdry_function a function of the boundary facet midpoints, shape 

58 (dim, n), returning a boolean mask: True for the wall. Used if given. 

59 wall_tags the tags of the mesh's own boundary facet tags that 

60 are the wall (a tag value or name, or a sequence of them), and 

61 symmetry_tags those that are symmetry planes. Used if there is no 

62 function. Every other boundary facet is far field. 

63 

64 The facet tags are `facet_tags` (a dolfinx MeshTags of the facets) or, 

65 if None, `inp_mesh.facet_tags`, which utils/mesh_io_utils.py:load_dolfinx_mesh 

66 sets (with `inp_mesh.facet_tag_names`, a {name: value} dict). If the 

67 function is used and the mesh also has tags, a UserWarning lists them. 

68 With the function, symmetry planes are not marked here (3D: y = 0 in 

69 DG_windtunnel_model); `symmetry_facet_mask` is None. 

70 """ 

71 # define input Dolfinx mesh object 

72 self.mesh = inp_mesh 

73 

74 self.baseline_nodes = self.mesh.geometry.x.copy() 

75 self.n_owned_nodes = int(self.mesh.geometry.index_map().size_local) 

76 self.node_match_tol = node_match_tol 

77 self.is_valid = True 

78 # design-parameter vector the current node coordinates were built for 

79 # (see apply_node_motions); None at baseline or when it is unknown 

80 self._deformation_design = None 

81 self._coordinate_space = None 

82 self._dof_to_node_map = None 

83 self._nodal_function = None 

84 self._quality = None 

85 self._node_sets = {} 

86 self.output_suffix = None 

87 self._deformation_writer = None 

88 self._quality_writers = None 

89 # callables run after every change of the node coordinates (see 

90 # add_geometry_listener) 

91 self._geometry_listeners = [] 

92 

93 self.n = ufl.FacetNormal(self.mesh) # cell boundary normal vector 

94 

95 # Volume-equivalent cell length |K|^(1/d). ufl.CellVolume is not 

96 # lowered for quads/hexes and FFCx cannot compile it, so it is written 

97 # pointwise from the Jacobian instead: exact on parallelogram cells, 

98 # the pointwise density elsewhere. Being a plain Jacobian expression, 

99 # it also differentiates w.r.t. the mesh coordinates, which an 

100 # assembled DG0 field would not. 

101 tdim = self.mesh.topology.dim 

102 self.cell_volume = (abs(ufl.JacobianDeterminant(self.mesh)) 

103 * ufl.classes.ReferenceCellVolume(self.mesh)) 

104 self.h_vol = self.cell_volume ** (1.0 / tdim) 

105 # Shortest edge, for the local pseudo-time step of PTC 

106 # (core/time_integration.py); never part of the residual. 

107 self.h_min = ufl.MinCellEdgeLength(self.mesh) 

108 

109 # The boundary queries below, and the wall-node lookups done later on 

110 # (entities_to_geometry) need these to exist already: the facet 

111 # entities, facet-to-cell connectivity (which also stores its 

112 # transpose), facet-to-vertex and vertex-to-cell connectivity 

113 self.mesh.topology.create_entities(tdim - 1) 

114 self.mesh.topology.create_connectivity(tdim - 1, tdim) 

115 self.mesh.topology.create_connectivity(tdim - 1, 0) 

116 self.mesh.topology.create_connectivity(0, tdim) 

117 

118 # we compute which facets form the mesh boundaries and their midpoints 

119 self.bdry_facets = dolfinx.mesh.exterior_facet_indices(self.mesh.topology) 

120 self.bdry_midpoints = dolfinx.mesh.compute_midpoints( 

121 self.mesh, self.mesh.topology.dim-1, self.bdry_facets) 

122 

123 # from all boundary facets, isolate the boundary of the wing/airfoil/etc. 

124 self.inner_bdry_facet_mask, self.symmetry_facet_mask = self._wall_and_symmetry_masks( 

125 inner_bdry_function, wall_tags, symmetry_tags, facet_tags) 

126 

127 # Initialize the meshtags object; used to define integration domains for boundary terms 

128 self.meshtags = None 

129 

130 # Instantiate FormManager for handling UFL forms and integration measures 

131 self.form_manager = FormManager(self) 

132 

133 # ------------------------------------------------------------------ 

134 # Wall and symmetry boundaries: from a function or from the mesh's tags 

135 # ------------------------------------------------------------------ 

136 def _tag_summary(self, facet_tags, names): 

137 """ 

138 '{value: (name, number of boundary facets)}' lines for the tags on the 

139 exterior facets, the same on every rank. Owned facets only, so a facet 

140 shared between ranks is counted once. 

141 """ 

142 comm = self.mesh.comm 

143 n_owned = self.mesh.topology.index_map(self.mesh.topology.dim - 1).size_local 

144 on_bdry = np.isin(facet_tags.indices, self.bdry_facets) & (facet_tags.indices < n_owned) 

145 values, counts = np.unique(facet_tags.values[on_bdry], return_counts=True) 

146 total = {} 

147 for v_c in comm.allgather(dict(zip(values.tolist(), counts.tolist()))): 

148 for v, c in v_c.items(): 

149 total[v] = total.get(v, 0) + c 

150 by_value = {v: n for n, v in names.items()} 

151 return {v: (by_value.get(v), total[v]) for v in sorted(total)} 

152 

153 @staticmethod 

154 def _describe(summary): 

155 return ", ".join("{}{} ({} facets)".format(v, " '{}'".format(n) if n else "", c) 

156 for v, (n, c) in summary.items()) 

157 

158 def _resolve_tags(self, tags, summary, names, what): 

159 """Tag values from a value, a name, or a sequence of them; checked to be present.""" 

160 if tags is None: 

161 return set() 

162 tags = [tags] if isinstance(tags, (int, np.integer, str)) else list(tags) 

163 values = set() 

164 for t in tags: 

165 if isinstance(t, str): 

166 if t not in names: 

167 raise ValueError("{}: the mesh has no boundary tag named '{}' (named tags: {})".format( 

168 what, t, sorted(names) or "none")) 

169 t = names[t] 

170 t = int(t) 

171 if t not in summary: 

172 raise ValueError("{}: tag {} is not on the mesh's boundary facets (tags present: {})".format( 

173 what, t, self._describe(summary) or "none")) 

174 values.add(t) 

175 return values 

176 

177 def _wall_and_symmetry_masks(self, inner_bdry_function, wall_tags, symmetry_tags, facet_tags): 

178 """(wall mask, symmetry mask or None) over self.bdry_facets; see __init__.""" 

179 facet_tags = facet_tags if facet_tags is not None else getattr(self.mesh, "facet_tags", None) 

180 names = getattr(self.mesh, "facet_tag_names", None) or {} 

181 summary = self._tag_summary(facet_tags, names) if facet_tags is not None else {} 

182 

183 if inner_bdry_function is not None: 

184 if self.mesh.comm.rank == 0 and (summary or wall_tags is not None 

185 or symmetry_tags is not None): 

186 if summary: 

187 message = ("The input mesh has boundary facet tags: {}. The boundary function " 

188 "is used to find the wall, so the boundary facet tags are " 

189 "currently unused.{} To use the tags " 

190 "instead, leave out the boundary function and pass wall_tags (the " 

191 "tags of the body) and, if there is a symmetry plane, " 

192 "symmetry_tags; all other boundaries are far field.").format( 

193 self._describe(summary), 

194 " wall_tags and symmetry_tags are ignored." if ( 

195 wall_tags is not None or symmetry_tags is not None) else "") 

196 else: 

197 message = ("wall_tags and symmetry_tags are ignored: the boundary function " 

198 "defines the wall.") 

199 warnings.warn(message, UserWarning, stacklevel=4) 

200 return np.asarray(inner_bdry_function(self.bdry_midpoints.T), dtype=bool), None 

201 

202 if wall_tags is None: 

203 if summary: 

204 raise ValueError( 

205 "No boundary function was given, and wall_tags was not set. The mesh has " 

206 "boundary facet tags: {}. Pass wall_tags (and symmetry_tags for symmetry " 

207 "planes); the other boundaries are far field.".format(self._describe(summary))) 

208 raise ValueError("Provide either a boundary function (facet midpoints -> wall mask) " 

209 "or wall_tags with a mesh that has boundary facet tags.") 

210 if facet_tags is None: 

211 raise ValueError("wall_tags was given but the mesh has no boundary facet tags (read it " 

212 "with load_dolfinx_mesh, or pass facet_tags)") 

213 wall = self._resolve_tags(wall_tags, summary, names, "wall_tags") 

214 sym = self._resolve_tags(symmetry_tags, summary, names, "symmetry_tags") 

215 if not wall: 

216 raise ValueError("wall_tags is empty") 

217 if wall & sym: 

218 raise ValueError("tags {} are both wall and symmetry".format(sorted(wall & sym))) 

219 

220 # tag value of every boundary facet of this rank (0: untagged) 

221 n_facets = self.mesh.topology.index_map(self.mesh.topology.dim - 1) 

222 value = np.zeros(n_facets.size_local + n_facets.num_ghosts, dtype=np.int64) 

223 value[facet_tags.indices] = facet_tags.values 

224 value = value[self.bdry_facets] 

225 wall_mask = np.isin(value, list(wall)) 

226 sym_mask = np.isin(value, list(sym)) if sym else None 

227 return wall_mask, sym_mask 

228 

229 def tag_boundary_facets(self, bdry_facet_value_tuples): 

230 # Construct the boundary mesh tags from a list of tuples (boundary facets, domain_boundary_ID) 

231 values = np.zeros(self.bdry_midpoints.shape[0], dtype=np.int32) 

232 for facets, value in bdry_facet_value_tuples: 

233 values[facets] = value 

234 

235 self.meshtags = dolfinx.mesh.meshtags( 

236 self.mesh, self.mesh.topology.dim - 1, self.bdry_facets, values) 

237 return self.meshtags 

238 

239 def _wall_boundary_node_coords(self, wall_tag_value): 

240 if self.meshtags is None: 

241 raise RuntimeError("Mesh tags not yet constructed. " 

242 "Call tag_boundary_facets first.") 

243 

244 wall_facet_indices = self.meshtags.find(wall_tag_value) 

245 

246 self.mesh.topology.create_connectivity(self.mesh.topology.dim - 1, 0) 

247 facet_to_node = self.mesh.topology.connectivity(self.mesh.topology.dim - 1, 0) 

248 

249 if len(wall_facet_indices) > 0: 

250 wall_vertex_topo_idxs = np.unique(np.concatenate( 

251 [facet_to_node.links(facet) for facet in wall_facet_indices])) 

252 wall_node_geom_idxs = dolfinx.cpp.mesh.entities_to_geometry( 

253 self.mesh._cpp_object, 0, wall_vertex_topo_idxs, False) 

254 return self.mesh.geometry.x[wall_node_geom_idxs[:, 0], :self.mesh.topology.dim] 

255 

256 return np.empty(shape=(0, self.mesh.topology.dim)) 

257 

258 def _global_bounding_box(self, coords): 

259 if coords.shape[0] > 0: 

260 local_mins = coords.min(axis=0) 

261 local_maxs = coords.max(axis=0) 

262 else: 

263 local_mins = np.full(self.mesh.topology.dim, np.inf) 

264 local_maxs = np.full(self.mesh.topology.dim, -np.inf) 

265 

266 return [(self.mesh.comm.allreduce(local_mins[i], op=MPI.MIN), 

267 self.mesh.comm.allreduce(local_maxs[i], op=MPI.MAX)) 

268 for i in range(self.mesh.topology.dim)] 

269 

270 def compute_wall_bounding_box(self, wall_tag_value): 

271 # Safe to call every iteration on the current (possibly deformed) mesh. 

272 coords = self._wall_boundary_node_coords(wall_tag_value) 

273 return self._global_bounding_box(coords) 

274 

275 

276 # ================================================================== 

277 # Geometry nodes and the CG1 coordinate space 

278 # ================================================================== 

279 

280 @property 

281 def coordinate_space(self): 

282 """Vector CG1 space whose dofs coincide with the geometry nodes. 

283 

284 The space mesh-coordinate derivatives (dR/dx_mesh) are taken in. 

285 Built on first use. 

286 """ 

287 if self._coordinate_space is None: 

288 element = basix.ufl.element( 

289 "Lagrange", self.mesh.basix_cell(), 1, 

290 shape=(self.mesh.geometry.dim,), dtype=dolfinx.default_real_type) 

291 self._coordinate_space = dolfinx.fem.functionspace(self.mesh, element) 

292 return self._coordinate_space 

293 

294 @property 

295 def dof_to_node_map(self): 

296 """(owned CG1 dofs) x (owned geometry nodes) 0/1 permutation matrix. 

297 

298 The CG1 vector dofs correspond 1:1 with the geometry nodes but are 

299 numbered differently: ``dof_to_node_map @ node_values`` reorders an 

300 (n_owned_nodes, gdim) nodal array into dof order, and its transpose 

301 maps dof-ordered derivatives back to geometry-node order. Matched by 

302 coordinates within ``node_match_tol`` -- an exact-match test with room 

303 for round-off. Built on first use. 

304 """ 

305 if self._dof_to_node_map is None: 

306 n = self.n_owned_nodes 

307 dof_tree = KDTree(self.coordinate_space.tabulate_dof_coordinates()[:n, :]) 

308 node_tree = KDTree(self.mesh.geometry.x[:n, :]) 

309 perm = dof_tree.sparse_distance_matrix( 

310 node_tree, max_distance=self.node_match_tol).tocsr() 

311 perm.data[:] = 1. 

312 self._dof_to_node_map = perm 

313 return self._dof_to_node_map 

314 

315 def nodal_function(self, node_values): 

316 """An (n_owned_nodes, gdim) nodal array as a Function on coordinate_space. 

317 

318 The Function is reused between calls. 

319 """ 

320 V = self.coordinate_space 

321 if self._nodal_function is None: 

322 self._nodal_function = dolfinx.fem.Function(V) 

323 n_local_dofs = V.dofmap.index_map.size_local * V.dofmap.index_map_bs 

324 self._nodal_function.x.array[:n_local_dofs] = \ 

325 (self.dof_to_node_map @ node_values).reshape(-1) 

326 self._nodal_function.x.scatter_forward() 

327 return self._nodal_function 

328 

329 # ================================================================== 

330 # Deformation 

331 # ================================================================== 

332 

333 def add_geometry_listener(self, callback): 

334 """Run `callback()` (collective) after every reset_nodes and every 

335 apply_node_motions that moves the nodes (not one skipped for an 

336 unchanged design), e.g. to refresh fields computed from the node 

337 coordinates (cell centres, wall distance).""" 

338 self._geometry_listeners.append(callback) 

339 

340 def _notify_geometry_listeners(self): 

341 for callback in self._geometry_listeners: 

342 callback() 

343 

344 @property 

345 def deformation_design(self): 

346 """Copy of the design vector the mesh is currently deformed for, or None.""" 

347 return None if self._deformation_design is None else self._deformation_design.copy() 

348 

349 def reset_nodes(self, notify=True): 

350 """Put every local node, ghosts included, back at its baseline position. 

351 

352 The callbacks no longer call this: the mesh stays deformed until 

353 apply_node_motions is given a different design. 

354 """ 

355 gdim = self.mesh.geometry.dim 

356 self.mesh.geometry.x[:, :gdim] = self.baseline_nodes[:, :gdim] 

357 self._deformation_design = None 

358 if notify: 

359 self._notify_geometry_listeners() 

360 

361 def _is_deformed_for(self, design): 

362 """True on every rank if the mesh already holds the deformation of ``design``. 

363 

364 ``design`` is replicated across ranks, but the comparison is reduced 

365 anyway: apply_node_motions is collective, so the ranks must agree on 

366 whether to skip it. 

367 """ 

368 stored = self._deformation_design 

369 same = (stored is not None and stored.shape == design.shape 

370 and bool(np.array_equal(stored, design))) 

371 return bool(self.mesh.comm.allreduce(same, op=MPI.MIN)) 

372 

373 def apply_node_motions(self, node_motions, design=None): 

374 """Move the mesh to baseline + ``node_motions`` and re-check its cells. 

375 

376 ``node_motions`` is (n_owned_nodes, gdim) in geometry-node order. 

377 Starting from the baseline makes the result depend on this call only, 

378 so repeated calls do not accumulate. Ghost nodes are updated from 

379 their owners. Sets ``is_valid`` (no inverted or zero-volume cell) and 

380 returns it. Collective. 

381 

382 ``design`` is the (replicated) design-parameter vector that 

383 ``node_motions`` was computed from. The mesh stays deformed after this 

384 call and remembers ``design``; a later call with an equal ``design`` 

385 returns the current ``is_valid`` without touching the mesh, so the 

386 deformation is built once per design. A different ``design`` resets 

387 the nodes and warps again. ``design=None`` always warps. 

388 """ 

389 if design is not None: 

390 design = np.array(design, dtype=np.float64, copy=True).ravel() 

391 if self._is_deformed_for(design): 

392 PETSc.Sys.Print("Mesh already deformed for this design: reusing it") 

393 return self.is_valid 

394 

395 # unknown until the new deformation is complete, so a failure part-way 

396 # never leaves a key that matches a half-moved mesh 

397 self._deformation_design = None 

398 self.reset_nodes(notify=False) 

399 n, gdim = self.n_owned_nodes, self.mesh.geometry.dim 

400 x = self.mesh.geometry.x 

401 x[:n, :gdim] += node_motions[:, :gdim] 

402 

403 # geometry.x has no scatter_forward of its own: route the owned 

404 # coordinates through a dolfinx.la.Vector on the geometry index map. 

405 geom_vec = dolfinx.la.vector(self.mesh.geometry.index_map(), bs=3) 

406 geom_vec.array[:n * 3] = x[:n, :].ravel() 

407 geom_vec.scatter_forward() 

408 x[:] = geom_vec.array.reshape(-1, 3) 

409 

410 self._update_quality() 

411 self._notify_geometry_listeners() 

412 self._deformation_design = design 

413 return self.is_valid 

414 

415 # ================================================================== 

416 # Mesh quality and output files 

417 # ================================================================== 

418 

419 def _setup_quality(self): 

420 """DG0 fields for the quality metrics apply_node_motions reports.""" 

421 if self._quality is not None: 

422 return self._quality 

423 msh = self.mesh 

424 V0 = dolfinx.fem.functionspace(msh, ("DG", 0)) 

425 pts = V0.element.interpolation_points 

426 q = {"V": V0, 

427 "vol_form": dolfinx.fem.form(ufl.TestFunction(V0) * ufl.dx), 

428 "vol": dolfinx.fem.Function(V0)} 

429 # Per-cell shape metric, as distinct from the whole-mesh max/min edge 

430 # ratio also reported: it is the per-cell value that says how badly an 

431 # isotropic length scale misrepresents a cell. 

432 for name, expr in (("jac", ufl.JacobianDeterminant(msh)), 

433 ("min_edge", ufl.MinCellEdgeLength(msh)), 

434 ("max_edge", ufl.MaxCellEdgeLength(msh)), 

435 ("cell_aspect_ratio", 

436 ufl.MaxCellEdgeLength(msh) / ufl.MinCellEdgeLength(msh))): 

437 q[name + "_expr"] = dolfinx.fem.Expression(expr, pts) 

438 q[name] = dolfinx.fem.Function(V0) 

439 self._quality = q 

440 return q 

441 

442 def _update_quality(self): 

443 """Evaluate, print and (if configured) write the quality metrics; set is_valid.""" 

444 q = self._setup_quality() 

445 comm = self.mesh.comm 

446 with q["vol"].x.petsc_vec.localForm() as loc: 

447 loc.set(0.0) 

448 dolfinx.fem.petsc.assemble_vector(q["vol"].x.petsc_vec, q["vol_form"]) 

449 q["vol"].x.scatter_forward() 

450 for name in ("jac", "min_edge", "max_edge", "cell_aspect_ratio"): 

451 q[name].interpolate(q[name + "_expr"]) 

452 

453 if self._quality_writers is not None: 

454 for name, writer in self._quality_writers.items(): 

455 writer.interpolate_and_write(q[name]) 

456 

457 def global_stat(name, fn, op, empty): 

458 a = q[name].x.array 

459 return comm.allreduce(fn(a) if a.size > 0 else empty, op=op) 

460 

461 min_vol = global_stat("vol", np.min, MPI.MIN, np.inf) 

462 max_vol = global_stat("vol", np.max, MPI.MAX, -np.inf) 

463 min_jac = global_stat("jac", np.min, MPI.MIN, np.inf) 

464 max_jac = global_stat("jac", np.max, MPI.MAX, -np.inf) 

465 min_edge = global_stat("min_edge", np.min, MPI.MIN, np.inf) 

466 max_edge = global_stat("max_edge", np.max, MPI.MAX, -np.inf) 

467 max_cell_ar = global_stat("cell_aspect_ratio", np.max, MPI.MAX, -np.inf) 

468 cell_ar_sum = global_stat("cell_aspect_ratio", np.sum, MPI.SUM, 0.0) 

469 n_cells = comm.allreduce(q["cell_aspect_ratio"].x.array.size, op=MPI.SUM) 

470 mean_cell_ar = cell_ar_sum / n_cells if n_cells > 0 else np.inf 

471 aspect_ratio = max_edge / min_edge if min_edge > 0 else np.inf 

472 

473 PETSc.Sys.Print("--- Mesh Quality Metrics ---") 

474 PETSc.Sys.Print("Min cell volume: {:.6e}".format(min_vol)) 

475 PETSc.Sys.Print("Max cell volume: {:.6e}".format(max_vol)) 

476 PETSc.Sys.Print("Min Jacobian determinant: {:.6e}".format(min_jac)) 

477 PETSc.Sys.Print("Max Jacobian determinant: {:.6e}".format(max_jac)) 

478 PETSc.Sys.Print("Min edge length: {:.6e}".format(min_edge)) 

479 PETSc.Sys.Print("Max edge length: {:.6e}".format(max_edge)) 

480 PETSc.Sys.Print("Max per-cell aspect ratio: {:.6e}".format(max_cell_ar)) 

481 PETSc.Sys.Print("Mean per-cell aspect ratio: {:.6e}".format(mean_cell_ar)) 

482 PETSc.Sys.Print("Max Global Aspect Ratio proxy (max_edge/min_edge over whole mesh): {:.6e}".format(aspect_ratio)) 

483 self.is_valid = bool(min_jac > 0) 

484 if not self.is_valid: 

485 PETSc.Sys.Print("WARNING: Mesh has inverted or zero-volume cells!") 

486 PETSc.Sys.Print("----------------------------") 

487 

488 def configure_output(self, suffix, deformation=True, quality=False): 

489 """Enable the per-evaluation mesh output files, named with ``suffix``. 

490 

491 ``deformation`` writes mesh_deformation_<suffix>.bp (through 

492 write_deformation_output) on the baseline mesh, so ParaView's Warp By 

493 Vector with it reproduces the deformed mesh; ``quality`` writes the DG0 

494 quality fields mesh_{vol,jac,min_edge,max_edge,cell_aspect_ratio}_<suffix>.bp 

495 on the deformed mesh on every apply_node_motions. ``suffix`` also names 

496 the memory logs of the classes that use this mesh. Collective. 

497 """ 

498 comm = self.mesh.comm 

499 self.output_suffix = suffix 

500 self._deformation_writer = ( 

501 FileWriter("mesh_deformation_{}".format(suffix), comm, self.coordinate_space, 

502 mesh_deformation=False) 

503 if deformation else None) 

504 self._quality_writers = None 

505 if quality: 

506 V0 = self._setup_quality()["V"] 

507 self._quality_writers = { 

508 name: FileWriter("mesh_{}_{}".format(name, suffix), comm, V0, mesh_deformation=True) 

509 for name in ("vol", "jac", "min_edge", "max_edge", "cell_aspect_ratio")} 

510 

511 def write_deformation_output(self, node_motions, write_counter): 

512 """Write (n_owned_nodes, gdim) node motions to mesh_deformation_<suffix>.bp. 

513 

514 A no-op unless configure_output enabled it. Called once per design 

515 evaluation with write_counter = eval_idx, so this file and the solution 

516 file share one index axis. 

517 """ 

518 if self._deformation_writer is None: 

519 return 

520 self._deformation_writer.interpolate_and_write( 

521 self.nodal_function(node_motions), write_counter=write_counter) 

522 

523 # ================================================================== 

524 # Tagged boundary node sets and facets 

525 # ================================================================== 

526 

527 def boundary_node_set(self, tags, exclude=None): 

528 """Owned nodes on facets tagged with any of ``tags``, numbered globally. 

529 

530 ``exclude`` (another BoundaryNodeSet) removes its nodes -- e.g. the 

531 wall from the farfield anchors, so a node shared by both stays with the 

532 wall. See BoundaryNodeSet for the numbering. Cached per (tags, 

533 exclude), so every caller gets the same numbering. Collective on first 

534 call; requires tag_boundary_facets. 

535 """ 

536 if self.meshtags is None: 

537 raise RuntimeError("Mesh tags not yet constructed. " 

538 "Call tag_boundary_facets first.") 

539 tags = tuple(int(t) for t in tags) 

540 key = (tags, None if exclude is None else exclude.tags) 

541 if key in self._node_sets: 

542 return self._node_sets[key] 

543 

544 msh = self.mesh 

545 comm = msh.comm 

546 gdim = msh.geometry.dim 

547 self._check_vertex_geometry_numbering() 

548 

549 parts = [find_local_surface_nodes(msh, self.meshtags, tag) for tag in tags] 

550 local_idx = (np.unique(np.concatenate(parts)) if parts 

551 else np.zeros((0,), dtype=np.int64)) 

552 if exclude is not None: 

553 local_idx = np.setdiff1d(local_idx, exclude.local_idx) 

554 local_idx = local_idx.astype(np.int64) 

555 

556 offset, n_global = assign_global_dof_offset(comm, local_idx.size) 

557 coords = allgather_numpy( 

558 comm, self.baseline_nodes[:self.n_owned_nodes, :gdim][local_idx].ravel(), 

559 n_global * gdim).reshape(n_global, gdim) 

560 # comm.allgather rather than allgather_numpy: these are int64 global 

561 # ids and allgather_numpy is float64-only. 

562 gids = msh.geometry.index_map().local_to_global( 

563 local_idx.astype(np.int32)).astype(np.int64) 

564 geom_gids = np.concatenate(comm.allgather(gids) + [np.zeros(0, dtype=np.int64)]) 

565 if np.unique(geom_gids).size != geom_gids.size: 

566 raise RuntimeError( 

567 "Mesh.boundary_node_set{}: a geometry node appears twice; each " 

568 "node should be owned by exactly one rank.".format(tags)) 

569 

570 node_set = BoundaryNodeSet(tags=tags, local_idx=local_idx, offset=int(offset), 

571 n_global=int(n_global), coords=coords, 

572 geom_gids=geom_gids) 

573 self._node_sets[key] = node_set 

574 return node_set 

575 

576 def _check_vertex_geometry_numbering(self): 

577 """Owned vertex indices must equal owned geometry-node indices. 

578 

579 find_local_surface_nodes works in topology vertex indices, the facet 

580 connectivity in geometry node indices. For first-order meshes the two 

581 coincide over the OWNED range, and only that range is used, so the node 

582 sets and the connectivity stay consistent. 

583 

584 They do NOT coincide over the ghost range -- every rank has a few 

585 hundred ghosts numbered differently in the two index spaces. That is 

586 fine: ghosts only enter through the facet connectivity, which is 

587 converted to global geometry ids via geometry.index_map().local_to_global 

588 before being used, and that map is valid for ghosts. 

589 """ 

590 msh = self.mesh 

591 n_owned_verts = msh.topology.index_map(0).size_local 

592 msh.topology.create_connectivity(0, msh.topology.dim) 

593 vertex_to_geom = dolfinx.mesh.entities_to_geometry( 

594 msh, 0, np.arange(n_owned_verts, dtype=np.int32)).reshape(-1) 

595 if not np.array_equal(vertex_to_geom, np.arange(n_owned_verts, dtype=np.int32)): 

596 raise RuntimeError( 

597 "Mesh: owned vertex indices and owned geometry node indices do " 

598 "not coincide (they do for first-order dolfinx meshes), so " 

599 "boundary node sets and facet connectivity would be indexed " 

600 "inconsistently.") 

601 

602 def tagged_facets(self, tags, node_gids): 

603 """This rank's OWNED facets tagged with any of ``tags``, as node-id rows. 

604 

605 A node id is the node's position in ``node_gids``, a replicated array 

606 of global geometry ids (a BoundaryNodeSet's ``geom_gids``, or several 

607 concatenated). Quad facets come back in cyclic vertex order. Only owned 

608 facets are taken: a facet shared between two ranks appears in both 

609 ranks' meshtags, and counting it twice would double it. 

610 

611 Returns ``(conn, facets)``: (n_facets, nodes_per_facet) int64 node ids 

612 and the local facet indices. 

613 """ 

614 msh = self.mesh 

615 fdim = msh.topology.dim - 1 

616 msh.topology.create_entities(fdim) 

617 msh.topology.create_connectivity(fdim, msh.topology.dim) 

618 msh.topology.create_connectivity(fdim, 0) 

619 empty = (np.zeros((0, 0), dtype=np.int64), np.zeros((0,), dtype=np.int32)) 

620 

621 parts = [self.meshtags.find(int(tag)) for tag in tags] 

622 if not parts: 

623 return empty 

624 facets = np.unique(np.concatenate(parts)).astype(np.int32) 

625 facets = facets[facets < msh.topology.index_map(fdim).size_local] 

626 if facets.size == 0: 

627 return empty 

628 

629 facet_geom = dolfinx.mesh.entities_to_geometry(msh, fdim, facets, False) 

630 gids = msh.geometry.index_map().local_to_global( 

631 facet_geom.reshape(-1).astype(np.int32)) 

632 

633 node_gids = np.asarray(node_gids, dtype=np.int64) 

634 order = np.argsort(node_gids) 

635 sorted_gids = node_gids[order] 

636 if np.unique(sorted_gids).size != sorted_gids.size: 

637 raise RuntimeError("Mesh.tagged_facets: node_gids contains duplicates.") 

638 pos = np.clip(np.searchsorted(sorted_gids, gids), 0, max(sorted_gids.size - 1, 0)) 

639 found = sorted_gids.size > 0 and np.all(sorted_gids[pos] == gids) 

640 if not found: 

641 raise RuntimeError( 

642 "Mesh.tagged_facets: {} facet node(s) are not in node_gids. This " 

643 "happens when a rank owns a tagged facet whose node is owned by a " 

644 "rank that saw no facet with that tag.".format( 

645 int((sorted_gids[pos] != gids).sum()) if sorted_gids.size else gids.size)) 

646 conn = order[pos].reshape(facet_geom.shape).astype(np.int64) 

647 if fdim == 2 and conn.shape[1] == 4: 

648 conn = conn[:, self._QUAD_CYCLIC_PERM] 

649 return conn, facets 

650 

651 def facet_cell_midpoints(self, facets): 

652 """(n_facets, 3) baseline midpoint of the single cell attached to each exterior facet.""" 

653 msh = self.mesh 

654 tdim = msh.topology.dim 

655 f_to_c = msh.topology.connectivity(tdim - 1, tdim) 

656 c_to_g = msh.geometry.dofmap 

657 x = self.baseline_nodes 

658 mids = np.empty((facets.size, 3), dtype=np.float64) 

659 for k, f in enumerate(facets): 

660 mids[k] = x[c_to_g[f_to_c.links(int(f))[0]]].mean(axis=0) 

661 return mids 

662 

663 @staticmethod 

664 def newell_normals(face_coords): 

665 """Newell's method: (n_faces, k, 3) -> (n_faces, 3), twice the area vector.""" 

666 rolled = np.roll(face_coords, -1, axis=1) 

667 return np.cross(face_coords, rolled).sum(axis=1) 

668 

669 @staticmethod 

670 def orient_facets(conn, coords, cell_mids): 

671 """Reverse the traversal of faces whose Newell normal points into the fluid. 

672 

673 After this every face normal points out of its adjacent fluid cell. 

674 DOLFINx orders facet vertices from the sorted global vertex numbering, 

675 which fixes neither the traversal order nor a consistent winding. 

676 Measured on the 3D wing mesh: without this, only ~51% of facets come 

677 out consistently oriented. Returns ``(conn, n_flipped)``. 

678 """ 

679 face_coords = coords[conn] 

680 normals = Mesh.newell_normals(face_coords) 

681 outward = np.einsum("ij,ij->i", normals, face_coords.mean(axis=1) - cell_mids) 

682 flip = outward < 0.0 

683 conn = conn.copy() 

684 conn[flip] = conn[flip][:, ::-1] 

685 return conn, int(flip.sum()) 

686 

687 def wall_facets(self, wall_tag): 

688 """Replicated wall facets as rows of wall node ids in ``[0, n_global)``. 

689 

690 Node ids are the ``boundary_node_set((wall_tag,))`` numbering. Segments 

691 in 2D, triangles or quads (cyclic vertex order) in 3D, in rank order, so 

692 every rank holds the same array. Consistently wound: the facet normal -- 

693 Newell's in 3D, ``(dy, -dx)`` of the traversal in 2D -- points out of 

694 the adjacent fluid cell, i.e. into the body. Collective. 

695 """ 

696 comm = self.mesh.comm 

697 gdim = self.mesh.geometry.dim 

698 wall = self.boundary_node_set((wall_tag,)) 

699 local_faces, local_facets = self.tagged_facets((wall_tag,), wall.geom_gids) 

700 if local_faces.size: 

701 coords = np.zeros((wall.n_global, 3), dtype=np.float64) 

702 coords[:, :gdim] = wall.coords 

703 cell_mids = self.facet_cell_midpoints(local_facets) 

704 if gdim == 2: 

705 a, b = local_faces[:, 0], local_faces[:, 1]