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