Skip to content

Commit a06497a

Browse files
committed
Add ElemExtrema
This is to be used in tracing routines in the future
1 parent 98d5001 commit a06497a

13 files changed

Lines changed: 667 additions & 110 deletions

File tree

Makefile.in

Lines changed: 111 additions & 40 deletions
Large diffs are not rendered by default.

include/Makefile.in

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -730,6 +730,7 @@ include_HEADERS = \
730730
geom/edge_inf_edge2.h \
731731
geom/elem.h \
732732
geom/elem_cutter.h \
733+
geom/elem_extrema.h \
733734
geom/elem_hash.h \
734735
geom/elem_internal.h \
735736
geom/elem_quality.h \

include/geom/elem.h

Lines changed: 7 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@
3939
#include <algorithm>
4040
#include <cstddef>
4141
#include <iostream>
42-
#include <limits.h> // CHAR_BIT
42+
#include <limits.h> // CHAR_BIT, std::numeric_limits
4343
#include <set>
4444
#include <vector>
4545
#include <memory>
@@ -224,18 +224,13 @@ class Elem : public ReferenceCountedObject<Elem>,
224224
/**
225225
* A static integral constant representing an invalid subdomain id.
226226
* See also DofObject::{invalid_id, invalid_unique_id, invalid_processor_id}.
227-
*
228-
* \note We don't use the static_cast(-1) trick here since
229-
* \p subdomain_id_type is sometimes a *signed* integer for
230-
* compatibility reasons (see libmesh/id_types.h).
231-
*
232-
* \note Normally you can declare static const integral types
233-
* directly in the header file (C++ standard, 9.4.2/4) but
234-
* std::numeric_limits<T>::max() is not considered a "constant
235-
* expression". This one is therefore defined in elem.C.
236-
* http://stackoverflow.com/questions/2738435/using-numeric-limitsmax-in-constant-expressions
237227
*/
238-
static const subdomain_id_type invalid_subdomain_id;
228+
static constexpr subdomain_id_type invalid_subdomain_id = std::numeric_limits<subdomain_id_type>::max();
229+
230+
/**
231+
* A static integral constant representing an invalid index to a vertex.
232+
*/
233+
static constexpr unsigned short invalid_vertex = std::numeric_limits<unsigned short>::max();
239234

240235
/**
241236
* \returns A pointer to the "reference element" associated

include/geom/elem_extrema.h

Lines changed: 176 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,176 @@
1+
// The libMesh Finite Element Library.
2+
// Copyright (C) 2002-2022 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3+
4+
// This library is free software; you can redistribute it and/or
5+
// modify it under the terms of the GNU Lesser General Public
6+
// License as published by the Free Software Foundation; either
7+
// version 2.1 of the License, or (at your option) any later version.
8+
9+
// This library is distributed in the hope that it will be useful,
10+
// but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12+
// Lesser General Public License for more details.
13+
14+
// You should have received a copy of the GNU Lesser General Public
15+
// License along with this library; if not, write to the Free Software
16+
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17+
18+
19+
20+
#ifndef LIBMESH_ELEMEXTREMA_H
21+
#define LIBMESH_ELEMEXTREMA_H
22+
23+
#include "libmesh/libmesh_common.h"
24+
25+
#if LIBMESH_DIM > 1
26+
27+
#include "libmesh/elem.h"
28+
29+
#include <utility> // std::pair
30+
31+
namespace libMesh
32+
{
33+
34+
/**
35+
* The \p ElemExtrema is a helper class for defining whether or not
36+
* _something_ is at an element vertex or edge (an _extrema_).
37+
*/
38+
class ElemExtrema : public std::pair<unsigned short, unsigned short>
39+
{
40+
public:
41+
/**
42+
* Default constructor: sets entires to invalid (not at vertex or edge)
43+
*/
44+
ElemExtrema()
45+
: std::pair<unsigned short, unsigned short>(Elem::invalid_vertex,
46+
Elem::invalid_vertex) { }
47+
48+
/**
49+
* Constructor with given vertex indices \p v1 and \p v2.
50+
*/
51+
ElemExtrema(const unsigned short v1, const unsigned short v2)
52+
: std::pair<unsigned short, unsigned short>(v1, v2) { }
53+
54+
/**
55+
* @returns true if at the extrema (edge or vertex)
56+
*/
57+
bool at_extrema() const { return first != Elem::invalid_vertex; }
58+
59+
/**
60+
* @returns true if the vertex index data is invalid
61+
*/
62+
bool is_invalid() const
63+
{ return first == Elem::invalid_vertex && second == Elem::invalid_vertex; }
64+
65+
/**
66+
* @returns true if at a vertex
67+
*/
68+
bool at_vertex() const
69+
{ return first != Elem::invalid_vertex && second == Elem::invalid_vertex; }
70+
/**
71+
* @returns true if at vertex \p v
72+
*/
73+
bool at_vertex(const unsigned short v) const
74+
{ return first == v && second == Elem::invalid_vertex; }
75+
76+
/**
77+
* Invalidates the current state
78+
*/
79+
void invalidate()
80+
{
81+
first = Elem::invalid_vertex;
82+
second = Elem::invalid_vertex;
83+
}
84+
85+
/**
86+
* @returns The vertex ID when at a vertex
87+
*/
88+
unsigned short vertex() const
89+
{ libmesh_assert(at_vertex()); return first; }
90+
91+
/**
92+
* Prints the current state (at edge, at vertex, not at either)
93+
*/
94+
std::string print() const;
95+
96+
/**
97+
* Sets the "at vertex" state
98+
*/
99+
void set_vertex(const unsigned short v)
100+
{
101+
libmesh_assert_not_equal_to(v, Elem::invalid_vertex);
102+
first = v;
103+
second = Elem::invalid_vertex;
104+
}
105+
106+
#if LIBMESH_DIM > 2
107+
/**
108+
* @returns true if at an edge
109+
*/
110+
bool at_edge() const
111+
{ return first != Elem::invalid_vertex && second != Elem::invalid_vertex; }
112+
/**
113+
* @returns true if at the edge defined by vertices \p v1 and \p v2
114+
*/
115+
bool at_edge(const unsigned short v1, const unsigned short v2) const
116+
{
117+
libmesh_assert_not_equal_to(v1, Elem::invalid_vertex);
118+
libmesh_assert_not_equal_to(v2, Elem::invalid_vertex);
119+
return (first == v1 && second == v2) || (first == v2 && second == v1);
120+
}
121+
/**
122+
* @returns true if at the edge with index \p e on element \p
123+
*/
124+
bool at_edge(const Elem & elem, const unsigned short e) const;
125+
126+
/**
127+
* @returns The vertices that contain the edge when at an edge
128+
*/
129+
const std::pair<unsigned short, unsigned short> & edge_vertices() const
130+
{ libmesh_assert(at_edge()); return *this; }
131+
132+
/**
133+
* Sets the "at edge" state
134+
*/
135+
void set_edge(const unsigned short v1, const unsigned short v2)
136+
{
137+
libmesh_assert_not_equal_to(v1, Elem::invalid_vertex);
138+
libmesh_assert_not_equal_to(v2, Elem::invalid_vertex);
139+
first = v1;
140+
second = v2;
141+
}
142+
/**
143+
* Sets the "at edge" state
144+
*/
145+
void set_edge(const std::pair<unsigned short, unsigned short> & vs)
146+
{ set_edge(vs.first, vs.second); }
147+
148+
/**
149+
* @returns The edge when at an edge
150+
*/
151+
std::unique_ptr<const Elem> build_edge(const Elem & elem) const;
152+
#endif
153+
154+
/**
155+
* @returns The vertex point when at a vertex
156+
*/
157+
const Point & vertex_point(const Elem & elem) const
158+
{ return elem.point(vertex()); }
159+
160+
/**
161+
* @returns Whether or not the current state (at vertex/edge) is valid
162+
* for the given \p elem and \p point. \p tol is used as the tolerance
163+
* to pass to the equality checks.
164+
*
165+
* This ONLY checks for validity when at_extrema().
166+
*/
167+
bool is_valid(const Elem & elem, const Point & point, const Real tol = TOLERANCE) const;
168+
};
169+
170+
} // namespace libMesh
171+
172+
std::ostream & operator<<(std::ostream & os, const libMesh::ElemExtrema & elem_extrema);
173+
174+
#endif // LIBMESH_DIM > 1
175+
176+
#endif // LIBMESH_ELEMEXTREMA_H

include/include_HEADERS

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -140,6 +140,7 @@ include_HEADERS = \
140140
geom/edge_inf_edge2.h \
141141
geom/elem.h \
142142
geom/elem_cutter.h \
143+
geom/elem_extrema.h \
143144
geom/elem_hash.h \
144145
geom/elem_internal.h \
145146
geom/elem_quality.h \

include/libmesh/Makefile.am

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -130,6 +130,7 @@ BUILT_SOURCES = \
130130
edge_inf_edge2.h \
131131
elem.h \
132132
elem_cutter.h \
133+
elem_extrema.h \
133134
elem_hash.h \
134135
elem_internal.h \
135136
elem_quality.h \
@@ -950,6 +951,9 @@ elem.h: $(top_srcdir)/include/geom/elem.h
950951
elem_cutter.h: $(top_srcdir)/include/geom/elem_cutter.h
951952
$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
952953

954+
elem_extrema.h: $(top_srcdir)/include/geom/elem_extrema.h
955+
$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
956+
953957
elem_hash.h: $(top_srcdir)/include/geom/elem_hash.h
954958
$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
955959

include/libmesh/Makefile.in

Lines changed: 16 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -547,19 +547,19 @@ BUILT_SOURCES = auto_ptr.h default_coupling.h dirichlet_boundaries.h \
547547
cell_pyramid5.h cell_tet.h cell_tet10.h cell_tet14.h \
548548
cell_tet4.h compare_elems_by_level.h edge.h edge_edge2.h \
549549
edge_edge3.h edge_edge4.h edge_inf_edge2.h elem.h \
550-
elem_cutter.h elem_hash.h elem_internal.h elem_quality.h \
551-
elem_range.h elem_side_builder.h face.h face_inf_quad.h \
552-
face_inf_quad4.h face_inf_quad6.h face_quad.h face_quad4.h \
553-
face_quad4_shell.h face_quad8.h face_quad8_shell.h \
554-
face_quad9.h face_tri.h face_tri3.h face_tri3_shell.h \
555-
face_tri3_subdivision.h face_tri6.h face_tri7.h node.h \
556-
node_elem.h node_range.h plane.h point.h reference_elem.h \
557-
remote_elem.h side.h sphere.h stored_range.h surface.h \
558-
abaqus_io.h boundary_info.h boundary_mesh.h checkpoint_io.h \
559-
distributed_mesh.h dyna_io.h ensight_io.h exodusII_io.h \
560-
exodusII_io_helper.h exodus_header_info.h fro_io.h gmsh_io.h \
561-
gmv_io.h gnuplot_io.h inf_elem_builder.h matlab_io.h \
562-
medit_io.h mesh.h mesh_base.h mesh_communication.h \
550+
elem_cutter.h elem_extrema.h elem_hash.h elem_internal.h \
551+
elem_quality.h elem_range.h elem_side_builder.h face.h \
552+
face_inf_quad.h face_inf_quad4.h face_inf_quad6.h face_quad.h \
553+
face_quad4.h face_quad4_shell.h face_quad8.h \
554+
face_quad8_shell.h face_quad9.h face_tri.h face_tri3.h \
555+
face_tri3_shell.h face_tri3_subdivision.h face_tri6.h \
556+
face_tri7.h node.h node_elem.h node_range.h plane.h point.h \
557+
reference_elem.h remote_elem.h side.h sphere.h stored_range.h \
558+
surface.h abaqus_io.h boundary_info.h boundary_mesh.h \
559+
checkpoint_io.h distributed_mesh.h dyna_io.h ensight_io.h \
560+
exodusII_io.h exodusII_io_helper.h exodus_header_info.h \
561+
fro_io.h gmsh_io.h gmv_io.h gnuplot_io.h inf_elem_builder.h \
562+
matlab_io.h medit_io.h mesh.h mesh_base.h mesh_communication.h \
563563
mesh_function.h mesh_generation.h mesh_input.h \
564564
mesh_inserter_iterator.h mesh_modification.h mesh_output.h \
565565
mesh_refinement.h mesh_serializer.h mesh_smoother.h \
@@ -1285,6 +1285,9 @@ elem.h: $(top_srcdir)/include/geom/elem.h
12851285
elem_cutter.h: $(top_srcdir)/include/geom/elem_cutter.h
12861286
$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
12871287

1288+
elem_extrema.h: $(top_srcdir)/include/geom/elem_extrema.h
1289+
$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
1290+
12881291
elem_hash.h: $(top_srcdir)/include/geom/elem_hash.h
12891292
$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
12901293

src/geom/elem.C

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -82,7 +82,6 @@
8282
#include <array>
8383
#include <iterator> // for std::ostream_iterator
8484
#include <sstream>
85-
#include <limits> // for std::numeric_limits<>
8685
#include <cmath> // for std::sqrt()
8786
#include <memory>
8887

@@ -93,8 +92,6 @@ namespace libMesh
9392
Threads::spin_mutex parent_indices_mutex;
9493
Threads::spin_mutex parent_bracketing_nodes_mutex;
9594

96-
const subdomain_id_type Elem::invalid_subdomain_id = std::numeric_limits<subdomain_id_type>::max();
97-
9895
// Initialize static member variables
9996
const unsigned int Elem::max_n_nodes;
10097

src/geom/elem_extrema.C

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
// The libMesh Finite Element Library.
2+
// Copyright (C) 2002-2022 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3+
4+
// This library is free software; you can redistribute it and/or
5+
// modify it under the terms of the GNU Lesser General Public
6+
// License as published by the Free Software Foundation; either
7+
// version 2.1 of the License, or (at your option) any later version.
8+
9+
// This library is distributed in the hope that it will be useful,
10+
// but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12+
// Lesser General Public License for more details.
13+
14+
// You should have received a copy of the GNU Lesser General Public
15+
// License along with this library; if not, write to the Free Software
16+
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17+
18+
19+
// Local includes
20+
#include "libmesh/elem_extrema.h"
21+
22+
#if LIBMESH_DIM > 1
23+
24+
namespace libMesh
25+
{
26+
27+
#if LIBMESH_DIM > 2
28+
bool ElemExtrema::at_edge(const Elem & elem, const unsigned short e) const
29+
{
30+
libmesh_assert_less(e, elem.n_edges());
31+
if (!at_edge())
32+
return false;
33+
const unsigned int * nodes_on_edge_ptr = elem.nodes_on_edge_ptr(e);
34+
return at_edge(nodes_on_edge_ptr[0], nodes_on_edge_ptr[1]);
35+
}
36+
37+
std::unique_ptr<const Elem> ElemExtrema::build_edge(const Elem & elem) const
38+
{
39+
libmesh_assert(at_edge());
40+
libmesh_assert_less(first, elem.n_vertices());
41+
libmesh_assert_less(second, elem.n_vertices());
42+
43+
for (const auto e : elem.edge_index_range())
44+
if (elem.is_node_on_edge(first, e) && elem.is_node_on_edge(second, e))
45+
return elem.build_edge_ptr(e);
46+
47+
libmesh_error_msg("Element does not contain vertices in ElemExtrema");
48+
}
49+
#endif
50+
51+
std::string ElemExtrema::print() const
52+
{
53+
std::stringstream oss;
54+
if (at_vertex())
55+
oss << "at vertex " << vertex();
56+
#if LIBMESH_DIM > 2
57+
else if (at_edge())
58+
oss << "at edge with vertices " << first << " and " << second;
59+
#endif
60+
else
61+
oss << "not at extrema";
62+
return oss.str();
63+
}
64+
65+
bool ElemExtrema::is_valid(const Elem & elem, const Point & point, const Real tol) const
66+
{
67+
libmesh_assert(first == Elem::invalid_vertex || first < elem.n_vertices());
68+
libmesh_assert(second == Elem::invalid_vertex || second < elem.n_vertices());
69+
70+
if (at_vertex())
71+
return vertex_point(elem).absolute_fuzzy_equals(point, tol);
72+
#if LIBMESH_DIM > 2
73+
if (elem.dim() == 3 && at_edge())
74+
return build_edge(elem)->contains_point(point, tol);
75+
#endif
76+
77+
return true;
78+
}
79+
80+
} // namespace libMesh
81+
82+
std::ostream &
83+
operator<<(std::ostream & os, const libMesh::ElemExtrema & elem_extrema)
84+
{
85+
os << elem_extrema.print();
86+
return os;
87+
}
88+
89+
#endif // LIBMESH_DIM > 1

0 commit comments

Comments
 (0)