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 : 193439028 : std::size_t operator()( const std::array< std::size_t, N >& p ) const {
83 : : using highwayhash::SipHash;
84 : : Shaper< N > shaper;
85 [ + + ]: 593796935 : for (std::size_t i=0; i<N; ++i) shaper.sizets[i] = p[i];
86 [ + - ]: 193439028 : std::sort( std::begin(shaper.sizets), std::end(shaper.sizets) );
87 [ + - ]: 386878056 : 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 : 120110788 : bool operator()( const std::array< std::size_t, N >& l,
103 : : const std::array< std::size_t, N >& r ) const
104 : : {
105 : 120110788 : std::array< std::size_t, N > s = l, p = r;
106 [ + - ]: 120110788 : std::sort( begin(s), end(s) );
107 [ + - ]: 120110788 : std::sort( begin(p), end(p) );
108 [ + - ]: 240221576 : 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 : 925 : explicit UnsMesh( const std::vector< std::size_t >& tetinp,
168 : 925 : const Coords& coord ) :
169 : 925 : m_graphsize( graphsize( tetinp ) ),
170 : : m_lininpoel(), m_triinpoel(),
171 : : m_tetinpoel( tetinp ),
172 : 925 : m_x( coord[0] ),
173 : 925 : m_y( coord[1] ),
174 [ + - ][ + - ]: 925 : m_z( coord[2] )
[ + - ][ + - ]
175 : : {
176 [ - + ][ - - ]: 925 : Assert( m_tetinpoel.size()%4 == 0,
[ - - ][ - - ]
177 : : "Size of tetinpoel must be divisible by 4" );
178 : 925 : }
179 : :
180 : : //! \brief Constructor copying over triangle element connectivity and array
181 : : //! of point coordinates
182 : 13 : 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< std::string >& elemvarnames = {},
187 : : const std::vector< tk::real >& vartimes = {},
188 : : const std::vector< std::vector< std::vector< tk::real > > >&
189 : : nodevars = {},
190 : : const std::vector< std::vector< std::vector< tk::real > > >&
191 : 13 : elemvars = {} ) :
192 : 13 : m_graphsize( graphsize( triinp ) ),
193 : : m_lininpoel(),
194 : : m_triinpoel( triinp ),
195 : : m_tetinpoel(),
196 : 13 : m_x( coord[0] ),
197 : 13 : m_y( coord[1] ),
198 : 13 : m_z( coord[2] ),
199 : : m_nodevarnames( nodevarnames ),
200 : : m_elemvarnames( elemvarnames ),
201 : : m_vartimes( vartimes ),
202 : : m_nodevars( nodevars ),
203 [ + - ][ + - ]: 13 : m_elemvars( elemvars )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
204 : : {
205 [ - + ][ - - ]: 13 : Assert( m_triinpoel.size()%3 == 0,
[ - - ][ - - ]
206 : : "Size of triinpoel must be divisible by 3" );
207 : 13 : }
208 : :
209 : : //! Constructor swallowing element connectivity and point coordinates
210 : : explicit UnsMesh( std::vector< std::size_t >&& tetinp,
211 : : std::vector< real >&& X,
212 : : std::vector< real >&& Y,
213 : : std::vector< real >&& Z ) :
214 : : m_graphsize( graphsize( tetinp ) ),
215 : : m_lininpoel(), m_triinpoel(),
216 : : m_tetinpoel( std::move(tetinp) ),
217 : : m_x( std::move(X) ),
218 : : m_y( std::move(Y) ),
219 : : m_z( std::move(Z) )
220 : : {
221 : : Assert( m_tetinpoel.size()%4 == 0,
222 : : "Size of tetinpoel must be divisible by 4" );
223 : : }
224 : :
225 : : //! \brief Constructor swallowing element connectivity and array of point
226 : : //! coordinates
227 : : explicit UnsMesh( std::vector< std::size_t >&& tetinp, Coords&& coord ) :
228 : : m_graphsize( graphsize( tetinp ) ),
229 : : m_lininpoel(), m_triinpoel(),
230 : : m_tetinpoel( std::move(tetinp) ),
231 : : m_x( std::move(coord[0]) ),
232 : : m_y( std::move(coord[1]) ),
233 : : m_z( std::move(coord[2]) )
234 : : {
235 : : Assert( m_tetinpoel.size()%4 == 0,
236 : : "Size of tetinpoel must be divisible by 4" );
237 : : }
238 : :
239 : : //! Constructor with connectivities and side set faces
240 : 53 : explicit UnsMesh(
241 : : const std::vector< std::size_t >& tetinp,
242 : : const Coords& coord,
243 : : const std::map< int, std::vector< std::size_t > >& bf,
244 : : const std::vector< std::size_t >& triinp,
245 : 53 : const std::map< int, std::vector< std::size_t > >& fid ) :
246 : 53 : m_graphsize( graphsize( tetinp ) ),
247 : : m_lininpoel(),
248 : : m_triinpoel( triinp ),
249 : : m_tetinpoel( tetinp ),
250 : 53 : m_x( coord[0] ),
251 : 53 : m_y( coord[1] ),
252 : 53 : m_z( coord[2] ),
253 : : m_bface( bf ),
254 [ + - ][ + - ]: 53 : m_faceid( fid )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
255 : : {
256 [ - + ][ - - ]: 53 : Assert( m_tetinpoel.size() % 4 == 0,
[ - - ][ - - ]
257 : : "Size of tetinpoel must be divisible by 4" );
258 [ - + ][ - - ]: 53 : Assert( m_triinpoel.size() % 3 == 0,
[ - - ][ - - ]
259 : : "Size of triinpoel must be divisible by 3" );
260 : 53 : }
261 : :
262 : : //! Constructor with connectivities and side set nodes
263 : 42 : explicit UnsMesh(
264 : : const std::vector< std::size_t >& tetinp,
265 : : const Coords& coord,
266 : 42 : const std::map< int, std::vector< std::size_t > >& bn ) :
267 : 42 : m_graphsize( graphsize( tetinp ) ),
268 : : m_lininpoel(),
269 : : m_triinpoel(),
270 : : m_tetinpoel( tetinp ),
271 : 42 : m_x( coord[0] ),
272 : 42 : m_y( coord[1] ),
273 : 42 : m_z( coord[2] ),
274 [ + - ][ + - ]: 42 : m_bnode( bn )
[ + - ][ + - ]
[ + - ]
275 : : {
276 [ - + ][ - - ]: 42 : Assert( m_tetinpoel.size() % 4 == 0,
[ - - ][ - - ]
277 : : "Size of tetinpoel must be divisible by 4" );
278 : 42 : }
279 : : ///@}
280 : :
281 : : /** @name Point coordinates accessors */
282 : : ///@{
283 : 37946 : const std::vector< real >& x() const noexcept { return m_x; }
284 : 37946 : const std::vector< real >& y() const noexcept { return m_y; }
285 : 37946 : const std::vector< real >& z() const noexcept { return m_z; }
286 : 80872 : std::vector< real >& x() noexcept { return m_x; }
287 : 80872 : std::vector< real >& y() noexcept { return m_y; }
288 : 80872 : std::vector< real >& z() noexcept { return m_z; }
289 : : ///@}
290 : :
291 : : /** @name Number of nodes accessors */
292 : : ///@{
293 : 37968 : std::size_t nnode() const noexcept { return m_x.size(); }
294 : : std::size_t nnode() noexcept { return m_x.size(); }
295 : : ///@}
296 : :
297 : : /** @name Graph size accessors */
298 : : ///@{
299 : : const std::size_t& size() const noexcept { return m_graphsize; }
300 : 36 : std::size_t& size() noexcept { return m_graphsize; }
301 : : ///@}
302 : :
303 : : //! Total number of elements accessor
304 : : std::size_t nelem() const noexcept {
305 : : return m_lininpoel.size()/2 + m_triinpoel.size()/3 + m_tetinpoel.size()/4;
306 : : }
307 : :
308 : : //! Number of element blocks accessor
309 : 1043 : std::size_t neblk() const noexcept {
310 : 1043 : return !m_triinpoel.empty() + !m_tetinpoel.empty();
311 : : }
312 : :
313 : : /** @name Line elements connectivity accessors */
314 : : ///@{
315 : 12 : const std::vector< std::size_t >& lininpoel() const noexcept
316 : 12 : { return m_lininpoel; }
317 : 7 : std::vector< std::size_t >& lininpoel() noexcept { return m_lininpoel; }
318 : : ///@}
319 : :
320 : : /** @name Triangle elements connectivity accessors */
321 : : ///@{
322 : 24841 : const std::vector< std::size_t >& triinpoel() const noexcept
323 : 24841 : { return m_triinpoel; }
324 : 202196 : std::vector< std::size_t >& triinpoel() noexcept { return m_triinpoel; }
325 : : ///@}
326 : :
327 : : /** @name Tetrahedra elements connectivity accessors */
328 : : ///@{
329 : 377502 : const std::vector< std::size_t >& tetinpoel() const noexcept
330 : 377502 : { return m_tetinpoel; }
331 : 1919372 : std::vector< std::size_t >& tetinpoel() noexcept { return m_tetinpoel; }
332 : : ///@}
333 : :
334 : : /** @name Side set face list accessors */
335 : : ///@{
336 : 2086 : const std::map< int, std::vector< std::size_t > >& bface() const noexcept
337 : 2086 : { return m_bface; }
338 : 901 : std::map< int, std::vector< std::size_t > >& bface() noexcept
339 : 901 : { return m_bface; }
340 : : ///@}
341 : :
342 : : /** @name Side set face id accessors */
343 : : ///@{
344 : 240 : const std::map< int, std::vector< std::size_t > >& faceid() const noexcept
345 : 240 : { return m_faceid; }
346 : 901 : std::map< int, std::vector< std::size_t > >& faceid() noexcept
347 : 901 : { return m_faceid; }
348 : : ///@}
349 : :
350 : : /** @name Side set node list accessors */
351 : : ///@{
352 : 2086 : const std::map< int, std::vector< std::size_t > >& bnode() const noexcept
353 : 2086 : { return m_bnode; }
354 : : std::map< int, std::vector< std::size_t > >& bnode() noexcept
355 : : { return m_bnode; }
356 : : ///@}
357 : :
358 : : /** @name Node variable names accessors */
359 : : ///@{
360 : 1043 : const std::vector< std::string >& nodevarnames() const noexcept
361 : 1043 : { return m_nodevarnames; }
362 : 81 : std::vector< std::string >& nodevarnames() noexcept
363 : 81 : { return m_nodevarnames; }
364 : : ///@}
365 : :
366 : : /** @name Element variable names accessors */
367 : : ///@{
368 : 1043 : const std::vector< std::string >& elemvarnames() const noexcept
369 : 1043 : { return m_elemvarnames; }
370 : 81 : std::vector< std::string >& elemvarnames() noexcept
371 : 81 : { return m_elemvarnames; }
372 : : ///@}
373 : :
374 : : /** @name Variable times accessors */
375 : : ///@{
376 : 1053 : const std::vector< tk::real >& vartimes() const noexcept
377 : 1053 : { return m_vartimes; }
378 : 117 : std::vector< tk::real >& vartimes() noexcept { return m_vartimes; }
379 : : ///@}
380 : :
381 : : /** @name Node variables accessors */
382 : : ///@{
383 : 1083 : const std::vector< std::vector< std::vector< tk::real > > >& nodevars()
384 : 1083 : const noexcept { return m_nodevars; }
385 : 31909 : std::vector< std::vector< std::vector< tk::real > > >& nodevars() noexcept
386 : 31909 : { return m_nodevars; }
387 : : ///@}
388 : :
389 : : /** @name Element variables accessors */
390 : : ///@{
391 : 1073 : const std::vector< std::vector< std::vector< tk::real > > >& elemvars()
392 : 1073 : const noexcept { return m_elemvars; }
393 : 63720 : std::vector< std::vector< std::vector< tk::real > > >& elemvars() noexcept
394 : 63720 : { return m_elemvars; }
395 : : ///@}
396 : :
397 : : private:
398 : : //! Number of nodes
399 : : //! \details Stores the size (number of nodes) of the mesh graph.
400 : : //! Used if only the graph is needed but not the node coordinates, e.g.,
401 : : //! for graph partitioning, in which case only the connectivity is
402 : : //! required. If the coordinates are also loaded, the member functions
403 : : //! nnode() and size() return the same.
404 : : std::size_t m_graphsize;
405 : :
406 : : //! Element connectivity
407 : : std::vector< std::size_t > m_lininpoel; //!< Line connectivity
408 : : std::vector< std::size_t > m_triinpoel; //!< Triangle connectivity
409 : : std::vector< std::size_t > m_tetinpoel; //!< Tetrahedron connectivity
410 : :
411 : : //! Node coordinates
412 : : std::vector< real > m_x;
413 : : std::vector< real > m_y;
414 : : std::vector< real > m_z;
415 : :
416 : : //! Side sets storing face ids adjacent to side sets
417 : : //! \details This stores lists of element IDs adjacent to faces associated
418 : : //! to side set IDs.
419 : : //! \note This is what ExodusII calls side set elem list.
420 : : std::map< int, std::vector< std::size_t > > m_bface;
421 : :
422 : : //! Side sets storing node ids adjacent to side sets
423 : : //! \details This stores lists of node IDs adjacent to faces associated
424 : : //! to side set IDs.
425 : : std::map< int, std::vector< std::size_t > > m_bnode;
426 : :
427 : : //! \brief Sides of faces used to define which face of an element is
428 : : //! adjacent to side set associated to side set id.
429 : : //! \note This is what ExodusII calls side set side list.
430 : : std::map< int, std::vector< std::size_t > > m_faceid;
431 : :
432 : : //! Node field data names
433 : : std::vector< std::string > m_nodevarnames;
434 : :
435 : : //! Element field data names
436 : : std::vector< std::string > m_elemvarnames;
437 : :
438 : : //! Time values for node field data
439 : : std::vector< tk::real > m_vartimes;
440 : :
441 : : //! Node field data
442 : : std::vector< std::vector< std::vector< tk::real > > > m_nodevars;
443 : :
444 : : //! Element field data
445 : : std::vector< std::vector< std::vector< tk::real > > > m_elemvars;
446 : :
447 : : //! Compute and return number of unique nodes in element connectivity
448 : : //! \param[in] inpoel Element connectivity
449 : : //! \return Number of unique node ids in connectivity, i.e., the graphsize
450 : : std::size_t
451 : 1037 : graphsize( const std::vector< std::size_t >& inpoel ) {
452 [ + - ]: 2074 : auto conn = inpoel;
453 [ + - ]: 1037 : unique( conn );
454 : 2074 : return conn.size();
455 : : }
456 : : };
457 : :
458 : : } // tk::
459 : :
460 : : #endif // UnsMesh_h
|
