libMesh
mesh_subdivision_support.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 // Local includes
23 #include "libmesh/mesh_tools.h"
24 #include "libmesh/mesh_subdivision_support.h"
25 #include "libmesh/boundary_info.h"
26 
27 namespace libMesh
28 {
29 
30 
31 void MeshTools::Subdivision::find_one_ring(const Tri3Subdivision * elem,
32  std::vector<const Node *> & nodes)
33 {
34  libmesh_assert(elem->is_subdivision_updated());
35  libmesh_assert(elem->get_ordered_node(0));
36 
37  unsigned int valence = elem->get_ordered_valence(0);
38  nodes.resize(valence + 6);
39 
40  // The first three vertices in the patch are the ones from the element triangle
41  nodes[0] = elem->get_ordered_node(0);
42  nodes[1] = elem->get_ordered_node(1);
43  nodes[valence] = elem->get_ordered_node(2);
44 
45  const unsigned int nn0 = elem->local_node_number(nodes[0]->id());
46 
47  const Tri3Subdivision * nb = dynamic_cast<const Tri3Subdivision *>(elem->neighbor_ptr(nn0));
48  libmesh_assert(nb);
49 
50  unsigned int j, i = 1;
51 
52  do
53  {
54  ++i;
55  j = nb->local_node_number(nodes[0]->id());
56  nodes[i] = nb->node_ptr(next[j]);
57  nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(j));
58  } while (nb != elem);
59 
60  /* for nodes connected with N (= valence[0]) */
61  nb = static_cast<const Tri3Subdivision *>(elem->neighbor_ptr(next[nn0]));
62  j = nb->local_node_number(nodes[1]->id());
63  nodes[valence+1] = nb->node_ptr(next[j]);
64 
65  nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(next[j]));
66  j = nb->local_node_number(nodes[valence+1]->id());
67  nodes[valence+4] = nb->node_ptr(next[j]);
68 
69  nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(next[j]));
70  j = nb->local_node_number(nodes[valence+4]->id());
71  nodes[valence+5] = nb->node_ptr(next[j]);
72 
73  /* for nodes connected with 1 */
74  nb = static_cast<const Tri3Subdivision *>(elem->neighbor_ptr(next[nn0]));
75  j = nb->local_node_number(nodes[1]->id());
76  // nodes[valence+1] has been determined already
77 
78  nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(j));
79  j = nb->local_node_number(nodes[1]->id());
80  nodes[valence+2] = nb->node_ptr(next[j]);
81 
82  nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(j));
83  j = nb->local_node_number(nodes[1]->id());
84  nodes[valence+3] = nb->node_ptr(next[j]);
85 
86  return;
87 }
88 
89 
91 {
92  std::vector<Elem *> new_elements;
93  new_elements.reserve(mesh.n_elem());
94  const bool mesh_has_boundary_data =
95  (mesh.get_boundary_info().n_boundary_ids() > 0);
96 
97  std::vector<Elem *> new_boundary_elements;
98  std::vector<short int> new_boundary_sides;
99  std::vector<boundary_id_type> new_boundary_ids;
100 
101  // Container to catch ids handed back from BoundaryInfo
102  std::vector<boundary_id_type> ids;
103 
104  for (const auto & elem : mesh.element_ptr_range())
105  {
106  libmesh_assert_equal_to(elem->type(), TRI3);
107 
108  Elem * tri = new Tri3Subdivision;
109  tri->set_id(elem->id());
110  tri->subdomain_id() = elem->subdomain_id();
111  tri->set_node(0) = elem->node_ptr(0);
112  tri->set_node(1) = elem->node_ptr(1);
113  tri->set_node(2) = elem->node_ptr(2);
114 
115  if (mesh_has_boundary_data)
116  {
117  for (auto side : elem->side_index_range())
118  {
119  mesh.get_boundary_info().boundary_ids(elem, side, ids);
120 
121  for (std::size_t id=0; id<ids.size(); ++id)
122  {
123  // add the boundary id to the list of new boundary ids
124  new_boundary_ids.push_back(ids[id]);
125  new_boundary_elements.push_back(tri);
126  new_boundary_sides.push_back(side);
127  }
128  }
129 
130  // remove the original element from the BoundaryInfo structure
131  mesh.get_boundary_info().remove(elem);
132  }
133 
134  new_elements.push_back(tri);
135  mesh.insert_elem(tri);
136  }
137  mesh.prepare_for_use();
138 
139  if (mesh_has_boundary_data)
140  {
141  // If the old mesh had boundary data, the new mesh better have some too.
142  libmesh_assert_greater(new_boundary_elements.size(), 0);
143 
144  // We should also be sure that the lengths of the new boundary data vectors
145  // are all the same.
146  libmesh_assert_equal_to(new_boundary_sides.size(), new_boundary_elements.size());
147  libmesh_assert_equal_to(new_boundary_sides.size(), new_boundary_ids.size());
148 
149  // Add the new boundary info to the mesh.
150  for (std::size_t s = 0; s < new_boundary_elements.size(); ++s)
151  mesh.get_boundary_info().add_side(new_boundary_elements[s],
152  new_boundary_sides[s],
153  new_boundary_ids[s]);
154  }
155 
156  mesh.prepare_for_use();
157 }
158 
159 
161 {
162  mesh.prepare_for_use();
163 
164  // convert all mesh elements to subdivision elements
165  all_subdivision(mesh);
166 
167  if (!ghosted)
168  {
169  // add the ghost elements for the boundaries
170  add_boundary_ghosts(mesh);
171  }
172  else
173  {
174  // This assumes that the mesh already has the ghosts. Only tagging them is required here.
175  tag_boundary_ghosts(mesh);
176  }
177 
178  mesh.prepare_for_use();
179 
180  std::vector<std::vector<const Elem *>> nodes_to_elem_map;
181  MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
182 
183  // compute the node valences
184  for (auto & node : mesh.node_ptr_range())
185  {
186  std::vector<const Node *> neighbors;
187  MeshTools::find_nodal_neighbors(mesh, *node, nodes_to_elem_map, neighbors);
188  const unsigned int valence =
189  cast_int<unsigned int>(neighbors.size());
190  libmesh_assert_greater(valence, 1);
191  node->set_valence(valence);
192  }
193 
194  for (auto & elem : mesh.element_ptr_range())
195  {
196  Tri3Subdivision * tri3s = dynamic_cast<Tri3Subdivision *>(elem);
197  libmesh_assert(tri3s);
198  if (!tri3s->is_ghost())
199  tri3s->prepare_subdivision_properties();
200  }
201 }
202 
203 
205 {
206  for (auto & elem : mesh.element_ptr_range())
207  {
208  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
209 
210  Tri3Subdivision * sd_elem = static_cast<Tri3Subdivision *>(elem);
211  for (auto i : elem->side_index_range())
212  {
213  if (elem->neighbor_ptr(i) == libmesh_nullptr)
214  {
215  sd_elem->set_ghost(true);
216  // set all other neighbors to ghosts as well
217  if (elem->neighbor_ptr(next[i]))
218  {
219  Tri3Subdivision * nb = static_cast<Tri3Subdivision *>(elem->neighbor_ptr(next[i]));
220  nb->set_ghost(true);
221  }
222  if (elem->neighbor_ptr(prev[i]))
223  {
224  Tri3Subdivision * nb = static_cast<Tri3Subdivision *>(elem->neighbor_ptr(prev[i]));
225  nb->set_ghost(true);
226  }
227  }
228  }
229  }
230 }
231 
232 
234 {
235  static const Real tol = 1e-5;
236 
237  // add the mirrored ghost elements (without using iterators, because the mesh is modified in the course)
238  std::vector<Tri3Subdivision *> ghost_elems;
239  std::vector<Node *> ghost_nodes;
240  const unsigned int n_elem = mesh.n_elem();
241  for (unsigned int eid = 0; eid < n_elem; ++eid)
242  {
243  Elem * elem = mesh.elem_ptr(eid);
244  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
245 
246  // If the triangle happens to be in a corner (two boundary
247  // edges), we perform a counter-clockwise loop by mirroring the
248  // previous triangle until we come back to the original
249  // triangle. This prevents degenerated triangles in the mesh
250  // corners and guarantees that the node in the middle of the
251  // loop is of valence=6.
252  for (auto i : elem->side_index_range())
253  {
254  libmesh_assert_not_equal_to(elem->neighbor_ptr(i), elem);
255 
256  if (elem->neighbor_ptr(i) == libmesh_nullptr &&
257  elem->neighbor_ptr(next[i]) == libmesh_nullptr)
258  {
259  Elem * nelem = elem;
260  unsigned int k = i;
261  for (unsigned int l=0;l<4;l++)
262  {
263  // this is the vertex to be mirrored
264  Point point = nelem->point(k) + nelem->point(next[k]) - nelem->point(prev[k]);
265 
266  // Check if the proposed vertex doesn't coincide
267  // with one of the existing vertices. This is
268  // necessary because for some triangulations, it can
269  // happen that two mirrored ghost vertices coincide,
270  // which would then lead to a zero size ghost
271  // element below.
272  Node * node = libmesh_nullptr;
273  for (std::size_t j = 0; j < ghost_nodes.size(); ++j)
274  {
275  if ((*ghost_nodes[j] - point).norm() < tol * (elem->point(k) - point).norm())
276  {
277  node = ghost_nodes[j];
278  break;
279  }
280  }
281 
282  // add the new vertex only if no other is nearby
283  if (node == libmesh_nullptr)
284  {
285  node = mesh.add_point(point);
286  ghost_nodes.push_back(node);
287  }
288 
289  Tri3Subdivision * newelem = new Tri3Subdivision();
290 
291  // add the first new ghost element to the list just as in the non-corner case
292  if (l == 0)
293  ghost_elems.push_back(newelem);
294 
295  newelem->set_node(0) = nelem->node_ptr(next[k]);
296  newelem->set_node(1) = nelem->node_ptr(k);
297  newelem->set_node(2) = node;
298  newelem->set_neighbor(0, nelem);
299  newelem->set_ghost(true);
300  if (l>0)
301  newelem->set_neighbor(2, libmesh_nullptr);
302  nelem->set_neighbor(k, newelem);
303 
304  mesh.add_elem(newelem);
305  mesh.get_boundary_info().add_node(nelem->node_ptr(k), 1);
306  mesh.get_boundary_info().add_node(nelem->node_ptr(next[k]), 1);
307  mesh.get_boundary_info().add_node(nelem->node_ptr(prev[k]), 1);
308  mesh.get_boundary_info().add_node(node, 1);
309 
310  nelem = newelem;
311  k = 2 ;
312  }
313 
314  Tri3Subdivision * newelem = new Tri3Subdivision();
315 
316  newelem->set_node(0) = elem->node_ptr(next[i]);
317  newelem->set_node(1) = nelem->node_ptr(2);
318  newelem->set_node(2) = elem->node_ptr(prev[i]);
319  newelem->set_neighbor(0, nelem);
320  nelem->set_neighbor(2, newelem);
321  newelem->set_ghost(true);
322  newelem->set_neighbor(2, elem);
323  elem->set_neighbor(next[i],newelem);
324 
325  mesh.add_elem(newelem);
326 
327  break;
328  }
329  }
330 
331  for (auto i : elem->side_index_range())
332  {
333  libmesh_assert_not_equal_to(elem->neighbor_ptr(i), elem);
334  if (elem->neighbor_ptr(i) == libmesh_nullptr)
335  {
336  // this is the vertex to be mirrored
337  Point point = elem->point(i) + elem->point(next[i]) - elem->point(prev[i]);
338 
339  // Check if the proposed vertex doesn't coincide with
340  // one of the existing vertices. This is necessary
341  // because for some triangulations, it can happen that
342  // two mirrored ghost vertices coincide, which would
343  // then lead to a zero size ghost element below.
344  Node * node = libmesh_nullptr;
345  for (std::size_t j = 0; j < ghost_nodes.size(); ++j)
346  {
347  if ((*ghost_nodes[j] - point).norm() < tol * (elem->point(i) - point).norm())
348  {
349  node = ghost_nodes[j];
350  break;
351  }
352  }
353 
354  // add the new vertex only if no other is nearby
355  if (node == libmesh_nullptr)
356  {
357  node = mesh.add_point(point);
358  ghost_nodes.push_back(node);
359  }
360 
361  Tri3Subdivision * newelem = new Tri3Subdivision();
362  ghost_elems.push_back(newelem);
363 
364  newelem->set_node(0) = elem->node_ptr(next[i]);
365  newelem->set_node(1) = elem->node_ptr(i);
366  newelem->set_node(2) = node;
367  newelem->set_neighbor(0, elem);
368  newelem->set_ghost(true);
369  elem->set_neighbor(i, newelem);
370 
371  mesh.add_elem(newelem);
372  mesh.get_boundary_info().add_node(elem->node_ptr(i), 1);
373  mesh.get_boundary_info().add_node(elem->node_ptr(next[i]), 1);
374  mesh.get_boundary_info().add_node(elem->node_ptr(prev[i]), 1);
375  mesh.get_boundary_info().add_node(node, 1);
376  }
377  }
378  }
379 
380  // add the missing ghost elements (connecting new ghost nodes)
381  std::vector<Tri3Subdivision *> missing_ghost_elems;
382  std::vector<Tri3Subdivision *>::iterator ghost_el = ghost_elems.begin();
383  const std::vector<Tri3Subdivision *>::iterator end_ghost_el = ghost_elems.end();
384  for (; ghost_el != end_ghost_el; ++ghost_el)
385  {
386  Tri3Subdivision * elem = *ghost_el;
387  libmesh_assert(elem->is_ghost());
388 
389  for (auto i : elem->side_index_range())
390  {
391  if (elem->neighbor_ptr(i) == libmesh_nullptr &&
392  elem->neighbor_ptr(prev[i]) != libmesh_nullptr)
393  {
394  // go around counter-clockwise
395  Tri3Subdivision * nb1 = static_cast<Tri3Subdivision *>(elem->neighbor_ptr(prev[i]));
396  Tri3Subdivision * nb2 = nb1;
397  unsigned int j = i;
398  unsigned int n_nb = 0;
399  while (nb1 != libmesh_nullptr && nb1->id() != elem->id())
400  {
401  j = nb1->local_node_number(elem->node_id(i));
402  nb2 = nb1;
403  nb1 = static_cast<Tri3Subdivision *>(nb1->neighbor_ptr(prev[j]));
404  libmesh_assert(nb1 == libmesh_nullptr || nb1->id() != nb2->id());
405  n_nb++;
406  }
407 
408  libmesh_assert_not_equal_to(nb2->id(), elem->id());
409 
410  // Above, we merged coinciding ghost vertices. Therefore, we need
411  // to exclude the case where there is no ghost element to add between
412  // these two (identical) ghost nodes.
413  if (elem->node_ptr(next[i])->id() == nb2->node_ptr(prev[j])->id())
414  break;
415 
416  // If the number of already present neighbors is less than 4, we add another extra element
417  // so that the node in the middle of the loop ends up being of valence=6.
418  // This case usually happens when the middle node corresponds to a corner of the original mesh,
419  // and the extra element below prevents degenerated triangles in the mesh corners.
420  if (n_nb < 4)
421  {
422  // this is the vertex to be mirrored
423  Point point = nb2->point(j) + nb2->point(prev[j]) - nb2->point(next[j]);
424 
425  // Check if the proposed vertex doesn't coincide with one of the existing vertices.
426  // This is necessary because for some triangulations, it can happen that two mirrored
427  // ghost vertices coincide, which would then lead to a zero size ghost element below.
428  Node * node = libmesh_nullptr;
429  for (std::size_t k = 0; k < ghost_nodes.size(); ++k)
430  {
431  if ((*ghost_nodes[k] - point).norm() < tol * (nb2->point(j) - point).norm())
432  {
433  node = ghost_nodes[k];
434  break;
435  }
436  }
437 
438  // add the new vertex only if no other is nearby
439  if (node == libmesh_nullptr)
440  {
441  node = mesh.add_point(point);
442  ghost_nodes.push_back(node);
443  }
444 
445  Tri3Subdivision * newelem = new Tri3Subdivision();
446 
447  newelem->set_node(0) = nb2->node_ptr(j);
448  newelem->set_node(1) = nb2->node_ptr(prev[j]);
449  newelem->set_node(2) = node;
450  newelem->set_neighbor(0, nb2);
451  newelem->set_neighbor(1, libmesh_nullptr);
452  newelem->set_ghost(true);
453  nb2->set_neighbor(prev[j], newelem);
454 
455  mesh.add_elem(newelem);
456  mesh.get_boundary_info().add_node(nb2->node_ptr(j), 1);
457  mesh.get_boundary_info().add_node(nb2->node_ptr(prev[j]), 1);
458  mesh.get_boundary_info().add_node(node, 1);
459 
460  nb2 = newelem;
461  j = nb2->local_node_number(elem->node_id(i));
462  }
463 
464  Tri3Subdivision * newelem = new Tri3Subdivision();
465  newelem->set_node(0) = elem->node_ptr(next[i]);
466  newelem->set_node(1) = elem->node_ptr(i);
467  newelem->set_node(2) = nb2->node_ptr(prev[j]);
468  newelem->set_neighbor(0, elem);
469  newelem->set_neighbor(1, nb2);
470  newelem->set_neighbor(2, libmesh_nullptr);
471  newelem->set_ghost(true);
472 
473  elem->set_neighbor(i, newelem);
474  nb2->set_neighbor(prev[j], newelem);
475 
476  missing_ghost_elems.push_back(newelem);
477  break;
478  }
479  } // end side loop
480  } // end ghost element loop
481 
482  // add the missing ghost elements to the mesh
483  std::vector<Tri3Subdivision *>::iterator missing_el = missing_ghost_elems.begin();
484  const std::vector<Tri3Subdivision *>::iterator end_missing_el = missing_ghost_elems.end();
485  for (; missing_el != end_missing_el; ++missing_el)
486  mesh.add_elem(*missing_el);
487 }
488 
489 } // namespace libMesh
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
A Node is like a Point, but with more information.
Definition: node.h:52
virtual ElemType type() const =0
void add_boundary_ghosts(MeshBase &mesh)
Adds a new layer of "ghost" elements along the domain boundaries.
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:656
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2083
void tag_boundary_ghosts(MeshBase &mesh)
Flags the outermost element layer along the domain boundaries as "ghost" elements.
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
unsigned short int side
Definition: xdr_io.C:49
static const unsigned int prev[3]
A lookup table for the decrement modulo 3 operation, for iterating through the three nodes per elemen...
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
void all_subdivision(MeshBase &mesh)
Turns a triangulated mesh into a subdivision mesh.
void build_nodes_to_elem_map(const MeshBase &mesh, std::vector< std::vector< dof_id_type >> &nodes_to_elem_map)
After calling this function the input vector nodes_to_elem_map will contain the node to element conne...
Definition: mesh_tools.C:257
The libMesh namespace provides an interface to certain functionality in the library.
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.
This is the MeshBase class.
Definition: mesh_base.h:68
dof_id_type & set_id()
Definition: dof_object.h:641
libmesh_assert(j)
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
std::vector< boundary_id_type > boundary_ids(const Node *node) const
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
void find_nodal_neighbors(const MeshBase &mesh, const Node &n, const std::vector< std::vector< const Elem * >> &nodes_to_elem_map, std::vector< const Node * > &neighbors)
Given a mesh and a node in the mesh, the vector will be filled with every node directly attached to t...
Definition: mesh_tools.C:699
virtual SimpleRange< element_iterator > element_ptr_range()=0
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1874
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
virtual SimpleRange< node_iterator > node_ptr_range()=0
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
void set_neighbor(const unsigned int i, Elem *n)
Assigns n as the neighbor.
Definition: elem.h:2000
virtual Elem * insert_elem(Elem *e)=0
Insert elem e to the element array, preserving its id and replacing/deleting any existing element wit...
static const unsigned int next[3]
A lookup table for the increment modulo 3 operation, for iterating through the three nodes per elemen...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
void find_one_ring(const Tri3Subdivision *elem, std::vector< const Node * > &nodes)
Determines the 1-ring of element elem, and writes it to the nodes vector.
std::size_t n_boundary_ids() const
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 prepare_subdivision_mesh(MeshBase &mesh, bool ghosted=false)
Prepares the mesh for use with subdivision elements.
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual dof_id_type n_elem() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38