mesh_modification.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for std::acos()
23 #include <algorithm>
24 #include <limits>
25 #include <map>
26 
27 // Local includes
28 #include "libmesh/boundary_info.h"
29 #include "libmesh/face_tri3.h"
30 #include "libmesh/face_tri6.h"
32 #include "libmesh/location_maps.h"
35 #include "libmesh/mesh_tools.h"
36 #include "libmesh/parallel.h"
37 #include "libmesh/remote_elem.h"
38 #include "libmesh/string_to_enum.h"
40 
41 namespace libMesh
42 {
43 
44 
45 // ------------------------------------------------------------
46 // MeshTools::Modification functions for mesh modification
48  const Real factor,
49  const bool perturb_boundary)
50 {
51  libmesh_assert (mesh.n_nodes());
52  libmesh_assert (mesh.n_elem());
53  libmesh_assert ((factor >= 0.) && (factor <= 1.));
54 
55  START_LOG("distort()", "MeshTools::Modification");
56 
57 
58 
59  // First find nodes on the boundary and flag them
60  // so that we don't move them
61  // on_boundary holds false (not on boundary) and true (on boundary)
62  std::vector<bool> on_boundary (mesh.n_nodes(), false);
63 
64  if (!perturb_boundary) MeshTools::find_boundary_nodes (mesh, on_boundary);
65 
66  // Now calculate the minimum distance to
67  // neighboring nodes for each node.
68  // hmin holds these distances.
69  std::vector<float> hmin (mesh.n_nodes(),
71 
74 
75  for (; el!=end; ++el)
76  for (unsigned int n=0; n<(*el)->n_nodes(); n++)
77  hmin[(*el)->node(n)] = std::min(hmin[(*el)->node(n)],
78  static_cast<float>((*el)->hmin()));
79 
80 
81  // Now actually move the nodes
82  {
83  const unsigned int seed = 123456;
84 
85  // seed the random number generator.
86  // We'll loop from 1 to n_nodes on every processor, even those
87  // that don't have a particular node, so that the pseudorandom
88  // numbers will be the same everywhere.
89  std::srand(seed);
90 
91  // If the node is on the boundary or
92  // the node is not used by any element (hmin[n]<1.e20)
93  // then we should not move it.
94  // [Note: Testing for (in)equality might be wrong
95  // (different types, namely float and double)]
96  for (unsigned int n=0; n<mesh.n_nodes(); n++)
97  if (!on_boundary[n] && (hmin[n] < 1.e20) )
98  {
99  // the direction, random but unit normalized
100 
101  Point dir( static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
102  (mesh.mesh_dimension() > 1) ?
103  static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
104  : 0.,
105  ((mesh.mesh_dimension() == 3) ?
106  static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
107  : 0.)
108  );
109 
110  dir(0) = (dir(0)-.5)*2.;
111  if (mesh.mesh_dimension() > 1)
112  dir(1) = (dir(1)-.5)*2.;
113  if (mesh.mesh_dimension() == 3)
114  dir(2) = (dir(2)-.5)*2.;
115 
116  dir = dir.unit();
117 
118  Node *node = mesh.node_ptr(n);
119  if (!node)
120  continue;
121 
122  (*node)(0) += dir(0)*factor*hmin[n];
123  if (mesh.mesh_dimension() > 1)
124  (*node)(1) += dir(1)*factor*hmin[n];
125  if (mesh.mesh_dimension() == 3)
126  (*node)(2) += dir(2)*factor*hmin[n];
127  }
128  }
129 
130 
131  // All done
132  STOP_LOG("distort()", "MeshTools::Modification");
133 }
134 
135 
136 
138  const Real xt,
139  const Real yt,
140  const Real zt)
141 {
142  const Point p(xt, yt, zt);
143 
144  const MeshBase::node_iterator nd_end = mesh.nodes_end();
145 
146  for (MeshBase::node_iterator nd = mesh.nodes_begin();
147  nd != nd_end; ++nd)
148  **nd += p;
149 }
150 
151 
152 // void MeshTools::Modification::rotate2D (MeshBase& mesh,
153 // const Real alpha)
154 // {
155 // libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
156 
157 // const Real pi = std::acos(-1);
158 // const Real a = alpha/180.*pi;
159 // for (unsigned int n=0; n<mesh.n_nodes(); n++)
160 // {
161 // const Point p = mesh.node(n);
162 // const Real x = p(0);
163 // const Real y = p(1);
164 // const Real z = p(2);
165 // mesh.node(n) = Point(std::cos(a)*x - std::sin(a)*y,
166 // std::sin(a)*x + std::cos(a)*y,
167 // z);
168 // }
169 
170 // }
171 
172 
173 
175  const Real phi,
176  const Real theta,
177  const Real psi)
178 {
179  libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
180 
181  const Real p = -phi/180.*libMesh::pi;
182  const Real t = -theta/180.*libMesh::pi;
183  const Real s = -psi/180.*libMesh::pi;
184  const Real sp = std::sin(p), cp = std::cos(p);
185  const Real st = std::sin(t), ct = std::cos(t);
186  const Real ss = std::sin(s), cs = std::cos(s);
187 
188  // We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
189  // (equations 6-14 give the entries of the composite transformation matrix).
190  // The rotations are performed sequentially about the z, x, and z axes, in that order.
191  // A positive angle yields a counter-clockwise rotation about the axis in question.
192  const MeshBase::node_iterator nd_end = mesh.nodes_end();
193 
194  for (MeshBase::node_iterator nd = mesh.nodes_begin();
195  nd != nd_end; ++nd)
196  {
197  const Point pt = **nd;
198  const Real x = pt(0);
199  const Real y = pt(1);
200  const Real z = pt(2);
201  **nd = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
202  (-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
203  ( sp*st)*x + (-cp*st)*y + (ct)*z );
204  }
205 }
206 
207 
209  const Real xs,
210  const Real ys,
211  const Real zs)
212 {
213  const Real x_scale = xs;
214  Real y_scale = ys;
215  Real z_scale = zs;
216 
217  if (ys == 0.)
218  {
219  libmesh_assert_equal_to (zs, 0.);
220 
221  y_scale = z_scale = x_scale;
222  }
223 
224  // Scale the x coordinate in all dimensions
225  const MeshBase::node_iterator nd_end = mesh.nodes_end();
226 
227  for (MeshBase::node_iterator nd = mesh.nodes_begin();
228  nd != nd_end; ++nd)
229  (**nd)(0) *= x_scale;
230 
231 
232  // Only scale the y coordinate in 2 and 3D
233  if (mesh.spatial_dimension() < 2)
234  return;
235 
236  for (MeshBase::node_iterator nd = mesh.nodes_begin();
237  nd != nd_end; ++nd)
238  (**nd)(1) *= y_scale;
239 
240  // Only scale the z coordinate in 3D
241  if (mesh.spatial_dimension() < 3)
242  return;
243 
244  for (MeshBase::node_iterator nd = mesh.nodes_begin();
245  nd != nd_end; ++nd)
246  (**nd)(2) *= z_scale;
247 }
248 
249 
250 
251 
252 // ------------------------------------------------------------
253 // UnstructuredMesh class member functions for mesh modification
255 {
256  /*
257  * when the mesh is not prepared,
258  * at least renumber the nodes and
259  * elements, so that the node ids
260  * are correct
261  */
262  if (!this->_is_prepared)
264 
265  START_LOG("all_first_order()", "Mesh");
266 
270  std::vector<bool> node_touched_by_me(this->max_node_id(), false);
271 
277  element_iterator endit = elements_end();
278  for (element_iterator it = elements_begin();
279  it != endit; ++it)
280  {
281  Elem* so_elem = *it;
282 
283  libmesh_assert(so_elem);
284 
285  /*
286  * build the first-order equivalent, add to
287  * the new_elements list.
288  */
289  Elem* lo_elem = Elem::build
291  (so_elem->type()), so_elem->parent()).release();
292 
293  for (unsigned int s=0; s != so_elem->n_sides(); ++s)
294  if (so_elem->neighbor(s) == remote_elem)
295  lo_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
296 
297 #ifdef LIBMESH_ENABLE_AMR
298  /*
299  * Reset the parent links of any child elements
300  */
301  if (so_elem->has_children())
302  for (unsigned int c=0; c != so_elem->n_children(); ++c)
303  {
304  so_elem->child(c)->set_parent(lo_elem);
305  lo_elem->add_child(so_elem->child(c), c);
306  }
307 
308  /*
309  * Reset the child link of any parent element
310  */
311  if (so_elem->parent())
312  {
313  unsigned int c =
314  so_elem->parent()->which_child_am_i(so_elem);
315  lo_elem->parent()->replace_child(lo_elem, c);
316  }
317 
318  /*
319  * Copy as much data to the new element as makes sense
320  */
321  lo_elem->set_p_level(so_elem->p_level());
322  lo_elem->set_refinement_flag(so_elem->refinement_flag());
323  lo_elem->set_p_refinement_flag(so_elem->p_refinement_flag());
324 #endif
325 
326  libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
327 
328  /*
329  * By definition the vertices of the linear and
330  * second order element are identically numbered.
331  * transfer these.
332  */
333  for (unsigned int v=0; v < so_elem->n_vertices(); v++)
334  {
335  lo_elem->set_node(v) = so_elem->get_node(v);
336  node_touched_by_me[lo_elem->node(v)] = true;
337  }
338 
345  libmesh_assert_equal_to (lo_elem->n_sides(), so_elem->n_sides());
346 
347  for (unsigned int s=0; s<so_elem->n_sides(); s++)
348  {
349  const std::vector<boundary_id_type> boundary_ids =
350  this->boundary_info->raw_boundary_ids (so_elem, s);
351 
352  this->boundary_info->add_side (lo_elem, s, boundary_ids);
353  }
354 
355  /*
356  * The new first-order element is ready.
357  * Inserting it into the mesh will replace and delete
358  * the second-order element.
359  */
360  lo_elem->set_id(so_elem->id());
361  lo_elem->processor_id() = so_elem->processor_id();
362  lo_elem->subdomain_id() = so_elem->subdomain_id();
363  this->insert_elem(lo_elem);
364  }
365 
366  const MeshBase::node_iterator nd_end = this->nodes_end();
368  while (nd != nd_end)
369  {
370  Node *the_node = *nd;
371  ++nd;
372  if (!node_touched_by_me[the_node->id()])
373  this->delete_node(the_node);
374  }
375 
376  STOP_LOG("all_first_order()", "Mesh");
377 
378  // On hanging nodes that used to also be second order nodes, we
379  // might now have an invalid nodal processor_id()
381 
382  // delete or renumber nodes, etc
383  this->prepare_for_use(/*skip_renumber =*/ false);
384 }
385 
386 
387 
388 void UnstructuredMesh::all_second_order (const bool full_ordered)
389 {
390  // This function must be run on all processors at once
391  parallel_object_only();
392 
393  /*
394  * when the mesh is not prepared,
395  * at least renumber the nodes and
396  * elements, so that the node ids
397  * are correct
398  */
399  if (!this->_is_prepared)
401 
402  /*
403  * If the mesh is empty
404  * then we have nothing to do
405  */
406  if (!this->n_elem())
407  return;
408 
409  /*
410  * If the mesh is already second order
411  * then we have nothing to do.
412  * We have to test for this in a round-about way to avoid
413  * a bug on distributed parallel meshes with more processors
414  * than elements.
415  */
416  bool already_second_order = false;
417  if (this->elements_begin() != this->elements_end() &&
418  (*(this->elements_begin()))->default_order() != FIRST)
419  already_second_order = true;
420  this->comm().max(already_second_order);
421  if (already_second_order)
422  return;
423 
424  START_LOG("all_second_order()", "Mesh");
425 
426  /*
427  * this map helps in identifying second order
428  * nodes. Namely, a second-order node:
429  * - edge node
430  * - face node
431  * - bubble node
432  * is uniquely defined through a set of adjacent
433  * vertices. This set of adjacent vertices is
434  * used to identify already added higher-order
435  * nodes. We are safe to use node id's since we
436  * make sure that these are correctly numbered.
437  */
438  std::map<std::vector<dof_id_type>, Node*> adj_vertices_to_so_nodes;
439 
440  /*
441  * for speed-up of the \p add_point() method, we
442  * can reserve memory. Guess the number of additional
443  * nodes for different dimensions
444  */
445  switch (this->mesh_dimension())
446  {
447  case 1:
448  /*
449  * in 1D, there can only be order-increase from Edge2
450  * to Edge3. Something like 1/2 of n_nodes() have
451  * to be added
452  */
453  this->reserve_nodes(static_cast<unsigned int>
454  (1.5*static_cast<double>(this->n_nodes())));
455  break;
456 
457  case 2:
458  /*
459  * in 2D, either refine from Tri3 to Tri6 (double the nodes)
460  * or from Quad4 to Quad8 (again, double) or Quad9 (2.25 that much)
461  */
462  this->reserve_nodes(static_cast<unsigned int>
463  (2*static_cast<double>(this->n_nodes())));
464  break;
465 
466 
467  case 3:
468  /*
469  * in 3D, either refine from Tet4 to Tet10 (factor = 2.5) up to
470  * Hex8 to Hex27 (something > 3). Since in 3D there _are_ already
471  * quite some nodes, and since we do not want to overburden the memory by
472  * a too conservative guess, use the lower bound
473  */
474  this->reserve_nodes(static_cast<unsigned int>
475  (2.5*static_cast<double>(this->n_nodes())));
476  break;
477 
478  default:
479  // Hm?
480  libmesh_error();
481  }
482 
483 
484 
485  /*
486  * form a vector that will hold the node id's of
487  * the vertices that are adjacent to the son-th
488  * second-order node. Pull this outside of the
489  * loop so that silly compilers don't repeatedly
490  * create and destroy the vector.
491  */
492  std::vector<dof_id_type> adjacent_vertices_ids;
493 
502  it != endit; ++it)
503  {
504  // the linear-order element
505  const Elem* lo_elem = *it;
506 
507  libmesh_assert(lo_elem);
508 
509  // make sure it is linear order
510  if (lo_elem->default_order() != FIRST)
511  {
512  libMesh::err << "ERROR: This is not a linear element: type="
513  << lo_elem->type() << std::endl;
514  libmesh_error();
515  }
516 
517  // this does _not_ work for refined elements
518  libmesh_assert_equal_to (lo_elem->level (), 0);
519 
520  /*
521  * build the second-order equivalent, add to
522  * the new_elements list. Note that this here
523  * is the only point where \p full_ordered
524  * is necessary. The remaining code works well
525  * for either type of seconrd-order equivalent, e.g.
526  * Hex20 or Hex27, as equivalents for Hex8
527  */
528  Elem* so_elem =
530  full_ordered) ).release();
531 
532  libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
533 
534 
535  /*
536  * By definition the vertices of the linear and
537  * second order element are identically numbered.
538  * transfer these.
539  */
540  for (unsigned int v=0; v < lo_elem->n_vertices(); v++)
541  so_elem->set_node(v) = lo_elem->get_node(v);
542 
543  /*
544  * Now handle the additional mid-side nodes. This
545  * is simply handled through a map that remembers
546  * the already-added nodes. This map maps the global
547  * ids of the vertices (that uniquely define this
548  * higher-order node) to the new node.
549  * Notation: son = second-order node
550  */
551  const unsigned int son_begin = so_elem->n_vertices();
552  const unsigned int son_end = so_elem->n_nodes();
553 
554 
555  for (unsigned int son=son_begin; son<son_end; son++)
556  {
557  const unsigned int n_adjacent_vertices =
558  so_elem->n_second_order_adjacent_vertices(son);
559 
560  adjacent_vertices_ids.resize(n_adjacent_vertices);
561 
562  for (unsigned int v=0; v<n_adjacent_vertices; v++)
563  adjacent_vertices_ids[v] =
564  so_elem->node( so_elem->second_order_adjacent_vertex(son,v) );
565 
566  /*
567  * \p adjacent_vertices_ids is now in order of the current
568  * side. sort it, so that comparisons with the
569  * \p adjacent_vertices_ids created through other elements'
570  * sides can match
571  */
572  std::sort(adjacent_vertices_ids.begin(),
573  adjacent_vertices_ids.end());
574 
575 
576  // does this set of vertices already has a mid-node added?
577  std::pair<std::map<std::vector<dof_id_type>, Node*>::iterator,
578  std::map<std::vector<dof_id_type>, Node*>::iterator>
579  pos = adj_vertices_to_so_nodes.equal_range (adjacent_vertices_ids);
580 
581  // no, not added yet
582  if (pos.first == pos.second)
583  {
584  /*
585  * for this set of vertices, there is no
586  * second_order node yet. Add it.
587  *
588  * compute the location of the new node as
589  * the average over the adjacent vertices.
590  */
591  Point new_location = this->point(adjacent_vertices_ids[0]);
592  for (unsigned int v=1; v<n_adjacent_vertices; v++)
593  new_location += this->point(adjacent_vertices_ids[v]);
594 
595  new_location /= static_cast<Real>(n_adjacent_vertices);
596 
597  /* Add the new point to the mesh, giving it a globally
598  * well-defined processor id.
599  */
600  Node* so_node = this->add_point
601  (new_location, DofObject::invalid_id,
602  this->node(adjacent_vertices_ids[0]).processor_id());
603 
604  /*
605  * insert the new node with its defining vertex
606  * set into the map, and relocate pos to this
607  * new entry, so that the so_elem can use
608  * \p pos for inserting the node
609  */
610  adj_vertices_to_so_nodes.insert(pos.first,
611  std::make_pair(adjacent_vertices_ids,
612  so_node));
613 
614  so_elem->set_node(son) = so_node;
615  }
616  // yes, already added.
617  else
618  {
619  libmesh_assert(pos.first->second);
620 
621  so_elem->set_node(son) = pos.first->second;
622  }
623  }
624 
625 
637  libmesh_assert_equal_to (lo_elem->n_sides(), so_elem->n_sides());
638 
639  for (unsigned int s=0; s<lo_elem->n_sides(); s++)
640  {
641  const std::vector<boundary_id_type> boundary_ids =
642  this->boundary_info->raw_boundary_ids (lo_elem, s);
643 
644  this->boundary_info->add_side (so_elem, s, boundary_ids);
645 
646  if (lo_elem->neighbor(s) == remote_elem)
647  so_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
648  }
649 
650  /*
651  * The new second-order element is ready.
652  * Inserting it into the mesh will replace and delete
653  * the first-order element.
654  */
655  so_elem->set_id(lo_elem->id());
656  so_elem->processor_id() = lo_elem->processor_id();
657  so_elem->subdomain_id() = lo_elem->subdomain_id();
658  this->insert_elem(so_elem);
659  }
660 
661  // we can clear the map
662  adj_vertices_to_so_nodes.clear();
663 
664 
665  STOP_LOG("all_second_order()", "Mesh");
666 
667  // In a ParallelMesh our ghost node processor ids may be bad and
668  // the ids of nodes touching remote elements may be inconsistent.
669  // Fix them.
670  if (!this->is_serial())
671  {
672  LocationMap<Node> loc_map;
674  (*this, loc_map);
675  }
676 
677  // renumber nodes, elements etc
678  this->prepare_for_use(/*skip_renumber =*/ false);
679 }
680 
681 
682 
683 
684 
685 
687 {
688  // The number of elements in the original mesh before any additions
689  // or deletions.
690  const dof_id_type n_orig_elem = mesh.n_elem();
691 
692  // We store pointers to the newly created elements in a vector
693  // until they are ready to be added to the mesh. This is because
694  // adding new elements on the fly can cause reallocation and invalidation
695  // of existing iterators.
696  std::vector<Elem*> new_elements;
697  new_elements.reserve (2*n_orig_elem);
698 
699  // If the original mesh has boundary data, we carry that over
700  // to the new mesh with triangular elements.
701  const bool mesh_has_boundary_data = (mesh.boundary_info->n_boundary_ids() > 0);
702 
703  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
704  std::vector<Elem*> new_bndry_elements;
705  std::vector<unsigned short int> new_bndry_sides;
706  std::vector<boundary_id_type> new_bndry_ids;
707 
708  // Iterate over the elements, splitting QUADS into
709  // pairs of conforming triangles.
710  // FIXME: This algorithm does not work on refined grids!
711  {
714 
715  for (; el!=end; ++el)
716  {
717  Elem* elem = *el;
718 
719  const ElemType etype = elem->type();
720 
721  // all_tri currently only works on coarse meshes
722  libmesh_assert (!elem->parent());
723 
724  // We split the quads using the shorter of the two diagonals
725  // to maintain the best angle properties.
726  bool edge_swap = false;
727 
728  // True if we actually split the current element.
729  bool split_elem = false;
730 
731  // The two new triangular elements we will split the quad into.
732  Elem* tri0 = NULL;
733  Elem* tri1 = NULL;
734 
735 
736  switch (etype)
737  {
738  case QUAD4:
739  {
740  split_elem = true;
741 
742  tri0 = new Tri3;
743  tri1 = new Tri3;
744 
745  // Check for possible edge swap
746  if ((elem->point(0) - elem->point(2)).size() <
747  (elem->point(1) - elem->point(3)).size())
748  {
749  tri0->set_node(0) = elem->get_node(0);
750  tri0->set_node(1) = elem->get_node(1);
751  tri0->set_node(2) = elem->get_node(2);
752 
753  tri1->set_node(0) = elem->get_node(0);
754  tri1->set_node(1) = elem->get_node(2);
755  tri1->set_node(2) = elem->get_node(3);
756  }
757 
758  else
759  {
760  edge_swap=true;
761 
762  tri0->set_node(0) = elem->get_node(0);
763  tri0->set_node(1) = elem->get_node(1);
764  tri0->set_node(2) = elem->get_node(3);
765 
766  tri1->set_node(0) = elem->get_node(1);
767  tri1->set_node(1) = elem->get_node(2);
768  tri1->set_node(2) = elem->get_node(3);
769  }
770 
771 
772  break;
773  }
774 
775  case QUAD8:
776  {
777  split_elem = true;
778 
779  tri0 = new Tri6;
780  tri1 = new Tri6;
781 
782  Node* new_node = mesh.add_point( (mesh.node(elem->node(0)) +
783  mesh.node(elem->node(1)) +
784  mesh.node(elem->node(2)) +
785  mesh.node(elem->node(3)) / 4)
786  );
787 
788  // Check for possible edge swap
789  if ((elem->point(0) - elem->point(2)).size() <
790  (elem->point(1) - elem->point(3)).size())
791  {
792  tri0->set_node(0) = elem->get_node(0);
793  tri0->set_node(1) = elem->get_node(1);
794  tri0->set_node(2) = elem->get_node(2);
795  tri0->set_node(3) = elem->get_node(4);
796  tri0->set_node(4) = elem->get_node(5);
797  tri0->set_node(5) = new_node;
798 
799  tri1->set_node(0) = elem->get_node(0);
800  tri1->set_node(1) = elem->get_node(2);
801  tri1->set_node(2) = elem->get_node(3);
802  tri1->set_node(3) = new_node;
803  tri1->set_node(4) = elem->get_node(6);
804  tri1->set_node(5) = elem->get_node(7);
805 
806  }
807 
808  else
809  {
810  edge_swap=true;
811 
812  tri0->set_node(0) = elem->get_node(3);
813  tri0->set_node(1) = elem->get_node(0);
814  tri0->set_node(2) = elem->get_node(1);
815  tri0->set_node(3) = elem->get_node(7);
816  tri0->set_node(4) = elem->get_node(4);
817  tri0->set_node(5) = new_node;
818 
819  tri1->set_node(0) = elem->get_node(1);
820  tri1->set_node(1) = elem->get_node(2);
821  tri1->set_node(2) = elem->get_node(3);
822  tri1->set_node(3) = elem->get_node(5);
823  tri1->set_node(4) = elem->get_node(6);
824  tri1->set_node(5) = new_node;
825  }
826 
827  break;
828  }
829 
830  case QUAD9:
831  {
832  split_elem = true;
833 
834  tri0 = new Tri6;
835  tri1 = new Tri6;
836 
837  // Check for possible edge swap
838  if ((elem->point(0) - elem->point(2)).size() <
839  (elem->point(1) - elem->point(3)).size())
840  {
841  tri0->set_node(0) = elem->get_node(0);
842  tri0->set_node(1) = elem->get_node(1);
843  tri0->set_node(2) = elem->get_node(2);
844  tri0->set_node(3) = elem->get_node(4);
845  tri0->set_node(4) = elem->get_node(5);
846  tri0->set_node(5) = elem->get_node(8);
847 
848  tri1->set_node(0) = elem->get_node(0);
849  tri1->set_node(1) = elem->get_node(2);
850  tri1->set_node(2) = elem->get_node(3);
851  tri1->set_node(3) = elem->get_node(8);
852  tri1->set_node(4) = elem->get_node(6);
853  tri1->set_node(5) = elem->get_node(7);
854  }
855 
856  else
857  {
858  edge_swap=true;
859 
860  tri0->set_node(0) = elem->get_node(0);
861  tri0->set_node(1) = elem->get_node(1);
862  tri0->set_node(2) = elem->get_node(3);
863  tri0->set_node(3) = elem->get_node(4);
864  tri0->set_node(4) = elem->get_node(8);
865  tri0->set_node(5) = elem->get_node(7);
866 
867  tri1->set_node(0) = elem->get_node(1);
868  tri1->set_node(1) = elem->get_node(2);
869  tri1->set_node(2) = elem->get_node(3);
870  tri1->set_node(3) = elem->get_node(5);
871  tri1->set_node(4) = elem->get_node(6);
872  tri1->set_node(5) = elem->get_node(8);
873  }
874 
875  break;
876  }
877  // No need to split elements that are already triangles
878  case TRI3:
879  case TRI6:
880  continue;
881  // Try to ignore non-2D elements for now
882  default:
883  {
884  libMesh::err << "Warning, encountered non-2D element "
885  << Utility::enum_to_string<ElemType>(etype)
886  << " in MeshTools::Modification::all_tri(), hope that's OK..."
887  << std::endl;
888  }
889  } // end switch (etype)
890 
891 
892 
893  if (split_elem)
894  {
895  // Be sure the correct ID's are also set for tri0 and
896  // tri1.
897  tri0->processor_id() = elem->processor_id();
898  tri0->subdomain_id() = elem->subdomain_id();
899  tri1->processor_id() = elem->processor_id();
900  tri1->subdomain_id() = elem->subdomain_id();
901 
902  if (mesh_has_boundary_data)
903  {
904  for (unsigned int sn=0; sn<elem->n_sides(); ++sn)
905  {
906  const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids(*el, sn);
907  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
908  {
909  const boundary_id_type b_id = *id_it;
910 
911  if (b_id != BoundaryInfo::invalid_id)
912  {
913  // Add the boundary ID to the list of new boundary ids
914  new_bndry_ids.push_back(b_id);
915 
916  // Convert the boundary side information of the old element to
917  // boundary side information for the new element.
918  if (!edge_swap)
919  {
920  switch (sn)
921  {
922  case 0:
923  {
924  // New boundary side is Tri 0, side 0
925  new_bndry_elements.push_back(tri0);
926  new_bndry_sides.push_back(0);
927  break;
928  }
929  case 1:
930  {
931  // New boundary side is Tri 0, side 1
932  new_bndry_elements.push_back(tri0);
933  new_bndry_sides.push_back(1);
934  break;
935  }
936  case 2:
937  {
938  // New boundary side is Tri 1, side 1
939  new_bndry_elements.push_back(tri1);
940  new_bndry_sides.push_back(1);
941  break;
942  }
943  case 3:
944  {
945  // New boundary side is Tri 1, side 2
946  new_bndry_elements.push_back(tri1);
947  new_bndry_sides.push_back(2);
948  break;
949  }
950 
951  default:
952  {
953  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
954  libmesh_error();
955  }
956  }
957  }
958 
959  else // edge_swap==true
960  {
961  switch (sn)
962  {
963  case 0:
964  {
965  // New boundary side is Tri 0, side 0
966  new_bndry_elements.push_back(tri0);
967  new_bndry_sides.push_back(0);
968  break;
969  }
970  case 1:
971  {
972  // New boundary side is Tri 1, side 0
973  new_bndry_elements.push_back(tri1);
974  new_bndry_sides.push_back(0);
975  break;
976  }
977  case 2:
978  {
979  // New boundary side is Tri 1, side 1
980  new_bndry_elements.push_back(tri1);
981  new_bndry_sides.push_back(1);
982  break;
983  }
984  case 3:
985  {
986  // New boundary side is Tri 0, side 2
987  new_bndry_elements.push_back(tri0);
988  new_bndry_sides.push_back(2);
989  break;
990  }
991 
992  default:
993  {
994  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
995  libmesh_error();
996  }
997  }
998  } // end edge_swap==true
999  } // end if (b_id != BoundaryInfo::invalid_id)
1000  } // end for loop over boundary IDs
1001  } // end for loop over sides
1002 
1003  // Remove the original element from the BoundaryInfo structure.
1004  mesh.boundary_info->remove(elem);
1005 
1006  } // end if (mesh_has_boundary_data)
1007 
1008 
1009  // On a distributed mesh, we need to preserve remote_elem
1010  // links, since prepare_for_use can't reconstruct them for
1011  // us.
1012  for (unsigned int sn=0; sn<elem->n_sides(); ++sn)
1013  {
1014  if (elem->neighbor(sn) == remote_elem)
1015  {
1016  // Create a remote_elem link on one of the new
1017  // elements corresponding to the link from the old
1018  // element.
1019  if (!edge_swap)
1020  {
1021  switch (sn)
1022  {
1023  case 0:
1024  {
1025  // New remote side is Tri 0, side 0
1026  tri0->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
1027  break;
1028  }
1029  case 1:
1030  {
1031  // New remote side is Tri 0, side 1
1032  tri0->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
1033  break;
1034  }
1035  case 2:
1036  {
1037  // New remote side is Tri 1, side 1
1038  tri1->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
1039  break;
1040  }
1041  case 3:
1042  {
1043  // New remote side is Tri 1, side 2
1044  tri1->set_neighbor(2, const_cast<RemoteElem*>(remote_elem));
1045  break;
1046  }
1047 
1048  default:
1049  {
1050  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
1051  libmesh_error();
1052  }
1053  }
1054  }
1055 
1056  else // edge_swap==true
1057  {
1058  switch (sn)
1059  {
1060  case 0:
1061  {
1062  // New remote side is Tri 0, side 0
1063  tri0->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
1064  break;
1065  }
1066  case 1:
1067  {
1068  // New remote side is Tri 1, side 0
1069  tri1->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
1070  break;
1071  }
1072  case 2:
1073  {
1074  // New remote side is Tri 1, side 1
1075  tri1->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
1076  break;
1077  }
1078  case 3:
1079  {
1080  // New remote side is Tri 0, side 2
1081  tri0->set_neighbor(2, const_cast<RemoteElem*>(remote_elem));
1082  break;
1083  }
1084 
1085  default:
1086  {
1087  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
1088  libmesh_error();
1089  }
1090  }
1091  } // end edge_swap==true
1092  } // end if (elem->neighbor(sn) == remote_elem)
1093  } // end for loop over sides
1094 
1095  // Determine new IDs for the split elements which will be
1096  // the same on all processors, therefore keeping the Mesh
1097  // in sync. Note: we offset the new IDs by n_orig_elem to
1098  // avoid overwriting any of the original IDs, this assumes
1099  // they were contiguously-numbered to begin with...
1100  tri0->set_id( n_orig_elem + 2*elem->id() + 0 );
1101  tri1->set_id( n_orig_elem + 2*elem->id() + 1 );
1102 
1103  // Add the newly-created triangles to the temporary vector of new elements.
1104  new_elements.push_back(tri0);
1105  new_elements.push_back(tri1);
1106 
1107  // Delete the original element
1108  mesh.delete_elem(elem);
1109  } // end if (split_elem)
1110  } // End for loop over elements
1111  } // end scope
1112 
1113 
1114  // Now, iterate over the new elements vector, and add them each to
1115  // the Mesh.
1116  {
1117  std::vector<Elem*>::iterator el = new_elements.begin();
1118  const std::vector<Elem*>::iterator end = new_elements.end();
1119  for (; el != end; ++el)
1120  mesh.add_elem(*el);
1121  }
1122 
1123  if (mesh_has_boundary_data)
1124  {
1125  // By this time, we should have removed all of the original boundary sides
1126  // - except on a hybrid mesh, where we can't "start from a blank slate"! - RHS
1127  // libmesh_assert_equal_to (mesh.boundary_info->n_boundary_conds(), 0);
1128 
1129  // Clear the boundary info, to be sure and start from a blank slate.
1130  // mesh.boundary_info->clear();
1131 
1132  // If the old mesh had boundary data, the new mesh better have some.
1133  libmesh_assert_greater (new_bndry_elements.size(), 0);
1134 
1135  // We should also be sure that the lengths of the new boundary data vectors
1136  // are all the same.
1137  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1138  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1139 
1140  // Add the new boundary info to the mesh
1141  for (unsigned int s=0; s<new_bndry_elements.size(); ++s)
1142  mesh.boundary_info->add_side(new_bndry_elements[s],
1143  new_bndry_sides[s],
1144  new_bndry_ids[s]);
1145  }
1146 
1147 
1148  // Prepare the newly created mesh for use.
1149  mesh.prepare_for_use(/*skip_renumber =*/ false);
1150 
1151  // Let the new_elements and new_bndry_elements vectors go out of scope.
1152 }
1153 
1154 
1156  const unsigned int n_iterations,
1157  const Real power)
1158 {
1162  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1163 
1164  /*
1165  * find the boundary nodes
1166  */
1167  std::vector<bool> on_boundary;
1168  MeshTools::find_boundary_nodes(mesh, on_boundary);
1169 
1170  for (unsigned int iter=0; iter<n_iterations; iter++)
1171 
1172  {
1173  /*
1174  * loop over the mesh refinement level
1175  */
1176  unsigned int n_levels = MeshTools::n_levels(mesh);
1177  for (unsigned int refinement_level=0; refinement_level != n_levels;
1178  refinement_level++)
1179  {
1180  // initialize the storage (have to do it on every level to get empty vectors
1181  std::vector<Point> new_positions;
1182  std::vector<Real> weight;
1183  new_positions.resize(mesh.n_nodes());
1184  weight.resize(mesh.n_nodes());
1185 
1186  {
1187  /*
1188  * Loop over the elements to calculate new node positions
1189  */
1190  MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
1191  const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
1192 
1193  for (; el != end; ++el)
1194  {
1195  /*
1196  * Constant handle for the element
1197  */
1198  const Elem* elem = *el;
1199 
1200  /*
1201  * We relax all nodes on level 0 first
1202  * If the element is refined (level > 0), we interpolate the
1203  * parents nodes with help of the embedding matrix
1204  */
1205  if (refinement_level == 0)
1206  {
1207  for (unsigned int s=0; s<elem->n_neighbors(); s++)
1208  {
1209  /*
1210  * Only operate on sides which are on the
1211  * boundary or for which the current element's
1212  * id is greater than its neighbor's.
1213  * Sides get only built once.
1214  */
1215  if ((elem->neighbor(s) != NULL) &&
1216  (elem->id() > elem->neighbor(s)->id()) )
1217  {
1218  AutoPtr<Elem> side(elem->build_side(s));
1219 
1220  Node* node0 = side->get_node(0);
1221  Node* node1 = side->get_node(1);
1222 
1223  Real node_weight = 1.;
1224  // calculate the weight of the nodes
1225  if (power > 0)
1226  {
1227  Point diff = (*node0)-(*node1);
1228  node_weight = std::pow( diff.size(), power );
1229  }
1230 
1231  const dof_id_type id0 = node0->id(), id1 = node1->id();
1232  new_positions[id0].add_scaled( *node1, node_weight );
1233  new_positions[id1].add_scaled( *node0, node_weight );
1234  weight[id0] += node_weight;
1235  weight[id1] += node_weight;
1236  }
1237  } // element neighbor loop
1238  }
1239 #ifdef LIBMESH_ENABLE_AMR
1240  else // refinement_level > 0
1241  {
1242  /*
1243  * Find the positions of the hanging nodes of refined elements.
1244  * We do this by calculating their position based on the parent
1245  * (one level less refined) element, and the embedding matrix
1246  */
1247 
1248  const Elem* parent = elem->parent();
1249 
1250  /*
1251  * find out which child I am
1252  */
1253  for (unsigned int c=0; c < parent->n_children(); c++)
1254  {
1255  if (parent->child(c) == elem)
1256  {
1257  /*
1258  *loop over the childs (that is, the current elements) nodes
1259  */
1260  for (unsigned int nc=0; nc < elem->n_nodes(); nc++)
1261  {
1262  /*
1263  * the new position of the node
1264  */
1265  Point point;
1266  for (unsigned int n=0; n<parent->n_nodes(); n++)
1267  {
1268  /*
1269  * The value from the embedding matrix
1270  */
1271  const float em_val = parent->embedding_matrix(c,nc,n);
1272 
1273  if (em_val != 0.)
1274  point.add_scaled (parent->point(n), em_val);
1275  }
1276 
1277  const dof_id_type id = elem->get_node(nc)->id();
1278  new_positions[id] = point;
1279  weight[id] = 1.;
1280  }
1281 
1282  } // if parent->child == elem
1283  } // for parent->n_children
1284  } // if element refinement_level
1285 #endif // #ifdef LIBMESH_ENABLE_AMR
1286 
1287  } // element loop
1288 
1289  /*
1290  * finally reposition the vertex nodes
1291  */
1292  for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
1293  if (!on_boundary[nid] && weight[nid] > 0.)
1294  mesh.node(nid) = new_positions[nid]/weight[nid];
1295  }
1296 
1297  {
1298  /*
1299  * Now handle the additional second_order nodes by calculating
1300  * their position based on the vertex postitions
1301  * we do a second loop over the level elements
1302  */
1303  MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
1304  const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
1305 
1306  for (; el != end; ++el)
1307  {
1308  /*
1309  * Constant handle for the element
1310  */
1311  const Elem* elem = *el;
1312  const unsigned int son_begin = elem->n_vertices();
1313  const unsigned int son_end = elem->n_nodes();
1314  for (unsigned int n=son_begin; n<son_end; n++)
1315  {
1316  const unsigned int n_adjacent_vertices =
1318 
1319  Point point;
1320  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1321  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1322 
1323  const dof_id_type id = elem->get_node(n)->id();
1324  mesh.node(id) = point/n_adjacent_vertices;
1325  }
1326  }
1327  }
1328 
1329  } // refinement_level loop
1330 
1331  } // end iteration
1332 }
1333 
1334 
1335 
1336 #ifdef LIBMESH_ENABLE_AMR
1338 {
1339  // Algorithm:
1340  // .) For each active element in the mesh: construct a
1341  // copy which is the same in every way *except* it is
1342  // a level 0 element. Store the pointers to these in
1343  // a separate vector. Save any boundary information as well.
1344  // Delete the active element from the mesh.
1345  // .) Loop over all (remaining) elements in the mesh, delete them.
1346  // .) Add the level-0 copies back to the mesh
1347 
1348  // Temporary storage for new element pointers
1349  std::vector<Elem*> new_elements;
1350 
1351  // BoundaryInfo Storage for element ids, sides, and BC ids
1352  std::vector<Elem*> saved_boundary_elements;
1353  std::vector<boundary_id_type> saved_bc_ids;
1354  std::vector<unsigned short int> saved_bc_sides;
1355 
1356  // Reserve a reasonable amt. of space for each
1357  new_elements.reserve(mesh.n_active_elem());
1358  saved_boundary_elements.reserve(mesh.boundary_info->n_boundary_conds());
1359  saved_bc_ids.reserve(mesh.boundary_info->n_boundary_conds());
1360  saved_bc_sides.reserve(mesh.boundary_info->n_boundary_conds());
1361  {
1364 
1365  for (; it != end; ++it)
1366  {
1367  Elem* elem = *it;
1368 
1369  // Make a new element of the same type
1370  Elem* copy = Elem::build(elem->type()).release();
1371 
1372  // Set node pointers (they still point to nodes in the original mesh)
1373  for(unsigned int n=0; n<elem->n_nodes(); n++)
1374  copy->set_node(n) = elem->get_node(n);
1375 
1376  // Copy over ids
1377  copy->processor_id() = elem->processor_id();
1378  copy->subdomain_id() = elem->subdomain_id();
1379 
1380  // Retain the original element's ID as well, otherwise ParallelMesh will
1381  // try to create one for you...
1382  copy->set_id( elem->id() );
1383 
1384  // This element could have boundary info or ParallelMesh
1385  // remote_elem links as well. We need to save the (elem,
1386  // side, bc_id) triples and those links
1387  for (unsigned int s=0; s<elem->n_sides(); s++)
1388  {
1389  if (elem->neighbor(s) == remote_elem)
1390  copy->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
1391 
1392  const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids(elem,s);
1393  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1394  {
1395  const boundary_id_type bc_id = *id_it;
1396 
1397  if (bc_id != BoundaryInfo::invalid_id)
1398  {
1399  saved_boundary_elements.push_back(copy);
1400  saved_bc_ids.push_back(bc_id);
1401  saved_bc_sides.push_back(s);
1402  }
1403  }
1404  }
1405 
1406 
1407  // We're done with this element
1408  mesh.delete_elem(elem);
1409 
1410  // But save the copy
1411  new_elements.push_back(copy);
1412  }
1413 
1414  // Make sure we saved the same number of boundary conditions
1415  // in each vector.
1416  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1417  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1418  }
1419 
1420 
1421  // Loop again, delete any remaining elements
1422  {
1425 
1426  for (; it != end; ++it)
1427  mesh.delete_elem( *it );
1428  }
1429 
1430 
1431  // Add the copied (now level-0) elements back to the mesh
1432  {
1433  for (std::vector<Elem*>::iterator it = new_elements.begin();
1434  it != new_elements.end();
1435  ++it)
1436  {
1437  dof_id_type orig_id = (*it)->id();
1438 
1439  Elem* added_elem = mesh.add_elem(*it);
1440 
1441  dof_id_type added_id = added_elem->id();
1442 
1443  // If the Elem, as it was re-added to the mesh, now has a
1444  // different ID (this is unlikely, so it's just an assert)
1445  // the boundary information will no longer be correct.
1446  libmesh_assert_equal_to (orig_id, added_id);
1447  }
1448  }
1449 
1450  // Finally, also add back the saved boundary information
1451  for (unsigned int e=0; e<saved_boundary_elements.size(); ++e)
1452  mesh.boundary_info->add_side(saved_boundary_elements[e],
1453  saved_bc_sides[e],
1454  saved_bc_ids[e]);
1455 
1456  // Trim unused and renumber nodes and elements
1457  mesh.prepare_for_use(/*skip_renumber =*/ false);
1458 }
1459 #endif // #ifdef LIBMESH_ENABLE_AMR
1460 
1461 
1462 
1464  const boundary_id_type old_id,
1465  const boundary_id_type new_id)
1466 {
1467  // Only level-0 elements store BCs. Loop over them.
1469  const MeshBase::element_iterator end_el = mesh.level_elements_end(0);
1470  for (; el != end_el; ++el)
1471  {
1472  Elem *elem = *el;
1473 
1474  unsigned int n_nodes = elem->n_nodes();
1475  for (unsigned int n=0; n != n_nodes; ++n)
1476  {
1477  const std::vector<boundary_id_type>& old_ids = mesh.boundary_info->boundary_ids(elem->get_node(n));
1478  if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
1479  {
1480  std::vector<boundary_id_type> new_ids(old_ids);
1481  std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
1482  mesh.boundary_info->remove(elem->get_node(n));
1483  mesh.boundary_info->add_node(elem->get_node(n), new_ids);
1484  }
1485  }
1486 
1487  unsigned int n_edges = elem->n_edges();
1488  for (unsigned int edge=0; edge != n_edges; ++edge)
1489  {
1490  const std::vector<boundary_id_type>& old_ids = mesh.boundary_info->edge_boundary_ids(elem, edge);
1491  if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
1492  {
1493  std::vector<boundary_id_type> new_ids(old_ids);
1494  std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
1495  mesh.boundary_info->remove_edge(elem, edge);
1496  mesh.boundary_info->add_edge(elem, edge, new_ids);
1497  }
1498  }
1499 
1500  unsigned int n_sides = elem->n_sides();
1501  for (unsigned int s=0; s != n_sides; ++s)
1502  {
1503  const std::vector<boundary_id_type>& old_ids = mesh.boundary_info->boundary_ids(elem, s);
1504  if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
1505  {
1506  std::vector<boundary_id_type> new_ids(old_ids);
1507  std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
1508  mesh.boundary_info->remove_side(elem, s);
1509  mesh.boundary_info->add_side(elem, s, new_ids);
1510  }
1511  }
1512  }
1513 }
1514 
1515 
1516 
1518  const subdomain_id_type old_id,
1519  const subdomain_id_type new_id)
1520 {
1522  const MeshBase::element_iterator end_el = mesh.elements_end();
1523 
1524  for (; el != end_el; ++el)
1525  {
1526  Elem *elem = *el;
1527 
1528  if (elem->subdomain_id() == old_id)
1529  elem->subdomain_id() = new_id;
1530  }
1531 }
1532 
1533 
1534 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo