libMesh
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | List of all members
DisjointNeighborTest Class Reference
Inheritance diagram for DisjointNeighborTest:
[legend]

Public Member Functions

 LIBMESH_CPPUNIT_TEST_SUITE (DisjointNeighborTest)
 
 CPPUNIT_TEST (testTempJump)
 
 CPPUNIT_TEST (testTempJumpRefine)
 
 CPPUNIT_TEST (testTempJumpLocalRefineFail)
 
 CPPUNIT_TEST (testCloneEquality)
 
 CPPUNIT_TEST (testStitchingDiscontinuousBoundaries)
 
 CPPUNIT_TEST (testPreserveDisjointNeighborPairsAfterStitch)
 
 CPPUNIT_TEST (testStitchCrossMesh)
 
 CPPUNIT_TEST (testDisjointNeighborConflictError)
 
 CPPUNIT_TEST_SUITE_END ()
 

Private Member Functions

void testTempJump ()
 
void testCloneEquality ()
 
void testStitchingDiscontinuousBoundaries ()
 
void testTempJumpRefine ()
 
void testTempJumpLocalRefineFail ()
 
void testPreserveDisjointNeighborPairsAfterStitch ()
 
void testStitchCrossMesh ()
 
void testDisjointNeighborConflictError ()
 
void build_two_disjoint_elems (Mesh &mesh, bool add_pair=true, const Point &translate_vec=Point(0.0, 0.0, 0.0))
 
void build_split_mesh_with_interface (Mesh &mesh, unsigned nx, unsigned ny, bool test_local_refinement=false)
 
void build_four_disjoint_elems (Mesh &mesh)
 

Detailed Description

Definition at line 193 of file disjoint_neighbor_test.C.

Member Function Documentation

◆ build_four_disjoint_elems()

void DisjointNeighborTest::build_four_disjoint_elems ( Mesh mesh)
inlineprivate

Definition at line 842 of file disjoint_neighbor_test.C.

843{
844 // Clear any existing data and initialize the mesh
845 mesh.clear();
847
848 BoundaryInfo &boundary = mesh.get_boundary_info();
849
850 // --------------------------------------------------------------------
851 // Domain layout:
852 // 4 elements aligned in a single row, each of width = 0.25
853 //
854 // y=1: x0---x1_L---x1_R---x2_L---x2_R---x3_L---x3_R---x4
855 //
856 // y=0: x0---x1_L---x1_R---x2_L---x2_R---x3_L---x3_R---x4
857 //
858 // The interfaces between adjacent elements are duplicated:
859 // for each internal interface, left and right sides have distinct nodes
860 // (geometrically coincident but topologically disconnected).
861 // --------------------------------------------------------------------
862 const Real dx = 0.25;
863
864 // ---- Add nodal coordinates ----
865 // Each interface has duplicated nodes (left and right) to represent discontinuity.
866 unsigned id = 0;
867 for (unsigned i = 0; i < 4; ++i)
868 {
869 Real xL = i * dx;
870 Real xR = (i + 1) * dx;
871
872 // Left nodes of the element
873 mesh.add_point(Point(xL, 0.0), id++); // bottom-left
874 mesh.add_point(Point(xL, 1.0), id++); // top-left
875
876 // Right nodes (duplicated interface nodes)
877 mesh.add_point(Point(xR, 0.0), id++); // bottom-right
878 mesh.add_point(Point(xR, 1.0), id++); // top-right
879 }
880
881 // Add one final pair of rightmost boundary nodes
882 mesh.add_point(Point(1.0, 0.0), id++);
883 mesh.add_point(Point(1.0, 1.0), id++);
884
885 // ---- Create four QUAD4 elements ----
886 // Node ordering: 0 = BL, 1 = BR, 2 = TR, 3 = TL
887 unsigned elem_id = 0;
888 for (unsigned i = 0; i < 4; ++i)
889 {
890 Elem *elem = mesh.add_elem(Elem::build_with_id(QUAD4, elem_id++));
891
892 unsigned base = i * 4;
893 elem->set_node(0, mesh.node_ptr(base)); // bottom-left
894 elem->set_node(1, mesh.node_ptr(base + 2)); // bottom-right
895 elem->set_node(2, mesh.node_ptr(base + 3)); // top-right
896 elem->set_node(3, mesh.node_ptr(base + 1)); // top-left
897
898 // Add outer boundaries
899 if (i == 0)
900 boundary.add_side(elem, 3, left_id); // left outer boundary
901 if (i == 3)
902 boundary.add_side(elem, 1, right_id); // right outer boundary
903 }
904
905 // ---- Define internal interface sides ----
906 // Element side 1 (right) connects to the next element’s side 3 (left)
907 // Each interface is treated as a disjoint (non-conforming) boundary pair
908 boundary.add_side(mesh.elem_ptr(0), 1, interface1_left_id);
909 boundary.add_side(mesh.elem_ptr(1), 3, interface1_right_id);
910
911 boundary.add_side(mesh.elem_ptr(1), 1, interface2_left_id);
912 boundary.add_side(mesh.elem_ptr(2), 3, interface2_right_id);
913
914 boundary.add_side(mesh.elem_ptr(2), 1, interface3_left_id);
915 boundary.add_side(mesh.elem_ptr(3), 3, interface3_right_id);
916
917 // ---- Register disjoint neighbor boundary pairs ----
918 // These pairs inform libMesh that the specified sides are geometrically
919 // adjacent but topologically disconnected.
923
924 // ---- Assign descriptive names to boundary IDs ----
925 boundary.sideset_name(left_id) = "left_boundary";
926 boundary.sideset_name(right_id) = "right_boundary";
927
928 boundary.sideset_name(interface1_left_id) = "interface1_left";
929 boundary.sideset_name(interface1_right_id) = "interface1_right";
930 boundary.sideset_name(interface2_left_id) = "interface2_left";
931 boundary.sideset_name(interface2_right_id) = "interface2_right";
932 boundary.sideset_name(interface3_left_id) = "interface3_left";
933 boundary.sideset_name(interface3_right_id) = "interface3_right";
934
935 // ---- Finalize mesh setup ----
936 mesh.allow_renumbering(false);
938}
The BoundaryInfo class contains information relevant to boundary conditions including storing faces,...
std::string & sideset_name(boundary_id_type id)
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure.
This is the base class from which all geometric element types are derived.
Definition elem.h:96
virtual Node *& set_node(const unsigned int i)
Definition elem.h:2567
static std::unique_ptr< Elem > build_with_id(const ElemType type, dof_id_type id)
Calls the build() method above with a nullptr parent, and additionally sets the newly-created Elem's ...
Definition elem.C:556
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition mesh_base.h:170
virtual const Node * node_ptr(const dof_id_type i) const =0
void allow_renumbering(bool allow)
If false is passed in then this mesh will no longer be renumbered when being prepared for use.
Definition mesh_base.h:1355
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly created (or read) mesh for use.
Definition mesh_base.C:824
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition mesh_base.h:423
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition mesh_base.C:1036
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
void add_disjoint_neighbor_boundary_pairs(const boundary_id_type b1, const boundary_id_type b2, const RealVectorValue &translation)
Register a pair of boundaries as disjoint neighbor boundary pairs.
Definition mesh_base.C:2272
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
static const boundary_id_type left_id
static const boundary_id_type interface1_left_id
static const boundary_id_type interface2_left_id
static const boundary_id_type interface3_left_id
static const boundary_id_type right_id
static const boundary_id_type interface3_right_id
static const boundary_id_type interface1_right_id
static const boundary_id_type interface2_right_id
MeshBase & mesh
VectorValue< Real > RealVectorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

References libMesh::MeshBase::add_disjoint_neighbor_boundary_pairs(), libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::MeshBase::allow_renumbering(), libMesh::Elem::build_with_id(), libMesh::MeshBase::clear(), libMesh::MeshBase::elem_ptr(), libMesh::MeshBase::get_boundary_info(), interface1_left_id, interface1_right_id, interface2_left_id, interface2_right_id, interface3_left_id, interface3_right_id, left_id, mesh, libMesh::MeshBase::node_ptr(), libMesh::MeshBase::prepare_for_use(), libMesh::QUAD4, libMesh::Real, right_id, libMesh::MeshBase::set_mesh_dimension(), libMesh::Elem::set_node(), and libMesh::BoundaryInfo::sideset_name().

Referenced by testPreserveDisjointNeighborPairsAfterStitch().

◆ build_split_mesh_with_interface()

void DisjointNeighborTest::build_split_mesh_with_interface ( Mesh mesh,
unsigned  nx,
unsigned  ny,
bool  test_local_refinement = false 
)
inlineprivate

Definition at line 716 of file disjoint_neighbor_test.C.

717 {
718 // Ensure nx is even so the interface aligns with element boundaries
719 libmesh_error_msg_if(nx % 2 != 0, "nx must be even!");
720
721 mesh.clear();
723
724 const unsigned nxL = nx / 2; // Number of elements in x-direction (left half)
725 const unsigned nxR = nx / 2; // Number of elements in x-direction (right half)
726 const double dx = 1.0 / static_cast<double>(nx);
727 const double dy = 1.0 / static_cast<double>(ny);
728
729 // --- Generate points for the left subdomain [0, 0.5] ---
730 // Each row of the left subdomain has (nxL + 1) nodes.
731 // Total nodes = (nxL + 1) * (ny + 1).
732 auto nidL = [&](unsigned i, unsigned j) {
733 return static_cast<dof_id_type>(j * (nxL + 1) + i);
734 };
735
736 for (unsigned j = 0; j <= ny; ++j)
737 {
738 const double y = j * dy;
739 for (unsigned i = 0; i <= nxL; ++i)
740 {
741 const double x = i * dx; // At i = nxL, x = 0.5 (interface)
742 mesh.add_point(Point(x, y), nidL(i, j));
743 }
744 }
745
746 // --- Generate points for the right subdomain [0.5, 1.0] ---
747 // Interface nodes at x = 0.5 are duplicated with new node IDs.
748 const dof_id_type baseR = static_cast<dof_id_type>((nxL + 1) * (ny + 1));
749
750 auto nidR = [&](unsigned i, unsigned j) {
751 return static_cast<dof_id_type>(baseR + j * (nxR + 1) + i);
752 };
753
754 for (unsigned j = 0; j <= ny; ++j)
755 {
756 const double y = j * dy;
757 for (unsigned i = 0; i <= nxR; ++i)
758 {
759 const double x = 0.5 + i * dx; // At i = 0, x = 0.5 (same coords as left interface, different ID)
760 mesh.add_point(Point(x, y), nidR(i, j));
761 }
762 }
763
764 BoundaryInfo &boundary = mesh.get_boundary_info();
765
766 // --- Create left subdomain elements (QUAD4) ---
767 // Node order: 0 = BL, 1 = BR, 2 = TR, 3 = TL
768 // Side order: 0 = bottom(0-1), 1 = right(1-2), 2 = top(2-3), 3 = left(3-0)
769 dof_id_type eid = 0;
770 for (unsigned j = 0; j < ny; ++j)
771 {
772 for (unsigned i = 0; i < nxL; ++i)
773 {
775 e->set_node(0, mesh.node_ptr(nidL(i, j)));
776 e->set_node(1, mesh.node_ptr(nidL(i+1, j)));
777 e->set_node(2, mesh.node_ptr(nidL(i+1, j+1)));
778 e->set_node(3, mesh.node_ptr(nidL(i, j+1)));
779
780 // Left outer boundary (i == 0 -> side 3)
781 if (i == 0)
782 boundary.add_side(e, 3, left_id);
783
784 // Interface (left side) boundary (i == nxL - 1 -> side 1)
785 if (i == nxL - 1)
786 boundary.add_side(e, 1, interface_left_id);
787 }
788 }
789
790 // --- Create right subdomain elements (QUAD4) ---
791 for (unsigned j = 0; j < ny; ++j)
792 {
793 for (unsigned i = 0; i < nxR; ++i)
794 {
796 e->set_node(0, mesh.node_ptr(nidR(i, j)));
797 e->set_node(1, mesh.node_ptr(nidR(i+1, j)));
798 e->set_node(2, mesh.node_ptr(nidR(i+1, j+1)));
799 e->set_node(3, mesh.node_ptr(nidR(i, j+1)));
800
801 // Interface (right side) boundary (i == 0 -> side 3)
802 if (i == 0)
803 boundary.add_side(e, 3, interface_right_id);
804
805 // Right outer boundary (i == nxR - 1 -> side 1)
806 if (i == nxR - 1)
807 boundary.add_side(e, 1, right_id);
808 }
809 }
810
811 // --- Assign human-readable boundary names ---
812 boundary.sideset_name(left_id) = "left_boundary";
813 boundary.sideset_name(right_id) = "right_boundary";
814 boundary.sideset_name(interface_left_id) = "interface_left";
815 boundary.sideset_name(interface_right_id) = "interface_right";
816
817
818 // This is the key testing step: inform libMesh about the disjoint neighbor boundary pairs
819 // And, in `prepare_for_use()`, libMesh will set up the disjoint neighbor relationships.
821
823
824#if LIBMESH_ENABLE_AMR
825 if (test_local_refinement)
826 {
827 // set refinement flags
828 for (auto & elem : mesh.element_ptr_range())
829 {
830 const Real xmid = elem->vertex_average()(0);
831 if ((xmid < 0.5 && xmid > 0.5 - dx) ||
832 (xmid > 0.5 && xmid < 0.5 + dx))
833 elem->set_refinement_flag(Elem::REFINE);
834 }
835
836 // refine
838 }
839#endif
840 }
Implements (adaptive) mesh refinement algorithms for a MeshBase.
bool refine_elements()
Only refines the user-requested elements.
static const boundary_id_type interface_left_id
static const boundary_id_type interface_right_id
uint8_t dof_id_type
Definition id_types.h:67

References libMesh::MeshBase::add_disjoint_neighbor_boundary_pairs(), libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::Elem::build_with_id(), libMesh::MeshBase::clear(), libMesh::MeshBase::get_boundary_info(), interface_left_id, interface_right_id, left_id, mesh, libMesh::MeshBase::node_ptr(), libMesh::MeshBase::prepare_for_use(), libMesh::QUAD4, libMesh::Real, libMesh::Elem::REFINE, libMesh::MeshRefinement::refine_elements(), right_id, libMesh::MeshBase::set_mesh_dimension(), libMesh::Elem::set_node(), and libMesh::BoundaryInfo::sideset_name().

Referenced by testTempJumpLocalRefineFail(), and testTempJumpRefine().

◆ build_two_disjoint_elems()

void DisjointNeighborTest::build_two_disjoint_elems ( Mesh mesh,
bool  add_pair = true,
const Point translate_vec = Point(0.0, 0.0, 0.0) 
)
inlineprivate

Definition at line 639 of file disjoint_neighbor_test.C.

640 {
641 // Domain: x in (0, 1), y in (0, 1)
642 // Split into two subdomains:
643 // Left subdomain: 0 <= x <= 0.5
644 // Right subdomain: 0.5 <= x <= 1
645 //
646 // Note: Points at x = 0.5 are duplicated (same coordinates but different node IDs)
647 // to represent an interface or discontinuity between the two subdomains.
648 //
649 // Coordinates layout:
650 //
651 // (0,1) (0.5,1)_L (0.5,1)_R (1,1)
652 // x--------x x--------x
653 // | | | |
654 // | Left | Interface | Right |
655 // | | | |
656 // x--------x x--------x
657 // (0,0) (0.5,0)_L (0.5,0)_R (1,0)
658
659
660 // ---- Define points ----
661
662 // Left subdomain nodes
663 mesh.add_point(Point(0.0, 0.0) + translate_vec, 0); // bottom-left corner
664 mesh.add_point(Point(0.5, 0.0) + translate_vec, 1); // bottom-right corner of left element (interface node)
665 mesh.add_point(Point(0.0, 1.0) + translate_vec, 2); // top-left corner
666 mesh.add_point(Point(0.5, 1.0) + translate_vec, 3); // top-right corner of left element (interface node)
667
668 // Right subdomain nodes (duplicated interface points)
669 mesh.add_point(Point(0.5, 0.0) + translate_vec, 4); // bottom-left corner of right element (same coords as node 1) (interface node)
670 mesh.add_point(Point(1.0, 0.0) + translate_vec, 5); // bottom-right corner
671 mesh.add_point(Point(0.5, 1.0) + translate_vec, 6); // top-left corner of right element (same coords as node 3) (interface node)
672 mesh.add_point(Point(1.0, 1.0) + translate_vec, 7); // top-right corner
673
674
675 // ---- Define elements ----
676 BoundaryInfo & boundary = mesh.get_boundary_info();
677 // Left element (element ID = 0)
678 {
680 elem->set_node(0, mesh.node_ptr(0)); // bottom-left (0,0)
681 elem->set_node(1, mesh.node_ptr(1)); // bottom-right (0.5,0)
682 elem->set_node(2, mesh.node_ptr(3)); // top-right (0.5,1)
683 elem->set_node(3, mesh.node_ptr(2)); // top-left (0,1)
684 boundary.add_side(elem, 3, left_id); // left boundary
685 boundary.add_side(elem, 1, interface_left_id);
686 boundary.sideset_name(left_id) = "left_boundary";
687 boundary.sideset_name(interface_left_id) = "interface_left";
688 }
689
690 // Right element (element ID = 1)
691 {
693 elem->set_node(0, mesh.node_ptr(4)); // bottom-left (0.5,0)_R
694 elem->set_node(1, mesh.node_ptr(5)); // bottom-right (1,0)
695 elem->set_node(2, mesh.node_ptr(7)); // top-right (1,1)
696 elem->set_node(3, mesh.node_ptr(6)); // top-left (0.5,1)_R
697 boundary.add_side(elem, 1, right_id); // right boundary
698 boundary.add_side(elem, 3, interface_right_id);
699 boundary.sideset_name(right_id) = "right_boundary";
700 boundary.sideset_name(interface_right_id) = "interface_right";
701 }
702
703 // This is the key testing step: inform libMesh about the disjoint neighbor boundary pairs
704 // And, in `prepare_for_use()`, libMesh will set up the disjoint neighbor relationships.
705 if (add_pair)
707
708 // libMesh shouldn't renumber, or our based-on-initial-id
709 // assertions later may fail.
710 mesh.allow_renumbering(false);
711
713 }

References libMesh::MeshBase::add_disjoint_neighbor_boundary_pairs(), libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::MeshBase::allow_renumbering(), libMesh::Elem::build_with_id(), libMesh::MeshBase::get_boundary_info(), interface_left_id, interface_right_id, left_id, mesh, libMesh::MeshBase::node_ptr(), libMesh::MeshBase::prepare_for_use(), libMesh::QUAD4, right_id, libMesh::Elem::set_node(), and libMesh::BoundaryInfo::sideset_name().

Referenced by testCloneEquality(), testDisjointNeighborConflictError(), testStitchCrossMesh(), testStitchingDiscontinuousBoundaries(), and testTempJump().

◆ CPPUNIT_TEST() [1/8]

DisjointNeighborTest::CPPUNIT_TEST ( testCloneEquality  )

◆ CPPUNIT_TEST() [2/8]

DisjointNeighborTest::CPPUNIT_TEST ( testDisjointNeighborConflictError  )

◆ CPPUNIT_TEST() [3/8]

DisjointNeighborTest::CPPUNIT_TEST ( testPreserveDisjointNeighborPairsAfterStitch  )

◆ CPPUNIT_TEST() [4/8]

DisjointNeighborTest::CPPUNIT_TEST ( testStitchCrossMesh  )

◆ CPPUNIT_TEST() [5/8]

DisjointNeighborTest::CPPUNIT_TEST ( testStitchingDiscontinuousBoundaries  )

◆ CPPUNIT_TEST() [6/8]

DisjointNeighborTest::CPPUNIT_TEST ( testTempJump  )

◆ CPPUNIT_TEST() [7/8]

DisjointNeighborTest::CPPUNIT_TEST ( testTempJumpLocalRefineFail  )

◆ CPPUNIT_TEST() [8/8]

DisjointNeighborTest::CPPUNIT_TEST ( testTempJumpRefine  )

◆ CPPUNIT_TEST_SUITE_END()

DisjointNeighborTest::CPPUNIT_TEST_SUITE_END ( )

◆ LIBMESH_CPPUNIT_TEST_SUITE()

DisjointNeighborTest::LIBMESH_CPPUNIT_TEST_SUITE ( DisjointNeighborTest  )

◆ testCloneEquality()

void DisjointNeighborTest::testCloneEquality ( )
inlineprivate

Definition at line 298 of file disjoint_neighbor_test.C.

299 {
300 LOG_UNIT_TEST;
301
302 Mesh mesh(*TestCommWorld, 2);
304
305 std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
306 CPPUNIT_ASSERT(*mesh_clone == mesh);
307 }
void build_two_disjoint_elems(Mesh &mesh, bool add_pair=true, const Point &translate_vec=Point(0.0, 0.0, 0.0))
virtual std::unique_ptr< MeshBase > clone() const =0
Virtual "copy constructor".
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition mesh.h:51

References build_two_disjoint_elems(), libMesh::MeshBase::clone(), mesh, and TestCommWorld.

◆ testDisjointNeighborConflictError()

void DisjointNeighborTest::testDisjointNeighborConflictError ( )
inlineprivate

Definition at line 552 of file disjoint_neighbor_test.C.

553 {
554 LOG_UNIT_TEST;
555
556 try
557 {
558 Mesh mesh_Y(*TestCommWorld, 2);
559 build_two_disjoint_elems(mesh_Y, false);
560 mesh_Y.add_disjoint_neighbor_boundary_pairs(right_id, left_id, RealVectorValue(-1.0, 0.0, 0.0));
561
562 // Create another mesh that defines the same pair properly
563 auto right_id_x = right_id + 10;
564 auto interface_left_id_x = interface_left_id + 10;
565 auto interface_right_id_x = interface_right_id + 10;
566
567 Mesh mesh_X(*TestCommWorld, 2);
568 const auto translate_vec = Point(1.0, 0.0, 0.0);
569 // Left subdomain nodes
570 mesh_X.add_point(Point(0.0, 0.0) + translate_vec, 0); // bottom-left corner
571 mesh_X.add_point(Point(0.5, 0.0) + translate_vec, 1); // bottom-right corner of left element (interface node)
572 mesh_X.add_point(Point(0.0, 1.0) + translate_vec, 2); // top-left corner
573 mesh_X.add_point(Point(0.5, 1.0) + translate_vec, 3); // top-right corner of left element (interface node)
574
575 // Right subdomain nodes (duplicated interface points)
576 mesh_X.add_point(Point(0.5, 0.0) + translate_vec, 4); // bottom-left corner of right element (same coords as node 1) (interface node)
577 mesh_X.add_point(Point(1.0, 0.0) + translate_vec, 5); // bottom-right corner
578 mesh_X.add_point(Point(0.5, 1.0) + translate_vec, 6); // top-left corner of right element (same coords as node 3) (interface node)
579 mesh_X.add_point(Point(1.0, 1.0) + translate_vec, 7); // top-right corner
580
581
582 // ---- Define elements ----
583 BoundaryInfo & boundary = mesh_X.get_boundary_info();
584 // Left element (element ID = 0)
585 {
586 Elem * elem = mesh_X.add_elem(Elem::build_with_id(QUAD4, 0));
587 elem->set_node(0, mesh_X.node_ptr(0)); // bottom-left (0,0)
588 elem->set_node(1, mesh_X.node_ptr(1)); // bottom-right (0.5,0)
589 elem->set_node(2, mesh_X.node_ptr(3)); // top-right (0.5,1)
590 elem->set_node(3, mesh_X.node_ptr(2)); // top-left (0,1)
591 boundary.add_side(elem, 3, right_id); // left boundary -> use the same ID as mesh_Y to trigger conflict
592 boundary.add_side(elem, 1, interface_left_id_x);
593 boundary.sideset_name(right_id) = "left_boundary_x";
594 boundary.sideset_name(interface_left_id_x) = "interface_left_x";
595 }
596
597 // Right element (element ID = 1)
598 {
599 Elem * elem = mesh_X.add_elem(Elem::build_with_id(QUAD4, 1));
600 elem->set_node(0, mesh_X.node_ptr(4)); // bottom-left (0.5,0)_R
601 elem->set_node(1, mesh_X.node_ptr(5)); // bottom-right (1,0)
602 elem->set_node(2, mesh_X.node_ptr(7)); // top-right (1,1)
603 elem->set_node(3, mesh_X.node_ptr(6)); // top-left (0.5,1)_R
604 boundary.add_side(elem, 1, right_id_x); // right boundary
605 boundary.add_side(elem, 3, interface_right_id_x);
606 boundary.sideset_name(right_id_x) = "right_boundary_x";
607 boundary.sideset_name(interface_right_id_x) = "interface_right_x";
608 }
609
610 mesh_X.add_disjoint_neighbor_boundary_pairs(right_id, right_id_x, RealVectorValue(1.0, 0.0, 0.0));
611
612 // libMesh shouldn't renumber, or our based-on-initial-id
613 // assertions later may fail.
614 mesh_X.allow_renumbering(false);
615
616 mesh_X.prepare_for_use();
617
618 // Stitching should trigger an error since Y already has the same
619 // boundary IDs but not registered as a pair
620 mesh_Y.stitch_meshes(mesh_X,
621 right_id, right_id_x,
622 1e-8,
623 true, // clear_stitched_boundary_ids
624 false, // verbose
625 false, // use_binary_search
626 false, // enforce_all_nodes_match_on_boundaries
627 false, // merge_boundary_nodes_all_or_nothing
628 false); // prepare_after_stitching
629 }
630 catch (const std::exception &e)
631 {
632 CPPUNIT_ASSERT(std::string(e.what()).find("Disjoint neighbor boundary pairing mismatch") != std::string::npos);
633 }
634 }

References libMesh::MeshBase::add_disjoint_neighbor_boundary_pairs(), libMesh::DistributedMesh::add_elem(), libMesh::DistributedMesh::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::MeshBase::allow_renumbering(), build_two_disjoint_elems(), libMesh::Elem::build_with_id(), libMesh::MeshBase::get_boundary_info(), interface_left_id, interface_right_id, left_id, libMesh::DistributedMesh::node_ptr(), libMesh::MeshBase::prepare_for_use(), libMesh::QUAD4, right_id, libMesh::Elem::set_node(), libMesh::BoundaryInfo::sideset_name(), libMesh::UnstructuredMesh::stitch_meshes(), and TestCommWorld.

◆ testPreserveDisjointNeighborPairsAfterStitch()

void DisjointNeighborTest::testPreserveDisjointNeighborPairsAfterStitch ( )
inlineprivate

Definition at line 410 of file disjoint_neighbor_test.C.

411 {
412 LOG_UNIT_TEST;
413 Mesh mesh(*TestCommWorld, 2);
415 // ExodusII_IO(mesh).write("four_elem_one_row.e");
416 CPPUNIT_ASSERT(mesh.parallel_n_nodes() == 16);
417
418 mesh.stitch_surfaces(interface2_left_id,
420 1e-8 /* tolerance */,
421 true /* clear_stitched_boundary_ids */,
422 false /* verbose */,
423 false /* use_binary_search */,
424 false /* enforce_all_nodes_match_on_boundaries */,
425 false /* merge_boundary_nodes_all_or_nothing */,
426 true /* prepare_after_stitching */);
427
428 // ExodusII_IO(mesh).write("four_elem_one_row_after_stitching.e");
429 CPPUNIT_ASSERT(mesh.parallel_n_nodes() == 14);
430 const auto * disjoint_pairs = mesh.get_disjoint_neighbor_boundary_pairs();
431
432 // make sure that both the interface1 and the interface3 pairing still exist afterward.
433 CPPUNIT_ASSERT(disjoint_pairs->size() == 4);
434 for (const auto & [bid, pb_ptr] : *disjoint_pairs)
435 {
436 const auto & pb = *pb_ptr;
437 const boundary_id_type a = pb.myboundary;
438 const boundary_id_type b = pb.pairedboundary;
439 if ((a == interface1_left_id && b == interface1_right_id) ||
443 continue;
444 else
445 CPPUNIT_FAIL("Disjoint neighbor boundary pair missing after stitching!");
446 }
447
448 }
void build_four_disjoint_elems(Mesh &mesh)
PeriodicBoundaries * get_disjoint_neighbor_boundary_pairs()
Definition mesh_base.C:2296
virtual dof_id_type parallel_n_nodes() const =0
static const Real b
int8_t boundary_id_type
Definition id_types.h:51

References b, build_four_disjoint_elems(), libMesh::MeshBase::get_disjoint_neighbor_boundary_pairs(), interface1_left_id, interface1_right_id, interface2_left_id, interface2_right_id, interface3_left_id, interface3_right_id, mesh, libMesh::MeshBase::parallel_n_nodes(), and TestCommWorld.

◆ testStitchCrossMesh()

void DisjointNeighborTest::testStitchCrossMesh ( )
inlineprivate

Definition at line 450 of file disjoint_neighbor_test.C.

451 {
452 LOG_UNIT_TEST;
453
454 // Create a mesh with interface boundaries but without pairing
455 Mesh mesh_Y(*TestCommWorld, 2);
457
458 // Create another mesh that defines the same pair properly
459 auto left_id_x = left_id + 10; // to avoid boundary ID conflict with mesh_Y
460 auto right_id_x = right_id + 10;
461 auto interface_left_id_x = interface_left_id + 10;
462 auto interface_right_id_x = interface_right_id + 10;
463
464 Mesh mesh_X(*TestCommWorld, 2);
465 const auto translate_vec = Point(1.0, 0.0, 0.0);
466 // Left subdomain nodes
467 mesh_X.add_point(Point(0.0, 0.0) + translate_vec, 0); // bottom-left corner
468 mesh_X.add_point(Point(0.5, 0.0) + translate_vec, 1); // bottom-right corner of left element (interface node)
469 mesh_X.add_point(Point(0.0, 1.0) + translate_vec, 2); // top-left corner
470 mesh_X.add_point(Point(0.5, 1.0) + translate_vec, 3); // top-right corner of left element (interface node)
471
472 // Right subdomain nodes (duplicated interface points)
473 mesh_X.add_point(Point(0.5, 0.0) + translate_vec, 4); // bottom-left corner of right element (same coords as node 1) (interface node)
474 mesh_X.add_point(Point(1.0, 0.0) + translate_vec, 5); // bottom-right corner
475 mesh_X.add_point(Point(0.5, 1.0) + translate_vec, 6); // top-left corner of right element (same coords as node 3) (interface node)
476 mesh_X.add_point(Point(1.0, 1.0) + translate_vec, 7); // top-right corner
477
478
479 // ---- Define elements ----
480 BoundaryInfo & boundary = mesh_X.get_boundary_info();
481 // Left element (element ID = 0)
482 {
483 Elem * elem = mesh_X.add_elem(Elem::build_with_id(QUAD4, 0));
484 elem->set_node(0, mesh_X.node_ptr(0)); // bottom-left (0,0)
485 elem->set_node(1, mesh_X.node_ptr(1)); // bottom-right (0.5,0)
486 elem->set_node(2, mesh_X.node_ptr(3)); // top-right (0.5,1)
487 elem->set_node(3, mesh_X.node_ptr(2)); // top-left (0,1)
488 boundary.add_side(elem, 3, left_id_x); // left boundary
489 boundary.add_side(elem, 1, interface_left_id_x);
490 boundary.sideset_name(left_id_x) = "left_boundary_x";
491 boundary.sideset_name(interface_left_id_x) = "interface_left_x";
492 }
493
494 // Right element (element ID = 1)
495 {
496 Elem * elem = mesh_X.add_elem(Elem::build_with_id(QUAD4, 1));
497 elem->set_node(0, mesh_X.node_ptr(4)); // bottom-left (0.5,0)_R
498 elem->set_node(1, mesh_X.node_ptr(5)); // bottom-right (1,0)
499 elem->set_node(2, mesh_X.node_ptr(7)); // top-right (1,1)
500 elem->set_node(3, mesh_X.node_ptr(6)); // top-left (0.5,1)_R
501 boundary.add_side(elem, 1, right_id_x); // right boundary
502 boundary.add_side(elem, 3, interface_right_id_x);
503 boundary.sideset_name(right_id_x) = "right_boundary_x";
504 boundary.sideset_name(interface_right_id_x) = "interface_right_x";
505 }
506
507 mesh_X.add_disjoint_neighbor_boundary_pairs(interface_left_id_x, interface_right_id_x, RealVectorValue(0.0, 0.0, 0.0));
508
509 // libMesh shouldn't renumber, or our based-on-initial-id
510 // assertions later may fail.
511 mesh_X.allow_renumbering(false);
512
513 mesh_X.prepare_for_use();
514
515 // Stitching should trigger an error since Y already has the same
516 // boundary IDs but not registered as a pair
517 mesh_Y.stitch_meshes(mesh_X,
518 right_id, left_id_x,
519 1e-8,
520 true, // clear_stitched_boundary_ids
521 false, // verbose
522 false, // use_binary_search
523 false, // enforce_all_nodes_match_on_boundaries
524 false, // merge_boundary_nodes_all_or_nothing
525 false); // prepare_after_stitching
526
527 CPPUNIT_ASSERT(mesh_Y.parallel_n_nodes() == 14);
528 const auto * disjoint_pairs = mesh_Y.get_disjoint_neighbor_boundary_pairs();
529
530 // ExodusII_IO(mesh_Y).write("mesh_Y_after_stitching.e");
531 // make sure that both the interface and the interface_x pairing still exist afterward.
532 CPPUNIT_ASSERT(disjoint_pairs->size() == 4);
533 for (const auto & [bid, pb_ptr] : *disjoint_pairs)
534 {
535 const auto & pb = *pb_ptr;
536 const boundary_id_type a = pb.myboundary;
537 const boundary_id_type b = pb.pairedboundary;
538 if ((a == interface_left_id_x && b == interface_right_id_x) ||
539 (a == interface_right_id_x && b == interface_left_id_x) ||
542 continue;
543 else
544 CPPUNIT_FAIL("Disjoint neighbor boundary pair missing after stitching!");
545 }
546
547
548 }

References libMesh::MeshBase::add_disjoint_neighbor_boundary_pairs(), libMesh::DistributedMesh::add_elem(), libMesh::DistributedMesh::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::MeshBase::allow_renumbering(), b, build_two_disjoint_elems(), libMesh::Elem::build_with_id(), libMesh::MeshBase::get_boundary_info(), libMesh::MeshBase::get_disjoint_neighbor_boundary_pairs(), interface_left_id, interface_right_id, left_id, libMesh::DistributedMesh::node_ptr(), libMesh::DistributedMesh::parallel_n_nodes(), libMesh::MeshBase::prepare_for_use(), libMesh::QUAD4, right_id, libMesh::Elem::set_node(), libMesh::BoundaryInfo::sideset_name(), libMesh::UnstructuredMesh::stitch_meshes(), and TestCommWorld.

◆ testStitchingDiscontinuousBoundaries()

void DisjointNeighborTest::testStitchingDiscontinuousBoundaries ( )
inlineprivate

Definition at line 309 of file disjoint_neighbor_test.C.

310 {
311 LOG_UNIT_TEST;
312
313 Mesh mesh(*TestCommWorld, 2);
315
316 CPPUNIT_ASSERT(mesh.parallel_n_nodes() == 8);
317
318 mesh.stitch_surfaces(interface_left_id, interface_right_id, 1e-8 /*tolerance*/,
319 true/*clear_stitched_boundary_ids*/, false /*verbose*/,false/*use_binary_search*/,
320 false /*enforce_all_nodes_match_on_boundaries*/, false /*merge_boundary_nodes_all_or_nothing*/,
321 false /*prepare_after_stitching*/);
322
323 // ExodusII_IO(mesh).write("stitched_disjoint_neighbor_boundary_pairs.e");
324 CPPUNIT_ASSERT(mesh.parallel_n_nodes() == 6);
325
326 Mesh mesh2(*TestCommWorld, 2);
328
329 CPPUNIT_ASSERT(mesh2.parallel_n_nodes() == 8);
330
331 mesh2.stitch_surfaces(interface_left_id, interface_right_id, 1e-8 /*tolerance*/,
332 true/*clear_stitched_boundary_ids*/, false /*verbose*/,false/*use_binary_search*/,
333 false /*enforce_all_nodes_match_on_boundaries*/, false /*merge_boundary_nodes_all_or_nothing*/,
334 true /*prepare_after_stitching*/);
335
336 CPPUNIT_ASSERT(mesh2.parallel_n_nodes() == 6);
337
338 // ExodusII_IO(mesh2).write("stitched_disjoint_neighbor_boundary_pairs_2.e");
339 }

References build_two_disjoint_elems(), interface_left_id, interface_right_id, mesh, libMesh::MeshBase::parallel_n_nodes(), libMesh::DistributedMesh::parallel_n_nodes(), libMesh::UnstructuredMesh::stitch_surfaces(), and TestCommWorld.

◆ testTempJump()

void DisjointNeighborTest::testTempJump ( )
inlineprivate

Definition at line 220 of file disjoint_neighbor_test.C.

221 {
222 LOG_UNIT_TEST;
223
224 Mesh mesh(*TestCommWorld, 2);
226
227 // Check elem 0 (the left element)
228 if (const Elem * elem_left = mesh.query_elem_ptr(0))
229 {
230 // This processor owns elem 0, so we can safely check its neighbor
231 const Elem * elem_right_neighbor = elem_left->neighbor_ptr(1);
232
233 // The neighbor relationship should have been set up by prepare_for_use()
234 libmesh_assert(elem_right_neighbor);
235
236 // Verify the neighbor is indeed elem 1
237 LIBMESH_ASSERT_NUMBERS_EQUAL(elem_right_neighbor->id(), 1, 1e-15);
238 }
239
240 // Check elem 1 (the right element)
241 if (const Elem * elem_right = mesh.query_elem_ptr(1))
242 {
243 // This processor owns elem 1, so we can safely check its neighbor
244 const Elem * elem_left_neighbor = elem_right->neighbor_ptr(3);
245
246 // The neighbor relationship should be valid
247 libmesh_assert(elem_left_neighbor);
248
249 // Verify the neighbor is indeed elem 0
250 LIBMESH_ASSERT_NUMBERS_EQUAL(elem_left_neighbor->id(), 0, 1e-15);
251 }
252
255 es.add_system<LinearImplicitSystem>("TempJump");
256
257 const unsigned int u_var = sys.add_variable("u", FEType(FIRST, LAGRANGE));
258
259 std::set<boundary_id_type> left_bdy { left_id };
260 std::set<boundary_id_type> right_bdy { right_id };
261 std::vector<unsigned int> vars (1, u_var);
262
263 Parameters params;
264 params.set<Real>("b") = b;
266 WrappedFunction<Number> right_val(sys, &right_solution_fn, &params);
267 DirichletBoundary bc_left (left_bdy, vars, left_val);
268 DirichletBoundary bc_right(right_bdy, vars, right_val);
269
270 sys.get_dof_map().add_dirichlet_boundary(bc_left);
271 sys.get_dof_map().add_dirichlet_boundary(bc_right);
272
274
275 // Ensure the DofMap creates sparsity entries between neighboring elements.
276 // Without this, PETSc may report "New nonzero at (a,b) caused a malloc."
278
279 es.init();
280 sys.solve();
281
282 // ExodusII_IO(mesh).write_equation_systems("temperature_jump.e", es);
283
284 for (Real x=0.; x<=1.; x+=0.05)
285 {
286 if (std::abs(x - 0.5) < 1e-12) continue; // skip interface
287 for (Real y=0.; y<=1.; y+=0.05)
288 {
289 Point p(x,y);
290 const Number exact = heat_exact(p,params,"","");
291 const Number approx = sys.point_value(0,p);
292 LIBMESH_ASSERT_NUMBERS_EQUAL(exact, approx, 1e-2);
293 }
294 }
295 }
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
void set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
Allow the implicit_neighbor_dofs flag to be set programmatically.
Definition dof_map.C:1891
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
dof_id_type id() const
Definition dof_object.h:819
const Elem * neighbor_ptr(unsigned int i) const
Definition elem.h:2615
This is the EquationSystems class.
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
virtual void solve() override
Assembles & solves the linear system A*x=b.
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition parameters.h:75
T & set(const std::string &)
Definition parameters.h:494
Number point_value(unsigned int var, const Point &p, const bool insist_on_success=true, const NumericVector< Number > *sol=nullptr) const
Definition system.C:2219
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition system.C:1959
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition system.C:1344
const DofMap & get_dof_map() const
Definition system.h:2417
Wrap a libMesh-style function pointer into a FunctionBase object.
void assemble_temperature_jump(EquationSystems &es, const std::string &)
Number right_solution_fn(const Point &, const Parameters &params, const std::string &, const std::string &)
Number heat_exact(const Point &p, const Parameters &, const std::string &, const std::string &)
Number left_solution_fn(const Point &, const Parameters &, const std::string &, const std::string &)
libmesh_assert(ctx)

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_temperature_jump(), libMesh::System::attach_assemble_function(), b, build_two_disjoint_elems(), libMesh::FIRST, libMesh::System::get_dof_map(), heat_exact(), libMesh::DofObject::id(), libMesh::EquationSystems::init(), libMesh::LAGRANGE, left_id, left_solution_fn(), libMesh::libmesh_assert(), mesh, libMesh::Elem::neighbor_ptr(), libMesh::System::point_value(), libMesh::MeshBase::query_elem_ptr(), libMesh::Real, right_id, right_solution_fn(), libMesh::Parameters::set(), libMesh::DofMap::set_implicit_neighbor_dofs(), libMesh::LinearImplicitSystem::solve(), and TestCommWorld.

◆ testTempJumpLocalRefineFail()

void DisjointNeighborTest::testTempJumpLocalRefineFail ( )
inlineprivate

Definition at line 394 of file disjoint_neighbor_test.C.

395 {
396 LOG_UNIT_TEST;
397
398 try
399 {
400 Mesh mesh(*TestCommWorld, 2);
402 }
403 catch (const std::exception &e)
404 {
405 CPPUNIT_ASSERT(std::string(e.what()).find("Mesh refinement is not yet implemented") != std::string::npos);
406 }
407 }
void build_split_mesh_with_interface(Mesh &mesh, unsigned nx, unsigned ny, bool test_local_refinement=false)
static const int Nx
static const int Ny

References build_split_mesh_with_interface(), mesh, Nx, Ny, and TestCommWorld.

◆ testTempJumpRefine()

void DisjointNeighborTest::testTempJumpRefine ( )
inlineprivate

Definition at line 341 of file disjoint_neighbor_test.C.

342 {
343 LOG_UNIT_TEST;
344
345 Mesh mesh(*TestCommWorld, 2);
347
350 es.add_system<LinearImplicitSystem>("TempJump");
351
352 const unsigned int u_var = sys.add_variable("u", FEType(FIRST, LAGRANGE));
353
354 std::set<boundary_id_type> left_bdy { left_id };
355 std::set<boundary_id_type> right_bdy { right_id };
356 std::vector<unsigned int> vars (1, u_var);
357
358 Parameters params;
359 params.set<Real>("b") = b;
361 WrappedFunction<Number> right_val(sys, &right_solution_fn, &params);
362 DirichletBoundary bc_left (left_bdy, vars, left_val);
363 DirichletBoundary bc_right(right_bdy, vars, right_val);
364
365 sys.get_dof_map().add_dirichlet_boundary(bc_left);
366 sys.get_dof_map().add_dirichlet_boundary(bc_right);
367
369
370 // Ensure the DofMap creates sparsity entries between neighboring elements.
371 // Without this, PETSc may report "New nonzero at (a,b) caused a malloc."
373
374 es.init();
375 sys.solve();
376
377 // ExodusII_IO(mesh).write_equation_systems("temperature_jump_refine.e", es);
378
379 for (Real x=0.; x<=1.; x+=0.05)
380 {
381 if (std::abs(x - 0.5) < 1e-12) continue; // skip interface
382 for (Real y=0.; y<=1.; y+=0.05)
383 {
384 Point p(x,y);
385 const Number exact = heat_exact(p,params,"","");
386 const Number approx = sys.point_value(0,p);
387 LIBMESH_ASSERT_NUMBERS_EQUAL(exact, approx, 1e-2);
388 }
389 }
390 }

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_temperature_jump(), libMesh::System::attach_assemble_function(), b, build_split_mesh_with_interface(), libMesh::FIRST, libMesh::System::get_dof_map(), heat_exact(), libMesh::EquationSystems::init(), libMesh::LAGRANGE, left_id, left_solution_fn(), mesh, Nx, Ny, libMesh::System::point_value(), libMesh::Real, right_id, right_solution_fn(), libMesh::Parameters::set(), libMesh::DofMap::set_implicit_neighbor_dofs(), libMesh::LinearImplicitSystem::solve(), and TestCommWorld.


The documentation for this class was generated from the following file: