Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Mesh/UnsMesh.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2025 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief 3D unstructured mesh class declaration
10 : : \details 3D unstructured mesh class declaration. This mesh class currently
11 : : supports line, triangle, and tetrahedron elements.
12 : : */
13 : : // *****************************************************************************
14 : : #ifndef UnsMesh_h
15 : : #define UnsMesh_h
16 : :
17 : : #include <vector>
18 : : #include <array>
19 : : #include <memory>
20 : : #include <tuple>
21 : : #include <map>
22 : : #include <unordered_set>
23 : : #include <unordered_map>
24 : :
25 : : #include "NoWarning/sip_hash.hpp"
26 : :
27 : : #include "Types.hpp"
28 : : #include "ContainerUtil.hpp"
29 : :
30 : : namespace tk {
31 : :
32 : : //! Highway hash "secret" key
33 : : //! \note No reason for these particular numbers, taken from highwayhash tests.
34 : : static constexpr highwayhash::HH_U64 hh_key[2] =
35 : : { 0x0706050403020100ULL, 0x0F0E0D0C0B0A0908ULL };
36 : :
37 : : //! 3D unstructured mesh class
38 : : class UnsMesh {
39 : :
40 : : private:
41 : : //! Union to access an C-style array of std::size_t as an array of bytes
42 : : //! \tparam N Number of entries to hold
43 : : //! \see UnsMesh::Hash
44 : : template< std::size_t N >
45 : : union Shaper {
46 : : char bytes[ N*sizeof(std::size_t) ];
47 : : std::size_t sizets[ N ];
48 : : };
49 : :
50 : : public:
51 : : using Coords = std::array< std::vector< real >, 3 >;
52 : : using Coord = std::array< real, 3 >;
53 : : using CoordMap = std::unordered_map< std::size_t, Coord >;
54 : :
55 : : //! Alias for storing a mesh chunk
56 : : //! \details The first vector is the element connectivity (local mesh node
57 : : //! IDs), the second vector is the global node IDs of owned elements,
58 : : //! while the third one is a map of global(key)->local(value) node IDs.
59 : : using Chunk = std::tuple< std::vector< std::size_t >,
60 : : std::vector< std::size_t >,
61 : : std::unordered_map< std::size_t, std::size_t > >;
62 : :
63 : : /** @name Aliases for element primitives */
64 : : ///@{
65 : : //! Edge: node IDs of two end-points
66 : : using Edge = std::array< std::size_t, 2 >;
67 : : //! Face: node IDs of a triangle (tetrahedron face)
68 : : using Face = std::array< std::size_t, 3 >;
69 : : //! Tet: node IDs of a tetrahedron
70 : : using Tet = std::array< std::size_t, 4 >;
71 : : ///@}
72 : :
73 : : //! Hash function class for element primitives, given by node IDs
74 : : //! \tparam N Number of nodes describing element primitive. E.g., Edge:2,
75 : : //! Face:3, Tet:4.
76 : : template< std::size_t N >
77 : : struct Hash {
78 : : //! Function call operator computing hash of node IDs
79 : : //! \param[in] p Array of node IDs of element primitive
80 : : //! \return Unique hash value for the same array of node IDs
81 : : //! \note The order of the nodes does not matter: the IDs are sorted
82 : : //! before the hash is computed.
83 : 35769898 : std::size_t operator()( const std::array< std::size_t, N >& p ) const {
84 : : using highwayhash::SipHash;
85 : : Shaper< N > shaper;
86 [ + + ]: 123153374 : for (std::size_t i=0; i<N; ++i) shaper.sizets[i] = p[i];
87 [ + - ]: 71539796 : std::sort( std::begin(shaper.sizets), std::end(shaper.sizets) );
88 [ + - ]: 71539796 : return SipHash( hh_key, shaper.bytes, N*sizeof(std::size_t) );
89 : : }
90 : : };
91 : :
92 : : //! Comparitor function class for element primitives, given by node IDs
93 : : //! \tparam N Number of nodes describing element primitive. E.g., Edge:2,
94 : : //! Face:3, Tet:4.
95 : : template< std::size_t N >
96 : : struct Eq {
97 : : //! Function call operator computing equality of element primitives
98 : : //! \param[in] l Left element primitive given by array of node IDs
99 : : //! \param[in] r Right element primitive given by array of node IDs
100 : : //! \return True if l = r, false otherwise
101 : : //! \note The order of the nodes does not matter: the IDs are sorted
102 : : //! before equality is determined.
103 : 18366859 : bool operator()( const std::array< std::size_t, N >& l,
104 : : const std::array< std::size_t, N >& r ) const
105 : : {
106 : 18366859 : std::array< std::size_t, N > s = l, p = r;
107 [ + - ]: 36733718 : std::sort( begin(s), end(s) );
108 [ + - ]: 36733718 : std::sort( begin(p), end(p) );
109 [ + - ]: 36733718 : return s == p;
110 : : }
111 : : };
112 : :
113 : : //! Unique set of edges
114 : : using EdgeSet = std::unordered_set< Edge, Hash<2>, Eq<2> >;
115 : :
116 : : //! Unique set of faces
117 : : using FaceSet = std::unordered_set< Face, Hash<3>, Eq<3> >;
118 : :
119 : : //! Unique set of tets
120 : : using TetSet = std::unordered_set< Tet, Hash<4>, Eq<4> >;
121 : :
122 : : /** @name Constructors */
123 : : ///@{
124 : : //! Constructor without initializing anything
125 : 61 : explicit UnsMesh() : m_graphsize(0), m_triinpoel(),
126 : 61 : m_tetinpoel(), m_x(), m_y(), m_z() {}
127 : :
128 : : //! Constructor copying over element connectivity
129 : : explicit UnsMesh( const std::vector< std::size_t >& tetinp ) :
130 : : m_graphsize( graphsize( tetinp ) ),
131 : : m_tetinpoel( tetinp )
132 : : {
133 : : Assert( m_tetinpoel.size()%4 == 0,
134 : : "Size of tetinpoel must be divisible by 4" );
135 : : }
136 : :
137 : : //! Constructor swallowing element connectivity
138 : 4 : explicit UnsMesh( std::vector< std::size_t >&& tetinp ) :
139 : 4 : m_graphsize( graphsize( tetinp ) ),
140 : 4 : m_tetinpoel( std::move(tetinp) )
141 : : {
142 [ - + ][ - - ]: 4 : Assert( m_tetinpoel.size()%4 == 0,
[ - - ][ - - ]
143 : : "Size of tetinpoel must be divisible by 4" );
144 : 4 : }
145 : :
146 : : //! Constructor copying over element connectivity and point coordinates
147 : : explicit UnsMesh( const std::vector< std::size_t >& tetinp,
148 : : const std::vector< real >& X,
149 : : const std::vector< real >& Y,
150 : : const std::vector< real >& Z ) :
151 : : m_graphsize( graphsize( tetinp ) ),
152 : : m_tetinpoel( tetinp ),
153 : : m_x( X ),
154 : : m_y( Y ),
155 : : m_z( Z )
156 : : {
157 : : Assert( m_tetinpoel.size()%4 == 0,
158 : : "Size of tetinpoel must be divisible by 4" );
159 : : }
160 : :
161 : : //! \brief Constructor copying over element connectivity and array of point
162 : : //! coordinates
163 : 838 : explicit UnsMesh( const std::vector< std::size_t >& tetinp,
164 : 838 : const Coords& coord ) :
165 : 838 : m_graphsize( graphsize( tetinp ) ),
166 [ + - ]: 838 : m_tetinpoel( tetinp ),
167 [ + - ]: 838 : m_x( coord[0] ),
168 [ + - ]: 838 : m_y( coord[1] ),
169 [ + - ]: 1676 : m_z( coord[2] )
170 : : {
171 [ - + ][ - - ]: 838 : Assert( m_tetinpoel.size()%4 == 0,
[ - - ][ - - ]
172 : : "Size of tetinpoel must be divisible by 4" );
173 : 838 : }
174 : :
175 : : //! \brief Constructor copying over triangle element connectivity and array
176 : : //! of point coordinates
177 : 50 : explicit UnsMesh(
178 : : const Coords& coord,
179 : : const std::vector< std::size_t >& triinp,
180 : : const std::vector< std::string >& nodevarnames = {},
181 : : const std::vector< tk::real >& vartimes = {},
182 : : const std::vector< std::vector< std::vector< tk::real > > >&
183 : 50 : nodevars = {} ) :
184 : 50 : m_graphsize( graphsize( triinp ) ),
185 : 50 : m_triinpoel( triinp ),
186 [ + - ]: 50 : m_x( coord[0] ),
187 [ + - ]: 50 : m_y( coord[1] ),
188 [ + - ]: 50 : m_z( coord[2] ),
189 [ + - ]: 50 : m_nodevarnames( nodevarnames ),
190 [ + - ]: 50 : m_vartimes( vartimes ),
191 [ + - ]: 150 : m_nodevars( nodevars )
192 : : {
193 [ - + ][ - - ]: 50 : Assert( m_triinpoel.size()%3 == 0,
[ - - ][ - - ]
194 : : "Size of triinpoel must be divisible by 3" );
195 : 50 : }
196 : :
197 : : //! Constructor swallowing element connectivity and point coordinates
198 : : explicit UnsMesh( std::vector< std::size_t >&& tetinp,
199 : : std::vector< real >&& X,
200 : : std::vector< real >&& Y,
201 : : std::vector< real >&& Z ) :
202 : : m_graphsize( graphsize( tetinp ) ),
203 : : m_tetinpoel( std::move(tetinp) ),
204 : : m_x( std::move(X) ),
205 : : m_y( std::move(Y) ),
206 : : m_z( std::move(Z) )
207 : : {
208 : : Assert( m_tetinpoel.size()%4 == 0,
209 : : "Size of tetinpoel must be divisible by 4" );
210 : : }
211 : :
212 : : //! \brief Constructor swallowing element connectivity and array of point
213 : : //! coordinates
214 : : explicit UnsMesh( std::vector< std::size_t >&& tetinp, Coords&& coord ) :
215 : : m_graphsize( graphsize( tetinp ) ),
216 : : m_tetinpoel( std::move(tetinp) ),
217 : : m_x( std::move(coord[0]) ),
218 : : m_y( std::move(coord[1]) ),
219 : : m_z( std::move(coord[2]) )
220 : : {
221 : : Assert( m_tetinpoel.size()%4 == 0,
222 : : "Size of tetinpoel must be divisible by 4" );
223 : : }
224 : :
225 : : //! Constructor with connectivities and side set faces
226 : : explicit UnsMesh(
227 : : const std::vector< std::size_t >& tetinp,
228 : : const Coords& coord,
229 : : const std::map< int, std::vector< std::size_t > >& bf,
230 : : const std::vector< std::size_t >& triinp,
231 : : const std::map< int, std::vector< std::size_t > >& fid ) :
232 : : m_graphsize( graphsize( tetinp ) ),
233 : : m_triinpoel( triinp ),
234 : : m_tetinpoel( tetinp ),
235 : : m_x( coord[0] ),
236 : : m_y( coord[1] ),
237 : : m_z( coord[2] ),
238 : : m_bface( bf ),
239 : : m_faceid( fid )
240 : : {
241 : : Assert( m_tetinpoel.size() % 4 == 0,
242 : : "Size of tetinpoel must be divisible by 4" );
243 : : Assert( m_triinpoel.size() % 3 == 0,
244 : : "Size of triinpoel must be divisible by 3" );
245 : : }
246 : :
247 : : //! Constructor with connectivities and side set nodes
248 : 84 : explicit UnsMesh(
249 : : const std::vector< std::size_t >& tetinp,
250 : : const Coords& coord,
251 : 84 : const std::map< int, std::vector< std::size_t > >& bn ) :
252 : 84 : m_graphsize( graphsize( tetinp ) ),
253 [ + - ]: 84 : m_tetinpoel( tetinp ),
254 [ + - ]: 84 : m_x( coord[0] ),
255 [ + - ]: 84 : m_y( coord[1] ),
256 [ + - ]: 84 : m_z( coord[2] ),
257 [ + - ]: 168 : m_bnode( bn )
258 : : {
259 [ - + ][ - - ]: 84 : Assert( m_tetinpoel.size() % 4 == 0,
[ - - ][ - - ]
260 : : "Size of tetinpoel must be divisible by 4" );
261 : 84 : }
262 : : ///@}
263 : :
264 : : /** @name Point coordinates accessors */
265 : : ///@{
266 : 37886 : const std::vector< real >& x() const noexcept { return m_x; }
267 : 37886 : const std::vector< real >& y() const noexcept { return m_y; }
268 : 37886 : const std::vector< real >& z() const noexcept { return m_z; }
269 : 80362 : std::vector< real >& x() noexcept { return m_x; }
270 : 80362 : std::vector< real >& y() noexcept { return m_y; }
271 : 80362 : std::vector< real >& z() noexcept { return m_z; }
272 : : ///@}
273 : :
274 : : /** @name Number of nodes accessors */
275 : : ///@{
276 : 37908 : std::size_t nnode() const noexcept { return m_x.size(); }
277 : : std::size_t nnode() noexcept { return m_x.size(); }
278 : : ///@}
279 : :
280 : : /** @name Graph size accessors */
281 : : ///@{
282 : : const std::size_t& size() const noexcept { return m_graphsize; }
283 : 36 : std::size_t& size() noexcept { return m_graphsize; }
284 : : ///@}
285 : :
286 : : //! Total number of elements accessor
287 : : std::size_t nelem() const noexcept {
288 : : return m_triinpoel.size()/3 + m_tetinpoel.size()/4;
289 : : }
290 : :
291 : : //! Number of element blocks accessor
292 : 983 : std::size_t neblk() const noexcept {
293 : 983 : return !m_triinpoel.empty() + !m_tetinpoel.empty();
294 : : }
295 : :
296 : : /** @name Triangle elements connectivity accessors */
297 : : ///@{
298 : 24721 : const std::vector< std::size_t >& triinpoel() const noexcept
299 : 24721 : { return m_triinpoel; }
300 : 202341 : std::vector< std::size_t >& triinpoel() noexcept { return m_triinpoel; }
301 : : ///@}
302 : :
303 : : /** @name Tetrahedra elements connectivity accessors */
304 : : ///@{
305 : 377382 : const std::vector< std::size_t >& tetinpoel() const noexcept
306 : 377382 : { return m_tetinpoel; }
307 : 1910318 : std::vector< std::size_t >& tetinpoel() noexcept { return m_tetinpoel; }
308 : : ///@}
309 : :
310 : : /** @name Side set face list accessors */
311 : : ///@{
312 : 1966 : const std::map< int, std::vector< std::size_t > >& bface() const noexcept
313 : 1966 : { return m_bface; }
314 : 1047 : std::map< int, std::vector< std::size_t > >& bface() noexcept
315 : 1047 : { return m_bface; }
316 : : ///@}
317 : :
318 : : /** @name Side set face id accessors */
319 : : ///@{
320 : 17 : const std::map< int, std::vector< std::size_t > >& faceid() const noexcept
321 : 17 : { return m_faceid; }
322 : 1045 : std::map< int, std::vector< std::size_t > >& faceid() noexcept
323 : 1045 : { return m_faceid; }
324 : : ///@}
325 : :
326 : : /** @name Side set node list accessors */
327 : : ///@{
328 : 1966 : const std::map< int, std::vector< std::size_t > >& bnode() const noexcept
329 : 1966 : { return m_bnode; }
330 : : std::map< int, std::vector< std::size_t > >& bnode() noexcept
331 : : { return m_bnode; }
332 : : ///@}
333 : :
334 : : /** @name Node variable names accessors */
335 : : ///@{
336 : 983 : const std::vector< std::string >& nodevarnames() const noexcept
337 : 983 : { return m_nodevarnames; }
338 : 81 : std::vector< std::string >& nodevarnames() noexcept
339 : 81 : { return m_nodevarnames; }
340 : : ///@}
341 : :
342 : : /** @name Variable times accessors */
343 : : ///@{
344 : 993 : const std::vector< tk::real >& vartimes() const noexcept
345 : 993 : { return m_vartimes; }
346 : 81 : std::vector< tk::real >& vartimes() noexcept { return m_vartimes; }
347 : : ///@}
348 : :
349 : : /** @name Node variables accessors */
350 : : ///@{
351 : 1023 : const std::vector< std::vector< std::vector< tk::real > > >& nodevars()
352 : 1023 : const noexcept { return m_nodevars; }
353 : 31909 : std::vector< std::vector< std::vector< tk::real > > >& nodevars() noexcept
354 : 31909 : { return m_nodevars; }
355 : : ///@}
356 : :
357 : : private:
358 : : //! Number of nodes
359 : : //! \details Stores the size (number of nodes) of the mesh graph.
360 : : //! Used if only the graph is needed but not the node coordinates, e.g.,
361 : : //! for graph partitioning, in which case only the connectivity is
362 : : //! required. If the coordinates are also loaded, the member functions
363 : : //! nnode() and size() return the same.
364 : : std::size_t m_graphsize;
365 : :
366 : : //! Element connectivity
367 : : std::vector< std::size_t > m_triinpoel; //!< Triangle connectivity
368 : : std::vector< std::size_t > m_tetinpoel; //!< Tetrahedron connectivity
369 : :
370 : : //! Node coordinates
371 : : std::vector< real > m_x;
372 : : std::vector< real > m_y;
373 : : std::vector< real > m_z;
374 : :
375 : : //! Side sets storing face ids adjacent to side sets
376 : : //! \details This stores lists of element IDs adjacent to faces associated
377 : : //! to side set IDs.
378 : : //! \note This is what ExodusII calls side set elem list.
379 : : std::map< int, std::vector< std::size_t > > m_bface;
380 : :
381 : : //! Side sets storing node ids adjacent to side sets
382 : : //! \details This stores lists of node IDs adjacent to faces associated
383 : : //! to side set IDs.
384 : : std::map< int, std::vector< std::size_t > > m_bnode;
385 : :
386 : : //! \brief Sides of faces used to define which face of an element is
387 : : //! adjacent to side set associated to side set id.
388 : : //! \note This is what ExodusII calls side set side list.
389 : : std::map< int, std::vector< std::size_t > > m_faceid;
390 : :
391 : : //! Node field data names
392 : : std::vector< std::string > m_nodevarnames;
393 : :
394 : : //! Time values for node field data
395 : : std::vector< tk::real > m_vartimes;
396 : :
397 : : //! Node field data
398 : : std::vector< std::vector< std::vector< tk::real > > > m_nodevars;
399 : :
400 : : //! Compute and return number of unique nodes in element connectivity
401 : : //! \param[in] inpoel Element connectivity
402 : : //! \return Number of unique node ids in connectivity, i.e., the graphsize
403 : : std::size_t
404 : 976 : graphsize( const std::vector< std::size_t >& inpoel ) {
405 [ + - ]: 976 : auto conn = inpoel;
406 [ + - ]: 976 : unique( conn );
407 : 1952 : return conn.size();
408 : 976 : }
409 : : };
410 : :
411 : : } // tk::
412 : :
413 : : #endif // UnsMesh_h
|
