libMesh
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | Private Attributes | List of all members
libMesh::Poly2TriTriangulator Class Reference

A C++ interface between LibMesh and the poly2tri library, with custom code for Steiner point insertion. More...

#include <poly2tri_triangulator.h>

Inheritance diagram for libMesh::Poly2TriTriangulator:
[legend]

Public Types

enum  TriangulationType { GENERATE_CONVEX_HULL = 0 , PSLG = 1 , INVALID_TRIANGULATION_TYPE }
 The TriangulationType is used with the general triangulate function defined below. More...
 

Public Member Functions

 Poly2TriTriangulator (UnstructuredMesh &mesh, dof_id_type n_boundary_nodes=DofObject::invalid_id)
 The constructor.
 
virtual ~Poly2TriTriangulator ()
 Empty destructor.
 
virtual void triangulate () override
 Internally, this calls the poly2tri triangulation code in a loop, inserting our owner Steiner points as necessary to promote mesh quality.
 
virtual void set_desired_area_function (FunctionBase< Real > *desired) override
 Set a function giving desired triangle area as a function of position.
 
virtual FunctionBase< Real > * get_desired_area_function () override
 Get the function giving desired triangle area as a function of position, or nullptr if no such function has been set.
 
virtual void set_refine_boundary_allowed (bool refine_bdy_allowed) override
 Set whether or not the triangulation is allowed to refine the mesh boundary when refining the interior.
 
virtual bool refine_boundary_allowed () const override
 Get whether or not the triangulation is allowed to refine the mesh boundary when refining the interior.
 
ElemTypeelem_type ()
 Sets and/or gets the desired element type.
 
Realdesired_area ()
 Sets and/or gets the desired triangle area.
 
Realminimum_angle ()
 Sets and/or gets the minimum desired angle.
 
TriangulationTypetriangulation_type ()
 Sets and/or gets the desired triangulation type.
 
bool & insert_extra_points ()
 Sets and/or gets the flag for inserting add'l points.
 
void set_interpolate_boundary_points (int n_points)
 Complicated setter, for compatibility with insert_extra_points()
 
int get_interpolate_boundary_points () const
 Complicated getter, for compatibility with insert_extra_points()
 
bool & smooth_after_generating ()
 Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the grid.
 
bool & quiet ()
 Whether not to silence internal messages to stdout.
 
void attach_hole_list (const std::vector< Hole * > *holes)
 Attaches a vector of Hole* pointers which will be meshed around.
 
void set_verify_hole_boundaries (bool v)
 Verifying that hole boundaries don't cross the outer boundary or each other is something like O(N_bdys^2*N_points_per_bdy^2), so we only do it if requested.
 
bool get_verify_hole_boundaries () const
 
void set_fixup_tri7_center_nodes (bool v)
 For TRI7 elements, after boundary mid-edge nodes have been snapped to a curved boundary in increase_triangle_order(), the interior node can be repositioned to the curved-mapping image of the reference centroid.
 
bool get_fixup_tri7_center_nodes () const
 
void attach_boundary_marker (const std::vector< int > *markers)
 Attaches boundary markers.
 
void attach_region_list (const std::vector< Region * > *regions)
 Attaches regions for using attribute to set subdomain IDs and better controlling the triangle sizes within the regions.
 
void set_auto_area_function (const Parallel::Communicator &comm, const unsigned int num_nearest_pts, const unsigned int power, const Real background_value, const Real background_eff_dist)
 Generate an auto area function based on spacing of boundary points.
 
bool has_auto_area_function ()
 Whether or not an auto area function has been set.
 
FunctionBase< Real > * get_auto_area_function ()
 Get the auto area function.
 
void calculate_auto_desired_area_samples (std::vector< Point > &function_points, std::vector< Real > &function_sizes, const Real &area_factor=1.5)
 The external boundary and all hole boundaries are collected.
 
void set_outer_boundary_ids (std::set< std::size_t > bdy_ids)
 A set of ids to allow on the outer boundary loop: interpreted as boundary ids of 2D elements and/or subdomain ids of 1D edges.
 
const std::set< std::size_t > & get_outer_boundary_ids () const
 

Public Attributes

std::vector< std::pair< unsigned int, unsigned int > > segments
 When constructing a PSLG, if the node numbers do not define the desired boundary segments implicitly through the ordering of the points, you can use the segments vector to specify the segments explicitly, Ex: unit square numbered counter-clockwise starting from origin segments[0] = (0,1) segments[1] = (1,2) segments[2] = (2,3) segments[3] = (3,0) (For the above case you could actually use the implicit ordering!)
 
std::vector< Pointsegment_midpoints
 When constructing a second-order triangulation from a second-order boundary, we may do the triangulation using first-order elements, in which case we need to save midpoint location data in order to reconstruct curvature along boundaries.
 
std::vector< Pointsegment_midpoints_keys
 When saving the midpoint location data, we need to save the corresponding segment information too.
 

Protected Member Functions

bool is_refine_boundary_allowed (const BoundaryInfo &boundary_info, const Elem &elem, unsigned int side)
 Is refining this element's boundary side allowed?
 
void triangulate_current_points ()
 Triangulate the current mesh and hole points.
 
bool insert_refinement_points ()
 Add Steiner points as new mesh nodes, as necessary to refine an existing trangulation.
 
bool should_refine_elem (Elem &elem)
 Returns true if the given element ought to be refined according to current criteria.
 
void elems_to_segments ()
 Helper function to create PSLG segments from our other boundary-defining options (1D mesh edges, 2D mesh boundary sides), if no segments already exist.
 
void nodes_to_segments (dof_id_type max_node_id)
 Helper function to create PSLG segments from our node ordering, up to the maximum node id, if no segments already exist.
 
void insert_any_extra_boundary_points ()
 Helper function to add extra points (midpoints of initial segments) to a PSLG triangulation.
 
void increase_triangle_order ()
 Helper function to upconvert Tri3 to any higher order triangle type if requested via _elem_type.
 
void fixup_tri7_center_nodes ()
 Helper called at the end of increase_triangle_order() for TRI7 to reposition the interior node to the curved-mapping image of the reference centroid after boundary midpoints have been snapped.
 
void verify_quadratic_elements ()
 Helper called at the end of increase_triangle_order() to verify that no boundary-midpoint snap has tangled a quadratic triangle, by sampling the signed element Jacobian at the reference nodes and centroid.
 
void verify_holes (const Hole &outer_bdy)
 Helper function to check holes for intersections if requested.
 
unsigned int total_hole_points ()
 Helper function to count points in and verify holes.
 

Protected Attributes

UnstructuredMesh_mesh
 Reference to the mesh which is to be created by triangle.
 
const std::vector< Hole * > * _holes
 A pointer to a vector of Hole*s.
 
const std::vector< int > * _markers
 Boundary markers.
 
const std::vector< Region * > * _regions
 A pointer to a vector of Regions*s.
 
std::set< std::size_t > _bdy_ids
 A set of ids to allow on the outer boundary loop.
 
ElemType _elem_type
 The type of elements to generate.
 
Real _desired_area
 The desired area for the elements in the resulting mesh.
 
Real _minimum_angle
 Minimum angle in triangles.
 
TriangulationType _triangulation_type
 The type of triangulation to perform: choices are: convex hull PSLG.
 
bool _insert_extra_points
 Flag which tells whether or not to insert additional nodes before triangulation.
 
int _interpolate_boundary_points
 Flag which tells how many additional nodes should be inserted between each pair of original mesh points.
 
bool _smooth_after_generating
 Flag which tells whether we should smooth the mesh after it is generated.
 
bool _quiet
 Flag which tells if we want to suppress stdout outputs.
 
bool _verify_hole_boundaries
 Flag which tells if we want to check hole geometry.
 
bool _fixup_tri7_center_nodes
 Flag which tells if we want to relocate TRI7 interior nodes to the curved-mapping centroid after snapping boundary midpoints.
 
std::unique_ptr< AutoAreaFunction_auto_area_function
 The auto area function based on the spacing of boundary points.
 

Private Attributes

std::map< const Hole *, std::unique_ptr< ArbitraryHole > > replaced_holes
 We might have to replace the user-provided holes with refined versions.
 
dof_id_type _n_boundary_nodes
 Keep track of how many mesh nodes are boundary nodes.
 
std::unique_ptr< FunctionBase< Real > > _desired_area_func
 Location-dependent area requirements.
 
bool _refine_bdy_allowed
 Whether to allow boundary refinement.
 

Detailed Description

A C++ interface between LibMesh and the poly2tri library, with custom code for Steiner point insertion.

Author
Roy H. Stogner
Date
2022

Definition at line 46 of file poly2tri_triangulator.h.

Member Enumeration Documentation

◆ TriangulationType

The TriangulationType is used with the general triangulate function defined below.

Enumerator
GENERATE_CONVEX_HULL 

First generate a convex hull from the set of points passed in, and then triangulate this set of points.

This is probably the most common type of usage.

PSLG 

Triangulate the interior of a Planar Straight Line Graph, which is defined implicitly by the order of the "points" vector: a straight line is assumed to lie between each successive pair of points, with an additional line joining the final and first points.

Explicitly telling the triangulator to add additional points may be important for this option.

INVALID_TRIANGULATION_TYPE 

Does nothing, used as a "null" value.

Definition at line 111 of file triangulator_interface.h.

112 {
119
130 PSLG = 1,
131
136 };
@ PSLG
Triangulate the interior of a Planar Straight Line Graph, which is defined implicitly by the order of...
@ INVALID_TRIANGULATION_TYPE
Does nothing, used as a "null" value.
@ GENERATE_CONVEX_HULL
First generate a convex hull from the set of points passed in, and then triangulate this set of point...

Constructor & Destructor Documentation

◆ Poly2TriTriangulator()

libMesh::Poly2TriTriangulator::Poly2TriTriangulator ( UnstructuredMesh mesh,
dof_id_type  n_boundary_nodes = DofObject::invalid_id 
)
explicit

The constructor.

A reference to the mesh containing the points which are to be triangulated must be provided. The first n_boundary_nodes are expected to form a closed loop around the mesh domain; any subsequent nodes are expected to be interior nodes or in the middle of (internal hole or external) boundary segments.

If n_boundary_nodes is not supplied or is invalid_id then all mesh points are expected to be boundary polyline points.

Definition at line 358 of file poly2tri_triangulator.C.

361 _n_boundary_nodes(n_boundary_nodes),
363{
364}
bool _refine_bdy_allowed
Whether to allow boundary refinement.
dof_id_type _n_boundary_nodes
Keep track of how many mesh nodes are boundary nodes.
MeshBase & mesh

◆ ~Poly2TriTriangulator()

libMesh::Poly2TriTriangulator::~Poly2TriTriangulator ( )
virtualdefault

Empty destructor.

Defaulted in the .C so we can forward declare unique_ptr contents.

Member Function Documentation

◆ attach_boundary_marker()

void libMesh::TriangulatorInterface::attach_boundary_marker ( const std::vector< int > *  markers)
inlineinherited

Attaches boundary markers.

If segments is set, the number of markers must be equal to the size of segments, otherwise, it is equal to the number of points.

Definition at line 313 of file triangulator_interface.h.

313{ _markers = markers; }
const std::vector< int > * _markers
Boundary markers.

References libMesh::TriangulatorInterface::_markers.

◆ attach_hole_list()

void libMesh::TriangulatorInterface::attach_hole_list ( const std::vector< Hole * > *  holes)
inlineinherited

◆ attach_region_list()

void libMesh::TriangulatorInterface::attach_region_list ( const std::vector< Region * > *  regions)
inlineinherited

Attaches regions for using attribute to set subdomain IDs and better controlling the triangle sizes within the regions.

Definition at line 319 of file triangulator_interface.h.

319{ _regions = regions; }
const std::vector< Region * > * _regions
A pointer to a vector of Regions*s.

References libMesh::TriangulatorInterface::_regions.

◆ calculate_auto_desired_area_samples()

void libMesh::TriangulatorInterface::calculate_auto_desired_area_samples ( std::vector< Point > &  function_points,
std::vector< Real > &  function_sizes,
const Real area_factor = 1.5 
)
inherited

The external boundary and all hole boundaries are collected.

The centroid of each EDGE element is used as the point position and the desired area is calculated as the area of the equilateral triangle with the edge length as the length of the EDGE element times an area_factor (default is 1.5).

Definition at line 631 of file triangulator_interface.C.

634{
635 // Get the hole mesh of the outer boundary
636 // Holes should already be attached if applicable when this function is called
637 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 for (unsigned int i = 0; i < bdry_mh.n_points(); i++)
641 {
642 function_points.push_back((bdry_mh.point(i) + bdry_mh.point((i + 1) % bdry_mh.n_points())) /
643 Real(2.0));
644 function_sizes.push_back(
645 (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 if(_holes)
649 for (const Hole * hole : *_holes)
650 {
651 for (unsigned int i = 0; i < hole->n_points(); i++)
652 {
653 function_points.push_back(
654 (hole->point(i) + hole->point((i + 1) % hole->n_points())) / Real(2.0));
655 function_sizes.push_back(
656 (hole->point(i) - hole->point((i + 1) % hole->n_points())).norm());
657 }
658 }
659
660 std::for_each(
661 function_sizes.begin(), function_sizes.end(), [&area_factor](Real & a) { a = a * a * area_factor * std::sqrt(3.0) / 4.0; });
662
663}
std::set< std::size_t > _bdy_ids
A set of ids to allow on the outer boundary loop.
UnstructuredMesh & _mesh
Reference to the mesh which is to be created by triangle.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

References libMesh::TriangulatorInterface::_bdy_ids, libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, and libMesh::Real.

Referenced by libMesh::TriangulatorInterface::get_auto_area_function().

◆ desired_area()

Real & libMesh::TriangulatorInterface::desired_area ( )
inlineinherited

Sets and/or gets the desired triangle area.

Set to zero to disable area constraint.

If a desired_area_function is set, then desired_area() should be used to set a minimum desired area; this will reduce "false negatives" by suggesting how finely to sample desired_area_function inside large triangles, where ideally the desired_area_function will be satisfied in the triangle interior and not just at the triangle vertices.

Definition at line 177 of file triangulator_interface.h.

177{return _desired_area;}
Real _desired_area
The desired area for the elements in the resulting mesh.

References libMesh::TriangulatorInterface::_desired_area.

Referenced by libMesh::MeshTools::Generation::build_delaunay_square(), MeshTriangulationTest::commonSettings(), insert_refinement_points(), main(), should_refine_elem(), MeshTriangulationTest::testPoly2TriRefinementBase(), MeshTriangulationTest::testTriangulatorHolePerturbed(), MeshTriangulationTest::testTriangulatorHoleTangentPerturbed(), MeshTriangulationTest::testTriangulatorInterp(), and triangulate_domain().

◆ elem_type()

ElemType & libMesh::TriangulatorInterface::elem_type ( )
inlineinherited

◆ elems_to_segments()

void libMesh::TriangulatorInterface::elems_to_segments ( )
protectedinherited

Helper function to create PSLG segments from our other boundary-defining options (1D mesh edges, 2D mesh boundary sides), if no segments already exist.

Definition at line 155 of file triangulator_interface.C.

156{
157 // Don't try to override manually specified segments
158 if (!this->segments.empty())
159 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 if (_mesh.n_elem())
166 {
167 // Mapping from points to node ids, to back those out from
168 // MeshedHole results later
169 std::map<Point, dof_id_type> point_id_map;
170
171 for (Node * node : _mesh.node_ptr_range())
172 {
173 // We're not going to support overlapping nodes on the boundary
174 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 point_id_map.emplace(*node, node->id());
180 }
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 for (Elem * elem : _mesh.element_ptr_range())
187 for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
188 point_id_map.erase(elem->point(n));
189
190 // We'll steal the ordering calculation from
191 // the MeshedHole code
192 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 std::unordered_set<Node *> nodes_to_delete;
204
205 for (Elem * elem : _mesh.element_ptr_range())
206 for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
207 nodes_to_delete.insert(elem->node_ptr(n));
208
209 if (!this->_bdy_ids.empty())
210 {
211 for (auto & node : _mesh.node_ptr_range())
212 if (!mh.contains(*node))
213 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.
219
220 // Make segments from boundary nodes; also make sure we don't
221 // delete them.
222 const std::size_t np = mh.n_points();
223 for (auto i : make_range(np))
224 {
225 const Point pt = mh.point(i);
226 const dof_id_type id0 = libmesh_map_find(point_id_map, pt);
227 nodes_to_delete.erase(_mesh.node_ptr(id0));
228 const Point next_pt = mh.point((np+i+1)%np);
229 const dof_id_type id1 = libmesh_map_find(point_id_map, next_pt);
230 this->segments.emplace_back(id0, id1);
231 for (auto m : make_range(mh.n_midpoints()))
232 {
233 this->segment_midpoints.emplace_back(mh.midpoint(m, i));
234 this->segment_midpoints_keys.emplace_back(pt);
235 }
236 }
237
238 for (Node * node : nodes_to_delete)
239 _mesh.delete_node(node);
240
241 if (this->_verify_hole_boundaries && _holes)
242 this->verify_holes(mh);
243 }
244}
virtual const Node * node_ptr(const dof_id_type i) const =0
virtual dof_id_type n_elem() const =0
virtual void clear_elems()=0
Deletes all the element data that is currently stored.
bool _verify_hole_boundaries
Flag which tells if we want to check hole geometry.
std::vector< std::pair< unsigned int, unsigned int > > segments
When constructing a PSLG, if the node numbers do not define the desired boundary segments implicitly ...
void verify_holes(const Hole &outer_bdy)
Helper function to check holes for intersections if requested.
std::vector< Point > segment_midpoints
When constructing a second-order triangulation from a second-order boundary, we may do the triangulat...
std::vector< Point > segment_midpoints_keys
When saving the midpoint location data, we need to save the corresponding segment information too.
bool contains(std::string_view superstring, std::string_view substring)
Look for a substring within a string.
Definition utility.C:205
uint8_t dof_id_type
Definition id_types.h:67
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition int_range.h:176
const dof_id_type n_nodes
Definition tecplot_io.C:67

References libMesh::TriangulatorInterface::_bdy_ids, libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, libMesh::TriangulatorInterface::_verify_hole_boundaries, libMesh::MeshBase::clear_elems(), libMesh::MeshBase::delete_node(), libMesh::make_range(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::node_ptr(), libMesh::TriangulatorInterface::segment_midpoints, libMesh::TriangulatorInterface::segment_midpoints_keys, libMesh::TriangulatorInterface::segments, and libMesh::TriangulatorInterface::verify_holes().

Referenced by libMesh::TriangleInterface::triangulate(), and triangulate().

◆ fixup_tri7_center_nodes()

void libMesh::TriangulatorInterface::fixup_tri7_center_nodes ( )
protectedinherited

Helper called at the end of increase_triangle_order() for TRI7 to reposition the interior node to the curved-mapping image of the reference centroid after boundary midpoints have been snapped.

Safe to call unconditionally for TRI7.

Definition at line 446 of file triangulator_interface.C.

447{
448 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 for (Elem * elem : _mesh.element_ptr_range())
459 {
460 libmesh_assert_equal_to(elem->n_vertices(), 3);
461 libmesh_assert_equal_to(elem->n_nodes(), 7u);
462
463 elem->point(6) = wv * (elem->point(0) +
464 elem->point(1) +
465 elem->point(2)) +
466 wm * (elem->point(3) +
467 elem->point(4) +
468 elem->point(5));
469 }
470}

References libMesh::TriangulatorInterface::_elem_type, libMesh::TriangulatorInterface::_mesh, libMesh::Real, and libMesh::TRI7.

Referenced by libMesh::TriangulatorInterface::increase_triangle_order().

◆ get_auto_area_function()

FunctionBase< Real > * libMesh::TriangulatorInterface::get_auto_area_function ( )
inherited

Get the auto area function.

Definition at line 618 of file triangulator_interface.C.

619{
620 if (!_auto_area_function->initialized())
621 {
622 // Points and target element sizes for the interpolation
623 std::vector<Point> function_points;
624 std::vector<Real> function_sizes;
625 calculate_auto_desired_area_samples(function_points, function_sizes);
626 _auto_area_function->init_mfi(function_points, function_sizes);
627 }
628 return _auto_area_function.get();
629}
void calculate_auto_desired_area_samples(std::vector< Point > &function_points, std::vector< Real > &function_sizes, const Real &area_factor=1.5)
The external boundary and all hole boundaries are collected.
std::unique_ptr< AutoAreaFunction > _auto_area_function
The auto area function based on the spacing of boundary points.

References libMesh::TriangulatorInterface::_auto_area_function, and libMesh::TriangulatorInterface::calculate_auto_desired_area_samples().

Referenced by should_refine_elem().

◆ get_desired_area_function()

FunctionBase< Real > * libMesh::Poly2TriTriangulator::get_desired_area_function ( )
overridevirtual

Get the function giving desired triangle area as a function of position, or nullptr if no such function has been set.

Reimplemented from libMesh::TriangulatorInterface.

Definition at line 460 of file poly2tri_triangulator.C.

461{
462 return _desired_area_func.get();
463}
std::unique_ptr< FunctionBase< Real > > _desired_area_func
Location-dependent area requirements.

References _desired_area_func.

Referenced by insert_refinement_points(), and should_refine_elem().

◆ get_fixup_tri7_center_nodes()

bool libMesh::TriangulatorInterface::get_fixup_tri7_center_nodes ( ) const
inlineinherited

Definition at line 275 of file triangulator_interface.h.

bool _fixup_tri7_center_nodes
Flag which tells if we want to relocate TRI7 interior nodes to the curved-mapping centroid after snap...

References libMesh::TriangulatorInterface::_fixup_tri7_center_nodes.

◆ get_interpolate_boundary_points()

int libMesh::TriangulatorInterface::get_interpolate_boundary_points ( ) const
inherited

Complicated getter, for compatibility with insert_extra_points()

Definition at line 143 of file triangulator_interface.C.

144{
145 // backwards compatibility - someone might have turned this off via
146 // the old API
148 return 0;
149
151}
int _interpolate_boundary_points
Flag which tells how many additional nodes should be inserted between each pair of original mesh poin...
bool _insert_extra_points
Flag which tells whether or not to insert additional nodes before triangulation.

References libMesh::TriangulatorInterface::_insert_extra_points, and libMesh::TriangulatorInterface::_interpolate_boundary_points.

Referenced by libMesh::TriangulatorInterface::insert_any_extra_boundary_points().

◆ get_outer_boundary_ids()

const std::set< std::size_t > & libMesh::TriangulatorInterface::get_outer_boundary_ids ( ) const
inlineinherited

Definition at line 371 of file triangulator_interface.h.

371{ return _bdy_ids; }

References libMesh::TriangulatorInterface::_bdy_ids.

◆ get_verify_hole_boundaries()

bool libMesh::TriangulatorInterface::get_verify_hole_boundaries ( ) const
inlineinherited

◆ has_auto_area_function()

bool libMesh::TriangulatorInterface::has_auto_area_function ( )
inlineinherited

Whether or not an auto area function has been set.

Definition at line 348 of file triangulator_interface.h.

348{return _auto_area_function != nullptr;}

References libMesh::TriangulatorInterface::_auto_area_function.

Referenced by insert_refinement_points(), and should_refine_elem().

◆ increase_triangle_order()

void libMesh::TriangulatorInterface::increase_triangle_order ( )
protectedinherited

Helper function to upconvert Tri3 to any higher order triangle type if requested via _elem_type.

Should be called at the end of triangulate()

Definition at line 341 of file triangulator_interface.C.

342{
343 switch (_elem_type)
344 {
345 case TRI3:
346 // Nothing to do if we're not requested to increase order
347 return;
348 case TRI6:
350 break;
351 case TRI7:
353 break;
354 default:
355 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 std::map<std::pair<Point, unsigned int>, Point> all_midpoints;
363 unsigned int n_midpoints =
364 this->segment_midpoints.size() / this->segments.size();
365 libmesh_assert_equal_to(this->segments.size() * n_midpoints,
366 this->segment_midpoints.size());
367 for (auto m : make_range(n_midpoints))
368 for (auto i : make_range(this->segments.size()))
369 {
370 const Point & p = segment_midpoints_keys[i*n_midpoints+m];
371 all_midpoints[{p,m}] =
372 this->segment_midpoints[i*n_midpoints+m];
373 }
374
375 if (_holes)
376 for (const Hole * hole : *_holes)
377 {
378 if (!hole->n_midpoints())
379 continue;
380 if (!n_midpoints)
381 n_midpoints = hole->n_midpoints();
382 else if (hole->n_midpoints() != n_midpoints)
383 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 const auto n_hole_points = hole->n_points();
393 libmesh_assert(n_hole_points);
394 for (auto m : make_range(n_midpoints))
395 {
396 for (auto i : make_range(n_hole_points-1))
397 {
398 const Point & p = hole->point(i+1);
399 all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, i);
400 }
401 const Point & p = hole->point(0);
402 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 if (n_midpoints > 1)
409 libmesh_not_implemented_msg
410 ("Cannot construct triangles with more than 1 midpoint per edge");
411
412 if (!n_midpoints)
413 return;
414
415 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 libmesh_assert_equal_to(elem->n_vertices(), 3);
420 libmesh_assert_not_equal_to(elem->default_order(), FIRST);
421
422 for (auto n : make_range(3))
423 {
424 // Only hole/outer boundary segments need adjusted midpoints
425 if (elem->neighbor_ptr(n))
426 continue;
427
428 const Point & p = elem->point(n);
429
430 if (const auto it = all_midpoints.find({p,0});
431 it != all_midpoints.end())
432 elem->point(n+3) = it->second;
433 }
434 }
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.
441
443}
void all_second_order(const bool full_ordered=true)
Calls the range-based version of this function with a range consisting of all elements in the mesh.
Definition mesh_base.C:1803
virtual void all_complete_order()
Calls the range-based version of this function with a range consisting of all elements in the mesh.
Definition mesh_base.C:1808
void fixup_tri7_center_nodes()
Helper called at the end of increase_triangle_order() for TRI7 to reposition the interior node to the...
void verify_quadratic_elements()
Helper called at the end of increase_triangle_order() to verify that no boundary-midpoint snap has ta...
libmesh_assert(ctx)

References libMesh::TriangulatorInterface::_elem_type, libMesh::TriangulatorInterface::_fixup_tri7_center_nodes, libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, libMesh::MeshBase::all_complete_order(), libMesh::MeshBase::all_second_order(), libMesh::FIRST, libMesh::TriangulatorInterface::fixup_tri7_center_nodes(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::TriangulatorInterface::segment_midpoints, libMesh::TriangulatorInterface::segment_midpoints_keys, libMesh::TriangulatorInterface::segments, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, and libMesh::TriangulatorInterface::verify_quadratic_elements().

Referenced by libMesh::TriangleInterface::triangulate(), and triangulate().

◆ insert_any_extra_boundary_points()

void libMesh::TriangulatorInterface::insert_any_extra_boundary_points ( )
protectedinherited

Helper function to add extra points (midpoints of initial segments) to a PSLG triangulation.

Definition at line 287 of file triangulator_interface.C.

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 const int n_interpolated = this->get_interpolate_boundary_points();
298 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.
303
304 std::vector<std::pair<unsigned int, unsigned int>> old_segments =
305 std::move(this->segments);
306
307 // We expect to have converted any elems and/or nodes into
308 // segments by now.
309 libmesh_assert(!old_segments.empty());
310
311 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 for (auto old_segment : old_segments)
318 {
319 Node * begin_node = _mesh.node_ptr(old_segment.first);
320 Node * end_node = _mesh.node_ptr(old_segment.second);
321 dof_id_type current_id = begin_node->id();
322 for (auto i : make_range(n_interpolated))
323 {
324 // new points are equispaced along the original segments
325 const Point new_point =
326 ((n_interpolated-i) * *(Point *)(begin_node) +
327 (i+1) * *(Point *)(end_node)) /
328 (n_interpolated + 1);
329 Node * next_node = _mesh.add_point(new_point, nn++);
330 this->segments.emplace_back(current_id,
331 next_node->id());
332 current_id = next_node->id();
333 }
334 this->segments.emplace_back(current_id,
335 end_node->id());
336 }
337 }
338}
dof_id_type id() const
Definition dof_object.h:819
virtual dof_id_type max_node_id() const =0
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
int get_interpolate_boundary_points() const
Complicated getter, for compatibility with insert_extra_points()
TriangulationType _triangulation_type
The type of triangulation to perform: choices are: convex hull PSLG.

References libMesh::TriangulatorInterface::_mesh, libMesh::TriangulatorInterface::_triangulation_type, libMesh::MeshBase::add_point(), libMesh::TriangulatorInterface::get_interpolate_boundary_points(), libMesh::DofObject::id(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::max_node_id(), libMesh::MeshBase::node_ptr(), libMesh::TriangulatorInterface::PSLG, and libMesh::TriangulatorInterface::segments.

Referenced by libMesh::TriangleInterface::triangulate(), and triangulate().

◆ insert_extra_points()

bool & libMesh::TriangulatorInterface::insert_extra_points ( )
inlineinherited

Sets and/or gets the flag for inserting add'l points.

Definition at line 210 of file triangulator_interface.h.

210{return _insert_extra_points;}

References libMesh::TriangulatorInterface::_insert_extra_points.

◆ insert_refinement_points()

bool libMesh::Poly2TriTriangulator::insert_refinement_points ( )
protected

Add Steiner points as new mesh nodes, as necessary to refine an existing trangulation.

Returns true iff new points were added.

Definition at line 786 of file poly2tri_triangulator.C.

787{
788 LOG_SCOPE("insert_refinement_points()", "Poly2TriTriangulator");
789
790 if (this->minimum_angle() != 0)
791 libmesh_not_implemented();
792
793 // We need neighbor pointers for ray casting and cavity finding
794 UnstructuredMesh & mesh = dynamic_cast<UnstructuredMesh &>(this->_mesh);
795 mesh.find_neighbors();
796
797 if (this->desired_area() == 0 &&
798 this->get_desired_area_function() == nullptr &&
799 !this->has_auto_area_function())
800 return false;
801
802 BoundaryInfo & boundary_info = _mesh.get_boundary_info();
803
804 // We won't immediately add these, lest we invalidate iterators on a
805 // ReplicatedMesh. They'll still be in the mesh neighbor topology
806 // for the purpose of doing Delaunay cavity stuff, so we need to
807 // manage memory here, but there's no point in adding them to the
808 // Mesh just to remove them again afterward when we hit up poly2tri.
809
810 // We'll need to be able to remove new elems from new_elems, in
811 // cases where a later refinement insertion has a not-yet-added
812 // element in its cavity, so we'll use a map here to make searching
813 // possible.
814 //
815 // For parallel consistency, we can't order a container we plan to
816 // iterate through based on Elem * or a hash of it. We'll be doing
817 // Delaunay swaps so we can't iterate based on geometry. These are
818 // not-yet-added elements so we can't iterate based on proper
819 // element ids ... but we can set a temporary element id to use for
820 // the purpose.
821 struct cmp {
822 bool operator()(Elem * a, Elem * b) const {
823 libmesh_assert(a == b || a->id() != b->id());
824 return (a->id() < b->id());
825 }
826 } comp;
827
828 std::map<Elem *, std::unique_ptr<Elem>, decltype(comp)> new_elems(comp);
829
830 // We should already be Delaunay when we get here, otherwise we
831 // won't be able to stay Delaunay later. But we're *not* always
832 // Delaunay when we get here? What the hell, poly2tri? Fixing this
833 // is expensive!
834 {
835 // restore_delaunay should get to the same Delaunay triangulation up to
836 // isomorphism regardless of ordering ... but we actually care
837 // about the isomorphisms! If a triangle's nodes are permuted on
838 // one processor vs another that's an issue. So sort our input
839 // carefully.
840 std::set<Elem *, decltype(comp)> all_elems
841 { mesh.elements_begin(), mesh.elements_end(), comp };
842
843 restore_delaunay(all_elems, boundary_info);
844
845 libmesh_assert_delaunay(mesh, new_elems);
846 }
847
848 // Map of which points follow which in the boundary polylines. If
849 // we have to add new boundary points, we'll use this to construct
850 // an updated this->segments to retriangulate with. If we have to
851 // add new hole points, we'll use this to insert points into an
852 // ArbitraryHole.
853 std::unordered_map<Point, Node *> next_boundary_node;
854
855 // In cases where we've been working with contiguous node id ranges;
856 // let's keep it that way.
859
860 // We can't handle duplicated nodes. We shouldn't ever create one,
861 // but let's make sure of that.
862#ifdef DEBUG
863 std::unordered_set<Point> mesh_points;
864 for (const Node * node : mesh.node_ptr_range())
865 {
866 libmesh_assert(!mesh_points.count(*node));
867 mesh_points.insert(*node);
868 }
869#endif
870
871 auto add_point = [&mesh,
872#ifdef DEBUG
873 &mesh_points,
874#endif
875 &nn](const Point & p)
876 {
877#ifdef DEBUG
878 libmesh_assert(!mesh_points.count(p));
879 mesh_points.insert(p);
880#endif
881 return mesh.add_point(p, nn++);
882 };
883
884 for (auto & elem : mesh.element_ptr_range())
885 {
886 // element_ptr_range skips deleted elements ... right?
887 libmesh_assert(elem);
888 libmesh_assert(elem->valid_id());
889
890 // We only handle triangles in our triangulation
891 libmesh_assert_equal_to(elem->level(), 0u);
892 libmesh_assert_equal_to(elem->type(), TRI3);
893
894 // If this triangle is as small as we desire, move along
895 if (!should_refine_elem(*elem))
896 continue;
897
898 // Otherwise add a Steiner point. We'd like to add the
899 // circumcenter ...
900 Point new_pt = elem->quasicircumcenter();
901
902 // And to give it a node;
903 Node * new_node = nullptr;
904
905 // But that might be outside our triangle, or even outside the
906 // boundary. We'll find a triangle that should contain our new
907 // point
908 Elem * cavity_elem = elem; // Start looking at elem anyway
909
910 // We'll refine a boundary later if necessary.
911 auto boundary_refine = [this, &next_boundary_node,
912 &cavity_elem, &new_node]
913 (unsigned int side)
914 {
915 libmesh_ignore(this); // Only used in dbg/devel
916 libmesh_assert(new_node);
917 libmesh_assert(new_node->valid_id());
918
919 Node * old_segment_start = cavity_elem->node_ptr(side),
920 * old_segment_end = cavity_elem->node_ptr((side+1)%3);
921 libmesh_assert(old_segment_start);
922 libmesh_assert(old_segment_start->valid_id());
923 libmesh_assert(old_segment_end);
924 libmesh_assert(old_segment_end->valid_id());
925
926 if (auto it = next_boundary_node.find(*old_segment_start);
927 it != next_boundary_node.end())
928 {
929 libmesh_assert(it->second == old_segment_end);
930 it->second = new_node;
931 }
932 else
933 {
934 // This would be an O(N) sanity check if we already
935 // have a segments vector or any holes. :-P
936 libmesh_assert(!this->segments.empty() ||
937 (_holes && !_holes->empty()) ||
938 (old_segment_end->id() ==
939 old_segment_start->id() + 1));
940 next_boundary_node[*old_segment_start] = new_node;
941 }
942
943 next_boundary_node[*new_node] = old_segment_end;
944 };
945
946 // Let's find a triangle containing our new point, or at least
947 // containing the end of a ray leading from our current triangle
948 // to the new point.
949 Point ray_start = elem->vertex_average();
950
951 // What side are we coming from, and what side are we going to?
952 unsigned int source_side = invalid_uint;
953 unsigned int side = invalid_uint;
954
955 while (!cavity_elem->contains_point(new_pt))
956 {
957 side = segment_intersection(*cavity_elem, ray_start, new_pt, source_side);
958
959 libmesh_assert_not_equal_to (side, invalid_uint);
960
961 Elem * neigh = cavity_elem->neighbor_ptr(side);
962 // If we're on a boundary, stop there. Refine the boundary
963 // if we're allowed, the boundary element otherwise.
964 if (!neigh)
965 {
966 if (this->is_refine_boundary_allowed(boundary_info,
967 *cavity_elem,
968 side))
969 {
970 new_pt = ray_start;
971 new_node = add_point(new_pt);
972 boundary_refine(side);
973 }
974 else
975 {
976 // Should we just add the vertex average of the
977 // boundary element, to minimize the number of
978 // slivers created?
979 //
980 // new_pt = cavity_elem->vertex_average();
981 //
982 // That works for a while, but it
983 // seems to be able to "run away" and leave us with
984 // crazy slivers on boundaries if we push interior
985 // refinement too far while disabling boundary
986 // refinement.
987 //
988 // Let's go back to refining our original problem
989 // element.
990 cavity_elem = elem;
991 new_pt = cavity_elem->vertex_average();
992 new_node = add_point(new_pt);
993
994 // This was going to be a side refinement but it's
995 // now an internal refinement
996 side = invalid_uint;
997 }
998
999 break;
1000 }
1001
1002 source_side = neigh->which_neighbor_am_i(cavity_elem);
1003 cavity_elem = neigh;
1004 side = invalid_uint;
1005 }
1006
1007 // If we're ready to create a new node and we're not on a
1008 // boundary ... should we be? We don't want to create any
1009 // sliver elements or confuse poly2tri or anything.
1010 if (side == invalid_uint && !new_node)
1011 {
1012 unsigned int worst_side = libMesh::invalid_uint;
1013 Real worst_cos = 1;
1014 for (auto s : make_range(3u))
1015 {
1016 // We never snap to a non-domain-boundary
1017 if (cavity_elem->neighbor_ptr(s))
1018 continue;
1019
1020 Real ax = cavity_elem->point(s)(0) - new_pt(0),
1021 ay = cavity_elem->point(s)(1) - new_pt(1),
1022 bx = cavity_elem->point((s+1)%3)(0) - new_pt(0),
1023 by = cavity_elem->point((s+1)%3)(1) - new_pt(1);
1024 const Real my_cos = (ax*bx+ay*by) /
1025 std::sqrt(ax*ax+ay*ay) /
1026 std::sqrt(bx*bx+by*by);
1027
1028 if (my_cos < worst_cos)
1029 {
1030 worst_side = s;
1031 worst_cos = my_cos;
1032 }
1033 }
1034
1035 // If we'd create a sliver element on the side, let's just
1036 // refine the side instead, if we're allowed.
1037 if (worst_cos < -0.6) // -0.5 is the best we could enforce?
1038 {
1039 side = worst_side;
1040
1041 if (this->is_refine_boundary_allowed(boundary_info,
1042 *cavity_elem,
1043 side))
1044 {
1045 // Let's just try bisecting for now
1046 new_pt = (cavity_elem->point(side) +
1047 cavity_elem->point((side+1)%3)) / 2;
1048 new_node = add_point(new_pt);
1049 boundary_refine(side);
1050 }
1051 else // Do the best we can under these restrictions
1052 {
1053 new_pt = cavity_elem->vertex_average();
1054 new_node = add_point(new_pt);
1055
1056 // This was going to be a side refinement but it's
1057 // now an internal refinement
1058 side = invalid_uint;
1059 }
1060 }
1061 else
1062 new_node = add_point(new_pt);
1063 }
1064 else
1065 libmesh_assert(new_node);
1066
1067 // Find the Delaunay cavity around the new point.
1068 std::set<Elem *, decltype(comp)> cavity(comp);
1069
1070 std::set<Elem *, decltype(comp)> unchecked_cavity ({cavity_elem}, comp);
1071 while (!unchecked_cavity.empty())
1072 {
1073 std::set<Elem *, decltype(comp)> checking_cavity(comp);
1074 checking_cavity.swap(unchecked_cavity);
1075 for (Elem * checking_elem : checking_cavity)
1076 {
1077 for (auto s : make_range(3u))
1078 {
1079 Elem * neigh = checking_elem->neighbor_ptr(s);
1080 if (!neigh || checking_cavity.count(neigh) || cavity.count(neigh))
1081 continue;
1082
1083 if (in_circumcircle(*neigh, new_pt, TOLERANCE*TOLERANCE))
1084 unchecked_cavity.insert(neigh);
1085 }
1086 }
1087
1088 libmesh_merge_move(cavity, checking_cavity);
1089 }
1090
1091 // Retriangulate the Delaunay cavity.
1092 // Each of our cavity triangle edges that are exterior to the
1093 // cavity will be a source of one new triangle.
1094
1095 // Set of elements that might need Delaunay swaps
1096 std::set<Elem *, decltype(comp)> check_delaunay_on(comp);
1097
1098 // Keep maps for doing neighbor pointer assignment. Not going
1099 // to iterate through these so hashing pointers is fine.
1100 std::unordered_map<Node *, std::pair<Elem *, boundary_id_type>>
1101 neighbors_CCW, neighbors_CW;
1102
1103 for (Elem * old_elem : cavity)
1104 {
1105 for (auto s : make_range(3u))
1106 {
1107 Elem * neigh = old_elem->neighbor_ptr(s);
1108 if (!neigh || !cavity.count(neigh))
1109 {
1110 Node * node_CW = old_elem->node_ptr(s),
1111 * node_CCW = old_elem->node_ptr((s+1)%3);
1112
1113 auto set_neighbors =
1114 [&neighbors_CW, &neighbors_CCW, &node_CW,
1115 &node_CCW, &boundary_info]
1116 (Elem * new_neigh, boundary_id_type bcid)
1117 {
1118 // Set clockwise neighbor and vice-versa if possible
1119 if (const auto CW_it = neighbors_CW.find(node_CW);
1120 CW_it == neighbors_CW.end())
1121 {
1122 libmesh_assert(!neighbors_CCW.count(node_CW));
1123 neighbors_CCW[node_CW] = std::make_pair(new_neigh, bcid);
1124 }
1125 else
1126 {
1127 Elem * neigh_CW = CW_it->second.first;
1128 if (new_neigh)
1129 {
1130 new_neigh->set_neighbor(0, neigh_CW);
1131 boundary_id_type bcid_CW = CW_it->second.second;
1132 if (bcid_CW != BoundaryInfo::invalid_id)
1133 boundary_info.add_side(new_neigh, 0, bcid_CW);
1134
1135 }
1136 if (neigh_CW)
1137 {
1138 neigh_CW->set_neighbor(2, new_neigh);
1139 if (bcid != BoundaryInfo::invalid_id)
1140 boundary_info.add_side(neigh_CW, 2, bcid);
1141 }
1142 neighbors_CW.erase(CW_it);
1143 }
1144
1145 // Set counter-CW neighbor and vice-versa if possible
1146 if (const auto CCW_it = neighbors_CCW.find(node_CCW);
1147 CCW_it == neighbors_CCW.end())
1148 {
1149 libmesh_assert(!neighbors_CW.count(node_CCW));
1150 neighbors_CW[node_CCW] = std::make_pair(new_neigh, bcid);
1151 }
1152 else
1153 {
1154 Elem * neigh_CCW = CCW_it->second.first;
1155 if (new_neigh)
1156 {
1157 boundary_id_type bcid_CCW = CCW_it->second.second;
1158 new_neigh->set_neighbor(2, neigh_CCW);
1159 if (bcid_CCW != BoundaryInfo::invalid_id)
1160 boundary_info.add_side(new_neigh, 2, bcid_CCW);
1161 }
1162 if (neigh_CCW)
1163 {
1164 neigh_CCW->set_neighbor(0, new_neigh);
1165 if (bcid != BoundaryInfo::invalid_id)
1166 boundary_info.add_side(neigh_CCW, 0, bcid);
1167 }
1168 neighbors_CCW.erase(CCW_it);
1169 }
1170 };
1171
1172 // We aren't going to try to add a sliver element if we
1173 // have a new boundary node here. We do need to
1174 // keep track of other elements' neighbors, though.
1175 if (old_elem == cavity_elem &&
1176 s == side)
1177 {
1178 std::vector<boundary_id_type> bcids;
1179 boundary_info.boundary_ids(old_elem, s, bcids);
1180 libmesh_assert_equal_to(bcids.size(), 1);
1181 set_neighbors(nullptr, bcids[0]);
1182 continue;
1183 }
1184
1185 auto new_elem = Elem::build_with_id(TRI3, ne++);
1186 new_elem->set_node(0, new_node);
1187 new_elem->set_node(1, node_CW);
1188 new_elem->set_node(2, node_CCW);
1189 libmesh_assert(!new_elem->is_flipped());
1190
1191 // Set in-and-out-of-cavity neighbor pointers
1192 new_elem->set_neighbor(1, neigh);
1193 if (neigh)
1194 {
1195 const unsigned int neigh_s =
1196 neigh->which_neighbor_am_i(old_elem);
1197 neigh->set_neighbor(neigh_s, new_elem.get());
1198 }
1199 else
1200 {
1201 std::vector<boundary_id_type> bcids;
1202 boundary_info.boundary_ids(old_elem, s, bcids);
1203 boundary_info.add_side(new_elem.get(), 1, bcids);
1204 }
1205
1206 // Set in-cavity neighbors' neighbor pointers
1207 set_neighbors(new_elem.get(), BoundaryInfo::invalid_id);
1208
1209 // C++ allows function argument evaluation in any
1210 // order, but we need get() to precede move
1211 Elem * new_elem_ptr = new_elem.get();
1212 new_elems.emplace(new_elem_ptr, std::move(new_elem));
1213
1214 check_delaunay_on.insert(new_elem_ptr);
1215 }
1216 }
1217
1218 boundary_info.remove(old_elem);
1219 }
1220
1221 // Now that we're done using our cavity elems (including with a
1222 // cavity.find() that used a comparator that dereferences the
1223 // elements!) it's safe to delete them.
1224 for (Elem * old_elem : cavity)
1225 {
1226 if (const auto it = new_elems.find(old_elem);
1227 it == new_elems.end())
1228 mesh.delete_elem(old_elem);
1229 else
1230 new_elems.erase(it);
1231 }
1232
1233 // Everybody found their match?
1234 libmesh_assert(neighbors_CW.empty());
1235 libmesh_assert(neighbors_CCW.empty());
1236
1237 // Because we're preserving boundaries here, our naive cavity
1238 // triangulation might not be a Delaunay triangulation. Let's
1239 // check and if necessary fix that; we depend on it when doing
1240 // future point insertions.
1241 restore_delaunay(check_delaunay_on, boundary_info);
1242
1243 // This is too expensive to do on every cavity in devel mode
1244#ifdef DEBUG
1245 libmesh_assert_delaunay(mesh, new_elems);
1246#endif
1247 }
1248
1249 // If we added any new boundary nodes, we're going to need to keep
1250 // track of the changes they made to the outer polyline and/or to
1251 // any holes.
1252 if (!next_boundary_node.empty())
1253 {
1254 auto checked_emplace = [this](dof_id_type new_first,
1255 dof_id_type new_second)
1256 {
1257#ifdef DEBUG
1258 for (auto [first, second] : this->segments)
1259 {
1260 libmesh_assert_not_equal_to(first, new_first);
1261 libmesh_assert_not_equal_to(second, new_second);
1262 }
1263 if (!this->segments.empty())
1264 libmesh_assert_equal_to(this->segments.back().second, new_first);
1265#endif
1266 libmesh_assert_not_equal_to(new_first, new_second);
1267
1268 this->segments.emplace_back (new_first, new_second);
1269 };
1270
1271 // Keep track of the outer polyline
1272 if (this->segments.empty())
1273 {
1275
1276 // Custom loop because we increment node_it 1+ times inside
1277 for (auto node_it = _mesh.nodes_begin(),
1278 node_end = _mesh.nodes_end();
1279 node_it != node_end;)
1280 {
1281 Node & node = **node_it;
1282 ++node_it;
1283
1284 const dof_id_type node_id = node.id();
1285
1286 // Don't add Steiner points
1287 if (node_id >= _n_boundary_nodes)
1288 break;
1289
1290 // Connect up the previous node, if we didn't already
1291 // connect it after some newly inserted nodes
1292 if (!this->segments.empty())
1293 last_id = this->segments.back().second;
1294
1295 if (last_id != DofObject::invalid_id &&
1296 last_id != node_id)
1297 checked_emplace(last_id, node_id);
1298
1299 last_id = node_id;
1300
1301 // Connect to any newly inserted nodes
1302 Node * this_node = &node;
1303 auto it = next_boundary_node.find(*this_node);
1304 while (it != next_boundary_node.end())
1305 {
1306 libmesh_assert(this_node->valid_id());
1307 Node * next_node = it->second;
1308 libmesh_assert(next_node->valid_id());
1309
1310 if (node_it != node_end &&
1311 next_node == *node_it)
1312 ++node_it;
1313
1314 checked_emplace(this_node->id(), next_node->id());
1315
1316 this_node = next_node;
1317 if (this_node->id() == this->segments.front().first)
1318 break;
1319
1320 it = next_boundary_node.find(*this_node);
1321 }
1322 }
1323
1324 // We expect a closed loop here
1325 if (this->segments.back().second != this->segments.front().first)
1326 checked_emplace(this->segments.back().second,
1327 this->segments.front().first);
1328 }
1329 else
1330 {
1331 std::vector<std::pair<unsigned int, unsigned int>> old_segments;
1332 old_segments.swap(this->segments);
1333
1334 auto old_it = old_segments.begin();
1335
1336 const Node * node = _mesh.node_ptr(old_it->first);
1337 const Node * const first_node = node;
1338
1339 do
1340 {
1341 const dof_id_type node_id = node->id();
1342 if (const auto it = next_boundary_node.find(*node);
1343 it == next_boundary_node.end())
1344 {
1345 while (node_id != old_it->first)
1346 {
1347 ++old_it;
1348 libmesh_assert(old_it != old_segments.end());
1349 }
1350 node = mesh.node_ptr(old_it->second);
1351 }
1352 else
1353 {
1354 node = it->second;
1355 }
1356
1357 checked_emplace(node_id, node->id());
1358 }
1359 while (node != first_node);
1360 }
1361
1362 // Keep track of any holes
1363 if (this->_holes)
1364 {
1365 // Do we have any holes that need to be newly replaced?
1366 for (const Hole * hole : *this->_holes)
1367 {
1368 if (this->replaced_holes.count(hole))
1369 continue;
1370
1371 bool hole_point_insertion = false;
1372 for (auto p : make_range(hole->n_points()))
1373 if (next_boundary_node.count(hole->point(p)))
1374 {
1375 hole_point_insertion = true;
1376 break;
1377 }
1378 if (hole_point_insertion)
1379 this->replaced_holes.emplace
1380 (hole, std::make_unique<ArbitraryHole>(*hole));
1381 }
1382
1383 // If we have any holes that are being replaced, make sure
1384 // their replacements are up to date.
1385 for (const Hole * hole : *this->_holes)
1386 {
1387 auto hole_it = replaced_holes.find(hole);
1388 if (hole_it == replaced_holes.end())
1389 continue;
1390
1391 ArbitraryHole & arb = *hole_it->second;
1392
1393 // We only need to update a replacement that's just had
1394 // new points inserted
1395 bool point_inserted = false;
1396 for (const Point & point : arb.get_points())
1397 if (next_boundary_node.count(point))
1398 {
1399 point_inserted = true;
1400 break;
1401 }
1402
1403 if (!point_inserted)
1404 continue;
1405
1406 // Find all points in the replacement hole
1407 std::vector<Point> new_points;
1408
1409 // Our outer polyline is expected to have points in
1410 // counter-clockwise order, so it proceeds "to the left"
1411 // from the point of view of rays inside the domain
1412 // pointing outward, and our next_boundary_node ordering
1413 // was filled accordingly.
1414 //
1415 // Our inner holes are expected to have points in
1416 // counter-clockwise order, but for holes "to the left
1417 // as viewed from the hole interior is the *opposite* of
1418 // "to the left as viewed from the domain interior". We
1419 // need to build the updated hole ordering "backwards".
1420
1421 // We should never see duplicate points when we add one
1422 // to a hole; if we do then we did something wrong.
1423 auto push_back_new_point = [&new_points](const Point & p) {
1424 // O(1) assert in devel
1425 libmesh_assert(new_points.empty() ||
1426 new_points.back() != p);
1427#ifdef DEBUG
1428 // O(N) asserts in dbg
1429 for (auto old_p : new_points)
1430 libmesh_assert_not_equal_to(old_p, p);
1431#endif
1432 new_points.push_back(p);
1433 };
1434
1435 for (auto point_it = arb.get_points().rbegin(),
1436 point_end = arb.get_points().rend();
1437 point_it != point_end;)
1438 {
1439 Point point = *point_it;
1440 ++point_it;
1441
1442 if (new_points.empty() ||
1443 (point != new_points.back() &&
1444 point != new_points.front()))
1445 push_back_new_point(point);
1446
1447 auto it = next_boundary_node.find(point);
1448 while (it != next_boundary_node.end())
1449 {
1450 point = *it->second;
1451 if (point == new_points.front())
1452 break;
1453 if (point_it != point_end &&
1454 point == *point_it)
1455 ++point_it;
1456 push_back_new_point(point);
1457 it = next_boundary_node.find(point);
1458 }
1459 }
1460
1461 std::reverse(new_points.begin(), new_points.end());
1462
1463 arb.set_points(std::move(new_points));
1464 }
1465 }
1466 }
1467
1468 // Okay, *now* we can add the new elements.
1469 for (auto & [raw_elem, unique_elem] : new_elems)
1470 {
1471 libmesh_assert_equal_to(raw_elem, unique_elem.get());
1472 libmesh_assert(!raw_elem->is_flipped());
1473 libmesh_ignore(raw_elem); // Old gcc warns "unused variable"
1474 mesh.add_elem(std::move(unique_elem));
1475 }
1476
1477 // Did we add anything?
1478 return !new_elems.empty();
1479}
The BoundaryInfo class contains information relevant to boundary conditions including storing faces,...
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
static const boundary_id_type invalid_id
Number used for internal use.
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure.
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
bool valid_id() const
Definition dof_object.h:861
static constexpr dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition dof_object.h:473
This is the base class from which all geometric element types are derived.
Definition elem.h:96
const Point & point(const unsigned int i) const
Definition elem.h:2462
static std::unique_ptr< Elem > build_with_id(const ElemType type, dof_id_type id)
Calls the build() method above with a nullptr parent, and additionally sets the newly-created Elem's ...
Definition elem.C:556
void set_neighbor(const unsigned int i, Elem *n)
Assigns n as the neighbor.
Definition elem.h:2635
unsigned int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition elem.h:2936
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
Definition elem.C:2784
const Node * node_ptr(const unsigned int i) const
Definition elem.h:2516
const Elem * neighbor_ptr(unsigned int i) const
Definition elem.h:2615
Point vertex_average() const
Definition elem.C:669
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition mesh_base.h:170
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
virtual dof_id_type max_elem_id() const =0
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
A Node is like a Point, but with more information.
Definition node.h:55
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
virtual FunctionBase< Real > * get_desired_area_function() override
Get the function giving desired triangle area as a function of position, or nullptr if no such functi...
bool is_refine_boundary_allowed(const BoundaryInfo &boundary_info, const Elem &elem, unsigned int side)
Is refining this element's boundary side allowed?
std::map< const Hole *, std::unique_ptr< ArbitraryHole > > replaced_holes
We might have to replace the user-provided holes with refined versions.
bool should_refine_elem(Elem &elem)
Returns true if the given element ought to be refined according to current criteria.
Real & minimum_angle()
Sets and/or gets the minimum desired angle.
bool has_auto_area_function()
Whether or not an auto area function has been set.
Real & desired_area()
Sets and/or gets the desired triangle area.
The UnstructuredMesh class is derived from the MeshBase class.
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true, const bool assert_valid=true) override
Other functions from MeshBase requiring re-definition.
static const Real b
int8_t boundary_id_type
Definition id_types.h:51
void libmesh_ignore(const Args &...)
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition libmesh.h:303
void libmesh_merge_move(T &target, T &source)
static constexpr Real TOLERANCE

References libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, _n_boundary_nodes, libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), b, libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build_with_id(), libMesh::Elem::contains_point(), libMesh::MeshBase::delete_elem(), libMesh::TriangulatorInterface::desired_area(), libMesh::UnstructuredMesh::find_neighbors(), libMesh::MeshBase::get_boundary_info(), get_desired_area_function(), libMesh::TriangulatorInterface::ArbitraryHole::get_points(), libMesh::TriangulatorInterface::has_auto_area_function(), libMesh::DofObject::id(), libMesh::DofObject::invalid_id, libMesh::BoundaryInfo::invalid_id, libMesh::invalid_uint, is_refine_boundary_allowed(), libMesh::Elem::level(), libMesh::libmesh_assert(), libMesh::libmesh_ignore(), libMesh::libmesh_merge_move(), libMesh::make_range(), libMesh::MeshBase::max_elem_id(), libMesh::MeshBase::max_node_id(), mesh, libMesh::TriangulatorInterface::minimum_angle(), libMesh::Elem::neighbor_ptr(), libMesh::MeshBase::node_ptr(), libMesh::Elem::node_ptr(), libMesh::Elem::point(), libMesh::Elem::quasicircumcenter(), libMesh::Real, libMesh::BoundaryInfo::remove(), replaced_holes, libMesh::TriangulatorInterface::segments, libMesh::Elem::set_neighbor(), libMesh::TriangulatorInterface::ArbitraryHole::set_points(), should_refine_elem(), libMesh::TOLERANCE, libMesh::TRI3, libMesh::Elem::type(), libMesh::DofObject::valid_id(), libMesh::Elem::vertex_average(), and libMesh::Elem::which_neighbor_am_i().

Referenced by triangulate().

◆ is_refine_boundary_allowed()

bool libMesh::Poly2TriTriangulator::is_refine_boundary_allowed ( const BoundaryInfo boundary_info,
const Elem elem,
unsigned int  side 
)
protected

Is refining this element's boundary side allowed?

Definition at line 466 of file poly2tri_triangulator.C.

470{
471 // We should only be calling this on a boundary side
472 libmesh_assert(!elem.neighbor_ptr(side));
473
474 std::vector<boundary_id_type> bcids;
475 boundary_info.boundary_ids(&elem, side, bcids);
476
477 // We should have one bcid on every boundary side.
478 libmesh_assert_equal_to(bcids.size(), 1);
479
480 if (bcids[0] == 0)
481 return this->refine_boundary_allowed();
482
483 // If we're not on an outer boundary side we'd better be on a hole
484 // side
485 libmesh_assert(this->_holes);
486
487 const boundary_id_type hole_num = bcids[0]-1;
488 libmesh_assert_less(hole_num, this->_holes->size());
489 const Hole * hole = (*this->_holes)[hole_num];
490 return hole->refine_boundary_allowed();
491}
virtual bool refine_boundary_allowed() const override
Get whether or not the triangulation is allowed to refine the mesh boundary when refining the interio...

References libMesh::TriangulatorInterface::_holes, libMesh::BoundaryInfo::boundary_ids(), libMesh::libmesh_assert(), libMesh::Elem::neighbor_ptr(), libMesh::TriangulatorInterface::Hole::refine_boundary_allowed(), and refine_boundary_allowed().

Referenced by insert_refinement_points().

◆ minimum_angle()

Real & libMesh::TriangulatorInterface::minimum_angle ( )
inlineinherited

Sets and/or gets the minimum desired angle.

Set to zero to disable angle constraint.

Definition at line 200 of file triangulator_interface.h.

200{return _minimum_angle;}
Real _minimum_angle
Minimum angle in triangles.

References libMesh::TriangulatorInterface::_minimum_angle.

Referenced by MeshTriangulationTest::commonSettings(), insert_refinement_points(), and main().

◆ nodes_to_segments()

void libMesh::TriangulatorInterface::nodes_to_segments ( dof_id_type  max_node_id)
protectedinherited

Helper function to create PSLG segments from our node ordering, up to the maximum node id, if no segments already exist.

Definition at line 248 of file triangulator_interface.C.

249{
250 // Don't try to override manually specified segments, or try to add
251 // segments if we're doing a convex hull
252 if (!this->segments.empty() || _triangulation_type != PSLG)
253 return;
254
255 for (auto node_it = _mesh.nodes_begin(),
256 node_end = _mesh.nodes_end();
257 node_it != node_end;)
258 {
259 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 if (node->id() >= max_node_id)
264 break;
265
266 ++node_it;
267
268 Node * next_node = (node_it == node_end) ?
269 *_mesh.nodes_begin() : *node_it;
270
271 this->segments.emplace_back(node->id(), next_node->id());
272 }
273
274 if (this->_verify_hole_boundaries && _holes)
275 {
276 std::vector<Point> outer_pts;
277 for (auto segment : this->segments)
278 outer_pts.push_back(_mesh.point(segment.first));
279
280 ArbitraryHole ah(outer_pts);
281 this->verify_holes(ah);
282 }
283}

References libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, libMesh::TriangulatorInterface::_triangulation_type, libMesh::TriangulatorInterface::_verify_hole_boundaries, libMesh::DofObject::id(), libMesh::MeshBase::point(), libMesh::TriangulatorInterface::PSLG, libMesh::TriangulatorInterface::segments, and libMesh::TriangulatorInterface::verify_holes().

Referenced by libMesh::TriangleInterface::triangulate(), and triangulate().

◆ quiet()

bool & libMesh::TriangulatorInterface::quiet ( )
inlineinherited

Whether not to silence internal messages to stdout.

Definition at line 249 of file triangulator_interface.h.

249{return _quiet;}
bool _quiet
Flag which tells if we want to suppress stdout outputs.

References libMesh::TriangulatorInterface::_quiet.

◆ refine_boundary_allowed()

virtual bool libMesh::Poly2TriTriangulator::refine_boundary_allowed ( ) const
inlineoverridevirtual

Get whether or not the triangulation is allowed to refine the mesh boundary when refining the interior.

True by default.

Reimplemented from libMesh::TriangulatorInterface.

Definition at line 104 of file poly2tri_triangulator.h.

105 { return _refine_bdy_allowed; }

References _refine_bdy_allowed.

Referenced by is_refine_boundary_allowed().

◆ set_auto_area_function()

void libMesh::TriangulatorInterface::set_auto_area_function ( const Parallel::Communicator comm,
const unsigned int  num_nearest_pts,
const unsigned int  power,
const Real  background_value,
const Real  background_eff_dist 
)
inherited

Generate an auto area function based on spacing of boundary points.

The external boundary as well as the hole boundaries are taken into consideration to generate the auto area function based on inverse distance interpolation. For each EDGE element on these boundaries, its centroid (midpoint) is used as the point position and the desired area is calculated as 1.5 times of the area of the equilateral triangle with the edge length as the length of the EDGE element. For a given position, the inverse distance interpolation only considers a number of nearest points (set by num_nearest_pts) to calculate the desired area. The weight of the value at each point is calculated as 1/distance^power.

In addition to these conventional inverse distance interpolation features, a concept of "background value" and "background effective distance" is introduced. The background value belongs to a virtual point located at a constant distance (background effective distance) from the given position. The weight of the value at this virtual point is calculated as 1/background_effective_distance^power. Effectively, the background value is the value when the given position is far away from the boundary points.

Definition at line 609 of file triangulator_interface.C.

614{
615 _auto_area_function = std::make_unique<AutoAreaFunction>(comm, num_nearest_pts, power, background_value, background_eff_dist);
616}

References libMesh::TriangulatorInterface::_auto_area_function.

◆ set_desired_area_function()

void libMesh::Poly2TriTriangulator::set_desired_area_function ( FunctionBase< Real > *  desired)
overridevirtual

Set a function giving desired triangle area as a function of position.

Set this to nullptr to disable position-dependent area constraint (falling back on desired_area()).

Reimplemented from libMesh::TriangulatorInterface.

Definition at line 450 of file poly2tri_triangulator.C.

452{
453 if (desired)
454 _desired_area_func = desired->clone();
455 else
456 _desired_area_func.reset();
457}
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0

References _desired_area_func, and libMesh::FunctionBase< Output >::clone().

Referenced by MeshTriangulationTest::testPoly2TriRefinementBase().

◆ set_fixup_tri7_center_nodes()

void libMesh::TriangulatorInterface::set_fixup_tri7_center_nodes ( bool  v)
inlineinherited

For TRI7 elements, after boundary mid-edge nodes have been snapped to a curved boundary in increase_triangle_order(), the interior node can be repositioned to the curved-mapping image of the reference centroid.

This is off by default; enable it only if the consuming code wants the relocated interior node.

Definition at line 273 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_fixup_tri7_center_nodes.

Referenced by MeshTriangulationTest::testPoly2TriEdge3ToTri7CenterFixup().

◆ set_interpolate_boundary_points()

void libMesh::TriangulatorInterface::set_interpolate_boundary_points ( int  n_points)
inherited

Complicated setter, for compatibility with insert_extra_points()

Definition at line 129 of file triangulator_interface.C.

130{
131 // Maybe we'll reserve a meaning for negatives later?
132 libmesh_assert(n_points >= 0);
133
135
136 // backwards compatibility - someone (including us) might want to
137 // query this via the old API.
138 _insert_extra_points = n_points;
139}

References libMesh::TriangulatorInterface::_insert_extra_points, libMesh::TriangulatorInterface::_interpolate_boundary_points, and libMesh::libmesh_assert().

Referenced by MeshTriangulationTest::testTriangulatorInterp().

◆ set_outer_boundary_ids()

void libMesh::TriangulatorInterface::set_outer_boundary_ids ( std::set< std::size_t >  bdy_ids)
inlineinherited

A set of ids to allow on the outer boundary loop: interpreted as boundary ids of 2D elements and/or subdomain ids of 1D edges.

If this is empty, then the outer boundary may be constructed from boundary edges of any id!

Definition at line 370 of file triangulator_interface.h.

370{ _bdy_ids = std::move(bdy_ids); }

References libMesh::TriangulatorInterface::_bdy_ids.

Referenced by MeshTriangulationTest::testHalfDomain().

◆ set_refine_boundary_allowed()

virtual void libMesh::Poly2TriTriangulator::set_refine_boundary_allowed ( bool  refine_bdy_allowed)
inlineoverridevirtual

Set whether or not the triangulation is allowed to refine the mesh boundary when refining the interior.

This is true by default, but may be set to false to make the mesh boundary more predictable (and so easier to stitch to other meshes) later.

Reimplemented from libMesh::TriangulatorInterface.

Definition at line 97 of file poly2tri_triangulator.h.

98 { _refine_bdy_allowed = refine_bdy_allowed; }

References _refine_bdy_allowed.

Referenced by MeshTriangulationTest::testPoly2TriHolesInterpRefined().

◆ set_verify_hole_boundaries()

void libMesh::TriangulatorInterface::set_verify_hole_boundaries ( bool  v)
inlineinherited

Verifying that hole boundaries don't cross the outer boundary or each other is something like O(N_bdys^2*N_points_per_bdy^2), so we only do it if requested.

Definition at line 262 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_verify_hole_boundaries.

Referenced by MeshTriangulationTest::commonSettings().

◆ should_refine_elem()

bool libMesh::Poly2TriTriangulator::should_refine_elem ( Elem elem)
protected

Returns true if the given element ought to be refined according to current criteria.

Definition at line 1482 of file poly2tri_triangulator.C.

1483{
1484 const Real min_area_target = this->desired_area();
1486
1487 // If this isn't a question, why are we here?
1488 libmesh_assert(min_area_target > 0 ||
1489 area_func != nullptr ||
1490 this->has_auto_area_function());
1491
1492 const Real area = elem.volume();
1493
1494 // If we don't have position-dependent area targets we can make a
1495 // decision quickly
1496 if (!area_func && !this->has_auto_area_function())
1497 return (area > min_area_target);
1498 else if(area_func && this->has_auto_area_function())
1499 libmesh_warning("WARNING: both desired are function and automatic area function are set. Using automatic area function.");
1500
1501 // If we do?
1502 //
1503 // See if we're meeting the local area target at all the elem
1504 // vertices first
1505 for (auto v : make_range(elem.n_vertices()))
1506 {
1507 // If we have an auto area function, we'll use it and override other area options
1508 const Real local_area_target = (*area_func)(elem.point(v));
1509 libmesh_error_msg_if
1510 (local_area_target <= 0,
1511 "Non-positive desired element areas are unachievable");
1512 if (area > local_area_target)
1513 return true;
1514 }
1515
1516 // If our vertices are happy, it's still possible that our interior
1517 // isn't. Are we allowed not to bother checking it?
1518 if (!min_area_target)
1519 return false;
1520
1521 libmesh_not_implemented_msg
1522 ("Combining a minimum desired_area with an area function isn't yet supported.");
1523}
virtual Real volume() const
Definition elem.C:3462
Base class for functors that can be evaluated at a point and (optionally) time.
FunctionBase< Real > * get_auto_area_function()
Get the auto area function.

References libMesh::TriangulatorInterface::desired_area(), libMesh::TriangulatorInterface::get_auto_area_function(), get_desired_area_function(), libMesh::TriangulatorInterface::has_auto_area_function(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::Elem::n_vertices(), libMesh::Elem::point(), libMesh::Real, and libMesh::Elem::volume().

Referenced by insert_refinement_points().

◆ smooth_after_generating()

bool & libMesh::TriangulatorInterface::smooth_after_generating ( )
inlineinherited

Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the grid.

Definition at line 244 of file triangulator_interface.h.

bool _smooth_after_generating
Flag which tells whether we should smooth the mesh after it is generated.

References libMesh::TriangulatorInterface::_smooth_after_generating.

Referenced by MeshTriangulationTest::commonSettings(), and triangulate_domain().

◆ total_hole_points()

unsigned int libMesh::TriangulatorInterface::total_hole_points ( )
protectedinherited

Helper function to count points in and verify holes.

Definition at line 585 of file triangulator_interface.C.

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 unsigned int n_hole_points = 0;
593
594 if (_holes)
595 for (const auto & hole : *_holes)
596 {
597 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 libmesh_assert_greater(hole->segment_indices().size(), 1);
603 libmesh_assert_equal_to(hole->segment_indices().back(), hole->n_points());
604 }
605
606 return n_hole_points;
607}

References libMesh::TriangulatorInterface::_holes.

Referenced by libMesh::TriangleInterface::triangulate().

◆ triangulate()

void libMesh::Poly2TriTriangulator::triangulate ( )
overridevirtual

Internally, this calls the poly2tri triangulation code in a loop, inserting our owner Steiner points as necessary to promote mesh quality.

Implements libMesh::TriangulatorInterface.

Definition at line 371 of file poly2tri_triangulator.C.

372{
373 LOG_SCOPE("triangulate()", "Poly2TriTriangulator");
374
375 // We only operate on serialized meshes. And it's not safe to
376 // serialize earlier, because it would then be possible for the user
377 // to re-parallelize the mesh in between there and here.
378 MeshSerializer serializer(_mesh);
379
380 // We don't yet support every set of Triangulator options in the
381 // poly2tri implementation
382
383 // We don't support convex hull triangulation, only triangulation of
384 // (implicitly defined, by node ordering) polygons (with holes if
385 // requested)
387 libmesh_not_implemented();
388
389 // We currently don't handle region specifications
390 if (_regions)
391 libmesh_not_implemented();
392
393 // We won't support quads any time soon, or 1D/3D in this interface
394 // ever.
395 if (_elem_type != TRI3 &&
396 _elem_type != TRI6 &&
397 _elem_type != TRI7)
398 libmesh_not_implemented();
399
400 // If we have no explicit segments defined, we may get them from
401 // mesh elements
402 this->elems_to_segments();
403
404 // If we *still* have no explicit segments defined, we get them from
405 // the order of nodes.
407
408 // Insert additional new points in between existing boundary points,
409 // if that is requested and reasonable
411
412 // Triangulate the points we have, then see if we need to add more;
413 // repeat until we don't need to add more.
414 //
415 // This is currently done redundantly in parallel; make sure no
416 // processor quits before the others.
417 do
418 {
419 libmesh_parallel_only(_mesh.comm());
421 }
422 while (this->insert_refinement_points());
423
424 libmesh_parallel_only(_mesh.comm());
425
426 // Okay, we really do need to support boundary ids soon, but we
427 // don't yet
428 if (_markers)
429 libmesh_not_implemented();
430
432
433 // To the naked eye, a few smoothing iterations usually looks better,
434 // so we do this by default unless the user says not to.
435 if (this->_smooth_after_generating)
437
438 // The user might have requested TRI6 or higher instead of TRI3. We
439 // can do this before prepare_for_use() because all we need for it
440 // is find_neighbors(), which we did in insert_refinement_points()
442
443 // Prepare the mesh for use before returning. This ensures (among
444 // other things) that it is partitioned and therefore users can
445 // iterate over local elements, etc.
447}
This class defines the data structures necessary for Laplace smoothing.
virtual void smooth() override
Redefinition of the smooth function from the base class.
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly created (or read) mesh for use.
Definition mesh_base.C:824
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition mesh_base.h:423
Temporarily serialize a DistributedMesh for non-distributed-mesh capable code paths.
const Parallel::Communicator & comm() const
bool insert_refinement_points()
Add Steiner points as new mesh nodes, as necessary to refine an existing trangulation.
void triangulate_current_points()
Triangulate the current mesh and hole points.
void increase_triangle_order()
Helper function to upconvert Tri3 to any higher order triangle type if requested via _elem_type.
void nodes_to_segments(dof_id_type max_node_id)
Helper function to create PSLG segments from our node ordering, up to the maximum node id,...
void elems_to_segments()
Helper function to create PSLG segments from our other boundary-defining options (1D mesh edges,...
void insert_any_extra_boundary_points()
Helper function to add extra points (midpoints of initial segments) to a PSLG triangulation.

References libMesh::TriangulatorInterface::_elem_type, libMesh::TriangulatorInterface::_markers, libMesh::TriangulatorInterface::_mesh, _n_boundary_nodes, libMesh::TriangulatorInterface::_regions, libMesh::TriangulatorInterface::_smooth_after_generating, libMesh::TriangulatorInterface::_triangulation_type, libMesh::ParallelObject::comm(), libMesh::TriangulatorInterface::elems_to_segments(), libMesh::TriangulatorInterface::increase_triangle_order(), libMesh::TriangulatorInterface::insert_any_extra_boundary_points(), insert_refinement_points(), libMesh::TriangulatorInterface::nodes_to_segments(), libMesh::MeshBase::prepare_for_use(), libMesh::TriangulatorInterface::PSLG, libMesh::MeshBase::set_mesh_dimension(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, and triangulate_current_points().

Referenced by main(), MeshTriangulationTest::testPoly2TriEdge3ToTri6FlipThrows(), MeshTriangulationTest::testPoly2TriEdge3ToTri7CenterFixup(), and MeshTriangulationTest::testPoly2TriRefinementBase().

◆ triangulate_current_points()

void libMesh::Poly2TriTriangulator::triangulate_current_points ( )
protected

Triangulate the current mesh and hole points.

Definition at line 494 of file poly2tri_triangulator.C.

495{
496 LOG_SCOPE("triangulate_current_points()", "Poly2TriTriangulator");
497
498 // Will the triangulation have holes?
499 const std::size_t n_holes = _holes != nullptr ? _holes->size() : 0;
500
501 // Mapping from Poly2Tri points to libMesh nodes, so we can get the
502 // connectivity translated back later.
503 std::map<const p2t::Point, Node *, P2TPointCompare> point_node_map;
504
505 // Poly2Tri data structures
506 // Poly2Tri takes vectors of pointers-to-Point for some reason, but
507 // we'll just make those shims to vectors of Point rather than
508 // individually/manually heap allocating everything.
509 std::vector<p2t::Point> outer_boundary_points;
510 std::vector<std::vector<p2t::Point>> inner_hole_points(n_holes);
511
513 libmesh_error_msg_if
514 (!nn, "Poly2TriTriangulator cannot triangulate an empty mesh!");
515
516 // Unless we're using an explicit segments list, we assume node ids
517 // are contiguous here.
518 if (this->segments.empty())
519 libmesh_error_msg_if
520 (_mesh.n_nodes() != nn,
521 "Poly2TriTriangulator needs contiguous node ids or explicit segments!");
522
523 // And if we have more nodes than outer boundary points, the rest
524 // may be interior "Steiner points". We use a set here so we can
525 // cheaply search and erase from it later, when we're identifying
526 // hole points.
527 std::set<p2t::Point, P2TPointCompare> steiner_points;
528
529 // If we were asked to use all mesh nodes as boundary nodes, now's
530 // the time to see how many that is.
532 {
534 libmesh_assert_equal_to(std::ptrdiff_t(_n_boundary_nodes),
535 std::distance(_mesh.nodes_begin(),
536 _mesh.nodes_end()));
537
538 }
539 else
540 libmesh_assert_less_equal(_n_boundary_nodes,
541 _mesh.n_nodes());
542
543 // Prepare poly2tri points for our nodes, sorted into outer boundary
544 // points and interior Steiner points.
545
546 if (this->segments.empty())
547 {
548 // If we have no segments even after taking elems into account,
549 // the nodal id ordering defines our outer polyline ordering
550 for (auto & node : _mesh.node_ptr_range())
551 {
552 const p2t::Point pt = to_p2t(*node);
553
554 // If we're out of boundary nodes, the rest are going to be
555 // Steiner points or hole points
556 if (node->id() < _n_boundary_nodes)
557 outer_boundary_points.push_back(pt);
558 else
559 steiner_points.insert(pt);
560
561 // We're not going to support overlapping nodes on the boundary
562 if (point_node_map.count(pt))
563 libmesh_not_implemented();
564
565 point_node_map.emplace(pt, node);
566 }
567 }
568 // If we have explicit segments defined, that's our outer polyline
569 // ordering:
570 else
571 {
572 // Let's make sure our segments are in order
574
575 // Add nodes from every segment, in order, to the outer polyline
576 for (auto [segment_start, segment_end] : this->segments)
577 {
578 if (last_id != DofObject::invalid_id)
579 libmesh_error_msg_if(segment_start != last_id,
580 "Disconnected triangulator segments");
581 last_id = segment_end;
582
583 Node * node = _mesh.query_node_ptr(segment_start);
584
585 libmesh_error_msg_if(!node,
586 "Triangulator segments reference nonexistent node id " <<
587 segment_start);
588
589 outer_boundary_points.emplace_back(double((*node)(0)), double((*node)(1)));
590 p2t::Point * pt = &outer_boundary_points.back();
591
592 // We're not going to support overlapping nodes on the boundary
593 if (point_node_map.count(*pt))
594 libmesh_not_implemented_msg
595 ("Triangulating overlapping boundary nodes is unsupported");
596
597 point_node_map.emplace(*pt, node);
598 }
599
600 libmesh_error_msg_if(last_id != this->segments[0].first,
601 "Non-closed-loop triangulator segments");
602
603 // If we have points that aren't in any segments, those will be
604 // Steiner points
605 for (auto & node : _mesh.node_ptr_range())
606 {
607 const p2t::Point pt = to_p2t(*node);
608 if (const auto it = point_node_map.find(pt);
609 it == point_node_map.end())
610 {
611 steiner_points.insert(pt);
612 point_node_map.emplace(pt, node);
613 }
614 else
615 libmesh_assert_equal_to(it->second, node);
616 }
617 }
618
619 // If we have any elements from a previous triangulation, we're
620 // going to replace them with our own. If we have any elements that
621 // were used to create our segments, we're done creating and we no
622 // longer need them.
624
625 // Keep track of what boundary ids we want to assign to each new
626 // triangle. We'll give the outer boundary BC 0, and give holes ids
627 // starting from 1. We've already got the point_node_map to find
628 // nodes, so we can just key on pairs of node ids to identify a side.
629 std::unordered_map<std::pair<dof_id_type,dof_id_type>,
630 boundary_id_type, libMesh::hash> side_boundary_id;
631
632 const boundary_id_type outer_bcid = 0;
633 const std::size_t n_outer = outer_boundary_points.size();
634
635 for (auto i : make_range(n_outer))
636 {
637 const Node * node1 =
638 libmesh_map_find(point_node_map, outer_boundary_points[i]),
639 * node2 =
640 libmesh_map_find(point_node_map, outer_boundary_points[(i+1)%n_outer]);
641
642 side_boundary_id.emplace(std::make_pair(node1->id(),
643 node2->id()),
644 outer_bcid);
645 }
646
647 // Create poly2tri triangulator with our mesh points
648 std::vector<p2t::Point *> outer_boundary_pointers(n_outer);
649 std::transform(outer_boundary_points.begin(),
650 outer_boundary_points.end(),
651 outer_boundary_pointers.begin(),
652 [](p2t::Point & p) { return &p; });
653
654
655 // Make sure shims for holes last as long as the CDT does; the
656 // poly2tri headers don't make clear whether or not they're hanging
657 // on to references to these vectors, and it would be reasonable for
658 // them to do so.
659 std::vector<std::vector<p2t::Point *>> inner_hole_pointers(n_holes);
660
661 p2t::CDT cdt{outer_boundary_pointers};
662
663 // Add any holes
664 for (auto h : make_range(n_holes))
665 {
666 const Hole * initial_hole = (*_holes)[h];
667 auto it = replaced_holes.find(initial_hole);
668 const Hole & our_hole =
669 (it == replaced_holes.end()) ?
670 *initial_hole : *it->second;
671 auto & poly2tri_hole = inner_hole_points[h];
672
673 for (auto i : make_range(our_hole.n_points()))
674 {
675 Point p = our_hole.point(i);
676 poly2tri_hole.emplace_back(to_p2t(p));
677
678 const auto & pt = poly2tri_hole.back();
679
680 // This won't be a steiner point.
681 steiner_points.erase(pt);
682
683 // If we see a hole point already in the mesh, we'll share
684 // that node. This might be a problem if it's a boundary
685 // node, but it might just be the same hole point already
686 // added during a previous triangulation refinement step.
687 if (point_node_map.count(pt))
688 {
689 libmesh_assert_equal_to
690 (point_node_map[pt],
691 _mesh.query_node_ptr(point_node_map[pt]->id()));
692 }
693 else
694 {
695 Node * node = _mesh.add_point(p, nn++);
696 point_node_map[pt] = node;
697 }
698 }
699
700 const boundary_id_type inner_bcid = h+1;
701 const std::size_t n_inner = poly2tri_hole.size();
702
703 for (auto i : make_range(n_inner))
704 {
705 const Node * node1 =
706 libmesh_map_find(point_node_map, poly2tri_hole[i]),
707 * node2 =
708 libmesh_map_find(point_node_map, poly2tri_hole[(i+1)%n_inner]);
709
710 side_boundary_id.emplace(std::make_pair(node1->id(),
711 node2->id()),
712 inner_bcid);
713 }
714
715 auto & poly2tri_ptrs = inner_hole_pointers[h];
716 poly2tri_ptrs.resize(n_inner);
717
718 std::transform(poly2tri_hole.begin(),
719 poly2tri_hole.end(),
720 poly2tri_ptrs.begin(),
721 [](p2t::Point & p) { return &p; });
722
723 cdt.AddHole(poly2tri_ptrs);
724 }
725
726 // Add any steiner points. We had them in a set, but post-C++11
727 // that won't give us non-const element access (even if we
728 // pinky-promise not to change the elements in any way that affects
729 // our Comparator), and Poly2Tri wants non-const elements (to store
730 // edge data?), so we need to move them here.
731 std::vector<p2t::Point> steiner_vector(steiner_points.begin(), steiner_points.end());
732 steiner_points.clear();
733 for (auto & p : steiner_vector)
734 cdt.AddPoint(&p);
735
736 // Triangulate!
737 cdt.Triangulate();
738
739 // Get poly2tri triangles, turn them into libMesh triangles
740 std::vector<p2t::Triangle *> triangles = cdt.GetTriangles();
741
742 // Do our own numbering, even on DistributedMesh
743 dof_id_type next_id = 0;
744
745 BoundaryInfo & boundary_info = _mesh.get_boundary_info();
746 boundary_info.clear();
747
748 // Add the triangles to our Mesh data structure.
749 for (auto ptri_ptr : triangles)
750 {
751 p2t::Triangle & ptri = *ptri_ptr;
752
753 // We always use TRI3 here, since that's what we have nodes for;
754 // if we need a higher order we can convert at the end.
755 auto elem = Elem::build_with_id(TRI3, next_id++);
756 for (auto v : make_range(3))
757 {
758 const p2t::Point & vertex = *ptri.GetPoint(v);
759
760 Node * node = libmesh_map_find(point_node_map, vertex);
761 libmesh_assert(node);
762 elem->set_node(v, node);
763 }
764
765 // We expect a consistent triangle orientation
766 libmesh_assert(!elem->is_flipped());
767
768 Elem * added_elem = _mesh.add_elem(std::move(elem));
769
770 for (auto v : make_range(3))
771 {
772 const Node & node1 = added_elem->node_ref(v),
773 & node2 = added_elem->node_ref((v+1)%3);
774
775 auto it = side_boundary_id.find(std::make_pair(node1.id(), node2.id()));
776 if (it == side_boundary_id.end())
777 it = side_boundary_id.find(std::make_pair(node2.id(), node1.id()));
778 if (it != side_boundary_id.end())
779 boundary_info.add_side(added_elem, v, it->second);
780 }
781 }
782}
void clear()
Clears the underlying data structures and restores the object to a pristine state with no data stored...
const Node & node_ref(const unsigned int i) const
Definition elem.h:2538
virtual dof_id_type n_nodes() const =0
virtual const Node * query_node_ptr(const dof_id_type i) const =0

References libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, _n_boundary_nodes, libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::Elem::build_with_id(), libMesh::BoundaryInfo::clear(), libMesh::MeshBase::clear_elems(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::DofObject::invalid_id, libMesh::Elem::is_flipped(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::max_node_id(), libMesh::MeshBase::n_nodes(), libMesh::TriangulatorInterface::Hole::n_points(), libMesh::Elem::node_ref(), libMesh::TriangulatorInterface::Hole::point(), libMesh::MeshBase::query_node_ptr(), replaced_holes, libMesh::TriangulatorInterface::segments, libMesh::Elem::set_node(), and libMesh::TRI3.

Referenced by triangulate().

◆ triangulation_type()

TriangulationType & libMesh::TriangulatorInterface::triangulation_type ( )
inlineinherited

◆ verify_holes()

void libMesh::TriangulatorInterface::verify_holes ( const Hole outer_bdy)
protectedinherited

Helper function to check holes for intersections if requested.

Definition at line 560 of file triangulator_interface.C.

561{
562 for (const Hole * hole : *_holes)
563 {
564 for (const Hole * hole2 : *_holes)
565 {
566 if (hole == hole2)
567 continue;
568
569 for (auto i : make_range(hole2->n_points()))
570 if (hole->contains(hole2->point(i)))
571 libmesh_error_msg
572 ("Found point " << hole2->point(i) <<
573 " on one hole boundary and another's interior");
574 }
575
576 for (auto i : make_range(hole->n_points()))
577 if (!outer_bdy.contains(hole->point(i)))
578 libmesh_error_msg
579 ("Found point " << hole->point(i) <<
580 " on hole boundary but outside outer boundary");
581 }
582}

References libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::Hole::contains(), and libMesh::make_range().

Referenced by libMesh::TriangulatorInterface::elems_to_segments(), and libMesh::TriangulatorInterface::nodes_to_segments().

◆ verify_quadratic_elements()

void libMesh::TriangulatorInterface::verify_quadratic_elements ( )
protectedinherited

Helper called at the end of increase_triangle_order() to verify that no boundary-midpoint snap has tangled a quadratic triangle, by sampling the signed element Jacobian at the reference nodes and centroid.

Errors out with a diagnostic naming the bad element and side if the Jacobian is non-positive.

Definition at line 473 of file triangulator_interface.C.

474{
475 if (_elem_type != TRI6 && _elem_type != TRI7)
476 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 for (Elem * elem : _mesh.element_ptr_range())
488 {
489 libmesh_assert_equal_to(elem->n_vertices(), 3);
490 libmesh_assert_greater_equal(elem->n_nodes(), 6u);
491
492 const Point & x0 = elem->point(0);
493 const Point & x1 = elem->point(1);
494 const Point & x2 = elem->point(2);
495 const Point & x3 = elem->point(3);
496 const Point & x4 = elem->point(4);
497 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 const Point a1 = 4*x0 + 4*x1 - 8*x3;
502 const Point b1 = 4*x0 - 4*x3 + 4*x4 - 4*x5;
503 const Point c1 = -3*x0 - 1*x1 + 4*x3;
504 const Point b2 = 4*x0 + 4*x2 - 8*x5;
505 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 const Real ref_area = 0.5 * cross_norm(x1 - x0, x2 - x0);
510 const Real jac_tol = TOLERANCE * ref_area;
511
512 Real min_jac = std::numeric_limits<Real>::max();
513 unsigned int worst_sample = 0;
514 for (unsigned int s = 0; s != 7; ++s)
515 {
516 const Real xi = xi_samples[s];
517 const Real eta = eta_samples[s];
518 const Point dxi = xi*a1 + eta*b1 + c1;
519 const Point deta = xi*b1 + eta*b2 + c2;
520 // z-component of the cross product; the elements are planar.
521 const Real jac = dxi(0)*deta(1) - dxi(1)*deta(0);
522 if (jac < min_jac)
523 {
524 min_jac = jac;
525 worst_sample = s;
526 }
527 }
528
529 if (min_jac > jac_tol)
530 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 std::ostringstream sides;
536 for (unsigned int n = 0; n != 3; ++n)
537 if (!elem->neighbor_ptr(n))
538 {
539 const Point straight =
540 0.5 * (elem->point(n) + elem->point((n+1) % 3));
541 sides << " (boundary side " << n
542 << ": straight midpoint " << straight
543 << ", snapped midpoint " << elem->point(n+3) << ")";
544 }
545
546 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 }
557}
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Calls cross_norm_sq() and takes the square root of the result.

References libMesh::TriangulatorInterface::_elem_type, libMesh::TriangulatorInterface::_mesh, libMesh::cross_norm(), libMesh::Real, libMesh::TOLERANCE, libMesh::TRI6, and libMesh::TRI7.

Referenced by libMesh::TriangulatorInterface::increase_triangle_order().

Member Data Documentation

◆ _auto_area_function

std::unique_ptr<AutoAreaFunction> libMesh::TriangulatorInterface::_auto_area_function
protectedinherited

◆ _bdy_ids

std::set<std::size_t> libMesh::TriangulatorInterface::_bdy_ids
protectedinherited

◆ _desired_area

Real libMesh::TriangulatorInterface::_desired_area
protectedinherited

The desired area for the elements in the resulting mesh.

Definition at line 463 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::desired_area(), and libMesh::TriangleInterface::triangulate().

◆ _desired_area_func

std::unique_ptr<FunctionBase<Real> > libMesh::Poly2TriTriangulator::_desired_area_func
private

Location-dependent area requirements.

Definition at line 149 of file poly2tri_triangulator.h.

Referenced by get_desired_area_function(), and set_desired_area_function().

◆ _elem_type

ElemType libMesh::TriangulatorInterface::_elem_type
protectedinherited

◆ _fixup_tri7_center_nodes

bool libMesh::TriangulatorInterface::_fixup_tri7_center_nodes
protectedinherited

Flag which tells if we want to relocate TRI7 interior nodes to the curved-mapping centroid after snapping boundary midpoints.

Definition at line 513 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::get_fixup_tri7_center_nodes(), libMesh::TriangulatorInterface::increase_triangle_order(), and libMesh::TriangulatorInterface::set_fixup_tri7_center_nodes().

◆ _holes

const std::vector<Hole*>* libMesh::TriangulatorInterface::_holes
protectedinherited

◆ _insert_extra_points

bool libMesh::TriangulatorInterface::_insert_extra_points
protectedinherited

Flag which tells whether or not to insert additional nodes before triangulation.

This can sometimes be used to "de-regularize" the resulting triangulation.

This flag is supported for backwards compatibility; setting _interpolate_boundary_points = 1 is equivalent.

Definition at line 485 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::get_interpolate_boundary_points(), libMesh::TriangulatorInterface::insert_extra_points(), and libMesh::TriangulatorInterface::set_interpolate_boundary_points().

◆ _interpolate_boundary_points

int libMesh::TriangulatorInterface::_interpolate_boundary_points
protectedinherited

Flag which tells how many additional nodes should be inserted between each pair of original mesh points.

Definition at line 491 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::get_interpolate_boundary_points(), and libMesh::TriangulatorInterface::set_interpolate_boundary_points().

◆ _markers

const std::vector<int>* libMesh::TriangulatorInterface::_markers
protectedinherited

◆ _mesh

UnstructuredMesh& libMesh::TriangulatorInterface::_mesh
protectedinherited

◆ _minimum_angle

Real libMesh::TriangulatorInterface::_minimum_angle
protectedinherited

Minimum angle in triangles.

Definition at line 468 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::minimum_angle(), and libMesh::TriangleInterface::triangulate().

◆ _n_boundary_nodes

dof_id_type libMesh::Poly2TriTriangulator::_n_boundary_nodes
private

Keep track of how many mesh nodes are boundary nodes.

Definition at line 144 of file poly2tri_triangulator.h.

Referenced by insert_refinement_points(), triangulate(), and triangulate_current_points().

◆ _quiet

bool libMesh::TriangulatorInterface::_quiet
protectedinherited

Flag which tells if we want to suppress stdout outputs.

Definition at line 502 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::quiet(), and libMesh::TriangleInterface::triangulate().

◆ _refine_bdy_allowed

bool libMesh::Poly2TriTriangulator::_refine_bdy_allowed
private

Whether to allow boundary refinement.

Definition at line 154 of file poly2tri_triangulator.h.

Referenced by refine_boundary_allowed(), and set_refine_boundary_allowed().

◆ _regions

const std::vector<Region*>* libMesh::TriangulatorInterface::_regions
protectedinherited

A pointer to a vector of Regions*s.

If this is nullptr, there are no regions!

Definition at line 447 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::attach_region_list(), libMesh::TriangleInterface::triangulate(), and triangulate().

◆ _smooth_after_generating

bool libMesh::TriangulatorInterface::_smooth_after_generating
protectedinherited

Flag which tells whether we should smooth the mesh after it is generated.

True by default.

Definition at line 497 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::smooth_after_generating(), libMesh::TriangleInterface::triangulate(), and triangulate().

◆ _triangulation_type

TriangulationType libMesh::TriangulatorInterface::_triangulation_type
protectedinherited

◆ _verify_hole_boundaries

bool libMesh::TriangulatorInterface::_verify_hole_boundaries
protectedinherited

◆ replaced_holes

std::map<const Hole *, std::unique_ptr<ArbitraryHole> > libMesh::Poly2TriTriangulator::replaced_holes
private

We might have to replace the user-provided holes with refined versions.

Definition at line 139 of file poly2tri_triangulator.h.

Referenced by insert_refinement_points(), and triangulate_current_points().

◆ segment_midpoints

std::vector<Point> libMesh::TriangulatorInterface::segment_midpoints
inherited

When constructing a second-order triangulation from a second-order boundary, we may do the triangulation using first-order elements, in which case we need to save midpoint location data in order to reconstruct curvature along boundaries.

Definition at line 298 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::elems_to_segments(), and libMesh::TriangulatorInterface::increase_triangle_order().

◆ segment_midpoints_keys

std::vector<Point> libMesh::TriangulatorInterface::segment_midpoints_keys
inherited

When saving the midpoint location data, we need to save the corresponding segment information too.

Here the first point of the segment is saved so that it can be used as a key to find the corresponding segment midpoint.

Definition at line 306 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::elems_to_segments(), and libMesh::TriangulatorInterface::increase_triangle_order().

◆ segments

std::vector<std::pair<unsigned int, unsigned int> > libMesh::TriangulatorInterface::segments
inherited

When constructing a PSLG, if the node numbers do not define the desired boundary segments implicitly through the ordering of the points, you can use the segments vector to specify the segments explicitly, Ex: unit square numbered counter-clockwise starting from origin segments[0] = (0,1) segments[1] = (1,2) segments[2] = (2,3) segments[3] = (3,0) (For the above case you could actually use the implicit ordering!)

Definition at line 290 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::elems_to_segments(), libMesh::TriangulatorInterface::increase_triangle_order(), libMesh::TriangulatorInterface::insert_any_extra_boundary_points(), insert_refinement_points(), libMesh::TriangulatorInterface::nodes_to_segments(), MeshTriangulationTest::testTriangulatorSegments(), libMesh::TriangleInterface::triangulate(), and triangulate_current_points().


The documentation for this class was generated from the following files: