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