410 {
411 LOG_UNIT_TEST;
412
413 if (multiple_subdomains && tangle_mesh)
414 libmesh_not_implemented_msg(
415 "Arbitrary mesh tangling with multiple subdomains is not supported.");
416
417
419 const auto dim = ref_elem->dim();
424
425
426
428 const auto side_length = 1.0 / n_elems_per_side;
429
431 {
432 case 1:
434 break;
435 case 2:
437 mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
438 break;
439
440 case 3:
442 n_elems_per_side,
443 n_elems_per_side,
444 n_elems_per_side,
445 0.,
446 1.,
447 0.,
448 1.,
449 0.,
450 1.,
451 type);
452 break;
453
454 default:
455 libmesh_error_msg(
"Unsupported dimension " <<
dim);
456 }
457
458
459 DistortHyperCube dh(
dim);
461
464
465 if (tangle_mesh)
466 {
467
468
469
470
471
472
473
477
478
479
480
481
482 std::map<dof_id_type, Real> dist1_map;
483
484 std::map<dof_id_type, Real> dist2_map;
485
486 for (
const auto * node :
mesh.local_node_ptr_range())
487 {
488 dist1_map[node->id()] = (*node - p1).
norm();
489 dist2_map[node->id()] = (*node - p2).
norm();
490 }
491
492
493
494
495 auto get_closet_point_accross_all_procs =
496 [&
mesh, &proc_id](
const std::map<dof_id_type, Real> &dist_map) {
497
499 Real d_min_local = std::numeric_limits<Real>::max();
501
502
503
504 if (dist_map.size())
505 {
506
507 auto min_it =
508 std::min_element(dist_map.begin(), dist_map.end(),
509 [](
const auto &a,
const auto &
b) {
510 return a.second < b.second;
511 });
512
513 node_id = min_it->first;
514 d_min_local = min_it->second;
515 }
516
517
518 auto d_min_global = d_min_local;
520
521
522
524 if (d_min_local == d_min_global) {
525 broadcasting_proc = proc_id;
527 }
528
530
533
534 return std::pair(node_id, node);
535 };
536
537 const auto [node_id1, node1] =
538 get_closet_point_accross_all_procs(dist1_map);
539 const auto [node_id2, node2] =
540 get_closet_point_accross_all_procs(dist2_map);
541
542
543 const auto displacement = tangle_damping_factor * (node1 - node2);
548
551 mesh.nodes_end(), sync_object);
552
553
554
557 CPPUNIT_ASSERT(unsmoothed_info.mesh_is_tangled);
558 }
559
560
561 std::unordered_map<subdomain_id_type, Real> distorted_subdomain_volumes;
563 if (multiple_subdomains)
564 {
565
566 for (
auto * elem :
mesh.active_element_ptr_range())
567 {
568 unsigned int subdomain_id = 0;
570 if (elem->vertex_average()(d) > 0.5)
571 ++subdomain_id;
572 elem->subdomain_id() += subdomain_id;
573 }
574
575
576
577
578 for (
auto *elem :
mesh.active_element_ptr_range()) {
579 const auto sub_id = elem->subdomain_id();
580
581 if (elem->processor_id() != proc_id)
582 continue;
583 distorted_subdomain_volumes[sub_id] += elem->volume();
584 highest_subdomain_id = std::max(sub_id, highest_subdomain_id);
585 }
586
587
589 for (
const auto sub_id :
make_range(highest_subdomain_id + 1)) {
590
591 if (distorted_subdomain_volumes.find(sub_id) ==
592 distorted_subdomain_volumes.end())
593 distorted_subdomain_volumes[sub_id] = 0.;
594
595 mesh.
comm().
sum(distorted_subdomain_volumes[sub_id]);
596 }
597 }
598
599
600
601
602
603
605
606
608 libmesh_error_msg_if(elem_orders.size() != 1,
609 "The variational smoother cannot be used for mixed-order meshes!");
610
611
612
613
614 const auto scale_factor = *elem_orders.begin() * ((type_is_pyramid || type_is_tet) ? 2 * 4 : 1);
615
616
617 auto node_distortion_is = [&n_elems_per_side, &
dim, &boundary_info, &scale_factor, &type_is_prism](
619
620 std::vector<boundary_id_type> boundary_ids;
621 boundary_info.boundary_ids(&node, boundary_ids);
622
623
624 const auto num_dofs =
dim - boundary_ids.size();
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644 std::size_t num_zero_or_one = 0;
645
646 bool distorted = false;
648 {
649 const Real r = node(d);
650 const Real R = r * n_elems_per_side * scale_factor;
651 CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, r);
652 CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, 1 - r);
653
654 bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
655 if (type_is_prism && (scale_factor == 3))
656 {
657
658
659
660 const Real R_prism = R / scale_factor * 2;
661 const bool d_distorted_prism =
662 std::abs(R_prism - std::round(R_prism)) > distortion_tol;
663 d_distorted &= d_distorted_prism;
664 }
665 distorted |= d_distorted;
667 }
668
669 CPPUNIT_ASSERT_GREATEREQUAL(
dim - num_dofs, num_zero_or_one);
670
671
672 if (num_dofs == 0)
673
674 return true;
675 return distorted == distortion;
676 };
677
678
679 for (
auto node :
mesh.node_ptr_range())
680 CPPUNIT_ASSERT(node_distortion_is(*node, true));
681
682
683
684
685 if (type_is_tri)
686 {
687
688
689
690 SquareToParallelogram stp;
692 }
693 else if (type_is_prism)
694 {
695
696
697
698 CubeToParallelepiped ctp;
700 }
701
703
704 if (type_is_tri)
705 {
706
707
708
709 ParallelogramToSquare pts;
711 }
712 else if (type_is_prism)
713 {
714
715
716
717 ParallelepipedToCube ptc;
719 }
720
721 if (multiple_subdomains)
722 {
723
724
725
726 std::unordered_map<subdomain_id_type, Real> smoothed_subdomain_volumes;
727 for (
auto *elem :
mesh.active_element_ptr_range()) {
728
729 if (elem->processor_id() != proc_id)
730 continue;
731 smoothed_subdomain_volumes[elem->subdomain_id()] += elem->volume();
732 }
733
734
735 for (
const auto sub_id :
make_range(highest_subdomain_id + 1)) {
736
737 if (smoothed_subdomain_volumes.find(sub_id) ==
738 smoothed_subdomain_volumes.end())
739 smoothed_subdomain_volumes[sub_id] = 0.;
740
741 mesh.
comm().
sum(smoothed_subdomain_volumes[sub_id]);
742 }
743
744 for (
const auto sub_id :
make_range(highest_subdomain_id + 1))
746 libmesh_map_find(distorted_subdomain_volumes, sub_id),
747 libmesh_map_find(smoothed_subdomain_volumes, sub_id),
TOLERANCE));
748 }
749 else
750 {
751 std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
753
754
755
756 std::set<dof_id_type> nodes_checked;
758
759 for (
const auto * elem :
mesh.active_element_ptr_range())
760 {
762 {
763 const auto & node = elem->node_ref(local_node_id);
764 if (nodes_checked.find(node.id()) != nodes_checked.end())
765 continue;
766
767 nodes_checked.insert(node.id());
768
769
770
771 if (type_is_pyramid)
772 {
773 if (local_node_id > 8 && local_node_id < 13)
774 {
775
776
777
778
779
780
781
782
783
784
785 const auto & base = elem->node_ref(local_node_id - 9);
786 const auto & apex = elem->node_ref(4);
787 const Real x = (type ==
PYRAMID18) ? 0.56646084 : 0.54985875;
788
789 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(base + x * (apex - base), tol));
790 continue;
791 }
792 else if (local_node_id > 13)
793 {
794
795
796
797
798
799
800
801
802 const auto & base1 = elem->node_ref(local_node_id - 14);
803 const auto & base2 = elem->node_ref((local_node_id - 13) % 4);
804 const auto & apex = elem->node_ref(4);
805
806 const auto node_approx =
807 (0.31401599 * base1 + 0.31401599 * base2 + 0.37196802 * apex);
808 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(node_approx, tol));
809 continue;
810 }
811 }
812
813
814
815 else if (type_is_tet && !elem->is_vertex(local_node_id))
816 {
817
818
819
820 std::vector<const Node *> neighbors;
822 mesh, node, nodes_to_elem_map, neighbors);
823
824 switch (neighbors.size())
825 {
826 case 2: {
827
828
829
830
831 const auto is_0_cube_center =
833 const auto is_1_cube_center =
836
837 if (is_0_cube_center || is_1_cube_center)
838 {
839 const auto & cube_center =
840 is_0_cube_center ? *neighbors[0] : *neighbors[1];
841 const auto & other =
842 is_0_cube_center ? *neighbors[1] : *neighbors[0];
843
844
845
846
848 {
849 const Real x = (type ==
TET10) ? 0.42895041 : 0.41486385;
850 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(
851 other + x * (cube_center - other), tol));
852 }
853
855 {
856 const Real x = (type ==
TET10) ? 0.55388920 : 0.58093516;
857 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(
858 other + x * (cube_center - other), tol));
859 }
860 }
861
862 else
863 {
864
865
866
867 const auto is_0_cube_vertex =
869 const auto is_1_cube_vertex =
871 const auto is_0_cube_face_center =
873 const auto is_1_cube_face_center =
875
876 if (is_0_cube_vertex && is_1_cube_vertex)
877
878
879
880 CPPUNIT_ASSERT(node_distortion_is(node, false));
881
882 else
883 {
884 libmesh_error_msg_if(
885 (is_0_cube_center || is_0_cube_face_center) &&
886 (is_1_cube_center || is_1_cube_face_center),
887 "We should never get here!");
888 const auto & cube_vertex =
889 is_0_cube_vertex ? *neighbors[0] : *neighbors[1];
890 const auto & cube_face_center =
891 is_0_cube_face_center ? *neighbors[0] : *neighbors[1];
892 const Real x = (type ==
TET10) ? 0.61299101 : 0.65125580;
893 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(
894 cube_vertex + x * (cube_face_center - cube_vertex), tol));
895 }
896 }
897
898 continue;
899 break;
900 }
901
902 case 6: {
903
904
905 neighbors.erase(std::remove_if(neighbors.begin(),
906 neighbors.end(),
907 [&elem](
const Node * n) {
908 return !elem->is_vertex(
909 elem->local_node(n->id()));
910 }),
911 neighbors.end());
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926 unsigned int num_vertices_at_cube_center = 0;
927 unsigned int num_vertices_at_cube_vertices = 0;
928 unsigned int num_vertices_at_cube_face_centers = 0;
929 for (const auto * neighbor : neighbors)
930 {
932 num_vertices_at_cube_center += 1;
934 num_vertices_at_cube_vertices += 1;
936 num_vertices_at_cube_face_centers += 1;
937 }
938
939
940
942
943
944 if (num_vertices_at_cube_vertices == 2)
945 {
946
947 if (num_vertices_at_cube_center)
948 for (const auto * neighbor : neighbors)
949 {
955 else
956 libmesh_error_msg("We should never get here!");
957
958 node_approx +=
weight * (*neighbor);
959 }
960
961
962 else if (num_vertices_at_cube_face_centers)
963 for (const auto * neighbor : neighbors)
964 {
970 else
971 libmesh_error_msg("We should never get here!");
972
973 node_approx +=
weight * (*neighbor);
974 }
975
976 else
977 libmesh_error_msg("We should never get here!");
978 }
979
980
981 else if (num_vertices_at_cube_center && num_vertices_at_cube_vertices &&
982 num_vertices_at_cube_face_centers)
983 for (const auto * neighbor : neighbors)
984 {
992 else
993 libmesh_error_msg("We should never get here!");
994
995 node_approx +=
weight * (*neighbor);
996 }
997
998 else
999 libmesh_error_msg("We should never get here!");
1000
1001 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(node_approx, tol));
1002
1003 continue;
1004 break;
1005 }
1006 default: {
1007 libmesh_error_msg(node << " has unexpected number of neighbors ("
1008 << neighbors.size() << ")");
1009 break;
1010 }
1011 }
1012 }
1013
1014 CPPUNIT_ASSERT(node_distortion_is(node, false));
1015
1016 }
1017 }
1018 }
1019 }
bool pointIsCubeFaceCenter(const Point &point, const Real &side_length, const Real &tol=TOLERANCE)
bool pointIsCubeCenter(const Point &point, const Real &side_length, const Real &tol=TOLERANCE)
bool pointIsCubeVertex(const Point &point, const Real &side_length, const Real &tol=TOLERANCE)
void max(const T &r, T &o, Request &req) const
void min(const T &r, T &o, Request &req) const
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
static constexpr dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
static const Order type_to_default_order_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the default approximation order of...
virtual const Node & node_ref(const dof_id_type i) const
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
virtual const Node * query_node_ptr(const dof_id_type i) const =0
const std::set< Order > & elem_default_orders() const
processor_id_type processor_id() const
const Parallel::Communicator & comm() const
A Point defines a location in LIBMESH_DIM dimensional Real space.
virtual void smooth() override
Redefinition of the smooth function from the base class.
const MeshQualityInfo & get_mesh_info() const
Getter for the _system's _mesh_info attribute.
virtual void setup()
Setup method that creates equation systems, system, and constraints, to be called just prior to smoot...
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
const Elem & get(const ElemType type_in)
std::string enum_to_string(const T e)
bool relative_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
bool absolute_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within an absolute tolerance.
uint8_t processor_id_type
const dof_id_type n_nodes
