Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Transfer/NodeSearch.cpp
4 : : \copyright 2020-2021 Charmworks, Inc.,
5 : : 2022-2025 J. Bakosi
6 : : All rights reserved. See the LICENSE file for details.
7 : : \brief Definitions pertaining to node search between 3d meshes
8 : : */
9 : : // *****************************************************************************
10 : :
11 : : #include <cmath>
12 : :
13 : : #include "NodeSearch.hpp"
14 : : #include "Transfer.hpp"
15 : : #include "DerivedData.hpp"
16 : : #include "InciterConfig.hpp"
17 : :
18 : : #include "NoWarning/collidecharm.h"
19 : :
20 : : #if defined(__clang__)
21 : : #pragma clang diagnostic push
22 : : #pragma clang diagnostic ignored "-Wold-style-cast"
23 : : #endif
24 : :
25 : 0 : PUPbytes( Collision );
26 : :
27 : : #if defined(__clang__)
28 : : #pragma clang diagnostic pop
29 : : #endif
30 : :
31 : : #include "NoWarning/transfer.decl.h"
32 : :
33 : : namespace transfer {
34 : : extern CollideHandle g_collideHandle;
35 : : extern CProxy_Transfer g_transferProxy;
36 : : }
37 : :
38 : : namespace inciter {
39 : : extern ctr::Config g_cfg;
40 : : }
41 : :
42 : : using transfer::NodeSearch;
43 : :
44 : 0 : NodeSearch::NodeSearch( CkArrayID p, MeshData mesh, CkCallback cb )
45 [ - - ]: 0 : : m_firstchunk( mesh.firstchunk )
46 : : // *****************************************************************************
47 : : // Constructor
48 : : //! \param[in] firstchunk Chunk ID used for the collision detection library
49 : : //! \param[in] cb Callback to inform application that the library is ready
50 : : // *****************************************************************************
51 : : {
52 : : //std::cout << "NodeSearch: " << mesh.nchare << '\n';
53 : : mesh.proxy = thisProxy;
54 [ - - ]: 0 : CollideRegister( g_collideHandle, /*ignored*/0 );
55 [ - - ]: 0 : g_transferProxy.ckLocalBranch()->setMesh( p, mesh );
56 [ - - ]: 0 : contribute( cb );
57 : 0 : }
58 : :
59 : : void
60 : 0 : NodeSearch::setSourceTets( const std::vector< std::size_t >& inpoel,
61 : : const std::array< std::vector< double >, 3 >& coord,
62 : : const tk::Fields& u,
63 : : const std::vector< double >& flag,
64 : : bool dir,
65 : : CkCallback cb )
66 : : // *****************************************************************************
67 : : // Set the data for the source tetrahedrons to be collided
68 : : //! \param[in] inpoel Pointer to the connectivity data for the source mesh
69 : : //! \param[in] coord Pointer to the coordinate data for the source mesh
70 : : //! \param[in] u Pointer to the solution data for the source mesh
71 : : //! \param[in] flag Transfer flags
72 : : //! \param[in] dir Transfer direction: 0: bg to overset, 1: overset to bg
73 : : //! \param[in] cb Callback to call when src side of transfer is done
74 : : // *****************************************************************************
75 : : {
76 : 0 : m_inpoel = const_cast< std::vector< std::size_t >* >( &inpoel );
77 : 0 : m_coord = const_cast< std::array< std::vector< double >, 3 >* >( &coord );
78 : 0 : m_u = const_cast< tk::Fields* >( &u );
79 : 0 : m_flag = const_cast< std::vector< double >* >( &flag );
80 : 0 : m_dir = dir;
81 : 0 : m_done = cb;
82 : 0 : m_srcnotified = 0;
83 : :
84 : : // Send tetrahedron data to the collision detection library
85 : 0 : collideTets();
86 : 0 : }
87 : :
88 : : void
89 : 0 : NodeSearch::setDestPoints( const std::array< std::vector< double >, 3 >& coord,
90 : : tk::Fields& u,
91 : : std::vector< double >& flag,
92 : : bool trflag,
93 : : bool dir,
94 : : CkCallback cb )
95 : : // *****************************************************************************
96 : : // Set the data for the destination points to be collided
97 : : //! \param[in] coord Pointer to the coordinate data for the destination mesh
98 : : //! \param[in,out] u Pointer to the solution data for the destination mesh
99 : : //! \param[in,out] flag Transfer flags
100 : : //! \param[in] trflag Transfer flags if true
101 : : //! \param[in] dir Transfer direction: 0: bg to overset, 1: overset to bg
102 : : //! \param[in] cb Callback to call when dst side of transfer is done
103 : : // *****************************************************************************
104 : : {
105 : 0 : m_coord = const_cast< std::array< std::vector< double >, 3 >* >( &coord );
106 : 0 : m_u = static_cast< tk::Fields* >( &u );
107 : 0 : m_flag = const_cast< std::vector< double >* >( &flag );
108 : 0 : m_trflag = trflag;
109 : 0 : m_dir = dir;
110 : 0 : m_done = cb;
111 : :
112 : : // Initialize msg counters, callback, and background solution data
113 : 0 : m_numsent = 0;
114 : 0 : m_numreceived = 0;
115 : : //background();
116 : :
117 : : // Send vertex data to the collision detection library
118 : 0 : collideVertices();
119 : 0 : }
120 : :
121 : : void
122 : 0 : NodeSearch::background()
123 : : // *****************************************************************************
124 : : // Initialize dest mesh solution with background data
125 : : //! \details This is useful to see what points did not receive solution.
126 : : // *****************************************************************************
127 : : {
128 : 0 : tk::Fields& u = *m_u;
129 [ - - ]: 0 : for (std::size_t i = 0; i < u.nunk(); ++i) u(i,0) = -1.0;
130 : 0 : }
131 : :
132 : : void
133 : 0 : NodeSearch::collideVertices()
134 : : // *****************************************************************************
135 : : // Pass vertex information to the collision detection library
136 : : // *****************************************************************************
137 : : {
138 : 0 : const std::array< std::vector< double >, 3 >& coord = *m_coord;
139 : : auto nVertices = coord[0].size();
140 : : std::size_t nBoxes = 0;
141 : 0 : std::vector< bbox3d > boxes( nVertices );
142 [ - - ][ - - ]: 0 : std::vector< int > prio( nVertices );
143 : 0 : auto firstchunk = static_cast< int >( m_firstchunk );
144 [ - - ]: 0 : for (std::size_t i=0; i<nVertices; ++i) {
145 : : boxes[nBoxes].empty();
146 : 0 : boxes[nBoxes].add(CkVector3d(coord[0][i], coord[1][i], coord[2][i]));
147 : 0 : prio[nBoxes] = firstchunk;
148 : 0 : ++nBoxes;
149 : : }
150 : :
151 [ - - ]: 0 : CollideBoxesPrio( g_collideHandle, firstchunk + thisIndex,
152 : : static_cast<int>(nBoxes), boxes.data(), prio.data() );
153 : 0 : }
154 : :
155 : : void
156 : 0 : NodeSearch::collideTets() const
157 : : // *****************************************************************************
158 : : // Pass tet information to the collision detection library
159 : : // *****************************************************************************
160 : : {
161 : 0 : const std::vector< std::size_t >& inpoel = *m_inpoel;
162 : 0 : const std::array< std::vector< double >, 3 >& coord = *m_coord;
163 : 0 : auto nBoxes = inpoel.size() / 4;
164 : 0 : std::vector< bbox3d > boxes( nBoxes );
165 [ - - ][ - - ]: 0 : std::vector< int > prio( nBoxes );
166 : 0 : auto firstchunk = static_cast< int >( m_firstchunk );
167 [ - - ]: 0 : for (std::size_t i=0; i<nBoxes; ++i) {
168 : : boxes[i].empty();
169 : 0 : prio[i] = firstchunk;
170 [ - - ]: 0 : for (std::size_t j=0; j<4; ++j) {
171 : : // Get index of the jth point of the ith tet
172 : 0 : auto p = inpoel[i*4+j];
173 : : // Add that point to the tets bounding box
174 : 0 : boxes[i].add( CkVector3d( coord[0][p], coord[1][p], coord[2][p] ) );
175 : : }
176 : : }
177 : :
178 [ - - ]: 0 : CollideBoxesPrio( g_collideHandle, firstchunk + thisIndex,
179 : : static_cast<int>(nBoxes), boxes.data(), prio.data() );
180 : 0 : }
181 : :
182 : : void
183 : 0 : NodeSearch::processCollisions( const MeshData& src,
184 : : int nColl,
185 : : Collision* colls )
186 : : // *****************************************************************************
187 : : // Process potential collisions by sending my points to the source mesh chares
188 : : // that they potentially collide with.
189 : : //! \param[in] src Source mesh config data
190 : : //! \param[in] nColl Number of potential collisions to process
191 : : //! \param[in] colls List of potential collisions
192 : : // *****************************************************************************
193 : : {
194 : 0 : const std::array< std::vector< double >, 3 >& coord = *m_coord;
195 : 0 : int mychunk = m_firstchunk + thisIndex;
196 : :
197 : : std::vector< std::vector< PotentialCollision > >
198 : 0 : pColls( static_cast<std::size_t>(src.nchare) );
199 : :
200 : : // Separate potential collisions into lists based on the source mesh chare
201 : : // that is involved in the potential collision
202 [ - - ]: 0 : for (int i=0; i<nColl; ++i) {
203 : : int chareindex;
204 : : PotentialCollision pColl;
205 [ - - ]: 0 : if (colls[i].A.chunk == mychunk) {
206 : 0 : chareindex = colls[i].B.chunk - src.firstchunk;
207 : 0 : pColl.dst = static_cast<std::size_t>(colls[i].A.number);
208 : 0 : pColl.src = static_cast<std::size_t>(colls[i].B.number);
209 : : } else {
210 : 0 : chareindex = colls[i].A.chunk - src.firstchunk;
211 : 0 : pColl.dst = static_cast<std::size_t>(colls[i].B.number);
212 : 0 : pColl.src = static_cast<std::size_t>(colls[i].A.number);
213 : : }
214 : :
215 : : #if defined(STRICT_GNUC)
216 : : #pragma GCC diagnostic push
217 : : #pragma GCC diagnostic ignored "-Wdeprecated-copy"
218 : : #endif
219 [ - - ]: 0 : auto dst = pColl.dst;
220 : 0 : pColl.point = { coord[0][dst], coord[1][dst], coord[2][dst] };
221 : : #if defined(STRICT_GNUC)
222 : : #pragma GCC diagnostic pop
223 : : #endif
224 : :
225 [ - - ]: 0 : pColls[ static_cast<std::size_t>(chareindex) ].push_back( pColl );
226 : : }
227 : :
228 : : // Send out the lists of potential collisions to the source mesh chares
229 [ - - ]: 0 : for (int i=0; i<src.nchare; ++i) {
230 : 0 : auto I = static_cast< std::size_t >( i );
231 : 0 : m_numsent++;
232 [ - - ][ - - ]: 0 : src.proxy[i].determineActualCollisions( thisProxy, thisIndex,
233 : : static_cast<int>(pColls[I].size()), pColls[I].data() );
234 : : }
235 : 0 : }
236 : :
237 : : void
238 : 0 : NodeSearch::determineActualCollisions( CProxy_NodeSearch proxy,
239 : : int index,
240 : : int nColls,
241 : : PotentialCollision* colls )
242 : : // *****************************************************************************
243 : : // Identify actual collisions by calling intet on all possible collisions, and
244 : : // interpolate solution values to send back to the destination mesh.
245 : : //! \param[in] proxy The proxy of the destination mesh chare array
246 : : //! \param[in] index The index in proxy to return the solution data to
247 : : //! \param[in] nColls Number of collisions to be checked
248 : : //! \param[in] colls List of potential collisions
249 : : // *****************************************************************************
250 : : {
251 : 0 : const std::vector< std::size_t >& inpoel = *m_inpoel;
252 : 0 : const tk::Fields& u = *m_u;
253 : 0 : const std::vector< double >& f = *m_flag;
254 : : //CkPrintf( "Source chare %i received data for %i potential collisions\n",
255 : : // thisIndex, nColls);
256 : :
257 : : // Slightly shift dest mesh if symmetry is specified
258 : : auto eps = std::numeric_limits< tk::real >::epsilon();
259 : : const auto& sym = inciter::g_cfg.get< tag::overset, tag::sym_ >()[ 1 ];
260 [ - - ]: 0 : std::array< double, 3 > shift{ 0.0, 0.0, sym=="z" ? eps : 0.0 };
261 : :
262 : : // Iterate over my potential collisions and call intet() to determine
263 : : // if an actual collision occurred, and if so interpolate solution to dest
264 : : int numInTet = 0;
265 : : std::vector< SolutionData > exp;
266 [ - - ]: 0 : for (int i=0; i<nColls; ++i) {
267 : 0 : const auto& p = colls[i].point;
268 [ - - ]: 0 : std::vector< double > point{ p.x, p.y, p.z };
269 : : std::array< double, 4 > N;
270 [ - - ][ - - ]: 0 : if (tk::intet( *m_coord, *m_inpoel, point, colls[i].src, N, shift )) {
271 : : ++numInTet;
272 : : SolutionData data;
273 : 0 : data.dst = colls[i].dst;
274 : 0 : auto e = colls[i].src;
275 [ - - ]: 0 : const auto A = inpoel[e*4+0];
276 : 0 : const auto B = inpoel[e*4+1];
277 : 0 : const auto C = inpoel[e*4+2];
278 [ - - ]: 0 : const auto D = inpoel[e*4+3];
279 [ - - ]: 0 : data.sol.resize( u.nprop() + 1 );
280 [ - - ]: 0 : for (std::size_t c=0; c<u.nprop(); ++c) {
281 : 0 : data.sol[c] = N[0]*u(A,c) + N[1]*u(B,c) + N[2]*u(C,c) + N[3]*u(D,c);
282 : : }
283 [ - - ]: 0 : data.sol.back() = N[0]*f[A] + N[1]*f[B] + N[2]*f[C] + N[3]*f[D];
284 [ - - ]: 0 : exp.push_back( data );
285 : : }
286 : : }
287 : :
288 : : //if (numInTet) {
289 : : // CkPrintf( "Source chare %i found %i/%i actual collisions\n",
290 : : // thisIndex, numInTet, nColls );
291 : : //}
292 : :
293 : : // Send the solution data for the actual collisions back to the dest mesh
294 [ - - ][ - - ]: 0 : proxy[index].transferSolution( exp );
295 [ - - ]: 0 : if (not m_srcnotified) {
296 [ - - ]: 0 : m_done.send();
297 : 0 : m_srcnotified = 1;
298 : : }
299 : 0 : }
300 : :
301 : : void
302 : 0 : NodeSearch::transferSolution( const std::vector< SolutionData >& data )
303 : : // *****************************************************************************
304 : : // Transfer the interpolated solution data to destination mesh
305 : : //! \param[in] data List of solutions at nodes (multiple scalars transferred)
306 : : // *****************************************************************************
307 : : {
308 : 0 : tk::Fields& u = *m_u;
309 : 0 : std::vector< double >& f = *m_flag;
310 : : //auto eps = std::numeric_limits< tk::real >::epsilon();
311 : :
312 [ - - ]: 0 : for (std::size_t i=0; i<data.size(); ++i) {
313 : : const auto& d = data[i];
314 [ - - ]: 0 : if (m_trflag) { // transfer flags only
315 [ - - ]: 0 : if (m_dir) f[d.dst] = d.sol.back(); // overset to background only
316 : : }
317 : : else { // transfer solution only
318 [ - - ]: 0 : if (m_dir) { // overset to background
319 [ - - ]: 0 : if (f[d.dst] > 0.5) {
320 [ - - ]: 0 : for (std::size_t c=0; c<u.nprop(); ++c) u(d.dst,c) = d.sol[c];
321 : : }
322 : : }
323 : : else { // background to overset
324 [ - - ][ - - ]: 0 : if (f[d.dst] > -0.5 and f[d.dst] < 0.5) {
325 [ - - ]: 0 : for (std::size_t c=0; c<u.nprop(); ++c) u(d.dst,c) = d.sol[c];
326 : : }
327 : : }
328 : : }
329 : : }
330 : :
331 : : // Inform the caller if we've received all solution data
332 : 0 : m_numreceived++;
333 [ - - ]: 0 : if (m_numreceived == m_numsent) {
334 : 0 : m_done.send();
335 : : }
336 : 0 : }
337 : :
338 : : #include "NoWarning/nodesearch.def.h"
|
