Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/libmesh_config.h"
20 :
21 : // Avoid anything here that might break precise IEEE754 compatibility,
22 : // to give us more reproduceable results.
23 : //
24 : // We can't disable excess x87 precision from pragmas, but hopefully
25 : // anyone optimizing will be using SSE instead anyway.
26 : #if defined(__clang__)
27 : # pragma float_control(precise, on)
28 : # pragma clang fp contract(off) reassociate(off)
29 : #elif defined(__NVCOMPILER)
30 : // We can't get -Kieee from pragmas, but so far FMA contractions are
31 : // the only thing we've caught breaking us, and nvc++ inherits a
32 : // pragma for those from LLVM.
33 : # pragma clang fp contract(off) reassociate(off)
34 : #elif defined(__GNUC__)
35 : // GCC goes last, because other compilers define __GNUC__, because
36 : // they're liars.
37 : # pragma GCC optimize("-fno-unsafe-math-optimizations")
38 : # pragma GCC optimize("-ffp-contract=off")
39 : #endif
40 :
41 : // libmesh includes
42 : #include "libmesh/mesh_triangle_interface.h"
43 : #include "libmesh/unstructured_mesh.h"
44 : #include "libmesh/face_tri3.h"
45 : #include "libmesh/face_tri6.h"
46 : #include "libmesh/mesh_generation.h"
47 : #include "libmesh/mesh_smoother_laplace.h"
48 : #include "libmesh/boundary_info.h"
49 : #include "libmesh/mesh_triangle_holes.h"
50 : #include "libmesh/mesh_triangle_wrapper.h"
51 : #include "libmesh/enum_elem_type.h"
52 : #include "libmesh/enum_order.h"
53 : #include "libmesh/enum_to_string.h"
54 : #include "libmesh/utility.h"
55 :
56 : #include "libmesh/meshfree_interpolation.h"
57 :
58 : // C/C++ includes
59 : #include <sstream>
60 :
61 :
62 : namespace libMesh
63 : {
64 : //
65 : // Function definitions for the AutoAreaFunction class
66 : //
67 :
68 : // Constructor
69 0 : AutoAreaFunction::AutoAreaFunction (const Parallel::Communicator &comm,
70 : const unsigned int num_nearest_pts,
71 : const unsigned int power,
72 : const Real background_value,
73 0 : const Real background_eff_dist):
74 0 : _comm(comm),
75 0 : _num_nearest_pts(num_nearest_pts),
76 0 : _power(power),
77 0 : _background_value(background_value),
78 0 : _background_eff_dist(background_eff_dist),
79 0 : _auto_area_mfi(std::make_unique<InverseDistanceInterpolation<3>>(_comm, _num_nearest_pts, _power, _background_value, _background_eff_dist))
80 : {
81 0 : this->_initialized = false;
82 0 : this->_is_time_dependent = false;
83 0 : }
84 :
85 : // Destructor
86 0 : AutoAreaFunction::~AutoAreaFunction () = default;
87 :
88 0 : void AutoAreaFunction::init_mfi (const std::vector<Point> & input_pts,
89 : const std::vector<Real> & input_vals)
90 : {
91 0 : std::vector<std::string> field_vars{"f"};
92 0 : _auto_area_mfi->set_field_variables(field_vars);
93 0 : _auto_area_mfi->get_source_points() = input_pts;
94 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
95 : std::vector<Number> input_complex_vals;
96 0 : for (const auto & input_val : input_vals)
97 0 : input_complex_vals.push_back(Complex (input_val, 0.0));
98 0 : _auto_area_mfi->get_source_vals() = input_complex_vals;
99 : #else
100 0 : _auto_area_mfi->get_source_vals() = input_vals;
101 : #endif
102 0 : _auto_area_mfi->prepare_for_use();
103 0 : this->_initialized = true;
104 0 : }
105 :
106 0 : Real AutoAreaFunction::operator() (const Point & p,
107 : const Real /*time*/)
108 : {
109 0 : libmesh_assert(this->_initialized);
110 :
111 0 : std::vector<Point> target_pts;
112 0 : std::vector<Number> target_vals;
113 :
114 0 : target_pts.push_back(p);
115 0 : target_vals.resize(1);
116 :
117 0 : _auto_area_mfi->interpolate_field_data(_auto_area_mfi->field_variables(), target_pts, target_vals);
118 :
119 0 : return libmesh_real(target_vals.front());
120 : }
121 :
122 : //
123 : // Function definitions for the TriangulatorInterface class
124 : //
125 :
126 : // Constructor
127 1977 : TriangulatorInterface::TriangulatorInterface(UnstructuredMesh & mesh)
128 1809 : : _mesh(mesh),
129 1809 : _holes(nullptr),
130 1809 : _markers(nullptr),
131 1809 : _regions(nullptr),
132 1809 : _elem_type(TRI3),
133 1809 : _desired_area(0.1),
134 1809 : _minimum_angle(20.0),
135 1809 : _triangulation_type(GENERATE_CONVEX_HULL),
136 1809 : _insert_extra_points(false),
137 1809 : _smooth_after_generating(true),
138 1809 : _quiet(true),
139 2145 : _auto_area_function(nullptr)
140 1977 : {}
141 :
142 :
143 221 : void TriangulatorInterface::set_interpolate_boundary_points (int n_points)
144 : {
145 : // Maybe we'll reserve a meaning for negatives later?
146 10 : libmesh_assert(n_points >= 0);
147 :
148 221 : _interpolate_boundary_points = n_points;
149 :
150 : // backwards compatibility - someone (including us) might want to
151 : // query this via the old API.
152 221 : _insert_extra_points = n_points;
153 221 : }
154 :
155 :
156 :
157 2137 : int TriangulatorInterface::get_interpolate_boundary_points () const
158 : {
159 : // backwards compatibility - someone might have turned this off via
160 : // the old API
161 2137 : if (!_insert_extra_points)
162 288 : return 0;
163 :
164 221 : return _interpolate_boundary_points;
165 : }
166 :
167 :
168 :
169 2350 : void TriangulatorInterface::elems_to_segments()
170 : {
171 : // Don't try to override manually specified segments
172 2350 : if (!this->segments.empty())
173 4 : return;
174 :
175 : // If we have edges, they should form the polyline with the ordering
176 : // we want. Let's turn them into segments for later use, because
177 : // we're going to delete the original elements to replace with our
178 : // triangulation.
179 2275 : if (_mesh.n_elem())
180 : {
181 : // Mapping from points to node ids, to back those out from
182 : // MeshedHole results later
183 48 : std::map<Point, dof_id_type> point_id_map;
184 :
185 6119 : for (Node * node : _mesh.node_ptr_range())
186 : {
187 : // We're not going to support overlapping nodes on the boundary
188 4841 : libmesh_error_msg_if
189 : (point_id_map.count(*node),
190 : "TriangulatorInterface does not support overlapping nodes found at "
191 : << static_cast<Point&>(*node));
192 :
193 4841 : point_id_map.emplace(*node, node->id());
194 603 : }
195 :
196 : // We don't support directly generating Tri6, so for
197 : // compatibility with future stitching we need to be working
198 : // with first-order elements. Let's get rid of any non-vertex
199 : // nodes we just added.
200 7556 : for (Elem * elem : _mesh.element_ptr_range())
201 4486 : for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
202 687 : point_id_map.erase(elem->point(n));
203 :
204 : // We'll steal the ordering calculation from
205 : // the MeshedHole code
206 1320 : const TriangulatorInterface::MeshedHole mh { _mesh, this->_bdy_ids };
207 :
208 : // If we've specified only a subset of the mesh as our outer
209 : // boundary, then we may have nodes that don't actually fall
210 : // inside that boundary. Triangulator code doesn't like Steiner
211 : // points that aren't inside the triangulation domain, so we
212 : // need to get rid of them.
213 : //
214 : // Also, if we're using Edge3 elements to define our outer
215 : // boundary, we're only dealing with their 2 end nodes and we'll
216 : // need to get rid of their central nodes.
217 36 : std::unordered_set<Node *> nodes_to_delete;
218 :
219 5618 : for (Elem * elem : _mesh.element_ptr_range())
220 3705 : for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
221 1667 : nodes_to_delete.insert(elem->node_ptr(n));
222 :
223 438 : if (!this->_bdy_ids.empty())
224 : {
225 4048 : for (auto & node : _mesh.node_ptr_range())
226 1953 : if (!mh.contains(*node))
227 204 : nodes_to_delete.insert(node);
228 : }
229 :
230 : // And now we're done with elements. Delete them lest they have
231 : // dangling pointers to nodes we'll be deleting.
232 438 : _mesh.clear_elems();
233 :
234 : // Make segments from boundary nodes; also make sure we don't
235 : // delete them.
236 438 : const std::size_t np = mh.n_points();
237 2190 : for (auto i : make_range(np))
238 : {
239 1752 : const Point pt = mh.point(i);
240 1752 : const dof_id_type id0 = libmesh_map_find(point_id_map, pt);
241 1752 : nodes_to_delete.erase(_mesh.node_ptr(id0));
242 1752 : const Point next_pt = mh.point((np+i+1)%np);
243 1752 : const dof_id_type id1 = libmesh_map_find(point_id_map, next_pt);
244 1752 : this->segments.emplace_back(id0, id1);
245 2620 : for (auto m : make_range(mh.n_midpoints()))
246 : {
247 868 : this->segment_midpoints.emplace_back(mh.midpoint(m, i));
248 868 : this->segment_midpoints_keys.emplace_back(pt);
249 : }
250 : }
251 :
252 2391 : for (Node * node : nodes_to_delete)
253 1953 : _mesh.delete_node(node);
254 :
255 438 : if (this->_verify_hole_boundaries && _holes)
256 0 : this->verify_holes(mh);
257 402 : }
258 : }
259 :
260 :
261 :
262 2137 : void TriangulatorInterface::nodes_to_segments(dof_id_type max_node_id)
263 : {
264 : // Don't try to override manually specified segments, or try to add
265 : // segments if we're doing a convex hull
266 2137 : if (!this->segments.empty() || _triangulation_type != PSLG)
267 252 : return;
268 :
269 1210 : for (auto node_it = _mesh.nodes_begin(),
270 1210 : node_end = _mesh.nodes_end();
271 5828 : node_it != node_end;)
272 : {
273 4664 : Node * node = *node_it;
274 :
275 : // If we're out of boundary nodes, the rest are going to be
276 : // Steiner points or hole points
277 4664 : if (node->id() >= max_node_id)
278 0 : break;
279 :
280 4476 : ++node_it;
281 :
282 9140 : Node * next_node = (node_it == node_end) ?
283 1306 : *_mesh.nodes_begin() : *node_it;
284 :
285 4664 : this->segments.emplace_back(node->id(), next_node->id());
286 : }
287 :
288 1164 : if (this->_verify_hole_boundaries && _holes)
289 : {
290 48 : std::vector<Point> outer_pts;
291 3262 : for (auto segment : this->segments)
292 2610 : outer_pts.push_back(_mesh.point(segment.first));
293 :
294 701 : ArbitraryHole ah(outer_pts);
295 652 : this->verify_holes(ah);
296 603 : }
297 : }
298 :
299 :
300 :
301 2137 : void TriangulatorInterface::insert_any_extra_boundary_points()
302 : {
303 : // If the initial PSLG is really simple, e.g. an L-shaped domain or
304 : // a square/rectangle, the resulting triangulation may be very
305 : // "structured" looking. Sometimes this is a problem if your
306 : // intention is to work with an "unstructured" looking grid. We can
307 : // attempt to work around this limitation by inserting midpoints
308 : // into the original PSLG. Inserting additional points into a
309 : // set of points meant to be a convex hull usually makes less sense.
310 :
311 2137 : const int n_interpolated = this->get_interpolate_boundary_points();
312 2137 : if ((_triangulation_type==PSLG) && n_interpolated)
313 : {
314 : // If we were lucky enough to start with contiguous node ids,
315 : // let's keep them that way.
316 221 : dof_id_type nn = _mesh.max_node_id();
317 :
318 : std::vector<std::pair<unsigned int, unsigned int>> old_segments =
319 231 : std::move(this->segments);
320 :
321 : // We expect to have converted any elems and/or nodes into
322 : // segments by now.
323 10 : libmesh_assert(!old_segments.empty());
324 :
325 10 : this->segments.clear();
326 :
327 : // Insert a new point on each segment at evenly spaced locations
328 : // between existing boundary points.
329 : // np=index into new points vector
330 : // n =index into original points vector
331 1105 : for (auto old_segment : old_segments)
332 : {
333 884 : Node * begin_node = _mesh.node_ptr(old_segment.first);
334 884 : Node * end_node = _mesh.node_ptr(old_segment.second);
335 884 : dof_id_type current_id = begin_node->id();
336 2920 : for (auto i : make_range(n_interpolated))
337 : {
338 : // new points are equispaced along the original segments
339 : const Point new_point =
340 2036 : ((n_interpolated-i) * *(Point *)(begin_node) +
341 2036 : (i+1) * *(Point *)(end_node)) /
342 2116 : (n_interpolated + 1);
343 2036 : Node * next_node = _mesh.add_point(new_point, nn++);
344 1876 : this->segments.emplace_back(current_id,
345 2036 : next_node->id());
346 2036 : current_id = next_node->id();
347 : }
348 804 : this->segments.emplace_back(current_id,
349 884 : end_node->id());
350 : }
351 : }
352 2137 : }
353 :
354 :
355 1641 : void TriangulatorInterface::increase_triangle_order()
356 : {
357 1641 : switch (_elem_type)
358 : {
359 42 : case TRI3:
360 : // Nothing to do if we're not requested to increase order
361 1491 : return;
362 150 : case TRI6:
363 150 : _mesh.all_second_order();
364 8 : break;
365 0 : case TRI7:
366 0 : _mesh.all_complete_order();
367 0 : break;
368 0 : default:
369 0 : libmesh_not_implemented();
370 : }
371 :
372 : // If we have any midpoint location data, we'll want to look it up
373 : // by point. all_midpoints[{p, m}] will be the mth midpoint
374 : // location following after point p (when traversing a triangle
375 : // counter-clockwise)
376 8 : std::map<std::pair<Point, unsigned int>, Point> all_midpoints;
377 : unsigned int n_midpoints =
378 166 : this->segment_midpoints.size() / this->segments.size();
379 8 : libmesh_assert_equal_to(this->segments.size() * n_midpoints,
380 : this->segment_midpoints.size());
381 225 : for (auto m : make_range(n_midpoints))
382 375 : for (auto i : make_range(this->segments.size()))
383 : {
384 300 : const Point & p = segment_midpoints_keys[i*n_midpoints+m];
385 300 : all_midpoints[{p,m}] =
386 32 : this->segment_midpoints[i*n_midpoints+m];
387 : }
388 :
389 150 : if (_holes)
390 150 : for (const Hole * hole : *_holes)
391 : {
392 75 : if (!hole->n_midpoints())
393 0 : continue;
394 75 : if (!n_midpoints)
395 75 : n_midpoints = hole->n_midpoints();
396 0 : else if (hole->n_midpoints() != n_midpoints)
397 0 : libmesh_not_implemented_msg
398 : ("Differing boundary midpoint counts " <<
399 : hole->n_midpoints() << " and " << n_midpoints);
400 :
401 : // Our inner holes are expected to have points in
402 : // counter-clockwise order, which is backwards from how we
403 : // want to traverse them when iterating in counter-clockwise
404 : // order over a triangle, so we'll need to reverse our maps
405 : // carefully here.
406 75 : const auto n_hole_points = hole->n_points();
407 4 : libmesh_assert(n_hole_points);
408 150 : for (auto m : make_range(n_midpoints))
409 : {
410 600 : for (auto i : make_range(n_hole_points-1))
411 : {
412 525 : const Point & p = hole->point(i+1);
413 525 : all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, i);
414 : }
415 75 : const Point & p = hole->point(0);
416 75 : all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, n_hole_points-1);
417 : }
418 : }
419 :
420 : // The n_midpoints > 1 case is for future proofing, but in the
421 : // present we have EDGE4 and no TRI10 yet.
422 150 : if (n_midpoints > 1)
423 0 : libmesh_not_implemented_msg
424 : ("Cannot construct triangles with more than 1 midpoint per edge");
425 :
426 150 : if (!n_midpoints)
427 0 : return;
428 :
429 2336 : for (Elem * elem : _mesh.element_ptr_range())
430 : {
431 : // This should only be called right after we've finished
432 : // converting a triangulation to higher order
433 56 : libmesh_assert_equal_to(elem->n_vertices(), 3);
434 56 : libmesh_assert_not_equal_to(elem->default_order(), FIRST);
435 :
436 4200 : for (auto n : make_range(3))
437 : {
438 : // Only hole/outer boundary segments need adjusted midpoints
439 3318 : if (elem->neighbor_ptr(n))
440 1846 : continue;
441 :
442 128 : const Point & p = elem->point(n);
443 :
444 1200 : if (const auto it = all_midpoints.find({p,0});
445 64 : it != all_midpoints.end())
446 948 : elem->point(n+3) = it->second;
447 : }
448 134 : }
449 : }
450 :
451 :
452 652 : void TriangulatorInterface::verify_holes(const Hole & outer_bdy)
453 : {
454 1874 : for (const Hole * hole : *_holes)
455 : {
456 7562 : for (const Hole * hole2 : *_holes)
457 : {
458 6340 : if (hole == hole2)
459 1179 : continue;
460 :
461 25872 : for (auto i : make_range(hole2->n_points()))
462 20754 : if (hole->contains(hole2->point(i)))
463 0 : libmesh_error_msg
464 : ("Found point " << hole2->point(i) <<
465 : " on one hole boundary and another's interior");
466 : }
467 :
468 6851 : for (auto i : make_range(hole->n_points()))
469 5629 : if (!outer_bdy.contains(hole->point(i)))
470 0 : libmesh_error_msg
471 : ("Found point " << hole->point(i) <<
472 : " on hole boundary but outside outer boundary");
473 : }
474 652 : }
475 :
476 :
477 504 : unsigned int TriangulatorInterface::total_hole_points()
478 : {
479 : // If the holes vector is non-nullptr (and non-empty) we need to determine
480 : // the number of additional points which the holes will add to the
481 : // triangulation.
482 : // Note that the number of points is always equal to the number of segments
483 : // that form the holes.
484 252 : unsigned int n_hole_points = 0;
485 :
486 504 : if (_holes)
487 40 : for (const auto & hole : *_holes)
488 : {
489 24 : n_hole_points += hole->n_points();
490 : // A hole at least has one enclosure.
491 : // Points on enclosures are ordered so that we can add segments implicitly.
492 : // Elements in segment_indices() indicates the starting points of all enclosures.
493 : // The last element in segment_indices() is the number of total points.
494 12 : libmesh_assert_greater(hole->segment_indices().size(), 1);
495 12 : libmesh_assert_equal_to(hole->segment_indices().back(), hole->n_points());
496 : }
497 :
498 504 : return n_hole_points;
499 : }
500 :
501 0 : void TriangulatorInterface::set_auto_area_function(const Parallel::Communicator &comm,
502 : const unsigned int num_nearest_pts,
503 : const unsigned int power,
504 : const Real background_value,
505 : const Real background_eff_dist)
506 : {
507 0 : _auto_area_function = std::make_unique<AutoAreaFunction>(comm, num_nearest_pts, power, background_value, background_eff_dist);
508 0 : }
509 :
510 0 : FunctionBase<Real> * TriangulatorInterface::get_auto_area_function()
511 : {
512 0 : if (!_auto_area_function->initialized())
513 : {
514 : // Points and target element sizes for the interpolation
515 0 : std::vector<Point> function_points;
516 0 : std::vector<Real> function_sizes;
517 0 : calculate_auto_desired_area_samples(function_points, function_sizes);
518 0 : _auto_area_function->init_mfi(function_points, function_sizes);
519 : }
520 0 : return _auto_area_function.get();
521 : }
522 :
523 0 : void TriangulatorInterface::calculate_auto_desired_area_samples(std::vector<Point> & function_points,
524 : std::vector<Real> & function_sizes,
525 : const Real & area_factor)
526 : {
527 : // Get the hole mesh of the outer boundary
528 : // Holes should already be attached if applicable when this function is called
529 0 : const TriangulatorInterface::MeshedHole bdry_mh { _mesh, this->_bdy_ids };
530 : // Collect all the centroid points of the outer boundary segments
531 : // and the corresponding element sizes
532 0 : for (unsigned int i = 0; i < bdry_mh.n_points(); i++)
533 : {
534 0 : function_points.push_back((bdry_mh.point(i) + bdry_mh.point((i + 1) % bdry_mh.n_points())) /
535 0 : Real(2.0));
536 0 : function_sizes.push_back(
537 0 : (bdry_mh.point(i) - bdry_mh.point((i + 1) % bdry_mh.n_points())).norm());
538 : }
539 : // If holes are present, do the same for the hole boundaries
540 0 : if(_holes)
541 0 : for (const Hole * hole : *_holes)
542 : {
543 0 : for (unsigned int i = 0; i < hole->n_points(); i++)
544 : {
545 0 : function_points.push_back(
546 0 : (hole->point(i) + hole->point((i + 1) % hole->n_points())) / Real(2.0));
547 0 : function_sizes.push_back(
548 0 : (hole->point(i) - hole->point((i + 1) % hole->n_points())).norm());
549 : }
550 : }
551 :
552 0 : std::for_each(
553 0 : function_sizes.begin(), function_sizes.end(), [&area_factor](Real & a) { a = a * a * area_factor * std::sqrt(3.0) / 4.0; });
554 :
555 0 : }
556 : } // namespace libMesh
557 :
|
