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 1902 : TriangulatorInterface::TriangulatorInterface(UnstructuredMesh & mesh)
108 1742 : : _mesh(mesh),
109 1742 : _holes(nullptr),
110 1742 : _markers(nullptr),
111 1742 : _regions(nullptr),
112 1742 : _elem_type(TRI3),
113 1742 : _desired_area(0.1),
114 1742 : _minimum_angle(20.0),
115 1742 : _triangulation_type(GENERATE_CONVEX_HULL),
116 1742 : _insert_extra_points(false),
117 1742 : _smooth_after_generating(true),
118 1742 : _quiet(true),
119 2062 : _auto_area_function(nullptr)
120 1902 : {}
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 2062 : int TriangulatorInterface::get_interpolate_boundary_points () const
138 : {
139 : // backwards compatibility - someone might have turned this off via
140 : // the old API
141 2062 : if (!_insert_extra_points)
142 284 : return 0;
143 :
144 221 : return _interpolate_boundary_points;
145 : }
146 :
147 :
148 :
149 2275 : void TriangulatorInterface::elems_to_segments()
150 : {
151 : // Don't try to override manually specified segments
152 2275 : 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 2200 : if (_mesh.n_elem())
160 : {
161 : // Mapping from points to node ids, to back those out from
162 : // MeshedHole results later
163 40 : std::map<Point, dof_id_type> point_id_map;
164 :
165 5373 : for (Node * node : _mesh.node_ptr_range())
166 : {
167 : // We're not going to support overlapping nodes on the boundary
168 4241 : 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 4241 : point_id_map.emplace(*node, node->id());
174 536 : }
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 6842 : for (Elem * elem : _mesh.element_ptr_range())
181 3886 : for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
182 588 : point_id_map.erase(elem->point(n));
183 :
184 : // We'll steal the ordering calculation from
185 : // the MeshedHole code
186 1166 : 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 28 : std::unordered_set<Node *> nodes_to_delete;
198 :
199 4904 : for (Elem * elem : _mesh.element_ptr_range())
200 3105 : for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
201 1284 : nodes_to_delete.insert(elem->node_ptr(n));
202 :
203 363 : 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 363 : _mesh.clear_elems();
213 :
214 : // Make segments from boundary nodes; also make sure we don't
215 : // delete them.
216 363 : const std::size_t np = mh.n_points();
217 1815 : for (auto i : make_range(np))
218 : {
219 1452 : const Point pt = mh.point(i);
220 1452 : const dof_id_type id0 = libmesh_map_find(point_id_map, pt);
221 1452 : nodes_to_delete.erase(_mesh.node_ptr(id0));
222 1452 : const Point next_pt = mh.point((np+i+1)%np);
223 1452 : const dof_id_type id1 = libmesh_map_find(point_id_map, next_pt);
224 1452 : this->segments.emplace_back(id0, id1);
225 2020 : for (auto m : make_range(mh.n_midpoints()))
226 568 : this->segment_midpoints.emplace_back(mh.midpoint(m, i));
227 : }
228 :
229 2016 : for (Node * node : nodes_to_delete)
230 1653 : _mesh.delete_node(node);
231 :
232 363 : if (this->_verify_hole_boundaries && _holes)
233 0 : this->verify_holes(mh);
234 335 : }
235 : }
236 :
237 :
238 :
239 2062 : void TriangulatorInterface::nodes_to_segments(dof_id_type max_node_id)
240 : {
241 : // Don't try to override manually specified segments, or try to add
242 : // segments if we're doing a convex hull
243 2062 : if (!this->segments.empty() || _triangulation_type != PSLG)
244 248 : return;
245 :
246 1210 : for (auto node_it = _mesh.nodes_begin(),
247 1210 : node_end = _mesh.nodes_end();
248 5828 : node_it != node_end;)
249 : {
250 4664 : Node * node = *node_it;
251 :
252 : // If we're out of boundary nodes, the rest are going to be
253 : // Steiner points or hole points
254 4664 : if (node->id() >= max_node_id)
255 0 : break;
256 :
257 4476 : ++node_it;
258 :
259 9140 : Node * next_node = (node_it == node_end) ?
260 1306 : *_mesh.nodes_begin() : *node_it;
261 :
262 4664 : this->segments.emplace_back(node->id(), next_node->id());
263 : }
264 :
265 1164 : if (this->_verify_hole_boundaries && _holes)
266 : {
267 48 : std::vector<Point> outer_pts;
268 3255 : for (auto segment : this->segments)
269 2604 : outer_pts.push_back(_mesh.point(segment.first));
270 :
271 699 : ArbitraryHole ah(outer_pts);
272 651 : this->verify_holes(ah);
273 603 : }
274 : }
275 :
276 :
277 :
278 2062 : void TriangulatorInterface::insert_any_extra_boundary_points()
279 : {
280 : // If the initial PSLG is really simple, e.g. an L-shaped domain or
281 : // a square/rectangle, the resulting triangulation may be very
282 : // "structured" looking. Sometimes this is a problem if your
283 : // intention is to work with an "unstructured" looking grid. We can
284 : // attempt to work around this limitation by inserting midpoints
285 : // into the original PSLG. Inserting additional points into a
286 : // set of points meant to be a convex hull usually makes less sense.
287 :
288 2062 : const int n_interpolated = this->get_interpolate_boundary_points();
289 2062 : if ((_triangulation_type==PSLG) && n_interpolated)
290 : {
291 : // If we were lucky enough to start with contiguous node ids,
292 : // let's keep them that way.
293 221 : dof_id_type nn = _mesh.max_node_id();
294 :
295 : std::vector<std::pair<unsigned int, unsigned int>> old_segments =
296 231 : std::move(this->segments);
297 :
298 : // We expect to have converted any elems and/or nodes into
299 : // segments by now.
300 10 : libmesh_assert(!old_segments.empty());
301 :
302 10 : this->segments.clear();
303 :
304 : // Insert a new point on each segment at evenly spaced locations
305 : // between existing boundary points.
306 : // np=index into new points vector
307 : // n =index into original points vector
308 1105 : for (auto old_segment : old_segments)
309 : {
310 884 : Node * begin_node = _mesh.node_ptr(old_segment.first);
311 884 : Node * end_node = _mesh.node_ptr(old_segment.second);
312 884 : dof_id_type current_id = begin_node->id();
313 2920 : for (auto i : make_range(n_interpolated))
314 : {
315 : // new points are equispaced along the original segments
316 : const Point new_point =
317 2036 : ((n_interpolated-i) * *(Point *)(begin_node) +
318 2036 : (i+1) * *(Point *)(end_node)) /
319 2116 : (n_interpolated + 1);
320 2036 : Node * next_node = _mesh.add_point(new_point, nn++);
321 1876 : this->segments.emplace_back(current_id,
322 2036 : next_node->id());
323 2036 : current_id = next_node->id();
324 : }
325 804 : this->segments.emplace_back(current_id,
326 884 : end_node->id());
327 : }
328 : }
329 2062 : }
330 :
331 :
332 1566 : void TriangulatorInterface::increase_triangle_order()
333 : {
334 1566 : switch (_elem_type)
335 : {
336 42 : case TRI3:
337 : // Nothing to do if we're not requested to increase order
338 1491 : return;
339 75 : case TRI6:
340 75 : _mesh.all_second_order();
341 4 : break;
342 0 : case TRI7:
343 0 : _mesh.all_complete_order();
344 0 : break;
345 0 : default:
346 0 : libmesh_not_implemented();
347 : }
348 :
349 : // If we have any midpoint location data, we'll want to look it up
350 : // by point. all_midpoints[{p, m}] will be the mth midpoint
351 : // location following after point p (when traversing a triangle
352 : // counter-clockwise)
353 4 : std::map<std::pair<Point, unsigned int>, Point> all_midpoints;
354 : unsigned int n_midpoints =
355 83 : this->segment_midpoints.size() / this->segments.size();
356 4 : libmesh_assert_equal_to(this->segments.size() * n_midpoints,
357 : this->segment_midpoints.size());
358 75 : for (auto m : make_range(n_midpoints))
359 0 : for (auto i : make_range(this->segments.size()))
360 : {
361 0 : const Point & p = _mesh.point(this->segments[i].first);
362 0 : all_midpoints[{p,m}] =
363 0 : this->segment_midpoints[i*n_midpoints+m];
364 : }
365 :
366 75 : if (_holes)
367 150 : for (const Hole * hole : *_holes)
368 : {
369 75 : if (!hole->n_midpoints())
370 0 : continue;
371 75 : if (!n_midpoints)
372 75 : n_midpoints = hole->n_midpoints();
373 0 : else if (hole->n_midpoints() != n_midpoints)
374 0 : libmesh_not_implemented_msg
375 : ("Differing boundary midpoint counts " <<
376 : hole->n_midpoints() << " and " << n_midpoints);
377 :
378 : // Our inner holes are expected to have points in
379 : // counter-clockwise order, which is backwards from how we
380 : // want to traverse them when iterating in counter-clockwise
381 : // order over a triangle, so we'll need to reverse our maps
382 : // carefully here.
383 75 : const auto n_hole_points = hole->n_points();
384 4 : libmesh_assert(n_hole_points);
385 150 : for (auto m : make_range(n_midpoints))
386 : {
387 600 : for (auto i : make_range(n_hole_points-1))
388 : {
389 525 : const Point & p = hole->point(i+1);
390 525 : all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, i);
391 : }
392 75 : const Point & p = hole->point(0);
393 75 : all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, n_hole_points-1);
394 : }
395 : }
396 :
397 : // The n_midpoints > 1 case is for future proofing, but in the
398 : // present we have EDGE4 and no TRI10 yet.
399 75 : if (n_midpoints > 1)
400 0 : libmesh_not_implemented_msg
401 : ("Cannot construct triangles with more than 1 midpoint per edge");
402 :
403 75 : if (!n_midpoints)
404 0 : return;
405 :
406 1898 : for (Elem * elem : _mesh.element_ptr_range())
407 : {
408 : // This should only be called right after we've finished
409 : // converting a triangulation to higher order
410 48 : libmesh_assert_equal_to(elem->n_vertices(), 3);
411 48 : libmesh_assert_not_equal_to(elem->default_order(), FIRST);
412 :
413 3600 : for (auto n : make_range(3))
414 : {
415 : // Only hole/outer boundary segments need adjusted midpoints
416 2844 : if (elem->neighbor_ptr(n))
417 1704 : continue;
418 :
419 96 : const Point & p = elem->point(n);
420 :
421 900 : if (const auto it = all_midpoints.find({p,0});
422 48 : it != all_midpoints.end())
423 632 : elem->point(n+3) = it->second;
424 : }
425 67 : }
426 : }
427 :
428 :
429 651 : void TriangulatorInterface::verify_holes(const Hole & outer_bdy)
430 : {
431 1870 : for (const Hole * hole : *_holes)
432 : {
433 7550 : for (const Hole * hole2 : *_holes)
434 : {
435 6331 : if (hole == hole2)
436 1179 : continue;
437 :
438 25560 : for (auto i : make_range(hole2->n_points()))
439 20448 : if (hole->contains(hole2->point(i)))
440 0 : libmesh_error_msg
441 : ("Found point " << hole2->point(i) <<
442 : " on one hole boundary and another's interior");
443 : }
444 :
445 6695 : for (auto i : make_range(hole->n_points()))
446 5476 : if (!outer_bdy.contains(hole->point(i)))
447 0 : libmesh_error_msg
448 : ("Found point " << hole->point(i) <<
449 : " on hole boundary but outside outer boundary");
450 : }
451 651 : }
452 :
453 :
454 500 : unsigned int TriangulatorInterface::total_hole_points()
455 : {
456 : // If the holes vector is non-nullptr (and non-empty) we need to determine
457 : // the number of additional points which the holes will add to the
458 : // triangulation.
459 : // Note that the number of points is always equal to the number of segments
460 : // that form the holes.
461 250 : unsigned int n_hole_points = 0;
462 :
463 500 : if (_holes)
464 40 : for (const auto & hole : *_holes)
465 : {
466 24 : n_hole_points += hole->n_points();
467 : // A hole at least has one enclosure.
468 : // Points on enclosures are ordered so that we can add segments implicitly.
469 : // Elements in segment_indices() indicates the starting points of all enclosures.
470 : // The last element in segment_indices() is the number of total points.
471 12 : libmesh_assert_greater(hole->segment_indices().size(), 1);
472 12 : libmesh_assert_equal_to(hole->segment_indices().back(), hole->n_points());
473 : }
474 :
475 500 : return n_hole_points;
476 : }
477 :
478 0 : void TriangulatorInterface::set_auto_area_function(const Parallel::Communicator &comm,
479 : const unsigned int num_nearest_pts,
480 : const unsigned int power,
481 : const Real background_value,
482 : const Real background_eff_dist)
483 : {
484 0 : _auto_area_function = std::make_unique<AutoAreaFunction>(comm, num_nearest_pts, power, background_value, background_eff_dist);
485 0 : }
486 :
487 0 : FunctionBase<Real> * TriangulatorInterface::get_auto_area_function()
488 : {
489 0 : if (!_auto_area_function->initialized())
490 : {
491 : // Points and target element sizes for the interpolation
492 0 : std::vector<Point> function_points;
493 0 : std::vector<Real> function_sizes;
494 0 : calculate_auto_desired_area_samples(function_points, function_sizes);
495 0 : _auto_area_function->init_mfi(function_points, function_sizes);
496 : }
497 0 : return _auto_area_function.get();
498 : }
499 :
500 0 : void TriangulatorInterface::calculate_auto_desired_area_samples(std::vector<Point> & function_points,
501 : std::vector<Real> & function_sizes,
502 : const Real & area_factor)
503 : {
504 : // Get the hole mesh of the outer boundary
505 : // Holes should already be attached if applicable when this function is called
506 0 : const TriangulatorInterface::MeshedHole bdry_mh { _mesh, this->_bdy_ids };
507 : // Collect all the centroid points of the outer boundary segments
508 : // and the corresponding element sizes
509 0 : for (unsigned int i = 0; i < bdry_mh.n_points(); i++)
510 : {
511 0 : function_points.push_back((bdry_mh.point(i) + bdry_mh.point((i + 1) % bdry_mh.n_points())) /
512 0 : Real(2.0));
513 0 : function_sizes.push_back(
514 0 : (bdry_mh.point(i) - bdry_mh.point((i + 1) % bdry_mh.n_points())).norm());
515 : }
516 : // If holes are present, do the same for the hole boundaries
517 0 : if(_holes)
518 0 : for (const Hole * hole : *_holes)
519 : {
520 0 : for (unsigned int i = 0; i < hole->n_points(); i++)
521 : {
522 0 : function_points.push_back(
523 0 : (hole->point(i) + hole->point((i + 1) % hole->n_points())) / Real(2.0));
524 0 : function_sizes.push_back(
525 0 : (hole->point(i) - hole->point((i + 1) % hole->n_points())).norm());
526 : }
527 : }
528 :
529 0 : std::for_each(
530 0 : function_sizes.begin(), function_sizes.end(), [&area_factor](Real & a) { a = a * a * area_factor * std::sqrt(3.0) / 4.0; });
531 :
532 0 : }
533 : } // namespace libMesh
534 :
|
