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 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Generate data structures derived from unstructured mesh
10 : : \details Generate data structures derived from the connectivity information
11 : : of an unstructured mesh.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <set>
16 : : #include <map>
17 : : #include <iterator>
18 : : #include <numeric>
19 : : #include <algorithm>
20 : : #include <type_traits>
21 : : #include <cstddef>
22 : : #include <array>
23 : : #include <unordered_set>
24 : : #include <unordered_map>
25 : : #include <iostream>
26 : : #include <cfenv>
27 : :
28 : : #include "Exception.hpp"
29 : : #include "DerivedData.hpp"
30 : : #include "ContainerUtil.hpp"
31 : : #include "Vector.hpp"
32 : :
33 : : namespace tk {
34 : :
35 : : std::size_t
36 : 2582 : npoin_in_graph( const std::vector< std::size_t >& inpoel )
37 : : // *****************************************************************************
38 : : // Compute number of points (nodes) in mesh from connectivity
39 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
40 : : //! \return Number of mesh points (nodes)
41 : : // *****************************************************************************
42 : : {
43 [ + - ]: 2582 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
44 [ - + ][ - - ]: 2582 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
45 : 2582 : return *minmax.second + 1;
46 : : }
47 : :
48 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
49 : 10283 : genEsup( const std::vector< std::size_t >& inpoel, std::size_t nnpe )
50 : : // *****************************************************************************
51 : : // Generate derived data structure, elements surrounding points
52 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
53 : : //! node ids of each element of an unstructured mesh. Example:
54 : : //! \code{.cpp}
55 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
56 : : //! 10, 14, 13, 12 };
57 : : //! \endcode
58 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
59 : : //! and { 10, 14, 13, 12 }.
60 : : //! \param[in] nnpe Number of nodes per element
61 : : //! \return Linked lists storing elements surrounding points
62 : : //! \warning It is not okay to call this function with an empty container or a
63 : : //! non-positive number of nodes per element; it will throw an exception.
64 : : //! \details The data generated here is stored in a linked list, more precisely,
65 : : //! two linked arrays (vectors), _esup1_ and _esup2_, where _esup2_ holds the
66 : : //! indices at which _esup1_ holds the element ids surrounding points. Looping
67 : : //! over all elements surrounding all points can then be accomplished by the
68 : : //! following loop:
69 : : //! \code{.cpp}
70 : : //! for (std::size_t p=0; p<npoin; ++p)
71 : : //! for (auto i=esup.second[p]+1; i<=esup.second[p+1]; ++i)
72 : : //! use element id esup.first[i]
73 : : //! \endcode
74 : : //! To find out the number of points, _npoin_, the mesh connectivity,
75 : : //! _inpoel_, can be queried:
76 : : //! \code{.cpp}
77 : : //! auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
78 : : //! Assert( *minmax.first == 0, "node ids should start from zero" );
79 : : //! auto npoin = *minmax.second + 1;
80 : : //! \endcode
81 : : //! \note In principle, this function *should* work for any positive nnpe,
82 : : //! however, only nnpe = 4 (tetrahedra) and nnpe = 3 (triangles) are tested.
83 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
84 : : // *****************************************************************************
85 : : {
86 [ + + ][ + - ]: 10283 : Assert( !inpoel.empty(), "Attempt to call genEsup() on empty container" );
[ + - ][ + - ]
87 [ + + ][ + - ]: 10281 : Assert( nnpe > 0, "Attempt to call genEsup() with zero nodes per element" );
[ + - ][ + - ]
88 [ + + ][ + - ]: 10279 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ + - ][ + - ]
89 : :
90 : : // find out number of points in mesh connectivity
91 [ + - ]: 10265 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
92 [ - + ][ - - ]: 10265 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
93 : 10265 : auto npoin = *minmax.second + 1;
94 : :
95 : : // allocate one of the linked lists storing elements surrounding points: esup2
96 : : // fill with zeros
97 [ + - ]: 10265 : std::vector< std::size_t > esup2( npoin+1, 0 );
98 : :
99 : : // element pass 1: count number of elements connected to each point
100 [ + + ]: 19389825 : for (auto n : inpoel) ++esup2[ n + 1 ];
101 : :
102 : : // storage/reshuffling pass 1: update storage counter and store
103 : : // also find out the maximum size of esup1 (mesup)
104 : 10265 : auto mesup = esup2[0]+1;
105 [ + + ]: 1379951 : for (std::size_t i=1; i<npoin+1; ++i) {
106 : 1369686 : esup2[i] += esup2[i-1];
107 [ + - ]: 1369686 : if (esup2[i]+1 > mesup) mesup = esup2[i]+1;
108 : : }
109 : :
110 : : // now we know mesup, so allocate the other one of the linked lists storing
111 : : // elements surrounding points: esup1
112 [ + - ]: 10265 : std::vector< std::size_t > esup1( mesup );
113 : :
114 : : // store the elements in esup1
115 : 10265 : std::size_t e = 0;
116 [ + + ]: 19389825 : for (auto n : inpoel) {
117 : 19379560 : auto j = esup2[n]+1;
118 : 19379560 : esup2[n] = j;
119 : 19379560 : esup1[j] = e/nnpe;
120 : 19379560 : ++e;
121 : : }
122 : :
123 : : // storage/reshuffling pass 2
124 [ + + ]: 1379951 : for (auto i=npoin; i>0; --i) esup2[i] = esup2[i-1];
125 : 10265 : esup2[0] = 0;
126 : :
127 : : // Return (move out) linked lists
128 [ + - ]: 20530 : return std::make_pair( std::move(esup1), std::move(esup2) );
129 : 10265 : }
130 : :
131 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
132 : 4959 : genPsup( const std::vector< std::size_t >& inpoel,
133 : : std::size_t nnpe,
134 : : const std::pair< std::vector< std::size_t >,
135 : : std::vector< std::size_t > >& esup )
136 : : // *****************************************************************************
137 : : // Generate derived data structure, points surrounding points
138 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
139 : : //! node ids of each element of an unstructured mesh. Example:
140 : : //! \code{.cpp}
141 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
142 : : //! 10, 14, 13, 12 };
143 : : //! \endcode
144 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
145 : : //! and { 10, 14, 13, 12 }.
146 : : //! \param[in] nnpe Number of nodes per element
147 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
148 : : //! \return Linked lists storing points surrounding points
149 : : //! \warning It is not okay to call this function with an empty container for
150 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
151 : : //! element; it will throw an exception.
152 : : //! \details The data generated here is stored in a linked list, more precisely,
153 : : //! two linked arrays (vectors), _psup1_ and _psup2_, where _psup2_ holds the
154 : : //! indices at which _psup1_ holds the point ids surrounding points. Looping
155 : : //! over all points surrounding all points can then be accomplished by the
156 : : //! following loop:
157 : : //! \code{.cpp}
158 : : //! for (std::size_t p=0; p<npoin; ++p)
159 : : //! for (auto i=psup.second[p]+1; i<=psup.second[p+1]; ++i)
160 : : //! use point id psup.first[i]
161 : : //! \endcode
162 : : //! To find out the number of points, _npoin_, the mesh connectivity,
163 : : //! _inpoel_, can be queried:
164 : : //! \code{.cpp}
165 : : //! auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
166 : : //! Assert( *minmax.first == 0, "node ids should start from zero" );
167 : : //! auto npoin = *minmax.second + 1;
168 : : //! \endcode
169 : : //! or the length-1 of the generated index list:
170 : : //! \code{.cpp}
171 : : //! auto npoin = psup.second.size()-1;
172 : : //! \endcode
173 : : //! \note In principle, this function *should* work for any positive nnpe,
174 : : //! however, only nnpe = 4 (tetrahedra) and nnpe = 3 (triangles) are tested.
175 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
176 : : // *****************************************************************************
177 : : {
178 [ + + ][ + - ]: 4959 : Assert( !inpoel.empty(), "Attempt to call genPsup() on empty container" );
[ + - ][ + - ]
179 [ + + ][ + - ]: 4957 : Assert( nnpe > 0, "Attempt to call genPsup() with zero nodes per element" );
[ + - ][ + - ]
180 [ - + ][ - - ]: 4955 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
181 [ + + ][ + - ]: 4955 : Assert( !esup.first.empty(), "Attempt to call genPsup() with empty esup1" );
[ + - ][ + - ]
182 [ + + ][ + - ]: 4953 : Assert( !esup.second.empty(), "Attempt to call genPsup() with empty esup2" );
[ + - ][ + - ]
183 : :
184 : : // find out number of points in mesh connectivity
185 [ + - ]: 4951 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
186 [ - + ][ - - ]: 4951 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
187 : 4951 : auto npoin = *minmax.second + 1;
188 : :
189 : 4951 : auto& esup1 = esup.first;
190 : 4951 : auto& esup2 = esup.second;
191 : :
192 : : // allocate both of the linked lists storing points surrounding points, we
193 : : // only know the size of psup2, put in a single zero in psup1
194 [ + - ][ + - ]: 4951 : std::vector< std::size_t > psup2( npoin+1 ), psup1( 1, 0 );
195 : :
196 : : // allocate and fill with zeros a temporary array, only used locally
197 [ + - ]: 4951 : std::vector< std::size_t > lpoin( npoin, 0 );
198 : :
199 : : // fill both psup1 and psup2
200 : 4951 : psup2[0] = 0;
201 : 4951 : std::size_t j = 0;
202 [ + + ]: 636817 : for (std::size_t p=0; p<npoin; ++p) {
203 [ + + ]: 8907142 : for (std::size_t i=esup2[p]+1; i<=esup2[p+1]; ++i ) {
204 [ + + ]: 41376092 : for (std::size_t n=0; n<nnpe; ++n) {
205 : 33100816 : auto q = inpoel[ esup1[i] * nnpe + n ];
206 [ + + ][ + + ]: 33100816 : if (q != p && lpoin[q] != p+1) {
[ + + ]
207 : 6396398 : ++j;
208 [ + - ]: 6396398 : psup1.push_back( q );
209 : 6396398 : lpoin[q] = p+1;
210 : : }
211 : : }
212 : : }
213 : 631866 : psup2[p+1] = j;
214 : : }
215 : :
216 : : // sort point ids for each point in psup1
217 [ + + ]: 636817 : for (std::size_t p=0; p<npoin; ++p)
218 [ + - ][ + - ]: 1263732 : std::sort(
[ + - ]
219 : 631866 : std::next( begin(psup1), static_cast<std::ptrdiff_t>(psup2[p]+1) ),
220 : 631866 : std::next( begin(psup1), static_cast<std::ptrdiff_t>(psup2[p+1]+1) ) );
221 : :
222 : : // Return (move out) linked lists
223 [ + - ]: 9902 : return std::make_pair( std::move(psup1), std::move(psup2) );
224 : 4951 : }
225 : :
226 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
227 : 16 : genEdsup( const std::vector< std::size_t >& inpoel,
228 : : std::size_t nnpe,
229 : : const std::pair< std::vector< std::size_t >,
230 : : std::vector< std::size_t > >& esup )
231 : : // *****************************************************************************
232 : : // Generate derived data structure, edges surrounding points
233 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
234 : : //! node ids of each element of an unstructured mesh. Example:
235 : : //! \code{.cpp}
236 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
237 : : //! 10, 14, 13, 12 };
238 : : //! \endcode
239 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
240 : : //! and { 10, 14, 13, 12 }.
241 : : //! \param[in] nnpe Number of nodes per element (3 or 4)
242 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
243 : : //! \return Linked lists storing edges (point ids p < q) emanating from points
244 : : //! \warning It is not okay to call this function with an empty container for
245 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
246 : : //! element; it will throw an exception.
247 : : //! \details The data generated here is stored in a linked list, more precisely,
248 : : //! two linked arrays (vectors), _edsup1_ and _edsup2_, where _edsup2_ holds
249 : : //! the indices at which _edsup1_ holds the edge-end point ids emanating from
250 : : //! points for all points. The generated data structure, linked lists edsup1
251 : : //! and edsup2, are very similar to psup1 and psup2, generated by genPsup(),
252 : : //! except here only unique edges are stored, i.e., for edges with point ids
253 : : //! p < q, only ids q are stored that are still associated to point p. Looping
254 : : //! over all unique edges can then be accomplished by the following loop:
255 : : //! \code{.cpp}
256 : : //! for (std::size_t p=0; p<npoin; ++p)
257 : : //! for (auto i=edsup.second[p]+1; i<=edsup.second[p+1]; ++i)
258 : : //! use edge with point ids p < edsup.first[i]
259 : : //! \endcode
260 : : //! To find out the number of points, _npoin_, the mesh connectivity,
261 : : //! _inpoel_, can be queried:
262 : : //! \code{.cpp}
263 : : //! auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
264 : : //! Assert( *minmax.first == 0, "node ids should start from zero" );
265 : : //! auto npoin = *minmax.second + 1;
266 : : //! \endcode
267 : : //! \note At first sight, this function seems to work for elements with more
268 : : //! vertices than that of tetrahedra. However, that is not the case since the
269 : : //! algorithm for nnpe > 4 would erronously identify any two combination of
270 : : //! vertices as a valid edge of an element. Since only triangles and
271 : : //! tetrahedra have no internal edges, this algorithm only works for triangle
272 : : //! and tetrahedra element connectivity.
273 : : //! \see tk::genInpoed for similar data that sometimes may be more advantageous
274 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
275 : : // *****************************************************************************
276 : : {
277 [ + + ][ + - ]: 16 : Assert( !inpoel.empty(), "Attempt to call genEdsup() on empty container" );
[ + - ][ + - ]
278 [ + + ][ + - ]: 14 : Assert( nnpe > 0, "Attempt to call genEdsup() with zero nodes per element" );
[ + - ][ + - ]
279 [ + + ][ + + ]: 12 : Assert( nnpe == 3 || nnpe == 4,
[ + - ][ + - ]
[ + - ]
280 : : "Attempt to call genEdsup() with nodes per element, nnpe, that is "
281 : : "neither 4 (tetrahedra) nor 3 (triangles)." );
282 [ - + ][ - - ]: 10 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
283 [ + + ][ + - ]: 10 : Assert( !esup.first.empty(), "Attempt to call genEdsup() with empty esup1" );
[ + - ][ + - ]
284 [ + + ][ + - ]: 8 : Assert( !esup.second.empty(), "Attempt to call genEdsup() with empty esup2" );
[ + - ][ + - ]
285 : :
286 : : // find out number of points in mesh connectivity
287 [ + - ]: 6 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
288 [ - + ][ - - ]: 6 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
289 : 6 : auto npoin = *minmax.second + 1;
290 : :
291 : 6 : const auto& esup1 = esup.first;
292 : 6 : const auto& esup2 = esup.second;
293 : :
294 : : // allocate and fill with zeros a temporary array, only used locally
295 [ + - ]: 6 : std::vector< std::size_t > lpoin( npoin, 0 );
296 : :
297 : : // map to contain stars, a point associated to points connected with edges
298 : : // storing only the end-point id, q, of point ids p < q
299 : 6 : std::map< std::size_t, std::vector< std::size_t > > star;
300 : :
301 : : // generate edge connectivity and store as stars where center id < spike id
302 [ + + ]: 90 : for (std::size_t p=0; p<npoin; ++p)
303 [ + + ]: 588 : for (std::size_t i=esup2[p]+1; i<=esup2[p+1]; ++i )
304 [ + + ]: 2304 : for (std::size_t n=0; n<nnpe; ++n) {
305 : 1800 : auto q = inpoel[ esup1[i] * nnpe + n ];
306 [ + + ][ + + ]: 1800 : if (q != p && lpoin[q] != p+1) {
[ + + ]
307 [ + + ][ + - ]: 510 : if (p < q) star[p].push_back(q);
[ + - ]
308 : 510 : lpoin[q] = p+1;
309 : : }
310 : : }
311 : :
312 : : // linked lists (vectors) to store edges surrounding points and their indices
313 [ + - ][ + - ]: 6 : std::vector< std::size_t > edsup1( 1, 0 ), edsup2( 1, 0 );
314 : :
315 : : // sort non-center points of each star and store nodes and indices in vectors
316 [ + + ]: 69 : for (auto& p : star) {
317 [ + - ]: 63 : std::sort( begin(p.second), end(p.second) );
318 [ + - ]: 63 : edsup2.push_back( edsup2.back() + p.second.size() );
319 [ + - ]: 63 : edsup1.insert( end(edsup1), begin(p.second), end(p.second) );
320 : : }
321 : : // fill up index array with the last index for points with no new edges
322 [ + + ]: 27 : for (std::size_t i=0; i<npoin-star.size(); ++i)
323 [ + - ]: 21 : edsup2.push_back( edsup2.back() );
324 : :
325 : : // Return (move out) linked lists
326 [ + - ]: 12 : return std::make_pair( std::move(edsup1), std::move(edsup2) );
327 : 6 : }
328 : :
329 : : std::vector< std::size_t >
330 : 30 : genInpoed( const std::vector< std::size_t >& inpoel,
331 : : std::size_t nnpe,
332 : : const std::pair< std::vector< std::size_t >,
333 : : std::vector< std::size_t > >& esup )
334 : : // *****************************************************************************
335 : : // Generate derived data structure, edge connectivity
336 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
337 : : //! node ids of each element of an unstructured mesh. Example:
338 : : //! \code{.cpp}
339 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
340 : : //! 10, 14, 13, 12 };
341 : : //! \endcode
342 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
343 : : //! and { 10, 14, 13, 12 }.
344 : : //! \param[in] nnpe Number of nodes per element (3 or 4)
345 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
346 : : //! \return Linear vector storing edge connectivity (point ids p < q)
347 : : //! \warning It is not okay to call this function with an empty container for
348 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
349 : : //! element; it will throw an exception.
350 : : //! \details The data generated here is stored in a linear vector and is very
351 : : //! similar to the linked lists, _edsup1_ and _edsup2, generated by
352 : : //! genEdsup(). The difference is that in the linear vector, inpoed, generated
353 : : //! here, both edge point ids are stored as a pair, p < q, as opposed to the
354 : : //! linked lists edsup1 and edsup2, in which edsup1 only stores the edge-end
355 : : //! point ids (still associated to edge-start point ids when used together
356 : : //! with edsup2). The rationale is that while inpoed is larger in memory, it
357 : : //! allows direct access to edges (pair of point ids making up an edge),
358 : : //! edsup1 and edsup2 are smaller in memory, still allow accessing the same
359 : : //! data (edge point id pairs) but only in a linear fashion, not by direct
360 : : //! access to particular edges. Accessing all unique edges using the edge
361 : : //! connectivity data structure, inpoed, generated here can be accomplished by
362 : : //! \code{.cpp}
363 : : //! for (std::size_t e=0; e<inpoed.size()/2; ++e) {
364 : : //! use point id p of edge e = inpoed[e*2];
365 : : //! use point id q of edge e = inpoed[e*2+1];
366 : : //! }
367 : : //! \endcode
368 : : //! \note At first sight, this function seems to work for elements with more
369 : : //! vertices than that of tetrahedra. However, that is not the case since the
370 : : //! algorithm for nnpe > 4 would erronously identify any two combination of
371 : : //! vertices as a valid edge of an element. Since only triangles and
372 : : //! tetrahedra have no internal edges, this algorithm only works for triangle
373 : : //! and tetrahedra element connectivity.
374 : : //! \see tk::genEdsup for similar data that sometimes may be more advantageous
375 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
376 : : // *****************************************************************************
377 : : {
378 [ + + ][ + - ]: 30 : Assert( !inpoel.empty(), "Attempt to call genInpoed() on empty container" );
[ + - ][ + - ]
379 [ + + ][ + - ]: 28 : Assert( nnpe > 0, "Attempt to call genInpoed() with zero nodes per element" );
[ + - ][ + - ]
380 [ + + ][ + + ]: 26 : Assert( nnpe == 3 || nnpe == 4,
[ + - ][ + - ]
[ + - ]
381 : : "Attempt to call genInpoed() with nodes per element, nnpe, that is "
382 : : "neither 4 (tetrahedra) nor 3 (triangles)." );
383 [ - + ][ - - ]: 24 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
384 [ + + ][ + - ]: 24 : Assert( !esup.first.empty(), "Attempt to call genInpoed() with empty esup1" );
[ + - ][ + - ]
385 [ + + ][ + - ]: 22 : Assert( !esup.second.empty(),
[ + - ][ + - ]
386 : : "Attempt to call genInpoed() with empty esup2" );
387 : :
388 : : // find out number of points in mesh connectivity
389 [ + - ]: 20 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
390 [ - + ][ - - ]: 20 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
391 : 20 : auto npoin = *minmax.second + 1;
392 : :
393 : 20 : const auto& esup1 = esup.first;
394 : 20 : const auto& esup2 = esup.second;
395 : :
396 : : // allocate and fill with zeros a temporary array, only used locally
397 [ + - ]: 20 : std::vector< std::size_t > lpoin( npoin, 0 );
398 : :
399 : : // map to contain stars, a point associated to points connected with edges,
400 : : // storing only the end-point id, q, of point ids p < q
401 : 20 : std::map< std::size_t, std::vector< std::size_t > > star;
402 : :
403 : : // generate edge connectivity and store as stars where center id < spike id
404 [ + + ]: 240 : for (std::size_t p=0; p<npoin; ++p)
405 [ + + ]: 1492 : for (std::size_t i=esup2[p]+1; i<=esup2[p+1]; ++i )
406 [ + + ]: 6072 : for (std::size_t n=0; n<nnpe; ++n) {
407 : 4800 : auto q = inpoel[ esup1[i] * nnpe + n ];
408 [ + + ][ + + ]: 4800 : if (q != p && lpoin[q] != p+1) {
[ + + ]
409 [ + + ][ + - ]: 1340 : if (p < q) star[p].push_back( q );
[ + - ]
410 : 1340 : lpoin[q] = p+1;
411 : : }
412 : : }
413 : :
414 : : // linear vector to store edge connectivity and their indices
415 : 20 : std::vector< std::size_t > inpoed;
416 : :
417 : : // sort non-center points of each star and store both start and end points of
418 : : // each star in linear vector
419 [ + + ]: 200 : for (auto& p : star) {
420 [ + - ]: 180 : std::sort( begin(p.second), end(p.second) );
421 [ + + ]: 850 : for (auto e : p.second) {
422 [ + - ]: 670 : inpoed.push_back( p.first );
423 [ + - ]: 670 : inpoed.push_back( e );
424 : : }
425 : : }
426 : :
427 : : // Return (move out) linear vector
428 : 40 : return inpoed;
429 : 20 : }
430 : :
431 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
432 : 14 : genEsupel( const std::vector< std::size_t >& inpoel,
433 : : std::size_t nnpe,
434 : : const std::pair< std::vector< std::size_t >,
435 : : std::vector< std::size_t > >& esup )
436 : : // *****************************************************************************
437 : : // Generate derived data structure, elements surrounding points of elements
438 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
439 : : //! node ids of each element of an unstructured mesh. Example:
440 : : //! \code{.cpp}
441 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
442 : : //! 10, 14, 13, 12 };
443 : : //! \endcode
444 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
445 : : //! and { 10, 14, 13, 12 }.
446 : : //! \param[in] nnpe Number of nodes per element
447 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
448 : : //! \return Linked lists storing elements surrounding points of elements
449 : : //! \warning It is not okay to call this function with an empty container for
450 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
451 : : //! element; it will throw an exception.
452 : : //! \details The data generated here is stored in a linked list, more precisely,
453 : : //! two linked arrays (vectors), _esupel1_ and _esupel2_, where _esupel2_
454 : : //! holds the indices at which _esupel1_ holds the element ids surrounding
455 : : //! points of elements. Looping over all elements surrounding the points of
456 : : //! all elements can then be accomplished by the following loop:
457 : : //! \code{.cpp}
458 : : //! for (std::size_t e=0; e<nelem; ++e)
459 : : //! for (auto i=esupel.second[e]+1; i<=esupel.second[e+1]; ++i)
460 : : //! use element id esupel.first[i]
461 : : //! \endcode
462 : : //! To find out the number of elements, _nelem_, the size of the mesh
463 : : //! connectivity vector, _inpoel_, can be devided by the number of nodes per
464 : : //! elements, _nnpe_:
465 : : //! \code{.cpp}
466 : : //! auto nelem = inpoel.size()/nnpe;
467 : : //! \endcode
468 : : //! \note In principle, this function *should* work for any positive nnpe,
469 : : //! however, only nnpe = 4 (tetrahedra) and nnpe = 3 (triangles) are tested.
470 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
471 : : // *****************************************************************************
472 : : {
473 [ + + ][ + - ]: 14 : Assert( !inpoel.empty(), "Attempt to call genEsupel() on empty container" );
[ + - ][ + - ]
474 [ + + ][ + - ]: 12 : Assert( nnpe > 0, "Attempt to call genEsupel() with zero nodes per element" );
[ + - ][ + - ]
475 [ - + ][ - - ]: 10 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
476 [ + + ][ + - ]: 10 : Assert( !esup.first.empty(), "Attempt to call genEsupel() with empty esup1" );
[ + - ][ + - ]
477 [ + + ][ + - ]: 8 : Assert( !esup.second.empty(),
[ + - ][ + - ]
478 : : "Attempt to call genEsupel() with empty esup2" );
479 : :
480 : 6 : const auto& esup1 = esup.first;
481 : 6 : const auto& esup2 = esup.second;
482 : :
483 : : // linked lists storing elements surrounding points of elements, put in a
484 : : // single zero in both
485 [ + - ][ + - ]: 6 : std::vector< std::size_t > esupel2( 1, 0 ), esupel1( 1, 0 );
486 : :
487 : 6 : std::size_t e = 0;
488 : 6 : std::set< std::size_t > esuel;
489 [ + + ]: 510 : for (auto p : inpoel) { // loop over all points of all elements
490 : : // collect unique element ids of elements surrounding points of element
491 [ + - ][ + + ]: 4104 : for (auto i=esup2[p]+1; i<=esup2[p+1]; ++i) esuel.insert( esup1[i] );
492 [ + + ]: 504 : if (++e%nnpe == 0) { // when finished checking all nodes of element
493 : : // erase element whose surrounding elements are considered
494 [ + - ]: 144 : esuel.erase( e/nnpe-1 );
495 : : // store unique element ids in esupel1
496 [ + - ]: 144 : esupel1.insert( end(esupel1), begin(esuel), end(esuel) );
497 : : // store end-index for element used to address into esupel1
498 [ + - ]: 144 : esupel2.push_back( esupel2.back() + esuel.size() );
499 : 144 : esuel.clear();
500 : : }
501 : : }
502 : :
503 : : // Return (move out) linked lists
504 [ + - ]: 12 : return std::make_pair( std::move(esupel1), std::move(esupel2) );
505 : 6 : }
506 : :
507 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
508 : 14 : genEsuel( const std::vector< std::size_t >& inpoel,
509 : : std::size_t nnpe,
510 : : const std::pair< std::vector< std::size_t >,
511 : : std::vector< std::size_t > >& esup )
512 : : // *****************************************************************************
513 : : // Generate derived data structure, elements surrounding elements
514 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
515 : : //! node ids of each element of an unstructured mesh. Example:
516 : : //! \code{.cpp}
517 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
518 : : //! 10, 14, 13, 12 };
519 : : //! \endcode
520 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
521 : : //! and { 10, 14, 13, 12 }.
522 : : //! \param[in] nnpe Number of nodes per element
523 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
524 : : //! \return Linked lists storing elements surrounding elements
525 : : //! \warning It is not okay to call this function with an empty container for
526 : : //! inpoel or esup.first or esup.second; it will throw an exception.
527 : : //! \details The data generated here is stored in a linked list, more precisely,
528 : : //! two linked arrays (vectors), _esuel1_ and _esuel2_, where _esuel2_ holds
529 : : //! the indices at which _esuel1_ holds the element ids surrounding elements.
530 : : //! Looping over elements surrounding elements can then be accomplished by the
531 : : //! following loop:
532 : : //! \code{.cpp}
533 : : //! for (std::size_t e=0; e<nelem; ++e)
534 : : //! for (auto i=esuel.second[e]+1; i<=esuel.second[e+1]; ++i)
535 : : //! use element id esuel.first[i]
536 : : //! \endcode
537 : : //! To find out the number of elements, _nelem_, the size of the mesh
538 : : //! connectivity vector, _inpoel_, can be devided by the number of nodes per
539 : : //! elements, _nnpe_:
540 : : //! \code{.cpp}
541 : : //! auto nelem = inpoel.size()/nnpe;
542 : : //! \endcode
543 : : //! \note In principle, this function *should* work for any positive nnpe,
544 : : //! however, only nnpe = 4 (tetrahedra) and nnpe = 3 (triangles) are tested.
545 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
546 : : // *****************************************************************************
547 : : {
548 [ + + ][ + - ]: 14 : Assert( !inpoel.empty(), "Attempt to call genEsuel() on empty container" );
[ + - ][ + - ]
549 [ + + ][ + - ]: 12 : Assert( nnpe > 0, "Attempt to call genEsuel() with zero nodes per element" );
[ + - ][ + - ]
550 [ - + ][ - - ]: 10 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by four" );
[ - - ][ - - ]
551 [ + + ][ + - ]: 10 : Assert( !esup.first.empty(), "Attempt to call genEsuel() with empty esuel1" );
[ + - ][ + - ]
552 [ + + ][ + - ]: 8 : Assert( !esup.second.empty(),
[ + - ][ + - ]
553 : : "Attempt to call genEsuel() with empty esuel2" );
554 : :
555 : 6 : const auto& esup1 = esup.first;
556 : 6 : const auto& esup2 = esup.second;
557 : :
558 : 6 : auto nelem = inpoel.size()/nnpe;
559 : :
560 : : // lambda that returns 1 if elements hel and gel share a face
561 : 185904 : auto adj = [ &inpoel, nnpe ]( std::size_t hel, std::size_t gel ) -> bool {
562 : 3600 : std::size_t sp = 0;
563 [ + + ]: 16848 : for (std::size_t h=0; h<nnpe; ++h)
564 [ + + ]: 62784 : for (std::size_t g=0; g<nnpe; ++g)
565 [ + + ]: 49536 : if (inpoel[hel*nnpe+h] == inpoel[gel*nnpe+g]) ++sp;
566 [ + + ]: 3600 : if (sp == nnpe-1) return true; else return false;
567 : 6 : };
568 : :
569 : : // map to associate unique elements and their surrounding elements
570 : 6 : std::map< std::size_t, std::vector< std::size_t > > es;
571 : :
572 [ + + ]: 150 : for (std::size_t e=0; e<nelem; ++e) {
573 : 144 : std::set< std::size_t > faces; // will collect elem ids of shared faces
574 [ + + ]: 648 : for (std::size_t n=0; n<nnpe; ++n) {
575 : 504 : auto i = inpoel[ e*nnpe+n ];
576 [ + + ]: 4104 : for (auto j=esup2[i]+1; j<=esup2[i+1]; ++j)
577 [ + + ][ + - ]: 3600 : if (adj( e, esup1[j] )) faces.insert( esup1[j] );
578 : : }
579 : : // store element ids of shared faces
580 [ + - ][ + - ]: 576 : for (auto j : faces) es[e].push_back(j);
[ + + ]
581 : 144 : }
582 : :
583 : : // storing elements surrounding elements
584 [ + - ][ + - ]: 6 : std::vector< std::size_t > esuel1( 1, 0 ), esuel2( 1, 0 );
585 : :
586 : : // store elements surrounding elements in linked lists
587 [ + + ]: 150 : for (const auto& e : es) {
588 [ + - ]: 144 : esuel2.push_back( esuel2.back() + e.second.size() );
589 [ + - ]: 144 : esuel1.insert( end(esuel1), begin(e.second), end(e.second) );
590 : : }
591 : :
592 : : // Return (move out) linked lists
593 [ + - ]: 12 : return std::make_pair( std::move(esuel1), std::move(esuel2) );
594 : 6 : }
595 : :
596 : : std::vector< std::size_t >
597 : 12 : genInedel( const std::vector< std::size_t >& inpoel,
598 : : std::size_t nnpe,
599 : : const std::vector< std::size_t >& inpoed )
600 : : // *****************************************************************************
601 : : // Generate derived data structure, edges of elements
602 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
603 : : //! node ids of each element of an unstructured mesh. Example:
604 : : //! \code{.cpp}
605 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
606 : : //! 10, 14, 13, 12 };
607 : : //! \endcode
608 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
609 : : //! and { 10, 14, 13, 12 }.
610 : : //! \param[in] nnpe Number of nodes per element
611 : : //! \param[in] inpoed Edge connectivity as linear vector, see tk::genInpoed
612 : : //! \return Linear vector storing all edge ids * 2 of all elements
613 : : //! \warning It is not okay to call this function with an empty container for
614 : : //! inpoel or inpoed or a non-positive number of nodes per element; it will
615 : : //! throw an exception.
616 : : //! \details The data generated here is stored in a linear vector with all
617 : : //! edge ids (as defined by inpoed) of all elements. The edge ids stored in
618 : : //! inedel can be directly used to index the vector inpoed. Because the
619 : : //! derived data structure generated here, inedel, is intended to be used in
620 : : //! conjunction with the linear vector inpoed and not with the linked lists
621 : : //! edsup1 and edsup2, this function takes inpoed as an argument. Accessing
622 : : //! the edges of element e using the edge of elements data structure, inedel,
623 : : //! generated here can be accomplished by
624 : : //! \code{.cpp}
625 : : //! for (std::size_t e=0; e<nelem; ++e) {
626 : : //! for (std::size_t i=0; i<nepe; ++i) {
627 : : //! use edge id inedel[e*nepe+i] of element e, or
628 : : //! use point ids p < q of edge id inedel[e*nepe+i] of element e as
629 : : //! p = inpoed[ inedel[e*nepe+i]*2 ]
630 : : //! q = inpoed[ inedel[e*nepe+i]*2+1 ]
631 : : //! }
632 : : //! }
633 : : //! \endcode
634 : : //! where _nepe_ denotes the number of edges per elements: 3 for triangles, 6
635 : : //! for tetrahedra. To find out the number of elements, _nelem_, the size of
636 : : //! the mesh connectivity vector, _inpoel_, can be devided by the number of
637 : : //! nodes per elements, _nnpe_:
638 : : //! \code{.cpp}
639 : : //! auto nelem = inpoel.size()/nnpe;
640 : : //! \endcode
641 : : //! \note At first sight, this function seems to work for elements with more
642 : : //! vertices than that of tetrahedra. However, that is not the case since the
643 : : //! algorithm for nnpe > 4 would erronously identify any two combination of
644 : : //! vertices as a valid edge of an element. Since only triangles and
645 : : //! tetrahedra have no internal edges, this algorithm only works for triangle
646 : : //! and tetrahedra element connectivity.
647 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
648 : : // *****************************************************************************
649 : : {
650 [ + + ][ + - ]: 12 : Assert( !inpoel.empty(), "Attempt to call genInedel() on empty container" );
[ + - ][ + - ]
651 [ + + ][ + - ]: 10 : Assert( nnpe > 0, "Attempt to call genInedel() with zero nodes per element" );
[ + - ][ + - ]
652 [ + + ][ + + ]: 6 : Assert( nnpe == 3 || nnpe == 4,
[ + - ][ + - ]
[ + - ]
653 : : "Attempt to call genInedel() with nodes per element, nnpe, that is "
654 : : "neither 4 (tetrahedra) nor 3 (triangles)." );
655 [ - + ][ - - ]: 4 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
656 [ - + ][ - - ]: 4 : Assert( !inpoed.empty(), "Attempt to call genInedel() with empty inpoed" );
[ - - ][ - - ]
657 : :
658 : : // find out number of points in mesh connectivity
659 [ + - ]: 4 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
660 [ - + ][ - - ]: 4 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
661 : 4 : auto npoin = *minmax.second + 1;
662 : :
663 : : // First, generate index of star centers. This is necessary to avoid a
664 : : // brute-force search for point ids of edges when searching for element edges.
665 : : // Note that this is the same as edsup2, generated by genEdsup(). However,
666 : : // because the derived data structure generated here, inedel, is intended to
667 : : // be used in conjunction with the linear vector inpoed and not with the
668 : : // linked lists edsup1 and edsup2, this function takes inpoed as an argument,
669 : : // and so edsup2 is temporarily generated here to avoid a brute-force search.
670 : :
671 : : // map to contain stars, a point associated to points connected with edges
672 : : // storing only the end-point id, q, of point ids p < q
673 : 4 : std::map< std::size_t, std::vector< std::size_t > > star;
674 : :
675 : : // generate stars from inpoed; starting with zero, every even is a star
676 : : // center, every odd is a spike
677 [ + + ]: 174 : for (std::size_t i=0; i<inpoed.size()/2; ++i)
678 [ + - ][ + - ]: 170 : star[ inpoed[i*2] ].push_back( inpoed[i*2+1] );
679 : :
680 : : // store index of star centers in vector; assume non-center points of each
681 : : // star have already been sorted
682 [ + - ]: 4 : std::vector< std::size_t > edsup2( 1, 0 );
683 [ + - ][ + + ]: 46 : for (const auto& p : star) edsup2.push_back(edsup2.back() + p.second.size());
684 : : // fill up index array with the last index for points with no new edges
685 [ + + ]: 18 : for (std::size_t i=0; i<npoin-star.size(); ++i)
686 [ + - ]: 14 : edsup2.push_back( edsup2.back() );
687 : 4 : star.clear();
688 : :
689 : : // Second, generate edges of elements
690 : :
691 : 4 : auto nelem = inpoel.size()/nnpe;
692 : :
693 : : // map associating elem id with vector of edge ids
694 : 4 : std::map< std::size_t, std::vector< std::size_t > > edges;
695 : :
696 : : // generate map of elements associated to edge ids
697 [ + + ]: 100 : for (std::size_t e=0; e<nelem; ++e)
698 [ + + ]: 432 : for (std::size_t n=0; n<nnpe; ++n) {
699 : 336 : auto p = inpoel[e*nnpe+n];
700 [ + + ]: 1324 : for (auto i=edsup2[p]+1; i<=edsup2[p+1]; ++i)
701 [ + + ]: 4508 : for (std::size_t j=0; j<nnpe; ++j)
702 [ + + ]: 3520 : if (inpoed[(i-1)*2+1] == inpoel[e*nnpe+j])
703 [ + - ][ + - ]: 432 : edges[e].push_back( i-1 );
704 : : }
705 : :
706 : : // linear vector to store the edge ids of all elements
707 [ + - ]: 4 : std::vector< std::size_t > inedel( sumvalsize(edges) );
708 : :
709 : : // store edge ids of elements in linear vector
710 : 4 : std::size_t j = 0;
711 [ + + ][ + + ]: 532 : for (const auto& e : edges) for (auto p : e.second) inedel[ j++ ] = p;
712 : :
713 : : // Return (move out) vector
714 : 8 : return inedel;
715 : 4 : }
716 : :
717 : : std::unordered_map< UnsMesh::Edge, std::vector< std::size_t >,
718 : : UnsMesh::Hash<2>, UnsMesh::Eq<2> >
719 : 15 : genEsued( const std::vector< std::size_t >& inpoel,
720 : : std::size_t nnpe,
721 : : const std::pair< std::vector< std::size_t >,
722 : : std::vector< std::size_t > >& esup )
723 : : // *****************************************************************************
724 : : // Generate derived data structure, elements surrounding edges
725 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
726 : : //! node ids of each element of an unstructured mesh. Example:
727 : : //! \code{.cpp}
728 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
729 : : //! 10, 14, 13, 12 };
730 : : //! \endcode
731 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
732 : : //! and { 10, 14, 13, 12 }.
733 : : //! \param[in] nnpe Number of nodes per element (3 or 4)
734 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
735 : : //! \return Associative container storing elements surrounding edges (value),
736 : : //! assigned to edge-end points (key)
737 : : //! \warning It is not okay to call this function with an empty container for
738 : : //! inpoel or esup.first or esup.second or a non-positive number of nodes per
739 : : //! element; it will throw an exception.
740 : : //! \details Looping over elements surrounding all edges can be accomplished by
741 : : //! the following loop:
742 : : //! \code{.cpp}
743 : : //! for (const auto& [edge,surr_elements] : esued) {
744 : : //! use element edge-end-point ids edge[0] and edge[1]
745 : : //! for (auto e : surr_elements) {
746 : : //! use element id e
747 : : //! }
748 : : //! }
749 : : //! \endcode
750 : : //! esued.size() equals the number of edges.
751 : : //! \note At first sight, this function seems to work for elements with more
752 : : //! vertices than that of tetrahedra. However, that is not the case since the
753 : : //! algorithm for nnpe > 4 would erronously identify any two combination of
754 : : //! vertices as a valid edge of an element. Since only triangles and
755 : : //! tetrahedra have no internal edges, this algorithm only works for triangle
756 : : //! and tetrahedra element connectivity.
757 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
758 : : // *****************************************************************************
759 : : {
760 [ + + ][ + - ]: 15 : Assert( !inpoel.empty(), "Attempt to call genEsued() on empty container" );
[ + - ][ + - ]
761 [ + + ][ + - ]: 13 : Assert( nnpe > 0, "Attempt to call genEsued() with zero nodes per element" );
[ + - ][ + - ]
762 [ + + ][ + + ]: 11 : Assert( nnpe == 3 || nnpe == 4,
[ + - ][ + - ]
[ + - ]
763 : : "Attempt to call genEsued() with nodes per element, nnpe, that is "
764 : : "neither 4 (tetrahedra) nor 3 (triangles)." );
765 [ - + ][ - - ]: 9 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by nnpe" );
[ - - ][ - - ]
766 [ + + ][ + - ]: 9 : Assert( !esup.first.empty(), "Attempt to call genEsued() with empty esup1" );
[ + - ][ + - ]
767 [ + + ][ + - ]: 7 : Assert( !esup.second.empty(), "Attempt to call genEsued() with empty esup2" );
[ + - ][ + - ]
768 : :
769 : : // find out number of points in mesh connectivity
770 [ + - ]: 5 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
771 [ - + ][ - - ]: 5 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
772 : 5 : auto npoin = *minmax.second + 1;
773 : :
774 : 5 : const auto& esup1 = esup.first;
775 : 5 : const auto& esup2 = esup.second;
776 : :
777 : : // allocate and fill with zeros a temporary array, only used locally
778 [ + - ]: 5 : std::vector< std::size_t > lpoin( npoin, 0 );
779 : :
780 : : // lambda that returns true if element e contains edge (p < q)
781 : 15162 : auto has = [ &inpoel, nnpe ]( std::size_t e, std::size_t p, std::size_t q ) {
782 : 1266 : int sp = 0;
783 [ + + ]: 5898 : for (std::size_t n=0; n<nnpe; ++n)
784 [ + + ][ + + ]: 4632 : if (inpoel[e*nnpe+n] == p || inpoel[e*nnpe+n] == q) ++sp;
[ + + ]
785 : 1266 : return sp == 2;
786 : 5 : };
787 : :
788 : : // map to associate edges to unique surrounding element ids
789 : : std::unordered_map< UnsMesh::Edge, std::vector< std::size_t >,
790 : 5 : UnsMesh::Hash<2>, UnsMesh::Eq<2> > esued;
791 : :
792 : : // generate edges and associated vector of unique surrounding element ids
793 [ + + ]: 75 : for (std::size_t p=0; p<npoin; ++p)
794 [ + + ]: 502 : for (std::size_t i=esup2[p]+1; i<=esup2[p+1]; ++i )
795 [ + + ]: 2016 : for (std::size_t n=0; n<nnpe; ++n) {
796 : 1584 : auto q = inpoel[ esup1[i] * nnpe + n ];
797 [ + + ][ + + ]: 1584 : if (q != p && lpoin[q] != p+1) {
[ + + ]
798 [ + + ]: 438 : if (p < q) { // for edge given point ids p < q
799 [ + + ]: 1485 : for (std::size_t j=esup2[p]+1; j<=esup2[p+1]; ++j ) {
800 : 1266 : auto e = esup1[j];
801 [ + + ][ + - ]: 1266 : if (has(e,p,q)) esued[{p,q}].push_back(e);
[ + - ]
802 : : }
803 : : }
804 : 438 : lpoin[q] = p+1;
805 : : }
806 : : }
807 : :
808 : : // sort element ids surrounding edges for each edge
809 [ + - ][ + + ]: 224 : for (auto& p : esued) std::sort( begin(p.second), end(p.second) );
810 : :
811 : : // Return elements surrounding edges data structure
812 : 10 : return esued;
813 : 5 : }
814 : :
815 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
816 : 26 : genEdpas( int mvecl, std::size_t nnpe, std::size_t npoin,
817 : : const std::vector< std::size_t >& inpoed )
818 : : // *****************************************************************************
819 : : // Generate vector-groups for edges
820 : : //! \param[in] mvecl Max vector length to target
821 : : //! \param[in] nnpe Number of nodes per (super-)edge
822 : : //! \param[in] npoin Number mesh points
823 : : //! \param[in] inpoed Edge connectivity as linear vector, see tk::genInpoed for
824 : : //! nnpe=2
825 : : //! \return Linked lists storing edge-groups so that any point of a group is
826 : : //! accessed only once within a group.
827 : : //! \warning It is not okay to call this function with an empty container or a
828 : : //! non-positive number of nodes per element; it will throw an exception.
829 : : //! \details The data generated here is stored in a linked list, more precisely,
830 : : //! two linked arrays (vectors), _edpas1_ and _edpas2_, where _edpas2_ holds
831 : : //! the indices at which _edpas1_ holds the edge ids of a vector group.
832 : : //! Looping over all groups can then be accomplished by the following loop:
833 : : //! \code{.cpp}
834 : : //! for (std::size_t w=0; w<edpas.second.size()-1; ++w)
835 : : //! for (auto i=edpas.second[w]+1; i<=edpas.second[w+1]; ++i)
836 : : //! use edge id edpas.first[i]
837 : : //! \endcode
838 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
839 : : // *****************************************************************************
840 : : {
841 [ + + ]: 26 : if (inpoed.empty()) return {};
842 : :
843 [ + + ][ + - ]: 24 : Assert( mvecl > 0, "Attempt to call genEdpas() with non-positive veclen" );
[ + - ][ + - ]
844 [ + + ][ + - ]: 22 : Assert( nnpe > 0, "Attempt to call genEdpas() with non-positive nnpe" );
[ + - ][ + - ]
845 [ + + ][ + - ]: 20 : Assert( npoin > 0, "Attempt to call genEdpas() with non-positive npoin" );
[ + - ][ + - ]
846 [ + + ][ + - ]: 18 : Assert( inpoed.size()%nnpe == 0, "Size of inpoed must be divisible by nnpe" );
[ + - ][ + - ]
847 : :
848 : 16 : auto nedge = inpoed.size() / nnpe;
849 : :
850 [ + - ]: 16 : std::vector< std::size_t > ledge( nedge, 0 );
851 [ + - ]: 16 : std::vector< std::size_t > lpoin( npoin, 0 );
852 : :
853 : 16 : std::pair< std::vector< std::size_t >, std::vector< std::size_t > > edpas;
854 [ + - ]: 16 : edpas.first.resize( nedge+1, 0 );
855 [ + - ]: 16 : edpas.second.push_back( 0 );
856 : :
857 [ + - ]: 16 : std::unordered_set< std::size_t > unedge( nedge );
858 [ + - ][ + + ]: 288 : for (std::size_t e=0; e<nedge; ++e) unedge.insert( e );
859 : :
860 : 16 : std::size_t nenew = 0, ngrou = 0;
861 : :
862 [ + + ]: 98 : while (nenew < nedge) {
863 : 82 : int nvecl = 0;
864 : 82 : ++ngrou;
865 [ + - ]: 82 : edpas.second.emplace_back();
866 [ + + ]: 910 : for (auto ie = begin(unedge); ie != end(unedge); ) {
867 : 852 : auto e = *ie;
868 : 852 : const auto N = inpoed.data() + e*nnpe;
869 : 852 : std::size_t nsw = 0;
870 [ + + ]: 1904 : for (std::size_t i=0; i<nnpe; ++i) {
871 [ + + ]: 1632 : if (lpoin[N[i]] == ngrou) break; else ++nsw;
872 : : }
873 [ + + ]: 852 : if (nsw == nnpe) {
874 [ + + ]: 868 : for (std::size_t i=0; i<nnpe; ++i) lpoin[N[i]] = ngrou;
875 : 272 : ledge[e] = ngrou;
876 : 272 : ++nenew;
877 : 272 : ++nvecl;
878 : 272 : edpas.first[nenew] = e;
879 : 272 : edpas.second[ngrou] = nenew;
880 [ + - ]: 272 : ie = unedge.erase( ie );
881 : : } else {
882 : 580 : ++ie;
883 : : }
884 [ + + ]: 852 : if (nvecl == mvecl) break;
885 : : }
886 : : }
887 : :
888 : : //std::size_t ne = 0;
889 : : //for (std::size_t i=0; i<edpas.second.size()-1; ++i) {
890 : : // std::cout << i+1 << ": " << edpas.second[i+1] - edpas.second[i] << '\n';
891 : : // ne += edpas.second[i+1] - edpas.second[i];
892 : : //}
893 : : //std::cout << "edges grouped: " << ne << " of " << nedge << '\n';
894 : : //std::cout << "edge groups:";
895 : : //for (std::size_t g=1; g<=edpas.second.size(); ++g) {
896 : : // std::cout << '\n';
897 : : // for (std::size_t e=0; e<nedge; ++e) {
898 : : // if (ledge[e] == g) {
899 : : // const auto N = inpoed.data() + e*nnpe;
900 : : // std::cout << e << ":\t";
901 : : // for (std::size_t i=0; i<nnpe-1; ++i) std::cout << N[i] << '-';
902 : : // std::cout << N[nnpe-1] << "\tgrp: " << ledge[e] << '\n';
903 : : // }
904 : : // }
905 : : //}
906 : : //std::cout << '\n';
907 : :
908 : : //std::cout << "\nnew access loop:\n";
909 : : //for (std::size_t g=0; g<edpas.second.size()-1; ++g) {
910 : : // //#pragma omp simd
911 : : // for (auto w=edpas.second[g]+1; w<=edpas.second[g+1]; ++w) {
912 : : // auto e = edpas.first[w];
913 : : // const auto N = inpoed.data() + e*nnpe;
914 : : // for (std::size_t i=0; i<nnpe-1; ++i) std::cout << N[i] << '-';
915 : : // std::cout << N[nnpe-1] << ' ';
916 : : // }
917 : : // std::cout << '\n';
918 : : //}
919 : : //std::cout << '\n';
920 : :
921 : : //std::cout << "old access loop:\n";
922 : : //for (std::size_t e=0; e<nedge; ++e) {
923 : : // const auto N = inpoed.data() + e*nnpe;
924 : : // for (std::size_t i=0; i<nnpe-1; ++i) std::cout << N[i] << '-';
925 : : // std::cout << N[nnpe-1] << ' ';
926 : : //}
927 : : //std::cout << "\n\n";
928 : :
929 : : // Return linked lists
930 : 16 : return edpas;
931 : 16 : }
932 : :
933 : : std::size_t
934 : 2 : genNbfacTet( std::size_t tnbfac,
935 : : const std::vector< std::size_t >& inpoel,
936 : : const std::vector< std::size_t >& triinpoel_complete,
937 : : const std::map< int, std::vector< std::size_t > >& bface_complete,
938 : : const std::unordered_map< std::size_t, std::size_t >& lid,
939 : : std::vector< std::size_t >& triinpoel,
940 : : std::map< int, std::vector< std::size_t > >& bface )
941 : : // *****************************************************************************
942 : : // Generate number of boundary-faces and triangle boundary-face connectivity
943 : : //! \param[in] tnbfac Total number of boundary faces in the entire mesh.
944 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
945 : : //! node ids of each element of an unstructured mesh.
946 : : //! \param[in] triinpoel_complete Interconnectivity of points and boundary-face
947 : : //! in the entire mesh.
948 : : //! \param[in] bface_complete Map of boundary-face lists mapped to corresponding
949 : : //! side set ids for the entire mesh.
950 : : //! \param[in] lid Mapping between the node indices used in the smaller inpoel
951 : : //! connectivity (a subset of the entire triinpoel_complete connectivity),
952 : : //! e.g., after mesh partitioning.
953 : : //! \param[inout] triinpoel Interconnectivity of points and boundary-face in
954 : : //! this mesh-partition.
955 : : //! \param[inout] bface Map of boundary-face lists mapped to corresponding
956 : : //! side set ids for this mesh-partition
957 : : //! \return Number of boundary-faces on this chare/mesh-partition.
958 : : //! \details This function takes a mesh by its domain-element
959 : : //! (tetrahedron-connectivity) in inpoel and a boundary-face (triangle)
960 : : //! connectivity in triinpoel_complete. Based on these two arrays, it
961 : : //! searches for those faces of triinpoel_complete that are also in inpoel
962 : : //! and as a result it generates (1) the number of boundary faces shared with
963 : : //! the mesh in inpoel and (2) the intersection of the triangle element
964 : : //! connectivity whose faces are shared with inpoel. An example use case is
965 : : //! where triinpoel_complete contains the connectivity for the boundary of the
966 : : //! full problem/mesh and inpoel contains the connectivity for only a chunk of
967 : : //! an already partitioned mesh. This function then intersects
968 : : //! triinpoel_complete with inpoel and returns only those faces that share
969 : : //! nodes with inpoel.
970 : : //! \warning This is for Triangular face-elements only.
971 : : // *****************************************************************************
972 : : {
973 : : // cppcheck-suppress unreadVariable
974 : 2 : std::size_t nbfac = 0, nnpf = 3;
975 : :
976 [ + - ]: 2 : if (tnbfac > 0)
977 : : {
978 : :
979 [ - + ][ - - ]: 2 : Assert( !inpoel.empty(),
[ - - ][ - - ]
980 : : "Attempt to call genNbfacTet() on empty inpoel container" );
981 [ - + ][ - - ]: 2 : Assert( !triinpoel_complete.empty(),
[ - - ][ - - ]
982 : : "Attempt to call genNbfacTet() on empty triinpoel_complete container" );
983 [ - + ][ - - ]: 2 : Assert( triinpoel_complete.size()/nnpf == tnbfac,
[ - - ][ - - ]
984 : : "Incorrect size of triinpoel in genNbfacTet()" );
985 : :
986 [ + - ]: 2 : auto nptet = inpoel;
987 [ + - ]: 2 : auto nptri = triinpoel_complete;
988 : :
989 [ + - ]: 2 : unique( nptet );
990 [ + - ]: 2 : unique( nptri );
991 : :
992 : 2 : std::unordered_set< std::size_t > snptet;
993 : :
994 : : // getting the reduced inpoel as a set for quick searches
995 [ + - ]: 2 : snptet.insert( begin(nptet), end(nptet));
996 : :
997 : : // vector to store boundary-face-nodes in this chunk
998 : 2 : std::vector< std::size_t > nptri_chunk;
999 : :
1000 : : // getting the nodes of the boundary-faces in this chunk
1001 [ + + ]: 30 : for (auto i : nptri)
1002 [ + - ][ + - ]: 28 : if (snptet.find(i) != end(snptet))
1003 [ + - ]: 28 : nptri_chunk.push_back(i);
1004 : :
1005 : : std::size_t tag, icoun;
1006 : :
1007 : : // matching nodes in nptri_chunk with nodes in inpoel and
1008 : : // triinpoel_complete to get the number of faces in this chunk
1009 [ + + ]: 4 : for (const auto& ss : bface_complete)
1010 : : {
1011 [ + + ]: 50 : for (auto f : ss.second)
1012 : : {
1013 : 48 : icoun = f*nnpf;
1014 : 48 : tag = 0;
1015 [ + + ]: 192 : for (std::size_t i=0; i<nnpf; ++i) {
1016 [ + + ]: 2160 : for (auto j : nptri_chunk) {
1017 : : // cppcheck-suppress useStlAlgorithm
1018 [ + + ]: 2016 : if (triinpoel_complete[icoun+i] == j) ++tag;
1019 : : }
1020 : : }
1021 [ + - ]: 48 : if (tag == nnpf)
1022 : : // this is a boundary face
1023 : : {
1024 [ + + ]: 192 : for (std::size_t i=0; i<nnpf; ++i)
1025 : : {
1026 : 144 : auto ip = triinpoel_complete[icoun+i];
1027 : :
1028 : : // find local renumbered node-id to store in triinpoel
1029 [ + - ][ + - ]: 144 : triinpoel.push_back( cref_find(lid,ip) );
1030 : : }
1031 : :
1032 [ + - ][ + - ]: 48 : bface[ss.first].push_back(nbfac);
1033 : 48 : ++nbfac;
1034 : : }
1035 : : }
1036 : : }
1037 : :
1038 : 2 : }
1039 : :
1040 : 2 : return nbfac;
1041 : : }
1042 : :
1043 : : std::vector< int >
1044 : 5215 : genEsuelTet( const std::vector< std::size_t >& inpoel,
1045 : : const std::pair< std::vector< std::size_t >,
1046 : : std::vector< std::size_t > >& esup )
1047 : : // *****************************************************************************
1048 : : // Generate derived data structure, elements surrounding elements
1049 : : // as a fixed length data structure as a full vector, including
1050 : : // boundary elements as -1.
1051 : : // \warning This is for Tetrahedra only.
1052 : : //! \param[in] inpoel Inteconnectivity of points and elements. These are the
1053 : : //! node ids of each element of an unstructured mesh. Example:
1054 : : //! \code{.cpp}
1055 : : //! std::vector< std::size_t > inpoel { 12, 14, 9, 11,
1056 : : //! 10, 14, 13, 12 };
1057 : : //! \endcode
1058 : : //! specifies two tetrahedra whose vertices (node ids) are { 12, 14, 9, 11 },
1059 : : //! and { 10, 14, 13, 12 }.
1060 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
1061 : : //! \return Vector storing elements surrounding elements
1062 : : //! \warning It is not okay to call this function with an empty container for
1063 : : //! inpoel or esup.first or esup.second; it will throw an exception.
1064 : : //! \details The data generated here is stored in a single vector, with length
1065 : : //! nfpe * nelem. Note however, that nelem is not explicitly provided, but
1066 : : //! calculated from inpoel. For boundary elements, at the boundary face, this
1067 : : //! esuelTet stores value -1 indicating that this is outside the domain. The
1068 : : //! convention for numbering the local face (triangle) connectivity is very
1069 : : //! important, e.g., in generating the inpofa array later. This node ordering
1070 : : //! convention is stored in tk::lpofa. Thus function is specific to
1071 : : //! tetrahedra, which is reflected in the fact that nnpe and nfpe are being
1072 : : //! set here in the function rather than being input arguments. To find out
1073 : : //! the number of elements, _nelem_, the size of the mesh connectivity vector,
1074 : : //! _inpoel_, can be devided by the number of nodes per elements, _nnpe_:
1075 : : //! \code{.cpp}
1076 : : //! auto nelem = inpoel.size()/nnpe;
1077 : : //! \endcode
1078 : : // *****************************************************************************
1079 : : {
1080 [ - + ][ - - ]: 5215 : Assert( !inpoel.empty(), "Attempt to call genEsuelTet() on empty container" );
[ - - ][ - - ]
1081 [ - + ][ - - ]: 5215 : Assert( !esup.first.empty(), "Attempt to call genEsuelTet() with empty esup1" );
[ - - ][ - - ]
1082 [ - + ][ - - ]: 5215 : Assert( !esup.second.empty(),
[ - - ][ - - ]
1083 : : "Attempt to call genEsuelTet() with empty esup2" );
1084 : :
1085 : 5215 : auto& esup1 = esup.first;
1086 : 5215 : auto& esup2 = esup.second;
1087 : :
1088 : : // set tetrahedron geometry
1089 : : // cppcheck-suppress unreadVariable
1090 : 5215 : std::size_t nnpe = 4, nfpe = 4, nnpf = 3;
1091 : :
1092 [ - + ][ - - ]: 5215 : Assert( inpoel.size()%nnpe == 0, "Size of inpoel must be divisible by four" );
[ - - ][ - - ]
1093 : :
1094 : : // get nelem and npoin
1095 : : // cppcheck-suppress unreadVariable
1096 : 5215 : auto nelem = inpoel.size()/nnpe;
1097 [ + - ]: 5215 : auto minmax = std::minmax_element( begin(inpoel), end(inpoel) );
1098 [ - + ][ - - ]: 5215 : Assert( *minmax.first == 0, "node ids should start from zero" );
[ - - ][ - - ]
1099 : 5215 : auto npoin = *minmax.second + 1;
1100 : :
1101 [ + - ]: 5215 : std::vector< int > esuelTet(nfpe*nelem, -1);
1102 [ + - ]: 5215 : std::vector< std::size_t > lhelp(nnpf,0),
1103 [ + - ]: 5215 : lpoin(npoin,0);
1104 : :
1105 [ + + ]: 2721138 : for (std::size_t e=0; e<nelem; ++e)
1106 : : {
1107 : 2715923 : auto mark = nnpe*e;
1108 [ + + ]: 13579615 : for (std::size_t fe=0; fe<nfpe; ++fe)
1109 : : {
1110 : : // array which stores points on this face
1111 : 10863692 : lhelp[0] = inpoel[mark+lpofa[fe][0]];
1112 : 10863692 : lhelp[1] = inpoel[mark+lpofa[fe][1]];
1113 : 10863692 : lhelp[2] = inpoel[mark+lpofa[fe][2]];
1114 : :
1115 : : // mark in this array
1116 : 10863692 : lpoin[lhelp[0]] = 1;
1117 : 10863692 : lpoin[lhelp[1]] = 1;
1118 : 10863692 : lpoin[lhelp[2]] = 1;
1119 : :
1120 : : // select a point on this face
1121 : 10863692 : auto ipoin = lhelp[0];
1122 : :
1123 : : // loop over elements around this point
1124 [ + + ]: 225538660 : for (std::size_t j=esup2[ipoin]+1; j<=esup2[ipoin+1]; ++j )
1125 : : {
1126 : 214674968 : auto jelem = esup1[j];
1127 : : // if this jelem is not e itself then proceed
1128 [ + + ]: 214674968 : if (jelem != e)
1129 : : {
1130 [ + + ]: 1019056380 : for (std::size_t fj=0; fj<nfpe; ++fj)
1131 : : {
1132 : 815245104 : std::size_t icoun(0);
1133 [ + + ]: 3260980416 : for (std::size_t jnofa=0; jnofa<nnpf; ++jnofa)
1134 : : {
1135 : 2445735312 : auto markj = jelem*nnpe;
1136 : 2445735312 : auto jpoin = inpoel[markj+lpofa[fj][jnofa]];
1137 [ + + ]: 2445735312 : if (lpoin[jpoin] == 1) { ++icoun; }
1138 : : }
1139 : : //store esuel if
1140 [ + + ]: 815245104 : if (icoun == nnpf)
1141 : : {
1142 : 9987708 : auto markf = nfpe*e;
1143 : 9987708 : esuelTet[markf+fe] = static_cast<int>(jelem);
1144 : :
1145 : 9987708 : markf = nfpe*jelem;
1146 : 9987708 : esuelTet[markf+fj] = static_cast<int>(e);
1147 : : }
1148 : : }
1149 : : }
1150 : : }
1151 : : // reset this array
1152 : 10863692 : lpoin[lhelp[0]] = 0;
1153 : 10863692 : lpoin[lhelp[1]] = 0;
1154 : 10863692 : lpoin[lhelp[2]] = 0;
1155 : : }
1156 : : }
1157 : :
1158 : 10430 : return esuelTet;
1159 : 5215 : }
1160 : :
1161 : : std::size_t
1162 : 8 : genNipfac( std::size_t nfpe,
1163 : : std::size_t nbfac,
1164 : : const std::vector< int >& esuelTet )
1165 : : // *****************************************************************************
1166 : : // Generate number of internal and physical-boundary faces
1167 : : //! \param[in] nfpe Number of faces per element.
1168 : : //! \param[in] nbfac Number of boundary faces.
1169 : : //! \param[in] esuelTet Elements surrounding elements.
1170 : : //! \return Total number of faces in the mesh
1171 : : //! \details The unsigned integer here gives the number of internal and
1172 : : //! physical-boundary faces in the mesh. The data structure does not include
1173 : : //! faces that are on partition/chare-boundaries.
1174 : : // *****************************************************************************
1175 : : {
1176 [ - + ][ - - ]: 8 : Assert( !esuelTet.empty(), "Attempt to call genNipfac() with empty esuelTet" );
[ - - ][ - - ]
1177 [ - + ][ - - ]: 8 : Assert( esuelTet.size()%nfpe == 0,
[ - - ][ - - ]
1178 : : "Size of esuelTet must be divisible by nfpe" );
1179 [ - + ][ - - ]: 8 : Assert( nfpe > 0, "Attempt to call genNipfac() with zero faces per element" );
[ - - ][ - - ]
1180 : :
1181 : 8 : auto nelem = esuelTet.size()/nfpe;
1182 : :
1183 : 8 : std::size_t nifac = 0;
1184 : :
1185 : : // loop through elements surrounding elements to find number of internal faces
1186 [ + + ]: 396 : for (std::size_t e=0; e<nelem; ++e)
1187 : : {
1188 [ + + ]: 1940 : for (std::size_t ip=nfpe*e; ip<nfpe*(e+1); ++ip)
1189 : : {
1190 [ + + ]: 1552 : if (esuelTet[ip] != -1)
1191 : : {
1192 [ + + ]: 1264 : if ( e<static_cast< std::size_t >(esuelTet[ip]) )
1193 : : {
1194 : 632 : ++nifac;
1195 : : }
1196 : : }
1197 : : }
1198 : : }
1199 : :
1200 : 8 : return nifac + nbfac;
1201 : : }
1202 : :
1203 : : std::vector< int >
1204 : 2 : genEsuf( std::size_t nfpe,
1205 : : std::size_t nipfac,
1206 : : std::size_t nbfac,
1207 : : const std::vector< std::size_t >& belem,
1208 : : const std::vector< int >& esuelTet )
1209 : : // *****************************************************************************
1210 : : // Generate derived data structure, elements surrounding faces
1211 : : //! \param[in] nfpe Number of faces per element.
1212 : : //! \param[in] nipfac Number of internal and physical-boundary faces.
1213 : : //! \param[in] nbfac Number of boundary faces.
1214 : : //! \param[in] belem Boundary element vector.
1215 : : //! \param[in] esuelTet Elements surrounding elements.
1216 : : //! \return Elements surrounding faces.
1217 : : //! \details The unsigned integer vector gives the IDs of the elements to the
1218 : : // left and the right of each face in the mesh. The convention followed
1219 : : // throughout is : The left element always has an ID smaller than the ID of
1220 : : // the right element.
1221 : : // *****************************************************************************
1222 : : {
1223 [ - + ][ - - ]: 2 : Assert( esuelTet.size()%nfpe == 0,
[ - - ][ - - ]
1224 : : "Size of esuelTet must be divisible by nfpe" );
1225 [ - + ][ - - ]: 2 : Assert( nfpe > 0, "Attempt to call genEsuf() with zero faces per element" );
[ - - ][ - - ]
1226 : :
1227 : 2 : auto nelem = esuelTet.size()/nfpe;
1228 : :
1229 [ + - ]: 2 : std::vector< int > esuf(2*nipfac);
1230 : :
1231 : : // counters for number of internal and boundary faces
1232 : 2 : std::size_t icoun(2*nbfac), bcoun(0);
1233 : :
1234 : : // loop to get face-element connectivity for internal faces
1235 [ + + ]: 148 : for (std::size_t e=0; e<nelem; ++e) {
1236 [ + + ]: 730 : for (std::size_t ip=nfpe*e; ip<nfpe*(e+1); ++ip) {
1237 : 584 : auto jelem = esuelTet[ip];
1238 [ + + ]: 584 : if (jelem != -1)
1239 : : {
1240 [ + + ]: 488 : if ( e < static_cast< std::size_t >(jelem) )
1241 : : {
1242 : 244 : esuf[icoun] = static_cast< int >(e);
1243 : 244 : esuf[icoun+1] = static_cast< int >(jelem);
1244 : 244 : icoun = icoun + 2;
1245 : : }
1246 : : }
1247 : : }
1248 : : }
1249 : :
1250 : : // loop to get face-element connectivity for physical-boundary faces
1251 : 2 : bcoun = 0;
1252 [ + + ]: 98 : for (auto ie : belem) {
1253 : 96 : esuf[bcoun] = static_cast< int >(ie);
1254 : 96 : esuf[bcoun+1] = -1; // outside domain
1255 : 96 : bcoun = bcoun + 2;
1256 : : }
1257 : :
1258 : 2 : return esuf;
1259 : : }
1260 : :
1261 : : std::vector< std::size_t >
1262 : 6 : genInpofaTet( std::size_t nipfac,
1263 : : std::size_t nbfac,
1264 : : const std::vector< std::size_t >& inpoel,
1265 : : const std::vector< std::size_t >& triinpoel,
1266 : : const std::vector< int >& esuelTet )
1267 : : // *****************************************************************************
1268 : : // Generate derived data structure, points on faces for tetrahedra only
1269 : : //! \param[in] nipfac Number of internal and physical-boundary faces.
1270 : : //! \param[in] nbfac Number of boundary faces.
1271 : : //! \param[in] inpoel Element-node connectivity.
1272 : : //! \param[in] triinpoel Face-node connectivity.
1273 : : //! \param[in] esuelTet Elements surrounding elements.
1274 : : //! \return Points surrounding faces. The unsigned integer vector gives the
1275 : : //! elements to the left and to the right of each face in the mesh.
1276 : : // *****************************************************************************
1277 : : {
1278 : 6 : std::vector< std::size_t > inpofa;
1279 : :
1280 : : // set tetrahedron geometry
1281 : : // cppcheck-suppress unreadVariable
1282 : 6 : std::size_t nnpe(4), nfpe(4), nnpf(3);
1283 : :
1284 [ - + ][ - - ]: 6 : Assert( esuelTet.size()%nfpe == 0,
[ - - ][ - - ]
1285 : : "Size of esuelTet must be divisible by nfpe" );
1286 [ - + ][ - - ]: 6 : Assert( inpoel.size()%nnpe == 0,
[ - - ][ - - ]
1287 : : "Size of inpoel must be divisible by nnpe" );
1288 : :
1289 [ + - ]: 6 : inpofa.resize(nnpf*nipfac);
1290 : :
1291 : : // counters for number of internal and boundary faces
1292 : 6 : std::size_t icoun(nnpf*nbfac);
1293 : :
1294 : : // loop over elems to get nodes on faces
1295 : : // this fills the interior face-node connectivity part
1296 [ + + ]: 248 : for (std::size_t e=0; e<inpoel.size()/nnpe; ++e)
1297 : : {
1298 : 242 : auto mark = nnpe*e;
1299 [ + + ]: 1210 : for (std::size_t f=0; f<nfpe ; ++f)
1300 : : {
1301 : 968 : auto ip = nfpe*e + f;
1302 : 968 : auto jelem = esuelTet[ip];
1303 [ + + ]: 968 : if (jelem != -1)
1304 : : {
1305 [ + + ]: 776 : if ( e < static_cast< std::size_t >(jelem) )
1306 : : {
1307 : 388 : inpofa[icoun] = inpoel[mark+lpofa[f][0]];
1308 : 388 : inpofa[icoun+1] = inpoel[mark+lpofa[f][1]];
1309 : 388 : inpofa[icoun+2] = inpoel[mark+lpofa[f][2]];
1310 : 388 : icoun = icoun + nnpf;
1311 : : }
1312 : : }
1313 : : }
1314 : : }
1315 : :
1316 : : // this fills the boundary face-node connectivity part
1317 : : // consistent with triinpoel
1318 [ + + ]: 198 : for (std::size_t f=0; f<nbfac; ++f)
1319 : : {
1320 : 192 : icoun = nnpf * f;
1321 : 192 : inpofa[icoun+0] = triinpoel[icoun+0];
1322 : 192 : inpofa[icoun+1] = triinpoel[icoun+1];
1323 : 192 : inpofa[icoun+2] = triinpoel[icoun+2];
1324 : : }
1325 : :
1326 : 6 : return inpofa;
1327 : 0 : }
1328 : :
1329 : : std::vector< std::size_t >
1330 : 4 : genBelemTet( std::size_t nbfac,
1331 : : const std::vector< std::size_t >& inpofa,
1332 : : const std::pair< std::vector< std::size_t >,
1333 : : std::vector< std::size_t > >& esup )
1334 : : // *****************************************************************************
1335 : : // Generate derived data, boundary elements
1336 : : //! \param[in] nbfac Number of boundary faces.
1337 : : //! \param[in] inpofa Face-node connectivity.
1338 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
1339 : : //! \return Host elements or boundary elements. The unsigned integer vector
1340 : : //! gives the elements to the left of each boundary face in the mesh.
1341 : : //! \details The data structure generated here contains an array of elements
1342 : : //! which share one or more of their faces with the physical boundary, i.e.,
1343 : : //! where exodus specifies a side-set for faces. Such elements are sometimes
1344 : : //! also called host or boundary elements.
1345 : : // *****************************************************************************
1346 : : {
1347 [ + - ]: 4 : std::vector< std::size_t > belem(nbfac);
1348 : :
1349 [ + - ]: 4 : if (nbfac > 0)
1350 : : {
1351 : :
1352 : : // set tetrahedron geometry
1353 : : // cppcheck-suppress unreadVariable
1354 : 4 : std::size_t nnpf = 3, tag = 0;
1355 : :
1356 : : // loop over all the boundary faces
1357 [ + + ]: 100 : for(std::size_t f=0; f<nbfac; ++f)
1358 : : {
1359 : 96 : belem[f] = 0;
1360 : :
1361 : : // array storing the element-cluster around face
1362 : 96 : std::vector< std::size_t > elemcluster;
1363 : :
1364 : : // loop over the nodes of this boundary face
1365 [ + + ]: 384 : for(std::size_t lp=0; lp<nnpf; ++lp)
1366 : : {
1367 : 288 : auto gp = inpofa[nnpf*f + lp];
1368 : :
1369 [ - + ][ - - ]: 288 : Assert( gp < esup.second.size(), "Indexing out of esup2" );
[ - - ][ - - ]
1370 : : // loop over elements surrounding this node
1371 [ + + ]: 2080 : for (auto i=esup.second[gp]+1; i<=esup.second[gp+1]; ++i)
1372 : : {
1373 : : // form element-cluster vector
1374 [ + - ]: 1792 : elemcluster.push_back(esup.first[i]);
1375 : : }
1376 : : }
1377 : :
1378 : : // loop over element cluster to find repeating elements
1379 [ + - ]: 272 : for(std::size_t i=0; i<elemcluster.size(); ++i)
1380 : : {
1381 : 272 : auto ge = elemcluster[i];
1382 : 272 : tag = 1;
1383 [ + + ]: 5384 : for(std::size_t j=0; j<elemcluster.size(); ++j)
1384 : : {
1385 [ + + ][ + + ]: 5112 : if ( i != j && elemcluster[j] == ge )
[ + + ]
1386 : : {
1387 : 248 : tag++;
1388 : : }
1389 : : }
1390 [ + + ]: 272 : if (tag == nnpf)
1391 : : {
1392 : : // this is the required boundary element
1393 : 96 : belem[f] = ge;
1394 : 96 : break;
1395 : : }
1396 : : }
1397 : 96 : }
1398 : : }
1399 : :
1400 : 4 : return belem;
1401 : 0 : }
1402 : :
1403 : : bool
1404 : 2601 : leakyPartition( const std::vector< int >& esueltet,
1405 : : const std::vector< std::size_t >& inpoel,
1406 : : const std::array< std::vector< real >, 3 >& coord )
1407 : : // *****************************************************************************
1408 : : // Perform leak-test on mesh (partition)
1409 : : //! \param[in] esueltet Elements surrounding elements for tetrahedra, see
1410 : : //! tk::genEsueltet()
1411 : : //! \param[in] inpoel Element connectivity
1412 : : //! \param[in] coord Node coordinates
1413 : : //! \details This function computes a surface integral over the boundary of the
1414 : : //! incoming mesh (partition). A non-zero vector result indicates a leak, e.g.,
1415 : : //! a hole in the mesh (partition), which indicates an error either in the
1416 : : // mesh geometry, mesh partitioning, or in the data structures that represent
1417 : : // faces.
1418 : : //! \return True if partition leaks.
1419 : : // *****************************************************************************
1420 : : {
1421 : 2601 : const auto& x = coord[0];
1422 : 2601 : const auto& y = coord[1];
1423 : 2601 : const auto& z = coord[2];
1424 : :
1425 : : // Storage for surface integral over our mesh partition
1426 : 2601 : std::array< real, 3 > s{{ 0.0, 0.0, 0.0}};
1427 : :
1428 [ + + ]: 1120652 : for (std::size_t e=0; e<esueltet.size()/4; ++e) { // for all our tets
1429 : 1118051 : auto mark = e*4;
1430 [ + + ]: 5590255 : for (std::size_t f=0; f<4; ++f) // for all tet faces
1431 [ + + ]: 4472204 : if (esueltet[mark+f] == -1) { // if face has no outside-neighbor tet
1432 : : // 3 local node IDs of face
1433 : 407252 : auto A = inpoel[ mark + lpofa[f][0] ];
1434 : 407252 : auto B = inpoel[ mark + lpofa[f][1] ];
1435 : 407252 : auto C = inpoel[ mark + lpofa[f][2] ];
1436 : : // Compute face area and normal
1437 : : real nx, ny, nz;
1438 : 407252 : auto a = normal( x[A],x[B],x[C], y[A],y[B],y[C], z[A],z[B],z[C],
1439 : : nx, ny, nz );
1440 : : // Sum up face area * face unit-normal
1441 : 407252 : s[0] += a * nx;
1442 : 407252 : s[1] += a * ny;
1443 : 407252 : s[2] += a * nz;
1444 : : }
1445 : : }
1446 : :
1447 : 2601 : auto eps = 1.0e-9;
1448 [ + - ][ + - ]: 2601 : return std::abs(s[0]) > eps || std::abs(s[1]) > eps || std::abs(s[2]) > eps;
[ - + ]
1449 : : }
1450 : :
1451 : : bool
1452 : 5007 : conforming( const std::vector< std::size_t >& inpoel,
1453 : : const std::array< std::vector< real >, 3 >& coord,
1454 : : bool cerr,
1455 : : const std::vector< std::size_t >& rid )
1456 : : // *****************************************************************************
1457 : : // Check if mesh (partition) is conforming
1458 : : //! \param[in] inpoel Element connectivity
1459 : : //! \param[in] coord Node coordinates
1460 : : //! \param[in] cerr True if hanging-node edge data should be output to
1461 : : //! std::cerr (true by default)
1462 : : //! \param[in] rid AMR Lib node id map
1463 : : //! std::cerr (true by default)
1464 : : //! \return True if mesh (partition) has no hanging nodes and thus the mesh is
1465 : : //! conforming, false if non-conforming.
1466 : : //! \details A conforming mesh by definition has no hanging nodes. A node is
1467 : : //! hanging if an edge of one element coincides with two (or more) edges (of
1468 : : //! two or more other elements). Thus, testing for conformity relies on
1469 : : //! checking the coordinates of all vertices: if any vertex coincides with
1470 : : //! that of a mid-point node of an edge, that is a hanging node. Note that
1471 : : //! this assumes that hanging nodes can only be at the mid-point of edges.
1472 : : //! This may happen after a mesh refinement step, due to a problem/bug,
1473 : : //! within the mesh refinement algorithm given by J. Waltz, Parallel adaptive
1474 : : //! refinement for unsteady flow calculations on 3D unstructured grids,
1475 : : //! International Journal for Numerical Methods in Fluids, 46: 37–57, 2004,
1476 : : //! which always adds/removes vertices at the mid-points of edges of a
1477 : : //! tetrahedron mesh within a single refinement step. Thus this algorithm is
1478 : : //! intended for this specific case, i.e., test for conformity after a
1479 : : //! single refinement step and not after multiple ones or for detecting
1480 : : //! hanging nodes in an arbitrary mesh.
1481 : : //*****************************************************************************
1482 : : {
1483 [ + + ][ + - ]: 5007 : Assert( !inpoel.empty(),
[ + - ][ + - ]
1484 : : "Attempt to call conforming() with empty mesh connectivity" );
1485 [ + + ][ + - ]: 5003 : Assert( inpoel.size() % 4 == 0,
[ + - ][ + - ]
1486 : : "Size of inpoel must be divisible by nnpe" );
1487 [ + - ][ + + ]: 5001 : Assert( *std::min_element( begin(inpoel), end(inpoel) ) == 0,
[ + - ][ + - ]
[ + - ]
1488 : : "Node ids should start from zero" );
1489 [ + - ][ + - ]: 4999 : Assert( !coord[0].empty() && !coord[1].empty() && !coord[2].empty(),
[ + - ][ - - ]
[ - - ][ - - ]
1490 : : "Attempt to call conforming() with empty coordinates container" );
1491 : :
1492 : : using Coord = UnsMesh::Coord;
1493 : : using Edge = UnsMesh::Edge;
1494 : : using Tet = UnsMesh::Tet;
1495 : :
1496 : : // Compare operator to be used as less-than for std::array< real, 3 >,
1497 : : // implemented as a lexicographic ordering.
1498 : : struct CoordLess {
1499 : : const real eps = std::numeric_limits< real >::epsilon();
1500 : 178128395 : bool operator() ( const Coord& lhs, const Coord& rhs ) const {
1501 [ + + ]: 178128395 : if (lhs[0] < rhs[0])
1502 : 87355407 : return true;
1503 [ + + ][ + + ]: 90772988 : else if (std::abs(lhs[0]-rhs[0]) < eps && lhs[1] < rhs[1])
[ + + ]
1504 : 4317931 : return true;
1505 : 86455057 : else if (std::abs(lhs[0]-rhs[0]) < eps &&
1506 [ + + ][ + + ]: 105885598 : std::abs(lhs[1]-rhs[1]) < eps &&
[ + + ]
1507 [ + + ]: 19430541 : lhs[2] < rhs[2])
1508 : 833701 : return true;
1509 : : else
1510 : 85621356 : return false;
1511 : : }
1512 : : };
1513 : :
1514 : : // Map associating data on potential hanging nodes. Key: coordinates of nodes
1515 : : // of edge-half points, value: tet id (local if in parallel), tet connectivity
1516 : : // (using local ids if in parallel), edge ids (local if in parallel).
1517 : : std::map< Coord, // edge-half node coordinates: x, y, z
1518 : : std::tuple< std::size_t, // element id of edge-half node
1519 : : Tet, // element node ids of edge-half node
1520 : : Edge >, // edge containing half-node
1521 : 4999 : CoordLess > edgeNodes;
1522 : :
1523 : 4999 : const auto& x = coord[0];
1524 : 4999 : const auto& y = coord[1];
1525 : 4999 : const auto& z = coord[2];
1526 : :
1527 : : fenv_t fe;
1528 : 4999 : feholdexcept( &fe );
1529 : :
1530 : : // Compute coordinates of nodes of mid-points of all edges
1531 [ + + ]: 2022933 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1532 : 2017934 : auto A = inpoel[e*4+0];
1533 : 2017934 : auto B = inpoel[e*4+1];
1534 : 2017934 : auto C = inpoel[e*4+2];
1535 : 2017934 : auto D = inpoel[e*4+3];
1536 : : std::array<Edge,6> edge{{ {{A,B}}, {{B,C}}, {{A,C}},
1537 : 2017934 : {{A,D}}, {{B,D}}, {{C,D}} }};
1538 [ + + ]: 14125538 : for (const auto& n : edge) {
1539 : 12107604 : Coord en{{ (x[n[0]] + x[n[1]]) / 2.0,
1540 : 12107604 : (y[n[0]] + y[n[1]]) / 2.0,
1541 : 24215208 : (z[n[0]] + z[n[1]]) / 2.0 }};
1542 [ + - ]: 12107604 : edgeNodes[ en ] = std::tuple<std::size_t,Tet,Edge>{ e, {{A,B,C,D}}, n };
1543 : : }
1544 : : }
1545 : :
1546 : 4999 : feclearexcept( FE_UNDERFLOW );
1547 : 4999 : feupdateenv( &fe );
1548 : :
1549 : : // Find hanging nodes. If the coordinates of an element vertex coincide with
1550 : : // that of a mid-point node of an edge, that is a hanging node. If we find one
1551 : : // such node we print out some info on it.
1552 : 4999 : auto ix = x.cbegin();
1553 : 4999 : auto iy = y.cbegin();
1554 : 4999 : auto iz = z.cbegin();
1555 : :
1556 : 4999 : bool hanging_node = false;
1557 : :
1558 [ + + ]: 575921 : while (ix != x.cend()) {
1559 : 570922 : Coord n{{ *ix, *iy, *iz }};
1560 [ + - ]: 570922 : auto i = edgeNodes.find( n );
1561 [ + + ]: 570922 : if (i != end(edgeNodes)) {
1562 : 2 : const auto& hanging_node_coord = i->first;
1563 : 2 : const auto& hanging_node_info = i->second;
1564 : 2 : auto tet_id = std::get< 0 >( hanging_node_info );
1565 : 2 : const auto& tet = std::get< 1 >( hanging_node_info );
1566 : 2 : const auto& edge = std::get< 2 >( hanging_node_info );
1567 [ - + ]: 2 : if (cerr) {
1568 : : std::cerr
1569 [ - - ]: 0 : << "Mesh conformity test found hanging node with coordinates"" ("
1570 [ - - ][ - - ]: 0 : << hanging_node_coord[0] << ", "
1571 [ - - ][ - - ]: 0 : << hanging_node_coord[1] << ", "
1572 [ - - ][ - - ]: 0 : << hanging_node_coord[2] << ") of tetrahedron element "
1573 [ - - ][ - - ]: 0 : << tet_id << " with connectivity (" << tet[0] << ','
[ - - ][ - - ]
1574 [ - - ][ - - ]: 0 : << tet[1] << ',' << tet[2] << ',' << tet[3] << ") on edge ("
[ - - ][ - - ]
[ - - ][ - - ]
1575 [ - - ][ - - ]: 0 : << edge[0] << ',' << edge[1] << ")"
[ - - ]
1576 [ - - ][ - - ]: 0 : << "AMR lib node ids for this edge: " << rid[edge[0]] << ','
[ - - ][ - - ]
1577 [ - - ][ - - ]: 0 : << rid[edge[1]] << std::endl;
1578 : : }
1579 : 2 : hanging_node = true;
1580 : : }
1581 : 570922 : ++ix; ++iy; ++iz;
1582 : : }
1583 : :
1584 [ + + ]: 4999 : if (hanging_node) return false;
1585 : :
1586 : 4997 : return true;
1587 : 4999 : }
1588 : :
1589 : : bool
1590 : 117636 : intet( const std::array< std::vector< real >, 3 >& coord,
1591 : : const std::vector< std::size_t >& inpoel,
1592 : : const std::vector< real >& p,
1593 : : std::size_t e,
1594 : : std::array< real, 4 >& N )
1595 : : // *****************************************************************************
1596 : : // Determine if a point is in a tetrahedron
1597 : : //! \param[in] coord Mesh node coordinates
1598 : : //! \param[in] inpoel Mesh element connectivity
1599 : : //! \param[in] p Point coordinates
1600 : : //! \param[in] e Mesh cell index
1601 : : //! \param[in,out] N Shapefunctions evaluated at the point
1602 : : //! \return True if ppoint is in mesh cell
1603 : : //! \see Lohner, An Introduction to Applied CFD Techniques, Wiley, 2008
1604 : : // *****************************************************************************
1605 : : {
1606 [ - + ][ - - ]: 117636 : Assert( p.size() == 3, "Size mismatch" );
[ - - ][ - - ]
1607 : :
1608 : : // Tetrahedron node indices
1609 : 117636 : const auto A = inpoel[e*4+0];
1610 : 117636 : const auto B = inpoel[e*4+1];
1611 : 117636 : const auto C = inpoel[e*4+2];
1612 : 117636 : const auto D = inpoel[e*4+3];
1613 : :
1614 : : // Tetrahedron node coordinates
1615 : 117636 : const auto& x = coord[0];
1616 : 117636 : const auto& y = coord[1];
1617 : 117636 : const auto& z = coord[2];
1618 : :
1619 : : // Point coordinates
1620 : 117636 : const auto& xp = p[0];
1621 : 117636 : const auto& yp = p[1];
1622 : 117636 : const auto& zp = p[2];
1623 : :
1624 : : // Evaluate linear shapefunctions at point locations using Cramer's Rule
1625 : : // | xp | | x1 x2 x3 x4 | | N1 |
1626 : : // | yp | = | y1 y2 y3 y4 | • | N2 |
1627 : : // | zp | | z1 z2 z3 z4 | | N3 |
1628 : : // | 1 | | 1 1 1 1 | | N4 |
1629 : :
1630 : 117636 : real DetX = (y[B]*z[C] - y[C]*z[B] - y[B]*z[D] + y[D]*z[B] +
1631 : 117636 : y[C]*z[D] - y[D]*z[C])*x[A] + x[B]*y[C]*z[A] - x[B]*y[A]*z[C] +
1632 : 117636 : x[C]*y[A]*z[B] - x[C]*y[B]*z[A] + x[B]*y[A]*z[D] - x[B]*y[D]*z[A] -
1633 : 117636 : x[D]*y[A]*z[B] + x[D]*y[B]*z[A] - x[C]*y[A]*z[D] + x[C]*y[D]*z[A] +
1634 : 117636 : x[D]*y[A]*z[C] - x[D]*y[C]*z[A] - x[B]*y[C]*z[D] + x[B]*y[D]*z[C] +
1635 : 117636 : x[C]*y[B]*z[D] - x[C]*y[D]*z[B] - x[D]*y[B]*z[C] + x[D]*y[C]*z[B];
1636 : :
1637 : 117636 : real DetX1 = (y[D]*z[C] - y[C]*z[D] + y[C]*zp - yp*z[C] -
1638 : 117636 : y[D]*zp + yp*z[D])*x[B] + x[C]*y[B]*z[D] - x[C]*y[D]*z[B] -
1639 : 117636 : x[D]*y[B]*z[C] + x[D]*y[C]*z[B] - x[C]*y[B]*zp + x[C]*yp*z[B] +
1640 : 117636 : xp*y[B]*z[C] - xp*y[C]*z[B] + x[D]*y[B]*zp - x[D]*yp*z[B] -
1641 : 117636 : xp*y[B]*z[D] + xp*y[D]*z[B] + x[C]*y[D]*zp - x[C]*yp*z[D] -
1642 : 117636 : x[D]*y[C]*zp + x[D]*yp*z[C] + xp*y[C]*z[D] - xp*y[D]*z[C];
1643 : :
1644 : 117636 : real DetX2 = (y[C]*z[D] - y[D]*z[C] - y[C]*zp + yp*z[C] +
1645 : 117636 : y[D]*zp - yp*z[D])*x[A] + x[C]*y[D]*z[A] - x[C]*y[A]*z[D] +
1646 : 117636 : x[D]*y[A]*z[C] - x[D]*y[C]*z[A] + x[C]*y[A]*zp - x[C]*yp*z[A] -
1647 : 117636 : xp*y[A]*z[C] + xp*y[C]*z[A] - x[D]*y[A]*zp + x[D]*yp*z[A] +
1648 : 117636 : xp*y[A]*z[D] - xp*y[D]*z[A] - x[C]*y[D]*zp + x[C]*yp*z[D] +
1649 : 117636 : x[D]*y[C]*zp - x[D]*yp*z[C] - xp*y[C]*z[D] + xp*y[D]*z[C];
1650 : :
1651 : 117636 : real DetX3 = (y[D]*z[B] - y[B]*z[D] + y[B]*zp - yp*z[B] -
1652 : 117636 : y[D]*zp + yp*z[D])*x[A] + x[B]*y[A]*z[D] - x[B]*y[D]*z[A] -
1653 : 117636 : x[D]*y[A]*z[B] + x[D]*y[B]*z[A] - x[B]*y[A]*zp + x[B]*yp*z[A] +
1654 : 117636 : xp*y[A]*z[B] - xp*y[B]*z[A] + x[D]*y[A]*zp - x[D]*yp*z[A] -
1655 : 117636 : xp*y[A]*z[D] + xp*y[D]*z[A] + x[B]*y[D]*zp - x[B]*yp*z[D] -
1656 : 117636 : x[D]*y[B]*zp + x[D]*yp*z[B] + xp*y[B]*z[D] - xp*y[D]*z[B];
1657 : :
1658 : 117636 : real DetX4 = (y[B]*z[C] - y[C]*z[B] - y[B]*zp + yp*z[B] +
1659 : 117636 : y[C]*zp - yp*z[C])*x[A] + x[B]*y[C]*z[A] - x[B]*y[A]*z[C] +
1660 : 117636 : x[C]*y[A]*z[B] - x[C]*y[B]*z[A] + x[B]*y[A]*zp - x[B]*yp*z[A] -
1661 : 117636 : xp*y[A]*z[B] + xp*y[B]*z[A] - x[C]*y[A]*zp + x[C]*yp*z[A] +
1662 : 117636 : xp*y[A]*z[C] - xp*y[C]*z[A] - x[B]*y[C]*zp + x[B]*yp*z[C] +
1663 : 117636 : x[C]*y[B]*zp - x[C]*yp*z[B] - xp*y[B]*z[C] + xp*y[C]*z[B];
1664 : :
1665 : : // Shape functions evaluated at point
1666 : 117636 : N[0] = DetX1/DetX;
1667 : 117636 : N[1] = DetX2/DetX;
1668 : 117636 : N[2] = DetX3/DetX;
1669 : 117636 : N[3] = DetX4/DetX;
1670 : :
1671 : : // if min( N^i, 1-N^i ) > 0 for all i, point is in cell
1672 [ + + ]: 127788 : if ( std::min(N[0],1.0-N[0]) > 0 && std::min(N[1],1.0-N[1]) > 0 &&
1673 [ + + ][ + + ]: 127788 : std::min(N[2],1.0-N[2]) > 0 && std::min(N[3],1.0-N[3]) > 0 )
[ + + ][ + + ]
1674 : : {
1675 : 75 : return true;
1676 : : } else {
1677 : 117561 : return false;
1678 : : }
1679 : : }
1680 : :
1681 : : } // tk::
|
