Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "geodata.hpp"
5 : #include "geodata_impl.hpp"
6 :
7 : #include <algorithm>
8 : #include <array>
9 : #include <limits>
10 : #include <numeric>
11 : #include <queue>
12 : #include <sstream>
13 : #include <string>
14 : #include <unordered_map>
15 : #include <unordered_set>
16 : #include <utility>
17 : #include <Eigen/Dense>
18 : #include <fmt/ranges.h>
19 : #include "fem/interpolator.hpp"
20 : #include "utils/communication.hpp"
21 : #include "utils/diagnostic.hpp"
22 : #include "utils/filesystem.hpp"
23 : #include "utils/meshio.hpp"
24 : #include "utils/omp.hpp"
25 : #include "utils/prettyprint.hpp"
26 : #include "utils/timer.hpp"
27 :
28 : namespace palace
29 : {
30 :
31 : using Vector3dMap = Eigen::Map<Eigen::Vector3d>;
32 : using CVector3dMap = Eigen::Map<const Eigen::Vector3d>;
33 :
34 : namespace
35 : {
36 :
37 : // Floating point precision for mesh IO. This precision is important, make sure nothing is
38 : // lost!
39 : constexpr auto MSH_FLT_PRECISION = std::numeric_limits<double>::max_digits10;
40 :
41 : // Load the serial mesh from disk.
42 : std::unique_ptr<mfem::Mesh> LoadMesh(const std::string &, bool,
43 : const config::BoundaryData &);
44 :
45 : // Clean the provided serial mesh by removing unused domain and boundary elements.
46 : void CleanMesh(std::unique_ptr<mfem::Mesh> &, const std::vector<int> &);
47 :
48 : // Create a new mesh by splitting all elements of the mesh into simplices or hexes
49 : // (using tet-to-hex). Optionally preserves curvature of the original mesh by interpolating
50 : // the high-order nodes with GSLIB.
51 : void SplitMeshElements(std::unique_ptr<mfem::Mesh> &, bool, bool);
52 :
53 : // Optionally reorder mesh elements based on MFEM's internal reordering tools for improved
54 : // cache usage.
55 : void ReorderMeshElements(mfem::Mesh &, bool = true);
56 :
57 : // Check that mesh boundary conditions are given for external boundaries.
58 : std::unordered_map<int, int> CheckMesh(const mfem::Mesh &, const config::BoundaryData &);
59 :
60 : // Adding boundary elements for material interfaces and exterior boundaries, and "crack"
61 : // desired internal boundary elements to disconnect the elements on either side.
62 : int AddInterfaceBdrElements(const IoData &, std::unique_ptr<mfem::Mesh> &,
63 : std::unordered_map<int, int> &, MPI_Comm comm);
64 :
65 : // Generate element-based mesh partitioning, using either a provided file or METIS.
66 : std::unique_ptr<int[]> GetMeshPartitioning(const mfem::Mesh &, int,
67 : const std::string & = "", bool = true);
68 :
69 : // Given a serial mesh on the root processor and element partitioning, create a parallel
70 : // mesh over the given communicator. The serial mesh is destroyed when no longer needed.
71 : std::unique_ptr<mfem::ParMesh> DistributeMesh(MPI_Comm, std::unique_ptr<mfem::Mesh> &,
72 : const int *, const std::string & = "");
73 :
74 : // Rebalance a conformal mesh across processor ranks, using the MeshPartitioner. Gathers the
75 : // mesh onto the root rank before scattering the partitioned mesh.
76 : void RebalanceConformalMesh(std::unique_ptr<mfem::ParMesh> &);
77 :
78 : } // namespace
79 :
80 : namespace mesh
81 : {
82 :
83 8 : std::unique_ptr<mfem::ParMesh> ReadMesh(const IoData &iodata, MPI_Comm comm)
84 : {
85 : // If possible on root, read the serial mesh (converting format if necessary), and do all
86 : // necessary serial preprocessing. When finished, distribute the mesh to all processes.
87 : // Count disk I/O time separately for the mesh read from file.
88 8 : BlockTimer bt0(Timer::MESH_PREPROCESS);
89 :
90 : // If not doing any local adaptation, or performing conformal adaptation, we can use the
91 : // mesh partitioner.
92 8 : std::unique_ptr<mfem::Mesh> smesh;
93 8 : const auto &refinement = iodata.model.refinement;
94 8 : const bool use_amr = (refinement.max_it > 0) || [&refinement]()
95 : {
96 8 : for (const auto &box : refinement.GetBoxes())
97 : {
98 0 : if (box.ref_levels > 0)
99 : {
100 : return true;
101 : }
102 : }
103 8 : for (const auto &sphere : refinement.GetSpheres())
104 : {
105 0 : if (sphere.ref_levels > 0)
106 : {
107 : return true;
108 : }
109 : }
110 : return false;
111 8 : }();
112 :
113 0 : const bool use_mesh_partitioner = [&]()
114 : {
115 : // Root must load the mesh to discover if nonconformal, as a previously adapted mesh
116 : // might be reused for nonadaptive simulations.
117 8 : BlockTimer bt(Timer::IO);
118 8 : bool use_mesh_partitioner = !use_amr || !refinement.nonconformal;
119 8 : if (Mpi::Root(comm))
120 : {
121 8 : smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature, iodata.boundaries);
122 8 : use_mesh_partitioner &= smesh->Conforming(); // The initial mesh must be conformal
123 : }
124 8 : Mpi::Broadcast(1, &use_mesh_partitioner, 0, comm);
125 8 : return use_mesh_partitioner;
126 16 : }();
127 :
128 : MPI_Comm node_comm;
129 8 : if (!use_mesh_partitioner)
130 : {
131 0 : MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, Mpi::Rank(comm), MPI_INFO_NULL,
132 : &node_comm);
133 : }
134 :
135 : {
136 8 : BlockTimer bt1(Timer::IO);
137 8 : if (!use_mesh_partitioner && Mpi::Root(node_comm) && !Mpi::Root(comm))
138 : {
139 : // Only one process per node reads the serial mesh, if not using mesh partitioner.
140 0 : smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature, iodata.boundaries);
141 : MFEM_VERIFY(!(smesh->Nonconforming() && use_mesh_partitioner),
142 : "Cannot use mesh partitioner on a nonconforming mesh!");
143 : }
144 8 : Mpi::Barrier(comm);
145 8 : }
146 :
147 : // Do some mesh preprocessing, and generate the partitioning.
148 : std::unique_ptr<int[]> partitioning;
149 8 : if (smesh)
150 : {
151 : // Check the the AMR specification and the mesh elements are compatible.
152 8 : const auto element_types = CheckElements(*smesh);
153 8 : MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_hexahedra ||
154 : refinement.nonconformal,
155 : "If there are tensor elements, AMR must be nonconformal!");
156 8 : MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_prisms ||
157 : refinement.nonconformal,
158 : "If there are wedge elements, AMR must be nonconformal!");
159 8 : MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_pyramids ||
160 : refinement.nonconformal,
161 : "If there are pyramid elements, AMR must be nonconformal!");
162 8 : MFEM_VERIFY(
163 : smesh->Conforming() || !use_amr || refinement.nonconformal,
164 : "The provided mesh is nonconformal, only nonconformal AMR can be performed!");
165 :
166 : // Clean up unused domain elements from the mesh.
167 8 : if (iodata.model.clean_unused_elements)
168 : {
169 : std::vector<int> attr_list;
170 8 : std::merge(iodata.domains.attributes.begin(), iodata.domains.attributes.end(),
171 : iodata.domains.postpro.attributes.begin(),
172 : iodata.domains.postpro.attributes.end(), std::back_inserter(attr_list));
173 8 : attr_list.erase(std::unique(attr_list.begin(), attr_list.end()), attr_list.end());
174 8 : CleanMesh(smesh, attr_list);
175 : }
176 :
177 : // Optionally convert mesh elements to simplices, for example in order to enable
178 : // conformal mesh refinement, or hexes.
179 8 : if (iodata.model.make_simplex || iodata.model.make_hex)
180 : {
181 0 : SplitMeshElements(smesh, iodata.model.make_simplex, iodata.model.make_hex);
182 : }
183 :
184 : // Optionally reorder elements (and vertices) based on spatial location after loading
185 : // the serial mesh.
186 8 : if (iodata.model.reorder_elements)
187 : {
188 0 : ReorderMeshElements(*smesh);
189 : }
190 :
191 : // Refine the serial mesh (not typically used, prefer parallel uniform refinement
192 : // instead).
193 : {
194 8 : int ne = smesh->GetNE();
195 8 : for (int l = 0; l < iodata.model.refinement.ser_uniform_ref_levels; l++)
196 : {
197 0 : smesh->UniformRefinement();
198 : }
199 8 : if (iodata.model.refinement.ser_uniform_ref_levels > 0)
200 : {
201 0 : Mpi::Print("Serial uniform mesh refinement levels added {:d} elements (initial = "
202 : "{:d}, final = {:d})\n",
203 0 : smesh->GetNE() - ne, ne, smesh->GetNE());
204 : }
205 : }
206 :
207 : // Check the final mesh, throwing warnings if there are exterior boundaries with no
208 : // associated boundary condition.
209 8 : if (smesh->Conforming())
210 : {
211 8 : auto face_to_be = CheckMesh(*smesh, iodata.boundaries);
212 :
213 : // Add new boundary elements for material interfaces if not present (with new unique
214 : // boundary attributes). Also duplicate internal boundary elements associated with
215 : // cracks if desired.
216 8 : if ((iodata.model.crack_bdr_elements || iodata.model.add_bdr_elements))
217 : {
218 : // Split all internal (non periodic) boundary elements for boundary attributes where
219 : // BC are applied (not just postprocessing).
220 8 : while (AddInterfaceBdrElements(iodata, smesh, face_to_be, comm) != 1)
221 : {
222 : // May require multiple calls due to early exit/retry approach.
223 : }
224 : }
225 : }
226 : else
227 : {
228 0 : Mpi::Warning("{} is a nonconformal mesh, assumed from previous AMR iteration.\n"
229 : "Skipping mesh modification preprocessing steps!\n\n",
230 0 : iodata.model.mesh);
231 : }
232 :
233 : // Finally, finalize the serial mesh. Mark tetrahedral meshes for refinement. There
234 : // should be no need to fix orientation as this was done during initial mesh loading
235 : // from disk.
236 : constexpr bool refine = true, fix_orientation = false;
237 8 : smesh->Finalize(refine, fix_orientation);
238 :
239 : // Generate the mesh partitioning.
240 8 : partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), iodata.model.partitioning);
241 : }
242 :
243 : // Distribute the mesh.
244 8 : std::unique_ptr<mfem::ParMesh> pmesh;
245 8 : if (use_mesh_partitioner)
246 : {
247 16 : pmesh = DistributeMesh(comm, smesh, partitioning.get(), iodata.problem.output);
248 : }
249 : else
250 : {
251 : // Send the preprocessed serial mesh and partitioning as a byte string.
252 0 : constexpr bool generate_edges = false, refine = true, fix_orientation = false;
253 : std::string so;
254 0 : int slen = 0;
255 0 : if (smesh)
256 : {
257 0 : std::ostringstream fo(std::stringstream::out);
258 : // fo << std::fixed;
259 : fo << std::scientific;
260 : fo.precision(MSH_FLT_PRECISION);
261 0 : smesh->Print(fo);
262 : smesh.reset(); // Root process needs to rebuild the mesh to ensure consistency with
263 : // the saved serial mesh (refinement marking, for example)
264 0 : so = fo.str();
265 : // so = zlib::CompressString(fo.str());
266 0 : slen = static_cast<int>(so.size());
267 0 : MFEM_VERIFY(so.size() == (std::size_t)slen, "Overflow in stringbuffer size!");
268 0 : }
269 0 : Mpi::Broadcast(1, &slen, 0, node_comm);
270 0 : if (so.empty())
271 : {
272 0 : so.resize(slen);
273 : }
274 0 : Mpi::Broadcast(slen, so.data(), 0, node_comm);
275 : {
276 0 : std::istringstream fi(so);
277 : // std::istringstream fi(zlib::DecompressString(so));
278 0 : smesh = std::make_unique<mfem::Mesh>(fi, generate_edges, refine, fix_orientation);
279 : so.clear();
280 0 : }
281 0 : if (refinement.nonconformal && use_amr)
282 : {
283 0 : smesh->EnsureNCMesh(true);
284 : }
285 0 : if (!partitioning)
286 : {
287 0 : partitioning = std::make_unique<int[]>(smesh->GetNE());
288 : }
289 0 : Mpi::Broadcast(smesh->GetNE(), partitioning.get(), 0, node_comm);
290 0 : MPI_Comm_free(&node_comm);
291 0 : pmesh = std::make_unique<mfem::ParMesh>(comm, *smesh, partitioning.get());
292 : smesh.reset();
293 : }
294 :
295 : if constexpr (false)
296 : {
297 : auto tmp = fs::path(iodata.problem.output) / "tmp";
298 : if (Mpi::Root(comm) && !fs::exists(tmp))
299 : {
300 : fs::create_directories(tmp);
301 : }
302 : int width = 1 + static_cast<int>(std::log10(Mpi::Size(comm) - 1));
303 : std::unique_ptr<mfem::Mesh> gsmesh =
304 : LoadMesh(iodata.model.mesh, iodata.model.remove_curvature, iodata.boundaries);
305 : std::unique_ptr<int[]> gpartitioning = GetMeshPartitioning(*gsmesh, Mpi::Size(comm));
306 : mfem::ParMesh gpmesh(comm, *gsmesh, gpartitioning.get(), 0);
307 : {
308 : std::string pfile =
309 : mfem::MakeParFilename(tmp.string() + "part.", Mpi::Rank(comm), ".mesh", width);
310 : std::ofstream fo(pfile);
311 : // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available
312 : fo.precision(MSH_FLT_PRECISION);
313 : gpmesh.ParPrint(fo);
314 : }
315 : {
316 : std::string pfile =
317 : mfem::MakeParFilename(tmp.string() + "final.", Mpi::Rank(comm), ".mesh", width);
318 : std::ofstream fo(pfile);
319 : // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available
320 : fo.precision(MSH_FLT_PRECISION);
321 : pmesh->ParPrint(fo);
322 : }
323 : }
324 :
325 8 : return pmesh;
326 8 : }
327 :
328 8 : void RefineMesh(const IoData &iodata, std::vector<std::unique_ptr<mfem::ParMesh>> &mesh)
329 : {
330 : // Prepare for uniform and region-based refinement.
331 8 : MFEM_VERIFY(mesh.size() == 1,
332 : "Input mesh vector before refinement has more than a single mesh!");
333 8 : int uniform_ref_levels = iodata.model.refinement.uniform_ref_levels;
334 : int max_region_ref_levels = 0;
335 8 : for (const auto &box : iodata.model.refinement.GetBoxes())
336 : {
337 0 : if (max_region_ref_levels < box.ref_levels)
338 : {
339 : max_region_ref_levels = box.ref_levels;
340 : }
341 : }
342 8 : for (const auto &sphere : iodata.model.refinement.GetSpheres())
343 : {
344 0 : if (max_region_ref_levels < sphere.ref_levels)
345 : {
346 : max_region_ref_levels = sphere.ref_levels;
347 : }
348 : }
349 8 : if (iodata.solver.linear.mg_use_mesh && iodata.solver.linear.mg_max_levels > 1)
350 : {
351 8 : mesh.reserve(1 + uniform_ref_levels + max_region_ref_levels);
352 : }
353 :
354 : // Prior to MFEM's PR #1046, the tetrahedral mesh required reorientation after all mesh
355 : // refinement in order to define higher-order Nedelec spaces on it. This is technically
356 : // not required after MFEM's PR #1046, but in case you want to be absolutely sure, we
357 : // reorient only the coarse mesh so that the refinements are still true refinements of
358 : // the original mesh (required for geometric multigrid). Otherwise, it happens after
359 : // refinement.
360 8 : if (iodata.model.reorient_tet_mesh && mesh.capacity() > 1)
361 : {
362 : PalacePragmaDiagnosticPush
363 : PalacePragmaDiagnosticDisableDeprecated
364 0 : mesh[0]->ReorientTetMesh();
365 : PalacePragmaDiagnosticPop
366 : }
367 :
368 : // Uniformly refine the mesh further in parallel, saving the level meshes for geometric
369 : // coarsening later on if desired.
370 8 : for (int l = 0; l < uniform_ref_levels; l++)
371 : {
372 0 : if (mesh.capacity() > 1)
373 : {
374 0 : mesh.emplace_back(std::make_unique<mfem::ParMesh>(*mesh.back()));
375 : }
376 0 : mesh.back()->UniformRefinement();
377 : }
378 :
379 : // Simplex meshes need to be re-finalized in order to use local refinement (see
380 : // the docstring for mfem::Mesh::UniformRefinement).
381 8 : const auto element_types = mesh::CheckElements(*mesh.back());
382 8 : if (element_types.has_simplices && uniform_ref_levels > 0 &&
383 0 : (max_region_ref_levels > 0 || iodata.model.refinement.max_it > 0))
384 : {
385 : constexpr bool refine = true, fix_orientation = false;
386 0 : Mpi::Print("\nFlattening mesh sequence:\n Local mesh refinement will start from the "
387 : "final uniformly-refined mesh\n");
388 : mesh.erase(mesh.begin(), mesh.end() - 1);
389 0 : mesh.back()->Finalize(refine, fix_orientation);
390 : }
391 :
392 : // Proceed with region-based refinement, level-by-level for all regions. Currently support
393 : // box and sphere region shapes. Any overlap between regions is ignored (take the union,
394 : // don't double-refine).
395 8 : MFEM_VERIFY(
396 : max_region_ref_levels == 0 ||
397 : !(element_types.has_hexahedra || element_types.has_prisms ||
398 : element_types.has_pyramids) ||
399 : mesh.back()->Nonconforming(),
400 : "Region-based refinement for non-simplex meshes requires a nonconformal mesh!");
401 : const bool use_nodes = (mesh.back()->GetNodes() != nullptr);
402 8 : const int ref = use_nodes ? mesh.back()->GetNodes()->FESpace()->GetMaxElementOrder() : 1;
403 : const int dim = mesh.back()->SpaceDimension();
404 : int region_ref_level = 0;
405 8 : while (region_ref_level < max_region_ref_levels)
406 : {
407 : // Mark elements for refinement in all regions. An element is marked for refinement if
408 : // any of its vertices are inside any refinement region for the given level.
409 : mfem::Array<mfem::Refinement> refinements;
410 0 : for (int i = 0; i < mesh.back()->GetNE(); i++)
411 : {
412 : bool refine = false;
413 0 : mfem::DenseMatrix pointmat;
414 0 : if (use_nodes)
415 : {
416 0 : mfem::ElementTransformation &T = *mesh.back()->GetElementTransformation(i);
417 : mfem::RefinedGeometry *RefG =
418 0 : mfem::GlobGeometryRefiner.Refine(T.GetGeometryType(), ref);
419 0 : T.Transform(RefG->RefPts, pointmat);
420 : }
421 : else
422 : {
423 0 : const int *verts = mesh.back()->GetElement(i)->GetVertices();
424 0 : const int nv = mesh.back()->GetElement(i)->GetNVertices();
425 0 : pointmat.SetSize(dim, nv);
426 0 : for (int j = 0; j < nv; j++)
427 : {
428 0 : const double *coord = mesh.back()->GetVertex(verts[j]);
429 0 : for (int d = 0; d < dim; d++)
430 : {
431 0 : pointmat(d, j) = coord[d];
432 : }
433 : }
434 : }
435 0 : for (const auto &box : iodata.model.refinement.GetBoxes())
436 : {
437 0 : if (region_ref_level < box.ref_levels)
438 : {
439 0 : for (int j = 0; j < pointmat.Width(); j++)
440 : {
441 : // Check if the point is inside the box.
442 : int d = 0;
443 0 : for (; d < pointmat.Height(); d++)
444 : {
445 0 : if (pointmat(d, j) < box.bbmin[d] || pointmat(d, j) > box.bbmax[d])
446 : {
447 : break;
448 : }
449 : }
450 0 : if (d == dim)
451 : {
452 : refine = true;
453 : break;
454 : }
455 : }
456 0 : if (refine)
457 : {
458 : break;
459 : }
460 : }
461 : }
462 0 : if (refine)
463 : {
464 0 : refinements.Append(mfem::Refinement(i));
465 : continue;
466 : }
467 0 : for (const auto &sphere : iodata.model.refinement.GetSpheres())
468 : {
469 0 : if (region_ref_level < sphere.ref_levels)
470 : {
471 0 : for (int j = 0; j < pointmat.Width(); j++)
472 : {
473 : // Check if the point is inside the sphere.
474 : double dist = 0.0;
475 0 : for (int d = 0; d < pointmat.Height(); d++)
476 : {
477 0 : double s = pointmat(d, j) - sphere.center[d];
478 0 : dist += s * s;
479 : }
480 0 : if (dist <= sphere.r * sphere.r)
481 : {
482 : refine = true;
483 : break;
484 : }
485 : }
486 0 : if (refine)
487 : {
488 : break;
489 : }
490 : }
491 : }
492 0 : if (refine)
493 : {
494 0 : refinements.Append(mfem::Refinement(i));
495 : }
496 0 : }
497 :
498 : // Do the refinement. For tensor element meshes, this may make the mesh nonconforming
499 : // (adds hanging nodes).
500 0 : if (mesh.capacity() > 1)
501 : {
502 0 : mesh.emplace_back(std::make_unique<mfem::ParMesh>(*mesh.back()));
503 : }
504 0 : mesh.back()->GeneralRefinement(refinements, -1);
505 0 : region_ref_level++;
506 : }
507 8 : if (max_region_ref_levels > 0 && mesh.capacity() == 1)
508 : {
509 0 : RebalanceMesh(iodata, mesh[0]);
510 : }
511 :
512 : // Prior to MFEM's PR #1046, the tetrahedral mesh required reorientation after all mesh
513 : // refinement in order to define higher-order Nedelec spaces on it. This is technically
514 : // not required after MFEM's PR #1046, but in case you want to be absolutely sure, we
515 : // reorient only the mesh after refinement if there is a single mesh (doesn't work with
516 : // h-refinement geometric multigrid).
517 8 : if (iodata.model.reorient_tet_mesh && mesh.capacity() == 1)
518 : {
519 : PalacePragmaDiagnosticPush
520 : PalacePragmaDiagnosticDisableDeprecated
521 0 : mesh[0]->ReorientTetMesh();
522 : PalacePragmaDiagnosticPop
523 : }
524 :
525 : // Print some mesh information.
526 : mfem::Vector bbmin, bbmax;
527 8 : GetAxisAlignedBoundingBox(*mesh[0], bbmin, bbmax);
528 : const double Lc = iodata.units.Dimensionalize<Units::ValueType::LENGTH>(1.0);
529 8 : Mpi::Print(mesh[0]->GetComm(), "\nMesh curvature order: {}\nMesh bounding box:\n",
530 : mesh[0]->GetNodes()
531 8 : ? std::to_string(mesh[0]->GetNodes()->FESpace()->GetMaxElementOrder())
532 8 : : "None");
533 8 : if (mesh[0]->SpaceDimension() == 3)
534 : {
535 8 : Mpi::Print(mesh[0]->GetComm(),
536 : " (Xmin, Ymin, Zmin) = ({:+.3e}, {:+.3e}, {:+.3e}) m\n"
537 : " (Xmax, Ymax, Zmax) = ({:+.3e}, {:+.3e}, {:+.3e}) m\n",
538 16 : bbmin[0] * Lc, bbmin[1] * Lc, bbmin[2] * Lc, bbmax[0] * Lc, bbmax[1] * Lc,
539 8 : bbmax[2] * Lc);
540 : }
541 : else
542 : {
543 0 : Mpi::Print(mesh[0]->GetComm(),
544 : " (Xmin, Ymin) = ({:+.3e}, {:+.3e}) m\n"
545 : " (Xmax, Ymax) = ({:+.3e}, {:+.3e}) m\n",
546 0 : bbmin[0] * Lc, bbmin[1] * Lc, bbmax[0] * Lc, bbmax[1] * Lc);
547 : }
548 16 : Mpi::Print(mesh[0]->GetComm(), "\n{}", (mesh.size() > 1) ? "Coarse " : "");
549 8 : mesh[0]->PrintInfo();
550 8 : if (mesh.size() > 1)
551 : {
552 0 : Mpi::Print(mesh[0]->GetComm(), "\nRefined ");
553 0 : mesh.back()->PrintInfo();
554 : }
555 8 : }
556 :
557 1 : mfem::Mesh MeshTetToHex(const mfem::Mesh &orig_mesh)
558 : {
559 : // Courtesy of https://gist.github.com/pazner/e9376f77055c0918d7c43e034e9e5888, only
560 : // supports tetrahedral elements for now. Eventually should be expanded to support prism
561 : // and pyramid elements but this mixed mesh support requires a bit more work.
562 1 : MFEM_VERIFY(orig_mesh.Dimension() == 3, "Tet-to-hex conversion only supports 3D meshes!");
563 : {
564 : // This checks the local mesh on each process, but the assertion failing on any single
565 : // process will terminate the program.
566 : mfem::Array<mfem::Geometry::Type> geoms;
567 1 : orig_mesh.GetGeometries(3, geoms);
568 1 : MFEM_VERIFY(geoms.Size() == 1 && geoms[0] == mfem::Geometry::TETRAHEDRON,
569 : "Tet-to-hex conversion only works for pure tetrahedral meshes!");
570 : }
571 :
572 : // Add new vertices in every edge, face, and volume. Each tet is subdivided into 4 hexes,
573 : // and each triangular face subdivided into 3 quads.
574 : const int nv_tet = orig_mesh.GetNV();
575 : const int nedge_tet = orig_mesh.GetNEdges();
576 : const int nface_tet = orig_mesh.GetNFaces();
577 : const int ne_tet = orig_mesh.GetNE();
578 : const int nbe_tet = orig_mesh.GetNBE();
579 1 : const int nv = nv_tet + nedge_tet + nface_tet + ne_tet;
580 1 : const int ne = 4 * ne_tet; // 4 hex per tet
581 1 : const int nbe = 3 * nbe_tet; // 3 square per tri
582 1 : mfem::Mesh hex_mesh(orig_mesh.Dimension(), nv, ne, nbe, orig_mesh.SpaceDimension());
583 :
584 : // Add original vertices.
585 5 : for (int v = 0; v < nv_tet; v++)
586 : {
587 4 : hex_mesh.AddVertex(orig_mesh.GetVertex(v));
588 : }
589 :
590 : // Add midpoints of edges, faces, and elements.
591 11 : auto AddCentroid = [&orig_mesh, &hex_mesh](const int *verts, int nv)
592 : {
593 11 : double coord[3] = {0.0, 0.0, 0.0};
594 39 : for (int i = 0; i < nv; i++)
595 : {
596 112 : for (int d = 0; d < orig_mesh.SpaceDimension(); d++)
597 : {
598 84 : coord[d] += orig_mesh.GetVertex(verts[i])[d] / nv;
599 : }
600 : }
601 11 : hex_mesh.AddVertex(coord);
602 1 : };
603 : {
604 : mfem::Array<int> verts;
605 7 : for (int e = 0; e < nedge_tet; ++e)
606 : {
607 6 : orig_mesh.GetEdgeVertices(e, verts);
608 6 : AddCentroid(verts.GetData(), verts.Size());
609 : }
610 : }
611 5 : for (int f = 0; f < nface_tet; ++f)
612 : {
613 8 : AddCentroid(orig_mesh.GetFace(f)->GetVertices(), orig_mesh.GetFace(f)->GetNVertices());
614 : }
615 2 : for (int e = 0; e < ne_tet; ++e)
616 : {
617 1 : AddCentroid(orig_mesh.GetElement(e)->GetVertices(),
618 1 : orig_mesh.GetElement(e)->GetNVertices());
619 : }
620 :
621 : // Connectivity of tetrahedron vertices to the edges. The vertices of the new mesh are
622 : // ordered so that the original tet vertices are first, then the vertices splitting each
623 : // edge, then the vertices at the center of each triangle face, then the center of the
624 : // tet. Thus the edge/face/element numbers can be used to index into the new array of
625 : // vertices, and the element local edge/face can be used to extract the global edge/face
626 : // index, and thus the corresponding vertex.
627 1 : constexpr int tet_vertex_edge_map[4 * 3] = {0, 1, 2, 3, 0, 4, 1, 3, 5, 5, 4, 2};
628 1 : constexpr int tet_vertex_face_map[4 * 3] = {3, 2, 1, 3, 0, 2, 3, 1, 0, 0, 1, 2};
629 1 : constexpr int tri_vertex_edge_map[3 * 2] = {0, 2, 1, 0, 2, 1};
630 :
631 : // Add four hexahedra for each tetrahedron.
632 : {
633 : mfem::Array<int> edges, faces, orients;
634 2 : for (int e = 0; e < ne_tet; ++e)
635 : {
636 : const int *verts = orig_mesh.GetElement(e)->GetVertices();
637 1 : orig_mesh.GetElementEdges(e, edges, orients);
638 1 : orig_mesh.GetElementFaces(e, faces, orients);
639 :
640 : // One hex for each vertex of the tet.
641 5 : for (int i = 0; i < 4; ++i)
642 : {
643 : int hex_v[8];
644 4 : hex_v[0] = verts[i];
645 4 : hex_v[1] = nv_tet + edges[tet_vertex_edge_map[3 * i + 0]];
646 4 : hex_v[2] = nv_tet + nedge_tet + faces[tet_vertex_face_map[3 * i + 0]];
647 4 : hex_v[3] = nv_tet + edges[tet_vertex_edge_map[3 * i + 1]];
648 4 : hex_v[4] = nv_tet + edges[tet_vertex_edge_map[3 * i + 2]];
649 4 : hex_v[5] = nv_tet + nedge_tet + faces[tet_vertex_face_map[3 * i + 1]];
650 4 : hex_v[6] = nv_tet + nedge_tet + nface_tet + e;
651 4 : hex_v[7] = nv_tet + nedge_tet + faces[tet_vertex_face_map[3 * i + 2]];
652 4 : hex_mesh.AddHex(hex_v, orig_mesh.GetAttribute(e));
653 : }
654 : }
655 : }
656 :
657 : // Add the boundary elements.
658 : {
659 : mfem::Array<int> edges, orients;
660 5 : for (int be = 0; be < nbe_tet; ++be)
661 : {
662 : int f, o;
663 : const int *verts = orig_mesh.GetBdrElement(be)->GetVertices();
664 4 : orig_mesh.GetBdrElementEdges(be, edges, orients);
665 4 : orig_mesh.GetBdrElementFace(be, &f, &o);
666 :
667 : // One quad for each vertex of the tri.
668 16 : for (int i = 0; i < 3; ++i)
669 : {
670 : int quad_v[4];
671 12 : quad_v[0] = verts[i];
672 12 : quad_v[1] = nv_tet + edges[tri_vertex_edge_map[2 * i + 0]];
673 12 : quad_v[2] = nv_tet + nedge_tet + f;
674 12 : quad_v[3] = nv_tet + edges[tri_vertex_edge_map[2 * i + 1]];
675 12 : hex_mesh.AddBdrQuad(quad_v, orig_mesh.GetBdrAttribute(be));
676 : }
677 : }
678 : }
679 :
680 : // Finalize the hex mesh topology. The mesh will be marked for refinement later on.
681 : constexpr bool generate_bdr = false;
682 1 : hex_mesh.FinalizeTopology(generate_bdr);
683 :
684 : // All elements have now been added, can construct the higher order field.
685 1 : if (orig_mesh.GetNodes())
686 : {
687 1 : hex_mesh.EnsureNodes();
688 : // Higher order associated to vertices are unchanged, and those for
689 : // previously existing edges. DOFs associated to new elements need to be set.
690 : const int sdim = orig_mesh.SpaceDimension();
691 : auto *orig_fespace = orig_mesh.GetNodes()->FESpace();
692 2 : hex_mesh.SetCurvature(orig_fespace->GetMaxElementOrder(), orig_fespace->IsDGSpace(),
693 : orig_mesh.SpaceDimension(), orig_fespace->GetOrdering());
694 :
695 : // Need to convert the hexahedra local coordinate system into the parent tetrahedra
696 : // system. Each hexahedra spans a different set of the tet's reference coordinates. To
697 : // convert between, define the reference coordinate locations of each of the vertices
698 : // the hexahedra will use, then perform trilinear interpolation in the reference space.
699 :
700 1 : auto [vert_loc, edge_loc, face_loc] = []()
701 : {
702 : std::array<std::array<double, 3>, 4> vert_loc{};
703 : vert_loc[0] = {0.0, 0.0, 0.0};
704 : vert_loc[1] = {1.0, 0.0, 0.0};
705 : vert_loc[2] = {0.0, 1.0, 0.0};
706 : vert_loc[3] = {0.0, 0.0, 1.0};
707 : std::array<std::array<double, 3>, 6> edge_loc{};
708 : edge_loc[0] = {0.5, 0.0, 0.0};
709 : edge_loc[1] = {0.0, 0.5, 0.0};
710 : edge_loc[2] = {0.0, 0.0, 0.5};
711 : edge_loc[3] = {0.5, 0.5, 0.0};
712 : edge_loc[4] = {0.5, 0.0, 0.5};
713 : edge_loc[5] = {0.0, 0.5, 0.5};
714 : std::array<std::array<double, 3>, 6> face_loc{};
715 : face_loc[0] = {1.0 / 3, 1.0 / 3, 1.0 / 3};
716 : face_loc[1] = {0.0, 1.0 / 3, 1.0 / 3};
717 : face_loc[2] = {1.0 / 3, 0.0, 1.0 / 3};
718 : face_loc[3] = {1.0 / 3, 1.0 / 3, 0.0};
719 : return std::make_tuple(vert_loc, edge_loc, face_loc);
720 : }();
721 : std::array<double, 3> centroid{{0.25, 0.25, 0.25}};
722 :
723 : // We assume the Nodes field is of a single order, and there is a single tet originally.
724 : // The nodes within the reference hex and parent tet are always the same, so we use the
725 : // typical FE. We then exploit the fact the map between reference spaces is always
726 : // linear, and construct the transformation explicitly.
727 1 : const auto *orig_FE = orig_mesh.GetNodes()->FESpace()->GetTypicalFE();
728 1 : const auto *child_FE = hex_mesh.GetNodes()->FESpace()->GetTypicalFE();
729 : // Original shape function (i), at new element nodes (j), for each new element (k).
730 1 : mfem::DenseTensor shape(orig_FE->GetDof(), child_FE->GetDof(), 4);
731 : mfem::Vector col; // For slicing into matrices within shape
732 5 : for (int i = 0; i < 4; i++)
733 : {
734 : // Collect the vertices of the new hex within the tet.
735 : std::array<std::array<double, 3>, 8> hex_verts;
736 4 : hex_verts[0] = vert_loc[i];
737 4 : hex_verts[1] = edge_loc[tet_vertex_edge_map[3 * i + 0]];
738 4 : hex_verts[2] = face_loc[tet_vertex_face_map[3 * i + 0]];
739 4 : hex_verts[3] = edge_loc[tet_vertex_edge_map[3 * i + 1]];
740 4 : hex_verts[4] = edge_loc[tet_vertex_edge_map[3 * i + 2]];
741 4 : hex_verts[5] = face_loc[tet_vertex_face_map[3 * i + 1]];
742 4 : hex_verts[6] = centroid;
743 4 : hex_verts[7] = face_loc[tet_vertex_face_map[3 * i + 2]];
744 36 : for (int j = 0; j < child_FE->GetNodes().Size(); j++)
745 : {
746 : const auto &cn = child_FE->GetNodes()[j];
747 : mfem::IntegrationPoint cn_in_orig;
748 :
749 : // Perform trilinear interpolation from (u,v,w) the unit ref coords in the new hex,
750 : // and the corresponding nodes in the containing tet.
751 : // clang-format off
752 : // x component
753 32 : cn_in_orig.x =
754 32 : hex_verts[0][0] * (1-cn.x) * (1-cn.y) * (1-cn.z) +
755 32 : hex_verts[1][0] * cn.x * (1-cn.y) * (1-cn.z) +
756 32 : hex_verts[2][0] * cn.x * cn.y * (1-cn.z) +
757 32 : hex_verts[3][0] * (1-cn.x) * cn.y * (1-cn.z) +
758 32 : hex_verts[4][0] * (1-cn.x) * (1-cn.y) * cn.z +
759 32 : hex_verts[5][0] * cn.x * (1-cn.y) * cn.z +
760 32 : hex_verts[6][0] * cn.x * cn.y * cn.z +
761 32 : hex_verts[7][0] * (1-cn.x) * cn.y * cn.z;
762 :
763 : // y component
764 32 : cn_in_orig.y =
765 32 : hex_verts[0][1] * (1-cn.x) * (1-cn.y) * (1-cn.z) +
766 32 : hex_verts[1][1] * cn.x * (1-cn.y) * (1-cn.z) +
767 32 : hex_verts[2][1] * cn.x * cn.y * (1-cn.z) +
768 32 : hex_verts[3][1] * (1-cn.x) * cn.y * (1-cn.z) +
769 32 : hex_verts[4][1] * (1-cn.x) * (1-cn.y) * cn.z +
770 32 : hex_verts[5][1] * cn.x * (1-cn.y) * cn.z +
771 32 : hex_verts[6][1] * cn.x * cn.y * cn.z +
772 32 : hex_verts[7][1] * (1-cn.x) * cn.y * cn.z;
773 :
774 : // z component
775 32 : cn_in_orig.z =
776 32 : hex_verts[0][2] * (1-cn.x) * (1-cn.y) * (1-cn.z) +
777 32 : hex_verts[1][2] * cn.x * (1-cn.y) * (1-cn.z) +
778 32 : hex_verts[2][2] * cn.x * cn.y * (1-cn.z) +
779 32 : hex_verts[3][2] * (1-cn.x) * cn.y * (1-cn.z) +
780 32 : hex_verts[4][2] * (1-cn.x) * (1-cn.y) * cn.z +
781 32 : hex_verts[5][2] * cn.x * (1-cn.y) * cn.z +
782 32 : hex_verts[6][2] * cn.x * cn.y * cn.z +
783 32 : hex_verts[7][2] * (1-cn.x) * cn.y * cn.z;
784 : // clang-format on
785 : shape(i).GetColumnReference(j, col);
786 32 : orig_FE->CalcShape(cn_in_orig, col);
787 : }
788 : }
789 :
790 : // Each submatrix of shape tensor now encodes the reference coordinates of each hex
791 : // within the containing tet. Extracting the specific element dof values, and applying
792 : // to the correct shape slice will now give the requisite higher order dofs evaluated at
793 : // the refined elements nodes.
794 : mfem::Array<int> hex_dofs;
795 1 : mfem::DenseMatrix point_matrix(child_FE->GetDof(), sdim); // nnode_child x sdim
796 1 : mfem::Vector dof_vals(orig_FE->GetDof() * sdim);
797 : mfem::DenseMatrix dof_vals_mat(dof_vals.GetData(), orig_FE->GetDof(), sdim);
798 2 : for (int e = 0; e < ne_tet; ++e)
799 : {
800 : // Returns byNODES no matter what, because FiniteElementSpace::GetElementVDofs does.
801 : // Matches the GetElementVDofs call below, which similarly always uses byNODES.
802 1 : orig_mesh.GetNodes()->GetElementDofValues(e, dof_vals);
803 5 : for (int i = 0; i < 4; i++)
804 : {
805 : // shape(i) : orig_FE->GetDof() x hex_FE->GetDof()
806 : // dof_vals_mat : orig_FE->GetDof() x sdim
807 : // point_matrix : child_FE->GetDof() x sdim
808 4 : MultAtB(shape(i), dof_vals_mat, point_matrix);
809 4 : hex_mesh.GetNodes()->FESpace()->GetElementVDofs(4 * e + i, hex_dofs);
810 4 : hex_mesh.GetNodes()->SetSubVector(hex_dofs, point_matrix.GetData());
811 : }
812 : }
813 2 : }
814 :
815 1 : return hex_mesh;
816 0 : }
817 :
818 : namespace
819 : {
820 :
821 86 : void ScaleMesh(mfem::Mesh &mesh, double L)
822 : {
823 86 : PalacePragmaOmp(parallel for schedule(static))
824 : for (int i = 0; i < mesh.GetNV(); i++)
825 : {
826 : double *v = mesh.GetVertex(i);
827 : std::transform(v, v + mesh.SpaceDimension(), v, [L](double val) { return val * L; });
828 : }
829 86 : if (auto *pmesh = dynamic_cast<mfem::ParMesh *>(&mesh))
830 : {
831 86 : PalacePragmaOmp(parallel for schedule(static))
832 : for (int i = 0; i < pmesh->face_nbr_vertices.Size(); i++)
833 : {
834 : double *v = pmesh->face_nbr_vertices[i]();
835 : std::transform(v, v + mesh.SpaceDimension(), v, [L](double val) { return val * L; });
836 : }
837 : }
838 86 : if (mesh.GetNodes())
839 : {
840 68 : *mesh.GetNodes() *= L;
841 68 : if (auto *pnodes = dynamic_cast<mfem::ParGridFunction *>(mesh.GetNodes()))
842 : {
843 68 : pnodes->FaceNbrData() *= L;
844 : }
845 : }
846 86 : }
847 :
848 : } // namespace
849 :
850 30 : void DimensionalizeMesh(mfem::Mesh &mesh, double L)
851 : {
852 30 : ScaleMesh(mesh, L);
853 30 : }
854 :
855 56 : void NondimensionalizeMesh(mfem::Mesh &mesh, double L)
856 : {
857 56 : ScaleMesh(mesh, 1.0 / L);
858 56 : }
859 :
860 0 : std::vector<mfem::Geometry::Type> ElementTypeInfo::GetGeomTypes() const
861 : {
862 : std::vector<mfem::Geometry::Type> geom_types;
863 0 : if (has_simplices)
864 : {
865 0 : geom_types.push_back(mfem::Geometry::TETRAHEDRON);
866 : }
867 0 : if (has_hexahedra)
868 : {
869 0 : geom_types.push_back(mfem::Geometry::CUBE);
870 : }
871 0 : if (has_prisms)
872 : {
873 0 : geom_types.push_back(mfem::Geometry::PRISM);
874 : }
875 0 : if (has_pyramids)
876 : {
877 0 : geom_types.push_back(mfem::Geometry::PYRAMID);
878 : }
879 0 : return geom_types;
880 : }
881 :
882 16 : ElementTypeInfo CheckElements(const mfem::Mesh &mesh)
883 : {
884 : // MeshGenerator is reduced over the communicator. This checks for geometries on any
885 : // processor.
886 : auto meshgen = mesh.MeshGenerator();
887 16 : return {bool(meshgen & 1), bool(meshgen & 2), bool(meshgen & 4), bool(meshgen & 8)};
888 : }
889 :
890 0 : bool CheckRefinementFlags(const mfem::Mesh &mesh)
891 : {
892 0 : bool marked = true;
893 0 : for (int e = 0; e < mesh.GetNE(); e++)
894 : {
895 : const mfem::Element *el = mesh.GetElement(e);
896 : const int geom = el->GetGeometryType();
897 0 : if (geom == mfem::Geometry::TETRAHEDRON)
898 : {
899 : const mfem::Tetrahedron *tet = static_cast<const mfem::Tetrahedron *>(el);
900 0 : if (tet->GetRefinementFlag() == 0)
901 : {
902 0 : marked = false;
903 : break;
904 : }
905 : }
906 : }
907 0 : if (const auto *pmesh = dynamic_cast<const mfem::ParMesh *>(&mesh))
908 : {
909 : Mpi::GlobalAnd(1, &marked, pmesh->GetComm());
910 : }
911 0 : return marked;
912 : }
913 :
914 317 : void AttrToMarker(int max_attr, const int *attr_list, int attr_list_size,
915 : mfem::Array<int> &marker, bool skip_invalid)
916 : {
917 457 : MFEM_VERIFY(skip_invalid || attr_list_size == 0 ||
918 : *std::max_element(attr_list, attr_list + attr_list_size) <= max_attr,
919 : "Invalid attribute number present ("
920 : << *std::max_element(attr_list, attr_list + attr_list_size) << ")!");
921 : marker.SetSize(max_attr);
922 317 : if (attr_list_size == 1 && attr_list[0] == -1)
923 : {
924 : marker = 1;
925 : }
926 : else
927 : {
928 : marker = 0;
929 725 : for (int i = 0; i < attr_list_size; i++)
930 : {
931 408 : int attr = attr_list[i];
932 408 : if ((attr <= 0 || attr > max_attr) && skip_invalid)
933 : {
934 0 : continue;
935 : }
936 0 : MFEM_VERIFY(attr > 0, "Attribute number less than one!");
937 408 : MFEM_VERIFY(marker[attr - 1] == 0, "Repeated attribute in attribute list!");
938 408 : marker[attr - 1] = 1;
939 : }
940 : }
941 317 : }
942 :
943 73 : void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
944 : bool bdr, mfem::Vector &min, mfem::Vector &max)
945 : {
946 73 : int dim = mesh.SpaceDimension();
947 73 : min.SetSize(dim);
948 73 : max.SetSize(dim);
949 292 : for (int d = 0; d < dim; d++)
950 : {
951 219 : min(d) = mfem::infinity();
952 219 : max(d) = -mfem::infinity();
953 : }
954 73 : if (!mesh.GetNodes())
955 : {
956 : auto BBUpdate =
957 97668 : [&mesh, &dim](const int *v, int nv, mfem::Vector &min, mfem::Vector &max)
958 : {
959 488340 : for (int j = 0; j < nv; j++)
960 : {
961 390672 : const double *coord = mesh.GetVertex(v[j]);
962 1562688 : for (int d = 0; d < dim; d++)
963 : {
964 1172016 : if (coord[d] < min(d))
965 : {
966 251 : min(d) = coord[d];
967 : }
968 1172016 : if (coord[d] > max(d))
969 : {
970 410 : max(d) = coord[d];
971 : }
972 : }
973 : }
974 97686 : };
975 18 : PalacePragmaOmp(parallel)
976 : {
977 : mfem::Vector loc_min(dim), loc_max(dim);
978 : for (int d = 0; d < dim; d++)
979 : {
980 : loc_min(d) = mfem::infinity();
981 : loc_max(d) = -mfem::infinity();
982 : }
983 : if (bdr)
984 : {
985 : PalacePragmaOmp(for schedule(static))
986 : for (int i = 0; i < mesh.GetNBE(); i++)
987 : {
988 : if (!marker[mesh.GetBdrAttribute(i) - 1])
989 : {
990 : continue;
991 : }
992 : const int *verts = mesh.GetBdrElement(i)->GetVertices();
993 : BBUpdate(verts, mesh.GetBdrElement(i)->GetNVertices(), loc_min, loc_max);
994 : }
995 : }
996 : else
997 : {
998 : PalacePragmaOmp(for schedule(static))
999 : for (int i = 0; i < mesh.GetNE(); i++)
1000 : {
1001 : if (!marker[mesh.GetAttribute(i) - 1])
1002 : {
1003 : continue;
1004 : }
1005 : const int *verts = mesh.GetElement(i)->GetVertices();
1006 : BBUpdate(verts, mesh.GetElement(i)->GetNVertices(), loc_min, loc_max);
1007 : }
1008 : }
1009 : PalacePragmaOmp(critical(BBUpdate))
1010 : {
1011 : for (int d = 0; d < dim; d++)
1012 : {
1013 : min(d) = std::min(min(d), loc_min(d));
1014 : max(d) = std::max(max(d), loc_max(d));
1015 : }
1016 : }
1017 : }
1018 : }
1019 : else
1020 : {
1021 55 : mesh.GetNodes()->HostRead();
1022 55 : const int ref = mesh.GetNodes()->FESpace()->GetMaxElementOrder();
1023 98731 : auto BBUpdate = [&ref](mfem::GeometryRefiner &refiner, mfem::ElementTransformation &T,
1024 : mfem::DenseMatrix &pointmat, mfem::Vector &min,
1025 : mfem::Vector &max)
1026 : {
1027 98731 : mfem::RefinedGeometry *RefG = refiner.Refine(T.GetGeometryType(), ref);
1028 98731 : T.Transform(RefG->RefPts, pointmat);
1029 493655 : for (int j = 0; j < pointmat.Width(); j++)
1030 : {
1031 1579696 : for (int d = 0; d < pointmat.Height(); d++)
1032 : {
1033 1184772 : if (pointmat(d, j) < min(d))
1034 : {
1035 711 : min(d) = pointmat(d, j);
1036 : }
1037 1184772 : if (pointmat(d, j) > max(d))
1038 : {
1039 715 : max(d) = pointmat(d, j);
1040 : }
1041 : }
1042 : }
1043 98786 : };
1044 55 : PalacePragmaOmp(parallel)
1045 : {
1046 : mfem::Vector loc_min(dim), loc_max(dim);
1047 : for (int d = 0; d < dim; d++)
1048 : {
1049 : loc_min(d) = mfem::infinity();
1050 : loc_max(d) = -mfem::infinity();
1051 : }
1052 : mfem::GeometryRefiner refiner;
1053 : mfem::IsoparametricTransformation T;
1054 : mfem::DenseMatrix pointmat;
1055 : if (bdr)
1056 : {
1057 : PalacePragmaOmp(for schedule(static))
1058 : for (int i = 0; i < mesh.GetNBE(); i++)
1059 : {
1060 : if (!marker[mesh.GetBdrAttribute(i) - 1])
1061 : {
1062 : continue;
1063 : }
1064 : mesh.GetBdrElementTransformation(i, &T);
1065 : BBUpdate(refiner, T, pointmat, loc_min, loc_max);
1066 : }
1067 : }
1068 : else
1069 : {
1070 : PalacePragmaOmp(for schedule(static))
1071 : for (int i = 0; i < mesh.GetNE(); i++)
1072 : {
1073 : if (!marker[mesh.GetAttribute(i) - 1])
1074 : {
1075 : continue;
1076 : }
1077 : mesh.GetElementTransformation(i, &T);
1078 : BBUpdate(refiner, T, pointmat, loc_min, loc_max);
1079 : }
1080 : }
1081 : PalacePragmaOmp(critical(BBUpdate))
1082 : {
1083 : for (int d = 0; d < dim; d++)
1084 : {
1085 : min(d) = std::min(min(d), loc_min(d));
1086 : max(d) = std::max(max(d), loc_max(d));
1087 : }
1088 : }
1089 : }
1090 : }
1091 73 : Mpi::GlobalMin(dim, min.HostReadWrite(), mesh.GetComm());
1092 73 : Mpi::GlobalMax(dim, max.HostReadWrite(), mesh.GetComm());
1093 73 : }
1094 :
1095 0 : double BoundingBox::Area() const
1096 : {
1097 0 : return 4.0 * CVector3dMap(axes[0].data()).cross(CVector3dMap(axes[1].data())).norm();
1098 : }
1099 :
1100 0 : double BoundingBox::Volume() const
1101 : {
1102 0 : return planar ? 0.0 : 2.0 * CVector3dMap(axes[2].data()).norm() * Area();
1103 : }
1104 :
1105 1 : std::array<std::array<double, 3>, 3> BoundingBox::Normals() const
1106 : {
1107 1 : std::array<std::array<double, 3>, 3> normals = {axes[0], axes[1], axes[2]};
1108 : Vector3dMap(normals[0].data()).normalize();
1109 : Vector3dMap(normals[1].data()).normalize();
1110 : Vector3dMap(normals[2].data()).normalize();
1111 1 : return normals;
1112 : }
1113 :
1114 31 : std::array<double, 3> BoundingBox::Lengths() const
1115 : {
1116 31 : return {2.0 * CVector3dMap(axes[0].data()).norm(),
1117 31 : 2.0 * CVector3dMap(axes[1].data()).norm(),
1118 31 : 2.0 * CVector3dMap(axes[2].data()).norm()};
1119 : }
1120 :
1121 30 : std::array<double, 3> BoundingBox::Deviations(const std::array<double, 3> &direction) const
1122 : {
1123 : const auto eig_dir = CVector3dMap(direction.data());
1124 : std::array<double, 3> deviation_deg;
1125 120 : for (std::size_t i = 0; i < 3; i++)
1126 : {
1127 90 : deviation_deg[i] =
1128 90 : std::acos(std::min(1.0, std::abs(eig_dir.normalized().dot(
1129 90 : CVector3dMap(axes[i].data()).normalized())))) *
1130 : (180.0 / M_PI);
1131 : }
1132 30 : return deviation_deg;
1133 : }
1134 :
1135 30 : BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
1136 : bool bdr)
1137 : {
1138 : std::vector<Eigen::Vector3d> vertices;
1139 30 : int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices);
1140 60 : return BoundingBoxFromPointCloud(mesh.GetComm(), vertices, dominant_rank);
1141 : }
1142 :
1143 0 : BoundingBox GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
1144 : bool bdr)
1145 : {
1146 : std::vector<Eigen::Vector3d> vertices;
1147 0 : int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices);
1148 0 : return BoundingBallFromPointCloud(mesh.GetComm(), vertices, dominant_rank);
1149 : }
1150 :
1151 30 : double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
1152 : bool bdr, const std::array<double, 3> &dir)
1153 : {
1154 : std::vector<Eigen::Vector3d> vertices;
1155 30 : int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices);
1156 : double length;
1157 30 : if (dominant_rank == Mpi::Rank(mesh.GetComm()))
1158 : {
1159 : CVector3dMap direction(dir.data());
1160 : auto Dot = [&](const auto &x, const auto &y)
1161 472 : { return direction.dot(x) < direction.dot(y); };
1162 28 : auto p_min = std::min_element(vertices.begin(), vertices.end(), Dot);
1163 28 : auto p_max = std::max_element(vertices.begin(), vertices.end(), Dot);
1164 28 : length = (*p_max - *p_min).dot(direction.normalized());
1165 : }
1166 : Mpi::Broadcast(1, &length, dominant_rank, mesh.GetComm());
1167 60 : return length;
1168 : }
1169 :
1170 0 : double GetDistanceFromPoint(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
1171 : bool bdr, const std::array<double, 3> &origin, bool max)
1172 : {
1173 : std::vector<Eigen::Vector3d> vertices;
1174 0 : int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices);
1175 : double dist;
1176 0 : if (dominant_rank == Mpi::Rank(mesh.GetComm()))
1177 : {
1178 : CVector3dMap x0(origin.data());
1179 : auto p =
1180 0 : max ? std::max_element(vertices.begin(), vertices.end(),
1181 0 : [&x0](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
1182 0 : { return (x - x0).norm() < (y - x0).norm(); })
1183 0 : : std::min_element(vertices.begin(), vertices.end(),
1184 0 : [&x0](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
1185 0 : { return (x - x0).norm() < (y - x0).norm(); });
1186 0 : dist = (*p - x0).norm();
1187 : }
1188 : Mpi::Broadcast(1, &dist, dominant_rank, mesh.GetComm());
1189 0 : return dist;
1190 : }
1191 :
1192 : // Given a mesh and boundary attribute marker array, compute a normal for the surface. If
1193 : // not averaging, use the first entry.
1194 78 : mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
1195 : bool average)
1196 : {
1197 : int dim = mesh.SpaceDimension();
1198 78 : mfem::IsoparametricTransformation T;
1199 : mfem::Vector loc_normal(dim), normal(dim);
1200 78 : normal = 0.0;
1201 78 : if (mesh.Dimension() == mesh.SpaceDimension())
1202 : {
1203 : // Loop over boundary elements. Exit early if not averaging and non-zero normal.
1204 3966 : for (int i = 0; i < mesh.GetNBE() && !(!average && normal.Norml2() > 0.0); i++)
1205 : {
1206 3888 : if (!marker[mesh.GetBdrAttribute(i) - 1])
1207 : {
1208 3672 : continue;
1209 : }
1210 216 : mesh.GetBdrElementTransformation(i, &T);
1211 216 : mesh::Normal(T, loc_normal, &normal);
1212 216 : normal += loc_normal;
1213 : }
1214 : }
1215 : else
1216 : {
1217 : // Loop over domain elements. Exit early if not averaging and non-zero normal.
1218 0 : for (int i = 0; i < mesh.GetNE() && !(!average && normal.Norml2() > 0.0); i++)
1219 : {
1220 0 : if (!marker[mesh.GetAttribute(i) - 1])
1221 : {
1222 0 : continue;
1223 : }
1224 0 : mesh.GetElementTransformation(i, &T);
1225 0 : mesh::Normal(T, loc_normal, &normal);
1226 0 : normal += loc_normal;
1227 : }
1228 : }
1229 :
1230 : // If different processors have different normal orientations, take that from the lowest
1231 : // rank processor.
1232 : MPI_Comm comm = mesh.GetComm();
1233 78 : int rank = Mpi::Size(comm);
1234 : mfem::Vector glob_normal(dim);
1235 78 : if (normal.Norml2() > 0.0)
1236 : {
1237 78 : rank = Mpi::Rank(comm);
1238 : }
1239 : Mpi::GlobalMin(1, &rank, comm);
1240 78 : if (rank == Mpi::Size(comm))
1241 : {
1242 : // No boundary elements are marked.
1243 0 : normal = 0.0;
1244 : return normal;
1245 : }
1246 78 : if (rank == Mpi::Rank(comm))
1247 : {
1248 76 : glob_normal = normal;
1249 : }
1250 78 : Mpi::Broadcast(dim, glob_normal.HostReadWrite(), rank, comm);
1251 78 : if (average)
1252 : {
1253 78 : if (normal * glob_normal < 0.0)
1254 : {
1255 0 : normal.Neg();
1256 : }
1257 : Mpi::GlobalSum(dim, normal.HostReadWrite(), comm);
1258 : }
1259 : else
1260 : {
1261 0 : normal = glob_normal;
1262 : }
1263 78 : normal /= normal.Norml2();
1264 :
1265 : if constexpr (false)
1266 : {
1267 : Mpi::Print(comm, " Surface normal = ({:+.3e})", fmt::join(normal, ", "));
1268 : }
1269 : return normal;
1270 78 : }
1271 :
1272 30 : double GetSurfaceArea(const mfem::ParMesh &mesh, const mfem::Array<int> &marker)
1273 : {
1274 30 : double area = 0.0;
1275 30 : PalacePragmaOmp(parallel reduction(+ : area))
1276 : {
1277 : mfem::IsoparametricTransformation T;
1278 : PalacePragmaOmp(for schedule(static))
1279 : for (int i = 0; i < mesh.GetNBE(); i++)
1280 : {
1281 : if (!marker[mesh.GetBdrAttribute(i) - 1])
1282 : {
1283 : continue;
1284 : }
1285 : mesh.GetBdrElementTransformation(i, &T);
1286 : const mfem::IntegrationRule &ir = mfem::IntRules.Get(T.GetGeometryType(), T.OrderJ());
1287 : for (int j = 0; j < ir.GetNPoints(); j++)
1288 : {
1289 : const mfem::IntegrationPoint &ip = ir.IntPoint(j);
1290 : T.SetIntPoint(&ip);
1291 : area += ip.weight * T.Weight();
1292 : }
1293 : }
1294 : }
1295 : Mpi::GlobalSum(1, &area, mesh.GetComm());
1296 30 : return area;
1297 : }
1298 :
1299 0 : double GetVolume(const mfem::ParMesh &mesh, const mfem::Array<int> &marker)
1300 : {
1301 0 : double volume = 0.0;
1302 0 : PalacePragmaOmp(parallel reduction(+ : volume))
1303 : {
1304 : mfem::IsoparametricTransformation T;
1305 : PalacePragmaOmp(for schedule(static))
1306 : for (int i = 0; i < mesh.GetNE(); i++)
1307 : {
1308 : if (!marker[mesh.GetAttribute(i) - 1])
1309 : {
1310 : continue;
1311 : }
1312 : mesh.GetElementTransformation(i, &T);
1313 : const mfem::IntegrationRule &ir = mfem::IntRules.Get(T.GetGeometryType(), T.OrderJ());
1314 : for (int j = 0; j < ir.GetNPoints(); j++)
1315 : {
1316 : const mfem::IntegrationPoint &ip = ir.IntPoint(j);
1317 : T.SetIntPoint(&ip);
1318 : volume += ip.weight * T.Weight();
1319 : }
1320 : }
1321 : }
1322 : Mpi::GlobalSum(1, &volume, mesh.GetComm());
1323 0 : return volume;
1324 : }
1325 :
1326 0 : double RebalanceMesh(const IoData &iodata, std::unique_ptr<mfem::ParMesh> &mesh)
1327 : {
1328 0 : BlockTimer bt0(Timer::REBALANCE);
1329 0 : MPI_Comm comm = mesh->GetComm();
1330 0 : if (iodata.model.refinement.save_adapt_mesh)
1331 : {
1332 : // Create a separate serial mesh to write to disk.
1333 0 : auto sfile = fs::path(iodata.problem.output) / fs::path(iodata.model.mesh).stem();
1334 0 : sfile += ".mesh";
1335 :
1336 0 : auto PrintSerial = [&](mfem::Mesh &smesh)
1337 : {
1338 0 : BlockTimer bt1(Timer::IO);
1339 0 : if (Mpi::Root(comm))
1340 : {
1341 0 : std::ofstream fo(sfile);
1342 : // mfem::ofgzstream fo(sfile, true); // Use zlib compression if available
1343 : // fo << std::fixed;
1344 : fo << std::scientific;
1345 : fo.precision(MSH_FLT_PRECISION);
1346 0 : mesh::DimensionalizeMesh(smesh, iodata.units.GetMeshLengthRelativeScale());
1347 0 : smesh.Mesh::Print(fo); // Do not need to nondimensionalize the temporary mesh
1348 0 : }
1349 0 : Mpi::Barrier(comm);
1350 0 : };
1351 :
1352 0 : if (mesh->Nonconforming())
1353 : {
1354 0 : mfem::ParMesh smesh(*mesh);
1355 : mfem::Array<int> serial_partition(mesh->GetNE());
1356 : serial_partition = 0;
1357 0 : smesh.Rebalance(serial_partition);
1358 0 : PrintSerial(smesh);
1359 0 : }
1360 : else
1361 : {
1362 0 : mfem::Mesh smesh = mesh->GetSerialMesh(0);
1363 0 : PrintSerial(smesh);
1364 0 : }
1365 0 : }
1366 :
1367 : // If there is more than one processor, may perform rebalancing.
1368 0 : if (Mpi::Size(comm) == 1)
1369 : {
1370 : return 1.0;
1371 : }
1372 : int min_elem, max_elem;
1373 0 : min_elem = max_elem = mesh->GetNE();
1374 0 : Mpi::GlobalMin(1, &min_elem, comm);
1375 0 : Mpi::GlobalMax(1, &max_elem, comm);
1376 0 : const double ratio = double(max_elem) / min_elem;
1377 0 : const double tol = iodata.model.refinement.maximum_imbalance;
1378 : if constexpr (false)
1379 : {
1380 : Mpi::Print("Rebalancing: max/min elements per processor = {:d}/{:d} (ratio = {:.3e}, "
1381 : "tol = {:.3e})\n",
1382 : max_elem, min_elem, ratio, tol);
1383 : }
1384 0 : if (ratio > tol)
1385 : {
1386 0 : if (mesh->Nonconforming())
1387 : {
1388 0 : mesh->Rebalance();
1389 : }
1390 : else
1391 : {
1392 : // Without access to a refinement tree, partitioning must be done on the root
1393 : // processor and then redistributed.
1394 0 : RebalanceConformalMesh(mesh);
1395 : }
1396 : }
1397 : return ratio;
1398 0 : }
1399 :
1400 : } // namespace mesh
1401 :
1402 : namespace
1403 : {
1404 :
1405 8 : std::unique_ptr<mfem::Mesh> LoadMesh(const std::string &mesh_file, bool remove_curvature,
1406 : const config::BoundaryData &boundaries)
1407 : {
1408 : // Read the (serial) mesh from the given mesh file. Handle preparation for refinement and
1409 : // orientations here to avoid possible reorientations and reordering later on. MFEM
1410 : // supports a native mesh format (.mesh), VTK/VTU, Gmsh, as well as some others. We use
1411 : // built-in converters for the types we know, otherwise rely on MFEM to do the conversion
1412 : // or error out if not supported.
1413 8 : constexpr bool generate_edges = false, refine = false, fix_orientation = true;
1414 8 : std::unique_ptr<mfem::Mesh> mesh;
1415 8 : fs::path mesh_path(mesh_file);
1416 40 : if (mesh_path.extension() == ".mphtxt" || mesh_path.extension() == ".mphbin" ||
1417 40 : mesh_path.extension() == ".nas" || mesh_path.extension() == ".bdf")
1418 : {
1419 : // Put translated mesh in temporary string buffer.
1420 0 : std::stringstream fi(std::stringstream::in | std::stringstream::out);
1421 : // fi << std::fixed;
1422 : fi << std::scientific;
1423 : fi.precision(MSH_FLT_PRECISION);
1424 0 : if (mesh_path.extension() == ".mphtxt" || mesh_path.extension() == ".mphbin")
1425 : {
1426 0 : mesh::ConvertMeshComsol(mesh_file, fi, remove_curvature);
1427 : // mesh::ConvertMeshComsol(mesh_file, fo, remove_curvature);
1428 : }
1429 : else
1430 : {
1431 0 : mesh::ConvertMeshNastran(mesh_file, fi, remove_curvature);
1432 : // mesh::ConvertMeshNastran(mesh_file, fo, remove_curvature);
1433 : }
1434 0 : mesh = std::make_unique<mfem::Mesh>(fi, generate_edges, refine, fix_orientation);
1435 0 : }
1436 : else
1437 : {
1438 : // Otherwise, just rely on MFEM load the mesh.
1439 8 : std::ifstream fi(mesh_file);
1440 8 : if (!fi.good())
1441 : {
1442 0 : MFEM_ABORT("Unable to open mesh file \"" << mesh_file << "\"!");
1443 : }
1444 8 : mesh = std::make_unique<mfem::Mesh>(fi, generate_edges, refine, fix_orientation);
1445 8 : }
1446 8 : if (remove_curvature)
1447 : {
1448 0 : mesh->SetCurvature(-1);
1449 : }
1450 : else
1451 : {
1452 8 : mesh->EnsureNodes();
1453 : }
1454 8 : if (!boundaries.periodic.boundary_pairs.empty())
1455 : {
1456 : auto periodic_mesh = std::move(mesh);
1457 :
1458 0 : for (const auto &data : boundaries.periodic.boundary_pairs)
1459 : {
1460 0 : auto periodic_mapping = mesh::DeterminePeriodicVertexMapping(periodic_mesh, data);
1461 0 : if (!periodic_mapping.empty())
1462 : {
1463 : auto p_mesh = std::make_unique<mfem::Mesh>(
1464 0 : mfem::Mesh::MakePeriodic(*periodic_mesh, periodic_mapping));
1465 : periodic_mesh = std::move(p_mesh);
1466 : }
1467 : }
1468 : mesh = std::move(periodic_mesh);
1469 : }
1470 :
1471 8 : return mesh;
1472 8 : }
1473 :
1474 : template <typename T = std::vector<int>>
1475 0 : void TransferHighOrderNodes(const mfem::Mesh &orig_mesh, mfem::Mesh &new_mesh,
1476 : const T *elem_delete_map = nullptr)
1477 : {
1478 : // This accounts for new boundary elements too since no new dofs are added. See the MFEM
1479 : // trimmer miniapp for reference.
1480 0 : MFEM_VERIFY(orig_mesh.GetNodes(), "No high-order nodes information to transfer!");
1481 : const mfem::GridFunction *nodes = orig_mesh.GetNodes();
1482 : const mfem::FiniteElementSpace *fespace = nodes->FESpace();
1483 : mfem::Ordering::Type ordering = fespace->GetOrdering();
1484 0 : int order = fespace->GetMaxElementOrder();
1485 : int sdim = orig_mesh.SpaceDimension();
1486 0 : bool discont =
1487 0 : (dynamic_cast<const mfem::L2_FECollection *>(fespace->FEColl()) != nullptr);
1488 0 : new_mesh.SetCurvature(order, discont, sdim, ordering);
1489 : mfem::GridFunction *new_nodes = new_mesh.GetNodes();
1490 : const mfem::FiniteElementSpace *new_fespace = new_nodes->FESpace();
1491 :
1492 : // Transfer dofs from the old mesh to the new ones. Either consider all elements (works
1493 : // for orientation or numbering changes), or use the prescribed old to new element index
1494 : // map.
1495 : mfem::Array<int> vdofs;
1496 : mfem::Vector loc_vec;
1497 0 : for (int e = 0; e < orig_mesh.GetNE(); e++)
1498 : {
1499 0 : if (!elem_delete_map || (*elem_delete_map)[e] >= 0)
1500 : {
1501 : // No need for DofTransformation here since spaces are H1 or L2.
1502 0 : fespace->GetElementVDofs(e, vdofs);
1503 0 : nodes->GetSubVector(vdofs, loc_vec);
1504 0 : new_fespace->GetElementVDofs(!elem_delete_map ? e : (*elem_delete_map)[e], vdofs);
1505 0 : new_nodes->SetSubVector(vdofs, loc_vec);
1506 : }
1507 : }
1508 0 : }
1509 :
1510 8 : void CleanMesh(std::unique_ptr<mfem::Mesh> &orig_mesh,
1511 : const std::vector<int> &mat_attr_list)
1512 : {
1513 : auto mat_marker = mesh::AttrToMarker(
1514 8 : orig_mesh->attributes.Size() ? orig_mesh->attributes.Max() : 0, mat_attr_list, true);
1515 8 : std::vector<int> elem_delete_map(orig_mesh->GetNE(), -1),
1516 8 : bdr_elem_delete_map(orig_mesh->GetNBE(), -1);
1517 :
1518 : // Delete domain and boundary elements which have no associated material or BC attribute
1519 : // from the mesh.
1520 8 : int new_ne = 0;
1521 344 : for (int e = 0; e < orig_mesh->GetNE(); e++)
1522 : {
1523 336 : if (mat_marker[orig_mesh->GetAttribute(e) - 1])
1524 : {
1525 336 : elem_delete_map[e] = new_ne++;
1526 : }
1527 : }
1528 :
1529 : // Make sure to remove any boundary elements which are no longer attached to elements in
1530 : // the domain.
1531 8 : int new_nbe = 0;
1532 392 : for (int be = 0; be < orig_mesh->GetNBE(); be++)
1533 : {
1534 : int f, o, e1, e2;
1535 384 : orig_mesh->GetBdrElementFace(be, &f, &o);
1536 384 : orig_mesh->GetFaceElements(f, &e1, &e2);
1537 384 : bool no_e1 = (e1 < 0 || elem_delete_map[e1] < 0);
1538 384 : bool no_e2 = (e2 < 0 || elem_delete_map[e2] < 0);
1539 384 : if (!no_e1 || !no_e2)
1540 : {
1541 384 : bdr_elem_delete_map[be] = new_nbe++;
1542 : }
1543 : else if constexpr (false)
1544 : {
1545 : Mpi::Print("Deleting an unattached boundary element!\n");
1546 : }
1547 : }
1548 8 : if (new_ne < orig_mesh->GetNE())
1549 : {
1550 0 : Mpi::Print("Removed {:d} unmarked domain elements from the mesh\n",
1551 0 : orig_mesh->GetNE() - new_ne);
1552 : }
1553 8 : if (new_nbe < orig_mesh->GetNBE())
1554 : {
1555 0 : Mpi::Print("Removed {:d} unattached boundary elements from the mesh\n",
1556 0 : orig_mesh->GetNBE() - new_nbe);
1557 : }
1558 :
1559 : // Create the new mesh.
1560 8 : if (new_ne == orig_mesh->GetNE() && new_nbe == orig_mesh->GetNBE())
1561 : {
1562 : return;
1563 : }
1564 0 : MFEM_VERIFY(!orig_mesh->Nonconforming(),
1565 : "Mesh element cleaning is not supported for nonconforming meshes!");
1566 : auto new_mesh =
1567 0 : std::make_unique<mfem::Mesh>(orig_mesh->Dimension(), orig_mesh->GetNV(), new_ne,
1568 0 : new_nbe, orig_mesh->SpaceDimension());
1569 :
1570 : // Copy vertices and non-deleted domain and boundary elements.
1571 0 : for (int v = 0; v < orig_mesh->GetNV(); v++)
1572 : {
1573 0 : new_mesh->AddVertex(orig_mesh->GetVertex(v));
1574 : }
1575 0 : for (int e = 0; e < orig_mesh->GetNE(); e++)
1576 : {
1577 0 : if (elem_delete_map[e] >= 0)
1578 : {
1579 0 : mfem::Element *el = orig_mesh->GetElement(e)->Duplicate(new_mesh.get());
1580 0 : new_mesh->AddElement(el);
1581 : }
1582 : }
1583 0 : for (int be = 0; be < orig_mesh->GetNBE(); be++)
1584 : {
1585 0 : if (bdr_elem_delete_map[be] >= 0)
1586 : {
1587 0 : mfem::Element *bdr_el = orig_mesh->GetBdrElement(be)->Duplicate(new_mesh.get());
1588 0 : new_mesh->AddBdrElement(bdr_el);
1589 : }
1590 : }
1591 :
1592 : // Finalize the new mesh topology and replace the old mesh. If a curved mesh, set up the
1593 : // new mesh by projecting nodes onto the new mesh for the non-trimmed vdofs. No need to
1594 : // mark for refinement or fix orientations, since everything is copied from the previous
1595 : // mesh.
1596 : constexpr bool generate_bdr = false;
1597 0 : new_mesh->FinalizeTopology(generate_bdr);
1598 0 : new_mesh->RemoveUnusedVertices(); // Remove vertices from the deleted elements
1599 0 : if (orig_mesh->GetNodes())
1600 : {
1601 0 : TransferHighOrderNodes(*orig_mesh, *new_mesh, &elem_delete_map);
1602 : }
1603 : orig_mesh = std::move(new_mesh);
1604 : }
1605 :
1606 0 : void SplitMeshElements(std::unique_ptr<mfem::Mesh> &orig_mesh, bool make_simplex,
1607 : bool make_hex)
1608 : {
1609 0 : if (!make_simplex && !make_hex)
1610 : {
1611 0 : return;
1612 : }
1613 : mfem::Mesh *mesh = orig_mesh.get();
1614 0 : mfem::Mesh new_mesh;
1615 :
1616 : // Convert all element types to simplices.
1617 0 : if (make_simplex)
1618 : {
1619 0 : const auto element_types = mesh::CheckElements(*mesh);
1620 0 : if (element_types.has_hexahedra || element_types.has_prisms ||
1621 0 : element_types.has_pyramids)
1622 : {
1623 0 : MFEM_VERIFY(!mesh->Nonconforming(),
1624 : "Mesh element splitting is not supported for nonconforming meshes!");
1625 0 : MFEM_VERIFY(
1626 : !element_types.has_pyramids,
1627 : "Splitting mesh elements to simplices does not support pyramid elements yet!");
1628 : int ne = mesh->GetNE();
1629 0 : new_mesh = mfem::Mesh::MakeSimplicial(*mesh);
1630 0 : Mpi::Print("Added {:d} elements to the mesh during conversion to simplices\n",
1631 0 : new_mesh.GetNE() - ne);
1632 : mesh = &new_mesh;
1633 : }
1634 : }
1635 :
1636 : // Convert all element types to hexahedra (currently only tet-to-hex).
1637 0 : if (make_hex)
1638 : {
1639 0 : const auto element_types = mesh::CheckElements(*mesh);
1640 0 : if (element_types.has_simplices || element_types.has_prisms ||
1641 0 : element_types.has_pyramids)
1642 : {
1643 0 : MFEM_VERIFY(!mesh->Nonconforming(),
1644 : "Mesh element splitting is not supported for nonconforming meshes!");
1645 0 : MFEM_VERIFY(!element_types.has_prisms && !element_types.has_pyramids,
1646 : "Splitting mesh elements to hexahedra only supports simplex elements "
1647 : "(tetrahedra) for now!");
1648 : int ne = mesh->GetNE();
1649 0 : new_mesh = mesh::MeshTetToHex(*mesh);
1650 0 : Mpi::Print("Added {:d} elements to the mesh during conversion to hexahedra\n",
1651 0 : new_mesh.GetNE() - ne);
1652 : mesh = &new_mesh;
1653 : }
1654 : }
1655 :
1656 : // Return if no modifications were made.
1657 0 : if (mesh == orig_mesh.get())
1658 : {
1659 : return;
1660 : }
1661 0 : orig_mesh = std::make_unique<mfem::Mesh>(std::move(new_mesh)); // Call move constructor
1662 0 : orig_mesh->FinalizeTopology();
1663 0 : }
1664 :
1665 0 : void ReorderMeshElements(mfem::Mesh &mesh, bool print)
1666 : {
1667 : mfem::Array<int> ordering;
1668 : if constexpr (false)
1669 : {
1670 : // Gecko reordering.
1671 : mfem::Array<int> tentative;
1672 : int outer = 3, inner = 3, window = 4, period = 2;
1673 : double best_cost = mfem::infinity();
1674 : for (int i = 0; i < outer; i++)
1675 : {
1676 : int seed = i + 1;
1677 : double cost =
1678 : mesh.GetGeckoElementOrdering(tentative, inner, window, period, seed, true);
1679 : if (cost < best_cost)
1680 : {
1681 : ordering = tentative;
1682 : best_cost = cost;
1683 : }
1684 : }
1685 : if (print)
1686 : {
1687 : Mpi::Print("Final cost: {:e}\n", best_cost);
1688 : }
1689 : }
1690 : else
1691 : {
1692 : // (Faster) Hilbert reordering.
1693 0 : mesh.GetHilbertElementOrdering(ordering);
1694 0 : mesh.ReorderElements(ordering);
1695 : }
1696 0 : }
1697 :
1698 8 : std::unordered_map<int, int> GetFaceToBdrElementMap(const mfem::Mesh &mesh,
1699 : const config::BoundaryData &boundaries)
1700 : {
1701 : std::unordered_map<int, int> face_to_be;
1702 8 : face_to_be.reserve(mesh.GetNBE());
1703 392 : for (int be = 0; be < mesh.GetNBE(); be++)
1704 : {
1705 384 : int f, o, e1 = -1, e2 = -1;
1706 384 : mesh.GetBdrElementFace(be, &f, &o);
1707 384 : int attr = mesh.GetBdrAttribute(be);
1708 384 : if (!boundaries.periodic.boundary_pairs.empty())
1709 : {
1710 0 : for (const auto &data : boundaries.periodic.boundary_pairs)
1711 : {
1712 : const auto &da = data.donor_attributes, &ra = data.receiver_attributes;
1713 0 : auto donor = std::find(da.begin(), da.end(), attr) != da.end();
1714 0 : auto receiver = std::find(ra.begin(), ra.end(), attr) != ra.end();
1715 0 : if (donor || receiver)
1716 : {
1717 0 : mesh.GetFaceElements(f, &e1, &e2);
1718 0 : MFEM_VERIFY(e1 >= 0 && e2 >= 0,
1719 : "Mesh is not periodic on attribute " << attr << "!");
1720 : }
1721 : }
1722 : }
1723 384 : MFEM_VERIFY((e1 >= 0 && e2 >= 0) || face_to_be.find(f) == face_to_be.end(),
1724 : "A non-periodic face ("
1725 : << f << ") cannot have multiple boundary elements! Attributes: " << attr
1726 : << ' ' << mesh.GetBdrAttribute(face_to_be[f]));
1727 384 : face_to_be[f] = be;
1728 : }
1729 8 : return face_to_be;
1730 : }
1731 :
1732 8 : std::unordered_map<int, int> CheckMesh(const mfem::Mesh &mesh,
1733 : const config::BoundaryData &boundaries)
1734 : {
1735 : // Check for:
1736 : // (1) Boundary elements with no prescribed boundary condition, and
1737 : // (2) Boundary faces which have no boundary element.
1738 : auto bdr_marker =
1739 8 : mesh::AttrToMarker(mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0,
1740 16 : boundaries.attributes, true);
1741 8 : std::unordered_map<int, int> face_to_be = GetFaceToBdrElementMap(mesh, boundaries);
1742 : std::unordered_set<int> bdr_warn_list;
1743 8 : int bdr_face_warn = 0;
1744 872 : for (int f = 0; f < mesh.GetNumFaces(); f++)
1745 : {
1746 : int e1, e2;
1747 864 : mesh.GetFaceElements(f, &e1, &e2);
1748 864 : if (e1 >= 0 && e2 >= 0)
1749 : {
1750 480 : continue; // Only consider true exterior faces
1751 : }
1752 : auto it = face_to_be.find(f);
1753 384 : if (it != face_to_be.end())
1754 : {
1755 384 : int attr = mesh.GetBdrAttribute(it->second);
1756 384 : if (!bdr_marker[attr - 1])
1757 : {
1758 : // Boundary element with no prescribed boundary condition.
1759 : bdr_warn_list.insert(attr);
1760 : }
1761 : }
1762 : else
1763 : {
1764 : // Boundary face with no attached boundary element.
1765 0 : bdr_face_warn++;
1766 : }
1767 : }
1768 8 : if (!bdr_warn_list.empty())
1769 : {
1770 8 : Mpi::Warning("One or more external boundary attributes has no associated boundary "
1771 : "condition!\n\"PMC\"/\"ZeroCharge\" condition is assumed!\n");
1772 16 : utils::PrettyPrint(bdr_warn_list, "Boundary attribute list:");
1773 8 : Mpi::Print("\n");
1774 : }
1775 8 : if (bdr_face_warn)
1776 : {
1777 0 : Mpi::Warning("{:d} mesh faces with no associated boundary element exist on the domain "
1778 : "boundary!\n",
1779 : bdr_face_warn);
1780 : }
1781 8 : return face_to_be;
1782 : }
1783 :
1784 : template <typename T>
1785 0 : class EdgeRefinementMesh : public mfem::Mesh
1786 : {
1787 : private:
1788 : // Container with keys being pairs of vertex indicies (match row/column indices of v_to_v)
1789 : // of edges which we desire to be refined.
1790 : const T &refinement_edges;
1791 :
1792 0 : void MarkTetMeshForRefinement(const mfem::DSTable &v_to_v) override
1793 : {
1794 : // The standard tetrahedral refinement in MFEM is to mark the longest edge of each tet
1795 : // for the first refinement. We hack this marking here in order to prioritize refinement
1796 : // of edges which are part of internal boundary faces being marked for refinement. This
1797 : // should hopefully limit the amount of extra refinement required to ensure conformity
1798 : // after the marked elements are refined. Marking will then discover only the longest
1799 : // edges, which are those within the boundary to be cracked.
1800 : mfem::Array<mfem::real_t> lengths;
1801 0 : GetEdgeLengths2(v_to_v, lengths);
1802 0 : const auto min_length = 0.01 * lengths.Min();
1803 0 : for (int i = 0; i < v_to_v.NumberOfRows(); i++)
1804 : {
1805 0 : for (mfem::DSTable::RowIterator it(v_to_v, i); !it; ++it)
1806 : {
1807 : int j = it.Column();
1808 0 : if (refinement_edges.find({i, j}) == refinement_edges.end())
1809 : {
1810 : // "Zero" the edge lengths which do not connect vertices on the interface. Avoid
1811 : // zero-length edges just in case.
1812 0 : lengths[it.Index()] = min_length;
1813 : }
1814 : }
1815 : }
1816 :
1817 : // Finish marking (see mfem::Mesh::MarkTetMeshForRefinement).
1818 0 : mfem::Array<int> indices(NumOfEdges);
1819 : std::iota(indices.begin(), indices.end(), 0);
1820 0 : for (int i = 0; i < NumOfElements; i++)
1821 : {
1822 0 : if (elements[i]->GetType() == mfem::Element::TETRAHEDRON)
1823 : {
1824 : MFEM_ASSERT(dynamic_cast<mfem::Tetrahedron *>(elements[i]),
1825 : "Unexpected non-Tetrahedron element type!");
1826 0 : static_cast<mfem::Tetrahedron *>(elements[i])->MarkEdge(v_to_v, lengths, indices);
1827 : }
1828 : }
1829 0 : for (int i = 0; i < NumOfBdrElements; i++)
1830 : {
1831 0 : if (boundary[i]->GetType() == mfem::Element::TRIANGLE)
1832 : {
1833 : MFEM_ASSERT(dynamic_cast<mfem::Triangle *>(boundary[i]),
1834 : "Unexpected non-Triangle element type!");
1835 0 : static_cast<mfem::Triangle *>(boundary[i])->MarkEdge(v_to_v, lengths, indices);
1836 : }
1837 : }
1838 0 : }
1839 :
1840 : public:
1841 0 : EdgeRefinementMesh(mfem::Mesh &&mesh, const T &refinement_edges)
1842 0 : : mfem::Mesh(std::move(mesh)), refinement_edges(refinement_edges)
1843 : {
1844 : }
1845 : };
1846 :
1847 : template <typename T>
1848 : struct UnorderedPair
1849 : {
1850 : T first, second;
1851 0 : UnorderedPair(T first, T second) : first(first), second(second) {}
1852 : bool operator==(const UnorderedPair &v) const
1853 : {
1854 0 : return ((v.first == first && v.second == second) ||
1855 0 : (v.first == second && v.second == first));
1856 : }
1857 : };
1858 :
1859 : template <typename T>
1860 : struct UnorderedPairHasher
1861 : {
1862 : std::size_t operator()(const UnorderedPair<T> &v) const
1863 : {
1864 : // Simple hash function for a pair, see https://tinyurl.com/2k4phapb.
1865 : return std::hash<T>()(std::min(v.first, v.second)) ^
1866 0 : std::hash<T>()(std::max(v.first, v.second)) << 1;
1867 : }
1868 : };
1869 :
1870 8 : int AddInterfaceBdrElements(const IoData &iodata, std::unique_ptr<mfem::Mesh> &orig_mesh,
1871 : std::unordered_map<int, int> &face_to_be, MPI_Comm comm)
1872 : {
1873 : // Exclude some internal boundary conditions for which cracking would give invalid
1874 : // results: lumpedports in particular.
1875 8 : const auto crack_boundary_attributes = [&iodata]()
1876 : {
1877 8 : auto cba = iodata.boundaries.attributes;
1878 32 : for (const auto &[idx, data] : iodata.boundaries.lumpedport)
1879 : {
1880 48 : for (const auto &e : data.elements)
1881 : {
1882 144 : auto attr_in_elem = [&](auto x)
1883 : {
1884 144 : return std::find(e.attributes.begin(), e.attributes.end(), x) !=
1885 144 : e.attributes.end();
1886 : };
1887 24 : cba.erase(std::remove_if(cba.begin(), cba.end(), attr_in_elem), cba.end());
1888 : }
1889 : }
1890 8 : return cba;
1891 8 : }();
1892 :
1893 : // Return if nothing to do. Otherwise, count vertices and boundary elements to add.
1894 8 : if (crack_boundary_attributes.empty() && !iodata.model.add_bdr_elements)
1895 : {
1896 : return 1; // Success
1897 : }
1898 :
1899 8 : if (face_to_be.size() != static_cast<std::size_t>(orig_mesh->GetNBE()))
1900 : {
1901 0 : face_to_be = GetFaceToBdrElementMap(*orig_mesh, iodata.boundaries);
1902 : }
1903 8 : int new_nv = orig_mesh->GetNV();
1904 8 : int new_nbe = orig_mesh->GetNBE();
1905 :
1906 : // Duplicate internal boundary elements from the given boundary attribute list, cracking
1907 : // the mesh such that domain elements on either side are no longer coupled. Correctly
1908 : // handles cracks where more than two domains intersect as well as seams where a crack
1909 : // ends and the vertices are not duplicated.
1910 : std::unordered_set<int> crack_bdr_elem;
1911 : std::unordered_map<int, std::vector<std::pair<int, std::unordered_set<int>>>>
1912 : crack_vert_duplicates;
1913 8 : std::unique_ptr<mfem::Table> vert_to_elem;
1914 8 : if (!crack_boundary_attributes.empty() && iodata.model.crack_bdr_elements)
1915 : {
1916 : auto crack_bdr_marker = mesh::AttrToMarker(
1917 0 : orig_mesh->bdr_attributes.Size() ? orig_mesh->bdr_attributes.Max() : 0,
1918 0 : crack_boundary_attributes, true);
1919 0 : for (int be = 0; be < orig_mesh->GetNBE(); be++)
1920 : {
1921 0 : if (crack_bdr_marker[orig_mesh->GetBdrAttribute(be) - 1])
1922 : {
1923 : int f, o, e1, e2;
1924 0 : orig_mesh->GetBdrElementFace(be, &f, &o);
1925 0 : orig_mesh->GetFaceElements(f, &e1, &e2);
1926 0 : if (e1 >= 0 && e2 >= 0)
1927 : {
1928 : crack_bdr_elem.insert(be);
1929 : }
1930 : }
1931 : }
1932 0 : MFEM_VERIFY(crack_bdr_elem.empty() || !orig_mesh->Nonconforming(),
1933 : "Duplicating internal boundary elements for interior boundaries is not "
1934 : "supported for nonconforming meshes!");
1935 :
1936 0 : vert_to_elem.reset(orig_mesh->GetVertexToElementTable()); // Owned by caller
1937 0 : const mfem::Table &elem_to_face = orig_mesh->ElementToFaceTable();
1938 0 : int new_nv_dups = 0;
1939 0 : for (auto be : crack_bdr_elem)
1940 : {
1941 : const mfem::Element *bdr_el = orig_mesh->GetBdrElement(be);
1942 : const int *verts = bdr_el->GetVertices();
1943 0 : for (int i = 0; i < bdr_el->GetNVertices(); i++)
1944 : {
1945 : // Skip vertices we have already processed.
1946 0 : const auto v = verts[i];
1947 0 : if (crack_vert_duplicates.find(v) != crack_vert_duplicates.end())
1948 : {
1949 0 : continue;
1950 : }
1951 :
1952 : // Collect connected components of elements connected to the vertex. Perform BFS
1953 : // on graph of all elements connected to this vertex, where adjacencies are
1954 : // determined by face connectivity excluding the crack faces.
1955 : std::vector<std::unordered_set<int>> components;
1956 : const int *elems = vert_to_elem->GetRow(v);
1957 0 : std::unordered_set<int> unvisited(elems, elems + vert_to_elem->RowSize(v));
1958 0 : while (!unvisited.empty())
1959 : {
1960 0 : auto &component = components.emplace_back();
1961 : component.reserve(unvisited.size());
1962 : std::queue<int> que;
1963 : que.push(*unvisited.begin());
1964 : unvisited.erase(unvisited.begin());
1965 0 : while (!que.empty())
1966 : {
1967 : // Process the current node.
1968 0 : int e = que.front();
1969 : que.pop();
1970 : component.insert(e);
1971 :
1972 : // Add neighbors.
1973 : const int *faces = elem_to_face.GetRow(e);
1974 0 : for (int j = 0; j < elem_to_face.RowSize(e); j++)
1975 : {
1976 0 : const auto f = faces[j];
1977 : {
1978 : auto it = face_to_be.find(f);
1979 0 : if (it != face_to_be.end() &&
1980 0 : crack_bdr_elem.find(it->second) != crack_bdr_elem.end())
1981 : {
1982 : // Skip element-element connectivities which cross the crack.
1983 : continue;
1984 : }
1985 : }
1986 : int e1, e2;
1987 0 : orig_mesh->GetFaceElements(f, &e1, &e2);
1988 0 : MFEM_VERIFY(
1989 : e == e1 || e == e2,
1990 : "Unexpected face-element connectivity in internal boundary cracking!");
1991 0 : int nbr = (e == e1) ? e2 : e1;
1992 0 : if (nbr >= 0)
1993 : {
1994 : auto it = unvisited.find(nbr);
1995 0 : if (it != unvisited.end())
1996 : {
1997 : que.push(nbr);
1998 : unvisited.erase(it);
1999 : }
2000 : }
2001 : }
2002 : }
2003 : }
2004 0 : MFEM_VERIFY(
2005 : !components.empty(),
2006 : "No connected components found for elements adjacent to a crack vertex!");
2007 : #if defined(MFEM_DEBUG)
2008 : {
2009 : std::size_t visited_size = 0;
2010 : for (const auto &component : components)
2011 : {
2012 : visited_size += component.size();
2013 : }
2014 : MFEM_VERIFY(visited_size == static_cast<std::size_t>(vert_to_elem->RowSize(v)),
2015 : "Failed to visit all elements in neighborhood of vertex when "
2016 : "counting connected components!");
2017 : }
2018 : #endif
2019 :
2020 : // Save mapping from original vertex to duplicate vertices, and the corresponding
2021 : // element groupings requiring renumbering. The first group doesn't need
2022 : // renumbering so is not saved. We still keep entries for non-duplicated crack
2023 : // vertices in the set to track them as processed, however.
2024 0 : auto &vert_components = crack_vert_duplicates.try_emplace(v).first->second;
2025 0 : for (auto it = components.begin() + 1; it != components.end(); ++it)
2026 : {
2027 0 : vert_components.emplace_back(-1, std::move(*it));
2028 0 : new_nv_dups++;
2029 : }
2030 0 : }
2031 : }
2032 :
2033 : // After processing all boundary elements, check if there are any elements which need
2034 : // refinement in order to successfully decouple both sides. This happens if we have an
2035 : // edge interior to the crack which connects to seam vertices (non-duplicated vertices
2036 : // attached to crack boundary elements). A previous implementation of refinement
2037 : // considered for refinement just all boundary elements with all attached vertices
2038 : // lying on the seam, but this doesn't catch all required cases.
2039 0 : if (iodata.model.refine_crack_elements)
2040 : {
2041 : std::unordered_map<UnorderedPair<int>, std::vector<int>, UnorderedPairHasher<int>>
2042 : coarse_crack_edge_to_be;
2043 0 : for (auto be : crack_bdr_elem)
2044 : {
2045 : const mfem::Element *bdr_el = orig_mesh->GetBdrElement(be);
2046 : const int *verts = bdr_el->GetVertices();
2047 0 : for (int i = 0; i < bdr_el->GetNEdges(); i++)
2048 : {
2049 0 : auto v0 = verts[bdr_el->GetEdgeVertices(i)[0]],
2050 0 : v1 = verts[bdr_el->GetEdgeVertices(i)[1]];
2051 : MFEM_ASSERT(crack_vert_duplicates.find(v0) != crack_vert_duplicates.end() &&
2052 : crack_vert_duplicates.find(v1) != crack_vert_duplicates.end(),
2053 : "Unable to locate crack vertices for an interior boundary element!");
2054 0 : if (crack_vert_duplicates[v0].empty() && crack_vert_duplicates[v1].empty())
2055 : {
2056 : // This is a seam edge, so add the attached boundary element to a list. The
2057 : // check for the edge being interior to the crack is indicated by visiting more
2058 : // than once.
2059 0 : auto it = coarse_crack_edge_to_be.find({v0, v1});
2060 : auto &adjacent_be =
2061 : (it == coarse_crack_edge_to_be.end())
2062 0 : ? coarse_crack_edge_to_be.try_emplace({v0, v1}).first->second
2063 0 : : it->second;
2064 : adjacent_be.push_back(be);
2065 : }
2066 : }
2067 : }
2068 0 : for (auto it = coarse_crack_edge_to_be.begin(); it != coarse_crack_edge_to_be.end();)
2069 : {
2070 : // Remove all seam edges which are on the "outside" of the crack (visited only
2071 : // once).
2072 0 : if (it->second.size() == 1)
2073 : {
2074 : it = coarse_crack_edge_to_be.erase(it);
2075 : }
2076 : else
2077 : {
2078 : ++it;
2079 : }
2080 : }
2081 : // Static reporting variables so can persist across retries.
2082 : static int new_ne_ref = 0;
2083 : static int new_ref_its = 0;
2084 0 : if (!coarse_crack_edge_to_be.empty())
2085 : {
2086 : // Locally refine the mesh using conformal refinement. If necessary, convert the
2087 : // mesh to simplices first to enable conforming refinement (this will do nothing
2088 : // if the mesh is already a simplex mesh).
2089 : // Note: Eventually we can implement manual conforming face refinement of pairs of
2090 : // elements sharing a face for all element types (insert a vertex at the boundary
2091 : // element center and connect it to all other element vertices). For now, this adds
2092 : // complexity and making use of conformal simplex refinement seems good enough for
2093 : // most use cases.
2094 : int ne = orig_mesh->GetNE();
2095 0 : SplitMeshElements(orig_mesh, true, false);
2096 0 : if (ne != orig_mesh->GetNE())
2097 : {
2098 : face_to_be.clear();
2099 0 : return 0; // Mesh was converted to simplices, start over
2100 : }
2101 : std::unordered_map<int, int> elem_to_refine;
2102 0 : for (const auto &[edge, adjacent_be] : coarse_crack_edge_to_be)
2103 : {
2104 0 : for (auto be : adjacent_be)
2105 : {
2106 : int f, o, e1, e2;
2107 0 : orig_mesh->GetBdrElementFace(be, &f, &o);
2108 0 : orig_mesh->GetFaceElements(f, &e1, &e2);
2109 : MFEM_ASSERT(e1 >= 0 && e2 >= 0,
2110 : "Invalid internal boundary element connectivity!");
2111 0 : elem_to_refine[e1]++; // Value-initialized to 0 at first access
2112 0 : elem_to_refine[e2]++;
2113 : }
2114 : }
2115 : mfem::Array<mfem::Refinement> refinements;
2116 0 : refinements.Reserve(elem_to_refine.size());
2117 0 : for (const auto &[e, count] : elem_to_refine)
2118 : {
2119 : // Tetrahedral bisection (vs. default octasection) will result in fewer added
2120 : // elements at the cost of a potential minor mesh quality degredation.
2121 0 : refinements.Append(mfem::Refinement(e, mfem::Refinement::X));
2122 : // refinements.Append(mfem::Refinement(e, (count > 1) ? mfem::Refinement::XY
2123 : // : mfem::Refinement::X));
2124 : }
2125 0 : if (mesh::CheckElements(*orig_mesh).has_simplices)
2126 : {
2127 : // Mark tetrahedral mesh for refinement before doing local refinement. This is a
2128 : // bit of a strange pattern to override the standard conforming refinement of the
2129 : // mfem::Mesh class. We want to implement our own edge marking of the tetrahedra,
2130 : // so we move the mesh to a constructed derived class object, mark it, and then
2131 : // move assign it to the original base class object before refining. All of these
2132 : // moves should be cheap without any extra memory allocation. Also, we mark the
2133 : // mesh every time to ensure multiple rounds of refinement target the interior
2134 : // boundary (we don't care about preserving the refinement hierarchy).
2135 : constexpr bool refine = true, fix_orientation = false;
2136 : EdgeRefinementMesh ref_mesh(std::move(*orig_mesh), coarse_crack_edge_to_be);
2137 0 : ref_mesh.Finalize(refine, fix_orientation);
2138 0 : *orig_mesh = std::move(ref_mesh);
2139 : }
2140 0 : orig_mesh->GeneralRefinement(refinements, 0);
2141 0 : new_ne_ref += orig_mesh->GetNE() - ne;
2142 0 : new_ref_its++;
2143 : face_to_be.clear();
2144 : return 0; // Mesh was refined (conformally), start over
2145 : }
2146 0 : else if (new_ne_ref > 0)
2147 : {
2148 0 : Mpi::Print(
2149 : "Added {:d} elements in {:d} iterations of local bisection for under-resolved "
2150 : "interior boundaries\n",
2151 : new_ne_ref, new_ref_its);
2152 : }
2153 : }
2154 :
2155 0 : new_nv += new_nv_dups;
2156 0 : new_nbe += crack_bdr_elem.size();
2157 0 : if (crack_bdr_elem.size() > 0)
2158 : {
2159 0 : Mpi::Print("Added {:d} duplicate vertices for interior boundaries in the mesh\n",
2160 : new_nv_dups);
2161 0 : Mpi::Print(
2162 : "Added {:d} duplicate boundary elements for interior boundaries in the mesh\n",
2163 0 : crack_bdr_elem.size());
2164 : }
2165 : }
2166 :
2167 : // Add new boundary elements at material interfaces or on the exterior boundary of the
2168 : // simulation domain, if there is not already a boundary element present.
2169 : std::unordered_map<int, int> new_face_bdr_elem;
2170 8 : if (iodata.model.add_bdr_elements)
2171 : {
2172 8 : int new_nbe_ext = 0, new_nbe_int = 0;
2173 872 : for (int f = 0; f < orig_mesh->GetNumFaces(); f++)
2174 : {
2175 : // Skip all faces which already have an associated boundary element (this includes
2176 : // any boundary elements which were duplicated during cracking in the previous step).
2177 864 : if (face_to_be.find(f) != face_to_be.end())
2178 : {
2179 384 : continue;
2180 : }
2181 : int e1, e2;
2182 480 : orig_mesh->GetFaceElements(f, &e1, &e2);
2183 480 : if (e1 < 0 || e2 < 0)
2184 : {
2185 : if constexpr (false)
2186 : {
2187 : Mpi::Print("Adding exterior boundary element!\n");
2188 : }
2189 0 : new_face_bdr_elem[f] = 1;
2190 0 : new_nbe_ext++;
2191 : }
2192 480 : else if (orig_mesh->GetAttribute(e1) != orig_mesh->GetAttribute(e2))
2193 : {
2194 : if constexpr (false)
2195 : {
2196 : Mpi::Print("Adding material interface boundary element!\n");
2197 : }
2198 0 : new_face_bdr_elem[f] = 1;
2199 0 : new_nbe_int++;
2200 : }
2201 : }
2202 8 : MFEM_VERIFY(new_nbe_ext + new_nbe_int == 0 || !orig_mesh->Nonconforming(),
2203 : "Adding material interface boundary elements is not supported for "
2204 : "nonconforming meshes!");
2205 :
2206 8 : new_nbe += (new_nbe_ext + new_nbe_int);
2207 8 : if (new_nbe_ext > 0)
2208 : {
2209 0 : Mpi::Print("Added {:d} boundary elements for exterior boundaries to the mesh\n",
2210 : new_nbe_ext);
2211 : }
2212 8 : if (new_nbe_int > 0)
2213 : {
2214 0 : Mpi::Print("Added {:d} boundary elements for material interfaces to the mesh\n",
2215 : new_nbe_int);
2216 : }
2217 : }
2218 :
2219 : // Export mesh after pre-processing, before cracking boundary elements.
2220 8 : if (iodata.model.export_prerefined_mesh && Mpi::Root(comm))
2221 : {
2222 0 : auto pos = iodata.model.mesh.find_last_of(".");
2223 0 : std::string meshfile = iodata.model.mesh.substr(0, pos) + "_preprocessed.mesh";
2224 0 : std::ofstream fo(meshfile);
2225 : fo.precision(MSH_FLT_PRECISION);
2226 0 : orig_mesh->Print(fo);
2227 0 : }
2228 :
2229 : // Create the new mesh. We can't just add the new vertices and boundary elements to the
2230 : // original mesh, because we need to keep it around in order to transfer the high-order
2231 : // nodes information to the new mesh.
2232 8 : if (new_nv == orig_mesh->GetNV() && new_nbe == orig_mesh->GetNBE())
2233 : {
2234 : return 1; // Success
2235 : }
2236 : auto new_mesh =
2237 0 : std::make_unique<mfem::Mesh>(orig_mesh->Dimension(), new_nv, orig_mesh->GetNE(),
2238 0 : new_nbe, orig_mesh->SpaceDimension());
2239 :
2240 : // Copy vertices and domain and boundary elements.
2241 0 : for (int v = 0; v < orig_mesh->GetNV(); v++)
2242 : {
2243 0 : new_mesh->AddVertex(orig_mesh->GetVertex(v));
2244 : }
2245 0 : for (int e = 0; e < orig_mesh->GetNE(); e++)
2246 : {
2247 0 : mfem::Element *el = orig_mesh->GetElement(e)->Duplicate(new_mesh.get());
2248 0 : new_mesh->AddElement(el);
2249 : }
2250 0 : for (int be = 0; be < orig_mesh->GetNBE(); be++)
2251 : {
2252 0 : mfem::Element *bdr_el = orig_mesh->GetBdrElement(be)->Duplicate(new_mesh.get());
2253 0 : new_mesh->AddBdrElement(bdr_el);
2254 : }
2255 :
2256 : // Add duplicated vertices from interior boundary cracking, renumber the vertices of
2257 : // domain and boundary elements to tear the mesh, and add new crack boundary elements.
2258 0 : if (!crack_boundary_attributes.empty() && !crack_bdr_elem.empty())
2259 : {
2260 : // Add duplicate vertices. We assign the vertex number of the duplicated vertex in order
2261 : // to update the element connectivities in the next step.
2262 0 : for (auto &[orig_v, vert_components] : crack_vert_duplicates)
2263 : {
2264 0 : for (auto &[dup_v, component] : vert_components)
2265 : {
2266 0 : dup_v = new_mesh->AddVertex(orig_mesh->GetVertex(orig_v));
2267 : }
2268 : }
2269 :
2270 : // Renumber the duplicated vertex in the domain elements.
2271 0 : for (const auto &[orig_v, vert_components] : crack_vert_duplicates)
2272 : {
2273 0 : if (vert_components.empty())
2274 : {
2275 0 : continue; // Can skip vertices which were not duplicated
2276 : }
2277 0 : const int *elems = vert_to_elem->GetRow(orig_v);
2278 0 : for (int i = 0; i < vert_to_elem->RowSize(orig_v); i++)
2279 : {
2280 : // Find vertex in the element.
2281 0 : const auto e = elems[i];
2282 : mfem::Element *el = new_mesh->GetElement(e);
2283 0 : int *verts = el->GetVertices(), j;
2284 0 : for (j = 0; j < el->GetNVertices(); j++)
2285 : {
2286 0 : if (verts[j] == orig_v)
2287 : {
2288 : break;
2289 : }
2290 : }
2291 0 : MFEM_VERIFY(j < el->GetNVertices(), "Unable to locate vertex in element!");
2292 :
2293 : // Find the correct duplicate for this vertex. It's OK if the element is not in
2294 : // any of the connected components, this indicates that it keeps the original
2295 : // vertex and its connectivity is unmodified.
2296 0 : for (const auto &[dup_v, component] : vert_components)
2297 : {
2298 0 : if (component.find(e) != component.end())
2299 : {
2300 0 : verts[j] = dup_v;
2301 0 : break;
2302 : }
2303 : }
2304 : }
2305 : }
2306 :
2307 : // Finally, we insert the new duplicate boundary elements for the crack interface and
2308 : // also renumber the original boundary elements. To renumber the original boundary
2309 : // elements in the mesh, we use the updated vertex connectivity from the torn elements
2310 : // in the new mesh (done above).
2311 0 : const mfem::Table &elem_to_face = orig_mesh->ElementToFaceTable();
2312 0 : for (int be = 0; be < orig_mesh->GetNBE(); be++)
2313 : {
2314 : // Whether on the crack or not, we renumber the boundary element vertices as needed
2315 : // based on the neighboring element. For non-crack boundary elements, both
2316 : // neighboring elements must be part of the same connected component. First we find
2317 : // the index of the face in the old element, which should match the new element.
2318 : int f, o, e1, e2;
2319 0 : orig_mesh->GetBdrElementFace(be, &f, &o);
2320 0 : orig_mesh->GetFaceElements(f, &e1, &e2);
2321 0 : MFEM_VERIFY(e1 >= 0, "Boundary element with no attached elements!");
2322 : const int *faces = elem_to_face.GetRow(e1);
2323 : int i;
2324 0 : for (i = 0; i < elem_to_face.RowSize(e1); i++)
2325 : {
2326 0 : if (faces[i] == f)
2327 : {
2328 : break;
2329 : }
2330 : }
2331 0 : MFEM_VERIFY(i < elem_to_face.RowSize(e1), "Unable to locate face in element!");
2332 :
2333 : // Update the boundary element vertices.
2334 : mfem::Element *bdr_el = new_mesh->GetBdrElement(be);
2335 : const mfem::Element *el = new_mesh->GetElement(e1);
2336 0 : for (int j = 0; j < bdr_el->GetNVertices(); j++)
2337 : {
2338 0 : bdr_el->GetVertices()[j] = el->GetVertices()[el->GetFaceVertices(i)[j]];
2339 : }
2340 :
2341 : // Add the duplicate boundary element for boundary elements on the crack.
2342 0 : if (crack_bdr_elem.find(be) != crack_bdr_elem.end())
2343 : {
2344 0 : faces = elem_to_face.GetRow(e2);
2345 0 : for (i = 0; i < elem_to_face.RowSize(e2); i++)
2346 : {
2347 0 : if (faces[i] == f)
2348 : {
2349 : break;
2350 : }
2351 : }
2352 0 : MFEM_VERIFY(i < elem_to_face.RowSize(e2), "Unable to locate face in element!");
2353 :
2354 : // Add the interface boundary element attached to element 2 (the other part of the
2355 : // pair has been attached to element 1 in the previous step).
2356 0 : bdr_el = bdr_el->Duplicate(new_mesh.get());
2357 0 : el = new_mesh->GetElement(e2);
2358 0 : for (int j = 0; j < bdr_el->GetNVertices(); j++)
2359 : {
2360 0 : bdr_el->GetVertices()[j] = el->GetVertices()[el->GetFaceVertices(i)[j]];
2361 : }
2362 0 : new_mesh->AddBdrElement(bdr_el);
2363 : }
2364 : }
2365 : }
2366 :
2367 : // Add new boundary elements.
2368 0 : if (iodata.model.add_bdr_elements && !new_face_bdr_elem.empty())
2369 : {
2370 : // Some (1-based) boundary attributes may be empty since they were removed from the
2371 : // original mesh, but to keep attributes the same as config file we don't compress the
2372 : // list.
2373 0 : const mfem::Table &elem_to_face = orig_mesh->ElementToFaceTable();
2374 : int bdr_attr_max =
2375 0 : orig_mesh->bdr_attributes.Size() ? orig_mesh->bdr_attributes.Max() : 0;
2376 0 : for (int f = 0; f < orig_mesh->GetNumFaces(); f++)
2377 : {
2378 0 : if (new_face_bdr_elem[f] > 0)
2379 : {
2380 : // Assign new unique attribute based on attached elements. Save so that the
2381 : // attributes of e1 and e2 can be easily referenced using the new attribute. Since
2382 : // attributes are in 1-based indexing, a, b > 0. See also
2383 : // https://en.wikipedia.org/wiki/Pairing_function.
2384 : int e1, e2, a = 0, b = 0;
2385 0 : orig_mesh->GetFaceElements(f, &e1, &e2);
2386 0 : if (e1 >= 0 && e2 >= 0)
2387 : {
2388 : a = std::max(orig_mesh->GetAttribute(e1), orig_mesh->GetAttribute(e2));
2389 0 : b = (a == orig_mesh->GetAttribute(e1)) ? orig_mesh->GetAttribute(e2)
2390 : : orig_mesh->GetAttribute(e1);
2391 : }
2392 : else // e1 >= 0
2393 : {
2394 : a = orig_mesh->GetAttribute(e1);
2395 : b = 0;
2396 : }
2397 0 : MFEM_VERIFY(a + b > 0, "Invalid new boundary element attribute!");
2398 0 : long long int l_new_attr =
2399 0 : bdr_attr_max + (((a + b) * (long long int)(a + b + 1)) / 2) + a;
2400 : int new_attr = mfem::internal::to_int(l_new_attr); // At least bdr_attr_max + 1
2401 :
2402 : // Add the boundary elements with the new boundary attribute. The element vertices
2403 : // may have been renumbered in the new mesh, so the new face is not necessarily
2404 : // just a duplicate of the old one. First we find the index of the face in the old
2405 : // element, which should match the new element.
2406 : const int *faces = elem_to_face.GetRow(e1);
2407 : int i;
2408 0 : for (i = 0; i < elem_to_face.RowSize(e1); i++)
2409 : {
2410 0 : if (faces[i] == f)
2411 : {
2412 : break;
2413 : }
2414 : }
2415 0 : MFEM_VERIFY(i < elem_to_face.RowSize(e1), "Unable to locate face in element!");
2416 :
2417 : // Now add the boundary element(s).
2418 0 : mfem::Element *bdr_el = orig_mesh->GetFace(f)->Duplicate(new_mesh.get());
2419 : bdr_el->SetAttribute(new_attr);
2420 0 : const mfem::Element *el = new_mesh->GetElement(e1);
2421 0 : for (int j = 0; j < bdr_el->GetNVertices(); j++)
2422 : {
2423 0 : bdr_el->GetVertices()[j] = el->GetVertices()[el->GetFaceVertices(i)[j]];
2424 : }
2425 0 : new_mesh->AddBdrElement(bdr_el);
2426 : if constexpr (false)
2427 : {
2428 : Mpi::Print(
2429 : "Adding boundary element with attribute {:d} from elements {:d} and {:d}\n",
2430 : new_attr, a, b);
2431 : }
2432 0 : if (new_face_bdr_elem[f] > 1)
2433 : {
2434 : // Flip order of vertices to reverse normal direction of the second added element.
2435 0 : bdr_el = bdr_el->Duplicate(new_mesh.get());
2436 0 : std::reverse(bdr_el->GetVertices(),
2437 0 : bdr_el->GetVertices() + bdr_el->GetNVertices());
2438 0 : new_mesh->AddBdrElement(bdr_el);
2439 : if constexpr (false)
2440 : {
2441 : Mpi::Print("Adding second boundary element with attribute {:d} from elements "
2442 : "{:d} and {:d}\n",
2443 : new_attr, a, b);
2444 : }
2445 : }
2446 : }
2447 : }
2448 : }
2449 :
2450 : // Finalize the new mesh topology, and copy mesh curvature information if needed. This
2451 : // copies the nodes over correctly accounting for the element topology changes (the number
2452 : // of elements in the mesh has not changed, just their connectivity has).
2453 : constexpr bool generate_bdr = false;
2454 0 : new_mesh->FinalizeTopology(generate_bdr);
2455 0 : if (orig_mesh->GetNodes())
2456 : {
2457 0 : TransferHighOrderNodes(*orig_mesh, *new_mesh);
2458 : }
2459 :
2460 : // If we have added cracks for interior boundary elements, apply a very very small
2461 : // perturbation to separate the duplicated boundary elements on either side and prevent
2462 : // them from lying exactly on top of each other. This is mostly just for visualization
2463 : // and can be increased in magnitude for debugging.
2464 0 : if (!crack_boundary_attributes.empty() && !crack_bdr_elem.empty() &&
2465 0 : iodata.model.crack_displ_factor > 0.0)
2466 : {
2467 : // mfem::Mesh::MoveNodes expects byNODES ordering when using vertices.
2468 : mfem::GridFunction *nodes = new_mesh->GetNodes();
2469 : mfem::Ordering::Type ordering =
2470 0 : nodes ? nodes->FESpace()->GetOrdering() : mfem::Ordering::byNODES;
2471 0 : int sdim = new_mesh->SpaceDimension();
2472 0 : int nv = nodes ? nodes->Size() / sdim : new_mesh->GetNV();
2473 : auto Index = [ordering, sdim, nv](int v, int d)
2474 0 : { return (ordering == mfem::Ordering::byVDIM) ? sdim * v + d : d * nv + v; };
2475 : mfem::Vector normal(sdim);
2476 0 : mfem::IsoparametricTransformation T;
2477 : mfem::Array<int> dofs;
2478 :
2479 : // Compute the displacement as the average normal of the attached boundary elements.
2480 0 : mfem::Vector displacements(nv * sdim);
2481 0 : displacements = 0.0;
2482 0 : double h_min = mfem::infinity();
2483 0 : const mfem::Table &elem_to_face = orig_mesh->ElementToFaceTable();
2484 0 : const mfem::Table &new_elem_to_face = new_mesh->ElementToFaceTable();
2485 0 : for (auto be : crack_bdr_elem)
2486 : {
2487 : // Get the neighboring elements (same indices in the old and new mesh).
2488 : int f, o, e1, e2;
2489 0 : orig_mesh->GetBdrElementFace(be, &f, &o);
2490 0 : orig_mesh->GetFaceElements(f, &e1, &e2);
2491 :
2492 : // Perturb both new boundary elements in opposite directions.
2493 0 : for (auto e : {e1, e2})
2494 : {
2495 : // Find the index of the face in the old element, which matches the new element, so
2496 : // we can get the list of all vertices or nodes to perturb.
2497 : const int *faces = elem_to_face.GetRow(e);
2498 : int i;
2499 0 : for (i = 0; i < elem_to_face.RowSize(e); i++)
2500 : {
2501 0 : if (faces[i] == f)
2502 : {
2503 : break;
2504 : }
2505 : }
2506 0 : MFEM_VERIFY(i < elem_to_face.RowSize(e), "Unable to locate face in element!");
2507 :
2508 : // Compute the element normal, oriented to point outward from element 1 initially.
2509 0 : int new_f = new_elem_to_face.GetRow(e)[i];
2510 0 : if (e == e1)
2511 : {
2512 0 : new_mesh->GetFaceTransformation(new_f, &T);
2513 : const mfem::IntegrationPoint &ip =
2514 : mfem::Geometries.GetCenter(T.GetGeometryType());
2515 : T.SetIntPoint(&ip);
2516 0 : mfem::CalcOrtho(T.Jacobian(), normal);
2517 0 : double s = normal.Norml2();
2518 0 : h_min = std::min(h_min, std::sqrt(s));
2519 0 : normal /= -s; // We could also area-weight the average normal
2520 : }
2521 : else // e == e2
2522 : {
2523 0 : normal *= -1.0;
2524 : }
2525 :
2526 : // For all "nodes" associated with this crack face, update the direction of their
2527 : // displacements.
2528 0 : auto NodeUpdate = [&](int v)
2529 : {
2530 0 : for (int d = 0; d < sdim; d++)
2531 : {
2532 0 : const int idx = Index(v, d);
2533 0 : displacements(idx) += normal(d);
2534 : }
2535 0 : };
2536 0 : if (nodes)
2537 : {
2538 0 : nodes->FESpace()->GetFaceDofs(new_f, dofs);
2539 0 : for (int j = 0; j < dofs.Size(); j++)
2540 : {
2541 0 : NodeUpdate(dofs[j]);
2542 : }
2543 : }
2544 : else
2545 : {
2546 : const mfem::Element *el = new_mesh->GetElement(e);
2547 0 : for (int j = 0; j < el->GetNFaceVertices(i); j++)
2548 : {
2549 0 : NodeUpdate(el->GetVertices()[el->GetFaceVertices(i)[j]]);
2550 : }
2551 : }
2552 : }
2553 : }
2554 0 : for (int v = 0; v < nv; v++)
2555 : {
2556 : double s = 0.0;
2557 0 : for (int d = 0; d < sdim; d++)
2558 : {
2559 : const int idx = Index(v, d);
2560 0 : s += displacements(idx) * displacements(idx);
2561 : }
2562 0 : if (s > 0.0)
2563 : {
2564 0 : s = std::sqrt(s);
2565 0 : for (int d = 0; d < sdim; d++)
2566 : {
2567 : const int idx = Index(v, d);
2568 0 : displacements(idx) /= s;
2569 : }
2570 : }
2571 : }
2572 :
2573 : // Scale and apply the displacements. We don't need to do anything special to constrain
2574 : // the displacements at seam vertices (and associated high-order nodes on seam edges) to
2575 : // to zero, because the normals from both sides will average out to zero.
2576 0 : displacements *= (iodata.model.crack_displ_factor * h_min /
2577 0 : (nodes ? nodes->FESpace()->GetMaxElementOrder() : 1));
2578 0 : new_mesh->MoveNodes(displacements);
2579 0 : }
2580 :
2581 : orig_mesh = std::move(new_mesh);
2582 : return 1; // Success
2583 8 : }
2584 :
2585 8 : std::unique_ptr<int[]> GetMeshPartitioning(const mfem::Mesh &mesh, int size,
2586 : const std::string &part_file, bool print)
2587 : {
2588 8 : MFEM_VERIFY(size <= mesh.GetNE(), "Mesh partitioning must have parts <= mesh elements ("
2589 : << size << " vs. " << mesh.GetNE() << ")!");
2590 8 : if (part_file.length() == 0)
2591 : {
2592 : const int part_method = 1;
2593 : std::unique_ptr<int[]> partitioning(
2594 8 : const_cast<mfem::Mesh &>(mesh).GeneratePartitioning(size, part_method));
2595 8 : if (print)
2596 : {
2597 8 : Mpi::Print("Finished partitioning mesh into {:d} subdomain{}\n", size,
2598 24 : (size > 1) ? "s" : "");
2599 : }
2600 : return partitioning;
2601 : }
2602 : // User can optionally specify a mesh partitioning file as generated from the MFEM
2603 : // mesh-explorer miniapp, for example. It has the format:
2604 : //
2605 : // number_of_elements <NE>
2606 : // number_of_processors <NPART>
2607 : // <part[0]>
2608 : // ...
2609 : // <part[NE-1]>
2610 : //
2611 : int ne, np;
2612 0 : std::ifstream part_ifs(part_file);
2613 0 : part_ifs.ignore(std::numeric_limits<std::streamsize>::max(), ' ');
2614 0 : part_ifs >> ne;
2615 0 : if (ne != mesh.GetNE())
2616 : {
2617 0 : MFEM_ABORT("Invalid partitioning file (number of elements)!");
2618 : }
2619 0 : part_ifs.ignore(std::numeric_limits<std::streamsize>::max(), ' ');
2620 0 : part_ifs >> np;
2621 0 : if (np != size)
2622 : {
2623 0 : MFEM_ABORT("Invalid partitioning file (number of processors)!");
2624 : }
2625 0 : auto partitioning = std::make_unique<int[]>(mesh.GetNE());
2626 : int i = 0;
2627 0 : while (i < mesh.GetNE())
2628 : {
2629 0 : part_ifs >> partitioning[i++];
2630 : }
2631 0 : if (print)
2632 : {
2633 0 : Mpi::Print("Read mesh partitioning into {:d} subdomain{} from disk\n", size,
2634 0 : (size > 1) ? "s" : "");
2635 : }
2636 : return partitioning;
2637 0 : }
2638 :
2639 8 : std::unique_ptr<mfem::ParMesh> DistributeMesh(MPI_Comm comm,
2640 : std::unique_ptr<mfem::Mesh> &smesh,
2641 : const int *partitioning,
2642 : const std::string &output_dir)
2643 : {
2644 : // Take a serial mesh and partitioning on the root process and construct the global
2645 : // parallel mesh. For now, prefer the MPI-based version to the file IO one. When
2646 : // constructing the ParMesh, we mark for refinement since refinement flags are not copied
2647 : // from the serial mesh. Beware that mfem::ParMesh constructor argument order is not the
2648 : // same as mfem::Mesh! Each processor's component gets sent as a byte string.
2649 8 : constexpr bool generate_edges = false, refine = true, fix_orientation = false;
2650 8 : std::unique_ptr<mfem::ParMesh> pmesh;
2651 8 : if (Mpi::Root(comm))
2652 : {
2653 : mfem::MeshPartitioner partitioner(*smesh, Mpi::Size(comm),
2654 8 : const_cast<int *>(partitioning));
2655 8 : std::vector<MPI_Request> send_requests(Mpi::Size(comm) - 1, MPI_REQUEST_NULL);
2656 : std::vector<std::string> so;
2657 8 : so.reserve(Mpi::Size(comm));
2658 24 : for (int i = 0; i < Mpi::Size(comm); i++)
2659 : {
2660 8 : mfem::MeshPart part;
2661 8 : partitioner.ExtractPart(i, part);
2662 8 : std::ostringstream fo(std::stringstream::out);
2663 : // fo << std::fixed;
2664 : fo << std::scientific;
2665 : fo.precision(MSH_FLT_PRECISION);
2666 8 : part.Print(fo);
2667 8 : so.push_back(fo.str());
2668 : // so.push_back((i > 0) ? zlib::CompressString(fo.str()) : fo.str());
2669 8 : if (i > 0)
2670 : {
2671 0 : int slen = static_cast<int>(so[i].length());
2672 0 : MFEM_VERIFY(so[i].length() == (std::size_t)slen,
2673 : "Overflow error distributing parallel mesh!");
2674 0 : MPI_Isend(so[i].data(), slen, MPI_CHAR, i, i, comm, &send_requests[i - 1]);
2675 : }
2676 8 : }
2677 : smesh.reset();
2678 8 : std::istringstream fi(so[0]); // This is never compressed
2679 : pmesh =
2680 16 : std::make_unique<mfem::ParMesh>(comm, fi, refine, generate_edges, fix_orientation);
2681 8 : MPI_Waitall(static_cast<int>(send_requests.size()), send_requests.data(),
2682 : MPI_STATUSES_IGNORE);
2683 16 : }
2684 : else
2685 : {
2686 : MPI_Status status;
2687 : int rlen;
2688 : std::string si;
2689 0 : MPI_Probe(0, Mpi::Rank(comm), comm, &status);
2690 0 : MPI_Get_count(&status, MPI_CHAR, &rlen);
2691 0 : si.resize(rlen);
2692 0 : MPI_Recv(si.data(), rlen, MPI_CHAR, 0, Mpi::Rank(comm), comm, MPI_STATUS_IGNORE);
2693 0 : std::istringstream fi(si);
2694 : // std::istringstream fi(zlib::DecompressString(si));
2695 : pmesh =
2696 0 : std::make_unique<mfem::ParMesh>(comm, fi, refine, generate_edges, fix_orientation);
2697 0 : }
2698 8 : return pmesh;
2699 : }
2700 :
2701 0 : void RebalanceConformalMesh(std::unique_ptr<mfem::ParMesh> &pmesh)
2702 : {
2703 : // Write the parallel mesh to a stream as a serial mesh, then read back in and partition
2704 : // using METIS.
2705 : MPI_Comm comm = pmesh->GetComm();
2706 : constexpr bool generate_edges = false, generate_bdr = false, refine = true,
2707 : fix_orientation = false;
2708 0 : std::unique_ptr<mfem::Mesh> smesh;
2709 : if constexpr (false)
2710 : {
2711 : // Write the serial mesh to a stream and read that through the Mesh constructor.
2712 : std::ostringstream fo(std::stringstream::out);
2713 : // fo << std::fixed;
2714 : fo << std::scientific;
2715 : fo.precision(MSH_FLT_PRECISION);
2716 : pmesh->PrintAsSerial(fo);
2717 : if (Mpi::Root(comm))
2718 : {
2719 : smesh = std::make_unique<mfem::Mesh>(fo, generate_edges, refine, fix_orientation);
2720 : }
2721 : }
2722 : else
2723 : {
2724 : // Directly ingest the generated Mesh and release the no longer needed memory.
2725 0 : smesh = std::make_unique<mfem::Mesh>(pmesh->GetSerialMesh(0));
2726 0 : if (Mpi::Root(comm))
2727 : {
2728 0 : smesh->FinalizeTopology(generate_bdr);
2729 0 : smesh->Finalize(refine, fix_orientation);
2730 : }
2731 : else
2732 : {
2733 : smesh.reset();
2734 : }
2735 : }
2736 :
2737 : // (Re)-construct the parallel mesh.
2738 : std::unique_ptr<int[]> partitioning;
2739 0 : if (Mpi::Root(comm))
2740 : {
2741 0 : partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), "", false);
2742 : }
2743 0 : pmesh = DistributeMesh(comm, smesh, partitioning.get());
2744 0 : }
2745 :
2746 : } // namespace
2747 :
2748 : } // namespace palace
|
