Functions
libMesh::MeshTools::Modification Namespace Reference

Functions

void distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
 
void redistribute (MeshBase &mesh, const FunctionBase< Real > &mapfunc)
 
void translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
 
void rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
 
void scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
 
void all_tri (MeshBase &mesh)
 
void smooth (MeshBase &, unsigned int, Real)
 
void flatten (MeshBase &mesh)
 
void change_boundary_id (MeshBase &mesh, const boundary_id_type old_id, const boundary_id_type new_id)
 
void change_subdomain_id (MeshBase &mesh, const subdomain_id_type old_id, const subdomain_id_type new_id)
 

Detailed Description

Tools for Mesh modification.

Author
Benjamin S. Kirk
Date
2004

Function Documentation

◆ all_tri()

void libMesh::MeshTools::Modification::all_tri ( MeshBase &  mesh)

Converts the 2D quadrilateral elements of a Mesh into triangular elements.

Note
Only works for 2D elements! 3D elements are ignored.
Probably won't do the right thing for meshes which have been refined previously.

Definition at line 271 of file mesh_modification.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build_side_ptr(), libMesh::ParallelObject::comm(), libMesh::MeshBase::delete_elem(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::MeshBase::element_ptr_range(), end, libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::INFEDGE2, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::DofObject::invalid_id, libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::is_serial(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::Parallel::Communicator::max(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::MeshBase::n_elem(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_id(), libMesh::Elem::node_ptr(), libMesh::MeshBase::parallel_max_unique_id(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::point(), libMesh::MeshBase::prepare_for_use(), libMesh::PRISM18, libMesh::PRISM6, libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::remote_elem, libMesh::BoundaryInfo::remove(), libMesh::Elem::set_node(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::TET10, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, libMesh::Elem::type(), and libMesh::DofObject::unique_id().

272 {
273  // The number of elements in the original mesh before any additions
274  // or deletions.
275  const dof_id_type n_orig_elem = mesh.n_elem();
276  const dof_id_type max_orig_id = mesh.max_elem_id();
277 
278  // We store pointers to the newly created elements in a vector
279  // until they are ready to be added to the mesh. This is because
280  // adding new elements on the fly can cause reallocation and invalidation
281  // of existing mesh element_iterators.
282  std::vector<Elem *> new_elements;
283 
284  unsigned int max_subelems = 1; // in 1D nothing needs to change
285  if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
286  max_subelems = 2;
287  if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
288  max_subelems = 6;
289 
290  new_elements.reserve (max_subelems*n_orig_elem);
291 
292  // If the original mesh has *side* boundary data, we carry that over
293  // to the new mesh with triangular elements. We currently only
294  // support bringing over side-based BCs to the all-tri mesh, but
295  // that could probably be extended to node and edge-based BCs as
296  // well.
297  const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
298 
299  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
300  std::vector<Elem *> new_bndry_elements;
301  std::vector<unsigned short int> new_bndry_sides;
302  std::vector<boundary_id_type> new_bndry_ids;
303 
304  // We may need to add new points if we run into a 1.5th order
305  // element; if we do that on a DistributedMesh in a ghost element then
306  // we will need to fix their ids / unique_ids
307  bool added_new_ghost_point = false;
308 
309  // Iterate over the elements, splitting:
310  // QUADs into pairs of conforming triangles
311  // PYRAMIDs into pairs of conforming tets,
312  // PRISMs into triplets of conforming tets, and
313  // HEXs into quintets or sextets of conforming tets.
314  // We split on the shortest diagonal to give us better
315  // triangle quality in 2D, and we split based on node ids
316  // to guarantee consistency in 3D.
317 
318  // FIXME: This algorithm does not work on refined grids!
319  {
320 #ifdef LIBMESH_ENABLE_UNIQUE_ID
321  unique_id_type max_unique_id = mesh.parallel_max_unique_id();
322 #endif
323 
324  for (auto & elem : mesh.element_ptr_range())
325  {
326  const ElemType etype = elem->type();
327 
328  // all_tri currently only works on coarse meshes
329  libmesh_assert (!elem->parent());
330 
331  // The new elements we will split the quad into.
332  // In 3D we may need as many as 6 tets per hex
333  Elem * subelem[6];
334 
335  for (unsigned int i = 0; i != max_subelems; ++i)
336  subelem[i] = nullptr;
337 
338  switch (etype)
339  {
340  case QUAD4:
341  {
342  subelem[0] = new Tri3;
343  subelem[1] = new Tri3;
344 
345  // Check for possible edge swap
346  if ((elem->point(0) - elem->point(2)).norm() <
347  (elem->point(1) - elem->point(3)).norm())
348  {
349  subelem[0]->set_node(0) = elem->node_ptr(0);
350  subelem[0]->set_node(1) = elem->node_ptr(1);
351  subelem[0]->set_node(2) = elem->node_ptr(2);
352 
353  subelem[1]->set_node(0) = elem->node_ptr(0);
354  subelem[1]->set_node(1) = elem->node_ptr(2);
355  subelem[1]->set_node(2) = elem->node_ptr(3);
356  }
357 
358  else
359  {
360  subelem[0]->set_node(0) = elem->node_ptr(0);
361  subelem[0]->set_node(1) = elem->node_ptr(1);
362  subelem[0]->set_node(2) = elem->node_ptr(3);
363 
364  subelem[1]->set_node(0) = elem->node_ptr(1);
365  subelem[1]->set_node(1) = elem->node_ptr(2);
366  subelem[1]->set_node(2) = elem->node_ptr(3);
367  }
368 
369 
370  break;
371  }
372 
373  case QUAD8:
374  {
375  if (elem->processor_id() != mesh.processor_id())
376  added_new_ghost_point = true;
377 
378  subelem[0] = new Tri6;
379  subelem[1] = new Tri6;
380 
381  // Add a new node at the center (vertex average) of the element.
382  Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
383  mesh.point(elem->node_id(1)) +
384  mesh.point(elem->node_id(2)) +
385  mesh.point(elem->node_id(3)))/4,
386  DofObject::invalid_id,
387  elem->processor_id());
388 
389  // Check for possible edge swap
390  if ((elem->point(0) - elem->point(2)).norm() <
391  (elem->point(1) - elem->point(3)).norm())
392  {
393  subelem[0]->set_node(0) = elem->node_ptr(0);
394  subelem[0]->set_node(1) = elem->node_ptr(1);
395  subelem[0]->set_node(2) = elem->node_ptr(2);
396  subelem[0]->set_node(3) = elem->node_ptr(4);
397  subelem[0]->set_node(4) = elem->node_ptr(5);
398  subelem[0]->set_node(5) = new_node;
399 
400  subelem[1]->set_node(0) = elem->node_ptr(0);
401  subelem[1]->set_node(1) = elem->node_ptr(2);
402  subelem[1]->set_node(2) = elem->node_ptr(3);
403  subelem[1]->set_node(3) = new_node;
404  subelem[1]->set_node(4) = elem->node_ptr(6);
405  subelem[1]->set_node(5) = elem->node_ptr(7);
406 
407  }
408 
409  else
410  {
411  subelem[0]->set_node(0) = elem->node_ptr(3);
412  subelem[0]->set_node(1) = elem->node_ptr(0);
413  subelem[0]->set_node(2) = elem->node_ptr(1);
414  subelem[0]->set_node(3) = elem->node_ptr(7);
415  subelem[0]->set_node(4) = elem->node_ptr(4);
416  subelem[0]->set_node(5) = new_node;
417 
418  subelem[1]->set_node(0) = elem->node_ptr(1);
419  subelem[1]->set_node(1) = elem->node_ptr(2);
420  subelem[1]->set_node(2) = elem->node_ptr(3);
421  subelem[1]->set_node(3) = elem->node_ptr(5);
422  subelem[1]->set_node(4) = elem->node_ptr(6);
423  subelem[1]->set_node(5) = new_node;
424  }
425 
426  break;
427  }
428 
429  case QUAD9:
430  {
431  subelem[0] = new Tri6;
432  subelem[1] = new Tri6;
433 
434  // Check for possible edge swap
435  if ((elem->point(0) - elem->point(2)).norm() <
436  (elem->point(1) - elem->point(3)).norm())
437  {
438  subelem[0]->set_node(0) = elem->node_ptr(0);
439  subelem[0]->set_node(1) = elem->node_ptr(1);
440  subelem[0]->set_node(2) = elem->node_ptr(2);
441  subelem[0]->set_node(3) = elem->node_ptr(4);
442  subelem[0]->set_node(4) = elem->node_ptr(5);
443  subelem[0]->set_node(5) = elem->node_ptr(8);
444 
445  subelem[1]->set_node(0) = elem->node_ptr(0);
446  subelem[1]->set_node(1) = elem->node_ptr(2);
447  subelem[1]->set_node(2) = elem->node_ptr(3);
448  subelem[1]->set_node(3) = elem->node_ptr(8);
449  subelem[1]->set_node(4) = elem->node_ptr(6);
450  subelem[1]->set_node(5) = elem->node_ptr(7);
451  }
452 
453  else
454  {
455  subelem[0]->set_node(0) = elem->node_ptr(0);
456  subelem[0]->set_node(1) = elem->node_ptr(1);
457  subelem[0]->set_node(2) = elem->node_ptr(3);
458  subelem[0]->set_node(3) = elem->node_ptr(4);
459  subelem[0]->set_node(4) = elem->node_ptr(8);
460  subelem[0]->set_node(5) = elem->node_ptr(7);
461 
462  subelem[1]->set_node(0) = elem->node_ptr(1);
463  subelem[1]->set_node(1) = elem->node_ptr(2);
464  subelem[1]->set_node(2) = elem->node_ptr(3);
465  subelem[1]->set_node(3) = elem->node_ptr(5);
466  subelem[1]->set_node(4) = elem->node_ptr(6);
467  subelem[1]->set_node(5) = elem->node_ptr(8);
468  }
469 
470  break;
471  }
472 
473  case PRISM6:
474  {
475  // Prisms all split into three tetrahedra
476  subelem[0] = new Tet4;
477  subelem[1] = new Tet4;
478  subelem[2] = new Tet4;
479 
480  // Triangular faces are not split.
481 
482  // On quad faces, we choose the node with the highest
483  // global id, and we split on the diagonal which
484  // includes that node. This ensures that (even in
485  // parallel, even on distributed meshes) the same
486  // diagonal split will be chosen for elements on either
487  // side of the same quad face. It also ensures that we
488  // always have a mix of "clockwise" and
489  // "counterclockwise" split faces (two of one and one
490  // of the other on each prism; this is useful since the
491  // alternative all-clockwise or all-counterclockwise
492  // face splittings can't be turned into tets without
493  // adding more nodes
494 
495  // Split on 0-4 diagonal
496  if (split_first_diagonal(elem, 0,4, 1,3))
497  {
498  // Split on 0-5 diagonal
499  if (split_first_diagonal(elem, 0,5, 2,3))
500  {
501  // Split on 1-5 diagonal
502  if (split_first_diagonal(elem, 1,5, 2,4))
503  {
504  subelem[0]->set_node(0) = elem->node_ptr(0);
505  subelem[0]->set_node(1) = elem->node_ptr(4);
506  subelem[0]->set_node(2) = elem->node_ptr(5);
507  subelem[0]->set_node(3) = elem->node_ptr(3);
508 
509  subelem[1]->set_node(0) = elem->node_ptr(0);
510  subelem[1]->set_node(1) = elem->node_ptr(4);
511  subelem[1]->set_node(2) = elem->node_ptr(1);
512  subelem[1]->set_node(3) = elem->node_ptr(5);
513 
514  subelem[2]->set_node(0) = elem->node_ptr(0);
515  subelem[2]->set_node(1) = elem->node_ptr(1);
516  subelem[2]->set_node(2) = elem->node_ptr(2);
517  subelem[2]->set_node(3) = elem->node_ptr(5);
518  }
519  else // Split on 2-4 diagonal
520  {
521  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
522 
523  subelem[0]->set_node(0) = elem->node_ptr(0);
524  subelem[0]->set_node(1) = elem->node_ptr(4);
525  subelem[0]->set_node(2) = elem->node_ptr(5);
526  subelem[0]->set_node(3) = elem->node_ptr(3);
527 
528  subelem[1]->set_node(0) = elem->node_ptr(0);
529  subelem[1]->set_node(1) = elem->node_ptr(4);
530  subelem[1]->set_node(2) = elem->node_ptr(2);
531  subelem[1]->set_node(3) = elem->node_ptr(5);
532 
533  subelem[2]->set_node(0) = elem->node_ptr(0);
534  subelem[2]->set_node(1) = elem->node_ptr(1);
535  subelem[2]->set_node(2) = elem->node_ptr(2);
536  subelem[2]->set_node(3) = elem->node_ptr(4);
537  }
538  }
539  else // Split on 2-3 diagonal
540  {
541  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
542 
543  // 0-4 and 2-3 split implies 2-4 split
544  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
545 
546  subelem[0]->set_node(0) = elem->node_ptr(0);
547  subelem[0]->set_node(1) = elem->node_ptr(4);
548  subelem[0]->set_node(2) = elem->node_ptr(2);
549  subelem[0]->set_node(3) = elem->node_ptr(3);
550 
551  subelem[1]->set_node(0) = elem->node_ptr(3);
552  subelem[1]->set_node(1) = elem->node_ptr(4);
553  subelem[1]->set_node(2) = elem->node_ptr(2);
554  subelem[1]->set_node(3) = elem->node_ptr(5);
555 
556  subelem[2]->set_node(0) = elem->node_ptr(0);
557  subelem[2]->set_node(1) = elem->node_ptr(1);
558  subelem[2]->set_node(2) = elem->node_ptr(2);
559  subelem[2]->set_node(3) = elem->node_ptr(4);
560  }
561  }
562  else // Split on 1-3 diagonal
563  {
564  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
565 
566  // Split on 0-5 diagonal
567  if (split_first_diagonal(elem, 0,5, 2,3))
568  {
569  // 1-3 and 0-5 split implies 1-5 split
570  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
571 
572  subelem[0]->set_node(0) = elem->node_ptr(1);
573  subelem[0]->set_node(1) = elem->node_ptr(3);
574  subelem[0]->set_node(2) = elem->node_ptr(4);
575  subelem[0]->set_node(3) = elem->node_ptr(5);
576 
577  subelem[1]->set_node(0) = elem->node_ptr(1);
578  subelem[1]->set_node(1) = elem->node_ptr(0);
579  subelem[1]->set_node(2) = elem->node_ptr(3);
580  subelem[1]->set_node(3) = elem->node_ptr(5);
581 
582  subelem[2]->set_node(0) = elem->node_ptr(0);
583  subelem[2]->set_node(1) = elem->node_ptr(1);
584  subelem[2]->set_node(2) = elem->node_ptr(2);
585  subelem[2]->set_node(3) = elem->node_ptr(5);
586  }
587  else // Split on 2-3 diagonal
588  {
589  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
590 
591  // Split on 1-5 diagonal
592  if (split_first_diagonal(elem, 1,5, 2,4))
593  {
594  subelem[0]->set_node(0) = elem->node_ptr(0);
595  subelem[0]->set_node(1) = elem->node_ptr(1);
596  subelem[0]->set_node(2) = elem->node_ptr(2);
597  subelem[0]->set_node(3) = elem->node_ptr(3);
598 
599  subelem[1]->set_node(0) = elem->node_ptr(3);
600  subelem[1]->set_node(1) = elem->node_ptr(1);
601  subelem[1]->set_node(2) = elem->node_ptr(2);
602  subelem[1]->set_node(3) = elem->node_ptr(5);
603 
604  subelem[2]->set_node(0) = elem->node_ptr(1);
605  subelem[2]->set_node(1) = elem->node_ptr(3);
606  subelem[2]->set_node(2) = elem->node_ptr(4);
607  subelem[2]->set_node(3) = elem->node_ptr(5);
608  }
609  else // Split on 2-4 diagonal
610  {
611  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
612 
613  subelem[0]->set_node(0) = elem->node_ptr(0);
614  subelem[0]->set_node(1) = elem->node_ptr(1);
615  subelem[0]->set_node(2) = elem->node_ptr(2);
616  subelem[0]->set_node(3) = elem->node_ptr(3);
617 
618  subelem[1]->set_node(0) = elem->node_ptr(2);
619  subelem[1]->set_node(1) = elem->node_ptr(3);
620  subelem[1]->set_node(2) = elem->node_ptr(4);
621  subelem[1]->set_node(3) = elem->node_ptr(5);
622 
623  subelem[2]->set_node(0) = elem->node_ptr(3);
624  subelem[2]->set_node(1) = elem->node_ptr(1);
625  subelem[2]->set_node(2) = elem->node_ptr(2);
626  subelem[2]->set_node(3) = elem->node_ptr(4);
627  }
628  }
629  }
630 
631  break;
632  }
633 
634  case PRISM18:
635  {
636  subelem[0] = new Tet10;
637  subelem[1] = new Tet10;
638  subelem[2] = new Tet10;
639 
640  // Split on 0-4 diagonal
641  if (split_first_diagonal(elem, 0,4, 1,3))
642  {
643  // Split on 0-5 diagonal
644  if (split_first_diagonal(elem, 0,5, 2,3))
645  {
646  // Split on 1-5 diagonal
647  if (split_first_diagonal(elem, 1,5, 2,4))
648  {
649  subelem[0]->set_node(0) = elem->node_ptr(0);
650  subelem[0]->set_node(1) = elem->node_ptr(4);
651  subelem[0]->set_node(2) = elem->node_ptr(5);
652  subelem[0]->set_node(3) = elem->node_ptr(3);
653 
654  subelem[0]->set_node(4) = elem->node_ptr(15);
655  subelem[0]->set_node(5) = elem->node_ptr(13);
656  subelem[0]->set_node(6) = elem->node_ptr(17);
657  subelem[0]->set_node(7) = elem->node_ptr(9);
658  subelem[0]->set_node(8) = elem->node_ptr(12);
659  subelem[0]->set_node(9) = elem->node_ptr(14);
660 
661  subelem[1]->set_node(0) = elem->node_ptr(0);
662  subelem[1]->set_node(1) = elem->node_ptr(4);
663  subelem[1]->set_node(2) = elem->node_ptr(1);
664  subelem[1]->set_node(3) = elem->node_ptr(5);
665 
666  subelem[1]->set_node(4) = elem->node_ptr(15);
667  subelem[1]->set_node(5) = elem->node_ptr(10);
668  subelem[1]->set_node(6) = elem->node_ptr(6);
669  subelem[1]->set_node(7) = elem->node_ptr(17);
670  subelem[1]->set_node(8) = elem->node_ptr(13);
671  subelem[1]->set_node(9) = elem->node_ptr(16);
672 
673  subelem[2]->set_node(0) = elem->node_ptr(0);
674  subelem[2]->set_node(1) = elem->node_ptr(1);
675  subelem[2]->set_node(2) = elem->node_ptr(2);
676  subelem[2]->set_node(3) = elem->node_ptr(5);
677 
678  subelem[2]->set_node(4) = elem->node_ptr(6);
679  subelem[2]->set_node(5) = elem->node_ptr(7);
680  subelem[2]->set_node(6) = elem->node_ptr(8);
681  subelem[2]->set_node(7) = elem->node_ptr(17);
682  subelem[2]->set_node(8) = elem->node_ptr(16);
683  subelem[2]->set_node(9) = elem->node_ptr(11);
684  }
685  else // Split on 2-4 diagonal
686  {
687  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
688 
689  subelem[0]->set_node(0) = elem->node_ptr(0);
690  subelem[0]->set_node(1) = elem->node_ptr(4);
691  subelem[0]->set_node(2) = elem->node_ptr(5);
692  subelem[0]->set_node(3) = elem->node_ptr(3);
693 
694  subelem[0]->set_node(4) = elem->node_ptr(15);
695  subelem[0]->set_node(5) = elem->node_ptr(13);
696  subelem[0]->set_node(6) = elem->node_ptr(17);
697  subelem[0]->set_node(7) = elem->node_ptr(9);
698  subelem[0]->set_node(8) = elem->node_ptr(12);
699  subelem[0]->set_node(9) = elem->node_ptr(14);
700 
701  subelem[1]->set_node(0) = elem->node_ptr(0);
702  subelem[1]->set_node(1) = elem->node_ptr(4);
703  subelem[1]->set_node(2) = elem->node_ptr(2);
704  subelem[1]->set_node(3) = elem->node_ptr(5);
705 
706  subelem[1]->set_node(4) = elem->node_ptr(15);
707  subelem[1]->set_node(5) = elem->node_ptr(16);
708  subelem[1]->set_node(6) = elem->node_ptr(8);
709  subelem[1]->set_node(7) = elem->node_ptr(17);
710  subelem[1]->set_node(8) = elem->node_ptr(13);
711  subelem[1]->set_node(9) = elem->node_ptr(11);
712 
713  subelem[2]->set_node(0) = elem->node_ptr(0);
714  subelem[2]->set_node(1) = elem->node_ptr(1);
715  subelem[2]->set_node(2) = elem->node_ptr(2);
716  subelem[2]->set_node(3) = elem->node_ptr(4);
717 
718  subelem[2]->set_node(4) = elem->node_ptr(6);
719  subelem[2]->set_node(5) = elem->node_ptr(7);
720  subelem[2]->set_node(6) = elem->node_ptr(8);
721  subelem[2]->set_node(7) = elem->node_ptr(15);
722  subelem[2]->set_node(8) = elem->node_ptr(10);
723  subelem[2]->set_node(9) = elem->node_ptr(16);
724  }
725  }
726  else // Split on 2-3 diagonal
727  {
728  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
729 
730  // 0-4 and 2-3 split implies 2-4 split
731  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
732 
733  subelem[0]->set_node(0) = elem->node_ptr(0);
734  subelem[0]->set_node(1) = elem->node_ptr(4);
735  subelem[0]->set_node(2) = elem->node_ptr(2);
736  subelem[0]->set_node(3) = elem->node_ptr(3);
737 
738  subelem[0]->set_node(4) = elem->node_ptr(15);
739  subelem[0]->set_node(5) = elem->node_ptr(16);
740  subelem[0]->set_node(6) = elem->node_ptr(8);
741  subelem[0]->set_node(7) = elem->node_ptr(9);
742  subelem[0]->set_node(8) = elem->node_ptr(12);
743  subelem[0]->set_node(9) = elem->node_ptr(17);
744 
745  subelem[1]->set_node(0) = elem->node_ptr(3);
746  subelem[1]->set_node(1) = elem->node_ptr(4);
747  subelem[1]->set_node(2) = elem->node_ptr(2);
748  subelem[1]->set_node(3) = elem->node_ptr(5);
749 
750  subelem[1]->set_node(4) = elem->node_ptr(12);
751  subelem[1]->set_node(5) = elem->node_ptr(16);
752  subelem[1]->set_node(6) = elem->node_ptr(17);
753  subelem[1]->set_node(7) = elem->node_ptr(14);
754  subelem[1]->set_node(8) = elem->node_ptr(13);
755  subelem[1]->set_node(9) = elem->node_ptr(11);
756 
757  subelem[2]->set_node(0) = elem->node_ptr(0);
758  subelem[2]->set_node(1) = elem->node_ptr(1);
759  subelem[2]->set_node(2) = elem->node_ptr(2);
760  subelem[2]->set_node(3) = elem->node_ptr(4);
761 
762  subelem[2]->set_node(4) = elem->node_ptr(6);
763  subelem[2]->set_node(5) = elem->node_ptr(7);
764  subelem[2]->set_node(6) = elem->node_ptr(8);
765  subelem[2]->set_node(7) = elem->node_ptr(15);
766  subelem[2]->set_node(8) = elem->node_ptr(10);
767  subelem[2]->set_node(9) = elem->node_ptr(16);
768  }
769  }
770  else // Split on 1-3 diagonal
771  {
772  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
773 
774  // Split on 0-5 diagonal
775  if (split_first_diagonal(elem, 0,5, 2,3))
776  {
777  // 1-3 and 0-5 split implies 1-5 split
778  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
779 
780  subelem[0]->set_node(0) = elem->node_ptr(1);
781  subelem[0]->set_node(1) = elem->node_ptr(3);
782  subelem[0]->set_node(2) = elem->node_ptr(4);
783  subelem[0]->set_node(3) = elem->node_ptr(5);
784 
785  subelem[0]->set_node(4) = elem->node_ptr(15);
786  subelem[0]->set_node(5) = elem->node_ptr(12);
787  subelem[0]->set_node(6) = elem->node_ptr(10);
788  subelem[0]->set_node(7) = elem->node_ptr(16);
789  subelem[0]->set_node(8) = elem->node_ptr(14);
790  subelem[0]->set_node(9) = elem->node_ptr(13);
791 
792  subelem[1]->set_node(0) = elem->node_ptr(1);
793  subelem[1]->set_node(1) = elem->node_ptr(0);
794  subelem[1]->set_node(2) = elem->node_ptr(3);
795  subelem[1]->set_node(3) = elem->node_ptr(5);
796 
797  subelem[1]->set_node(4) = elem->node_ptr(6);
798  subelem[1]->set_node(5) = elem->node_ptr(9);
799  subelem[1]->set_node(6) = elem->node_ptr(15);
800  subelem[1]->set_node(7) = elem->node_ptr(16);
801  subelem[1]->set_node(8) = elem->node_ptr(17);
802  subelem[1]->set_node(9) = elem->node_ptr(14);
803 
804  subelem[2]->set_node(0) = elem->node_ptr(0);
805  subelem[2]->set_node(1) = elem->node_ptr(1);
806  subelem[2]->set_node(2) = elem->node_ptr(2);
807  subelem[2]->set_node(3) = elem->node_ptr(5);
808 
809  subelem[2]->set_node(4) = elem->node_ptr(6);
810  subelem[2]->set_node(5) = elem->node_ptr(7);
811  subelem[2]->set_node(6) = elem->node_ptr(8);
812  subelem[2]->set_node(7) = elem->node_ptr(17);
813  subelem[2]->set_node(8) = elem->node_ptr(16);
814  subelem[2]->set_node(9) = elem->node_ptr(11);
815  }
816  else // Split on 2-3 diagonal
817  {
818  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
819 
820  // Split on 1-5 diagonal
821  if (split_first_diagonal(elem, 1,5, 2,4))
822  {
823  subelem[0]->set_node(0) = elem->node_ptr(0);
824  subelem[0]->set_node(1) = elem->node_ptr(1);
825  subelem[0]->set_node(2) = elem->node_ptr(2);
826  subelem[0]->set_node(3) = elem->node_ptr(3);
827 
828  subelem[0]->set_node(4) = elem->node_ptr(6);
829  subelem[0]->set_node(5) = elem->node_ptr(7);
830  subelem[0]->set_node(6) = elem->node_ptr(8);
831  subelem[0]->set_node(7) = elem->node_ptr(9);
832  subelem[0]->set_node(8) = elem->node_ptr(15);
833  subelem[0]->set_node(9) = elem->node_ptr(17);
834 
835  subelem[1]->set_node(0) = elem->node_ptr(3);
836  subelem[1]->set_node(1) = elem->node_ptr(1);
837  subelem[1]->set_node(2) = elem->node_ptr(2);
838  subelem[1]->set_node(3) = elem->node_ptr(5);
839 
840  subelem[1]->set_node(4) = elem->node_ptr(15);
841  subelem[1]->set_node(5) = elem->node_ptr(7);
842  subelem[1]->set_node(6) = elem->node_ptr(17);
843  subelem[1]->set_node(7) = elem->node_ptr(14);
844  subelem[1]->set_node(8) = elem->node_ptr(16);
845  subelem[1]->set_node(9) = elem->node_ptr(11);
846 
847  subelem[2]->set_node(0) = elem->node_ptr(1);
848  subelem[2]->set_node(1) = elem->node_ptr(3);
849  subelem[2]->set_node(2) = elem->node_ptr(4);
850  subelem[2]->set_node(3) = elem->node_ptr(5);
851 
852  subelem[2]->set_node(4) = elem->node_ptr(15);
853  subelem[2]->set_node(5) = elem->node_ptr(12);
854  subelem[2]->set_node(6) = elem->node_ptr(10);
855  subelem[2]->set_node(7) = elem->node_ptr(16);
856  subelem[2]->set_node(8) = elem->node_ptr(14);
857  subelem[2]->set_node(9) = elem->node_ptr(13);
858  }
859  else // Split on 2-4 diagonal
860  {
861  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
862 
863  subelem[0]->set_node(0) = elem->node_ptr(0);
864  subelem[0]->set_node(1) = elem->node_ptr(1);
865  subelem[0]->set_node(2) = elem->node_ptr(2);
866  subelem[0]->set_node(3) = elem->node_ptr(3);
867 
868  subelem[0]->set_node(4) = elem->node_ptr(6);
869  subelem[0]->set_node(5) = elem->node_ptr(7);
870  subelem[0]->set_node(6) = elem->node_ptr(8);
871  subelem[0]->set_node(7) = elem->node_ptr(9);
872  subelem[0]->set_node(8) = elem->node_ptr(15);
873  subelem[0]->set_node(9) = elem->node_ptr(17);
874 
875  subelem[1]->set_node(0) = elem->node_ptr(2);
876  subelem[1]->set_node(1) = elem->node_ptr(3);
877  subelem[1]->set_node(2) = elem->node_ptr(4);
878  subelem[1]->set_node(3) = elem->node_ptr(5);
879 
880  subelem[1]->set_node(4) = elem->node_ptr(17);
881  subelem[1]->set_node(5) = elem->node_ptr(12);
882  subelem[1]->set_node(6) = elem->node_ptr(16);
883  subelem[1]->set_node(7) = elem->node_ptr(11);
884  subelem[1]->set_node(8) = elem->node_ptr(14);
885  subelem[1]->set_node(9) = elem->node_ptr(13);
886 
887  subelem[2]->set_node(0) = elem->node_ptr(3);
888  subelem[2]->set_node(1) = elem->node_ptr(1);
889  subelem[2]->set_node(2) = elem->node_ptr(2);
890  subelem[2]->set_node(3) = elem->node_ptr(4);
891 
892  subelem[2]->set_node(4) = elem->node_ptr(15);
893  subelem[2]->set_node(5) = elem->node_ptr(7);
894  subelem[2]->set_node(6) = elem->node_ptr(17);
895  subelem[2]->set_node(7) = elem->node_ptr(12);
896  subelem[2]->set_node(8) = elem->node_ptr(10);
897  subelem[2]->set_node(9) = elem->node_ptr(16);
898  }
899  }
900  }
901 
902  break;
903  }
904 
905  // No need to split elements that are already simplicial:
906  case EDGE2:
907  case EDGE3:
908  case EDGE4:
909  case TRI3:
910  case TRI6:
911  case TET4:
912  case TET10:
913  case INFEDGE2:
914  // No way to split infinite quad/prism elements, so
915  // hopefully no need to
916  case INFQUAD4:
917  case INFQUAD6:
918  case INFPRISM6:
919  case INFPRISM12:
920  continue;
921  // If we're left with an unimplemented hex we're probably
922  // out of luck. TODO: implement hexes
923  default:
924  {
925  libMesh::err << "Error, encountered unimplemented element "
926  << Utility::enum_to_string<ElemType>(etype)
927  << " in MeshTools::Modification::all_tri()..."
928  << std::endl;
929  libmesh_not_implemented();
930  }
931  } // end switch (etype)
932 
933 
934 
935  // Be sure the correct IDs are also set for all subelems.
936  for (unsigned int i=0; i != max_subelems; ++i)
937  if (subelem[i]) {
938  subelem[i]->processor_id() = elem->processor_id();
939  subelem[i]->subdomain_id() = elem->subdomain_id();
940  }
941 
942  // On a mesh with boundary data, we need to move that data to
943  // the new elements.
944 
945  // On a mesh which is distributed, we need to move
946  // remote_elem links to the new elements.
947  bool mesh_is_serial = mesh.is_serial();
948 
949  if (mesh_has_boundary_data || mesh_is_serial)
950  {
951  // Container to key boundary IDs handed back by the BoundaryInfo object.
952  std::vector<boundary_id_type> bc_ids;
953 
954  for (auto sn : elem->side_index_range())
955  {
956  mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
957  for (const auto & b_id : bc_ids)
958  {
959  if (mesh_is_serial && b_id == BoundaryInfo::invalid_id)
960  continue;
961 
962  // Make a sorted list of node ids for elem->side(sn)
963  std::unique_ptr<Elem> elem_side = elem->build_side_ptr(sn);
964  std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
965  for (unsigned int esn=0,
966  n_esn = cast_int<unsigned int>(elem_side_nodes.size());
967  esn != n_esn; ++esn)
968  elem_side_nodes[esn] = elem_side->node_id(esn);
969  std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
970 
971  for (unsigned int i=0; i != max_subelems; ++i)
972  if (subelem[i])
973  {
974  for (auto subside : subelem[i]->side_index_range())
975  {
976  std::unique_ptr<Elem> subside_elem = subelem[i]->build_side_ptr(subside);
977 
978  // Make a list of *vertex* node ids for this subside, see if they are all present
979  // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
980  // subelem[i]->key(subside) in the Prism cases, since the new side is
981  // a different type. Note 2: we only use vertex nodes since, in the future,
982  // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
983  // original face will not contain the mid-edge node.
984  std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
985  for (unsigned int ssn=0,
986  n_ssn = cast_int<unsigned int>(subside_nodes.size());
987  ssn != n_ssn; ++ssn)
988  subside_nodes[ssn] = subside_elem->node_id(ssn);
989  std::sort(subside_nodes.begin(), subside_nodes.end());
990 
991  // std::includes returns true if every element of the second sorted range is
992  // contained in the first sorted range.
993  if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
994  subside_nodes.begin(), subside_nodes.end()))
995  {
996  if (b_id != BoundaryInfo::invalid_id)
997  {
998  new_bndry_ids.push_back(b_id);
999  new_bndry_elements.push_back(subelem[i]);
1000  new_bndry_sides.push_back(subside);
1001  }
1002 
1003  // If the original element had a RemoteElem neighbor on side 'sn',
1004  // then the subelem has one on side 'subside'.
1005  if (elem->neighbor_ptr(sn) == remote_elem)
1006  subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1007  }
1008  }
1009  }
1010  } // end for loop over boundary IDs
1011  } // end for loop over sides
1012 
1013  // Remove the original element from the BoundaryInfo structure.
1014  mesh.get_boundary_info().remove(elem);
1015 
1016  } // end if (mesh_has_boundary_data)
1017 
1018  // Determine new IDs for the split elements which will be
1019  // the same on all processors, therefore keeping the Mesh
1020  // in sync. Note: we offset the new IDs by max_orig_id to
1021  // avoid overwriting any of the original IDs.
1022  for (unsigned int i=0; i != max_subelems; ++i)
1023  if (subelem[i])
1024  {
1025  // Determine new IDs for the split elements which will be
1026  // the same on all processors, therefore keeping the Mesh
1027  // in sync. Note: we offset the new IDs by the max of the
1028  // pre-existing ids to avoid conflicting with originals.
1029  subelem[i]->set_id( max_orig_id + 6*elem->id() + i );
1030 
1031 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1032  subelem[i]->set_unique_id() = max_unique_id + 2*elem->unique_id() + i;
1033 #endif
1034 
1035  // Prepare to add the newly-created simplices
1036  new_elements.push_back(subelem[i]);
1037  }
1038 
1039  // Delete the original element
1040  mesh.delete_elem(elem);
1041  } // End for loop over elements
1042  } // end scope
1043 
1044 
1045  // Now, iterate over the new elements vector, and add them each to
1046  // the Mesh.
1047  {
1048  std::vector<Elem *>::iterator el = new_elements.begin();
1049  const std::vector<Elem *>::iterator end = new_elements.end();
1050  for (; el != end; ++el)
1051  mesh.add_elem(*el);
1052  }
1053 
1054  if (mesh_has_boundary_data)
1055  {
1056  // If the old mesh had boundary data, the new mesh better have
1057  // some. However, we can't assert that the size of
1058  // new_bndry_elements vector is > 0, since we may not have split
1059  // any elements actually on the boundary. We also can't assert
1060  // that the original number of boundary sides is equal to the
1061  // sum of the boundary sides currently in the mesh and the
1062  // newly-added boundary sides, since in 3D, we may have split a
1063  // boundary QUAD into two boundary TRIs. Therefore, we won't be
1064  // too picky about the actual number of BCs, and just assert that
1065  // there are some, somewhere.
1066 #ifndef NDEBUG
1067  bool nbe_nonempty = new_bndry_elements.size();
1068  mesh.comm().max(nbe_nonempty);
1069  libmesh_assert(nbe_nonempty ||
1070  mesh.get_boundary_info().n_boundary_conds()>0);
1071 #endif
1072 
1073  // We should also be sure that the lengths of the new boundary data vectors
1074  // are all the same.
1075  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1076  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1077 
1078  // Add the new boundary info to the mesh
1079  for (std::size_t s=0; s<new_bndry_elements.size(); ++s)
1080  mesh.get_boundary_info().add_side(new_bndry_elements[s],
1081  new_bndry_sides[s],
1082  new_bndry_ids[s]);
1083  }
1084 
1085  // In a DistributedMesh any newly added ghost node ids may be
1086  // inconsistent, and unique_ids of newly added ghost nodes remain
1087  // unset.
1088  // make_nodes_parallel_consistent() will fix all this.
1089  if (!mesh.is_serial())
1090  {
1091  mesh.comm().max(added_new_ghost_point);
1092 
1093  if (added_new_ghost_point)
1094  MeshCommunication().make_nodes_parallel_consistent (mesh);
1095  }
1096 
1097 
1098 
1099  // Prepare the newly created mesh for use.
1100  mesh.prepare_for_use(/*skip_renumber =*/ false);
1101 
1102  // Let the new_elements and new_bndry_elements vectors go out of scope.
1103 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
end
IterBase * end
Definition: variant_filter_iterator.h:343
libMesh::err
OStreamProxy err(std::cerr)
Definition: libmesh_common.h:229
libMesh::unique_id_type
uint8_t unique_id_type
Definition: id_types.h:79
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:64
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57

◆ change_boundary_id()

void libMesh::MeshTools::Modification::change_boundary_id ( MeshBase &  mesh,
const boundary_id_type  old_id,
const boundary_id_type  new_id 
)

Finds any boundary ids that are currently old_id, changes them to new_id

Definition at line 1374 of file mesh_modification.C.

References libMesh::BoundaryInfo::add_edge(), libMesh::BoundaryInfo::add_node(), libMesh::BoundaryInfo::add_shellface(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::BoundaryInfo::build_edge_list(), libMesh::BoundaryInfo::build_node_list(), libMesh::BoundaryInfo::build_shellface_list(), libMesh::BoundaryInfo::build_side_list(), libMesh::BoundaryInfo::edge_boundary_ids(), libMesh::MeshBase::elem_ptr(), libMesh::MeshBase::get_boundary_info(), mesh, libMesh::MeshBase::node_ptr(), libMesh::BoundaryInfo::remove(), libMesh::BoundaryInfo::remove_edge(), libMesh::BoundaryInfo::remove_id(), libMesh::BoundaryInfo::remove_shellface(), libMesh::BoundaryInfo::remove_side(), libMesh::BoundaryInfo::shellface_boundary_ids(), and side.

1377 {
1378  if (old_id == new_id)
1379  {
1380  // If the IDs are the same, this is a no-op.
1381  return;
1382  }
1383 
1384  // A reference to the Mesh's BoundaryInfo object, for convenience.
1385  BoundaryInfo & bi = mesh.get_boundary_info();
1386 
1387  {
1388  // Temporary vector to hold ids
1389  std::vector<boundary_id_type> bndry_ids;
1390 
1391  // build_node_list returns a vector of (node, bc) tuples.
1392  for (const auto & t : bi.build_node_list())
1393  if (std::get<1>(t) == old_id)
1394  {
1395  // Get the node in question
1396  const Node * node = mesh.node_ptr(std::get<0>(t));
1397 
1398  // Get all the current IDs for this node.
1399  bi.boundary_ids(node, bndry_ids);
1400 
1401  // Update the IDs accordingly
1402  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1403 
1404  // Remove all traces of that node from the BoundaryInfo object
1405  bi.remove(node);
1406 
1407  // Add it back with the updated IDs
1408  bi.add_node(node, bndry_ids);
1409  }
1410  }
1411 
1412  {
1413  // Temporary vector to hold ids
1414  std::vector<boundary_id_type> bndry_ids;
1415 
1416  // build_edge_list returns a vector of (elem, side, bc) tuples.
1417  for (const auto & t : bi.build_edge_list())
1418  if (std::get<2>(t) == old_id)
1419  {
1420  // Get the elem in question
1421  const Elem * elem = mesh.elem_ptr(std::get<0>(t));
1422 
1423  // The edge of the elem in question
1424  unsigned short int edge = std::get<1>(t);
1425 
1426  // Get all the current IDs for the edge in question.
1427  bi.edge_boundary_ids(elem, edge, bndry_ids);
1428 
1429  // Update the IDs accordingly
1430  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1431 
1432  // Remove all traces of that edge from the BoundaryInfo object
1433  bi.remove_edge(elem, edge);
1434 
1435  // Add it back with the updated IDs
1436  bi.add_edge(elem, edge, bndry_ids);
1437  }
1438  }
1439 
1440  {
1441  // Temporary vector to hold ids
1442  std::vector<boundary_id_type> bndry_ids;
1443 
1444  // build_shellface_list returns a vector of (elem, side, bc) tuples.
1445  for (const auto & t : bi.build_shellface_list())
1446  if (std::get<2>(t) == old_id)
1447  {
1448  // Get the elem in question
1449  const Elem * elem = mesh.elem_ptr(std::get<0>(t));
1450 
1451  // The shellface of the elem in question
1452  unsigned short int shellface = std::get<1>(t);
1453 
1454  // Get all the current IDs for the shellface in question.
1455  bi.shellface_boundary_ids(elem, shellface, bndry_ids);
1456 
1457  // Update the IDs accordingly
1458  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1459 
1460  // Remove all traces of that shellface from the BoundaryInfo object
1461  bi.remove_shellface(elem, shellface);
1462 
1463  // Add it back with the updated IDs
1464  bi.add_shellface(elem, shellface, bndry_ids);
1465  }
1466  }
1467 
1468  {
1469  // Temporary vector to hold ids
1470  std::vector<boundary_id_type> bndry_ids;
1471 
1472  // build_side_list returns a vector of (elem, side, bc) tuples.
1473  for (const auto & t : bi.build_side_list())
1474  if (std::get<2>(t) == old_id)
1475  {
1476  // Get the elem in question
1477  const Elem * elem = mesh.elem_ptr(std::get<0>(t));
1478 
1479  // The side of the elem in question
1480  unsigned short int side = std::get<1>(t);
1481 
1482  // Get all the current IDs for the side in question.
1483  bi.boundary_ids(elem, side, bndry_ids);
1484 
1485  // Update the IDs accordingly
1486  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1487 
1488  // Remove all traces of that side from the BoundaryInfo object
1489  bi.remove_side(elem, side);
1490 
1491  // Add it back with the updated IDs
1492  bi.add_side(elem, side, bndry_ids);
1493  }
1494  }
1495 
1496  // Remove any remaining references to the old_id so it does not show
1497  // up in lists of boundary ids, etc.
1498  bi.remove_id(old_id);
1499 }
side
unsigned short int side
Definition: xdr_io.C:50
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245

◆ change_subdomain_id()

void libMesh::MeshTools::Modification::change_subdomain_id ( MeshBase &  mesh,
const subdomain_id_type  old_id,
const subdomain_id_type  new_id 
)

Finds any subdomain ids that are currently old_id, changes them to new_id

Definition at line 1503 of file mesh_modification.C.

References libMesh::MeshBase::element_ptr_range(), mesh, and libMesh::Elem::subdomain_id().

1506 {
1507  if (old_id == new_id)
1508  {
1509  // If the IDs are the same, this is a no-op.
1510  return;
1511  }
1512 
1513  for (auto & elem : mesh.element_ptr_range())
1514  {
1515  if (elem->subdomain_id() == old_id)
1516  elem->subdomain_id() = new_id;
1517  }
1518 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245

◆ distort()

void libMesh::MeshTools::Modification::distort ( MeshBase &  mesh,
const Real  factor,
const bool  perturb_boundary = false 
)

Randomly perturb the nodal locations. This function will move each node factor fraction of its minimum neighboring node separation distance. Nodes on the boundary are not moved by default, however they may be by setting the flag perturb_boundary true.

Definition at line 65 of file mesh_modification.C.

References libMesh::MeshBase::active_element_ptr_range(), libMesh::MeshTools::find_boundary_nodes(), libMesh::Elem::hmin(), std::max(), libMesh::MeshBase::max_node_id(), mesh, libMesh::MeshBase::mesh_dimension(), std::min(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::Elem::node_ref_range(), libMesh::MeshBase::query_node_ptr(), and libMesh::TypeVector< T >::unit().

68 {
69  libmesh_assert (mesh.n_nodes());
70  libmesh_assert (mesh.n_elem());
71  libmesh_assert ((factor >= 0.) && (factor <= 1.));
72 
73  LOG_SCOPE("distort()", "MeshTools::Modification");
74 
75  // If we are not perturbing boundary nodes, make a
76  // quickly-searchable list of node ids we can check against.
77  std::unordered_set<dof_id_type> boundary_node_ids;
78  if (!perturb_boundary)
79  boundary_node_ids = MeshTools::find_boundary_nodes (mesh);
80 
81  // Now calculate the minimum distance to
82  // neighboring nodes for each node.
83  // hmin holds these distances.
84  std::vector<float> hmin (mesh.max_node_id(),
85  std::numeric_limits<float>::max());
86 
87  for (const auto & elem : mesh.active_element_ptr_range())
88  for (auto & n : elem->node_ref_range())
89  hmin[n.id()] = std::min(hmin[n.id()],
90  static_cast<float>(elem->hmin()));
91 
92  // Now actually move the nodes
93  {
94  const unsigned int seed = 123456;
95 
96  // seed the random number generator.
97  // We'll loop from 1 to n_nodes on every processor, even those
98  // that don't have a particular node, so that the pseudorandom
99  // numbers will be the same everywhere.
100  std::srand(seed);
101 
102  // If the node is on the boundary or
103  // the node is not used by any element (hmin[n]<1.e20)
104  // then we should not move it.
105  // [Note: Testing for (in)equality might be wrong
106  // (different types, namely float and double)]
107  for (unsigned int n=0; n<mesh.max_node_id(); n++)
108  if ((perturb_boundary || !boundary_node_ids.count(n)) && hmin[n] < 1.e20)
109  {
110  // the direction, random but unit normalized
111  Point dir (static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
112  (mesh.mesh_dimension() > 1) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.,
113  ((mesh.mesh_dimension() == 3) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.));
114 
115  dir(0) = (dir(0)-.5)*2.;
116  if (mesh.mesh_dimension() > 1)
117  dir(1) = (dir(1)-.5)*2.;
118  if (mesh.mesh_dimension() == 3)
119  dir(2) = (dir(2)-.5)*2.;
120 
121  dir = dir.unit();
122 
123  Node * node = mesh.query_node_ptr(n);
124  if (!node)
125  continue;
126 
127  (*node)(0) += dir(0)*factor*hmin[n];
128  if (mesh.mesh_dimension() > 1)
129  (*node)(1) += dir(1)*factor*hmin[n];
130  if (mesh.mesh_dimension() == 3)
131  (*node)(2) += dir(2)*factor*hmin[n];
132  }
133  }
134 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
std::max
long double max(long double a, double b)
Definition: libmesh_augment_std_namespace.h:27
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
std::min
long double min(long double a, double b)
Definition: libmesh_augment_std_namespace.h:29
libMesh::MeshTools::find_boundary_nodes
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Definition: mesh_tools.C:303

◆ flatten()

void libMesh::MeshTools::Modification::flatten ( MeshBase &  mesh)

Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elements. This is useful when you want to write out a uniformly-refined grid to be treated later as an initial mesh.

Note
Many functions in LibMesh assume a conforming (with no hanging nodes) grid exists at some level, so you probably only want to do this on meshes which have been uniformly refined.

Definition at line 1260 of file mesh_modification.C.

References libMesh::MeshBase::active_element_ptr_range(), libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_side(), bc_id, libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::element_ptr_range(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::BoundaryInfo::invalid_id, libMesh::libmesh_ignore(), mesh, libMesh::MeshBase::n_active_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_index_range(), libMesh::Elem::node_ptr(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMesh::remote_elem, libMesh::DofObject::set_id(), libMesh::Elem::set_neighbor(), libMesh::Elem::set_node(), libMesh::DofObject::set_unique_id(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::Elem::type(), and libMesh::DofObject::unique_id().

1261 {
1262  // Algorithm:
1263  // .) For each active element in the mesh: construct a
1264  // copy which is the same in every way *except* it is
1265  // a level 0 element. Store the pointers to these in
1266  // a separate vector. Save any boundary information as well.
1267  // Delete the active element from the mesh.
1268  // .) Loop over all (remaining) elements in the mesh, delete them.
1269  // .) Add the level-0 copies back to the mesh
1270 
1271  // Temporary storage for new element pointers
1272  std::vector<Elem *> new_elements;
1273 
1274  // BoundaryInfo Storage for element ids, sides, and BC ids
1275  std::vector<Elem *> saved_boundary_elements;
1276  std::vector<boundary_id_type> saved_bc_ids;
1277  std::vector<unsigned short int> saved_bc_sides;
1278 
1279  // Container to catch boundary ids passed back by BoundaryInfo
1280  std::vector<boundary_id_type> bc_ids;
1281 
1282  // Reserve a reasonable amt. of space for each
1283  new_elements.reserve(mesh.n_active_elem());
1284  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1285  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1286  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1287 
1288  for (auto & elem : mesh.active_element_ptr_range())
1289  {
1290  // Make a new element of the same type
1291  Elem * copy = Elem::build(elem->type()).release();
1292 
1293  // Set node pointers (they still point to nodes in the original mesh)
1294  for (auto n : elem->node_index_range())
1295  copy->set_node(n) = elem->node_ptr(n);
1296 
1297  // Copy over ids
1298  copy->processor_id() = elem->processor_id();
1299  copy->subdomain_id() = elem->subdomain_id();
1300 
1301  // Retain the original element's ID(s) as well, otherwise
1302  // the Mesh may try to create them for you...
1303  copy->set_id( elem->id() );
1304 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1305  copy->set_unique_id() = elem->unique_id();
1306 #endif
1307 
1308  // This element could have boundary info or DistributedMesh
1309  // remote_elem links as well. We need to save the (elem,
1310  // side, bc_id) triples and those links
1311  for (auto s : elem->side_index_range())
1312  {
1313  if (elem->neighbor_ptr(s) == remote_elem)
1314  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1315 
1316  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1317  for (const auto & bc_id : bc_ids)
1318  if (bc_id != BoundaryInfo::invalid_id)
1319  {
1320  saved_boundary_elements.push_back(copy);
1321  saved_bc_ids.push_back(bc_id);
1322  saved_bc_sides.push_back(s);
1323  }
1324  }
1325 
1326 
1327  // We're done with this element
1328  mesh.delete_elem(elem);
1329 
1330  // But save the copy
1331  new_elements.push_back(copy);
1332  }
1333 
1334  // Make sure we saved the same number of boundary conditions
1335  // in each vector.
1336  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1337  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1338 
1339  // Loop again, delete any remaining elements
1340  for (auto & elem : mesh.element_ptr_range())
1341  mesh.delete_elem(elem);
1342 
1343  // Add the copied (now level-0) elements back to the mesh
1344  for (auto & new_elem : new_elements)
1345  {
1346  // Save the original ID, because the act of adding the Elem can
1347  // change new_elem's id!
1348  dof_id_type orig_id = new_elem->id();
1349 
1350  Elem * added_elem = mesh.add_elem(new_elem);
1351 
1352  // If the Elem, as it was re-added to the mesh, now has a
1353  // different ID (this is unlikely, so it's just an assert)
1354  // the boundary information will no longer be correct.
1355  libmesh_assert_equal_to (orig_id, added_elem->id());
1356 
1357  // Avoid compiler warnings in opt mode.
1358  libmesh_ignore(added_elem, orig_id);
1359  }
1360 
1361  // Finally, also add back the saved boundary information
1362  for (std::size_t e=0; e<saved_boundary_elements.size(); ++e)
1363  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1364  saved_bc_sides[e],
1365  saved_bc_ids[e]);
1366 
1367  // Trim unused and renumber nodes and elements
1368  mesh.prepare_for_use(/*skip_renumber =*/ false);
1369 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
bc_id
boundary_id_type bc_id
Definition: xdr_io.C:51
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:461
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:64
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57

◆ redistribute()

void libMesh::MeshTools::Modification::redistribute ( MeshBase &  mesh,
const FunctionBase< Real > &  mapfunc 
)

Deterministically perturb the nodal locations. This function will move each node from it's current x/y/z coordinates to a new x/y/z coordinate given by the first LIBMESH_DIM components of the specified function mapfunc

Nodes on the boundary are also moved.

Currently, non-vertex nodes are moved in the same way as vertex nodes, according to (newx,newy,newz) = mapfunc(x,y,z). This behavior is often suboptimal for higher order geometries and may be subject to change in future libMesh versions.

Definition at line 138 of file mesh_modification.C.

References libMesh::FunctionBase< Output >::clone(), mesh, libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), and libMesh::MeshBase::node_ptr_range().

Referenced by libMesh::MeshTools::Generation::build_cube().

140 {
141  libmesh_assert (mesh.n_nodes());
142  libmesh_assert (mesh.n_elem());
143 
144  LOG_SCOPE("redistribute()", "MeshTools::Modification");
145 
146  DenseVector<Real> output_vec(LIBMESH_DIM);
147 
148  // FIXME - we should thread this later.
149  std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
150 
151  for (auto & node : mesh.node_ptr_range())
152  {
153  (*myfunc)(*node, output_vec);
154 
155  (*node)(0) = output_vec(0);
156 #if LIBMESH_DIM > 1
157  (*node)(1) = output_vec(1);
158 #endif
159 #if LIBMESH_DIM > 2
160  (*node)(2) = output_vec(2);
161 #endif
162  }
163 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::FunctionBase::clone
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0

◆ rotate()

void libMesh::MeshTools::Modification::rotate ( MeshBase &  mesh,
const Real  phi,
const Real  theta = 0.,
const Real  psi = 0. 
)

Rotates the mesh in the xy plane. The rotation is counter-clock-wise (mathematical definition). The angle is in degrees (360 make a full circle) Rotates the mesh in 3D space. Here the standard Euler angles are adopted (http://mathworld.wolfram.com/EulerAngles.html) The angles are in degrees (360 make a full circle)

Definition at line 201 of file mesh_modification.C.

References mesh, libMesh::MeshBase::node_ptr_range(), libMesh::pi, and libMesh::Real.

205 {
206 #if LIBMESH_DIM == 3
207  const Real p = -phi/180.*libMesh::pi;
208  const Real t = -theta/180.*libMesh::pi;
209  const Real s = -psi/180.*libMesh::pi;
210  const Real sp = std::sin(p), cp = std::cos(p);
211  const Real st = std::sin(t), ct = std::cos(t);
212  const Real ss = std::sin(s), cs = std::cos(s);
213 
214  // We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
215  // (equations 6-14 give the entries of the composite transformation matrix).
216  // The rotations are performed sequentially about the z, x, and z axes, in that order.
217  // A positive angle yields a counter-clockwise rotation about the axis in question.
218  for (auto & node : mesh.node_ptr_range())
219  {
220  const Point pt = *node;
221  const Real x = pt(0);
222  const Real y = pt(1);
223  const Real z = pt(2);
224  *node = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
225  (-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
226  ( sp*st)*x + (-cp*st)*y + (ct)*z );
227  }
228 #else
229  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
230 #endif
231 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::pi
const Real pi
Definition: libmesh.h:233

◆ scale()

void libMesh::MeshTools::Modification::scale ( MeshBase &  mesh,
const Real  xs,
const Real  ys = 0.,
const Real  zs = 0. 
)

Scales the mesh. The grid points are scaled in the x direction by xs, in the y direction by ys, etc... If only xs is specified then the scaling is assumed uniform in all directions.

Definition at line 234 of file mesh_modification.C.

References mesh, libMesh::MeshBase::node_ptr_range(), and libMesh::Real.

Referenced by libMesh::DenseVector< Output >::operator*=(), and libMesh::DenseMatrix< Number >::operator*=().

238 {
239  const Real x_scale = xs;
240  Real y_scale = ys;
241  Real z_scale = zs;
242 
243  if (ys == 0.)
244  {
245  libmesh_assert_equal_to (zs, 0.);
246 
247  y_scale = z_scale = x_scale;
248  }
249 
250  // Scale the x coordinate in all dimensions
251  for (auto & node : mesh.node_ptr_range())
252  (*node)(0) *= x_scale;
253 
254  // Only scale the y coordinate in 2 and 3D
255  if (LIBMESH_DIM < 2)
256  return;
257 
258  for (auto & node : mesh.node_ptr_range())
259  (*node)(1) *= y_scale;
260 
261  // Only scale the z coordinate in 3D
262  if (LIBMESH_DIM < 3)
263  return;
264 
265  for (auto & node : mesh.node_ptr_range())
266  (*node)(2) *= z_scale;
267 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131

◆ smooth()

void libMesh::MeshTools::Modification::smooth ( MeshBase &  mesh,
unsigned int  n_iterations,
Real  power 
)

Smooth the mesh with a simple Laplace smoothing algorithm. The mesh is smoothed n_iterations times. If the parameter power is 0, each node is moved to the average position of the neighboring connected nodes. If power > 0, the node positions are weighted by their distance. The positions of higher order nodes, and nodes living in refined elements, are calculated from the vertex positions of their parent nodes. Only works in 2D.

Author
Martin Luthi (luthi.nosp@m.@gi..nosp@m.alask.nosp@m.a.ed.nosp@m.u)
Date
2005

This implementation assumes every element "side" has only 2 nodes.

Definition at line 1106 of file mesh_modification.C.

References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::as_range(), libMesh::Elem::build_side_ptr(), libMesh::Elem::embedding_matrix(), libMesh::MeshTools::find_boundary_nodes(), libMesh::DofObject::id(), libMesh::MeshBase::level_elements_begin(), libMesh::MeshBase::level_elements_end(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::MeshTools::n_levels(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_nodes(), libMesh::Elem::n_second_order_adjacent_vertices(), libMesh::Elem::n_vertices(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_index_range(), libMesh::Elem::node_ptr(), libMesh::MeshBase::node_ref(), libMesh::TypeVector< T >::norm(), libMesh::Elem::parent(), libMesh::Elem::point(), std::pow(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), side, libMesh::Elem::side_index_range(), libMesh::MeshTools::weight(), and libMesh::Elem::which_child_am_i().

1109 {
1113  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1114 
1115  /*
1116  * Create a quickly-searchable list of boundary nodes.
1117  */
1118  std::unordered_set<dof_id_type> boundary_node_ids =
1119  MeshTools::find_boundary_nodes (mesh);
1120 
1121  for (unsigned int iter=0; iter<n_iterations; iter++)
1122  {
1123  /*
1124  * loop over the mesh refinement level
1125  */
1126  unsigned int n_levels = MeshTools::n_levels(mesh);
1127  for (unsigned int refinement_level=0; refinement_level != n_levels;
1128  refinement_level++)
1129  {
1130  // initialize the storage (have to do it on every level to get empty vectors
1131  std::vector<Point> new_positions;
1132  std::vector<Real> weight;
1133  new_positions.resize(mesh.n_nodes());
1134  weight.resize(mesh.n_nodes());
1135 
1136  {
1137  // Loop over the elements to calculate new node positions
1138  for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1139  mesh.level_elements_end(refinement_level)))
1140  {
1141  /*
1142  * We relax all nodes on level 0 first
1143  * If the element is refined (level > 0), we interpolate the
1144  * parents nodes with help of the embedding matrix
1145  */
1146  if (refinement_level == 0)
1147  {
1148  for (auto s : elem->side_index_range())
1149  {
1150  /*
1151  * Only operate on sides which are on the
1152  * boundary or for which the current element's
1153  * id is greater than its neighbor's.
1154  * Sides get only built once.
1155  */
1156  if ((elem->neighbor_ptr(s) != nullptr) &&
1157  (elem->id() > elem->neighbor_ptr(s)->id()))
1158  {
1159  std::unique_ptr<const Elem> side(elem->build_side_ptr(s));
1160 
1161  const Node & node0 = side->node_ref(0);
1162  const Node & node1 = side->node_ref(1);
1163 
1164  Real node_weight = 1.;
1165  // calculate the weight of the nodes
1166  if (power > 0)
1167  {
1168  Point diff = node0-node1;
1169  node_weight = std::pow(diff.norm(), power);
1170  }
1171 
1172  const dof_id_type id0 = node0.id(), id1 = node1.id();
1173  new_positions[id0].add_scaled( node1, node_weight );
1174  new_positions[id1].add_scaled( node0, node_weight );
1175  weight[id0] += node_weight;
1176  weight[id1] += node_weight;
1177  }
1178  } // element neighbor loop
1179  }
1180 #ifdef LIBMESH_ENABLE_AMR
1181  else // refinement_level > 0
1182  {
1183  /*
1184  * Find the positions of the hanging nodes of refined elements.
1185  * We do this by calculating their position based on the parent
1186  * (one level less refined) element, and the embedding matrix
1187  */
1188 
1189  const Elem * parent = elem->parent();
1190 
1191  /*
1192  * find out which child I am
1193  */
1194  unsigned int c = parent->which_child_am_i(elem);
1195  /*
1196  *loop over the childs (that is, the current elements) nodes
1197  */
1198  for (auto nc : elem->node_index_range())
1199  {
1200  /*
1201  * the new position of the node
1202  */
1203  Point point;
1204  for (auto n : parent->node_index_range())
1205  {
1206  /*
1207  * The value from the embedding matrix
1208  */
1209  const float em_val = parent->embedding_matrix(c,nc,n);
1210 
1211  if (em_val != 0.)
1212  point.add_scaled (parent->point(n), em_val);
1213  }
1214 
1215  const dof_id_type id = elem->node_ptr(nc)->id();
1216  new_positions[id] = point;
1217  weight[id] = 1.;
1218  }
1219  } // if element refinement_level
1220 #endif // #ifdef LIBMESH_ENABLE_AMR
1221 
1222  } // element loop
1223 
1224  /*
1225  * finally reposition the vertex nodes
1226  */
1227  for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
1228  if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1229  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1230  }
1231 
1232  // Now handle the additional second_order nodes by calculating
1233  // their position based on the vertex positions
1234  // we do a second loop over the level elements
1235  for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1236  mesh.level_elements_end(refinement_level)))
1237  {
1238  const unsigned int son_begin = elem->n_vertices();
1239  const unsigned int son_end = elem->n_nodes();
1240  for (unsigned int n=son_begin; n<son_end; n++)
1241  {
1242  const unsigned int n_adjacent_vertices =
1243  elem->n_second_order_adjacent_vertices(n);
1244 
1245  Point point;
1246  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1247  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1248 
1249  const dof_id_type id = elem->node_ptr(n)->id();
1250  mesh.node_ref(id) = point/n_adjacent_vertices;
1251  }
1252  }
1253  } // refinement_level loop
1254  } // end iteration
1255 }
side
unsigned short int side
Definition: xdr_io.C:50
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::MeshTools::weight
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:233
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:653
libMesh::as_range
SimpleRange< I > as_range(const std::pair< I, I > &p)
Definition: simple_range.h:57
std::pow
double pow(double a, int b)
Definition: libmesh_augment_std_namespace.h:66
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::MeshTools::find_boundary_nodes
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Definition: mesh_tools.C:303
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:64

◆ translate()

void libMesh::MeshTools::Modification::translate ( MeshBase &  mesh,
const Real  xt = 0.,
const Real  yt = 0.,
const Real  zt = 0. 
)

Translates the mesh. The grid points are translated in the x direction by xt, in the y direction by yt, etc...

Definition at line 167 of file mesh_modification.C.

References mesh, and libMesh::MeshBase::node_ptr_range().

171 {
172  const Point p(xt, yt, zt);
173 
174  for (auto & node : mesh.node_ptr_range())
175  *node += p;
176 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
generated by doxygen