Skip to content

Commit 9f54211

Browse files
committed
BoundaryInfo: fix internal boundary lookup without AMR
Restore previous behavior in `side_with_boundary_id()`. When AMR is disabled, internal boundaries incorrectly returned invalid_uint. Add fallback return for non-external sides. Add `testInternalBoundary()` to verify correct side detection on a QUAD4 mesh. Update `disconnected_neighbor_test` to use default API.
1 parent 0bcdbb4 commit 9f54211

File tree

4 files changed

+63
-44
lines changed

4 files changed

+63
-44
lines changed

include/mesh/boundary_info.h

Lines changed: 1 addition & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -616,22 +616,12 @@ class BoundaryInfo : public ParallelObject
616616
/**
617617
* \returns A side of element \p elem whose associated boundary id is
618618
* \p boundary_id if such a side exists, and \p invalid_uint otherwise.
619-
* \p include_internal_boundary
620-
* If true, sides on internal logical boundaries are also considered.
621-
* If false (default), only sides on the geometric boundary
622-
* (i.e. sides without a valid neighbor element) are inspected.
623-
*
624-
* Setting this to true is useful for intentionally disconnected
625-
* interfaces (e.g. cohesive zone interfaces) where
626-
* a neighboring element exists but the side should still be treated
627-
* as a boundary for enforcement of interface conditions.
628619
*
629620
* \note If multiple sides of \p elem have the same id, only the lowest numbered
630621
* such side is returned.
631622
*/
632623
unsigned int side_with_boundary_id(const Elem * const elem,
633-
const boundary_id_type boundary_id,
634-
bool include_internal_boundary = false) const;
624+
const boundary_id_type boundary_id) const;
635625

636626
/**
637627
* \returns All sides of element \p elem whose associated boundary id is

src/mesh/boundary_info.C

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -2298,8 +2298,7 @@ void BoundaryInfo::renumber_node_id (boundary_id_type old_id,
22982298

22992299

23002300
unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
2301-
const boundary_id_type boundary_id_in,
2302-
bool include_internal_boundary) const
2301+
const boundary_id_type boundary_id_in) const
23032302
{
23042303
const Elem * searched_elem = elem;
23052304

@@ -2326,24 +2325,25 @@ unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
23262325
if (elem->neighbor_ptr(side) == nullptr)
23272326
return side;
23282327

2328+
// Internal boundary case
2329+
const Elem * p = elem;
2330+
2331+
#ifdef LIBMESH_ENABLE_AMR
23292332
// If we're on an internal boundary then we need to be sure
23302333
// it's the same internal boundary as our top_parent
2331-
bool aligned_with_parent = false;
2332-
#ifdef LIBMESH_ENABLE_AMR
2333-
const Elem * p = elem;
23342334
while (p != nullptr)
23352335
{
23362336
const Elem * parent = p->parent();
23372337
if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side))
23382338
break;
23392339
p = parent;
23402340
}
2341-
aligned_with_parent = (p == nullptr);
2341+
#else
2342+
// do not forget to return the internal boundary when AMR is disabled
2343+
return side;
23422344
#endif
2343-
// Two valid cases:
2344-
// 1. The internal face is geometrically aligned with the parent's boundary side (original behavior)
2345-
// 2. The caller explicitly allows logical/internal boundaries (e.g. CZM interfaces)
2346-
if (aligned_with_parent || include_internal_boundary)
2345+
// We're on that side of our top_parent; return it
2346+
if (!p)
23472347
return side;
23482348
}
23492349
// Otherwise we need to check if the child's ancestors have something on

tests/mesh/boundary_info.C

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ public:
2929
#if LIBMESH_DIM > 1
3030
CPPUNIT_TEST( testMesh );
3131
CPPUNIT_TEST( testRenumber );
32+
CPPUNIT_TEST( testInternalBoundary );
3233
# if LIBMESH_DIM > 2
3334
CPPUNIT_TEST( testSelectiveRenumber );
3435
# endif
@@ -1110,6 +1111,49 @@ public:
11101111
}
11111112
}
11121113

1114+
1115+
void testInternalBoundary()
1116+
{
1117+
LOG_UNIT_TEST;
1118+
1119+
Mesh mesh(*TestCommWorld);
1120+
1121+
// Build a 2x2 QUAD4 structured mesh on [0,1]x[0,1].
1122+
MeshTools::Generation::build_square(mesh,
1123+
2, 2,
1124+
0., 1.,
1125+
0., 1.,
1126+
QUAD4);
1127+
1128+
BoundaryInfo & bi = mesh.get_boundary_info();
1129+
1130+
// Select an element on the left side of the interior boundary
1131+
// (average x-coordinate < 0.5).
1132+
Elem * left_elem = nullptr;
1133+
for (auto & elem : mesh.active_element_ptr_range())
1134+
{
1135+
if (elem->vertex_average()(0) < 0.5)
1136+
{
1137+
left_elem = elem;
1138+
break;
1139+
}
1140+
}
1141+
1142+
CPPUNIT_ASSERT(left_elem);
1143+
unsigned int internal_side = 1;
1144+
CPPUNIT_ASSERT(left_elem->neighbor_ptr(internal_side) != nullptr);
1145+
1146+
// Assign a custom boundary ID on this internal side.
1147+
boundary_id_type BID = 7;
1148+
bi.add_side(left_elem, internal_side, BID);
1149+
1150+
mesh.prepare_for_use();
1151+
1152+
unsigned int found_side = bi.side_with_boundary_id(left_elem, BID);
1153+
CPPUNIT_ASSERT_EQUAL(internal_side, found_side);
1154+
}
1155+
1156+
11131157
};
11141158

11151159
CPPUNIT_TEST_SUITE_REGISTRATION( BoundaryInfoTest );

tests/systems/disconnected_neighbor_test.C

Lines changed: 8 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,12 @@ void assemble_temperature_jump(EquationSystems &es,
7878
fe_face_L->attach_quadrature_rule(&qface);
7979
fe_face_R->attach_quadrature_rule(&qface);
8080

81+
const auto &JxW_vol = fe->get_JxW();
82+
const auto &dphi_vol = fe->get_dphi();
83+
const auto &phi_L = fe_face_L->get_phi();
84+
const auto &phi_R = fe_face_R->get_phi();
85+
const auto &JxW_face = fe_face_L->get_JxW();
86+
8187
DenseMatrix<Number> Ke;
8288
DenseVector<Number> Fe;
8389
std::vector<dof_id_type> dof_indices;
@@ -96,34 +102,21 @@ void assemble_temperature_jump(EquationSystems &es,
96102
Fe.zero();
97103

98104
// --- Volume contribution ---
99-
fe->get_dphi(); // sets calculate_dphi = true
100-
fe->get_JxW(); // sets calculate_map = true
101105
fe->reinit(elem);
102-
const auto &JxW_vol = fe->get_JxW();
103-
const auto &dphi_vol = fe->get_dphi();
104-
105106
for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
106107
for (unsigned int i = 0; i < n_dofs; i++)
107108
for (unsigned int j = 0; j < n_dofs; j++)
108109
Ke(i, j) += conductance * JxW_vol[qp] * (dphi_vol[i][qp] * dphi_vol[j][qp]);
109110

110111
// --- Left-side interface ---
111-
unsigned int side = boundary.side_with_boundary_id(elem, interface_left_id, true /*include_internal_boundary*/);
112+
unsigned int side = boundary.side_with_boundary_id(elem, interface_left_id);
112113
if (side != libMesh::invalid_uint)
113114
{
114-
fe_face_L->get_phi(); // sets calculate_phi = true
115-
fe_face_L->get_JxW(); // sets calculate_map = true
116115
fe_face_L->reinit(elem, side);
117-
const auto &phi_L = fe_face_L->get_phi();
118-
const auto &JxW_face = fe_face_L->get_JxW();
119-
120116
const Elem *neighbor = elem->neighbor_ptr(side);
121117
libmesh_assert(neighbor);
122118
unsigned int n_side = neighbor->which_neighbor_am_i(elem);
123-
fe_face_R->get_phi();
124-
fe_face_R->get_JxW();
125119
fe_face_R->reinit(neighbor, n_side);
126-
const auto &phi_R = fe_face_R->get_phi();
127120

128121
std::vector<dof_id_type> dofs_L, dofs_R;
129122
dof_map.dof_indices(elem, dofs_L);
@@ -148,22 +141,14 @@ void assemble_temperature_jump(EquationSystems &es,
148141
}
149142

150143
// --- Right-side interface ---
151-
side = boundary.side_with_boundary_id(elem, interface_right_id, true /*include_internal_boundary*/);
144+
side = boundary.side_with_boundary_id(elem, interface_right_id);
152145
if (side != libMesh::invalid_uint)
153146
{
154-
fe_face_R->get_phi(); // sets calculate_phi = true
155-
fe_face_R->get_JxW(); // sets calculate_map = true
156147
fe_face_R->reinit(elem, side);
157-
const auto &phi_R = fe_face_R->get_phi();
158-
const auto &JxW_face = fe_face_R->get_JxW();
159-
160148
const Elem *neighbor = elem->neighbor_ptr(side);
161149
libmesh_assert(neighbor);
162150
unsigned int n_side = neighbor->which_neighbor_am_i(elem);
163-
fe_face_L->get_phi();
164-
fe_face_L->get_JxW();
165151
fe_face_L->reinit(neighbor, n_side);
166-
const auto &phi_L = fe_face_L->get_phi();
167152

168153
std::vector<dof_id_type> dofs_R, dofs_L;
169154
dof_map.dof_indices(elem, dofs_R);

0 commit comments

Comments
 (0)