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 : // 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 <limits>
40 : #include <sstream>
41 :
42 :
43 : namespace libMesh
44 : {
45 : //
46 : // Function definitions for the AutoAreaFunction class
47 : //
48 :
49 : // Constructor
50 0 : AutoAreaFunction::AutoAreaFunction (const Parallel::Communicator &comm,
51 : const unsigned int num_nearest_pts,
52 : const unsigned int power,
53 : const Real background_value,
54 0 : const Real background_eff_dist):
55 0 : _comm(comm),
56 0 : _num_nearest_pts(num_nearest_pts),
57 0 : _power(power),
58 0 : _background_value(background_value),
59 0 : _background_eff_dist(background_eff_dist),
60 0 : _auto_area_mfi(std::make_unique<InverseDistanceInterpolation<3>>(_comm, _num_nearest_pts, _power, _background_value, _background_eff_dist))
61 : {
62 0 : this->_initialized = false;
63 0 : this->_is_time_dependent = false;
64 0 : }
65 :
66 : // Destructor
67 0 : AutoAreaFunction::~AutoAreaFunction () = default;
68 :
69 0 : void AutoAreaFunction::init_mfi (const std::vector<Point> & input_pts,
70 : const std::vector<Real> & input_vals)
71 : {
72 0 : std::vector<std::string> field_vars{"f"};
73 0 : _auto_area_mfi->set_field_variables(field_vars);
74 0 : _auto_area_mfi->get_source_points() = input_pts;
75 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
76 : std::vector<Number> input_complex_vals;
77 0 : for (const auto & input_val : input_vals)
78 0 : input_complex_vals.push_back(Complex (input_val, 0.0));
79 0 : _auto_area_mfi->get_source_vals() = input_complex_vals;
80 : #else
81 0 : _auto_area_mfi->get_source_vals() = input_vals;
82 : #endif
83 0 : _auto_area_mfi->prepare_for_use();
84 0 : this->_initialized = true;
85 0 : }
86 :
87 0 : Real AutoAreaFunction::operator() (const Point & p,
88 : const Real /*time*/)
89 : {
90 0 : libmesh_assert(this->_initialized);
91 :
92 0 : std::vector<Point> target_pts;
93 0 : std::vector<Number> target_vals;
94 :
95 0 : target_pts.push_back(p);
96 0 : target_vals.resize(1);
97 :
98 0 : _auto_area_mfi->interpolate_field_data(_auto_area_mfi->field_variables(), target_pts, target_vals);
99 :
100 0 : return libmesh_real(target_vals.front());
101 : }
102 :
103 : //
104 : // Function definitions for the TriangulatorInterface class
105 : //
106 :
107 : // Constructor
108 2119 : TriangulatorInterface::TriangulatorInterface(UnstructuredMesh & mesh)
109 1943 : : _mesh(mesh),
110 1943 : _holes(nullptr),
111 1943 : _markers(nullptr),
112 1943 : _regions(nullptr),
113 1943 : _elem_type(TRI3),
114 1943 : _desired_area(0.1),
115 1943 : _minimum_angle(20.0),
116 1943 : _triangulation_type(GENERATE_CONVEX_HULL),
117 1943 : _insert_extra_points(false),
118 1943 : _smooth_after_generating(true),
119 1943 : _quiet(true),
120 1943 : _fixup_tri7_center_nodes(false),
121 2295 : _auto_area_function(nullptr)
122 2119 : {}
123 :
124 :
125 221 : void TriangulatorInterface::set_interpolate_boundary_points (int n_points)
126 : {
127 : // Maybe we'll reserve a meaning for negatives later?
128 10 : libmesh_assert(n_points >= 0);
129 :
130 221 : _interpolate_boundary_points = n_points;
131 :
132 : // backwards compatibility - someone (including us) might want to
133 : // query this via the old API.
134 221 : _insert_extra_points = n_points;
135 221 : }
136 :
137 :
138 :
139 2279 : int TriangulatorInterface::get_interpolate_boundary_points () const
140 : {
141 : // backwards compatibility - someone might have turned this off via
142 : // the old API
143 2279 : if (!_insert_extra_points)
144 292 : return 0;
145 :
146 221 : return _interpolate_boundary_points;
147 : }
148 :
149 :
150 :
151 2492 : void TriangulatorInterface::elems_to_segments()
152 : {
153 : // Don't try to override manually specified segments
154 2492 : if (!this->segments.empty())
155 4 : return;
156 :
157 : // If we have edges, they should form the polyline with the ordering
158 : // we want. Let's turn them into segments for later use, because
159 : // we're going to delete the original elements to replace with our
160 : // triangulation.
161 2417 : if (_mesh.n_elem())
162 : {
163 : // Mapping from points to node ids, to back those out from
164 : // MeshedHole results later
165 56 : std::map<Point, dof_id_type> point_id_map;
166 :
167 7393 : for (Node * node : _mesh.node_ptr_range())
168 : {
169 : // We're not going to support overlapping nodes on the boundary
170 5835 : libmesh_error_msg_if
171 : (point_id_map.count(*node),
172 : "TriangulatorInterface does not support overlapping nodes found at "
173 : << static_cast<Point&>(*node));
174 :
175 5835 : point_id_map.emplace(*node, node->id());
176 737 : }
177 :
178 : // We don't support directly generating Tri6, so for
179 : // compatibility with future stitching we need to be working
180 : // with first-order elements. Let's get rid of any non-vertex
181 : // nodes we just added.
182 8802 : for (Elem * elem : _mesh.element_ptr_range())
183 5480 : for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
184 849 : point_id_map.erase(elem->point(n));
185 :
186 : // We'll steal the ordering calculation from
187 : // the MeshedHole code
188 1608 : const TriangulatorInterface::MeshedHole mh { _mesh, this->_bdy_ids };
189 :
190 : // If we've specified only a subset of the mesh as our outer
191 : // boundary, then we may have nodes that don't actually fall
192 : // inside that boundary. Triangulator code doesn't like Steiner
193 : // points that aren't inside the triangulation domain, so we
194 : // need to get rid of them.
195 : //
196 : // Also, if we're using Edge3 elements to define our outer
197 : // boundary, we're only dealing with their 2 end nodes and we'll
198 : // need to get rid of their central nodes.
199 44 : std::unordered_set<Node *> nodes_to_delete;
200 :
201 6864 : for (Elem * elem : _mesh.element_ptr_range())
202 4699 : for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
203 2312 : nodes_to_delete.insert(elem->node_ptr(n));
204 :
205 580 : if (!this->_bdy_ids.empty())
206 : {
207 4048 : for (auto & node : _mesh.node_ptr_range())
208 1953 : if (!mh.contains(*node))
209 204 : nodes_to_delete.insert(node);
210 : }
211 :
212 : // And now we're done with elements. Delete them lest they have
213 : // dangling pointers to nodes we'll be deleting.
214 580 : _mesh.clear_elems();
215 :
216 : // Make segments from boundary nodes; also make sure we don't
217 : // delete them.
218 580 : const std::size_t np = mh.n_points();
219 2829 : for (auto i : make_range(np))
220 : {
221 2249 : const Point pt = mh.point(i);
222 2249 : const dof_id_type id0 = libmesh_map_find(point_id_map, pt);
223 2249 : nodes_to_delete.erase(_mesh.node_ptr(id0));
224 2249 : const Point next_pt = mh.point((np+i+1)%np);
225 2249 : const dof_id_type id1 = libmesh_map_find(point_id_map, next_pt);
226 2249 : this->segments.emplace_back(id0, id1);
227 3614 : for (auto m : make_range(mh.n_midpoints()))
228 : {
229 1365 : this->segment_midpoints.emplace_back(mh.midpoint(m, i));
230 1365 : this->segment_midpoints_keys.emplace_back(pt);
231 : }
232 : }
233 :
234 3030 : for (Node * node : nodes_to_delete)
235 2450 : _mesh.delete_node(node);
236 :
237 580 : if (this->_verify_hole_boundaries && _holes)
238 0 : this->verify_holes(mh);
239 536 : }
240 : }
241 :
242 :
243 :
244 2279 : void TriangulatorInterface::nodes_to_segments(dof_id_type max_node_id)
245 : {
246 : // Don't try to override manually specified segments, or try to add
247 : // segments if we're doing a convex hull
248 2279 : if (!this->segments.empty() || _triangulation_type != PSLG)
249 256 : return;
250 :
251 1210 : for (auto node_it = _mesh.nodes_begin(),
252 1210 : node_end = _mesh.nodes_end();
253 5828 : node_it != node_end;)
254 : {
255 4664 : Node * node = *node_it;
256 :
257 : // If we're out of boundary nodes, the rest are going to be
258 : // Steiner points or hole points
259 4664 : if (node->id() >= max_node_id)
260 0 : break;
261 :
262 4476 : ++node_it;
263 :
264 9140 : Node * next_node = (node_it == node_end) ?
265 1306 : *_mesh.nodes_begin() : *node_it;
266 :
267 4664 : this->segments.emplace_back(node->id(), next_node->id());
268 : }
269 :
270 1164 : if (this->_verify_hole_boundaries && _holes)
271 : {
272 48 : std::vector<Point> outer_pts;
273 3255 : for (auto segment : this->segments)
274 2604 : outer_pts.push_back(_mesh.point(segment.first));
275 :
276 699 : ArbitraryHole ah(outer_pts);
277 651 : this->verify_holes(ah);
278 603 : }
279 : }
280 :
281 :
282 :
283 2279 : void TriangulatorInterface::insert_any_extra_boundary_points()
284 : {
285 : // If the initial PSLG is really simple, e.g. an L-shaped domain or
286 : // a square/rectangle, the resulting triangulation may be very
287 : // "structured" looking. Sometimes this is a problem if your
288 : // intention is to work with an "unstructured" looking grid. We can
289 : // attempt to work around this limitation by inserting midpoints
290 : // into the original PSLG. Inserting additional points into a
291 : // set of points meant to be a convex hull usually makes less sense.
292 :
293 2279 : const int n_interpolated = this->get_interpolate_boundary_points();
294 2279 : if ((_triangulation_type==PSLG) && n_interpolated)
295 : {
296 : // If we were lucky enough to start with contiguous node ids,
297 : // let's keep them that way.
298 221 : dof_id_type nn = _mesh.max_node_id();
299 :
300 : std::vector<std::pair<unsigned int, unsigned int>> old_segments =
301 231 : std::move(this->segments);
302 :
303 : // We expect to have converted any elems and/or nodes into
304 : // segments by now.
305 10 : libmesh_assert(!old_segments.empty());
306 :
307 10 : this->segments.clear();
308 :
309 : // Insert a new point on each segment at evenly spaced locations
310 : // between existing boundary points.
311 : // np=index into new points vector
312 : // n =index into original points vector
313 1105 : for (auto old_segment : old_segments)
314 : {
315 884 : Node * begin_node = _mesh.node_ptr(old_segment.first);
316 884 : Node * end_node = _mesh.node_ptr(old_segment.second);
317 884 : dof_id_type current_id = begin_node->id();
318 2920 : for (auto i : make_range(n_interpolated))
319 : {
320 : // new points are equispaced along the original segments
321 : const Point new_point =
322 2036 : ((n_interpolated-i) * *(Point *)(begin_node) +
323 2036 : (i+1) * *(Point *)(end_node)) /
324 2116 : (n_interpolated + 1);
325 2036 : Node * next_node = _mesh.add_point(new_point, nn++);
326 1876 : this->segments.emplace_back(current_id,
327 2036 : next_node->id());
328 2036 : current_id = next_node->id();
329 : }
330 804 : this->segments.emplace_back(current_id,
331 884 : end_node->id());
332 : }
333 : }
334 2279 : }
335 :
336 :
337 1783 : void TriangulatorInterface::increase_triangle_order()
338 : {
339 1783 : switch (_elem_type)
340 : {
341 42 : case TRI3:
342 : // Nothing to do if we're not requested to increase order
343 1491 : return;
344 221 : case TRI6:
345 221 : _mesh.all_second_order();
346 10 : break;
347 71 : case TRI7:
348 71 : _mesh.all_complete_order();
349 2 : break;
350 0 : default:
351 0 : libmesh_not_implemented();
352 : }
353 :
354 : // If we have any midpoint location data, we'll want to look it up
355 : // by point. all_midpoints[{p, m}] will be the mth midpoint
356 : // location following after point p (when traversing a triangle
357 : // counter-clockwise)
358 14 : std::map<std::pair<Point, unsigned int>, Point> all_midpoints;
359 : unsigned int n_midpoints =
360 316 : this->segment_midpoints.size() / this->segments.size();
361 12 : libmesh_assert_equal_to(this->segments.size() * n_midpoints,
362 : this->segment_midpoints.size());
363 509 : for (auto m : make_range(n_midpoints))
364 1014 : for (auto i : make_range(this->segments.size()))
365 : {
366 797 : const Point & p = segment_midpoints_keys[i*n_midpoints+m];
367 864 : all_midpoints[{p,m}] =
368 60 : this->segment_midpoints[i*n_midpoints+m];
369 : }
370 :
371 292 : if (_holes)
372 150 : for (const Hole * hole : *_holes)
373 : {
374 75 : if (!hole->n_midpoints())
375 0 : continue;
376 75 : if (!n_midpoints)
377 75 : n_midpoints = hole->n_midpoints();
378 0 : else if (hole->n_midpoints() != n_midpoints)
379 0 : libmesh_not_implemented_msg
380 : ("Differing boundary midpoint counts " <<
381 : hole->n_midpoints() << " and " << n_midpoints);
382 :
383 : // Our inner holes are expected to have points in
384 : // counter-clockwise order, which is backwards from how we
385 : // want to traverse them when iterating in counter-clockwise
386 : // order over a triangle, so we'll need to reverse our maps
387 : // carefully here.
388 75 : const auto n_hole_points = hole->n_points();
389 4 : libmesh_assert(n_hole_points);
390 150 : for (auto m : make_range(n_midpoints))
391 : {
392 600 : for (auto i : make_range(n_hole_points-1))
393 : {
394 525 : const Point & p = hole->point(i+1);
395 525 : all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, i);
396 : }
397 75 : const Point & p = hole->point(0);
398 77 : all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, n_hole_points-1);
399 : }
400 : }
401 :
402 : // The n_midpoints > 1 case is for future proofing, but in the
403 : // present we have EDGE4 and no TRI10 yet.
404 292 : if (n_midpoints > 1)
405 0 : libmesh_not_implemented_msg
406 : ("Cannot construct triangles with more than 1 midpoint per edge");
407 :
408 292 : if (!n_midpoints)
409 0 : return;
410 :
411 3036 : for (Elem * elem : _mesh.element_ptr_range())
412 : {
413 : // This should only be called right after we've finished
414 : // converting a triangulation to higher order
415 62 : libmesh_assert_equal_to(elem->n_vertices(), 3);
416 62 : libmesh_assert_not_equal_to(elem->default_order(), FIRST);
417 :
418 5052 : for (auto n : make_range(3))
419 : {
420 : // Only hole/outer boundary segments need adjusted midpoints
421 3975 : if (elem->neighbor_ptr(n))
422 1984 : continue;
423 :
424 156 : const Point & p = elem->point(n);
425 :
426 1697 : if (const auto it = all_midpoints.find({p,0});
427 78 : it != all_midpoints.end())
428 1459 : elem->point(n+3) = it->second;
429 : }
430 268 : }
431 :
432 : // Moving boundary mid-edge nodes can displace the TRI7 interior node
433 : // and tangle the element map. Repositioning the interior node is
434 : // opt-in (off by default); the validity check always runs.
435 292 : if (_elem_type == TRI7 && _fixup_tri7_center_nodes)
436 71 : this->fixup_tri7_center_nodes();
437 :
438 292 : this->verify_quadratic_elements();
439 : }
440 :
441 :
442 71 : void TriangulatorInterface::fixup_tri7_center_nodes()
443 : {
444 2 : libmesh_assert_equal_to(_elem_type, TRI7);
445 :
446 : // Place the interior node at the image of the reference centroid
447 : // (xi, eta) = (1/3, 1/3) under the curved Tri6 map, using the Tri6
448 : // shape function values there as weights: -1/9 on the vertices and
449 : // 4/9 on the mid-edges. This reduces to the straight-edge centroid
450 : // when no boundary midpoint has moved.
451 : static const Real wv = -Real(1)/9;
452 : static const Real wm = Real(4)/9;
453 :
454 416 : for (Elem * elem : _mesh.element_ptr_range())
455 : {
456 4 : libmesh_assert_equal_to(elem->n_vertices(), 3);
457 4 : libmesh_assert_equal_to(elem->n_nodes(), 7u);
458 :
459 142 : elem->point(6) = wv * (elem->point(0) +
460 4 : elem->point(1) +
461 8 : elem->point(2)) +
462 0 : wm * (elem->point(3) +
463 8 : elem->point(4) +
464 16 : elem->point(5));
465 67 : }
466 71 : }
467 :
468 :
469 292 : void TriangulatorInterface::verify_quadratic_elements()
470 : {
471 292 : if (_elem_type != TRI6 && _elem_type != TRI7)
472 0 : return;
473 :
474 : // Once fixup_tri7_center_nodes() has placed node 6, the TRI6 and TRI7
475 : // mappings coincide and this Tri6 formula serves both.
476 : static const Real xi_samples[7] = {Real(0), Real(1), Real(0),
477 : Real(1)/2, Real(1)/2, Real(0),
478 : Real(1)/3};
479 : static const Real eta_samples[7] = {Real(0), Real(0), Real(1),
480 : Real(0), Real(1)/2, Real(1)/2,
481 : Real(1)/3};
482 :
483 2909 : for (Elem * elem : _mesh.element_ptr_range())
484 : {
485 62 : libmesh_assert_equal_to(elem->n_vertices(), 3);
486 62 : libmesh_assert_greater_equal(elem->n_nodes(), 6u);
487 :
488 124 : const Point & x0 = elem->point(0);
489 62 : const Point & x1 = elem->point(1);
490 62 : const Point & x2 = elem->point(2);
491 62 : const Point & x3 = elem->point(3);
492 62 : const Point & x4 = elem->point(4);
493 62 : const Point & x5 = elem->point(5);
494 :
495 : // Tri6 mapping derivative coefficients (see Tri6::volume()):
496 : // dx/dxi = xi*a1 + eta*b1 + c1, dx/deta = xi*b1 + eta*b2 + c2.
497 62 : const Point a1 = 4*x0 + 4*x1 - 8*x3;
498 62 : const Point b1 = 4*x0 - 4*x3 + 4*x4 - 4*x5;
499 62 : const Point c1 = -3*x0 - 1*x1 + 4*x3;
500 62 : const Point b2 = 4*x0 + 4*x2 - 8*x5;
501 62 : const Point c2 = -3*x0 - 1*x2 + 4*x5;
502 :
503 : // Scale the tolerance by the straight-edge triangle area, which
504 : // is strictly positive for the valid TRI3 poly2tri input.
505 1263 : const Real ref_area = 0.5 * cross_norm(x1 - x0, x2 - x0);
506 1263 : const Real jac_tol = TOLERANCE * ref_area;
507 :
508 62 : Real min_jac = std::numeric_limits<Real>::max();
509 62 : unsigned int worst_sample = 0;
510 10104 : for (unsigned int s = 0; s != 7; ++s)
511 : {
512 8841 : const Real xi = xi_samples[s];
513 8841 : const Real eta = eta_samples[s];
514 434 : const Point dxi = xi*a1 + eta*b1 + c1;
515 434 : const Point deta = xi*b1 + eta*b2 + c2;
516 : // z-component of the cross product; the elements are planar.
517 8841 : const Real jac = dxi(0)*deta(1) - dxi(1)*deta(0);
518 8841 : if (jac < min_jac)
519 : {
520 100 : min_jac = jac;
521 100 : worst_sample = s;
522 : }
523 : }
524 :
525 1263 : if (min_jac > jac_tol)
526 60 : continue;
527 :
528 : // Build a diagnostic naming every snapped boundary side on this
529 : // element so the user can immediately see which curved-boundary
530 : // input caused the tangle.
531 75 : std::ostringstream sides;
532 284 : for (unsigned int n = 0; n != 3; ++n)
533 219 : if (!elem->neighbor_ptr(n))
534 : {
535 : const Point straight =
536 213 : 0.5 * (elem->point(n) + elem->point((n+1) % 3));
537 207 : sides << " (boundary side " << n
538 207 : << ": straight midpoint " << straight
539 219 : << ", snapped midpoint " << elem->point(n+3) << ")";
540 : }
541 :
542 485 : libmesh_error_msg(
543 : "TriangulatorInterface: snapping a boundary midpoint produced a "
544 : "tangled quadratic triangle (element " << elem->id()
545 : << ", non-positive Jacobian " << min_jac
546 : << " at reference sample (" << xi_samples[worst_sample] << ", "
547 : << eta_samples[worst_sample] << "); reference triangle area "
548 : << ref_area << ")." << sides.str()
549 : << " Refine the boundary discretization so that recorded "
550 : "midpoints lie closer to their straight-line midpoints, "
551 : "then retry.");
552 335 : }
553 : }
554 :
555 :
556 651 : void TriangulatorInterface::verify_holes(const Hole & outer_bdy)
557 : {
558 1870 : for (const Hole * hole : *_holes)
559 : {
560 7550 : for (const Hole * hole2 : *_holes)
561 : {
562 6331 : if (hole == hole2)
563 1179 : continue;
564 :
565 25560 : for (auto i : make_range(hole2->n_points()))
566 20448 : if (hole->contains(hole2->point(i)))
567 0 : libmesh_error_msg
568 : ("Found point " << hole2->point(i) <<
569 : " on one hole boundary and another's interior");
570 : }
571 :
572 6695 : for (auto i : make_range(hole->n_points()))
573 5476 : if (!outer_bdy.contains(hole->point(i)))
574 0 : libmesh_error_msg
575 : ("Found point " << hole->point(i) <<
576 : " on hole boundary but outside outer boundary");
577 : }
578 651 : }
579 :
580 :
581 504 : unsigned int TriangulatorInterface::total_hole_points()
582 : {
583 : // If the holes vector is non-nullptr (and non-empty) we need to determine
584 : // the number of additional points which the holes will add to the
585 : // triangulation.
586 : // Note that the number of points is always equal to the number of segments
587 : // that form the holes.
588 252 : unsigned int n_hole_points = 0;
589 :
590 504 : if (_holes)
591 40 : for (const auto & hole : *_holes)
592 : {
593 24 : n_hole_points += hole->n_points();
594 : // A hole at least has one enclosure.
595 : // Points on enclosures are ordered so that we can add segments implicitly.
596 : // Elements in segment_indices() indicates the starting points of all enclosures.
597 : // The last element in segment_indices() is the number of total points.
598 12 : libmesh_assert_greater(hole->segment_indices().size(), 1);
599 12 : libmesh_assert_equal_to(hole->segment_indices().back(), hole->n_points());
600 : }
601 :
602 504 : return n_hole_points;
603 : }
604 :
605 0 : void TriangulatorInterface::set_auto_area_function(const Parallel::Communicator &comm,
606 : const unsigned int num_nearest_pts,
607 : const unsigned int power,
608 : const Real background_value,
609 : const Real background_eff_dist)
610 : {
611 0 : _auto_area_function = std::make_unique<AutoAreaFunction>(comm, num_nearest_pts, power, background_value, background_eff_dist);
612 0 : }
613 :
614 0 : FunctionBase<Real> * TriangulatorInterface::get_auto_area_function()
615 : {
616 0 : if (!_auto_area_function->initialized())
617 : {
618 : // Points and target element sizes for the interpolation
619 0 : std::vector<Point> function_points;
620 0 : std::vector<Real> function_sizes;
621 0 : calculate_auto_desired_area_samples(function_points, function_sizes);
622 0 : _auto_area_function->init_mfi(function_points, function_sizes);
623 : }
624 0 : return _auto_area_function.get();
625 : }
626 :
627 0 : void TriangulatorInterface::calculate_auto_desired_area_samples(std::vector<Point> & function_points,
628 : std::vector<Real> & function_sizes,
629 : const Real & area_factor)
630 : {
631 : // Get the hole mesh of the outer boundary
632 : // Holes should already be attached if applicable when this function is called
633 0 : const TriangulatorInterface::MeshedHole bdry_mh { _mesh, this->_bdy_ids };
634 : // Collect all the centroid points of the outer boundary segments
635 : // and the corresponding element sizes
636 0 : for (unsigned int i = 0; i < bdry_mh.n_points(); i++)
637 : {
638 0 : function_points.push_back((bdry_mh.point(i) + bdry_mh.point((i + 1) % bdry_mh.n_points())) /
639 0 : Real(2.0));
640 0 : function_sizes.push_back(
641 0 : (bdry_mh.point(i) - bdry_mh.point((i + 1) % bdry_mh.n_points())).norm());
642 : }
643 : // If holes are present, do the same for the hole boundaries
644 0 : if(_holes)
645 0 : for (const Hole * hole : *_holes)
646 : {
647 0 : for (unsigned int i = 0; i < hole->n_points(); i++)
648 : {
649 0 : function_points.push_back(
650 0 : (hole->point(i) + hole->point((i + 1) % hole->n_points())) / Real(2.0));
651 0 : function_sizes.push_back(
652 0 : (hole->point(i) - hole->point((i + 1) % hole->n_points())).norm());
653 : }
654 : }
655 :
656 0 : std::for_each(
657 0 : function_sizes.begin(), function_sizes.end(), [&area_factor](Real & a) { a = a * a * area_factor * std::sqrt(3.0) / 4.0; });
658 :
659 0 : }
660 : } // namespace libMesh
661 :
|
