Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Mesh/DerivedData.cpp
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 Generate data structures derived from unstructured mesh
9 : : \details Generate data structures derived from the connectivity information
10 : : of an unstructured mesh.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include <set>
15 : : #include <map>
16 : : #include <iterator>
17 : : #include <numeric>
18 : : #include <algorithm>
19 : : #include <type_traits>
20 : : #include <cstddef>
21 : : #include <array>
22 : : #include <unordered_set>
23 : : #include <unordered_map>
24 : : #include <iostream>
25 : : #include <cfenv>
26 : :
27 : : #include "Exception.hpp"
28 : : #include "DerivedData.hpp"
29 : : #include "ContainerUtil.hpp"
30 : : #include "Vector.hpp"
31 : :
32 : : namespace tk {
33 : :
34 : : int
35 : 491137308 : orient( const UnsMesh::Edge& t, const UnsMesh::Edge& e )
36 : : // *****************************************************************************
37 : : // Determine edge orientation
38 : : //! \return -1.0 if edge t is oriented the same as edge e, +1.0 opposite
39 : : // *****************************************************************************
40 : : {
41 [ + + ][ + + ]: 491137308 : if (t[0] == e[0] && t[1] == e[1])
[ + + ]
42 : 45682003 : return -1;
43 [ + + ][ + + ]: 445455305 : else if (t[0] == e[1] && t[1] == e[0])
[ + + ]
44 : 36174215 : return 1;
45 : : else
46 : 409281090 : return 0;
47 : : }
48 : :
49 : : std::size_t
50 : 2128 : npoin_in_graph( const std::vector< std::size_t >& inpoel )
51 : : // *****************************************************************************
52 : : // Compute number of points (nodes) in mesh from connectivity
53 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
54 : : //! \return Number of mesh points (nodes)
55 : : // *****************************************************************************
56 : : {
57 [ + - ]: 2128 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
58 [ - + ][ - - ]: 2128 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
59 : 2128 : return *minmax.second + 1;
60 : : }
61 : :
62 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
63 : 46970 : genEsup( const std::vector< std::size_t >& inpoel, std::size_t nnpe )
64 : : // *****************************************************************************
65 : : // Generate derived data structure, elements surrounding points
66 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
67 : : //! node ids of each element of an unstructured mesh. Example:
68 : : //! \code{.cpp}
69 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
70 : : //! 10, 14, 13, 12 };
71 : : //! \endcode
72 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
73 : : //! and { 10, 14, 13, 12 }.
74 : : //! \param[in] nnpe Number of nodes per element
75 : : //! \return Linked lists storing elements surrounding points
76 : : //! \warning It is not okay to call this function with an empty container or a
77 : : //! non-positive number of nodes per element; it will throw an exception.
78 : : //! \details The data generated here is stored in a linked list, more precisely,
79 : : //! two linked arrays (vectors), _esup1_ and _esup2_, where _esup2_ holds the
80 : : //! indices at which _esup1_ holds the element ids surrounding points. Looping
81 : : //! over all elements surrounding all points can then be accomplished by the
82 : : //! following loop:
83 : : //! \code{.cpp}
84 : : //! for (std::size_t p=0; p<npoin; ++p)
85 : : //! for (auto i=esup.second[p]+1; i<=esup.second[p+1]; ++i)
86 : : //! use element id esup.first[i]
87 : : //! \endcode
88 : : //! To find out the number of points, _npoin_, the mesh connectivity,
89 : : //! _inpoel_, can be queried:
90 : : //! \code{.cpp}
91 : : //! auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
92 : : //! Assert( *minmax.first == 0, "node ids should start from zero" );
93 : : //! auto npoin = *minmax.second + 1;
94 : : //! \endcode
95 : : //! \note In principle, this function *should* work for any positive nnpe,
96 : : //! however, only nnpe = 4 (tetrahedra) and nnpe = 3 (triangles) are tested.
97 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
98 : : // *****************************************************************************
99 : : {
100 [ + + ][ + - ]: 46970 : Assert( !inpoel.empty(), "Attempt to call genEsup() on empty container" );
[ + - ][ + - ]
101 [ + + ][ + - ]: 46969 : Assert( nnpe > 0, "Attempt to call genEsup() with zero nodes per element" );
[ + - ][ + - ]
102 [ + + ][ + - ]: 46968 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ + - ][ + - ]
103 : :
104 : : // find out number of points in mesh connectivity
105 [ + - ]: 46961 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
106 [ - + ][ - - ]: 46961 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
107 : 46961 : auto npoin = *minmax.second + 1;
108 : :
109 : : // allocate one of the linked lists storing elements surrounding points: esup2
110 : : // fill with zeros
111 [ + - ]: 93922 : std::vector< std::size_t > esup2( npoin+1, 0 );
112 : :
113 : : // element pass 1: count number of elements connected to each point
114 [ + + ]: 67817065 : for (auto n : inpoel) ++esup2[ n + 1 ];
115 : :
116 : : // storage/reshuffling pass 1: update storage counter and store
117 : : // also find out the maximum size of esup1 (mesup)
118 : 46961 : auto mesup = esup2[0]+1;
119 [ + + ]: 5893621 : for (std::size_t i=1; i<npoin+1; ++i) {
120 : 5846660 : esup2[i] += esup2[i-1];
121 [ + - ]: 5846660 : if (esup2[i]+1 > mesup) mesup = esup2[i]+1;
122 : : }
123 : :
124 : : // now we know mesup, so allocate the other one of the linked lists storing
125 : : // elements surrounding points: esup1
126 [ + - ]: 93922 : std::vector< std::size_t > esup1( mesup );
127 : :
128 : : // store the elements in esup1
129 : 46961 : std::size_t e = 0;
130 [ + + ]: 67817065 : for (auto n : inpoel) {
131 : 67770104 : auto j = esup2[n]+1;
132 : 67770104 : esup2[n] = j;
133 : 67770104 : esup1[j] = e/nnpe;
134 : 67770104 : ++e;
135 : : }
136 : :
137 : : // storage/reshuffling pass 2
138 [ + + ]: 5893621 : for (auto i=npoin; i>0; --i) esup2[i] = esup2[i-1];
139 : 46961 : esup2[0] = 0;
140 : :
141 : : // Return (move out) linked lists
142 [ + - ]: 93922 : return std::make_pair( std::move(esup1), std::move(esup2) );
143 : : }
144 : :
145 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
146 : 3318 : genPsup( const std::vector< std::size_t >& inpoel,
147 : : std::size_t nnpe,
148 : : const std::pair< std::vector< std::size_t >,
149 : : std::vector< std::size_t > >& esup )
150 : : // *****************************************************************************
151 : : // Generate derived data structure, points surrounding points
152 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
153 : : //! node ids of each element of an unstructured mesh. Example:
154 : : //! \code{.cpp}
155 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
156 : : //! 10, 14, 13, 12 };
157 : : //! \endcode
158 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
159 : : //! and { 10, 14, 13, 12 }.
160 : : //! \param[in] nnpe Number of nodes per element
161 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
162 : : //! \return Linked lists storing points surrounding points
163 : : //! \warning It is not okay to call this function with an empty container for
164 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
165 : : //! element; it will throw an exception.
166 : : //! \details The data generated here is stored in a linked list, more precisely,
167 : : //! two linked arrays (vectors), _psup1_ and _psup2_, where _psup2_ holds the
168 : : //! indices at which _psup1_ holds the point ids surrounding points. Looping
169 : : //! over all points surrounding all points can then be accomplished by the
170 : : //! following loop:
171 : : //! \code{.cpp}
172 : : //! for (std::size_t p=0; p<npoin; ++p)
173 : : //! for (auto i=psup.second[p]+1; i<=psup.second[p+1]; ++i)
174 : : //! use point id psup.first[i]
175 : : //! \endcode
176 : : //! To find out the number of points, _npoin_, the mesh connectivity,
177 : : //! _inpoel_, can be queried:
178 : : //! \code{.cpp}
179 : : //! auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
180 : : //! Assert( *minmax.first == 0, "node ids should start from zero" );
181 : : //! auto npoin = *minmax.second + 1;
182 : : //! \endcode
183 : : //! or the length-1 of the generated index list:
184 : : //! \code{.cpp}
185 : : //! auto npoin = psup.second.size()-1;
186 : : //! \endcode
187 : : //! \note In principle, this function *should* work for any positive nnpe,
188 : : //! however, only nnpe = 4 (tetrahedra) and nnpe = 3 (triangles) are tested.
189 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
190 : : // *****************************************************************************
191 : : {
192 [ + + ][ + - ]: 3318 : Assert( !inpoel.empty(), "Attempt to call genPsup() on empty container" );
[ + - ][ + - ]
193 [ + + ][ + - ]: 3317 : Assert( nnpe > 0, "Attempt to call genPsup() with zero nodes per element" );
[ + - ][ + - ]
194 [ - + ][ - - ]: 3316 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
195 [ + + ][ + - ]: 3316 : Assert( !esup.first.empty(), "Attempt to call genPsup() with empty esup1" );
[ + - ][ + - ]
196 [ + + ][ + - ]: 3315 : Assert( !esup.second.empty(), "Attempt to call genPsup() with empty esup2" );
[ + - ][ + - ]
197 : :
198 : : // find out number of points in mesh connectivity
199 [ + - ]: 3314 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
200 [ - + ][ - - ]: 3314 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
201 : 3314 : auto npoin = *minmax.second + 1;
202 : :
203 : 3314 : auto& esup1 = esup.first;
204 : 3314 : auto& esup2 = esup.second;
205 : :
206 : : // allocate both of the linked lists storing points surrounding points, we
207 : : // only know the size of psup2, put in a single zero in psup1
208 [ + - ][ + - ]: 9942 : std::vector< std::size_t > psup2( npoin+1 ), psup1( 1, 0 );
209 : :
210 : : // allocate and fill with zeros a temporary array, only used locally
211 [ + - ]: 6628 : std::vector< std::size_t > lpoin( npoin, 0 );
212 : :
213 : : // fill both psup1 and psup2
214 : 3314 : psup2[0] = 0;
215 : 3314 : std::size_t j = 0;
216 [ + + ]: 458859 : for (std::size_t p=0; p<npoin; ++p) {
217 [ + + ]: 6389357 : for (std::size_t i=esup2[p]+1; i<=esup2[p+1]; ++i ) {
218 [ + + ]: 29668844 : for (std::size_t n=0; n<nnpe; ++n) {
219 : 23735032 : auto q = inpoel[ esup1[i] * nnpe + n ];
220 [ + + ][ + + ]: 23735032 : if (q != p && lpoin[q] != p+1) {
[ + + ]
221 : 4564524 : ++j;
222 [ + - ]: 4564524 : psup1.push_back( q );
223 : 4564524 : lpoin[q] = p+1;
224 : : }
225 : : }
226 : : }
227 : 455545 : psup2[p+1] = j;
228 : : }
229 : :
230 : : // sort point ids for each point in psup1
231 [ + + ]: 458859 : for (std::size_t p=0; p<npoin; ++p)
232 [ + - ][ + - ]: 911090 : std::sort(
[ + - ]
233 : 455545 : std::next( begin(psup1), static_cast<std::ptrdiff_t>(psup2[p]+1) ),
234 : 455545 : std::next( begin(psup1), static_cast<std::ptrdiff_t>(psup2[p+1]+1) ) );
235 : :
236 : : // Return (move out) linked lists
237 [ + - ]: 6628 : return std::make_pair( std::move(psup1), std::move(psup2) );
238 : : }
239 : :
240 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
241 : 9 : genEdsup( const std::vector< std::size_t >& inpoel,
242 : : std::size_t nnpe,
243 : : const std::pair< std::vector< std::size_t >,
244 : : std::vector< std::size_t > >& esup )
245 : : // *****************************************************************************
246 : : // Generate derived data structure, edges surrounding points
247 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
248 : : //! node ids of each element of an unstructured mesh. Example:
249 : : //! \code{.cpp}
250 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
251 : : //! 10, 14, 13, 12 };
252 : : //! \endcode
253 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
254 : : //! and { 10, 14, 13, 12 }.
255 : : //! \param[in] nnpe Number of nodes per element (3 or 4)
256 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
257 : : //! \return Linked lists storing edges (point ids p < q) emanating from points
258 : : //! \warning It is not okay to call this function with an empty container for
259 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
260 : : //! element; it will throw an exception.
261 : : //! \details The data generated here is stored in a linked list, more precisely,
262 : : //! two linked arrays (vectors), _edsup1_ and _edsup2_, where _edsup2_ holds
263 : : //! the indices at which _edsup1_ holds the edge-end point ids emanating from
264 : : //! points for all points. The generated data structure, linked lists edsup1
265 : : //! and edsup2, are very similar to psup1 and psup2, generated by genPsup(),
266 : : //! except here only unique edges are stored, i.e., for edges with point ids
267 : : //! p < q, only ids q are stored that are still associated to point p. Looping
268 : : //! over all unique edges can then be accomplished by the following loop:
269 : : //! \code{.cpp}
270 : : //! for (std::size_t p=0; p<npoin; ++p)
271 : : //! for (auto i=edsup.second[p]+1; i<=edsup.second[p+1]; ++i)
272 : : //! use edge with point ids p < edsup.first[i]
273 : : //! \endcode
274 : : //! To find out the number of points, _npoin_, the mesh connectivity,
275 : : //! _inpoel_, can be queried:
276 : : //! \code{.cpp}
277 : : //! auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
278 : : //! Assert( *minmax.first == 0, "node ids should start from zero" );
279 : : //! auto npoin = *minmax.second + 1;
280 : : //! \endcode
281 : : //! \note At first sight, this function seems to work for elements with more
282 : : //! vertices than that of tetrahedra. However, that is not the case since the
283 : : //! algorithm for nnpe > 4 would erronously identify any two combination of
284 : : //! vertices as a valid edge of an element. Since only triangles and
285 : : //! tetrahedra have no internal edges, this algorithm only works for triangle
286 : : //! and tetrahedra element connectivity.
287 : : //! \see tk::genInpoed for similar data that sometimes may be more advantageous
288 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
289 : : // *****************************************************************************
290 : : {
291 [ + + ][ + - ]: 9 : Assert( !inpoel.empty(), "Attempt to call genEdsup() on empty container" );
[ + - ][ + - ]
292 [ + + ][ + - ]: 8 : Assert( nnpe > 0, "Attempt to call genEdsup() with zero nodes per element" );
[ + - ][ + - ]
293 [ + + ][ + + ]: 7 : Assert( nnpe == 3 || nnpe == 4,
[ + - ][ + - ]
[ + - ]
294 : : "Attempt to call genEdsup() with nodes per element, nnpe, that is "
295 : : "neither 4 (tetrahedra) nor 3 (triangles)." );
296 [ - + ][ - - ]: 6 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
297 [ + + ][ + - ]: 6 : Assert( !esup.first.empty(), "Attempt to call genEdsup() with empty esup1" );
[ + - ][ + - ]
298 [ + + ][ + - ]: 5 : Assert( !esup.second.empty(), "Attempt to call genEdsup() with empty esup2" );
[ + - ][ + - ]
299 : :
300 : : // find out number of points in mesh connectivity
301 [ + - ]: 4 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
302 [ - + ][ - - ]: 4 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
303 : 4 : auto npoin = *minmax.second + 1;
304 : :
305 : 4 : auto& esup1 = esup.first;
306 : 4 : auto& esup2 = esup.second;
307 : :
308 : : // allocate and fill with zeros a temporary array, only used locally
309 [ + - ]: 8 : std::vector< std::size_t > lpoin( npoin, 0 );
310 : :
311 : : // map to contain stars, a point associated to points connected with edges
312 : : // storing only the end-point id, q, of point ids p < q
313 : 8 : std::map< std::size_t, std::vector< std::size_t > > star;
314 : :
315 : : // generate edge connectivity and store as stars where center id < spike id
316 [ + + ]: 60 : for (std::size_t p=0; p<npoin; ++p)
317 [ + + ]: 392 : for (std::size_t i=esup2[p]+1; i<=esup2[p+1]; ++i )
318 [ + + ]: 1536 : for (std::size_t n=0; n<nnpe; ++n) {
319 : 1200 : auto q = inpoel[ esup1[i] * nnpe + n ];
320 [ + + ][ + + ]: 1200 : if (q != p && lpoin[q] != p+1) {
[ + + ]
321 [ + + ][ + - ]: 340 : if (p < q) star[p].push_back(q);
[ + - ]
322 : 340 : lpoin[q] = p+1;
323 : : }
324 : : }
325 : :
326 : : // linked lists (vectors) to store edges surrounding points and their indices
327 [ + - ][ + - ]: 8 : std::vector< std::size_t > edsup1( 1, 0 ), edsup2( 1, 0 );
328 : :
329 : : // sort non-center points of each star and store nodes and indices in vectors
330 [ + + ]: 46 : for (auto& p : star) {
331 [ + - ]: 42 : std::sort( begin(p.second), end(p.second) );
332 [ + - ]: 42 : edsup2.push_back( edsup2.back() + p.second.size() );
333 [ + - ]: 42 : edsup1.insert( end(edsup1), begin(p.second), end(p.second) );
334 : : }
335 : : // fill up index array with the last index for points with no new edges
336 [ + + ]: 18 : for (std::size_t i=0; i<npoin-star.size(); ++i)
337 [ + - ]: 14 : edsup2.push_back( edsup2.back() );
338 : :
339 : : // Return (move out) linked lists
340 [ + - ]: 8 : return std::make_pair( std::move(edsup1), std::move(edsup2) );
341 : : }
342 : :
343 : : std::vector< std::size_t >
344 : 14 : genInpoed( const std::vector< std::size_t >& inpoel,
345 : : std::size_t nnpe,
346 : : const std::pair< std::vector< std::size_t >,
347 : : std::vector< std::size_t > >& esup )
348 : : // *****************************************************************************
349 : : // Generate derived data structure, edge connectivity
350 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
351 : : //! node ids of each element of an unstructured mesh. Example:
352 : : //! \code{.cpp}
353 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
354 : : //! 10, 14, 13, 12 };
355 : : //! \endcode
356 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
357 : : //! and { 10, 14, 13, 12 }.
358 : : //! \param[in] nnpe Number of nodes per element (3 or 4)
359 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
360 : : //! \return Linear vector storing edge connectivity (point ids p < q)
361 : : //! \warning It is not okay to call this function with an empty container for
362 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
363 : : //! element; it will throw an exception.
364 : : //! \details The data generated here is stored in a linear vector and is very
365 : : //! similar to the linked lists, _edsup1_ and _edsup2, generated by
366 : : //! genEdsup(). The difference is that in the linear vector, inpoed, generated
367 : : //! here, both edge point ids are stored as a pair, p < q, as opposed to the
368 : : //! linked lists edsup1 and edsup2, in which edsup1 only stores the edge-end
369 : : //! point ids (still associated to edge-start point ids when used together
370 : : //! with edsup2). The rationale is that while inpoed is larger in memory, it
371 : : //! allows direct access to edges (pair of point ids making up an edge),
372 : : //! edsup1 and edsup2 are smaller in memory, still allow accessing the same
373 : : //! data (edge point id pairs) but only in a linear fashion, not by direct
374 : : //! access to particular edges. Accessing all unique edges using the edge
375 : : //! connectivity data structure, inpoed, generated here can be accomplished by
376 : : //! \code{.cpp}
377 : : //! for (std::size_t e=0; e<inpoed.size()/2; ++e) {
378 : : //! use point id p of edge e = inpoed[e*2];
379 : : //! use point id q of edge e = inpoed[e*2+1];
380 : : //! }
381 : : //! \endcode
382 : : //! \note At first sight, this function seems to work for elements with more
383 : : //! vertices than that of tetrahedra. However, that is not the case since the
384 : : //! algorithm for nnpe > 4 would erronously identify any two combination of
385 : : //! vertices as a valid edge of an element. Since only triangles and
386 : : //! tetrahedra have no internal edges, this algorithm only works for triangle
387 : : //! and tetrahedra element connectivity.
388 : : //! \see tk::genEdsup for similar data that sometimes may be more advantageous
389 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
390 : : // *****************************************************************************
391 : : {
392 [ + + ][ + - ]: 14 : Assert( !inpoel.empty(), "Attempt to call genInpoed() on empty container" );
[ + - ][ + - ]
393 [ + + ][ + - ]: 13 : Assert( nnpe > 0, "Attempt to call genInpoed() with zero nodes per element" );
[ + - ][ + - ]
394 [ + + ][ + + ]: 12 : Assert( nnpe == 3 || nnpe == 4,
[ + - ][ + - ]
[ + - ]
395 : : "Attempt to call genInpoed() with nodes per element, nnpe, that is "
396 : : "neither 4 (tetrahedra) nor 3 (triangles)." );
397 [ - + ][ - - ]: 11 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
398 [ + + ][ + - ]: 11 : Assert( !esup.first.empty(), "Attempt to call genInpoed() with empty esup1" );
[ + - ][ + - ]
399 [ + + ][ + - ]: 10 : Assert( !esup.second.empty(),
[ + - ][ + - ]
400 : : "Attempt to call genInpoed() with empty esup2" );
401 : :
402 : : // find out number of points in mesh connectivity
403 [ + - ]: 9 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
404 [ - + ][ - - ]: 9 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
405 : 9 : auto npoin = *minmax.second + 1;
406 : :
407 : 9 : auto& esup1 = esup.first;
408 : 9 : auto& esup2 = esup.second;
409 : :
410 : : // allocate and fill with zeros a temporary array, only used locally
411 [ + - ]: 18 : std::vector< std::size_t > lpoin( npoin, 0 );
412 : :
413 : : // map to contain stars, a point associated to points connected with edges,
414 : : // storing only the end-point id, q, of point ids p < q
415 : 18 : std::map< std::size_t, std::vector< std::size_t > > star;
416 : :
417 : : // generate edge connectivity and store as stars where center id < spike id
418 [ + + ]: 105 : for (std::size_t p=0; p<npoin; ++p)
419 [ + + ]: 636 : for (std::size_t i=esup2[p]+1; i<=esup2[p+1]; ++i )
420 [ + + ]: 2556 : for (std::size_t n=0; n<nnpe; ++n) {
421 : 2016 : auto q = inpoel[ esup1[i] * nnpe + n ];
422 [ + + ][ + + ]: 2016 : if (q != p && lpoin[q] != p+1) {
[ + + ]
423 [ + + ][ + - ]: 572 : if (p < q) star[p].push_back( q );
[ + - ]
424 : 572 : lpoin[q] = p+1;
425 : : }
426 : : }
427 : :
428 : : // linear vector to store edge connectivity and their indices
429 : 9 : std::vector< std::size_t > inpoed;
430 : :
431 : : // sort non-center points of each star and store both start and end points of
432 : : // each star in linear vector
433 [ + + ]: 86 : for (auto& p : star) {
434 [ + - ]: 77 : std::sort( begin(p.second), end(p.second) );
435 [ + + ]: 363 : for (auto e : p.second) {
436 [ + - ]: 286 : inpoed.push_back( p.first );
437 [ + - ]: 286 : inpoed.push_back( e );
438 : : }
439 : : }
440 : :
441 : : // Return (move out) linear vector
442 : 18 : return inpoed;
443 : : }
444 : :
445 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
446 : 8 : genEsupel( const std::vector< std::size_t >& inpoel,
447 : : std::size_t nnpe,
448 : : const std::pair< std::vector< std::size_t >,
449 : : std::vector< std::size_t > >& esup )
450 : : // *****************************************************************************
451 : : // Generate derived data structure, elements surrounding points of elements
452 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
453 : : //! node ids of each element of an unstructured mesh. Example:
454 : : //! \code{.cpp}
455 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
456 : : //! 10, 14, 13, 12 };
457 : : //! \endcode
458 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
459 : : //! and { 10, 14, 13, 12 }.
460 : : //! \param[in] nnpe Number of nodes per element
461 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
462 : : //! \return Linked lists storing elements surrounding points of elements
463 : : //! \warning It is not okay to call this function with an empty container for
464 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
465 : : //! element; it will throw an exception.
466 : : //! \details The data generated here is stored in a linked list, more precisely,
467 : : //! two linked arrays (vectors), _esupel1_ and _esupel2_, where _esupel2_
468 : : //! holds the indices at which _esupel1_ holds the element ids surrounding
469 : : //! points of elements. Looping over all elements surrounding the points of
470 : : //! all elements can then be accomplished by the following loop:
471 : : //! \code{.cpp}
472 : : //! for (std::size_t e=0; e<nelem; ++e)
473 : : //! for (auto i=esupel.second[e]+1; i<=esupel.second[e+1]; ++i)
474 : : //! use element id esupel.first[i]
475 : : //! \endcode
476 : : //! To find out the number of elements, _nelem_, the size of the mesh
477 : : //! connectivity vector, _inpoel_, can be devided by the number of nodes per
478 : : //! elements, _nnpe_:
479 : : //! \code{.cpp}
480 : : //! auto nelem = inpoel.size()/nnpe;
481 : : //! \endcode
482 : : //! \note In principle, this function *should* work for any positive nnpe,
483 : : //! however, only nnpe = 4 (tetrahedra) and nnpe = 3 (triangles) are tested.
484 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
485 : : // *****************************************************************************
486 : : {
487 [ + + ][ + - ]: 8 : Assert( !inpoel.empty(), "Attempt to call genEsupel() on empty container" );
[ + - ][ + - ]
488 [ + + ][ + - ]: 7 : Assert( nnpe > 0, "Attempt to call genEsupel() with zero nodes per element" );
[ + - ][ + - ]
489 [ - + ][ - - ]: 6 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
490 [ + + ][ + - ]: 6 : Assert( !esup.first.empty(), "Attempt to call genEsupel() with empty esup1" );
[ + - ][ + - ]
491 [ + + ][ + - ]: 5 : Assert( !esup.second.empty(),
[ + - ][ + - ]
492 : : "Attempt to call genEsupel() with empty esup2" );
493 : :
494 : 4 : auto& esup1 = esup.first;
495 : 4 : auto& esup2 = esup.second;
496 : :
497 : : // linked lists storing elements surrounding points of elements, put in a
498 : : // single zero in both
499 [ + - ][ + - ]: 8 : std::vector< std::size_t > esupel2( 1, 0 ), esupel1( 1, 0 );
500 : :
501 : 4 : std::size_t e = 0;
502 : 8 : std::set< std::size_t > esuel;
503 [ + + ]: 340 : for (auto p : inpoel) { // loop over all points of all elements
504 : : // collect unique element ids of elements surrounding points of element
505 [ + + ][ + - ]: 2736 : for (auto i=esup2[p]+1; i<=esup2[p+1]; ++i) esuel.insert( esup1[i] );
506 [ + + ]: 336 : if (++e%nnpe == 0) { // when finished checking all nodes of element
507 : : // erase element whose surrounding elements are considered
508 [ + - ]: 96 : esuel.erase( e/nnpe-1 );
509 : : // store unique element ids in esupel1
510 [ + - ]: 96 : esupel1.insert( end(esupel1), begin(esuel), end(esuel) );
511 : : // store end-index for element used to address into esupel1
512 [ + - ]: 96 : esupel2.push_back( esupel2.back() + esuel.size() );
513 : 96 : esuel.clear();
514 : : }
515 : : }
516 : :
517 : : // Return (move out) linked lists
518 [ + - ]: 8 : return std::make_pair( std::move(esupel1), std::move(esupel2) );
519 : : }
520 : :
521 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
522 : 8 : genEsuel( const std::vector< std::size_t >& inpoel,
523 : : std::size_t nnpe,
524 : : const std::pair< std::vector< std::size_t >,
525 : : std::vector< std::size_t > >& esup )
526 : : // *****************************************************************************
527 : : // Generate derived data structure, elements surrounding elements
528 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
529 : : //! node ids of each element of an unstructured mesh. Example:
530 : : //! \code{.cpp}
531 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
532 : : //! 10, 14, 13, 12 };
533 : : //! \endcode
534 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
535 : : //! and { 10, 14, 13, 12 }.
536 : : //! \param[in] nnpe Number of nodes per element
537 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
538 : : //! \return Linked lists storing elements surrounding elements
539 : : //! \warning It is not okay to call this function with an empty container for
540 : : //! inpoel or esup.first or esup.second; it will throw an exception.
541 : : //! \details The data generated here is stored in a linked list, more precisely,
542 : : //! two linked arrays (vectors), _esuel1_ and _esuel2_, where _esuel2_ holds
543 : : //! the indices at which _esuel1_ holds the element ids surrounding elements.
544 : : //! Looping over elements surrounding elements can then be accomplished by the
545 : : //! following loop:
546 : : //! \code{.cpp}
547 : : //! for (std::size_t e=0; e<nelem; ++e)
548 : : //! for (auto i=esuel.second[e]+1; i<=esuel.second[e+1]; ++i)
549 : : //! use element id esuel.first[i]
550 : : //! \endcode
551 : : //! To find out the number of elements, _nelem_, the size of the mesh
552 : : //! connectivity vector, _inpoel_, can be devided by the number of nodes per
553 : : //! elements, _nnpe_:
554 : : //! \code{.cpp}
555 : : //! auto nelem = inpoel.size()/nnpe;
556 : : //! \endcode
557 : : //! \note In principle, this function *should* work for any positive nnpe,
558 : : //! however, only nnpe = 4 (tetrahedra) and nnpe = 3 (triangles) are tested.
559 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
560 : : // *****************************************************************************
561 : : {
562 [ + + ][ + - ]: 8 : Assert( !inpoel.empty(), "Attempt to call genEsuel() on empty container" );
[ + - ][ + - ]
563 [ + + ][ + - ]: 7 : Assert( nnpe > 0, "Attempt to call genEsuel() with zero nodes per element" );
[ + - ][ + - ]
564 [ - + ][ - - ]: 6 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by four" );
[ - - ][ - - ]
565 [ + + ][ + - ]: 6 : Assert( !esup.first.empty(), "Attempt to call genEsuel() with empty esuel1" );
[ + - ][ + - ]
566 [ + + ][ + - ]: 5 : Assert( !esup.second.empty(),
[ + - ][ + - ]
567 : : "Attempt to call genEsuel() with empty esuel2" );
568 : :
569 : 4 : auto& esup1 = esup.first;
570 : 4 : auto& esup2 = esup.second;
571 : :
572 : 4 : auto nelem = inpoel.size()/nnpe;
573 : :
574 : : // lambda that returns 1 if elements hel and gel share a face
575 : 123936 : auto adj = [ &inpoel, nnpe ]( std::size_t hel, std::size_t gel ) -> bool {
576 : 2400 : std::size_t sp = 0;
577 [ + + ]: 11232 : for (std::size_t h=0; h<nnpe; ++h)
578 [ + + ]: 41856 : for (std::size_t g=0; g<nnpe; ++g)
579 [ + + ]: 33024 : if (inpoel[hel*nnpe+h] == inpoel[gel*nnpe+g]) ++sp;
580 [ + + ]: 2400 : if (sp == nnpe-1) return true; else return false;
581 : 4 : };
582 : :
583 : : // map to associate unique elements and their surrounding elements
584 : 8 : std::map< std::size_t, std::vector< std::size_t > > es;
585 : :
586 [ + + ]: 100 : for (std::size_t e=0; e<nelem; ++e) {
587 : 192 : std::set< std::size_t > faces; // will collect elem ids of shared faces
588 [ + + ]: 432 : for (std::size_t n=0; n<nnpe; ++n) {
589 : 336 : auto i = inpoel[ e*nnpe+n ];
590 [ + + ]: 2736 : for (auto j=esup2[i]+1; j<=esup2[i+1]; ++j)
591 [ + + ][ + - ]: 2400 : if (adj( e, esup1[j] )) faces.insert( esup1[j] );
592 : : }
593 : : // store element ids of shared faces
594 [ + + ][ + - ]: 384 : for (auto j : faces) es[e].push_back(j);
[ + - ]
595 : : }
596 : :
597 : : // storing elements surrounding elements
598 [ + - ][ + - ]: 8 : std::vector< std::size_t > esuel1( 1, 0 ), esuel2( 1, 0 );
599 : :
600 : : // store elements surrounding elements in linked lists
601 [ + + ]: 100 : for (const auto& e : es) {
602 [ + - ]: 96 : esuel2.push_back( esuel2.back() + e.second.size() );
603 [ + - ]: 96 : esuel1.insert( end(esuel1), begin(e.second), end(e.second) );
604 : : }
605 : :
606 : : // Return (move out) linked lists
607 [ + - ]: 8 : return std::make_pair( std::move(esuel1), std::move(esuel2) );
608 : : }
609 : :
610 : : std::vector< std::size_t >
611 : 6 : genInedel( const std::vector< std::size_t >& inpoel,
612 : : std::size_t nnpe,
613 : : const std::vector< std::size_t >& inpoed )
614 : : // *****************************************************************************
615 : : // Generate derived data structure, edges of elements
616 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
617 : : //! node ids of each element of an unstructured mesh. Example:
618 : : //! \code{.cpp}
619 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
620 : : //! 10, 14, 13, 12 };
621 : : //! \endcode
622 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
623 : : //! and { 10, 14, 13, 12 }.
624 : : //! \param[in] nnpe Number of nodes per element
625 : : //! \param[in] inpoed Edge connectivity as linear vector, see tk::genInpoed
626 : : //! \return Linear vector storing all edge ids * 2 of all elements
627 : : //! \warning It is not okay to call this function with an empty container for
628 : : //! inpoel or inpoed or a non-positive number of nodes per element; it will
629 : : //! throw an exception.
630 : : //! \details The data generated here is stored in a linear vector with all
631 : : //! edge ids (as defined by inpoed) of all elements. The edge ids stored in
632 : : //! inedel can be directly used to index the vector inpoed. Because the
633 : : //! derived data structure generated here, inedel, is intended to be used in
634 : : //! conjunction with the linear vector inpoed and not with the linked lists
635 : : //! edsup1 and edsup2, this function takes inpoed as an argument. Accessing
636 : : //! the edges of element e using the edge of elements data structure, inedel,
637 : : //! generated here can be accomplished by
638 : : //! \code{.cpp}
639 : : //! for (std::size_t e=0; e<nelem; ++e) {
640 : : //! for (std::size_t i=0; i<nepe; ++i) {
641 : : //! use edge id inedel[e*nepe+i] of element e, or
642 : : //! use point ids p < q of edge id inedel[e*nepe+i] of element e as
643 : : //! p = inpoed[ inedel[e*nepe+i]*2 ]
644 : : //! q = inpoed[ inedel[e*nepe+i]*2+1 ]
645 : : //! }
646 : : //! }
647 : : //! \endcode
648 : : //! where _nepe_ denotes the number of edges per elements: 3 for triangles, 6
649 : : //! for tetrahedra. To find out the number of elements, _nelem_, the size of
650 : : //! the mesh connectivity vector, _inpoel_, can be devided by the number of
651 : : //! nodes per elements, _nnpe_:
652 : : //! \code{.cpp}
653 : : //! auto nelem = inpoel.size()/nnpe;
654 : : //! \endcode
655 : : //! \note At first sight, this function seems to work for elements with more
656 : : //! vertices than that of tetrahedra. However, that is not the case since the
657 : : //! algorithm for nnpe > 4 would erronously identify any two combination of
658 : : //! vertices as a valid edge of an element. Since only triangles and
659 : : //! tetrahedra have no internal edges, this algorithm only works for triangle
660 : : //! and tetrahedra element connectivity.
661 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
662 : : // *****************************************************************************
663 : : {
664 [ + + ][ + - ]: 6 : Assert( !inpoel.empty(), "Attempt to call genInedel() on empty container" );
[ + - ][ + - ]
665 [ + + ][ + - ]: 5 : Assert( nnpe > 0, "Attempt to call genInedel() with zero nodes per element" );
[ + - ][ + - ]
666 [ + + ][ + + ]: 3 : Assert( nnpe == 3 || nnpe == 4,
[ + - ][ + - ]
[ + - ]
667 : : "Attempt to call genInedel() with nodes per element, nnpe, that is "
668 : : "neither 4 (tetrahedra) nor 3 (triangles)." );
669 [ - + ][ - - ]: 2 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
670 [ - + ][ - - ]: 2 : Assert( !inpoed.empty(), "Attempt to call genInedel() with empty inpoed" );
[ - - ][ - - ]
671 : :
672 : : // find out number of points in mesh connectivity
673 [ + - ]: 2 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
674 [ - + ][ - - ]: 2 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
675 : 2 : auto npoin = *minmax.second + 1;
676 : :
677 : : // First, generate index of star centers. This is necessary to avoid a
678 : : // brute-force search for point ids of edges when searching for element edges.
679 : : // Note that this is the same as edsup2, generated by genEdsup(). However,
680 : : // because the derived data structure generated here, inedel, is intended to
681 : : // be used in conjunction with the linear vector inpoed and not with the
682 : : // linked lists edsup1 and edsup2, this function takes inpoed as an argument,
683 : : // and so edsup2 is temporarily generated here to avoid a brute-force search.
684 : :
685 : : // map to contain stars, a point associated to points connected with edges
686 : : // storing only the end-point id, q, of point ids p < q
687 : 4 : std::map< std::size_t, std::vector< std::size_t > > star;
688 : :
689 : : // generate stars from inpoed; starting with zero, every even is a star
690 : : // center, every odd is a spike
691 [ + + ]: 87 : for (std::size_t i=0; i<inpoed.size()/2; ++i)
692 [ + - ][ + - ]: 85 : star[ inpoed[i*2] ].push_back( inpoed[i*2+1] );
693 : :
694 : : // store index of star centers in vector; assume non-center points of each
695 : : // star have already been sorted
696 [ + - ]: 4 : std::vector< std::size_t > edsup2( 1, 0 );
697 [ + + ][ + - ]: 23 : for (const auto& p : star) edsup2.push_back(edsup2.back() + p.second.size());
698 : : // fill up index array with the last index for points with no new edges
699 [ + + ]: 9 : for (std::size_t i=0; i<npoin-star.size(); ++i)
700 [ + - ]: 7 : edsup2.push_back( edsup2.back() );
701 : 2 : star.clear();
702 : :
703 : : // Second, generate edges of elements
704 : :
705 : 2 : auto nelem = inpoel.size()/nnpe;
706 : :
707 : : // map associating elem id with vector of edge ids
708 : 4 : std::map< std::size_t, std::vector< std::size_t > > edges;
709 : :
710 : : // generate map of elements associated to edge ids
711 [ + + ]: 50 : for (std::size_t e=0; e<nelem; ++e)
712 [ + + ]: 216 : for (std::size_t n=0; n<nnpe; ++n) {
713 : 168 : auto p = inpoel[e*nnpe+n];
714 [ + + ]: 662 : for (auto i=edsup2[p]+1; i<=edsup2[p+1]; ++i)
715 [ + + ]: 2254 : for (std::size_t j=0; j<nnpe; ++j)
716 [ + + ]: 1760 : if (inpoed[(i-1)*2+1] == inpoel[e*nnpe+j])
717 [ + - ][ + - ]: 216 : edges[e].push_back( i-1 );
718 : : }
719 : :
720 : : // linear vector to store the edge ids of all elements
721 [ + - ]: 2 : std::vector< std::size_t > inedel( sumvalsize(edges) );
722 : :
723 : : // store edge ids of elements in linear vector
724 : 2 : std::size_t j = 0;
725 [ + + ][ + + ]: 266 : for (const auto& e : edges) for (auto p : e.second) inedel[ j++ ] = p;
726 : :
727 : : // Return (move out) vector
728 : 4 : return inedel;
729 : : }
730 : :
731 : : std::unordered_map< UnsMesh::Edge, std::vector< std::size_t >,
732 : : UnsMesh::Hash<2>, UnsMesh::Eq<2> >
733 : 32883 : genEsued( const std::vector< std::size_t >& inpoel,
734 : : std::size_t nnpe,
735 : : const std::pair< std::vector< std::size_t >,
736 : : std::vector< std::size_t > >& esup )
737 : : // *****************************************************************************
738 : : // Generate derived data structure, elements surrounding edges
739 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
740 : : //! node ids of each element of an unstructured mesh. Example:
741 : : //! \code{.cpp}
742 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
743 : : //! 10, 14, 13, 12 };
744 : : //! \endcode
745 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
746 : : //! and { 10, 14, 13, 12 }.
747 : : //! \param[in] nnpe Number of nodes per element (3 or 4)
748 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
749 : : //! \return Associative container storing elements surrounding edges (value),
750 : : //! assigned to edge-end points (key)
751 : : //! \warning It is not okay to call this function with an empty container for
752 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
753 : : //! element; it will throw an exception.
754 : : //! \details Looping over elements surrounding all edges can be accomplished by
755 : : //! the following loop:
756 : : //! \code{.cpp}
757 : : //! for (const auto& [edge,surr_elements] : esued) {
758 : : //! use element edge-end-point ids edge[0] and edge[1]
759 : : //! for (auto e : surr_elements) {
760 : : //! use element id e
761 : : //! }
762 : : //! }
763 : : //! \endcode
764 : : //! esued.size() equals the number of edges.
765 : : //! \note At first sight, this function seems to work for elements with more
766 : : //! vertices than that of tetrahedra. However, that is not the case since the
767 : : //! algorithm for nnpe > 4 would erronously identify any two combination of
768 : : //! vertices as a valid edge of an element. Since only triangles and
769 : : //! tetrahedra have no internal edges, this algorithm only works for triangle
770 : : //! and tetrahedra element connectivity.
771 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
772 : : // *****************************************************************************
773 : : {
774 [ + + ][ + - ]: 32883 : Assert( !inpoel.empty(), "Attempt to call genEsued() on empty container" );
[ + - ][ + - ]
775 [ + + ][ + - ]: 32882 : Assert( nnpe > 0, "Attempt to call genEsued() with zero nodes per element" );
[ + - ][ + - ]
776 [ + + ][ + + ]: 32881 : Assert( nnpe == 3 || nnpe == 4,
[ + - ][ + - ]
[ + - ]
777 : : "Attempt to call genEsued() with nodes per element, nnpe, that is "
778 : : "neither 4 (tetrahedra) nor 3 (triangles)." );
779 [ - + ][ - - ]: 32880 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
780 [ + + ][ + - ]: 32880 : Assert( !esup.first.empty(), "Attempt to call genEsued() with empty esup1" );
[ + - ][ + - ]
781 [ + + ][ + - ]: 32879 : Assert( !esup.second.empty(), "Attempt to call genEsued() with empty esup2" );
[ + - ][ + - ]
782 : :
783 : : // find out number of points in mesh connectivity
784 [ + - ]: 32878 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
785 [ - + ][ - - ]: 32878 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
786 : 32878 : auto npoin = *minmax.second + 1;
787 : :
788 : 32878 : auto& esup1 = esup.first;
789 : 32878 : auto& esup2 = esup.second;
790 : :
791 : : // allocate and fill with zeros a temporary array, only used locally
792 [ + - ]: 65756 : std::vector< std::size_t > lpoin( npoin, 0 );
793 : :
794 : : // lambda that returns true if element e contains edge (p < q)
795 : 3113379141 : auto has = [ &inpoel, nnpe ]( std::size_t e, std::size_t p, std::size_t q ) {
796 : 239490753 : int sp = 0;
797 [ + + ]: 1197453549 : for (std::size_t n=0; n<nnpe; ++n)
798 [ + + ][ + + ]: 957962796 : if (inpoel[e*nnpe+n] == p || inpoel[e*nnpe+n] == q) ++sp;
[ + + ]
799 : 239490753 : return sp == 2;
800 : 32878 : };
801 : :
802 : : // map to associate edges to unique surrounding element ids
803 : : std::unordered_map< UnsMesh::Edge, std::vector< std::size_t >,
804 : 32878 : UnsMesh::Hash<2>, UnsMesh::Eq<2> > esued;
805 : :
806 : : // generate edges and associated vector of unique surrounding element ids
807 [ + + ]: 3801408 : for (std::size_t p=0; p<npoin; ++p)
808 [ + + ]: 44277622 : for (std::size_t i=esup2[p]+1; i<=esup2[p+1]; ++i )
809 [ + + ]: 202545388 : for (std::size_t n=0; n<nnpe; ++n) {
810 : 162036296 : auto q = inpoel[ esup1[i] * nnpe + n ];
811 [ + + ][ + + ]: 162036296 : if (q != p && lpoin[q] != p+1) {
[ + + ]
812 [ + + ]: 34539934 : if (p < q) { // for edge given point ids p < q
813 [ + + ]: 256760720 : for (std::size_t j=esup2[p]+1; j<=esup2[p+1]; ++j ) {
814 : 239490753 : auto e = esup1[j];
815 [ + + ][ + - ]: 239490753 : if (has(e,p,q)) esued[{p,q}].push_back(e);
[ + - ]
816 : : }
817 : : }
818 : 34539934 : lpoin[q] = p+1;
819 : : }
820 : : }
821 : :
822 : : // sort element ids surrounding edges for each edge
823 [ + + ][ + - ]: 17302845 : for (auto& p : esued) std::sort( begin(p.second), end(p.second) );
824 : :
825 : : // Return elements surrounding edges data structure
826 : 65756 : return esued;
827 : : }
828 : :
829 : : std::size_t
830 : 1 : genNbfacTet( std::size_t tnbfac,
831 : : const std::vector< std::size_t >& inpoel,
832 : : const std::vector< std::size_t >& triinpoel_complete,
833 : : const std::map< int, std::vector< std::size_t > >& bface_complete,
834 : : const std::unordered_map< std::size_t, std::size_t >& lid,
835 : : std::vector< std::size_t >& triinpoel,
836 : : std::map< int, std::vector< std::size_t > >& bface )
837 : : // *****************************************************************************
838 : : // Generate the number of boundary-faces and the triangle boundary-face
839 : : // connectivity for a chunk of a full mesh.
840 : : // \warning This is for Triangular face-elements only.
841 : : //! \param[in] tnbfac Total number of boundary faces in the entire mesh.
842 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
843 : : //! node ids of each element of an unstructured mesh.
844 : : //! \param[in] triinpoel_complete Interconnectivity of points and boundary-face
845 : : //! in the entire mesh.
846 : : //! \param[in] bface_complete Map of boundary-face lists mapped to corresponding
847 : : //! side set ids for the entire mesh.
848 : : //! \param[in] lid Mapping between the node indices used in the smaller inpoel
849 : : //! connectivity (a subset of the entire triinpoel_complete connectivity),
850 : : //! e.g., after mesh partitioning.
851 : : //! \param[inout] triinpoel Interconnectivity of points and boundary-face in
852 : : //! this mesh-partition.
853 : : //! \param[inout] bface Map of boundary-face lists mapped to corresponding
854 : : //! side set ids for this mesh-partition
855 : : //! \return Number of boundary-faces on this chare/mesh-partition.
856 : : //! \details This function takes a mesh by its domain-element
857 : : //! (tetrahedron-connectivity) in inpoel and a boundary-face (triangle)
858 : : //! connectivity in triinpoel_complete. Based on these two arrays, it
859 : : //! searches for those faces of triinpoel_complete that are also in inpoel
860 : : //! and as a result it generates (1) the number of boundary faces shared with
861 : : //! the mesh in inpoel and (2) the intersection of the triangle element
862 : : //! connectivity whose faces are shared with inpoel. An example use case is
863 : : //! where triinpoel_complete contains the connectivity for the boundary of the
864 : : //! full problem/mesh and inpoel contains the connectivity for only a chunk of
865 : : //! an already partitioned mesh. This function then intersects
866 : : //! triinpoel_complete with inpoel and returns only those faces that share
867 : : //! nodes with inpoel.
868 : : // *****************************************************************************
869 : : {
870 : 1 : std::size_t nbfac(0), nnpf(3);
871 : :
872 [ + - ]: 1 : if (tnbfac > 0)
873 : : {
874 : :
875 [ - + ][ - - ]: 1 : Assert( !inpoel.empty(),
[ - - ][ - - ]
876 : : "Attempt to call genNbfacTet() on empty inpoel container" );
877 [ - + ][ - - ]: 1 : Assert( !triinpoel_complete.empty(),
[ - - ][ - - ]
878 : : "Attempt to call genNbfacTet() on empty triinpoel_complete container" );
879 [ - + ][ - - ]: 1 : Assert( triinpoel_complete.size()/nnpf == tnbfac,
[ - - ][ - - ]
880 : : "Incorrect size of triinpoel in genNbfacTet()" );
881 : :
882 [ + - ]: 2 : auto nptet = inpoel;
883 [ + - ]: 2 : auto nptri = triinpoel_complete;
884 : :
885 [ + - ]: 1 : unique( nptet );
886 [ + - ]: 1 : unique( nptri );
887 : :
888 : 2 : std::unordered_set< std::size_t > snptet;
889 : :
890 : : // getting the reduced inpoel as a set for quick searches
891 [ + - ]: 1 : snptet.insert( begin(nptet), end(nptet));
892 : :
893 : : // vector to store boundary-face-nodes in this chunk
894 : 2 : std::vector< std::size_t > nptri_chunk;
895 : :
896 : : // getting the nodes of the boundary-faces in this chunk
897 [ + + ]: 15 : for (auto i : nptri)
898 [ + - ][ + - ]: 14 : if (snptet.find(i) != end(snptet))
899 [ + - ]: 14 : nptri_chunk.push_back(i);
900 : :
901 : : std::size_t tag, icoun;
902 : :
903 : : // matching nodes in nptri_chunk with nodes in inpoel and
904 : : // triinpoel_complete to get the number of faces in this chunk
905 [ + + ]: 2 : for (const auto& ss : bface_complete)
906 : : {
907 [ + + ]: 25 : for (auto f : ss.second)
908 : : {
909 : 24 : icoun = f*nnpf;
910 : 24 : tag = 0;
911 [ + + ]: 96 : for (std::size_t i=0; i<nnpf; ++i) {
912 [ + + ]: 1080 : for (auto j : nptri_chunk) {
913 : : // cppcheck-suppress useStlAlgorithm
914 [ + + ]: 1008 : if (triinpoel_complete[icoun+i] == j) ++tag;
915 : : }
916 : : }
917 [ + - ]: 24 : if (tag == nnpf)
918 : : // this is a boundary face
919 : : {
920 [ + + ]: 96 : for (std::size_t i=0; i<nnpf; ++i)
921 : : {
922 : 72 : auto ip = triinpoel_complete[icoun+i];
923 : :
924 : : // find local renumbered node-id to store in triinpoel
925 [ + - ][ + - ]: 72 : triinpoel.push_back( cref_find(lid,ip) );
926 : : }
927 : :
928 [ + - ][ + - ]: 24 : bface[ss.first].push_back(nbfac);
929 : 24 : ++nbfac;
930 : : }
931 : : }
932 : : }
933 : :
934 : : }
935 : :
936 : 1 : return nbfac;
937 : : }
938 : :
939 : : std::vector< int >
940 : 5769 : genEsuelTet( const std::vector< std::size_t >& inpoel,
941 : : const std::pair< std::vector< std::size_t >,
942 : : std::vector< std::size_t > >& esup )
943 : : // *****************************************************************************
944 : : // Generate derived data structure, elements surrounding elements
945 : : // as a fixed length data structure as a full vector, including
946 : : // boundary elements as -1.
947 : : // \warning This is for Tetrahedra only.
948 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
949 : : //! node ids of each element of an unstructured mesh. Example:
950 : : //! \code{.cpp}
951 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
952 : : //! 10, 14, 13, 12 };
953 : : //! \endcode
954 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
955 : : //! and { 10, 14, 13, 12 }.
956 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
957 : : //! \return Vector storing elements surrounding elements
958 : : //! \warning It is not okay to call this function with an empty container for
959 : : //! inpoel or esup.first or esup.second; it will throw an exception.
960 : : //! \details The data generated here is stored in a single vector, with length
961 : : //! nfpe * nelem. Note however, that nelem is not explicitly provided, but
962 : : //! calculated from inpoel. For boundary elements, at the boundary face, this
963 : : //! esuelTet stores value -1 indicating that this is outside the domain. The
964 : : //! convention for numbering the local face (triangle) connectivity is very
965 : : //! important, e.g., in generating the inpofa array later. This node ordering
966 : : //! convention is stored in tk::lpofa. Thus function is specific to
967 : : //! tetrahedra, which is reflected in the fact that nnpe and nfpe are being
968 : : //! set here in the function rather than being input arguments. To find out
969 : : //! the number of elements, _nelem_, the size of the mesh connectivity vector,
970 : : //! _inpoel_, can be devided by the number of nodes per elements, _nnpe_:
971 : : //! \code{.cpp}
972 : : //! auto nelem = inpoel.size()/nnpe;
973 : : //! \endcode
974 : : // *****************************************************************************
975 : : {
976 [ - + ][ - - ]: 5769 : Assert( !inpoel.empty(), "Attempt to call genEsuelTet() on empty container" );
[ - - ][ - - ]
977 [ - + ][ - - ]: 5769 : Assert( !esup.first.empty(), "Attempt to call genEsuelTet() with empty esup1" );
[ - - ][ - - ]
978 [ - + ][ - - ]: 5769 : Assert( !esup.second.empty(),
[ - - ][ - - ]
979 : : "Attempt to call genEsuelTet() with empty esup2" );
980 : :
981 : 5769 : auto& esup1 = esup.first;
982 : 5769 : auto& esup2 = esup.second;
983 : :
984 : : // set tetrahedron geometry
985 : 5769 : std::size_t nnpe(4), nfpe(4), nnpf(3);
986 : :
987 [ - + ][ - - ]: 5769 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by four" );
[ - - ][ - - ]
988 : :
989 : : // get nelem and npoin
990 : 5769 : auto nelem = inpoel.size()/nnpe;
991 [ + - ]: 5769 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
992 [ - + ][ - - ]: 5769 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
993 : 5769 : auto npoin = *minmax.second + 1;
994 : :
995 [ + - ]: 5769 : std::vector< int > esuelTet(nfpe*nelem, -1);
996 [ + - ]: 11538 : std::vector< std::size_t > lhelp(nnpf,0),
997 [ + - ]: 11538 : lpoin(npoin,0);
998 : :
999 [ + + ]: 3195178 : for (std::size_t e=0; e<nelem; ++e)
1000 : : {
1001 : 3189409 : auto mark = nnpe*e;
1002 [ + + ]: 15947045 : for (std::size_t fe=0; fe<nfpe; ++fe)
1003 : : {
1004 : : // array which stores points on this face
1005 : 12757636 : lhelp[0] = inpoel[mark+lpofa[fe][0]];
1006 : 12757636 : lhelp[1] = inpoel[mark+lpofa[fe][1]];
1007 : 12757636 : lhelp[2] = inpoel[mark+lpofa[fe][2]];
1008 : :
1009 : : // mark in this array
1010 : 12757636 : lpoin[lhelp[0]] = 1;
1011 : 12757636 : lpoin[lhelp[1]] = 1;
1012 : 12757636 : lpoin[lhelp[2]] = 1;
1013 : :
1014 : : // select a point on this face
1015 : 12757636 : auto ipoin = lhelp[0];
1016 : :
1017 : : // loop over elements around this point
1018 [ + + ]: 249127538 : for (std::size_t j=esup2[ipoin]+1; j<=esup2[ipoin+1]; ++j )
1019 : : {
1020 : 236369902 : auto jelem = esup1[j];
1021 : : // if this jelem is not e itself then proceed
1022 [ + + ]: 236369902 : if (jelem != e)
1023 : : {
1024 [ + + ]: 1118061330 : for (std::size_t fj=0; fj<nfpe; ++fj)
1025 : : {
1026 : 894449064 : std::size_t icoun(0);
1027 [ + + ]: 3577796256 : for (std::size_t jnofa=0; jnofa<nnpf; ++jnofa)
1028 : : {
1029 : 2683347192 : auto markj = jelem*nnpe;
1030 : 2683347192 : auto jpoin = inpoel[markj+lpofa[fj][jnofa]];
1031 [ + + ]: 2683347192 : if (lpoin[jpoin] == 1) { ++icoun; }
1032 : : }
1033 : : //store esuel if
1034 [ + + ]: 894449064 : if (icoun == nnpf)
1035 : : {
1036 : 11414994 : auto markf = nfpe*e;
1037 : 11414994 : esuelTet[markf+fe] = static_cast<int>(jelem);
1038 : :
1039 : 11414994 : markf = nfpe*jelem;
1040 : 11414994 : esuelTet[markf+fj] = static_cast<int>(e);
1041 : : }
1042 : : }
1043 : : }
1044 : : }
1045 : : // reset this array
1046 : 12757636 : lpoin[lhelp[0]] = 0;
1047 : 12757636 : lpoin[lhelp[1]] = 0;
1048 : 12757636 : lpoin[lhelp[2]] = 0;
1049 : : }
1050 : : }
1051 : :
1052 : 11538 : return esuelTet;
1053 : : }
1054 : :
1055 : : std::size_t
1056 : 1103 : genNipfac( std::size_t nfpe,
1057 : : std::size_t nbfac,
1058 : : const std::vector< int >& esuelTet )
1059 : : // *****************************************************************************
1060 : : // Generate derived data, number of internal and physical-boundary faces in the
1061 : : // mesh. This does not include faces that are on partition/chare-boundaries.
1062 : : //! \param[in] nfpe Number of faces per element.
1063 : : //! \param[in] nbfac Number of boundary faces.
1064 : : //! \param[in] esuelTet Elements surrounding elements.
1065 : : //! \return Total number of faces in the mesh
1066 : : //! \details The unsigned integer here gives the number of internal and
1067 : : // physical-boundary faces in the mesh.
1068 : : // *****************************************************************************
1069 : : {
1070 [ - + ][ - - ]: 1103 : Assert( !esuelTet.empty(), "Attempt to call genNipfac() with empty esuelTet" );
[ - - ][ - - ]
1071 [ - + ][ - - ]: 1103 : Assert( esuelTet.size()%nfpe == 0,
[ - - ][ - - ]
1072 : : "Size of esuelTet must be divisible by nfpe" );
1073 [ - + ][ - - ]: 1103 : Assert( nfpe > 0, "Attempt to call genNipfac() with zero faces per element" );
[ - - ][ - - ]
1074 : :
1075 : 1103 : auto nelem = esuelTet.size()/nfpe;
1076 : :
1077 : 1103 : std::size_t nifac = 0;
1078 : :
1079 : : // loop through elements surrounding elements to find number of internal faces
1080 [ + + ]: 274653 : for (std::size_t e=0; e<nelem; ++e)
1081 : : {
1082 [ + + ]: 1367750 : for (std::size_t ip=nfpe*e; ip<nfpe*(e+1); ++ip)
1083 : : {
1084 [ + + ]: 1094200 : if (esuelTet[ip] != -1)
1085 : : {
1086 [ + + ]: 933710 : if ( e<static_cast< std::size_t >(esuelTet[ip]) )
1087 : : {
1088 : 466855 : ++nifac;
1089 : : }
1090 : : }
1091 : : }
1092 : : }
1093 : :
1094 : 1103 : return nifac + nbfac;
1095 : : }
1096 : :
1097 : : std::vector< int >
1098 : 1100 : genEsuf( std::size_t nfpe,
1099 : : std::size_t nipfac,
1100 : : std::size_t nbfac,
1101 : : const std::vector< std::size_t >& belem,
1102 : : const std::vector< int >& esuelTet )
1103 : : // *****************************************************************************
1104 : : // Generate derived data, elements surrounding faces
1105 : : //! \param[in] nfpe Number of faces per element.
1106 : : //! \param[in] nipfac Number of internal and physical-boundary faces.
1107 : : //! \param[in] nbfac Number of boundary faces.
1108 : : //! \param[in] belem Boundary element vector.
1109 : : //! \param[in] esuelTet Elements surrounding elements.
1110 : : //! \return Elements surrounding faces.
1111 : : //! \details The unsigned integer vector gives the IDs of the elements to the
1112 : : // left and the right of each face in the mesh. The convention followed
1113 : : // throughout is : The left element always has an ID smaller than the ID of
1114 : : // the right element.
1115 : : // *****************************************************************************
1116 : : {
1117 [ - + ][ - - ]: 1100 : Assert( esuelTet.size()%nfpe == 0,
[ - - ][ - - ]
1118 : : "Size of esuelTet must be divisible by nfpe" );
1119 [ - + ][ - - ]: 1100 : Assert( nfpe > 0, "Attempt to call genEsuf() with zero faces per element" );
[ - - ][ - - ]
1120 : :
1121 : 1100 : auto nelem = esuelTet.size()/nfpe;
1122 : :
1123 [ + - ]: 1100 : std::vector< int > esuf(2*nipfac);
1124 : :
1125 : : // counters for number of internal and boundary faces
1126 : 1100 : std::size_t icoun(2*nbfac), bcoun(0);
1127 : :
1128 : : // loop to get face-element connectivity for internal faces
1129 [ + + ]: 274529 : for (std::size_t e=0; e<nelem; ++e) {
1130 [ + + ]: 1367145 : for (std::size_t ip=nfpe*e; ip<nfpe*(e+1); ++ip) {
1131 : 1093716 : auto jelem = esuelTet[ip];
1132 [ + + ]: 1093716 : if (jelem != -1)
1133 : : {
1134 [ + + ]: 933322 : if ( e < static_cast< std::size_t >(jelem) )
1135 : : {
1136 : 466661 : esuf[icoun] = static_cast< int >(e);
1137 : 466661 : esuf[icoun+1] = static_cast< int >(jelem);
1138 : 466661 : icoun = icoun + 2;
1139 : : }
1140 : : }
1141 : : }
1142 : : }
1143 : :
1144 : : // loop to get face-element connectivity for physical-boundary faces
1145 : 1100 : bcoun = 0;
1146 [ + + ]: 123838 : for (auto ie : belem) {
1147 : 122738 : esuf[bcoun] = static_cast< int >(ie);
1148 : 122738 : esuf[bcoun+1] = -1; // outside domain
1149 : 122738 : bcoun = bcoun + 2;
1150 : : }
1151 : :
1152 : 1100 : return esuf;
1153 : : }
1154 : :
1155 : : std::vector< std::size_t >
1156 : 1102 : genInpofaTet( std::size_t nipfac,
1157 : : std::size_t nbfac,
1158 : : const std::vector< std::size_t >& inpoel,
1159 : : const std::vector< std::size_t >& triinpoel,
1160 : : const std::vector< int >& esuelTet )
1161 : : // *****************************************************************************
1162 : : // Generate derived data, points on faces for tetrahedra only
1163 : : //! \param[in] nipfac Number of internal and physical-boundary faces.
1164 : : //! \param[in] nbfac Number of boundary faces.
1165 : : //! \param[in] inpoel Element-node connectivity.
1166 : : //! \param[in] triinpoel Face-node connectivity.
1167 : : //! \param[in] esuelTet Elements surrounding elements.
1168 : : //! \return Points surrounding faces. The unsigned integer vector gives the
1169 : : //! elements to the left and to the right of each face in the mesh.
1170 : : // *****************************************************************************
1171 : : {
1172 : 1102 : std::vector< std::size_t > inpofa;
1173 : :
1174 : : // set tetrahedron geometry
1175 : 1102 : std::size_t nnpe(4), nfpe(4), nnpf(3);
1176 : :
1177 [ - + ][ - - ]: 1102 : Assert( esuelTet.size()%nfpe == 0,
[ - - ][ - - ]
1178 : : "Size of esuelTet must be divisible by nfpe" );
1179 [ - + ][ - - ]: 1102 : Assert( inpoel.size()%nnpe == 0,
[ - - ][ - - ]
1180 : : "Size of inpoel must be divisible by nnpe" );
1181 : :
1182 [ + - ]: 1102 : inpofa.resize(nnpf*nipfac);
1183 : :
1184 : : // counters for number of internal and boundary faces
1185 : 1102 : std::size_t icoun(nnpf*nbfac);
1186 : :
1187 : : // loop over elems to get nodes on faces
1188 : : // this fills the interior face-node connectivity part
1189 [ + + ]: 274579 : for (std::size_t e=0; e<inpoel.size()/nnpe; ++e)
1190 : : {
1191 : 273477 : auto mark = nnpe*e;
1192 [ + + ]: 1367385 : for (std::size_t f=0; f<nfpe ; ++f)
1193 : : {
1194 : 1093908 : auto ip = nfpe*e + f;
1195 : 1093908 : auto jelem = esuelTet[ip];
1196 [ + + ]: 1093908 : if (jelem != -1)
1197 : : {
1198 [ + + ]: 933466 : if ( e < static_cast< std::size_t >(jelem) )
1199 : : {
1200 : 466733 : inpofa[icoun] = inpoel[mark+lpofa[f][0]];
1201 : 466733 : inpofa[icoun+1] = inpoel[mark+lpofa[f][1]];
1202 : 466733 : inpofa[icoun+2] = inpoel[mark+lpofa[f][2]];
1203 : 466733 : icoun = icoun + nnpf;
1204 : : }
1205 : : }
1206 : : }
1207 : : }
1208 : :
1209 : : // this fills the boundary face-node connectivity part
1210 : : // consistent with triinpoel
1211 [ + + ]: 123888 : for (std::size_t f=0; f<nbfac; ++f)
1212 : : {
1213 : 122786 : icoun = nnpf * f;
1214 : 122786 : inpofa[icoun+0] = triinpoel[icoun+0];
1215 : 122786 : inpofa[icoun+1] = triinpoel[icoun+1];
1216 : 122786 : inpofa[icoun+2] = triinpoel[icoun+2];
1217 : : }
1218 : :
1219 : 1102 : return inpofa;
1220 : : }
1221 : :
1222 : : std::vector< std::size_t >
1223 : 1101 : genBelemTet( std::size_t nbfac,
1224 : : const std::vector< std::size_t >& inpofa,
1225 : : const std::pair< std::vector< std::size_t >,
1226 : : std::vector< std::size_t > >& esup )
1227 : : // *****************************************************************************
1228 : : // Generate derived data, and array of elements which share one or more of
1229 : : // their faces with the physical boundary, i.e. where exodus specifies a
1230 : : // side-set for faces. Such elements are sometimes also called host or
1231 : : // boundary elements.
1232 : : //! \param[in] nbfac Number of boundary faces.
1233 : : //! \param[in] inpofa Face-node connectivity.
1234 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
1235 : : //! \return Host elements or boundary elements. The unsigned integer vector
1236 : : //! gives the elements to the left of each boundary face in the mesh.
1237 : : // *****************************************************************************
1238 : : {
1239 [ + - ]: 1101 : std::vector< std::size_t > belem(nbfac);
1240 : :
1241 [ + + ]: 1101 : if (nbfac > 0)
1242 : : {
1243 : :
1244 : : // set tetrahedron geometry
1245 : 1059 : std::size_t nnpf(3), tag(0);
1246 : :
1247 : : // loop over all the boundary faces
1248 [ + + ]: 123797 : for(std::size_t f=0; f<nbfac; ++f)
1249 : : {
1250 : 122738 : belem[f] = 0;
1251 : :
1252 : : // array storing the element-cluster around face
1253 : 245476 : std::vector< std::size_t > elemcluster;
1254 : :
1255 : : // loop over the nodes of this boundary face
1256 [ + + ]: 490952 : for(std::size_t lp=0; lp<nnpf; ++lp)
1257 : : {
1258 : 368214 : auto gp = inpofa[nnpf*f + lp];
1259 : :
1260 [ - + ][ - - ]: 368214 : Assert( gp < esup.second.size(), "Indexing out of esup2" );
[ - - ][ - - ]
1261 : : // loop over elements surrounding this node
1262 [ + + ]: 4175754 : for (auto i=esup.second[gp]+1; i<=esup.second[gp+1]; ++i)
1263 : : {
1264 : : // form element-cluster vector
1265 [ + - ]: 3807540 : elemcluster.push_back(esup.first[i]);
1266 : : }
1267 : : }
1268 : :
1269 : : // loop over element cluster to find repeating elements
1270 [ + - ]: 747875 : for(std::size_t i=0; i<elemcluster.size(); ++i)
1271 : : {
1272 : 747875 : auto ge = elemcluster[i];
1273 : 747875 : tag = 1;
1274 [ + + ]: 25552991 : for(std::size_t j=0; j<elemcluster.size(); ++j)
1275 : : {
1276 [ + + ][ + + ]: 24805116 : if ( i != j && elemcluster[j] == ge )
[ + + ]
1277 : : {
1278 : 453019 : tag++;
1279 : : }
1280 : : }
1281 [ + + ]: 747875 : if (tag == nnpf)
1282 : : {
1283 : : // this is the required boundary element
1284 : 122738 : belem[f] = ge;
1285 : 122738 : break;
1286 : : }
1287 : : }
1288 : : }
1289 : : }
1290 : :
1291 : 1101 : return belem;
1292 : : }
1293 : :
1294 : : Fields
1295 : 1100 : genGeoFaceTri( std::size_t nipfac,
1296 : : const std::vector< std::size_t >& inpofa,
1297 : : const UnsMesh::Coords& coord )
1298 : : // *****************************************************************************
1299 : : // Generate derived data, which stores the geometry details both internal and
1300 : : // boundary triangular faces in the mesh.
1301 : : //! \param[in] nipfac Number of internal and physical-boundary faces.
1302 : : //! \param[in] inpofa Face-node connectivity.
1303 : : //! \param[in] coord Co-ordinates of nodes in this mesh-chunk.
1304 : : //! \return Face geometry information. This includes face area, unit normal
1305 : : //! pointing outward of the element to the left of the face, and face
1306 : : //! centroid coordinates. Use the following examples to access this
1307 : : //! information for face-f.
1308 : : //! \details
1309 : : //! face area: geoFace(f,0,0),
1310 : : //! unit-normal x-component: geoFace(f,1,0),
1311 : : //! y-component: geoFace(f,2,0),
1312 : : //! z-component: geoFace(f,3,0),
1313 : : //! centroid x-coordinate: geoFace(f,4,0),
1314 : : //! y-coordinate: geoFace(f,5,0),
1315 : : //! z-coordinate: geoFace(f,6,0).
1316 : : // *****************************************************************************
1317 : : {
1318 : 1100 : Fields geoFace( nipfac, 7 );
1319 : :
1320 : : // set triangle geometry
1321 : 1100 : std::size_t nnpf(3);
1322 : :
1323 [ - + ][ - - ]: 1100 : Assert( inpofa.size()%nnpf == 0,
[ - - ][ - - ]
1324 : : "Size of inpofa must be divisible by nnpf" );
1325 : :
1326 [ + + ]: 590333 : for(std::size_t f=0; f<nipfac; ++f)
1327 : : {
1328 : : std::size_t ip1, ip2, ip3;
1329 : : real xp1, yp1, zp1,
1330 : : xp2, yp2, zp2,
1331 : : xp3, yp3, zp3;
1332 : :
1333 : : // get area
1334 : 589233 : ip1 = inpofa[nnpf*f];
1335 : 589233 : ip2 = inpofa[nnpf*f + 1];
1336 : 589233 : ip3 = inpofa[nnpf*f + 2];
1337 : :
1338 : 589233 : xp1 = coord[0][ip1];
1339 : 589233 : yp1 = coord[1][ip1];
1340 : 589233 : zp1 = coord[2][ip1];
1341 : :
1342 : 589233 : xp2 = coord[0][ip2];
1343 : 589233 : yp2 = coord[1][ip2];
1344 : 589233 : zp2 = coord[2][ip2];
1345 : :
1346 : 589233 : xp3 = coord[0][ip3];
1347 : 589233 : yp3 = coord[1][ip3];
1348 : 589233 : zp3 = coord[2][ip3];
1349 : :
1350 : : auto geoif = geoFaceTri( {{xp1, xp2, xp3}},
1351 : : {{yp1, yp2, yp3}},
1352 [ + - ]: 1178466 : {{zp1, zp2, zp3}} );
1353 : :
1354 [ + + ]: 4713864 : for (std::size_t i=0; i<7; ++i)
1355 [ + - ][ + - ]: 4124631 : geoFace(f,i) = geoif(0,i);
1356 : : }
1357 : :
1358 : 1100 : return geoFace;
1359 : : }
1360 : :
1361 : : std::array< real, 3 >
1362 : 2996317 : normal( const std::array< real, 3 >& x,
1363 : : const std::array< real, 3 >& y,
1364 : : const std::array< real, 3 >& z )
1365 : : // *****************************************************************************
1366 : : //! Compute the unit normal vector of a triangle
1367 : : //! \param[in] x x-coordinates of the three vertices of the triangle
1368 : : //! \param[in] y y-coordinates of the three vertices of the triangle
1369 : : //! \param[in] z z-coordinates of the three vertices of the triangle
1370 : : //! \return Unit normal
1371 : : // *****************************************************************************
1372 : : {
1373 : : real nx, ny, nz;
1374 : 2996317 : normal( x[0],x[1],x[2], y[0],y[1],y[2], z[0],z[1],z[2], nx, ny, nz );
1375 : 2996317 : return { std::move(nx), std::move(ny), std::move(nz) };
1376 : : }
1377 : :
1378 : : real
1379 : 2996317 : area( const std::array< real, 3 >& x,
1380 : : const std::array< real, 3 >& y,
1381 : : const std::array< real, 3 >& z )
1382 : : // *****************************************************************************
1383 : : //! Compute the are of a triangle
1384 : : //! \param[in] x x-coordinates of the three vertices of the triangle
1385 : : //! \param[in] y y-coordinates of the three vertices of the triangle
1386 : : //! \param[in] z z-coordinates of the three vertices of the triangle
1387 : : //! \return Area
1388 : : // *****************************************************************************
1389 : : {
1390 : 2996317 : return area( x[0],x[1],x[2], y[0],y[1],y[2], z[0],z[1],z[2] );
1391 : : }
1392 : :
1393 : : Fields
1394 : 2996317 : geoFaceTri( const std::array< real, 3 >& x,
1395 : : const std::array< real, 3 >& y,
1396 : : const std::array< real, 3 >& z )
1397 : : // *****************************************************************************
1398 : : //! Compute geometry of the face given by three vertices
1399 : : //! \param[in] x x-coordinates of the three vertices of the triangular face.
1400 : : //! \param[in] y y-coordinates of the three vertices of the triangular face.
1401 : : //! \param[in] z z-coordinates of the three vertices of the triangular face.
1402 : : //! \return Face geometry information. This includes face area, unit normal
1403 : : //! pointing outward of the element to the left of the face, and face
1404 : : //! centroid coordinates.
1405 : : //! \details
1406 : : //! face area: geoFace(f,0,0),
1407 : : //! unit-normal x-component: geoFace(f,1,0),
1408 : : //! y-component: geoFace(f,2,0),
1409 : : //! z-component: geoFace(f,3,0),
1410 : : //! centroid x-coordinate: geoFace(f,4,0),
1411 : : //! y-coordinate: geoFace(f,5,0),
1412 : : //! z-coordinate: geoFace(f,6,0).
1413 : : // *****************************************************************************
1414 : : {
1415 [ + - ]: 2996317 : Fields geoiFace( 1, 7 );
1416 : :
1417 : : // compute area
1418 [ + - ][ + - ]: 2996317 : geoiFace(0,0) = area( x, y, z );
1419 : :
1420 : : // compute unit normal to face
1421 [ + - ]: 2996317 : auto n = normal( x, y, z );
1422 [ + - ]: 2996317 : geoiFace(0,1) = n[0];
1423 [ + - ]: 2996317 : geoiFace(0,2) = n[1];
1424 [ + - ]: 2996317 : geoiFace(0,3) = n[2];
1425 : :
1426 : : // compute centroid
1427 [ + - ]: 2996317 : geoiFace(0,4) = (x[0]+x[1]+x[2])/3.0;
1428 [ + - ]: 2996317 : geoiFace(0,5) = (y[0]+y[1]+y[2])/3.0;
1429 [ + - ]: 2996317 : geoiFace(0,6) = (z[0]+z[1]+z[2])/3.0;
1430 : :
1431 : 5992634 : return geoiFace;
1432 : : }
1433 : :
1434 : : Fields
1435 : 3282 : genGeoElemTet( const std::vector< std::size_t >& inpoel,
1436 : : const UnsMesh::Coords& coord )
1437 : : // *****************************************************************************
1438 : : // Generate derived data, which stores the geometry details of tetrahedral
1439 : : // elements.
1440 : : //! \param[in] inpoel Element-node connectivity.
1441 : : //! \param[in] coord Co-ordinates of nodes in this mesh-chunk.
1442 : : //! \return Element geometry information. This includes element volume,
1443 : : //! element centroid coordinates, and minimum edge length. Use the following
1444 : : //! examples to access this information for element-e.
1445 : : //! volume: geoElem(e,0,0),
1446 : : //! centroid x-coordinate: geoElem(e,1,0),
1447 : : //! y-coordinate: geoElem(e,2,0),
1448 : : //! z-coordinate: geoElem(e,3,0).
1449 : : //! minimum edge-length: geoElem(e,4,0).
1450 : : // *****************************************************************************
1451 : : {
1452 : : // set tetrahedron geometry
1453 : 3282 : std::size_t nnpe(4);
1454 : :
1455 [ - + ][ - - ]: 3282 : Assert( inpoel.size()%nnpe == 0,
[ - - ][ - - ]
1456 : : "Size of inpoel must be divisible by nnpe" );
1457 : :
1458 : 3282 : auto nelem = inpoel.size()/nnpe;
1459 : :
1460 : 3282 : Fields geoElem( nelem, 5 );
1461 : :
1462 : 3282 : const auto& x = coord[0];
1463 : 3282 : const auto& y = coord[1];
1464 : 3282 : const auto& z = coord[2];
1465 : :
1466 [ + + ]: 1874069 : for(std::size_t e=0; e<nelem; ++e)
1467 : : {
1468 : : // get volume
1469 : 1870787 : const auto A = inpoel[nnpe*e+0];
1470 : 1870787 : const auto B = inpoel[nnpe*e+1];
1471 : 1870787 : const auto C = inpoel[nnpe*e+2];
1472 : 1870787 : const auto D = inpoel[nnpe*e+3];
1473 : 1870787 : std::array< real, 3 > ba{{ x[B]-x[A], y[B]-y[A], z[B]-z[A] }},
1474 : 1870787 : ca{{ x[C]-x[A], y[C]-y[A], z[C]-z[A] }},
1475 : 1870787 : da{{ x[D]-x[A], y[D]-y[A], z[D]-z[A] }};
1476 : :
1477 : 1870787 : const auto vole = triple( ba, ca, da ) / 6.0;
1478 : :
1479 [ - + ][ - - ]: 1870787 : Assert( vole > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
1480 : :
1481 [ + - ]: 1870787 : geoElem(e,0) = vole;
1482 : :
1483 : : // get centroid
1484 [ + - ]: 1870787 : geoElem(e,1) = (x[A]+x[B]+x[C]+x[D])/4.0;
1485 [ + - ]: 1870787 : geoElem(e,2) = (y[A]+y[B]+y[C]+y[D])/4.0;
1486 [ + - ]: 1870787 : geoElem(e,3) = (z[A]+z[B]+z[C]+z[D])/4.0;
1487 : :
1488 : : // calculate minimum edge-length
1489 : 1870787 : tk::real edgelen = std::numeric_limits< tk::real >::max();
1490 [ + + ]: 7483148 : for (std::size_t i=0; i<nnpe-1; ++i)
1491 : : {
1492 [ + + ]: 16837083 : for (std::size_t j=i+1; j<nnpe; ++j)
1493 : : {
1494 : 11224722 : auto ni(inpoel[nnpe*e+i]), nj(inpoel[nnpe*e+j]);
1495 : 22449444 : edgelen = std::min( edgelen, tk::length( x[ni]-x[nj], y[ni]-y[nj],
1496 : 22449444 : z[ni]-z[nj] ) );
1497 : : }
1498 : : }
1499 [ + - ]: 1870787 : geoElem(e,4) = edgelen;
1500 : : }
1501 : :
1502 : 3282 : return geoElem;
1503 : : }
1504 : :
1505 : : bool
1506 : 4495 : leakyPartition( const std::vector< int >& esueltet,
1507 : : const std::vector< std::size_t >& inpoel,
1508 : : const UnsMesh::Coords& coord )
1509 : : // *****************************************************************************
1510 : : // Perform leak-test on mesh (partition)
1511 : : //! \param[in] esueltet Elements surrounding elements for tetrahedra, see
1512 : : //! tk::genEsueltet()
1513 : : //! \param[in] inpoel Element connectivity
1514 : : //! \param[in] coord Node coordinates
1515 : : //! \details This function computes a surface integral over the boundary of the
1516 : : //! incoming mesh (partition). A non-zero vector result indicates a leak, e.g.,
1517 : : //! a hole in the mesh (partition), which indicates an error either in the
1518 : : // mesh geometry, mesh partitioning, or in the data structures that represent
1519 : : // faces.
1520 : : //! \return True if partition leaks.
1521 : : // *****************************************************************************
1522 : : {
1523 : 4495 : const auto& x = coord[0];
1524 : 4495 : const auto& y = coord[1];
1525 : 4495 : const auto& z = coord[2];
1526 : :
1527 : : // Storage for surface integral over our mesh partition
1528 : 4495 : std::array< real, 3 > s{{ 0.0, 0.0, 0.0}};
1529 : :
1530 [ + + ]: 1684244 : for (std::size_t e=0; e<esueltet.size()/4; ++e) { // for all our tets
1531 : 1679749 : auto mark = e*4;
1532 [ + + ]: 8398745 : for (std::size_t f=0; f<4; ++f) // for all tet faces
1533 [ + + ]: 6718996 : if (esueltet[mark+f] == -1) { // if face has no outside-neighbor tet
1534 : : // 3 local node IDs of face
1535 : 849798 : auto A = inpoel[ mark + lpofa[f][0] ];
1536 : 849798 : auto B = inpoel[ mark + lpofa[f][1] ];
1537 : 849798 : auto C = inpoel[ mark + lpofa[f][2] ];
1538 : : // Compute geometry data for face
1539 : 2549394 : auto geoface = geoFaceTri( {{x[A], x[B], x[C]}},
1540 : 2549394 : {{y[A], y[B], y[C]}},
1541 [ + - ]: 849798 : {{z[A], z[B], z[C]}} );
1542 : : // Sum up face area * face unit-normal
1543 [ + - ][ + - ]: 849798 : s[0] += geoface(0,0) * geoface(0,1);
1544 [ + - ][ + - ]: 849798 : s[1] += geoface(0,0) * geoface(0,2);
1545 [ + - ][ + - ]: 849798 : s[2] += geoface(0,0) * geoface(0,3);
1546 : : }
1547 : : }
1548 : :
1549 : 4495 : auto eps = 1.0e-9;
1550 [ + - ][ + - ]: 4495 : return std::abs(s[0]) > eps || std::abs(s[1]) > eps || std::abs(s[2]) > eps;
[ - + ]
1551 : : }
1552 : :
1553 : : bool
1554 : 4072 : conforming( const std::vector< std::size_t >& inpoel,
1555 : : const UnsMesh::Coords& coord,
1556 : : bool cerr,
1557 : : const std::vector< std::size_t >& rid )
1558 : : // *****************************************************************************
1559 : : // Check if mesh (partition) is conforming
1560 : : //! \param[in] inpoel Element connectivity
1561 : : //! \param[in] coord Node coordinates
1562 : : //! \param[in] cerr True if hanging-node edge data should be output to
1563 : : //! std::cerr (true by default)
1564 : : //! \param[in] rid AMR Lib node id map
1565 : : //! std::cerr (true by default)
1566 : : //! \return True if mesh (partition) has no hanging nodes and thus the mesh is
1567 : : //! conforming, false if non-conforming.
1568 : : //! \details A conforming mesh by definition has no hanging nodes. A node is
1569 : : //! hanging if an edge of one element coincides with two (or more) edges (of
1570 : : //! two or more other elements). Thus, testing for conformity relies on
1571 : : //! checking the coordinates of all vertices: if any vertex coincides with
1572 : : //! that of a mid-point node of an edge, that is a hanging node. Note that
1573 : : //! this assumes that hanging nodes can only be at the mid-point of edges.
1574 : : //! This may happen after a mesh refinement step, due to a problem/bug,
1575 : : //! within the mesh refinement algorithm given by J. Waltz, Parallel adaptive
1576 : : //! refinement for unsteady flow calculations on 3D unstructured grids,
1577 : : //! International Journal for Numerical Methods in Fluids, 46: 37–57, 2004,
1578 : : //! which always adds/removes vertices at the mid-points of edges of a
1579 : : //! tetrahedron mesh within a single refinement step. Thus this algorithm is
1580 : : //! intended for this specific case, i.e., test for conformity after a
1581 : : //! single refinement step and not after multiple ones or for detecting
1582 : : //! hanging nodes in an arbitrary mesh.
1583 : : //*****************************************************************************
1584 : : {
1585 [ + + ][ + - ]: 4072 : Assert( !inpoel.empty(),
[ + - ][ + - ]
1586 : : "Attempt to call conforming() with empty mesh connectivity" );
1587 [ + + ][ + - ]: 4070 : Assert( inpoel.size() % 4 == 0,
[ + - ][ + - ]
1588 : : "Size of inpoel must be divisible by nnpe" );
1589 [ + - ][ + + ]: 4069 : Assert( *std::min_element( begin(inpoel), end(inpoel) ) == 0,
[ + - ][ + - ]
[ + - ]
1590 : : "Node ids should start from zero" );
1591 [ + - ][ + - ]: 4068 : Assert( !coord[0].empty() && !coord[1].empty() && !coord[2].empty(),
[ + - ][ - - ]
[ - - ][ - - ]
1592 : : "Attempt to call conforming() with empty coordinates container" );
1593 : :
1594 : : using Coord = UnsMesh::Coord;
1595 : : using Edge = UnsMesh::Edge;
1596 : : using Tet = UnsMesh::Tet;
1597 : :
1598 : : // Compare operator to be used as less-than for std::array< tk::real, 3 >,
1599 : : // implemented as a lexicographic ordering.
1600 : : struct CoordLess {
1601 : : const real eps = std::numeric_limits< real >::epsilon();
1602 : 116919762 : bool operator() ( const Coord& lhs, const Coord& rhs ) const {
1603 [ + + ]: 116919762 : if (lhs[0] < rhs[0])
1604 : 50982417 : return true;
1605 [ + + ][ + + ]: 65937345 : else if (std::abs(lhs[0]-rhs[0]) < eps && lhs[1] < rhs[1])
[ + + ]
1606 : 5587267 : return true;
1607 : 60350078 : else if (std::abs(lhs[0]-rhs[0]) < eps &&
1608 [ + + ][ + + ]: 73754795 : std::abs(lhs[1]-rhs[1]) < eps &&
[ + + ]
1609 [ + + ]: 13404717 : lhs[2] < rhs[2])
1610 : 787183 : return true;
1611 : : else
1612 : 59562895 : return false;
1613 : : }
1614 : : };
1615 : :
1616 : : // Map associating data on potential hanging nodes. Key: coordinates of nodes
1617 : : // of edge-half points, value: tet id (local if in parallel), tet connectivity
1618 : : // (using local ids if in parallel), edge ids (local if in parallel).
1619 : : std::map< Coord, // edge-half node coordinates: x, y, z
1620 : : std::tuple< std::size_t, // element id of edge-half node
1621 : : Tet, // element node ids of edge-half node
1622 : : Edge >, // edge containing half-node
1623 : 8136 : CoordLess > edgeNodes;
1624 : :
1625 : 4068 : const auto& x = coord[0];
1626 : 4068 : const auto& y = coord[1];
1627 : 4068 : const auto& z = coord[2];
1628 : :
1629 : : fenv_t fe;
1630 : 4068 : feholdexcept( &fe );
1631 : :
1632 : : // Compute coordinates of nodes of mid-points of all edges
1633 [ + + ]: 1388279 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1634 : 1384211 : auto A = inpoel[e*4+0];
1635 : 1384211 : auto B = inpoel[e*4+1];
1636 : 1384211 : auto C = inpoel[e*4+2];
1637 : 1384211 : auto D = inpoel[e*4+3];
1638 : : std::array<Edge,6> edge{{ {{A,B}}, {{B,C}}, {{A,C}},
1639 : 1384211 : {{A,D}}, {{B,D}}, {{C,D}} }};
1640 [ + + ]: 9689477 : for (const auto& n : edge) {
1641 : 8305266 : Coord en{{ (x[n[0]] + x[n[1]]) / 2.0,
1642 : 8305266 : (y[n[0]] + y[n[1]]) / 2.0,
1643 : 16610532 : (z[n[0]] + z[n[1]]) / 2.0 }};
1644 [ + - ][ + - ]: 8305266 : edgeNodes[ en ] = std::tuple<std::size_t,Tet,Edge>{ e, {{A,B,C,D}}, n };
1645 : : }
1646 : : }
1647 : :
1648 : 4068 : feclearexcept( FE_UNDERFLOW );
1649 : 4068 : feupdateenv( &fe );
1650 : :
1651 : : // Find hanging nodes. If the coordinates of an element vertex coincide with
1652 : : // that of a mid-point node of an edge, that is a hanging node. If we find one
1653 : : // such node we print out some info on it.
1654 : 4068 : auto ix = x.cbegin();
1655 : 4068 : auto iy = y.cbegin();
1656 : 4068 : auto iz = z.cbegin();
1657 : :
1658 : 4068 : bool hanging_node = false;
1659 : :
1660 [ + + ]: 455214 : while (ix != x.cend()) {
1661 : 451146 : Coord n{{ *ix, *iy, *iz }};
1662 [ + - ]: 451146 : auto i = edgeNodes.find( n );
1663 [ + + ]: 451146 : if (i != end(edgeNodes)) {
1664 : 1 : const auto& hanging_node_coord = i->first;
1665 : 1 : const auto& hanging_node_info = i->second;
1666 : 1 : auto tet_id = std::get< 0 >( hanging_node_info );
1667 : 1 : const auto& tet = std::get< 1 >( hanging_node_info );
1668 : 1 : const auto& edge = std::get< 2 >( hanging_node_info );
1669 [ - + ]: 1 : if (cerr) {
1670 : : std::cerr
1671 [ - - ]: 0 : << "Mesh conformity test found hanging node with coordinates"" ("
1672 [ - - ][ - - ]: 0 : << hanging_node_coord[0] << ", "
1673 [ - - ][ - - ]: 0 : << hanging_node_coord[1] << ", "
1674 [ - - ][ - - ]: 0 : << hanging_node_coord[2] << ") of tetrahedron element "
1675 [ - - ][ - - ]: 0 : << tet_id << " with connectivity (" << tet[0] << ','
[ - - ][ - - ]
1676 [ - - ][ - - ]: 0 : << tet[1] << ',' << tet[2] << ',' << tet[3] << ") on edge ("
[ - - ][ - - ]
[ - - ][ - - ]
1677 [ - - ][ - - ]: 0 : << edge[0] << ',' << edge[1] << ")"
[ - - ]
1678 [ - - ][ - - ]: 0 : << "AMR lib node ids for this edge: " << rid[edge[0]] << ','
[ - - ][ - - ]
1679 [ - - ][ - - ]: 0 : << rid[edge[1]] << std::endl;
1680 : : }
1681 : 1 : hanging_node = true;
1682 : : }
1683 : 451146 : ++ix; ++iy; ++iz;
1684 : : }
1685 : :
1686 [ + + ]: 4068 : if (hanging_node) return false;
1687 : :
1688 : 4067 : return true;
1689 : : }
1690 : :
1691 : : bool
1692 : 3486022 : intet( const std::array< std::vector< real >, 3 >& coord,
1693 : : const std::vector< std::size_t >& inpoel,
1694 : : const std::vector< real >& p,
1695 : : std::size_t e,
1696 : : std::array< real, 4 >& N )
1697 : : // *****************************************************************************
1698 : : // Determine if a point is in a tetrahedron
1699 : : //! \param[in] coord Mesh node coordinates
1700 : : //! \param[in] inpoel Mesh element connectivity
1701 : : //! \param[in] p Point coordinates
1702 : : //! \param[in] e Mesh cell index
1703 : : //! \param[in,out] N Shapefunctions evaluated at the point
1704 : : //! \return True if ppoint is in mesh cell
1705 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
1706 : : // *****************************************************************************
1707 : : {
1708 [ - + ][ - - ]: 3486022 : Assert( p.size() == 3, "Size mismatch" );
[ - - ][ - - ]
1709 : :
1710 : : // Tetrahedron node indices
1711 : 3486022 : const auto A = inpoel[e*4+0];
1712 : 3486022 : const auto B = inpoel[e*4+1];
1713 : 3486022 : const auto C = inpoel[e*4+2];
1714 : 3486022 : const auto D = inpoel[e*4+3];
1715 : :
1716 : : // Tetrahedron node coordinates
1717 : 3486022 : const auto& x = coord[0];
1718 : 3486022 : const auto& y = coord[1];
1719 : 3486022 : const auto& z = coord[2];
1720 : :
1721 : : // Point coordinates
1722 : 3486022 : const auto& xp = p[0];
1723 : 3486022 : const auto& yp = p[1];
1724 : 3486022 : const auto& zp = p[2];
1725 : :
1726 : : // Evaluate linear shapefunctions at point locations using Cramer's Rule
1727 : : // | xp | | x1 x2 x3 x4 | | N1 |
1728 : : // | yp | = | y1 y2 y3 y4 | • | N2 |
1729 : : // | zp | | z1 z2 z3 z4 | | N3 |
1730 : : // | 1 | | 1 1 1 1 | | N4 |
1731 : :
1732 : 3486022 : real DetX = (y[B]*z[C] - y[C]*z[B] - y[B]*z[D] + y[D]*z[B] +
1733 : 3486022 : y[C]*z[D] - y[D]*z[C])*x[A] + x[B]*y[C]*z[A] - x[B]*y[A]*z[C] +
1734 : 3486022 : x[C]*y[A]*z[B] - x[C]*y[B]*z[A] + x[B]*y[A]*z[D] - x[B]*y[D]*z[A] -
1735 : 3486022 : x[D]*y[A]*z[B] + x[D]*y[B]*z[A] - x[C]*y[A]*z[D] + x[C]*y[D]*z[A] +
1736 : 3486022 : x[D]*y[A]*z[C] - x[D]*y[C]*z[A] - x[B]*y[C]*z[D] + x[B]*y[D]*z[C] +
1737 : 3486022 : x[C]*y[B]*z[D] - x[C]*y[D]*z[B] - x[D]*y[B]*z[C] + x[D]*y[C]*z[B];
1738 : :
1739 : 3486022 : real DetX1 = (y[D]*z[C] - y[C]*z[D] + y[C]*zp - yp*z[C] -
1740 : 3486022 : y[D]*zp + yp*z[D])*x[B] + x[C]*y[B]*z[D] - x[C]*y[D]*z[B] -
1741 : 3486022 : x[D]*y[B]*z[C] + x[D]*y[C]*z[B] - x[C]*y[B]*zp + x[C]*yp*z[B] +
1742 : 3486022 : xp*y[B]*z[C] - xp*y[C]*z[B] + x[D]*y[B]*zp - x[D]*yp*z[B] -
1743 : 3486022 : xp*y[B]*z[D] + xp*y[D]*z[B] + x[C]*y[D]*zp - x[C]*yp*z[D] -
1744 : 3486022 : x[D]*y[C]*zp + x[D]*yp*z[C] + xp*y[C]*z[D] - xp*y[D]*z[C];
1745 : :
1746 : 3486022 : real DetX2 = (y[C]*z[D] - y[D]*z[C] - y[C]*zp + yp*z[C] +
1747 : 3486022 : y[D]*zp - yp*z[D])*x[A] + x[C]*y[D]*z[A] - x[C]*y[A]*z[D] +
1748 : 3486022 : x[D]*y[A]*z[C] - x[D]*y[C]*z[A] + x[C]*y[A]*zp - x[C]*yp*z[A] -
1749 : 3486022 : xp*y[A]*z[C] + xp*y[C]*z[A] - x[D]*y[A]*zp + x[D]*yp*z[A] +
1750 : 3486022 : xp*y[A]*z[D] - xp*y[D]*z[A] - x[C]*y[D]*zp + x[C]*yp*z[D] +
1751 : 3486022 : x[D]*y[C]*zp - x[D]*yp*z[C] - xp*y[C]*z[D] + xp*y[D]*z[C];
1752 : :
1753 : 3486022 : real DetX3 = (y[D]*z[B] - y[B]*z[D] + y[B]*zp - yp*z[B] -
1754 : 3486022 : y[D]*zp + yp*z[D])*x[A] + x[B]*y[A]*z[D] - x[B]*y[D]*z[A] -
1755 : 3486022 : x[D]*y[A]*z[B] + x[D]*y[B]*z[A] - x[B]*y[A]*zp + x[B]*yp*z[A] +
1756 : 3486022 : xp*y[A]*z[B] - xp*y[B]*z[A] + x[D]*y[A]*zp - x[D]*yp*z[A] -
1757 : 3486022 : xp*y[A]*z[D] + xp*y[D]*z[A] + x[B]*y[D]*zp - x[B]*yp*z[D] -
1758 : 3486022 : x[D]*y[B]*zp + x[D]*yp*z[B] + xp*y[B]*z[D] - xp*y[D]*z[B];
1759 : :
1760 : 3486022 : real DetX4 = (y[B]*z[C] - y[C]*z[B] - y[B]*zp + yp*z[B] +
1761 : 3486022 : y[C]*zp - yp*z[C])*x[A] + x[B]*y[C]*z[A] - x[B]*y[A]*z[C] +
1762 : 3486022 : x[C]*y[A]*z[B] - x[C]*y[B]*z[A] + x[B]*y[A]*zp - x[B]*yp*z[A] -
1763 : 3486022 : xp*y[A]*z[B] + xp*y[B]*z[A] - x[C]*y[A]*zp + x[C]*yp*z[A] +
1764 : 3486022 : xp*y[A]*z[C] - xp*y[C]*z[A] - x[B]*y[C]*zp + x[B]*yp*z[C] +
1765 : 3486022 : x[C]*y[B]*zp - x[C]*yp*z[B] - xp*y[B]*z[C] + xp*y[C]*z[B];
1766 : :
1767 : : // Shape functions evaluated at point
1768 : 3486022 : N[0] = DetX1/DetX;
1769 : 3486022 : N[1] = DetX2/DetX;
1770 : 3486022 : N[2] = DetX3/DetX;
1771 : 3486022 : N[3] = DetX4/DetX;
1772 : :
1773 : : // if min( N^i, 1-N^i ) > 0 for all i, point is in cell
1774 [ + + ]: 6107165 : if ( std::min(N[0],1.0-N[0]) > 0 && std::min(N[1],1.0-N[1]) > 0 &&
1775 [ + + ][ + + ]: 6107165 : std::min(N[2],1.0-N[2]) > 0 && std::min(N[3],1.0-N[3]) > 0 )
[ + + ][ + + ]
1776 : : {
1777 : 460107 : return true;
1778 : : } else {
1779 : 3025915 : return false;
1780 : : }
1781 : : }
1782 : :
1783 : : tk::UnsMesh::Coords
1784 : 1113 : curl( const std::array< std::vector< tk::real >, 3 >& coord,
1785 : : const std::vector< std::size_t >& inpoel,
1786 : : const tk::UnsMesh::Coords& v )
1787 : : // *****************************************************************************
1788 : : // Compute curl of a vector field at nodes of unstructured tetrahedra mesh
1789 : : //! \param[in] coord Mesh node coordinates
1790 : : //! \param[in] inpoel Mesh element connectivity
1791 : : //! \param[in] v Vector field whose curl to compute
1792 : : //! \return Weak (partial) result of curl of v (partial beacuse it still needs
1793 : : //! a division by the nodal volumes.
1794 : : // *****************************************************************************
1795 : : {
1796 : : // access node cooordinates
1797 : 1113 : const auto& x = coord[0];
1798 : 1113 : const auto& y = coord[1];
1799 : 1113 : const auto& z = coord[2];
1800 : : // access vector field components
1801 : 1113 : const auto& vx = v[0];
1802 : 1113 : const auto& vy = v[1];
1803 : 1113 : const auto& vz = v[2];
1804 : :
1805 : 1113 : auto npoin = x.size();
1806 : 1113 : tk::UnsMesh::Coords curl;
1807 [ + - ]: 1113 : curl[0].resize( npoin, 0.0 );
1808 [ + - ]: 1113 : curl[1].resize( npoin, 0.0 );
1809 [ + - ]: 1113 : curl[2].resize( npoin, 0.0 );
1810 : :
1811 [ + + ]: 5107071 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1812 : : // access node IDs
1813 : : std::size_t N[4] =
1814 : 5105958 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
1815 : : // compute element Jacobi determinant, J = 6V
1816 : 5105958 : tk::real bax = x[N[1]]-x[N[0]];
1817 : 5105958 : tk::real bay = y[N[1]]-y[N[0]];
1818 : 5105958 : tk::real baz = z[N[1]]-z[N[0]];
1819 : 5105958 : tk::real cax = x[N[2]]-x[N[0]];
1820 : 5105958 : tk::real cay = y[N[2]]-y[N[0]];
1821 : 5105958 : tk::real caz = z[N[2]]-z[N[0]];
1822 : 5105958 : tk::real dax = x[N[3]]-x[N[0]];
1823 : 5105958 : tk::real day = y[N[3]]-y[N[0]];
1824 : 5105958 : tk::real daz = z[N[3]]-z[N[0]];
1825 : 5105958 : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
1826 : 5105958 : auto J24 = J/24.0;
1827 : : // shape function derivatives, nnode*ndim [4][3]
1828 : : tk::real g[4][3];
1829 : 5105958 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
1830 : : g[1][0], g[1][1], g[1][2] );
1831 : 5105958 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
1832 : : g[2][0], g[2][1], g[2][2] );
1833 : 5105958 : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
1834 : : g[3][0], g[3][1], g[3][2] );
1835 [ + + ]: 20423832 : for (std::size_t i=0; i<3; ++i)
1836 : 15317874 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
1837 [ + + ]: 25529790 : for (std::size_t b=0; b<4; ++b) {
1838 : : tk::real c[3];
1839 : 20423832 : tk::cross( g[b][0], g[b][1], g[b][2],
1840 : 20423832 : vx[N[b]], vy[N[b]], vz[N[b]],
1841 : : c[0], c[1], c[2] );
1842 [ + + ]: 102119160 : for (std::size_t a=0; a<4; ++a) {
1843 : 81695328 : curl[0][N[a]] += J24 * c[0];
1844 : 81695328 : curl[1][N[a]] += J24 * c[1];
1845 : 81695328 : curl[2][N[a]] += J24 * c[2];
1846 : : }
1847 : : }
1848 : : }
1849 : :
1850 : 1113 : return curl;
1851 : : }
1852 : :
1853 : : std::vector< tk::real >
1854 : 1111 : div( const std::array< std::vector< tk::real >, 3 >& coord,
1855 : : const std::vector< std::size_t >& inpoel,
1856 : : const tk::UnsMesh::Coords& v )
1857 : : // *****************************************************************************
1858 : : // Compute divergence of vector field at nodes of unstructured tetrahedra mesh
1859 : : //! \param[in] coord Mesh node coordinates
1860 : : //! \param[in] inpoel Mesh element connectivity
1861 : : //! \param[in] v Vector field whose divergence to compute
1862 : : //! \return Weak (partial) result of div v (partial beacuse it still needs
1863 : : //! a division by the nodal volumes.
1864 : : // *****************************************************************************
1865 : : {
1866 : : // access node cooordinates
1867 : 1111 : const auto& x = coord[0];
1868 : 1111 : const auto& y = coord[1];
1869 : 1111 : const auto& z = coord[2];
1870 : :
1871 : 1111 : auto npoin = x.size();
1872 [ + - ]: 1111 : std::vector< tk::real > div( npoin, 0.0 );
1873 : :
1874 [ + + ]: 5107021 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1875 : : // access node IDs
1876 : : std::size_t N[4] =
1877 : 5105910 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
1878 : : // compute element Jacobi determinant, J = 6V
1879 : 5105910 : tk::real bax = x[N[1]]-x[N[0]];
1880 : 5105910 : tk::real bay = y[N[1]]-y[N[0]];
1881 : 5105910 : tk::real baz = z[N[1]]-z[N[0]];
1882 : 5105910 : tk::real cax = x[N[2]]-x[N[0]];
1883 : 5105910 : tk::real cay = y[N[2]]-y[N[0]];
1884 : 5105910 : tk::real caz = z[N[2]]-z[N[0]];
1885 : 5105910 : tk::real dax = x[N[3]]-x[N[0]];
1886 : 5105910 : tk::real day = y[N[3]]-y[N[0]];
1887 : 5105910 : tk::real daz = z[N[3]]-z[N[0]];
1888 : 5105910 : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
1889 : 5105910 : auto J24 = J/24.0;
1890 : : // shape function derivatives, nnode*ndim [4][3]
1891 : : tk::real g[4][3];
1892 : 5105910 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
1893 : : g[1][0], g[1][1], g[1][2] );
1894 : 5105910 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
1895 : : g[2][0], g[2][1], g[2][2] );
1896 : 5105910 : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
1897 : : g[3][0], g[3][1], g[3][2] );
1898 [ + + ]: 20423640 : for (std::size_t i=0; i<3; ++i)
1899 : 15317730 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
1900 [ + + ]: 25529550 : for (std::size_t a=0; a<4; ++a)
1901 [ + + ]: 102118200 : for (std::size_t b=0; b<4; ++b)
1902 [ + + ]: 326778240 : for (std::size_t i=0; i<3; ++i)
1903 : 245083680 : div[N[a]] += J24 * g[b][i] * v[i][N[b]];
1904 : : }
1905 : :
1906 : 1111 : return div;
1907 : : }
1908 : :
1909 : : tk::UnsMesh::Coords
1910 : 124 : grad( const std::array< std::vector< tk::real >, 3 >& coord,
1911 : : const std::vector< std::size_t >& inpoel,
1912 : : const std::vector< tk::real >& phi )
1913 : : // *****************************************************************************
1914 : : // Compute gradient of a scalar field at nodes of unstructured tetrahedra mesh
1915 : : //! \param[in] coord Mesh node coordinates
1916 : : //! \param[in] inpoel Mesh element connectivity
1917 : : //! \param[in] phi Scalar field whose gradient to compute
1918 : : //! \return Weak (partial) result of grad phi (partial beacuse it still needs
1919 : : //! a division by the nodal volumes.
1920 : : // *****************************************************************************
1921 : : {
1922 : : // access node cooordinates
1923 : 124 : const auto& x = coord[0];
1924 : 124 : const auto& y = coord[1];
1925 : 124 : const auto& z = coord[2];
1926 : :
1927 : 124 : auto npoin = x.size();
1928 : 124 : tk::UnsMesh::Coords grad;
1929 [ + - ]: 124 : grad[0].resize( npoin, 0.0 );
1930 [ + - ]: 124 : grad[1].resize( npoin, 0.0 );
1931 [ + - ]: 124 : grad[2].resize( npoin, 0.0 );
1932 : :
1933 [ + + ]: 22754 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1934 : : // access node IDs
1935 : : std::size_t N[4] =
1936 : 22630 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
1937 : : // compute element Jacobi determinant, J = 6V
1938 : 22630 : tk::real bax = x[N[1]]-x[N[0]];
1939 : 22630 : tk::real bay = y[N[1]]-y[N[0]];
1940 : 22630 : tk::real baz = z[N[1]]-z[N[0]];
1941 : 22630 : tk::real cax = x[N[2]]-x[N[0]];
1942 : 22630 : tk::real cay = y[N[2]]-y[N[0]];
1943 : 22630 : tk::real caz = z[N[2]]-z[N[0]];
1944 : 22630 : tk::real dax = x[N[3]]-x[N[0]];
1945 : 22630 : tk::real day = y[N[3]]-y[N[0]];
1946 : 22630 : tk::real daz = z[N[3]]-z[N[0]];
1947 : 22630 : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
1948 : 22630 : auto J24 = J/24.0;
1949 : : // shape function derivatives, nnode*ndim [4][3]
1950 : : tk::real g[4][3];
1951 : 22630 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
1952 : : g[1][0], g[1][1], g[1][2] );
1953 : 22630 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
1954 : : g[2][0], g[2][1], g[2][2] );
1955 : 22630 : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
1956 : : g[3][0], g[3][1], g[3][2] );
1957 [ + + ]: 90520 : for (std::size_t i=0; i<3; ++i)
1958 : 67890 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
1959 [ + + ]: 113150 : for (std::size_t a=0; a<4; ++a)
1960 [ + + ]: 452600 : for (std::size_t b=0; b<4; ++b)
1961 [ + + ]: 1448320 : for (std::size_t i=0; i<3; ++i)
1962 : 1086240 : grad[i][N[a]] += J24 * g[b][i] * phi[N[b]];
1963 : : }
1964 : :
1965 : 124 : return grad;
1966 : : }
1967 : :
1968 : : } // tk::
|
