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 :
6 : #include "geodata_impl.hpp"
7 :
8 : #include <array>
9 : #include <limits>
10 : #include <random>
11 : #include <unordered_set>
12 : #include <Eigen/Dense>
13 : #include "utils/communication.hpp"
14 : #include "utils/omp.hpp"
15 :
16 : namespace palace::mesh
17 : {
18 :
19 : using Vector3dMap = Eigen::Map<Eigen::Vector3d>;
20 : using CVector3dMap = Eigen::Map<const Eigen::Vector3d>;
21 :
22 : // Compute a lexicographic comparison of Eigen Vector3d.
23 7176 : bool EigenLE(const Eigen::Vector3d &x, const Eigen::Vector3d &y)
24 : {
25 7176 : return std::lexicographical_compare(x.begin(), x.end(), y.begin(), y.end());
26 : }
27 :
28 : // Helper for collecting a point cloud from a mesh, used in calculating bounding boxes and
29 : // bounding balls. Returns the dominant rank, for which the vertices argument will be
30 : // filled, while all other ranks will have an empty vector. Vertices are de-duplicated to a
31 : // certain floating point precision.
32 60 : int CollectPointCloudOnRoot(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
33 : bool bdr, std::vector<Eigen::Vector3d> &vertices)
34 : {
35 60 : if (!mesh.GetNodes())
36 : {
37 : // Linear mesh, work with element vertices directly.
38 0 : PalacePragmaOmp(parallel)
39 : {
40 : std::unordered_set<int> vertex_indices;
41 : if (bdr)
42 : {
43 : PalacePragmaOmp(for schedule(static))
44 : for (int i = 0; i < mesh.GetNBE(); i++)
45 : {
46 : if (!marker[mesh.GetBdrAttribute(i) - 1])
47 : {
48 : continue;
49 : }
50 : const int *verts = mesh.GetBdrElement(i)->GetVertices();
51 : vertex_indices.insert(verts, verts + mesh.GetBdrElement(i)->GetNVertices());
52 : }
53 : }
54 : else
55 : {
56 : PalacePragmaOmp(for schedule(static))
57 : for (int i = 0; i < mesh.GetNE(); i++)
58 : {
59 : if (!marker[mesh.GetAttribute(i) - 1])
60 : {
61 : continue;
62 : }
63 : const int *verts = mesh.GetElement(i)->GetVertices();
64 : vertex_indices.insert(verts, verts + mesh.GetElement(i)->GetNVertices());
65 : }
66 : }
67 : PalacePragmaOmp(critical(PointCloud))
68 : {
69 : for (auto i : vertex_indices)
70 : {
71 : const auto &vx = mesh.GetVertex(i);
72 : vertices.emplace_back(vx[0], vx[1], vx[2]);
73 : }
74 : }
75 : }
76 : }
77 : else
78 : {
79 : // Curved mesh, need to process point matrices.
80 60 : const int ref = mesh.GetNodes()->FESpace()->GetMaxElementOrder();
81 432 : auto AddPoints = [&](mfem::GeometryRefiner &refiner, mfem::ElementTransformation &T,
82 : mfem::DenseMatrix &pointmat,
83 : std::vector<Eigen::Vector3d> &loc_vertices)
84 : {
85 432 : mfem::RefinedGeometry *RefG = refiner.Refine(T.GetGeometryType(), ref);
86 432 : T.Transform(RefG->RefPts, pointmat);
87 1728 : for (int j = 0; j < pointmat.Width(); j++)
88 : {
89 1296 : loc_vertices.emplace_back(pointmat(0, j), pointmat(1, j), pointmat(2, j));
90 : }
91 492 : };
92 60 : PalacePragmaOmp(parallel)
93 : {
94 : mfem::GeometryRefiner refiner;
95 : mfem::IsoparametricTransformation T;
96 : mfem::DenseMatrix pointmat; // 3 x N
97 : std::vector<Eigen::Vector3d> loc_vertices;
98 : if (bdr)
99 : {
100 : PalacePragmaOmp(for schedule(static))
101 : for (int i = 0; i < mesh.GetNBE(); i++)
102 : {
103 : if (!marker[mesh.GetBdrAttribute(i) - 1])
104 : {
105 : continue;
106 : }
107 : mesh.GetBdrElementTransformation(i, &T);
108 : AddPoints(refiner, T, pointmat, loc_vertices);
109 : }
110 : }
111 : else
112 : {
113 : PalacePragmaOmp(for schedule(static))
114 : for (int i = 0; i < mesh.GetNE(); i++)
115 : {
116 : if (!marker[mesh.GetAttribute(i) - 1])
117 : {
118 : continue;
119 : }
120 : mesh.GetElementTransformation(i, &T);
121 : AddPoints(refiner, T, pointmat, loc_vertices);
122 : }
123 : }
124 : PalacePragmaOmp(critical(PointCloud))
125 : {
126 : for (const auto &v : loc_vertices)
127 : {
128 : vertices.push_back(v);
129 : }
130 : }
131 : }
132 : }
133 :
134 : // dominant_rank will perform the calculation.
135 60 : MPI_Comm comm = mesh.GetComm();
136 60 : const auto num_vertices = int(vertices.size());
137 180 : const int dominant_rank = [&]()
138 : {
139 60 : int vert = num_vertices, rank = Mpi::Rank(comm);
140 60 : Mpi::GlobalMaxLoc(1, &vert, &rank, comm);
141 60 : return rank;
142 60 : }();
143 60 : std::vector<int> recv_counts(Mpi::Size(comm)), displacements;
144 : std::vector<Eigen::Vector3d> collected_vertices;
145 60 : MPI_Gather(&num_vertices, 1, MPI_INT, recv_counts.data(), 1, MPI_INT, dominant_rank,
146 : comm);
147 60 : if (dominant_rank == Mpi::Rank(comm))
148 : {
149 : // First displacement is zero, then after is the partial sum of recv_counts.
150 56 : displacements.resize(Mpi::Size(comm));
151 56 : displacements[0] = 0;
152 : std::partial_sum(recv_counts.begin(), recv_counts.end() - 1, displacements.begin() + 1);
153 :
154 : // Add on slots at the end of vertices for the incoming data.
155 56 : collected_vertices.resize(std::accumulate(recv_counts.begin(), recv_counts.end(), 0));
156 :
157 : // MPI transfer will be done with MPI_DOUBLE, so duplicate all these values.
158 116 : for (auto &x : displacements)
159 : {
160 60 : x *= 3;
161 : }
162 116 : for (auto &x : recv_counts)
163 : {
164 60 : x *= 3;
165 : }
166 : }
167 :
168 : // Gather the data to the dominant rank.
169 : static_assert(sizeof(Eigen::Vector3d) == 3 * sizeof(double));
170 60 : MPI_Gatherv(vertices.data(), 3 * num_vertices, MPI_DOUBLE, collected_vertices.data(),
171 : recv_counts.data(), displacements.data(), MPI_DOUBLE, dominant_rank, comm);
172 :
173 : // Deduplicate vertices. Given floating point precision, need a tolerance.
174 60 : if (dominant_rank == Mpi::Rank(comm))
175 : {
176 : auto vertex_equality = [](const auto &x, const auto &y)
177 : {
178 : constexpr double tolerance = 10.0 * std::numeric_limits<double>::epsilon();
179 1240 : return std::abs(x[0] - y[0]) < tolerance && std::abs(x[1] - y[1]) < tolerance &&
180 944 : std::abs(x[2] - y[2]) < tolerance;
181 : };
182 : vertices = std::move(collected_vertices);
183 56 : std::sort(vertices.begin(), vertices.end(), EigenLE);
184 56 : vertices.erase(std::unique(vertices.begin(), vertices.end(), vertex_equality),
185 : vertices.end());
186 : }
187 : else
188 : {
189 : vertices.clear();
190 : }
191 :
192 60 : return dominant_rank;
193 : }
194 :
195 : // Compute the distance from a point orthogonal to the list of normal axes, relative to
196 : // the given origin.
197 1365 : auto PerpendicularDistance(const std::initializer_list<Eigen::Vector3d> &normals,
198 : const Eigen::Vector3d &origin, const Eigen::Vector3d &v)
199 : {
200 : Eigen::Vector3d v0 = v - origin;
201 3463 : for (const auto &n : normals)
202 : {
203 : v0 -= n.dot(v0) * n;
204 : }
205 1365 : return v0.norm();
206 : }
207 :
208 : // Calculates a bounding box from a point cloud, result is broadcast across all processes.
209 31 : BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm,
210 : const std::vector<Eigen::Vector3d> &vertices,
211 : int dominant_rank)
212 : {
213 : BoundingBox box;
214 31 : if (dominant_rank == Mpi::Rank(comm))
215 : {
216 : // Pick a candidate 000 vertex using lexicographic sort. This can be vulnerable to
217 : // floating point precision if the box is axis aligned, but not floating point exact.
218 : // Pick candidate 111 as the furthest from this candidate, then reassign 000 as the
219 : // furthest from 111. Such a pair has to form the diagonal for a point cloud defining a
220 : // box. Verify that p_111 is also the maximum distance from p_000 -> a diagonal is
221 : // found.
222 29 : MFEM_VERIFY(vertices.size() >= 4,
223 : "A bounding box requires a minimum of four vertices for this algorithm!");
224 : auto p_000 = std::min_element(vertices.begin(), vertices.end(), EigenLE);
225 : auto p_111 =
226 29 : std::max_element(vertices.begin(), vertices.end(),
227 316 : [p_000](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
228 316 : { return (x - *p_000).norm() < (y - *p_000).norm(); });
229 29 : p_000 = std::max_element(vertices.begin(), vertices.end(),
230 316 : [p_111](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
231 316 : { return (x - *p_111).norm() < (y - *p_111).norm(); });
232 : MFEM_ASSERT(std::max_element(vertices.begin(), vertices.end(),
233 : [p_000](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
234 : { return (x - *p_000).norm() < (y - *p_000).norm(); }) ==
235 : p_111,
236 : "p_000 and p_111 must be mutually opposing points!");
237 :
238 : // Define a diagonal of the ASSUMED cuboid bounding box.
239 : const auto &v_000 = *p_000;
240 : const auto &v_111 = *p_111;
241 29 : MFEM_VERIFY(&v_000 != &v_111, "Minimum and maximum extents cannot be identical!");
242 : const Eigen::Vector3d origin = v_000;
243 29 : const Eigen::Vector3d n_1 = (v_111 - v_000).normalized();
244 :
245 : // Find the vertex furthest from the diagonal axis. We cannot know yet if this defines
246 : // (001) or (011).
247 29 : const auto &t_0 = *std::max_element(vertices.begin(), vertices.end(),
248 316 : [&](const auto &x, const auto &y)
249 : {
250 316 : return PerpendicularDistance({n_1}, origin, x) <
251 316 : PerpendicularDistance({n_1}, origin, y);
252 : });
253 29 : MFEM_VERIFY(&t_0 != &v_000 && &t_0 != &v_111, "Vertices are degenerate!");
254 :
255 : // Use the discovered vertex to define a second direction and thus a plane. n_1 and n_2
256 : // now define a planar coordinate system intersecting the main diagonal, and two
257 : // opposite edges of the cuboid.
258 : const Eigen::Vector3d n_2 =
259 29 : ((t_0 - origin) - ((t_0 - origin).dot(n_1) * n_1)).normalized();
260 :
261 : // Collect the furthest point from the plane to determine if the box is planar. Look for
262 : // a component that maximizes distance from the planar system: complete the axes with a
263 : // cross, then use a dot product to pick the greatest deviation.
264 29 : constexpr double rel_tol = 1.0e-6;
265 29 : auto max_distance = PerpendicularDistance(
266 : {n_1, n_2}, origin,
267 29 : *std::max_element(vertices.begin(), vertices.end(),
268 316 : [&](const auto &x, const auto &y)
269 : {
270 316 : return PerpendicularDistance({n_1, n_2}, origin, x) <
271 316 : PerpendicularDistance({n_1, n_2}, origin, y);
272 29 : }));
273 29 : box.planar = (max_distance < rel_tol * (v_111 - v_000).norm());
274 :
275 : // For the non-planar case, collect points furthest from the plane and choose the one
276 : // closest to the origin as the next vertex which might be (001) or (011).
277 29 : const auto &t_1 = [&]()
278 : {
279 29 : if (box.planar)
280 : {
281 21 : return t_0;
282 : }
283 : std::vector<Eigen::Vector3d> vertices_out_of_plane;
284 8 : std::copy_if(vertices.begin(), vertices.end(),
285 : std::back_inserter(vertices_out_of_plane),
286 72 : [&](const auto &v)
287 : {
288 72 : return std::abs(PerpendicularDistance({n_1, n_2}, origin, v) -
289 72 : max_distance) < rel_tol * max_distance;
290 : });
291 8 : return *std::min_element(vertices_out_of_plane.begin(), vertices_out_of_plane.end(),
292 0 : [&](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
293 0 : { return (x - origin).norm() < (y - origin).norm(); });
294 29 : }();
295 :
296 : // Given candidates t_0 and t_1, the closer to origin defines v_001.
297 : const bool t_0_gt_t_1 = (t_0 - origin).norm() > (t_1 - origin).norm();
298 29 : const auto &v_001 = t_0_gt_t_1 ? t_1 : t_0;
299 29 : const auto &v_011 = box.planar ? v_111 : (t_0_gt_t_1 ? t_0 : t_1);
300 :
301 : // Compute the center as halfway along the main diagonal.
302 29 : Vector3dMap(box.center.data()) = 0.5 * (v_000 + v_111);
303 :
304 : if constexpr (false)
305 : {
306 : fmt::print("box.center {}!\n", box.center);
307 : fmt::print("v_000 {}!\n", v_000);
308 : fmt::print("v_001 {}!\n", v_001);
309 : fmt::print("v_011 {}!\n", v_011);
310 : fmt::print("v_111 {}!\n", v_111);
311 : }
312 :
313 : // Compute the box axes. Using the 4 extremal points, we find the first two axes as the
314 : // edges which are closest to perpendicular. For a perfect rectangular prism point
315 : // cloud, we could instead compute the axes and length in each direction using the
316 : // found edges of the cuboid, but this does not work for non-rectangular prism
317 : // cross-sections or pyramid shapes.
318 : {
319 29 : const auto [e_0, e_1] = [&v_000, &v_001, &v_011, &v_111]()
320 : {
321 29 : std::array<const Eigen::Vector3d *, 4> verts = {&v_000, &v_001, &v_011, &v_111};
322 : Eigen::Vector3d e_0 = Eigen::Vector3d::Zero(), e_1 = Eigen::Vector3d::Zero();
323 : double dot_min = mfem::infinity();
324 33 : for (int i_0 = 0; i_0 < 4; i_0++)
325 : {
326 38 : for (int j_0 = i_0 + 1; j_0 < 4; j_0++)
327 : {
328 86 : for (int i_1 = 0; i_1 < 4; i_1++)
329 : {
330 200 : for (int j_1 = i_1 + 1; j_1 < 4; j_1++)
331 : {
332 148 : if ((i_1 == i_0 && j_1 == j_0) || verts[i_0] == verts[j_0] ||
333 109 : verts[i_1] == verts[j_1])
334 : {
335 44 : continue;
336 : }
337 104 : const auto e_ij_0 = (*verts[j_0] - *verts[i_0]).normalized();
338 104 : const auto e_ij_1 = (*verts[j_1] - *verts[i_1]).normalized();
339 : const auto dot = std::abs(e_ij_0.dot(e_ij_1));
340 104 : if (dot < dot_min)
341 : {
342 : if constexpr (false)
343 : {
344 : fmt::print("i_0 {} i_1 {} j_0 {} j_1 {}\n", i_0, i_1, j_0, j_1);
345 : fmt::print("e_ij_0 {}, e_ij_1 {}!\n", e_ij_0, e_ij_1);
346 : }
347 : dot_min = dot;
348 : e_0 = e_ij_0;
349 : e_1 = e_ij_1;
350 66 : if (dot_min < rel_tol)
351 : {
352 28 : return std::make_pair(e_0, e_1);
353 : }
354 : }
355 : }
356 : }
357 : }
358 : }
359 : return std::make_pair(e_0, e_1);
360 29 : }();
361 :
362 : if constexpr (false)
363 : {
364 : fmt::print("e_0 {}, e_1 {}!\n", e_0, e_1);
365 : }
366 :
367 : Vector3dMap(box.axes[0].data()) = e_0;
368 : Vector3dMap(box.axes[1].data()) = e_1;
369 : Vector3dMap(box.axes[2].data()) =
370 29 : box.planar ? Eigen::Vector3d::Zero() : e_0.cross(e_1);
371 : }
372 :
373 : // Scale axes by length of the box in each direction.
374 29 : std::array<double, 3> l = {0.0};
375 145 : for (const auto &v : {v_000, v_001, v_011, v_111})
376 : {
377 : const auto v_0 = v - Vector3dMap(box.center.data());
378 116 : l[0] = std::max(l[0], std::abs(v_0.dot(Vector3dMap(box.axes[0].data()))));
379 116 : l[1] = std::max(l[1], std::abs(v_0.dot(Vector3dMap(box.axes[1].data()))));
380 116 : l[2] = std::max(l[2], std::abs(v_0.dot(Vector3dMap(box.axes[2].data()))));
381 : }
382 : Vector3dMap(box.axes[0].data()) *= l[0];
383 : Vector3dMap(box.axes[1].data()) *= l[1];
384 : Vector3dMap(box.axes[2].data()) *= l[2];
385 :
386 : // Make sure the longest dimension comes first.
387 : std::sort(box.axes.begin(), box.axes.end(), [](const auto &x, const auto &y)
388 : { return CVector3dMap(x.data()).norm() > CVector3dMap(y.data()).norm(); });
389 : }
390 :
391 : // Broadcast result to all processors.
392 : Mpi::Broadcast(3, box.center.data(), dominant_rank, comm);
393 : Mpi::Broadcast(3 * 3, box.axes.data()->data(), dominant_rank, comm);
394 : Mpi::Broadcast(1, &box.planar, dominant_rank, comm);
395 :
396 31 : return box;
397 : }
398 :
399 : // Use 4 points to define a sphere in 3D. If the points are coplanar, 3 of them are used to
400 : // define a circle which is interpreted as the equator of the sphere. We assume the points
401 : // are unique and not collinear.
402 0 : BoundingBall SphereFromPoints(const std::vector<std::size_t> &indices,
403 : const std::vector<Eigen::Vector3d> &vertices)
404 : {
405 : // Given 0 or 1 points, just return a radius of 0.
406 0 : MFEM_VERIFY(
407 : indices.size() <= 4,
408 : "Determining a sphere in 3D requires 4 points (and a circle requires 3 points)!");
409 : BoundingBall ball;
410 0 : ball.planar = (indices.size() < 4);
411 0 : if (indices.size() < 2)
412 : {
413 : ball.origin = Eigen::Vector3d::Zero();
414 0 : ball.radius = 0.0;
415 0 : return ball;
416 : }
417 :
418 : // For two points, construct a circle with the segment as its diameter. This could also
419 : // handle the collinear case for more than 2 points.
420 0 : if (indices.size() == 2)
421 : {
422 0 : ball.origin = 0.5 * (vertices[indices[0]] + vertices[indices[1]]);
423 0 : ball.radius = (vertices[indices[0]] - ball.origin).norm();
424 0 : return ball;
425 : }
426 :
427 : // Check for coplanarity.
428 : constexpr double rel_tol = 1.0e-6;
429 0 : const Eigen::Vector3d AB = vertices[indices[1]] - vertices[indices[0]];
430 0 : const Eigen::Vector3d AC = vertices[indices[2]] - vertices[indices[0]];
431 : const Eigen::Vector3d ABAC = AB.cross(AC);
432 : Eigen::Vector3d AD = Eigen::Vector3d::Zero();
433 0 : if (!ball.planar)
434 : {
435 0 : AD = vertices[indices[3]] - vertices[indices[0]];
436 0 : ball.planar = (std::abs(AD.dot(ABAC)) < rel_tol * AD.norm() * ABAC.norm());
437 : }
438 :
439 : // Construct a circle passing through 3 points.
440 : // See: https://en.wikipedia.org/wiki/Circumcircle#Higher_dimensions.
441 0 : if (ball.planar)
442 : {
443 0 : ball.origin = (0.5 / ABAC.squaredNorm()) *
444 0 : ((AB.squaredNorm() * AC) - (AC.squaredNorm() * AB)).cross(ABAC);
445 0 : ball.radius = ball.origin.norm();
446 0 : ball.origin += vertices[indices[0]];
447 : #if defined(MFEM_DEBUG)
448 : const auto r1 = (vertices[indices[1]] - ball.origin).norm();
449 : const auto r2 = (vertices[indices[2]] - ball.origin).norm();
450 : MFEM_VERIFY((1.0 - rel_tol) * ball.radius < r1 && r1 < (1.0 + rel_tol) * ball.radius &&
451 : (1.0 - rel_tol) * ball.radius < r2 &&
452 : r2 < (1.0 + rel_tol) * ball.radius,
453 : "Invalid circle calculated from 3 points!");
454 : #endif
455 : return ball;
456 : }
457 :
458 : // Construct a sphere passing through 4 points.
459 : // See: https://steve.hollasch.net/cgindex/geometry/sphere4pts.html.
460 : Eigen::Matrix3d C;
461 : Eigen::Vector3d d;
462 0 : const auto s = vertices[indices[0]].squaredNorm();
463 : C.row(0) = AB.transpose();
464 : C.row(1) = AC.transpose();
465 : C.row(2) = AD.transpose();
466 0 : d(0) = 0.5 * (vertices[indices[1]].squaredNorm() - s);
467 0 : d(1) = 0.5 * (vertices[indices[2]].squaredNorm() - s);
468 0 : d(2) = 0.5 * (vertices[indices[3]].squaredNorm() - s);
469 0 : ball.origin = C.inverse() * d; // 3x3 matrix inverse might be faster than general LU
470 : // if Eigen uses the explicit closed-form solution
471 0 : ball.radius = (vertices[indices[0]] - ball.origin).norm();
472 : #if defined(MFEM_DEBUG)
473 : const auto r1 = (vertices[indices[1]] - ball.origin).norm();
474 : const auto r2 = (vertices[indices[2]] - ball.origin).norm();
475 : const auto r3 = (vertices[indices[3]] - ball.origin).norm();
476 : MFEM_VERIFY((1.0 - rel_tol) * ball.radius < r1 && r1 < (1.0 + rel_tol) * ball.radius &&
477 : (1.0 - rel_tol) * ball.radius < r2 &&
478 : r2 < (1.0 + rel_tol) * ball.radius &&
479 : (1.0 - rel_tol) * ball.radius < r3 && r3 < (1.0 + rel_tol) * ball.radius,
480 : "Invalid sphere calculated from 3 points!");
481 : #endif
482 0 : return ball;
483 : }
484 :
485 0 : BoundingBall Welzl(std::vector<std::size_t> P, std::vector<std::size_t> R,
486 : const std::vector<Eigen::Vector3d> &vertices)
487 : {
488 : // Base case.
489 0 : if (R.size() == 4 || P.empty())
490 : {
491 0 : return SphereFromPoints(R, vertices);
492 : }
493 :
494 : // Choose a p ∈ P randomly, and recurse for (P \ {p}, R). The set P has already been
495 : // randomized on input.
496 0 : const std::size_t p = P.back();
497 : P.pop_back();
498 0 : BoundingBall D = Welzl(P, R, vertices);
499 :
500 : // If p is outside the sphere, recurse for (P \ {p}, R U {p}).
501 : constexpr double rel_tol = 1.0e-6;
502 0 : if ((vertices[p] - D.origin).norm() >= (1.0 + rel_tol) * D.radius)
503 : {
504 : R.push_back(p);
505 0 : D = Welzl(P, R, vertices);
506 : }
507 :
508 : return D;
509 : }
510 :
511 : // Calculates a bounding ball from a point cloud using Welzl's algorithm, result is
512 : // broadcast across all processes. We don't operate on the convex hull, since the number of
513 : // points should be small enough that operating on the full set should be OK. If only three
514 : // points are provided, the bounding circle is computed (likewise for if the points are
515 : // coplanar).
516 0 : BoundingBox BoundingBallFromPointCloud(MPI_Comm comm,
517 : const std::vector<Eigen::Vector3d> &vertices,
518 : int dominant_rank)
519 : {
520 : BoundingBox ball;
521 0 : if (dominant_rank == Mpi::Rank(comm))
522 : {
523 0 : MFEM_VERIFY(vertices.size() >= 3,
524 : "A bounding ball requires a minimum of three vertices for this algorithm!");
525 0 : std::vector<std::size_t> indices(vertices.size());
526 : std::iota(indices.begin(), indices.end(), 0);
527 :
528 : // Acceleration from https://informatica.vu.lt/journal/INFORMATICA/article/1251. Allow
529 : // for duplicate points and just add the 4 points to the end of the indicies list to be
530 : // considered first. The two points are not necessarily the maximizer of the distance
531 : // between all pairs, but they should be a good estimate.
532 : {
533 0 : auto p_1 = std::min_element(vertices.begin(), vertices.end(), EigenLE);
534 0 : auto p_2 = std::max_element(vertices.begin(), vertices.end(),
535 0 : [p_1](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
536 0 : { return (x - *p_1).norm() < (y - *p_1).norm(); });
537 0 : p_1 = std::max_element(vertices.begin(), vertices.end(),
538 0 : [p_2](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
539 0 : { return (x - *p_2).norm() < (y - *p_2).norm(); });
540 :
541 : // Find the next point as the vertex furthest from the initial axis.
542 0 : const Eigen::Vector3d n_1 = (*p_2 - *p_1).normalized();
543 0 : auto p_3 = std::max_element(vertices.begin(), vertices.end(),
544 0 : [&](const auto &x, const auto &y)
545 : {
546 0 : return PerpendicularDistance({n_1}, *p_1, x) <
547 0 : PerpendicularDistance({n_1}, *p_1, y);
548 : });
549 0 : auto p_4 = std::max_element(vertices.begin(), vertices.end(),
550 0 : [p_3](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
551 0 : { return (x - *p_3).norm() < (y - *p_3).norm(); });
552 0 : MFEM_VERIFY(p_3 != p_1 && p_3 != p_2 && p_4 != p_1 && p_4 != p_2,
553 : "Vertices are degenerate!");
554 :
555 : // Start search with these points, which should be roughly extremal. With the search
556 : // for p_3 done in an orthogonal direction, p_1, p_2, p_3, and p_4 should all be
557 : // unique.
558 0 : std::swap(indices[indices.size() - 1], indices[p_1 - vertices.begin()]);
559 0 : std::swap(indices[indices.size() - 2], indices[p_2 - vertices.begin()]);
560 0 : std::swap(indices[indices.size() - 3], indices[p_3 - vertices.begin()]);
561 0 : std::swap(indices[indices.size() - 4], indices[p_4 - vertices.begin()]);
562 : }
563 :
564 : // Randomly permute the point set.
565 : {
566 0 : std::random_device rd;
567 0 : std::mt19937 g(rd());
568 0 : std::shuffle(indices.begin(), indices.end() - 4, g);
569 : }
570 :
571 : // Compute the bounding ball.
572 0 : BoundingBall min_ball = Welzl(indices, {}, vertices);
573 : Vector3dMap(ball.center.data()) = min_ball.origin;
574 0 : Vector3dMap(ball.axes[0].data()) = Eigen::Vector3d(min_ball.radius, 0.0, 0.0);
575 0 : Vector3dMap(ball.axes[1].data()) = Eigen::Vector3d(0.0, min_ball.radius, 0.0);
576 : Vector3dMap(ball.axes[2].data()) = Eigen::Vector3d(0.0, 0.0, min_ball.radius);
577 0 : ball.planar = min_ball.planar;
578 : }
579 :
580 : // Broadcast result to all processors.
581 : Mpi::Broadcast(3, ball.center.data(), dominant_rank, comm);
582 : Mpi::Broadcast(3 * 3, ball.axes.data()->data(), dominant_rank, comm);
583 : Mpi::Broadcast(1, &ball.planar, dominant_rank, comm);
584 :
585 0 : return ball;
586 : }
587 :
588 : // Compute a normal vector from an element transformation, optionally ensure aligned
589 : // (| normal ⋅ align | > 0)
590 240 : void Normal(mfem::ElementTransformation &T, mfem::Vector &normal,
591 : const mfem::Vector *const align)
592 : {
593 : const mfem::IntegrationPoint &ip = mfem::Geometries.GetCenter(T.GetGeometryType());
594 : T.SetIntPoint(&ip);
595 240 : mfem::CalcOrtho(T.Jacobian(), normal);
596 240 : normal /= normal.Norml2();
597 240 : if (align && (normal * (*align) < 0))
598 : {
599 12 : normal *= -1;
600 : }
601 240 : }
602 :
603 : // Compute the centroid of a container of vertices.
604 : template <typename T>
605 26 : mfem::Vector ComputeCentroid(const std::unique_ptr<mfem::Mesh> &mesh, const T &vertidxs)
606 : {
607 : mfem::Vector centroid(mesh->SpaceDimension());
608 26 : centroid = 0.0;
609 : int c = 0;
610 156 : for (const int v : vertidxs)
611 : {
612 130 : centroid += mfem::Vector(mesh->GetVertex(v), 3);
613 130 : c++;
614 : }
615 52 : return centroid /= c;
616 : }
617 :
618 : // Compute the normal vector for a set of elements. If "inside" is true, normal will
619 : // point inside the mesh, otherwise it will point outside the mesh.
620 : template <typename T>
621 2 : mfem::Vector ComputeNormal(const std::unique_ptr<mfem::Mesh> &mesh, const T &elem_set,
622 : bool inside, bool check_planar = true)
623 : {
624 : const int sdim = mesh->SpaceDimension();
625 2 : mfem::IsoparametricTransformation trans;
626 : mfem::Vector normal(sdim), last_normal(sdim), align(sdim);
627 2 : normal = 0.0;
628 2 : last_normal = 0.0;
629 : mfem::Array<int> vert_bdr;
630 :
631 : // Ensure that the computed normal points "inside" or "outside".
632 26 : auto Alignment = [&](int el, auto &align)
633 : {
634 : int eladj, info;
635 24 : mesh->GetBdrElementAdjacentElement(el, eladj, info);
636 24 : mesh->GetElementCenter(eladj, align);
637 24 : mesh->GetBdrElementVertices(el, vert_bdr);
638 24 : align -= ComputeCentroid(mesh, vert_bdr);
639 24 : if (!inside) // align points inwards
640 : {
641 12 : align *= -1;
642 : }
643 : };
644 :
645 26 : for (auto elem : elem_set)
646 : {
647 24 : Alignment(elem, align);
648 24 : mesh->GetBdrElementTransformation(elem, &trans);
649 24 : mesh::Normal(trans, normal, &align);
650 24 : if (!check_planar)
651 : {
652 : break; // If not checking planar, use the first.
653 : }
654 24 : MFEM_VERIFY((last_normal * last_normal == 0.0) || ((last_normal * normal - 1) < 1e-8),
655 : "Periodic boundary mapping is only supported for planar boundaries!");
656 24 : last_normal = normal;
657 : }
658 2 : return normal;
659 2 : }
660 :
661 : struct Frame
662 : {
663 0 : Frame(const mfem::Vector &o) : origin(o)
664 : {
665 0 : for (auto &x : basis)
666 : {
667 0 : x.SetSize(3);
668 0 : x = 0.0;
669 : }
670 0 : }
671 : mfem::Vector origin;
672 : std::array<mfem::Vector, 3> basis;
673 : };
674 :
675 0 : Frame Find3DFrame(std::unique_ptr<mfem::Mesh> &mesh,
676 : const std::unordered_set<int> &vertidxs, const mfem::Vector ¢roid,
677 : const mfem::Vector &normal, double mesh_dim)
678 : {
679 0 : Frame frame(centroid);
680 0 : frame.basis[0] = normal;
681 :
682 : // For each point, compute its distance to the centroid.
683 : std::map<int, std::vector<int>, std::greater<int>> dist2points;
684 0 : for (const int v : vertidxs)
685 : {
686 : auto dist = centroid.DistanceTo(mesh->GetVertex(v));
687 : // Convert dist to integer to avoid floating point differences.
688 0 : dist2points[std::round(dist / mesh_dim * 1e8)].push_back(v);
689 : }
690 :
691 0 : for (const auto &[dist, verts] : dist2points)
692 : {
693 0 : if (verts.size() > 1 || dist == 0)
694 : {
695 : continue;
696 : }
697 0 : frame.basis[1] = mesh->GetVertex(verts.front());
698 0 : frame.basis[1] -= centroid;
699 0 : frame.basis[1] /= frame.basis[1].Norml2();
700 : break;
701 : }
702 :
703 : // Define final point by computing the cross product.
704 0 : frame.basis[0].cross3D(frame.basis[1], frame.basis[2]);
705 :
706 0 : return frame;
707 0 : }
708 :
709 : // Calculate the rotation matrix between two vectors.
710 0 : void ComputeRotation(const mfem::Vector &normal1, const mfem::Vector &normal2,
711 : mfem::DenseMatrix &transformation)
712 : {
713 0 : mfem::DenseMatrix R(3), vx(3), vx2(3);
714 :
715 : mfem::Vector v(normal1.Size());
716 0 : normal1.cross3D(normal2, v);
717 0 : double c = normal1 * normal2;
718 :
719 0 : vx(0, 1) = -v[2];
720 0 : vx(0, 2) = v[1];
721 0 : vx(1, 0) = v[2];
722 0 : vx(1, 2) = -v[0];
723 0 : vx(2, 0) = -v[1];
724 0 : vx(2, 1) = v[0];
725 :
726 0 : R(0, 0) = R(1, 1) = R(2, 2) = 1.0;
727 0 : R += vx;
728 0 : Mult(vx, vx, vx2);
729 0 : if (std::abs(1.0 + c) > 1e-8)
730 : {
731 0 : vx2.Set(1.0 / (1.0 + c), vx2);
732 : }
733 0 : R += vx2;
734 :
735 0 : for (int i = 0; i < 3; i++)
736 : {
737 0 : for (int j = 0; j < 3; j++)
738 : {
739 0 : transformation(i, j) = R(i, j);
740 : }
741 : }
742 0 : }
743 :
744 0 : mfem::DenseMatrix ComputeAffineTransformationMatrix(const Frame &donor,
745 : const Frame &receiver)
746 : {
747 :
748 0 : mfem::DenseMatrix A(4, 4);
749 0 : A = 0.0;
750 0 : if (donor.basis[1].Norml2() > 0.0 && receiver.basis[1].Norml2() > 0.0)
751 : {
752 : // Stably compute the rotation matrix from unit vectors.
753 : Eigen::Matrix3d source, target;
754 : Eigen::Vector3d source_centroid, target_centroid;
755 0 : for (int i = 0; i < 3; i++)
756 : {
757 0 : for (int j = 0; j < 3; j++)
758 : {
759 0 : source(i, j) = donor.basis[j](i);
760 0 : target(i, j) = receiver.basis[j](i);
761 : }
762 0 : source_centroid(i) = donor.origin(i);
763 0 : target_centroid(i) = receiver.origin(i);
764 : }
765 : Eigen::Matrix3d R = source * target.transpose();
766 : Eigen::JacobiSVD<Eigen::Matrix3d> svd(R, Eigen::ComputeFullU | Eigen::ComputeFullV);
767 0 : R = svd.matrixV() * svd.matrixU().transpose();
768 :
769 : // Account for possible reflection in R (det(R) = -1).
770 : Eigen::DiagonalMatrix<double, 3> Diag(1.0, 1.0, R.determinant());
771 0 : R = svd.matrixV() * Diag * svd.matrixU().transpose();
772 :
773 : // Compute translation and form transformation matrix.
774 0 : const Eigen::Vector3d translation = target_centroid - R * source_centroid;
775 0 : for (int i = 0; i < 3; i++)
776 : {
777 0 : for (int j = 0; j < 3; j++)
778 : {
779 0 : A(i, j) = R(i, j);
780 : }
781 0 : A(i, 3) = translation(i);
782 : }
783 0 : A(3, 3) = 1.0;
784 : }
785 : else
786 : {
787 : // If the donor or receiver basis is ambiguous, we assume no rotation around the
788 : // normals, and that the rotation comes only from realigning normal vectors.
789 0 : ComputeRotation(donor.basis[0], receiver.basis[0], A);
790 0 : for (int i = 0; i < 3; i++)
791 : {
792 0 : A(i, 3) = receiver.origin(i) - donor.origin(i);
793 : }
794 0 : A(3, 3) = 1.0;
795 : }
796 :
797 0 : return A;
798 0 : }
799 :
800 : // Create the vertex mapping between sets of donor and receiver pts related
801 : // by an affine transformation matrix.
802 0 : std::vector<int> CreatePeriodicVertexMapping(std::unique_ptr<mfem::Mesh> &mesh,
803 : const std::unordered_set<int> &donor_v,
804 : const std::unordered_set<int> &receiver_v,
805 : const mfem::DenseMatrix &transform,
806 : double tol = 1e-6)
807 : {
808 : // Similar to MFEM's CreatePeriodicVertexMapping, maps from replica to primary vertex.
809 : std::unordered_map<int, int> replica2primary;
810 :
811 : // KD-tree containing all the receiver points.
812 : mfem::KDTree3D kdtree;
813 0 : for (const int v : receiver_v)
814 : {
815 0 : kdtree.AddPoint(mesh->GetVertex(v), v);
816 : }
817 : kdtree.Sort();
818 :
819 : // Loop over donor points and find the corresponding receiver point.
820 : mfem::Vector from(4), to(4);
821 0 : for (int vi : donor_v)
822 : {
823 : // TODO: mfem patch to allow SetVector direct from pointer.
824 0 : std::copy(mesh->GetVertex(vi), mesh->GetVertex(vi) + 3, from.begin());
825 0 : from[3] = 1.0; // reset
826 0 : transform.Mult(from, to); // receiver = transform * donor
827 :
828 0 : const int vj = kdtree.FindClosestPoint(to.GetData());
829 0 : std::copy(mesh->GetVertex(vj), mesh->GetVertex(vj) + 3, from.begin());
830 0 : from -= to; // Check that the loaded vertex is identical to the transformed
831 0 : MFEM_VERIFY(from.Norml2() < tol,
832 : "Could not match points on periodic boundaries, transformed donor point "
833 : "does not correspond to a receiver point!");
834 0 : MFEM_VERIFY(replica2primary.find(vj) == replica2primary.end(),
835 : "Could not match points on periodic boundaries, multiple donor points map "
836 : "to the same receiver point!")
837 0 : replica2primary[vj] = vi;
838 : }
839 :
840 0 : std::vector<int> v2v(mesh->GetNV());
841 : std::iota(v2v.begin(), v2v.end(), 0);
842 0 : for (const auto &[r, p] : replica2primary)
843 : {
844 0 : v2v[r] = p;
845 : }
846 0 : return v2v;
847 : }
848 :
849 : // Determine the vertex mapping between donor and receiver boundary attributes.
850 : // Uses the translation vector or affine transformation matrix specified in the
851 : // configuration file. If not provided, attempts to automatically detect the
852 : // affine transformation between donor and receiver boundary vertices.
853 : std::vector<int>
854 1 : DeterminePeriodicVertexMapping(std::unique_ptr<mfem::Mesh> &mesh,
855 : const struct palace::config::PeriodicData &data,
856 : const double tol)
857 : {
858 : // Get mesh dimensions, will be used to define a reasonable tolerance in mesh units.
859 : mfem::Vector bbmin, bbmax;
860 1 : mesh->GetBoundingBox(bbmin, bbmax);
861 1 : bbmax -= bbmin;
862 1 : const double mesh_dim = bbmax.Norml2();
863 1 : const double mesh_tol = tol * mesh_dim;
864 :
865 : // Identify donor and receiver vertices and elements.
866 : const auto &da = data.donor_attributes, &ra = data.receiver_attributes;
867 : std::unordered_set<int> bdr_v_donor, bdr_v_receiver;
868 : std::unordered_set<int> bdr_e_donor, bdr_e_receiver;
869 49 : for (int be = 0; be < mesh->GetNBE(); be++)
870 : {
871 48 : int attr = mesh->GetBdrAttribute(be);
872 48 : auto donor = std::find(da.begin(), da.end(), attr) != da.end();
873 48 : auto receiver = std::find(ra.begin(), ra.end(), attr) != ra.end();
874 48 : if (donor || receiver)
875 : {
876 : int el, info;
877 24 : mesh->GetBdrElementAdjacentElement(be, el, info);
878 : mfem::Array<int> vertidxs;
879 : mesh->GetBdrElementVertices(be, vertidxs);
880 24 : (donor ? bdr_e_donor : bdr_e_receiver).insert(be);
881 24 : (donor ? bdr_v_donor : bdr_v_receiver).insert(vertidxs.begin(), vertidxs.end());
882 : }
883 : }
884 :
885 1 : MFEM_VERIFY(bdr_v_donor.size() == bdr_v_receiver.size(),
886 : "Different number of "
887 : "vertices on donor and receiver boundaries. Cannot create periodic mesh.");
888 :
889 : // Check if mesh has enough elements in periodic direction. MFEM's periodicity
890 : // fails for meshes with <=2 elements in the period direction.
891 : // Compare the number of mesh elements to the number of periodic boundary
892 : // elements.
893 1 : const int num_periodic_bc_elems = bdr_e_donor.size() + bdr_e_receiver.size();
894 : mfem::Array<mfem::Geometry::Type> geoms;
895 1 : mesh->GetGeometries(3, geoms);
896 1 : if (geoms.Size() == 1 && geoms[0] == mfem::Geometry::TETRAHEDRON)
897 : {
898 : // Pure tet mesh.
899 0 : MFEM_VERIFY(mesh->GetNE() > 3 * num_periodic_bc_elems,
900 : "Not enough mesh elements in periodic direction!");
901 : }
902 : else
903 : {
904 : // No tets.
905 1 : MFEM_VERIFY(mesh->GetNE() > num_periodic_bc_elems,
906 : "Not enough mesh elements in periodic direction!");
907 : }
908 :
909 : // Determine the affine transformation between donor and receiver points.
910 : // Use the translation vector or affine transformation matrix if provided
911 : // in the config file, otherwise automatically detect the transformation.
912 1 : mfem::DenseMatrix transformation(4);
913 1 : if (std::any_of(data.affine_transform.begin(), data.affine_transform.end(),
914 : [](auto y) { return std::abs(y) > 0.0; }))
915 : {
916 : // Use user-provided affine transformation matrix.
917 0 : for (int i = 0; i < 4; i++)
918 : {
919 0 : for (int j = 0; j < 4; j++)
920 : {
921 0 : transformation(i, j) = data.affine_transform[i * 4 + j]; // row major conversion
922 : }
923 : }
924 : }
925 : else
926 : {
927 : // Automatically detect transformation.
928 : // Compute the centroid for each boundary.
929 1 : auto donor_centroid = ComputeCentroid(mesh, bdr_v_donor);
930 1 : auto receiver_centroid = ComputeCentroid(mesh, bdr_v_receiver);
931 :
932 : // Compute the normal vector for each boundary.
933 1 : auto donor_normal = ComputeNormal(mesh, bdr_e_donor, true);
934 1 : auto receiver_normal = ComputeNormal(mesh, bdr_e_receiver, false);
935 :
936 : // Return empty mapping if centroids and normal vectors are the same (up to a sign).
937 1 : mfem::Vector diff = donor_centroid;
938 1 : diff -= receiver_centroid;
939 1 : double dot = donor_normal * receiver_normal;
940 1 : if (diff.Norml2() < mesh_tol && std::abs(std::abs(dot) - 1.0) < mesh_tol)
941 : {
942 1 : return {};
943 : }
944 :
945 : // Compute a frame (origin, normal, and two in plane points) for each boundary.
946 : auto donor_frame =
947 0 : Find3DFrame(mesh, bdr_v_donor, donor_centroid, donor_normal, mesh_dim);
948 : auto receiver_frame =
949 0 : Find3DFrame(mesh, bdr_v_receiver, receiver_centroid, receiver_normal, mesh_dim);
950 :
951 : // Compute the affine transformation matrix.
952 0 : transformation = ComputeAffineTransformationMatrix(donor_frame, receiver_frame);
953 0 : }
954 : return CreatePeriodicVertexMapping(mesh, bdr_v_donor, bdr_v_receiver, transformation,
955 0 : mesh_tol);
956 1 : }
957 :
958 : } // namespace palace::mesh
|
