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