libMesh
Functions
libMesh::MeshTools::Modification Namespace Reference

Tools for Mesh modification. More...

Functions

void distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
 Randomly perturb the nodal locations. More...
 
void permute_elements (MeshBase &mesh)
 Randomly permute the nodal ordering of each element (without twisting the element mapping). More...
 
void orient_elements (MeshBase &mesh)
 Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientation. More...
 
void redistribute (MeshBase &mesh, const FunctionBase< Real > &mapfunc)
 Deterministically perturb the nodal locations. More...
 
void translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
 Translates the mesh. More...
 
RealTensorValue rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
 Rotates the mesh in the xy plane. More...
 
void scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
 Scales the mesh. More...
 
void all_tri (MeshBase &mesh)
 Converts the 2D quadrilateral elements of a Mesh into triangular elements. More...
 
void smooth (MeshBase &, unsigned int, Real)
 Smooth the mesh with a simple Laplace smoothing algorithm. More...
 
void flatten (MeshBase &mesh)
 Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elements. More...
 
void 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. More...
 
void 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. More...
 

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 333 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(), libMesh::Elem::build_side_ptr(), libMesh::ParallelObject::comm(), libMesh::MeshBase::delete_elem(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::get_extra_integer(), libMesh::DofObject::id(), libMesh::index_range(), libMesh::INFEDGE2, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::DofObject::invalid_id, libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::is_prepared(), libMesh::MeshBase::is_replicated(), libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::Elem::mapping_data(), libMesh::Elem::mapping_type(), TIMPI::Communicator::max(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::MeshBase::n_elem(), libMesh::DofObject::n_extra_integers(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_id(), libMesh::Elem::node_ptr(), libMesh::TensorTools::norm(), libMesh::MeshBase::parallel_max_unique_id(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::point(), libMesh::MeshBase::prepare_for_use(), libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::remote_elem, libMesh::BoundaryInfo::remove(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, libMesh::Elem::type(), and libMesh::DofObject::unique_id().

Referenced by OverlappingFunctorTest::checkCouplingFunctorTri(), OverlappingFunctorTest::checkCouplingFunctorTriUnifRef(), main(), AllTriTest::test_helper_2D(), and AllTriTest::test_helper_3D().

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

◆ 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 1464 of file mesh_modification.C.

References libMesh::MeshBase::get_boundary_info(), mesh, and libMesh::BoundaryInfo::renumber_id().

Referenced by MeshStitchTest::renameAndShift().

1467 {
1468  // This is just a shim around the member implementation, now
1469  mesh.get_boundary_info().renumber_id(old_id, new_id);
1470 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308

◆ 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 1474 of file mesh_modification.C.

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

1477 {
1478  if (old_id == new_id)
1479  {
1480  // If the IDs are the same, this is a no-op.
1481  return;
1482  }
1483 
1484  for (auto & elem : mesh.element_ptr_range())
1485  {
1486  if (elem->subdomain_id() == old_id)
1487  elem->subdomain_id() = new_id;
1488  }
1489 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308

◆ 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 68 of file mesh_modification.C.

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

Referenced by main(), DistortTest::perturb_and_check(), VolumeTest::test_true_centroid_and_volume(), and VolumeTest::testQuad4TrueCentroid().

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

◆ 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 1338 of file mesh_modification.C.

References libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::get_extra_integer(), libMesh::DofObject::id(), libMesh::index_range(), libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::is_prepared(), libMesh::MeshBase::is_replicated(), libMesh::libmesh_assert(), libMesh::libmesh_ignore(), libMesh::Elem::mapping_data(), libMesh::Elem::mapping_type(), mesh, libMesh::MeshBase::n_active_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::DofObject::n_extra_integers(), 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::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::Elem::type(), and libMesh::DofObject::unique_id().

Referenced by main().

1339 {
1340  libmesh_assert(mesh.is_prepared() || mesh.is_replicated());
1341 
1342  // Algorithm:
1343  // .) For each active element in the mesh: construct a
1344  // copy which is the same in every way *except* it is
1345  // a level 0 element. Store the pointers to these in
1346  // a separate vector. Save any boundary information as well.
1347  // Delete the active element from the mesh.
1348  // .) Loop over all (remaining) elements in the mesh, delete them.
1349  // .) Add the level-0 copies back to the mesh
1350 
1351  // Temporary storage for new element pointers
1352  std::vector<std::unique_ptr<Elem>> new_elements;
1353 
1354  // BoundaryInfo Storage for element ids, sides, and BC ids
1355  std::vector<Elem *> saved_boundary_elements;
1356  std::vector<boundary_id_type> saved_bc_ids;
1357  std::vector<unsigned short int> saved_bc_sides;
1358 
1359  // Container to catch boundary ids passed back by BoundaryInfo
1360  std::vector<boundary_id_type> bc_ids;
1361 
1362  // Reserve a reasonable amt. of space for each
1363  new_elements.reserve(mesh.n_active_elem());
1364  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1365  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1366  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1367 
1368  for (auto & elem : mesh.active_element_ptr_range())
1369  {
1370  // Make a new element of the same type
1371  auto copy = Elem::build(elem->type());
1372 
1373  // Set node pointers (they still point to nodes in the original mesh)
1374  for (auto n : elem->node_index_range())
1375  copy->set_node(n) = elem->node_ptr(n);
1376 
1377  // Copy over ids
1378  copy->processor_id() = elem->processor_id();
1379  copy->subdomain_id() = elem->subdomain_id();
1380 
1381  // Retain the original element's ID(s) as well, otherwise
1382  // the Mesh may try to create them for you...
1383  copy->set_id( elem->id() );
1384 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1385  copy->set_unique_id(elem->unique_id());
1386 #endif
1387 
1388  // This element could have boundary info or DistributedMesh
1389  // remote_elem links as well. We need to save the (elem,
1390  // side, bc_id) triples and those links
1391  for (auto s : elem->side_index_range())
1392  {
1393  if (elem->neighbor_ptr(s) == remote_elem)
1394  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1395 
1396  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1397  for (const auto & bc_id : bc_ids)
1398  if (bc_id != BoundaryInfo::invalid_id)
1399  {
1400  saved_boundary_elements.push_back(copy.get());
1401  saved_bc_ids.push_back(bc_id);
1402  saved_bc_sides.push_back(s);
1403  }
1404  }
1405 
1406  // Copy any extra element data. Since the copy hasn't been
1407  // added to the mesh yet any allocation has to be done manually.
1408  const unsigned int nei = elem->n_extra_integers();
1409  copy->add_extra_integers(nei);
1410  for (unsigned int i=0; i != nei; ++i)
1411  copy->set_extra_integer(i, elem->get_extra_integer(i));
1412 
1413  // Copy any mapping data.
1414  copy->set_mapping_type(elem->mapping_type());
1415  copy->set_mapping_data(elem->mapping_data());
1416 
1417  // We're done with this element
1418  mesh.delete_elem(elem);
1419 
1420  // But save the copy
1421  new_elements.push_back(std::move(copy));
1422  }
1423 
1424  // Make sure we saved the same number of boundary conditions
1425  // in each vector.
1426  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1427  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1428 
1429  // Loop again, delete any remaining elements
1430  for (auto & elem : mesh.element_ptr_range())
1431  mesh.delete_elem(elem);
1432 
1433  // Add the copied (now level-0) elements back to the mesh
1434  for (auto & new_elem : new_elements)
1435  {
1436  // Save the original ID, because the act of adding the Elem can
1437  // change new_elem's id!
1438  dof_id_type orig_id = new_elem->id();
1439 
1440  Elem * added_elem = mesh.add_elem(std::move(new_elem));
1441 
1442  // If the Elem, as it was re-added to the mesh, now has a
1443  // different ID (this is unlikely, so it's just an assert)
1444  // the boundary information will no longer be correct.
1445  libmesh_assert_equal_to (orig_id, added_elem->id());
1446 
1447  // Avoid compiler warnings in opt mode.
1448  libmesh_ignore(added_elem, orig_id);
1449  }
1450 
1451  // Finally, also add back the saved boundary information
1452  for (auto e : index_range(saved_boundary_elements))
1453  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1454  saved_bc_sides[e],
1455  saved_bc_ids[e]);
1456 
1457  // Trim unused and renumber nodes and elements
1458  mesh.prepare_for_use();
1459 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:516
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:111
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:54

◆ orient_elements()

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

Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientation.

This function does not currently handle meshes with any element refinement.

Definition at line 189 of file mesh_modification.C.

References libMesh::MeshBase::get_boundary_info(), mesh, libMesh::MeshTools::n_levels(), and libMesh::Elem::orient().

Referenced by ElemTest< elem_type >::test_orient_elements().

190 {
191  LOG_SCOPE("orient_elements()", "MeshTools::Modification");
192 
193  // We don't yet support doing orient() on a parent element, which
194  // would require us to consistently orient all its children and
195  // give them different local child numbers.
196  unsigned int n_levels = MeshTools::n_levels(mesh);
197  if (n_levels > 1)
198  libmesh_not_implemented_msg("orient_elements() does not support refined meshes");
199 
200  BoundaryInfo & boundary_info = mesh.get_boundary_info();
201  for (auto elem : mesh.element_ptr_range())
202  elem->orient(&boundary_info);
203 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:801

◆ permute_elements()

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

Randomly permute the nodal ordering of each element (without twisting the element mapping).

This is useful for regression testing with a variety of element orientations.

This function does not currently handle meshes with any element refinement.

This function does not currently permute BoundaryInfo data associated with element sides, which will likely be scrambled by the permutation.

Definition at line 149 of file mesh_modification.C.

References libMesh::make_range(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshTools::n_levels(), libMesh::Elem::n_permutations(), libMesh::Elem::permute(), and libMesh::MeshBase::query_elem_ptr().

Referenced by main(), and FETestBase< order, family, elem_type, 1 >::setUp().

150 {
151  LOG_SCOPE("permute_elements()", "MeshTools::Modification");
152 
153  // We don't yet support doing permute() on a parent element, which
154  // would require us to consistently permute all its children and
155  // give them different local child numbers.
156  unsigned int n_levels = MeshTools::n_levels(mesh);
157  if (n_levels > 1)
158  libmesh_error();
159 
160  const unsigned int seed = 123456;
161 
162  // seed the random number generator.
163  // We'll loop from 1 to max_elem_id on every processor, even those
164  // that don't have a particular element, so that the pseudorandom
165  // numbers will be the same everywhere.
166  std::srand(seed);
167 
168 
169  for (auto e_id : make_range(mesh.max_elem_id()))
170  {
171  int my_rand = std::rand();
172 
173  Elem * elem = mesh.query_elem_ptr(e_id);
174 
175  if (!elem)
176  continue;
177 
178  const unsigned int max_permutation = elem->n_permutations();
179  if (!max_permutation)
180  continue;
181 
182  const unsigned int perm = my_rand % max_permutation;
183 
184  elem->permute(perm);
185  }
186 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:801
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134

◆ 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 207 of file mesh_modification.C.

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

Referenced by libMesh::MeshTools::Generation::build_cube(), and FETestBase< order, family, elem_type, 1 >::setUp().

209 {
210  libmesh_assert (mesh.n_nodes());
211  libmesh_assert (mesh.n_elem());
212 
213  LOG_SCOPE("redistribute()", "MeshTools::Modification");
214 
215  DenseVector<Real> output_vec(LIBMESH_DIM);
216 
217  // FIXME - we should thread this later.
218  std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
219 
220  for (auto & node : mesh.node_ptr_range())
221  {
222  (*myfunc)(*node, output_vec);
223 
224  (*node)(0) = output_vec(0);
225 #if LIBMESH_DIM > 1
226  (*node)(1) = output_vec(1);
227 #endif
228 #if LIBMESH_DIM > 2
229  (*node)(2) = output_vec(2);
230 #endif
231  }
232 }
libMesh::DenseVector< Real >
template class LIBMESH_EXPORT DenseVector< Real >
Definition: dense_vector.C:29
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FunctionBase::clone
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0

◆ rotate()

RealTensorValue 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)

Returns
the 3x3 rotation matrix implied by (phi, theta, psi)

Definition at line 271 of file mesh_modification.C.

References libMesh::TensorValue< T >::intrinsic_rotation_matrix(), libMesh::libmesh_ignore(), and mesh.

Referenced by libMesh::orientation(), libMesh::Prism6::permute(), libMesh::Hex8::permute(), libMesh::Hex20::permute(), libMesh::Tet10::permute(), libMesh::Tet14::permute(), libMesh::Prism15::permute(), libMesh::Tet4::permute(), libMesh::Hex27::permute(), libMesh::Prism18::permute(), libMesh::Prism20::permute(), libMesh::Prism21::permute(), FETestBase< order, family, elem_type, 1 >::setUp(), BBoxTest::test_no_degenerate(), and BBoxTest::test_two_degenerate().

275 {
276 #if LIBMESH_DIM == 3
277  const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
278 
279  for (auto & node : mesh.node_ptr_range())
280  {
281  Point & pt = *node;
282  pt = R * pt;
283  }
284 
285  return R;
286 
287 #else
288  libmesh_ignore(mesh, phi, theta, psi);
289  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
290  // We'll never get here
291  return RealTensorValue();
292 #endif
293 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308
libMesh::RealTensorValue
TensorValue< Real > RealTensorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
Definition: hp_coarsentest.h:48
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:516

◆ 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 296 of file mesh_modification.C.

References mesh, and libMesh::Real.

Referenced by libMesh::PetscVector< libMesh::Number >::add(), main(), libMesh::DenseVector< Output >::operator*=(), libMesh::DenseMatrix< Real >::operator*=(), and libMesh::RBEIMConstruction::train_eim_approximation_with_POD().

300 {
301  const Real x_scale = xs;
302  Real y_scale = ys;
303  Real z_scale = zs;
304 
305  if (ys == 0.)
306  {
307  libmesh_assert_equal_to (zs, 0.);
308 
309  y_scale = z_scale = x_scale;
310  }
311 
312  // Scale the x coordinate in all dimensions
313  for (auto & node : mesh.node_ptr_range())
314  (*node)(0) *= x_scale;
315 
316  // Only scale the y coordinate in 2 and 3D
317  if (LIBMESH_DIM < 2)
318  return;
319 
320  for (auto & node : mesh.node_ptr_range())
321  (*node)(1) *= y_scale;
322 
323  // Only scale the z coordinate in 3D
324  if (LIBMESH_DIM < 3)
325  return;
326 
327  for (auto & node : mesh.node_ptr_range())
328  (*node)(2) *= z_scale;
329 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:143

◆ 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 1182 of file mesh_modification.C.

References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::as_range(), libMesh::Elem::embedding_matrix(), libMesh::MeshTools::find_boundary_nodes(), libMesh::DofObject::id(), libMesh::make_range(), 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::Elem::node_ref(), libMesh::MeshBase::node_ref(), libMesh::TypeVector< T >::norm(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::Utility::pow(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), libMesh::Elem::side_index_range(), libMesh::MeshTools::weight(), and libMesh::Elem::which_child_am_i().

1185 {
1189  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1190 
1191  /*
1192  * Create a quickly-searchable list of boundary nodes.
1193  */
1194  std::unordered_set<dof_id_type> boundary_node_ids =
1195  MeshTools::find_boundary_nodes (mesh);
1196 
1197  // For avoiding extraneous element side allocation
1198  ElemSideBuilder side_builder;
1199 
1200  for (unsigned int iter=0; iter<n_iterations; iter++)
1201  {
1202  /*
1203  * loop over the mesh refinement level
1204  */
1205  unsigned int n_levels = MeshTools::n_levels(mesh);
1206  for (unsigned int refinement_level=0; refinement_level != n_levels;
1207  refinement_level++)
1208  {
1209  // initialize the storage (have to do it on every level to get empty vectors
1210  std::vector<Point> new_positions;
1211  std::vector<Real> weight;
1212  new_positions.resize(mesh.n_nodes());
1213  weight.resize(mesh.n_nodes());
1214 
1215  {
1216  // Loop over the elements to calculate new node positions
1217  for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1218  mesh.level_elements_end(refinement_level)))
1219  {
1220  /*
1221  * We relax all nodes on level 0 first
1222  * If the element is refined (level > 0), we interpolate the
1223  * parents nodes with help of the embedding matrix
1224  */
1225  if (refinement_level == 0)
1226  {
1227  for (auto s : elem->side_index_range())
1228  {
1229  /*
1230  * Only operate on sides which are on the
1231  * boundary or for which the current element's
1232  * id is greater than its neighbor's.
1233  * Sides get only built once.
1234  */
1235  if ((elem->neighbor_ptr(s) != nullptr) &&
1236  (elem->id() > elem->neighbor_ptr(s)->id()))
1237  {
1238  const Elem & side = side_builder(*elem, s);
1239  const Node & node0 = side.node_ref(0);
1240  const Node & node1 = side.node_ref(1);
1241 
1242  Real node_weight = 1.;
1243  // calculate the weight of the nodes
1244  if (power > 0)
1245  {
1246  Point diff = node0-node1;
1247  node_weight = std::pow(diff.norm(), power);
1248  }
1249 
1250  const dof_id_type id0 = node0.id(), id1 = node1.id();
1251  new_positions[id0].add_scaled( node1, node_weight );
1252  new_positions[id1].add_scaled( node0, node_weight );
1253  weight[id0] += node_weight;
1254  weight[id1] += node_weight;
1255  }
1256  } // element neighbor loop
1257  }
1258 #ifdef LIBMESH_ENABLE_AMR
1259  else // refinement_level > 0
1260  {
1261  /*
1262  * Find the positions of the hanging nodes of refined elements.
1263  * We do this by calculating their position based on the parent
1264  * (one level less refined) element, and the embedding matrix
1265  */
1266 
1267  const Elem * parent = elem->parent();
1268 
1269  /*
1270  * find out which child I am
1271  */
1272  unsigned int c = parent->which_child_am_i(elem);
1273  /*
1274  *loop over the childs (that is, the current elements) nodes
1275  */
1276  for (auto nc : elem->node_index_range())
1277  {
1278  /*
1279  * the new position of the node
1280  */
1281  Point point;
1282  for (auto n : parent->node_index_range())
1283  {
1284  /*
1285  * The value from the embedding matrix
1286  */
1287  const Real em_val = parent->embedding_matrix(c,nc,n);
1288 
1289  if (em_val != 0.)
1290  point.add_scaled (parent->point(n), em_val);
1291  }
1292 
1293  const dof_id_type id = elem->node_ptr(nc)->id();
1294  new_positions[id] = point;
1295  weight[id] = 1.;
1296  }
1297  } // if element refinement_level
1298 #endif // #ifdef LIBMESH_ENABLE_AMR
1299 
1300  } // element loop
1301 
1302  /*
1303  * finally reposition the vertex nodes
1304  */
1305  for (auto nid : make_range(mesh.n_nodes()))
1306  if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1307  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1308  }
1309 
1310  // Now handle the additional second_order nodes by calculating
1311  // their position based on the vertex positions
1312  // we do a second loop over the level elements
1313  for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1314  mesh.level_elements_end(refinement_level)))
1315  {
1316  const unsigned int son_begin = elem->n_vertices();
1317  const unsigned int son_end = elem->n_nodes();
1318  for (unsigned int n=son_begin; n<son_end; n++)
1319  {
1320  const unsigned int n_adjacent_vertices =
1321  elem->n_second_order_adjacent_vertices(n);
1322 
1323  Point point;
1324  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1325  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1326 
1327  const dof_id_type id = elem->node_ptr(n)->id();
1328  mesh.node_ref(id) = point/n_adjacent_vertices;
1329  }
1330  }
1331  } // refinement_level loop
1332  } // end iteration
1333 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308
libMesh::MeshTools::find_boundary_nodes
std::unordered_set< dof_id_type > find_boundary_nodes(const MeshBase &mesh)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:516
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:328
libMesh::MeshTools::weight
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:436
libMesh::as_range
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:801
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:143
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67

◆ 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 236 of file mesh_modification.C.

References mesh.

Referenced by MeshStitchTest::testMeshStitchElemsets(), and MeshTriangulationTest::testTriangulatorRoundHole().

240 {
241  const Point p(xt, yt, zt);
242 
243  for (auto & node : mesh.node_ptr_range())
244  *node += p;
245 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1308
Generated on Thu Apr 18 2024 15:19:21 for libMesh by   doxygen 1.8.14