408 {
409 LOG_UNIT_TEST;
410
411 if (multiple_subdomains && tangle_mesh)
412 libmesh_not_implemented_msg(
413 "Arbitrary mesh tangling with multiple subdomains is not supported.");
414
415
417 const auto dim = ref_elem->dim();
422
423
424
426 const auto side_length = 1.0 / n_elems_per_side;
427
429 {
430 case 1:
432 break;
433 case 2:
435 mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
436 break;
437
438 case 3:
440 n_elems_per_side,
441 n_elems_per_side,
442 n_elems_per_side,
443 0.,
444 1.,
445 0.,
446 1.,
447 0.,
448 1.,
449 type);
450 break;
451
452 default:
453 libmesh_error_msg(
"Unsupported dimension " <<
dim);
454 }
455
456
457 DistortHyperCube dh(
dim);
459
462
463 if (tangle_mesh)
464 {
465
466
467
468
469
470
471
475
476
477
478
479
480 std::map<dof_id_type, Real> dist1_map;
481
482 std::map<dof_id_type, Real> dist2_map;
483
484 for (
const auto * node :
mesh.local_node_ptr_range())
485 {
486 dist1_map[node->id()] = (*node - p1).
norm();
487 dist2_map[node->id()] = (*node - p2).
norm();
488 }
489
490
491
492
493 auto get_closet_point_accross_all_procs =
494 [&
mesh, &proc_id](
const std::map<dof_id_type, Real> &dist_map) {
495
497 Real d_min_local = std::numeric_limits<Real>::max();
499
500
501
502 if (dist_map.size())
503 {
504
505 auto min_it =
506 std::min_element(dist_map.begin(), dist_map.end(),
507 [](
const auto &a,
const auto &
b) {
508 return a.second < b.second;
509 });
510
511 node_id = min_it->first;
512 d_min_local = min_it->second;
513 }
514
515
516 auto d_min_global = d_min_local;
518
519
520
522 if (d_min_local == d_min_global) {
523 broadcasting_proc = proc_id;
525 }
526
528
531
532 return std::pair(node_id, node);
533 };
534
535 const auto [node_id1, node1] =
536 get_closet_point_accross_all_procs(dist1_map);
537 const auto [node_id2, node2] =
538 get_closet_point_accross_all_procs(dist2_map);
539
540
541 const auto displacement = tangle_damping_factor * (node1 - node2);
546
549 mesh.nodes_end(), sync_object);
550
551
552
555 CPPUNIT_ASSERT(unsmoothed_info.mesh_is_tangled);
556 }
557
558
559 std::unordered_map<subdomain_id_type, Real> distorted_subdomain_volumes;
561 if (multiple_subdomains)
562 {
563
564 for (
auto * elem :
mesh.active_element_ptr_range())
565 {
566 unsigned int subdomain_id = 0;
568 if (elem->vertex_average()(d) > 0.5)
569 ++subdomain_id;
570 elem->subdomain_id() += subdomain_id;
571 }
572
573
574
575
576 for (
auto *elem :
mesh.active_element_ptr_range()) {
577 const auto sub_id = elem->subdomain_id();
578
579 if (elem->processor_id() != proc_id)
580 continue;
581 distorted_subdomain_volumes[sub_id] += elem->volume();
582 highest_subdomain_id = std::max(sub_id, highest_subdomain_id);
583 }
584
585
587 for (
const auto sub_id :
make_range(highest_subdomain_id + 1)) {
588
589 if (distorted_subdomain_volumes.find(sub_id) ==
590 distorted_subdomain_volumes.end())
591 distorted_subdomain_volumes[sub_id] = 0.;
592
593 mesh.
comm().
sum(distorted_subdomain_volumes[sub_id]);
594 }
595
596
597
599 }
600
601
603 libmesh_error_msg_if(elem_orders.size() != 1,
604 "The variational smoother cannot be used for mixed-order meshes!");
605
606
607
608
609 const auto scale_factor = *elem_orders.begin() * ((type_is_pyramid || type_is_tet) ? 2 * 4 : 1);
610
611
612 auto node_distortion_is = [&n_elems_per_side, &
dim, &boundary_info, &scale_factor, &type_is_prism](
614
615 std::vector<boundary_id_type> boundary_ids;
616 boundary_info.boundary_ids(&node, boundary_ids);
617
618
619 const auto num_dofs =
dim - boundary_ids.size();
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639 std::size_t num_zero_or_one = 0;
640
641 bool distorted = false;
643 {
644 const Real r = node(d);
645 const Real R = r * n_elems_per_side * scale_factor;
646 CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, r);
647 CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, 1 - r);
648
649 bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
650 if (type_is_prism && (scale_factor == 3))
651 {
652
653
654
655 const Real R_prism = R / scale_factor * 2;
656 const bool d_distorted_prism =
657 std::abs(R_prism - std::round(R_prism)) > distortion_tol;
658 d_distorted &= d_distorted_prism;
659 }
660 distorted |= d_distorted;
662 }
663
664 CPPUNIT_ASSERT_GREATEREQUAL(
dim - num_dofs, num_zero_or_one);
665
666
667 if (num_dofs == 0)
668
669 return true;
670 return distorted == distortion;
671 };
672
673
674 for (
auto node :
mesh.node_ptr_range())
675 CPPUNIT_ASSERT(node_distortion_is(*node, true));
676
677
678
679
680 if (type_is_tri)
681 {
682
683
684
685 SquareToParallelogram stp;
687 }
688 else if (type_is_prism)
689 {
690
691
692
693 CubeToParallelepiped ctp;
695 }
696
698
699 if (type_is_tri)
700 {
701
702
703
704 ParallelogramToSquare pts;
706 }
707 else if (type_is_prism)
708 {
709
710
711
712 ParallelepipedToCube ptc;
714 }
715
716 if (multiple_subdomains)
717 {
718
719
720
721 std::unordered_map<subdomain_id_type, Real> smoothed_subdomain_volumes;
722 for (
auto *elem :
mesh.active_element_ptr_range()) {
723
724 if (elem->processor_id() != proc_id)
725 continue;
726 smoothed_subdomain_volumes[elem->subdomain_id()] += elem->volume();
727 }
728
729
730 for (
const auto sub_id :
make_range(highest_subdomain_id + 1)) {
731
732 if (smoothed_subdomain_volumes.find(sub_id) ==
733 smoothed_subdomain_volumes.end())
734 smoothed_subdomain_volumes[sub_id] = 0.;
735
736 mesh.
comm().
sum(smoothed_subdomain_volumes[sub_id]);
737 }
738
739 for (
const auto sub_id :
make_range(highest_subdomain_id + 1))
741 libmesh_map_find(distorted_subdomain_volumes, sub_id),
742 libmesh_map_find(smoothed_subdomain_volumes, sub_id),
TOLERANCE));
743 }
744 else
745 {
746 std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
748
749
750
751 std::set<dof_id_type> nodes_checked;
753
754 for (
const auto * elem :
mesh.active_element_ptr_range())
755 {
757 {
758 const auto & node = elem->node_ref(local_node_id);
759 if (nodes_checked.find(node.id()) != nodes_checked.end())
760 continue;
761
762 nodes_checked.insert(node.id());
763
764
765
766 if (type_is_pyramid)
767 {
768 if (local_node_id > 8 && local_node_id < 13)
769 {
770
771
772
773
774
775
776
777
778
779
780 const auto & base = elem->node_ref(local_node_id - 9);
781 const auto & apex = elem->node_ref(4);
782 const Real x = (type ==
PYRAMID18) ? 0.56646084 : 0.54985875;
783
784 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(base + x * (apex - base), tol));
785 continue;
786 }
787 else if (local_node_id > 13)
788 {
789
790
791
792
793
794
795
796
797 const auto & base1 = elem->node_ref(local_node_id - 14);
798 const auto & base2 = elem->node_ref((local_node_id - 13) % 4);
799 const auto & apex = elem->node_ref(4);
800
801 const auto node_approx =
802 (0.31401599 * base1 + 0.31401599 * base2 + 0.37196802 * apex);
803 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(node_approx, tol));
804 continue;
805 }
806 }
807
808
809
810 else if (type_is_tet && !elem->is_vertex(local_node_id))
811 {
812
813
814
815 std::vector<const Node *> neighbors;
817 mesh, node, nodes_to_elem_map, neighbors);
818
819 switch (neighbors.size())
820 {
821 case 2: {
822
823
824
825
826 const auto is_0_cube_center =
828 const auto is_1_cube_center =
831
832 if (is_0_cube_center || is_1_cube_center)
833 {
834 const auto & cube_center =
835 is_0_cube_center ? *neighbors[0] : *neighbors[1];
836 const auto & other =
837 is_0_cube_center ? *neighbors[1] : *neighbors[0];
838
839
840
841
843 {
844 const Real x = (type ==
TET10) ? 0.42895041 : 0.41486385;
845 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(
846 other + x * (cube_center - other), tol));
847 }
848
850 {
851 const Real x = (type ==
TET10) ? 0.55388920 : 0.58093516;
852 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(
853 other + x * (cube_center - other), tol));
854 }
855 }
856
857 else
858 {
859
860
861
862 const auto is_0_cube_vertex =
864 const auto is_1_cube_vertex =
866 const auto is_0_cube_face_center =
868 const auto is_1_cube_face_center =
870
871 if (is_0_cube_vertex && is_1_cube_vertex)
872
873
874
875 CPPUNIT_ASSERT(node_distortion_is(node, false));
876
877 else
878 {
879 libmesh_error_msg_if(
880 (is_0_cube_center || is_0_cube_face_center) &&
881 (is_1_cube_center || is_1_cube_face_center),
882 "We should never get here!");
883 const auto & cube_vertex =
884 is_0_cube_vertex ? *neighbors[0] : *neighbors[1];
885 const auto & cube_face_center =
886 is_0_cube_face_center ? *neighbors[0] : *neighbors[1];
887 const Real x = (type ==
TET10) ? 0.61299101 : 0.65125580;
888 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(
889 cube_vertex + x * (cube_face_center - cube_vertex), tol));
890 }
891 }
892
893 continue;
894 break;
895 }
896
897 case 6: {
898
899
900 neighbors.erase(std::remove_if(neighbors.begin(),
901 neighbors.end(),
902 [&elem](
const Node * n) {
903 return !elem->is_vertex(
904 elem->local_node(n->id()));
905 }),
906 neighbors.end());
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921 unsigned int num_vertices_at_cube_center = 0;
922 unsigned int num_vertices_at_cube_vertices = 0;
923 unsigned int num_vertices_at_cube_face_centers = 0;
924 for (const auto * neighbor : neighbors)
925 {
927 num_vertices_at_cube_center += 1;
929 num_vertices_at_cube_vertices += 1;
931 num_vertices_at_cube_face_centers += 1;
932 }
933
934
935
937
938
939 if (num_vertices_at_cube_vertices == 2)
940 {
941
942 if (num_vertices_at_cube_center)
943 for (const auto * neighbor : neighbors)
944 {
950 else
951 libmesh_error_msg("We should never get here!");
952
953 node_approx +=
weight * (*neighbor);
954 }
955
956
957 else if (num_vertices_at_cube_face_centers)
958 for (const auto * neighbor : neighbors)
959 {
965 else
966 libmesh_error_msg("We should never get here!");
967
968 node_approx +=
weight * (*neighbor);
969 }
970
971 else
972 libmesh_error_msg("We should never get here!");
973 }
974
975
976 else if (num_vertices_at_cube_center && num_vertices_at_cube_vertices &&
977 num_vertices_at_cube_face_centers)
978 for (const auto * neighbor : neighbors)
979 {
987 else
988 libmesh_error_msg("We should never get here!");
989
990 node_approx +=
weight * (*neighbor);
991 }
992
993 else
994 libmesh_error_msg("We should never get here!");
995
996 CPPUNIT_ASSERT(node.absolute_fuzzy_equals(node_approx, tol));
997
998 continue;
999 break;
1000 }
1001 default: {
1002 libmesh_error_msg(node << " has unexpected number of neighbors ("
1003 << neighbors.size() << ")");
1004 break;
1005 }
1006 }
1007 }
1008
1009 CPPUNIT_ASSERT(node_distortion_is(node, false));
1010
1011 }
1012 }
1013 }
1014 }
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
