Skip to content

Commit 08d51bb

Browse files
loganharbourloganharbour
committed
Add IntersectionTools::within_edge
1 parent 53ec530 commit 08d51bb

File tree

3 files changed

+214
-93
lines changed

3 files changed

+214
-93
lines changed

‎include/geom/intersection_tools.h‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -69,6 +69,34 @@ bool collinear(const Point & p1,
6969
const Point & p3,
7070
const Real tol = TOLERANCE);
7171

72+
/**
73+
* @returns Whether or not the edges on element \p element are
74+
* individually collinear.
75+
*/
76+
bool edges_are_collinear(const Elem & elem, const Real tol = TOLERANCE);
77+
78+
/**
79+
* \returns True if the given point is contained within an edge an element
80+
* @param elem The element
81+
* @param p The point
82+
* @param corner To be filled with the edge/vertex that the point
83+
* is within/at or, if any (must be initially invalid)
84+
* @param linearize Whether or not to "linearize" the check, if this
85+
* is set to false and edges are found to not be collinear, an error
86+
* is thrown
87+
*
88+
* \p corner will be set to an "at vertex" state if the point is
89+
* both within the edge _and_ at a vertex.
90+
*
91+
* This method is only implemented for three-dimensional, finite
92+
* elements.
93+
*/
94+
bool within_edge(const Elem & elem,
95+
const Point & p,
96+
ElemCorner & corner,
97+
const bool linearize = false,
98+
const Real tol = TOLERANCE);
99+
72100
/**
73101
* \returns True if the given point is contained within an edge
74102
* on the given side of an element
@@ -94,6 +122,20 @@ bool within_edge_on_side(const Elem & elem,
94122
const bool linearize = false,
95123
const Real tol = TOLERANCE);
96124

125+
namespace detail
126+
{
127+
/**
128+
* Internal method for checking whether or not the point \p
129+
* is within the edge defined by vertices \p v1 and \p v2
130+
* on element \p elem.
131+
*/
132+
bool _within_edge(const Elem & elem,
133+
const Point & p,
134+
ElemCorner & corner,
135+
const unsigned int v1,
136+
const unsigned int v2,
137+
const Real tol);
138+
} // namespace detail
97139
} // namespace IntersectionTools
98140
} // namespace libMesh
99141

‎src/geom/intersection_tools.C‎

Lines changed: 87 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,27 @@ bool collinear(const Point & p1,
7373
return (std::abs(p2_p1 * p3_p2) / denom) > ((Real)1 - tol);
7474
}
7575

76+
bool edges_are_collinear(const Elem & elem, const Real tol)
77+
{
78+
if (elem.default_order() < 2)
79+
return true;
80+
81+
for (const auto e : elem.edge_index_range())
82+
{
83+
// we should expect at least 3 nodes/edges for our higher order elems
84+
libmesh_assert_greater_equal(elem.n_nodes_on_edge(e), 3);
85+
86+
const unsigned int * edge_nodes_map = elem.nodes_on_edge_ptr(e);
87+
if (!collinear(elem.point(edge_nodes_map[0]),
88+
elem.point(edge_nodes_map[1]),
89+
elem.point(edge_nodes_map[2]),
90+
tol))
91+
return false;
92+
}
93+
94+
return true;
95+
}
96+
7697
bool within_edge_on_side(const Elem & elem,
7798
const Point & p,
7899
const unsigned short s,
@@ -84,21 +105,10 @@ bool within_edge_on_side(const Elem & elem,
84105
libmesh_assert(corner.is_invalid());
85106
libmesh_assert_equal_to(elem.dim(), 3);
86107

87-
// For higher order than linear without linearization, make sure
88-
// that our edges are collinear
89-
if (elem.default_order() > 1 && !linearize)
90-
for (const auto e : elem.edge_index_range())
91-
{
92-
// we should expect at least 3 nodes/edges for our higher order elems
93-
libmesh_assert_greater_equal(elem.n_nodes_on_edge(e), 3);
94-
95-
const unsigned int * edge_nodes_map = elem.nodes_on_edge_ptr(e);
96-
if (!collinear(elem.point(edge_nodes_map[0]),
97-
elem.point(edge_nodes_map[1]),
98-
elem.point(edge_nodes_map[2])))
99-
libmesh_error_msg("Failed to use Cell::without_edge_on_side without linearization "
100-
"because an edge was found that is not collinear.");
101-
}
108+
// Without linearization, make sure that our edges are collinear
109+
if (!linearize && !edges_are_collinear(elem, tol))
110+
libmesh_error_msg("Failed to use Cell::without_edge_on_side without linearization "
111+
"because an edge was found that is not collinear.");
102112

103113
const auto vs = elem.n_vertices_on_side(s);
104114
const unsigned int * side_nodes_map = elem.nodes_on_side_ptr(s);
@@ -109,35 +119,73 @@ bool within_edge_on_side(const Elem & elem,
109119
auto last_v = side_nodes_map[vs - 1];
110120
for (const auto side_v : make_range(vs))
111121
{
112-
const auto other_v = side_nodes_map[side_v];
113-
const auto within_result = within_segment(elem.point(last_v), elem.point(other_v), p, tol);
114-
if (within_result == NOT_WITHIN)
115-
{
116-
last_v = side_nodes_map[side_v];
117-
continue;
118-
}
119-
120-
if (within_result == BETWEEN)
121-
{
122-
corner.set_edge(last_v, other_v);
123-
libmesh_assert(corner.build_edge(elem)->contains_point(p, tol));
124-
}
125-
else if (within_result == AT_BEGINNING)
126-
{
127-
corner.set_vertex(last_v);
128-
libmesh_assert(elem.point(last_v).relative_fuzzy_equals(p, tol));
129-
}
130-
else
131-
{
132-
corner.set_vertex(other_v);
133-
libmesh_assert(elem.point(other_v).relative_fuzzy_equals(p, tol));
134-
}
135-
return true;
122+
if (detail::_within_edge(elem, p, corner, last_v, side_nodes_map[side_v], tol))
123+
return true;
124+
last_v = side_nodes_map[side_v];
125+
}
126+
127+
return false;
128+
}
129+
130+
bool within_edge(const Elem & elem,
131+
const Point & p,
132+
ElemCorner & corner,
133+
const bool linearize,
134+
const Real tol)
135+
{
136+
libmesh_assert(corner.is_invalid());
137+
138+
// Without linearization, make sure that our edges are collinear
139+
if (!linearize && !edges_are_collinear(elem, tol))
140+
libmesh_error_msg("Failed to use Cell::without_edge without linearization "
141+
"because an edge was found that is not collinear.");
142+
143+
for (const auto e : elem.edge_index_range())
144+
{
145+
const unsigned int * edge_nodes_map = elem.nodes_on_edge_ptr(e);
146+
if (detail::_within_edge(elem, p, corner, edge_nodes_map[0], edge_nodes_map[1], tol))
147+
return true;
136148
}
137149

138150
return false;
139151
}
140152

153+
namespace detail
154+
{
155+
bool _within_edge(const Elem & elem,
156+
const Point & p,
157+
ElemCorner & corner,
158+
const unsigned int v1,
159+
const unsigned int v2,
160+
const Real tol)
161+
{
162+
libmesh_assert(corner.is_invalid());
163+
libmesh_assert_less(v1, elem.n_vertices());
164+
libmesh_assert_less(v2, elem.n_vertices());
165+
166+
const auto within_result = within_segment(elem.point(v1), elem.point(v2), p, tol);
167+
if (within_result == NOT_WITHIN)
168+
return false;
169+
170+
if (within_result == BETWEEN)
171+
{
172+
corner.set_edge(v1, v2);
173+
libmesh_assert(corner.build_edge(elem)->contains_point(p, tol));
174+
}
175+
else if (within_result == AT_BEGINNING)
176+
{
177+
corner.set_vertex(v1);
178+
libmesh_assert(elem.point(v1).relative_fuzzy_equals(p, tol));
179+
}
180+
else
181+
{
182+
corner.set_vertex(v2);
183+
libmesh_assert(elem.point(v2).relative_fuzzy_equals(p, tol));
184+
}
185+
186+
return true;
187+
}
188+
} // namespace detail
141189
} // namespace libMesh::IntersectionTools
142190

143191

‎tests/geom/intersection_tools_test.C‎

Lines changed: 85 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -117,80 +117,102 @@ public:
117117
elem_type);
118118
}
119119

120-
void test_within_edge_on_side()
120+
void test_within_edge()
121121
{
122122
LOG_UNIT_TEST;
123123

124-
if (_mesh->mesh_dimension() != 3)
125-
return;
124+
const auto is_3d = _mesh->mesh_dimension() == 3;
126125

127126
// check locations at every node
128127
for (const auto & elem : _mesh->active_local_element_ptr_range())
129-
for (const auto s : elem->side_index_range())
130-
for (const auto e : elem->edge_index_range())
131-
for (const auto n : elem->nodes_on_edge(e))
132-
{
133-
ElemCorner corner;
134-
const auto within = IntersectionTools::within_edge_on_side(*elem,
135-
elem->point(n),
136-
s,
137-
corner);
138-
139-
CPPUNIT_ASSERT_EQUAL(within, elem->is_node_on_side(n, s));
140-
if (elem->is_node_on_side(n, s))
128+
for (const auto e : elem->edge_index_range())
129+
for (const auto n : elem->nodes_on_edge(e))
130+
{
131+
ElemCorner corner;
132+
const auto within = IntersectionTools::within_edge(*elem, elem->point(n), corner);
133+
CPPUNIT_ASSERT(within);
134+
135+
if (is_3d)
136+
for (const auto s : elem->side_index_range())
141137
{
142-
CPPUNIT_ASSERT_EQUAL(elem->is_vertex(n), corner.at_vertex(n));
143-
CPPUNIT_ASSERT_EQUAL(elem->is_vertex(n), !corner.at_edge(*elem, e));
138+
ElemCorner corner_on_side;
139+
const auto within_on_side = IntersectionTools::within_edge_on_side(*elem,
140+
elem->point(n),
141+
s,
142+
corner_on_side);
143+
144+
const auto node_is_on_side = elem->is_node_on_side(n, s);
145+
CPPUNIT_ASSERT_EQUAL(within_on_side, node_is_on_side);
146+
if (node_is_on_side)
147+
{
148+
CPPUNIT_ASSERT_EQUAL(corner.at_vertex(n), corner_on_side.at_vertex(n));
149+
CPPUNIT_ASSERT_EQUAL(corner.at_edge(*elem, e), corner_on_side.at_edge(*elem, e));
150+
CPPUNIT_ASSERT_EQUAL(elem->is_vertex(n), corner_on_side.at_vertex(n));
151+
CPPUNIT_ASSERT_EQUAL(elem->is_vertex(n), !corner_on_side.at_edge(*elem, e));
152+
}
144153
}
145-
}
154+
}
146155

147156
// cut edges into segments
148157
for (const auto & elem : _mesh->active_local_element_ptr_range())
149158
for (const auto e : elem->edge_index_range())
150-
for (const auto s : elem->side_index_range())
151-
if (elem->is_edge_on_side(e, s))
152-
{
153-
const auto nodes_on_edge = elem->nodes_on_edge(e);
154-
const auto & p1 = elem->point(nodes_on_edge[0]);
155-
const auto & p2 = elem->point(nodes_on_edge[1]);
156-
const auto length_vec = p2 - p1;
157-
const auto length = length_vec.norm();
158-
const auto p1_to_p2 = length_vec / length;
159-
160-
int segments = 5;
161-
Real dx = (Real)1 / segments * length;
162-
for (const auto i : make_range(-1, segments + 1))
163-
{
164-
const auto p = p1 + Real(i) * dx * p1_to_p2;
165-
ElemCorner corner;
166-
const auto within = IntersectionTools::within_edge_on_side(*elem,
167-
p,
168-
s,
169-
corner);
170-
171-
CPPUNIT_ASSERT_EQUAL(within, i >= 0 && i <= segments);
172-
CPPUNIT_ASSERT_EQUAL(corner.at_vertex(nodes_on_edge[0]), i == 0);
173-
CPPUNIT_ASSERT_EQUAL(corner.at_vertex(nodes_on_edge[1]), i == segments);
174-
CPPUNIT_ASSERT_EQUAL(corner.at_edge(*elem, e), i > 0 && i < segments);
175-
}
176-
}
159+
{
160+
const auto nodes_on_edge = elem->nodes_on_edge(e);
161+
const auto & p1 = elem->point(nodes_on_edge[0]);
162+
const auto & p2 = elem->point(nodes_on_edge[1]);
163+
const auto length_vec = p2 - p1;
164+
const auto length = length_vec.norm();
165+
const auto p1_to_p2 = length_vec / length;
166+
167+
int segments = 5;
168+
Real dx = (Real)1 / segments * length;
169+
for (const auto i : make_range(-1, segments + 1))
170+
{
171+
const auto p = p1 + Real(i) * dx * p1_to_p2;
172+
173+
ElemCorner corner;
174+
const auto within = IntersectionTools::within_edge(*elem, p, corner);
175+
CPPUNIT_ASSERT_EQUAL(within, i >= 0 && i <= segments);
176+
177+
if (is_3d)
178+
for (const auto s : elem->side_index_range())
179+
if (elem->is_edge_on_side(e, s))
180+
{
181+
ElemCorner corner_on_side;
182+
const auto within_on_side = IntersectionTools::within_edge_on_side(*elem,
183+
p,
184+
s,
185+
corner_on_side);
186+
187+
CPPUNIT_ASSERT_EQUAL(within_on_side, within);
188+
CPPUNIT_ASSERT_EQUAL(corner_on_side.at_vertex(nodes_on_edge[0]), i == 0);
189+
CPPUNIT_ASSERT_EQUAL(corner_on_side.at_vertex(nodes_on_edge[1]), i == segments);
190+
CPPUNIT_ASSERT_EQUAL(corner_on_side.at_edge(*elem, e), i > 0 && i < segments);
191+
}
192+
}
193+
}
177194

178195
// check elem centroids
179196
for (const auto & elem : _mesh->active_local_element_ptr_range())
180-
for (const auto s : elem->side_index_range())
181-
{
182-
ElemCorner corner;
183-
CPPUNIT_ASSERT(!IntersectionTools::within_edge_on_side(*elem,
184-
elem->vertex_average(),
185-
s,
186-
corner));
187-
}
197+
{
198+
const auto p = elem->vertex_average();
199+
200+
ElemCorner corner;
201+
CPPUNIT_ASSERT(!IntersectionTools::within_edge(*elem, p, corner));
202+
203+
if (is_3d)
204+
for (const auto s : elem->side_index_range())
205+
{
206+
ElemCorner corner_on_side;
207+
CPPUNIT_ASSERT(!IntersectionTools::within_edge_on_side(*elem, p, s, corner_on_side));
208+
}
209+
}
188210
}
189211

190212
};
191213

192214
#define MESHEDINTERSECTIONTOOLSTEST \
193-
CPPUNIT_TEST( test_within_edge_on_side );
215+
CPPUNIT_TEST( test_within_edge );
194216

195217
#define INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(elemtype) \
196218
class MeshedIntersectionToolsTest_##elemtype : public MeshedIntersectionToolsTest<elemtype> { \
@@ -209,6 +231,15 @@ public:
209231
\
210232
CPPUNIT_TEST_SUITE_REGISTRATION( MeshedIntersectionToolsTest_##elemtype )
211233

234+
#if LIBMESH_DIM > 1
235+
INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(TRI3);
236+
INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(TRI6);
237+
INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(TRI7);
238+
239+
INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(QUAD4);
240+
INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(QUAD8);
241+
INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(QUAD9);
242+
#endif // LIBMESH_DIM > 1
212243

213244
#if LIBMESH_DIM > 2
214245
INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(TET4);

0 commit comments

Comments
 (0)