Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Mesh/DerivedData.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief 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 : : #ifndef DerivedData_h
14 : : #define DerivedData_h
15 : :
16 : : #include <vector>
17 : : #include <map>
18 : : #include <utility>
19 : : #include <cstddef>
20 : : #include "Types.hpp"
21 : : #include "Fields.hpp"
22 : : #include "UnsMesh.hpp"
23 : :
24 : : namespace tk {
25 : :
26 : : //! Const array defining the node ordering convention for tetrahedron faces
27 : : //! \details This two-dimensional array stores the naming/ordering convention of
28 : : //! the node indices of a tetrahedron (tet) element. The dimensions are 4x3 as
29 : : //! a tetrahedron has a total of 4 nodes and each (triangle) face has 3 nodes.
30 : : //! Thus the array below associates tet node 0 with nodes {1,2,3}, tet node 1
31 : : //! with {2,0,3}, tet node 2 with {3,0,1}, and tet node 3 with {0,2,1}. Note
32 : : //! that not only these mappings are important, but also the order of the
33 : : //! nodes within the triplets as this specific order also defines the outwards
34 : : //! normal of each face.
35 : : const std::array< UnsMesh::Face, 4 >
36 : : lpofa{{ {{1,2,3}}, {{2,0,3}}, {{3,0,1}}, {{0,2,1}} }};
37 : :
38 : : //! Const array defining the node ordering convention for tetrahedron edges
39 : : const std::array< UnsMesh::Edge, 6 >
40 : : lpoed{{ {{0,1}}, {{1,2}}, {{0,2}}, {{0,3}}, {{1,3}}, {{2,3}} }};
41 : :
42 : : //! Const array defining the node ordering convention for triangle edges
43 : : const std::array< UnsMesh::Edge, 3 > lpoet{{ {{0,1}}, {{1,2}}, {{2,0}} }};
44 : :
45 : : //! Determine edge orientation
46 : : int
47 : : orient( const UnsMesh::Edge& t, const UnsMesh::Edge& e );
48 : :
49 : : //! Compute number of points (nodes) in mesh from connectivity
50 : : std::size_t
51 : : npoin_in_graph( const std::vector< std::size_t >& inpoel );
52 : :
53 : : //! Compute the unit normal vector of a triangle
54 : : //! \param[in] x1 x coordinate of the 1st vertex of the triangle
55 : : //! \param[in] x2 x coordinate of the 2nd vertex of the triangle
56 : : //! \param[in] x3 x coordinate of the 3rd vertex of the triangle
57 : : //! \param[in] y1 y coordinate of the 1st vertex of the triangle
58 : : //! \param[in] y2 y coordinate of the 2nd vertex of the triangle
59 : : //! \param[in] y3 y coordinate of the 3rd vertex of the triangle
60 : : //! \param[in] z1 z coordinate of the 1st vertex of the triangle
61 : : //! \param[in] z2 z coordinate of the 2nd vertex of the triangle
62 : : //! \param[in] z3 z coordinate of the 3rd vertex of the triangle
63 : : //! \param[out] nx x coordinate of the unit normal
64 : : //! \param[out] ny y coordinate of the unit normal
65 : : //! \param[out] nz z coordinate of the unit normal
66 : : #pragma omp declare simd
67 : : inline void
68 : 6327005 : normal( real x1, real x2, real x3,
69 : : real y1, real y2, real y3,
70 : : real z1, real z2, real z3,
71 : : real& nx, real& ny, real& nz )
72 : : {
73 : 6327005 : real ax = x2 - x1;
74 : 6327005 : real ay = y2 - y1;
75 : 6327005 : real az = z2 - z1;
76 : :
77 : 6327005 : real bx = x3 - x1;
78 : 6327005 : real by = y3 - y1;
79 : 6327005 : real bz = z3 - z1;
80 : :
81 : 6327005 : real n1 = ay*bz - az*by;
82 : 6327005 : real n2 = -(ax*bz - az*bx);
83 : 6327005 : real n3 = ax*by - ay*bx;
84 : :
85 : 6327005 : auto farea = std::sqrt( n1*n1 + n2*n2 + n3*n3 );
86 : :
87 : 6327005 : nx = n1/farea;
88 : 6327005 : ny = n2/farea;
89 : 6327005 : nz = n3/farea;
90 : 6327005 : }
91 : :
92 : : //! Compute the unit normal vector of a triangle
93 : : std::array< real, 3 >
94 : : normal( const std::array< real, 3 >& x,
95 : : const std::array< real, 3 >& y,
96 : : const std::array< real, 3 >& z );
97 : :
98 : : //! Compute the are of a triangle
99 : : //! \param[in] x1 x coordinate of the 1st vertex of the triangle
100 : : //! \param[in] x2 x coordinate of the 2nd vertex of the triangle
101 : : //! \param[in] x3 x coordinate of the 3rd vertex of the triangle
102 : : //! \param[in] y1 y coordinate of the 1st vertex of the triangle
103 : : //! \param[in] y2 y coordinate of the 2nd vertex of the triangle
104 : : //! \param[in] y3 y coordinate of the 3rd vertex of the triangle
105 : : //! \param[in] z1 z coordinate of the 1st vertex of the triangle
106 : : //! \param[in] z2 z coordinate of the 2nd vertex of the triangle
107 : : //! \param[in] z3 z coordinate of the 3rd vertex of the triangle
108 : : //! \return Area of the triangle
109 : : #pragma omp declare simd
110 : : inline real
111 : 6327005 : area( real x1, real x2, real x3,
112 : : real y1, real y2, real y3,
113 : : real z1, real z2, real z3 )
114 : : {
115 : 6327005 : auto sidea = sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1) );
116 : 6327005 : auto sideb = sqrt( (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) + (z3-z2)*(z3-z2) );
117 : 6327005 : auto sidec = sqrt( (x1-x3)*(x1-x3) + (y1-y3)*(y1-y3) + (z1-z3)*(z1-z3) );
118 : :
119 : 6327005 : auto semip = 0.5 * (sidea + sideb + sidec);
120 : :
121 : 6327005 : return sqrt( semip * (semip-sidea) * (semip-sideb) * (semip-sidec) );
122 : : }
123 : :
124 : : // Compute local face id wrt to tet based on tet and face connectivities
125 : 764086 : inline int opposite_vertex_of_tet(
126 : : const std::array< std::size_t,4 >& el,
127 : : const std::array< std::size_t,3 >& face) noexcept
128 : : {
129 : : // Bitmask trick
130 : : // mark the three face nodes
131 : 764086 : uint8_t mask = 0;
132 [ + + ]: 3056344 : for (auto fn : face) {
133 [ + - ][ + + ]: 5734538 : for (int i=0; i<4; ++i) if (el[static_cast<std::size_t>(i)]==fn)
134 : 2292258 : { mask |= (1u<<i); break; }
135 : : }
136 : : // the missing bit is the opposite vertex
137 [ + - ][ + + ]: 1906322 : for (int i=0; i<4; ++i) if ((mask & (1u<<i)) == 0) return i;
138 : 0 : return 0; // unreachable if topology is valid
139 : : }
140 : :
141 : : //! Compute the area of a triangle
142 : : real
143 : : area( const std::array< real, 3 >& x,
144 : : const std::array< real, 3 >& y,
145 : : const std::array< real, 3 >& z );
146 : :
147 : : //! Generate derived data structure, elements surrounding points
148 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
149 : : genEsup( const std::vector< std::size_t >& inpoel, std::size_t nnpe );
150 : :
151 : : //! Generate derived data structure, points surrounding points
152 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
153 : : genPsup( const std::vector< std::size_t >& inpoel,
154 : : std::size_t nnpe,
155 : : const std::pair< std::vector< std::size_t >,
156 : : std::vector< std::size_t > >& esup );
157 : :
158 : : //! Generate derived data structure, edges surrounding points
159 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
160 : : genEdsup( const std::vector< std::size_t >& inpoel,
161 : : std::size_t nnpe,
162 : : const std::pair< std::vector< std::size_t >,
163 : : std::vector< std::size_t > >& esup );
164 : :
165 : : //! Generate derived data structure, edge connectivity
166 : : std::vector< std::size_t >
167 : : genInpoed( const std::vector< std::size_t >& inpoel,
168 : : std::size_t nnpe,
169 : : const std::pair< std::vector< std::size_t >,
170 : : std::vector< std::size_t > >& esup );
171 : :
172 : : //! Generate derived data structure, elements surrounding points of elements
173 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
174 : : genEsupel( const std::vector< std::size_t >& inpoel,
175 : : std::size_t nnpe,
176 : : const std::pair< std::vector< std::size_t >,
177 : : std::vector< std::size_t > >& esup );
178 : :
179 : : //! Generate derived data structure, elements surrounding elements
180 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > >
181 : : genEsuel( const std::vector< std::size_t >& inpoel,
182 : : std::size_t nnpe,
183 : : const std::pair< std::vector< std::size_t >,
184 : : std::vector< std::size_t > >& esup );
185 : :
186 : : //! \brief Generate derived data structure, elements surrounding elements
187 : : //! as a fixed length data structure as a full vector, including boundary
188 : : //! elements as -1.
189 : : std::vector< int >
190 : : genEsuelTet( const std::vector< std::size_t >& inpoel,
191 : : const std::pair< std::vector< std::size_t >,
192 : : std::vector< std::size_t > >& esup );
193 : :
194 : : //! Generate derived data structure, edges of elements
195 : : std::vector< std::size_t >
196 : : genInedel( const std::vector< std::size_t >& inpoel,
197 : : std::size_t nnpe,
198 : : const std::vector< std::size_t >& inpoed );
199 : :
200 : : //! Generate derived data structure, elements surrounding edges
201 : : std::unordered_map< UnsMesh::Edge, std::vector< std::size_t >,
202 : : UnsMesh::Hash<2>, UnsMesh::Eq<2> >
203 : : genEsued( const std::vector< std::size_t >& inpoel,
204 : : std::size_t nnpe,
205 : : const std::pair< std::vector< std::size_t >,
206 : : std::vector< std::size_t > >& esup );
207 : :
208 : : //! Generate total number of boundary faces in this chunk
209 : : std::size_t
210 : : genNbfacTet( std::size_t tnbfac,
211 : : const std::vector< std::size_t >& inpoel,
212 : : const std::vector< std::size_t >& triinpoel_complete,
213 : : const std::map< int, std::vector< std::size_t > >& bface_complete,
214 : : const std::unordered_map< std::size_t, std::size_t >& lid,
215 : : std::vector< std::size_t >& triinpoel,
216 : : std::map< int, std::vector< std::size_t > >& bface );
217 : :
218 : : //! Generate number of internal and physical-boundary faces
219 : : std::size_t
220 : : genNipfac( std::size_t nfpe,
221 : : std::size_t nbfac,
222 : : const std::vector< int >& esuelTet );
223 : :
224 : : //! Generate derived data structure, elements surrounding faces
225 : : std::vector< int >
226 : : genEsuf( std::size_t nfpe,
227 : : std::size_t nipfac,
228 : : std::size_t nbfac,
229 : : const std::vector< std::size_t >& belem,
230 : : const std::vector< int >& esuelTet );
231 : :
232 : : //! Generate derived data structure, node-face connectivity
233 : : std::vector< std::size_t >
234 : : genInpofaTet( std::size_t nipfac,
235 : : std::size_t nbfac,
236 : : const std::vector< std::size_t >& inpoel,
237 : : const std::vector< std::size_t >& triinpoel,
238 : : const std::vector< int >& esuelTet );
239 : :
240 : : //! Generate derived data structure, host/boundary element
241 : : std::vector< std::size_t >
242 : : genBelemTet( std::size_t nbfac,
243 : : const std::vector< std::size_t >& inpofa,
244 : : const std::pair< std::vector< std::size_t >,
245 : : std::vector< std::size_t > >& esup );
246 : :
247 : : //! Generate derived data structure, face geometry
248 : : Fields
249 : : genGeoFaceTri( std::size_t nipfac,
250 : : const std::vector< std::size_t >& inpofa,
251 : : const UnsMesh::Coords& coord );
252 : :
253 : : //! Compute geometry of the face given by three vertices
254 : : Fields
255 : : geoFaceTri( const std::array< real, 3 >& x,
256 : : const std::array< real, 3 >& y,
257 : : const std::array< real, 3 >& z );
258 : :
259 : : //! Generate derived data structure, element geometry
260 : : Fields
261 : : genGeoElemTet( const std::vector< std::size_t >& inpoel,
262 : : const UnsMesh::Coords& coord );
263 : :
264 : : //! Perform leak-test on mesh (partition)
265 : : bool
266 : : leakyPartition( const std::vector< int >& esueltet,
267 : : const std::vector< std::size_t >& inpoel,
268 : : const UnsMesh::Coords& coord );
269 : :
270 : : //! Check if mesh (partition) is conforming
271 : : bool
272 : : conforming( const std::vector< std::size_t >& inpoel,
273 : : const UnsMesh::Coords& coord,
274 : : bool cerr = true,
275 : : const std::vector< std::size_t >& rid={} );
276 : :
277 : : //! Determine if a point is in a tetrahedron
278 : : bool
279 : : intet( const std::array< std::vector< real >, 3 >& coord,
280 : : const std::vector< std::size_t >& inpoel,
281 : : const std::vector< real >& p,
282 : : std::size_t e,
283 : : std::array< real, 4 >& N );
284 : :
285 : : //! Compute curl of a vector field at nodes of unstructured tetrahedra mesh
286 : : tk::UnsMesh::Coords
287 : : curl( const std::array< std::vector< tk::real >, 3 >& coord,
288 : : const std::vector< std::size_t >& inpoel,
289 : : const tk::UnsMesh::Coords& v );
290 : :
291 : : //! Compute divergence of vector field at nodes of unstructured tetrahedra mesh
292 : : std::vector< tk::real >
293 : : div( const std::array< std::vector< tk::real >, 3 >& coord,
294 : : const std::vector< std::size_t >& inpoel,
295 : : const tk::UnsMesh::Coords& v );
296 : :
297 : : //! Compute gradient of a scalar field at nodes of unstructured tetrahedra mesh
298 : : tk::UnsMesh::Coords
299 : : grad( const std::array< std::vector< tk::real >, 3 >& coord,
300 : : const std::vector< std::size_t >& inpoel,
301 : : const std::vector< tk::real >& phi );
302 : :
303 : : } // tk::
304 : :
305 : : #endif // DerivedData_h
|
