20 #include "libmesh/simplex_refiner.h" 22 #include "libmesh/boundary_info.h" 23 #include "libmesh/face_tri3.h" 24 #include "libmesh/int_range.h" 25 #include "libmesh/libmesh_logging.h" 26 #include "libmesh/mesh_base.h" 27 #include "libmesh/mesh_communication.h" 28 #include "libmesh/mesh_tools.h" 29 #include "libmesh/node.h" 30 #include "libmesh/parallel_algebra.h" 31 #include "libmesh/remote_elem.h" 32 #include "libmesh/unstructured_mesh.h" 33 #include "libmesh/utility.h" 72 volume_func !=
nullptr);
80 (
"Warning: trying to refine sliver element with volume " <<
84 return (
volume > min_volume_target);
93 const Real local_volume_target = (*volume_func)(elem.
point(v));
95 (local_volume_target <= 0,
96 "Non-positive desired element volumes are unachievable");
97 if (
volume > local_volume_target)
103 if (!min_volume_target)
106 libmesh_not_implemented_msg
107 (
"Combining a minimum desired_volume with an volume function isn't yet supported.");
120 libmesh_assert_equal_to(elem.
level(), 0u);
124 libmesh_assert_equal_to(elem.
type(),
TRI3);
126 const int n_sides = elem.
n_sides();
127 int side_to_refine = 0;
128 Real length_to_refine = 0;
129 std::pair<Node *, Node *> vertices_to_refine;
135 std::pair<Node *, Node *> vertices
139 if ((
Point&)(*vertices.first) >
140 (
Point&)(*vertices.second))
141 std::swap(vertices.first, vertices.second);
143 const auto sidevec = *(
Point*)vertices.first -
144 *(
Point*)vertices.second;
145 Real length = sidevec.norm();
153 elem.
id() == coarse_id)
155 vertices.second->id());
160 if (length > length_to_refine &&
164 side_to_refine = side;
165 length_to_refine = length;
166 vertices_to_refine = vertices;
171 if (!length_to_refine)
178 std::vector<Elem *> neighbors;
180 neighbors.push_back(neigh);
183 if (
auto it =
new_nodes.find(vertices_to_refine);
186 midedge_node = it->second;
194 *(
Point*)vertices_to_refine.second)/2,
200 for (
const Elem * neigh : neighbors)
203 neigh->processor_id());
205 new_nodes[vertices_to_refine] = midedge_node;
209 constexpr
int max_subelems = 2;
210 std::unique_ptr<Tri3> subelem[max_subelems];
213 subelem[0] = std::make_unique<Tri3>();
214 subelem[0]->set_node(0, elem.
node_ptr(side_to_refine));
215 subelem[0]->set_node(1, midedge_node);
216 subelem[0]->set_node(2, elem.
node_ptr((side_to_refine+2)%3));
218 subelem[1] = std::make_unique<Tri3>();
219 subelem[1]->set_node(0, elem.
node_ptr((side_to_refine+2)%3));
220 subelem[1]->set_node(1, midedge_node);
221 subelem[1]->set_node(2, elem.
node_ptr((side_to_refine+1)%3));
232 std::vector<boundary_id_type> bc_ids;
233 boundary_info.
boundary_ids(&elem, side_to_refine, bc_ids);
234 boundary_info.
add_side(subelem[0].
get(), 0, bc_ids);
235 boundary_info.
add_side(subelem[1].
get(), 1, bc_ids);
236 subelem[0]->set_neighbor(0, elem.
neighbor_ptr(side_to_refine));
237 subelem[1]->set_neighbor(1, elem.
neighbor_ptr(side_to_refine));
239 boundary_info.
boundary_ids(&elem, (side_to_refine+1)%3, bc_ids);
240 boundary_info.
add_side(subelem[1].
get(), 2, bc_ids);
241 subelem[1]->set_neighbor(2, elem.
neighbor_ptr((side_to_refine+1)%3));
243 boundary_info.
boundary_ids(&elem, (side_to_refine+2)%3, bc_ids);
244 boundary_info.
add_side(subelem[0].
get(), 2, bc_ids);
245 subelem[0]->set_neighbor(2, elem.
neighbor_ptr((side_to_refine+2)%3));
249 for (
unsigned int i=0; i != max_subelems; ++i)
260 subelem[i]->add_extra_integers(nei);
261 for (
unsigned int ei=0; ei != nei; ++ei)
264 subelem[i]->inherit_data_from(elem);
267 boundary_info.
remove(&elem);
272 auto & added = it->second;
274 sub_it != added.end())
276 if (&elem == sub_it->second)
294 std::size_t n_refined_elements = 1;
296 for (
unsigned int i=0; i != max_subelems; ++i)
298 Elem * add_elem = subelem[i].get();
303 n_refined_elements +=
307 return n_refined_elements;
320 for (
auto & elem :
_mesh.element_ptr_range())
321 for (
auto n :
make_range(elem->n_neighbors()))
329 if (!proxy_neighbor.get())
332 elem->set_neighbor(n, proxy_neighbor.get());
338 auto edge_gather_functor =
341 const std::vector<std::pair<dof_id_type, dof_id_type>> & edges,
342 std::vector<refinement_datum> & edge_refinements)
345 const std::size_t query_size = edges.size();
346 edge_refinements.resize(query_size);
347 for (std::size_t i=0; i != query_size; ++i)
349 auto vertex_ids = edges[i];
350 std::pair<Node *, Node *> vertices
357 auto edge_action_functor =
360 const std::vector<std::pair<dof_id_type, dof_id_type>> &,
361 const std::vector<refinement_datum> & edge_refinements
364 for (
const auto & one_edges_refinements : edge_refinements)
365 for (
const auto & refinement : one_edges_refinements)
367 std::pair<Node *, Node *> vertices
371 if ((
Point&)(*vertices.first) >
372 (
Point&)(*vertices.second))
373 std::swap(vertices.first, vertices.second);
379 std::get<2>(refinement));
380 it->second->processor_id() = std::get<3>(refinement);
386 *(
Point*)vertices.second)/2,
387 std::get<2>(refinement),
388 std::get<3>(refinement));
394 std::size_t refined_elements = 0;
395 std::size_t newly_refined_elements = 0;
398 newly_refined_elements = 0;
402 for (
auto & elem :
_mesh.element_ptr_range())
407 elem->id() : it->second;
409 newly_refined_elements +=
412 refined_elements += newly_refined_elements;
415 this->new_elements.clear();
425 if (have_edge_queries)
426 Parallel::pull_parallel_vector_data
428 edge_action_functor, refinement_data_ex);
431 while (newly_refined_elements);
439 std::unordered_map<processor_id_type, std::vector<dof_id_type>> elems_to_query;
442 if (added_elem_map.empty())
446 added_elem_map.begin()->second->processor_id();
450 elems_to_query[pid].push_back(coarse_id);
455 #ifdef LIBMESH_ENABLE_UNIQUE_ID 456 std::vector<std::tuple<Point, dof_id_type, unique_id_type>>
458 std::vector<std::tuple<Point, dof_id_type>>
460 elem_refinement_datum;
462 auto added_gather_functor =
465 const std::vector<dof_id_type> & coarse_elems,
466 std::vector<elem_refinement_datum> & coarse_refinements)
469 const std::size_t query_size = coarse_elems.size();
470 coarse_refinements.resize(query_size);
477 for (
auto [vertex_avg, elem] : added)
479 coarse_refinements[i].emplace_back
482 #ifdef LIBMESH_ENABLE_UNIQUE_ID 490 auto added_action_functor =
493 const std::vector<dof_id_type> & coarse_elems,
494 const std::vector<elem_refinement_datum> & coarse_refinements)
496 const std::size_t query_size = coarse_elems.size();
500 const auto & refinement_data = coarse_refinements[i];
501 const auto & our_added =
503 for (
auto [vertex_avg,
id 504 #ifdef LIBMESH_ENABLE_UNIQUE_ID
509 Elem & our_elem = *libmesh_map_find(our_added,
514 #ifdef LIBMESH_ENABLE_UNIQUE_ID 522 elem_refinement_datum * refinement_data_ex =
nullptr;
523 Parallel::pull_parallel_vector_data
524 (
_mesh.
comm(), elems_to_query, added_gather_functor,
525 added_action_functor, refinement_data_ex);
530 mc.make_node_ids_parallel_consistent (
_mesh);
531 mc.make_node_unique_ids_parallel_consistent (
_mesh);
540 # ifdef LIBMESH_ENABLE_UNIQUE_ID 545 return refined_elements;
554 _desired_volume_func = desired->
clone();
556 _desired_volume_func.reset();
574 vec.emplace_back(vertices.first->id(),
575 vertices.second->id(),
577 it->second->processor_id());
579 std::pair<Node *, Node *> subedge {vertices.first, it->second};
580 if ((
Point&)(*subedge.first) >
581 (
Point&)(*subedge.second))
582 std::swap(subedge.first, subedge.second);
585 subedge = std::make_pair(it->second, vertices.second);
586 if ((
Point&)(*subedge.first) >
587 (
Point&)(*subedge.second))
588 std::swap(subedge.first, subedge.second);
std::unordered_map< dof_id_type, std::unordered_map< Point, Elem * > > added_elements
Keep track of elements added within each original element, so we can answer queries about them from o...
virtual void set_desired_volume_function(FunctionBase< Real > *desired)
Set a function giving desired element volume as a function of position.
unique_id_type & set_unique_id()
std::unique_ptr< FunctionBase< Real > > _desired_volume_func
Location-dependent volume requirements.
A Node is like a Point, but with more information.
virtual void renumber_node(dof_id_type old_id, dof_id_type new_id)=0
Changes the id of node old_id, both by changing node(old_id)->id() and by moving node(old_id) in the ...
SimplexRefiner(UnstructuredMesh &mesh)
The constructor.
virtual FunctionBase< Real > * get_desired_volume_function()
Get the function giving desired element volume as a function of position, or nullptr if no such funct...
std::vector< std::tuple< dof_id_type, dof_id_type, dof_id_type, processor_id_type > > refinement_datum
static constexpr Real TOLERANCE
std::map< std::pair< Node *, Node * >, Node * > new_nodes
Keep track of new nodes on edges.
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
This is the base class from which all geometric element types are derived.
const Parallel::Communicator & comm() const
Real & desired_volume()
Sets and/or gets the desired element volume.
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.
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
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.
uint8_t processor_id_type
virtual bool is_serial() const
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
This is the MeshCommunication class.
UnstructuredMesh & _mesh
Reference to the mesh which is to be refined.
std::vector< std::unique_ptr< Elem > > new_elements
Keep track of elements to add so we don't invalidate iterators during an iteration over elements...
The UnstructuredMesh class is derived from the MeshBase class.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
bool refine_elements()
Finds elements which exceed the requested metric and refines them via subdivision into new simplices ...
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
bool should_refine_elem(Elem &elem)
Checks if an element exceeds the requested metric.
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
void fill_refinement_datum(std::pair< Node *, Node *> vertices, refinement_datum &vec)
Helper for responding to edge queries.
virtual unsigned int n_sides() const =0
std::unordered_map< Elem *, dof_id_type > coarse_parent
Keep track of added elements' original coarse ids, so we get subsequent splits assigned to the correc...
const Elem * neighbor_ptr(unsigned int i) const
unsigned int level() const
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void renumber_elem(dof_id_type old_id, dof_id_type new_id)=0
Changes the id of element old_id, both by changing elem(old_id)->id() and by moving elem(old_id) in t...
void max(const T &r, T &o, Request &req) const
const Node * node_ptr(const unsigned int i) const
virtual bool is_replicated() const
void make_node_proc_ids_parallel_consistent(MeshBase &)
Assuming all processor ids on nodes touching local elements are parallel consistent, this function makes all other processor ids parallel consistent as well.
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...
virtual Real volume() const
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...
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
unsigned int n_extra_integers() const
Returns how many extra integers are associated to the DofObject.
std::unordered_map< processor_id_type, std::unique_ptr< Elem > > proxy_elements
Keep track of what processor an element's neighbors came from.
std::size_t refine_via_edges()
Finds elements which exceed the requested metric and refines them via inserting new midedge nodes and...
std::unordered_map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > edge_queries
Keep track of coarse remote edges for which we should query about additional point splits...
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
processor_id_type processor_id() const
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
const Point & point(const unsigned int i) const
Point vertex_average() const
dof_id_type get_extra_integer(const unsigned int index) const
Gets the value on this object of the extra integer associated with index, which should have been obta...
const RemoteElem * remote_elem