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