Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/Refiner.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Mesh refiner for interfacing the mesh refinement library
9 : : \see Refiner.h for more info.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include <vector>
14 : : #include <algorithm>
15 : :
16 : : #include "Refiner.hpp"
17 : : #include "Reorder.hpp"
18 : : #include "AMR/mesh_adapter.hpp"
19 : : #include "AMR/Error.hpp"
20 : : #include "Inciter/InputDeck/InputDeck.hpp"
21 : : #include "CGPDE.hpp"
22 : : #include "DGPDE.hpp"
23 : : #include "DerivedData.hpp"
24 : : #include "UnsMesh.hpp"
25 : : #include "Centering.hpp"
26 : : #include "Around.hpp"
27 : : #include "Sorter.hpp"
28 : : #include "Discretization.hpp"
29 : :
30 : : namespace inciter {
31 : :
32 : : extern ctr::InputDeck g_inputdeck;
33 : : extern ctr::InputDeck g_inputdeck_defaults;
34 : : extern std::vector< CGPDE > g_cgpde;
35 : : extern std::vector< DGPDE > g_dgpde;
36 : :
37 : : } // inciter::
38 : :
39 : : using inciter::Refiner;
40 : :
41 : 2106 : Refiner::Refiner( std::size_t meshid,
42 : : const CProxy_Transporter& transporter,
43 : : const CProxy_Sorter& sorter,
44 : : const tk::CProxy_MeshWriter& meshwriter,
45 : : const std::vector< Scheme >& scheme,
46 : : const tk::RefinerCallback& cbr,
47 : : const tk::SorterCallback& cbs,
48 : : const std::vector< std::size_t >& ginpoel,
49 : : const tk::UnsMesh::CoordMap& coordmap,
50 : : const std::map< int, std::vector< std::size_t > >& bface,
51 : : const std::vector< std::size_t >& triinpoel,
52 : : const std::map< int, std::vector< std::size_t > >& bnode,
53 : 2106 : int nchare ) :
54 : : m_meshid( meshid ),
55 : : m_host( transporter ),
56 : : m_sorter( sorter ),
57 : : m_meshwriter( meshwriter ),
58 : : m_scheme( scheme ),
59 : : m_cbr( cbr ),
60 : : m_cbs( cbs ),
61 : : m_ginpoel( ginpoel ),
62 : : m_el( tk::global2local( ginpoel ) ), // fills m_inpoel, m_gid, m_lid
63 : : m_coordmap( coordmap ),
64 : : m_coord( flatcoord(coordmap) ),
65 : : m_bface( bface ),
66 : : m_bnode( bnode ),
67 : : m_triinpoel( triinpoel ),
68 : : m_nchare( nchare ),
69 : : m_mode( RefMode::T0REF ),
70 : : m_initref( g_inputdeck.get< tag::amr, tag::init >() ),
71 [ + - ]: 2106 : m_ninitref( g_inputdeck.get< tag::amr, tag::init >().size() ),
72 : 2106 : m_refiner( m_inpoel ),
73 : : m_nref( 0 ),
74 : : m_nbnd( 0 ),
75 : : m_extra( 0 ),
76 : : m_ch(),
77 : : m_edgech(),
78 : : m_chedge(),
79 : : m_localEdgeData(),
80 : : m_remoteEdgeData(),
81 : : m_nodeCommMap(),
82 : : m_oldTets(),
83 : : m_addedNodes(),
84 : : m_addedTets(),
85 : : m_removedNodes(),
86 : : m_oldntets( 0 ),
87 : : m_coarseBndFaces(),
88 : : m_coarseBndNodes(),
89 : : m_rid( m_coord[0].size() ),
90 : : m_oldrid(),
91 : : m_lref( m_rid.size() ),
92 : : m_parent(),
93 : : m_writeCallback(),
94 : : m_outref_ginpoel(),
95 : : m_outref_el(),
96 : : m_outref_coord(),
97 : : m_outref_addedNodes(),
98 : : m_outref_addedTets(),
99 : : m_outref_nodeCommMap(),
100 : : m_outref_bface(),
101 : : m_outref_bnode(),
102 [ + - ][ + - ]: 8424 : m_outref_triinpoel()
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - - ]
103 : : // *****************************************************************************
104 : : // Constructor
105 : : //! \param[in] meshid Mesh ID
106 : : //! \param[in] transporter Transporter (host) proxy
107 : : //! \param[in] sorter Mesh reordering (sorter) proxy
108 : : //! \param[in] meshwriter Mesh writer proxy
109 : : //! \param[in] scheme Discretization schemes (one per mesh)
110 : : //! \param[in] cbr Charm++ callbacks for Refiner
111 : : //! \param[in] cbs Charm++ callbacks for Sorter
112 : : //! \param[in] ginpoel Mesh connectivity (this chare) using global node IDs
113 : : //! \param[in] coordmap Mesh node coordinates (this chare) for global node IDs
114 : : //! \param[in] bface File-internal elem ids of side sets
115 : : //! \param[in] bnode Node lists of side sets
116 : : //! \param[in] triinpoel Triangle face connectivity with global IDs
117 : : //! \param[in] nchare Total number of refiner chares (chare array elements)
118 : : // *****************************************************************************
119 : : {
120 : : Assert( !m_ginpoel.empty(), "No elements assigned to refiner chare" );
121 : : Assert( tk::positiveJacobians( m_inpoel, m_coord ),
122 : : "Input mesh to Refiner Jacobian non-positive" );
123 : : Assert( !tk::leakyPartition(
124 : : tk::genEsuelTet( m_inpoel, tk::genEsup(m_inpoel,4) ),
125 : : m_inpoel, m_coord ),
126 : : "Input mesh to Refiner leaky" );
127 : :
128 : : #if not defined(__INTEL_COMPILER) || defined(NDEBUG)
129 : : // The above ifdef skips running the conformity test with the intel compiler
130 : : // in debug mode only. This is necessary because in tk::conforming(), filling
131 : : // up the map can fail with some meshes (only in parallel), e.g., tube.exo,
132 : : // used by some regression tests, due to the intel compiler generating some
133 : : // garbage incorrect code - only in debug, only in parallel, only with that
134 : : // mesh.
135 : : Assert( tk::conforming( m_inpoel, m_coord, true, m_rid ),
136 : : "Input mesh to Refiner not conforming" );
137 : : #endif
138 : :
139 : : // Generate local -> refiner lib node id map and its inverse
140 [ + - ]: 2106 : libmap();
141 : :
142 : : // Reverse initial mesh refinement type list (will pop from back)
143 : : std::reverse( begin(m_initref), end(m_initref) );
144 : :
145 : : // Generate boundary data structures for coarse mesh
146 [ + - ]: 2106 : coarseBnd();
147 : :
148 : : // If initial mesh refinement is configured, start initial mesh refinement.
149 : : // See also tk::grm::check_amr_errors in Control/Inciter/InputDeck/Ggrammar.h.
150 [ + + ][ + - ]: 2106 : if (g_inputdeck.get< tag::amr, tag::t0ref >() && m_ninitref > 0)
151 [ + - ]: 54 : t0ref();
152 : : else
153 [ + - ]: 2052 : endt0ref();
154 : 2106 : }
155 : :
156 : : void
157 : 2106 : Refiner::libmap()
158 : : // *****************************************************************************
159 : : // (Re-)generate local -> refiner lib node id map and its inverse
160 : : // *****************************************************************************
161 : : {
162 : : // Fill initial (matching) mapping between local and refiner node ids
163 : : std::iota( begin(m_rid), end(m_rid), 0 );
164 : :
165 : : // Fill in inverse, mapping refiner to local node ids
166 : : std::size_t i = 0;
167 [ + + ]: 647198 : for (auto r : m_rid) m_lref[r] = i++;
168 : 2106 : }
169 : :
170 : : void
171 : 2128 : Refiner::coarseBnd()
172 : : // *****************************************************************************
173 : : // (Re-)generate boundary data structures for coarse mesh
174 : : // *****************************************************************************
175 : : {
176 : : // Generate unique set of faces for each side set of the input (coarsest) mesh
177 : : m_coarseBndFaces.clear();
178 [ + + ]: 5018 : for (const auto& [ setid, faceids ] : m_bface) {
179 : : auto& faces = m_coarseBndFaces[ setid ];
180 [ + + ]: 315862 : for (auto f : faceids) {
181 : 312972 : faces.insert(
182 : 312972 : {{{ m_triinpoel[f*3+0], m_triinpoel[f*3+1], m_triinpoel[f*3+2] }}} );
183 : : }
184 : : }
185 : :
186 : : // Generate unique set of nodes for each side set of the input (coarsest) mesh
187 : : m_coarseBndNodes.clear();
188 [ + + ]: 3302 : for (const auto& [ setid, nodes ] : m_bnode) {
189 : : m_coarseBndNodes[ setid ].insert( begin(nodes), end(nodes) );
190 : : }
191 : 2128 : }
192 : :
193 : : void
194 : 2106 : Refiner::sendProxy()
195 : : // *****************************************************************************
196 : : // Send Refiner proxy to Discretization objects
197 : : //! \details This should be called when bound Discretization chare array
198 : : //! elements have already been created.
199 : : // *****************************************************************************
200 : : {
201 : : // Make sure (bound) Discretization chare is already created and accessible
202 : : Assert( m_scheme[m_meshid].disc()[thisIndex].ckLocal() != nullptr,
203 : : "About to dereference nullptr" );
204 : :
205 : : // Pass Refiner Charm++ chare proxy to fellow (bound) Discretization object
206 [ + - ]: 4212 : m_scheme[m_meshid].disc()[thisIndex].ckLocal()->setRefiner( thisProxy );
207 : 2106 : }
208 : :
209 : : void
210 : 22 : Refiner::reorder()
211 : : // *****************************************************************************
212 : : // Query Sorter and update local mesh with the reordered one
213 : : // *****************************************************************************
214 : : {
215 : 22 : m_sorter[thisIndex].ckLocal()->
216 [ + - ]: 44 : mesh( m_ginpoel, m_coordmap, m_triinpoel, m_bnode );
217 : :
218 : : // Update local mesh data based on data just received from Sorter
219 : 22 : m_el = tk::global2local( m_ginpoel ); // fills m_inpoel, m_gid, m_lid
220 : 22 : m_coord = flatcoord( m_coordmap );
221 : :
222 : : // Re-generate boundary data structures for coarse mesh
223 : 22 : coarseBnd();
224 : :
225 : : // WARNING: This re-creates the AMR lib which is probably not what we
226 : : // ultimately want, beacuse this deletes its history recorded during initial
227 : : // (t<0) refinement. However, this appears to correctly update the local mesh
228 : : // based on the reordered one (from Sorter) at least when t0ref is off.
229 : 22 : m_refiner = AMR::mesh_adapter_t( m_inpoel );
230 : 22 : }
231 : :
232 : : tk::UnsMesh::Coords
233 : 2128 : Refiner::flatcoord( const tk::UnsMesh::CoordMap& coordmap )
234 : : // *****************************************************************************
235 : : // Generate flat coordinate data from coordinate map
236 : : //! \param[in] coordmap Coordinates associated to global node IDs of mesh chunk
237 : : //! \return Flat coordinate data
238 : : // *****************************************************************************
239 : : {
240 : : tk::UnsMesh::Coords coord;
241 : :
242 : : // Convert node coordinates associated to global node IDs to a flat vector
243 : : auto npoin = coordmap.size();
244 : : Assert( m_gid.size() == npoin, "Size mismatch" );
245 [ + - ]: 2128 : coord[0].resize( npoin );
246 [ + - ]: 2128 : coord[1].resize( npoin );
247 [ + - ]: 2128 : coord[2].resize( npoin );
248 [ + + ]: 653930 : for (const auto& [ gid, coords ] : coordmap) {
249 : 651802 : auto i = tk::cref_find( m_lid, gid );
250 : : Assert( i < npoin, "Indexing out of coordinate map" );
251 : 651802 : coord[0][i] = coords[0];
252 : 651802 : coord[1][i] = coords[1];
253 : 651802 : coord[2][i] = coords[2];
254 : : }
255 : :
256 : 2128 : return coord;
257 : : }
258 : :
259 : : void
260 : 184 : Refiner::dtref( const std::map< int, std::vector< std::size_t > >& bface,
261 : : const std::map< int, std::vector< std::size_t > >& bnode,
262 : : const std::vector< std::size_t >& triinpoel )
263 : : // *****************************************************************************
264 : : // Start mesh refinement (during time stepping, t>0)
265 : : //! \param[in] bface Boundary-faces mapped to side set ids
266 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
267 : : //! \param[in] triinpoel Boundary-face connectivity
268 : : // *****************************************************************************
269 : : {
270 : 184 : m_mode = RefMode::DTREF;
271 : :
272 : : // Update boundary node lists
273 : : m_bface = bface;
274 : : m_bnode = bnode;
275 : 184 : m_triinpoel = triinpoel;
276 : :
277 : 184 : start();
278 : 184 : }
279 : :
280 : : void
281 : 66 : Refiner::outref( const std::map< int, std::vector< std::size_t > >& bface,
282 : : const std::map< int, std::vector< std::size_t > >& bnode,
283 : : const std::vector< std::size_t >& triinpoel,
284 : : CkCallback c,
285 : : RefMode mode )
286 : : // *****************************************************************************
287 : : // Start mesh refinement (for field output)
288 : : //! \param[in] bface Boundary-faces mapped to side set ids
289 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
290 : : //! \param[in] triinpoel Boundary-face connectivity
291 : : //! \param[in] c Function to continue with after the writing field output
292 : : //! \param[in] mode Refinement mode
293 : : // *****************************************************************************
294 : : {
295 : 66 : m_mode = mode;
296 : :
297 : 66 : m_writeCallback = c;
298 : :
299 : : // Update boundary node lists
300 : : m_bface = bface;
301 : : m_bnode = bnode;
302 : 66 : m_triinpoel = triinpoel;
303 : :
304 : 66 : start();
305 : 66 : }
306 : :
307 : : void
308 : 133 : Refiner::t0ref()
309 : : // *****************************************************************************
310 : : // Output mesh to file before a new step mesh refinement
311 : : // *****************************************************************************
312 : : {
313 : : Assert( m_ninitref > 0, "No initial mesh refinement steps configured" );
314 : : // Output initial mesh to file
315 [ + + ]: 133 : auto l = m_ninitref - m_initref.size(); // num initref steps completed
316 : 133 : auto t0 = g_inputdeck.get< tag::discr, tag::t0 >();
317 [ + + ]: 133 : if (l == 0) {
318 [ + - ][ + - ]: 108 : writeMesh( "t0ref", l, t0-1.0,
319 [ + - ][ - + ]: 162 : CkCallback( CkIndex_Refiner::start(), thisProxy[thisIndex] ) );
[ - - ]
320 : : } else {
321 : 79 : start();
322 : : }
323 : 133 : }
324 : :
325 : : void
326 : 383 : Refiner::start()
327 : : // *****************************************************************************
328 : : // Start new step of mesh refinement
329 : : // *****************************************************************************
330 : : {
331 : 383 : m_extra = 0;
332 : : m_ch.clear();
333 : : m_remoteEdgeData.clear();
334 : : m_remoteEdges.clear();
335 : :
336 : 383 : updateEdgeData();
337 : :
338 : : // Generate and communicate boundary edges
339 : 383 : bndEdges();
340 : 383 : }
341 : :
342 : : void
343 : 383 : Refiner::bndEdges()
344 : : // *****************************************************************************
345 : : // Generate boundary edges and send them to all chares
346 : : //! \details Extract edges on the boundary only. The boundary edges (shared by
347 : : //! multiple chares) will be agreed on a refinement that yields a conforming
348 : : //! mesh across chares boundaries.
349 : : // *****************************************************************************
350 : : {
351 : : // Compute the number of edges (chunksize) a chare will respond to when
352 : : // computing shared edges
353 : 383 : auto N = static_cast< std::size_t >( m_nchare );
354 [ + - ]: 383 : std::size_t chunksize = std::numeric_limits< std::size_t >::max() / N;
355 : :
356 : : // Generate boundary edges of our mesh chunk
357 : : std::unordered_map< int, EdgeSet > chbedges;
358 [ + - ]: 766 : auto esup = tk::genEsup( m_inpoel, 4 ); // elements surrounding points
359 [ + - ]: 383 : auto esuel = tk::genEsuelTet( m_inpoel, esup ); // elems surrounding elements
360 [ + + ]: 426073 : for (std::size_t e=0; e<esuel.size()/4; ++e) {
361 : 425690 : auto mark = e*4;
362 [ + + ]: 2128450 : for (std::size_t f=0; f<4; ++f) {
363 [ + + ]: 1702760 : if (esuel[mark+f] == -1) {
364 [ + - ]: 151770 : auto A = m_ginpoel[ mark+tk::lpofa[f][0] ];
365 : 151770 : auto B = m_ginpoel[ mark+tk::lpofa[f][1] ];
366 : 151770 : auto C = m_ginpoel[ mark+tk::lpofa[f][2] ];
367 : : Assert( m_lid.find( A ) != end(m_lid), "Local node ID not found" );
368 : : Assert( m_lid.find( B ) != end(m_lid), "Local node ID not found" );
369 : : Assert( m_lid.find( C ) != end(m_lid), "Local node ID not found" );
370 : : // assign edges to bins a single chare will respond to when computing
371 : : // shared edges
372 : 151770 : auto bin = A / chunksize;
373 : : Assert( bin < N, "Will index out of number of chares" );
374 [ + - ][ + - ]: 151770 : chbedges[ static_cast<int>(bin) ].insert( {A,B} );
375 : 151770 : bin = B / chunksize;
376 : : Assert( bin < N, "Will index out of number of chares" );
377 [ + - ][ + - ]: 151770 : chbedges[ static_cast<int>(bin) ].insert( {B,C} );
378 : 151770 : bin = C / chunksize;
379 : : Assert( bin < N, "Will index out of number of chares" );
380 [ + - ][ + - ]: 151770 : chbedges[ static_cast<int>(bin) ].insert( {C,A} );
381 : : }
382 : : }
383 : : }
384 : :
385 : : // Send edges in bins to chares that will compute shared edges
386 : 383 : m_nbnd = chbedges.size();
387 [ - + ]: 383 : if (m_nbnd == 0)
388 [ - - ]: 0 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
389 : : m_cbr.get< tag::queried >() );
390 : : else
391 [ + + ]: 1440 : for (const auto& [ targetchare, bndedges ] : chbedges)
392 [ + - ][ + - ]: 2114 : thisProxy[ targetchare ].query( thisIndex, bndedges );
[ - - ]
393 : 383 : }
394 : :
395 : : void
396 : 1057 : Refiner::query( int fromch, const EdgeSet& edges )
397 : : // *****************************************************************************
398 : : // Incoming query for a list boundary edges for which this chare compiles
399 : : // shared edges
400 : : //! \param[in] fromch Sender chare ID
401 : : //! \param[in] edges Chare-boundary edge list from another chare
402 : : // *****************************************************************************
403 : : {
404 : : // Store incoming edges in edge->chare and its inverse, chare->edge, maps
405 [ + + ]: 571221 : for (const auto& e : edges) m_edgech[ e ].push_back( fromch );
406 : : m_chedge[ fromch ].insert( begin(edges), end(edges) );
407 : : // Report back to chare message received from
408 [ + - ]: 1057 : thisProxy[ fromch ].recvquery();
409 : 1057 : }
410 : :
411 : : void
412 : 1057 : Refiner::recvquery()
413 : : // *****************************************************************************
414 : : // Receive receipt of boundary edge lists to query
415 : : // *****************************************************************************
416 : : {
417 [ + + ]: 1057 : if (--m_nbnd == 0)
418 : 383 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
419 : : m_cbr.get< tag::queried >() );
420 : 1057 : }
421 : :
422 : : void
423 : 383 : Refiner::response()
424 : : // *****************************************************************************
425 : : // Respond to boundary edge list queries
426 : : // *****************************************************************************
427 : : {
428 : : std::unordered_map< int, std::vector< int > > exp;
429 : :
430 : : // Compute shared edges whose chare ids will be sent back to querying chares
431 [ + + ]: 1440 : for (const auto& [ neighborchare, bndedges ] : m_chedge) {
432 : : auto& e = exp[ neighborchare ];
433 [ + + ]: 286139 : for (const auto& ed : bndedges)
434 [ + + ]: 630722 : for (auto d : tk::cref_find(m_edgech,ed))
435 [ + + ]: 345640 : if (d != neighborchare)
436 [ + - ]: 60558 : e.push_back( d );
437 : : }
438 : :
439 : : // Send chare ids of shared edges to chares that issued a query to us. Shared
440 : : // boundary edges assigned to chare ids sharing the boundary edge were
441 : : // computed above for those chares that queried this map from us. These
442 : : // boundary edges form a distributed table and we only work on a chunk of it.
443 : : // Note that we only send data back to those chares that have queried us. The
444 : : // receiving sides do not know in advance if they receive messages or not.
445 : : // Completion is detected by having the receiver respond back and counting
446 : : // the responses on the sender side, i.e., this chare.
447 : 383 : m_nbnd = exp.size();
448 [ + + ]: 383 : if (m_nbnd == 0)
449 [ + - ]: 157 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
450 : : m_cbr.get< tag::responded >() );
451 : : else
452 [ + + ]: 1283 : for (const auto& [ targetchare, bndedges ] : exp)
453 [ + - ][ + - ]: 2114 : thisProxy[ targetchare ].bnd( thisIndex, bndedges );
454 : 383 : }
455 : :
456 : : void
457 : 1057 : Refiner::bnd( int fromch, const std::vector< int >& chares )
458 : : // *****************************************************************************
459 : : // Receive shared boundary edges for our mesh chunk
460 : : //! \param[in] fromch Sender chare ID
461 : : //! \param[in] chares Chare ids we share edges with
462 : : // *****************************************************************************
463 : : {
464 : : // Store chare ids we share edges with
465 : : m_ch.insert( begin(chares), end(chares) );
466 : :
467 : : // Report back to chare message received from
468 [ + - ]: 1057 : thisProxy[ fromch ].recvbnd();
469 : 1057 : }
470 : :
471 : : void
472 : 1057 : Refiner::recvbnd()
473 : : // *****************************************************************************
474 : : // Receive receipt of shared boundary edges
475 : : // *****************************************************************************
476 : : {
477 [ + + ]: 1057 : if (--m_nbnd == 0)
478 : 226 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
479 : : m_cbr.get< tag::responded >() );
480 : 1057 : }
481 : :
482 : : void
483 : 383 : Refiner::refine()
484 : : // *****************************************************************************
485 : : // Do a single step of mesh refinement (really, only tag edges)
486 : : //! \details During initial (t<0, t0ref) mesh refinement, this is a single step
487 : : //! in a potentially multiple-entry list of initial adaptive mesh refinement
488 : : //! steps. Distribution of the chare-boundary edges must have preceded this
489 : : //! step, so that boundary edges (shared by multiple chares) can agree on a
490 : : //! refinement that yields a conforming mesh across chare boundaries.
491 : : //!
492 : : //! During-timestepping (t>0, dtref) mesh refinement this is called once, as
493 : : //! we only do a single step during time stepping.
494 : : //!
495 : : //! During field-output (outref) mesh refinement, this may be called multiple
496 : : //! times, depending the number of refinement levels needed for the field
497 : : //! output.
498 : : // *****************************************************************************
499 : : {
500 : : // Free memory used for computing shared boundary edges
501 : 383 : tk::destroy( m_edgech );
502 : 383 : tk::destroy( m_chedge );
503 : :
504 : : // Perform leak test on old mesh
505 : : Assert( !tk::leakyPartition(
506 : : tk::genEsuelTet( m_inpoel, tk::genEsup(m_inpoel,4) ),
507 : : m_inpoel, m_coord ),
508 : : "Mesh partition before refinement leaky" );
509 : :
510 [ + + ]: 383 : if (m_mode == RefMode::T0REF) {
511 : :
512 : : // Refine mesh based on next initial refinement type
513 [ + - ]: 133 : if (!m_initref.empty()) {
514 : 133 : auto r = m_initref.back(); // consume (reversed) list from its back
515 [ + + ]: 133 : if (r == ctr::AMRInitialType::UNIFORM)
516 : 70 : uniformRefine();
517 [ + + ]: 63 : else if (r == ctr::AMRInitialType::UNIFORM_DEREFINE)
518 : 39 : uniformDeRefine();
519 [ + + ]: 24 : else if (r == ctr::AMRInitialType::INITIAL_CONDITIONS)
520 : 22 : errorRefine();
521 [ + - ]: 2 : else if (r == ctr::AMRInitialType::COORDINATES)
522 : 2 : coordRefine();
523 [ - - ]: 0 : else if (r == ctr::AMRInitialType::EDGELIST)
524 : 0 : edgelistRefine();
525 [ - - ][ - - ]: 0 : else Throw( "Initial AMR type not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
526 : : }
527 : :
528 [ + + ]: 250 : } else if (m_mode == RefMode::DTREF) {
529 : :
530 [ + - ]: 184 : if (g_inputdeck.get< tag::amr, tag::dtref_uniform >())
531 : 184 : uniformRefine();
532 : : else
533 : 0 : errorRefine();
534 : :
535 [ + + ]: 66 : } else if (m_mode == RefMode::OUTREF) {
536 : :
537 : 33 : uniformRefine();
538 : :
539 [ + - ]: 33 : } else if (m_mode == RefMode::OUTDEREF) {
540 : :
541 : 33 : uniformDeRefine();
542 : :
543 [ - - ][ - - ]: 0 : } else Throw( "RefMode not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
544 : :
545 : : // Communicate extra edges
546 : 383 : comExtra();
547 : 383 : }
548 : :
549 : : void
550 [ + + ]: 391 : Refiner::comExtra()
551 : : // *****************************************************************************
552 : : // Communicate extra edges along chare boundaries
553 : : // *****************************************************************************
554 : : {
555 : : // Export extra added nodes on our mesh chunk boundary to other chares
556 [ + + ]: 391 : if (m_ch.empty()) {
557 : 51 : correctref();
558 : : } else {
559 [ + + ]: 1876 : for (auto c : m_ch) { // for all chares we share at least an edge with
560 [ + - ][ - + ]: 4608 : thisProxy[c].addRefBndEdges(thisIndex, m_localEdgeData, m_intermediates);
[ - - ]
561 : : }
562 : : }
563 : 391 : }
564 : :
565 : : void
566 : 1536 : Refiner::addRefBndEdges(
567 : : int fromch,
568 : : const AMR::EdgeData& ed,
569 : : const std::unordered_set< std::size_t >& intermediates )
570 : : // *****************************************************************************
571 : : //! Receive edges on our chare boundary from other chares
572 : : //! \param[in] fromch Chare call coming from
573 : : //! \param[in] ed Edges on chare boundary
574 : : //! \param[in] intermediates Intermediate nodes
575 : : // *****************************************************************************
576 : : {
577 : : // Save/augment buffers of edge data for each sender chare
578 : : auto& red = m_remoteEdgeData[ fromch ];
579 : : auto& re = m_remoteEdges[ fromch ];
580 : : using edge_data_t = std::tuple< Edge, int, AMR::Edge_Lock_Case >;
581 [ + + ]: 737211 : for (const auto& [ edge, data ] : ed) {
582 : 735675 : red.push_back( edge_data_t{ edge, data.first, data.second } );
583 : 735675 : re.push_back( edge );
584 : : }
585 : :
586 : : // Add intermediates to mesh refiner lib
587 [ + + ]: 1973 : for (const auto g : intermediates) {
588 : 437 : auto l = m_lid.find( g ); // convert to local node ids
589 [ + + ]: 437 : if (l != end(m_lid)) {
590 : 10 : m_refiner.tet_store.intermediate_list.insert( m_rid[l->second] );
591 : : }
592 : : }
593 : :
594 : : // Heard from every worker we share at least a single edge with
595 [ + + ]: 1536 : if (++m_nref == m_ch.size()) {
596 : 340 : m_nref = 0;
597 : : // Add intermediates to refiner lib
598 : : auto localedges_orig = m_localEdgeData;
599 : : //auto intermediates_orig = m_intermediates;
600 [ + - ]: 340 : m_refiner.lock_intermediates();
601 : : // Run compatibility algorithm
602 [ + - ]: 340 : m_refiner.mark_refinement();
603 : : // Update edge data from mesh refiner
604 [ + - ]: 340 : updateEdgeData();
605 : : // If refiner lib modified our edges, need to recommunicate
606 [ + + ]: 340 : auto modified = static_cast< std::size_t >(
607 : : localedges_orig != m_localEdgeData ? 1 : 0 );
608 : : //int modified = ( (localedges_orig != m_localEdgeData ||
609 : : // intermediates_orig != m_intermediates) ? 1 : 0 );
610 [ + - ]: 340 : std::vector< std::size_t > meshdata{ m_meshid, modified };
611 [ + - ]: 340 : contribute( meshdata, CkReduction::max_ulong,
612 : : m_cbr.get< tag::compatibility >() );
613 : : }
614 : 1536 : }
615 : :
616 : : void
617 : 387 : Refiner::correctref()
618 : : // *****************************************************************************
619 : : // Correct extra edges to arrive at conforming mesh across chare boundaries
620 : : //! \details This function is called repeatedly until there is not a a single
621 : : //! edge that needs correction for the whole distributed problem to arrive at
622 : : //! a conforming mesh across chare boundaries during a mesh refinement step.
623 : : // *****************************************************************************
624 : : {
625 : : auto unlocked = AMR::Edge_Lock_Case::unlocked;
626 : :
627 : : // Storage for edge data that need correction to yield a conforming mesh
628 : : AMR::EdgeData extra;
629 : :
630 : : // loop through all edges shared with other chares
631 [ + + ]: 1919 : for (const auto& [ neighborchare, edgedata ] : m_remoteEdgeData) {
632 [ + + ]: 737207 : for (const auto& [edge,remote_needs_refining,remote_lock_case] : edgedata) {
633 : : // find local data of remote edge
634 : : auto it = m_localEdgeData.find( edge );
635 [ + + ]: 735675 : if (it != end(m_localEdgeData)) {
636 : : auto& local = it->second;
637 : : auto& local_needs_refining = local.first;
638 : 42614 : auto& local_lock_case = local.second;
639 : :
640 : 42614 : auto local_needs_refining_orig = local_needs_refining;
641 [ - + ]: 42614 : auto local_lock_case_orig = local_lock_case;
642 : :
643 : : Assert( !(local_lock_case > unlocked && local_needs_refining),
644 : : "Invalid local edge: locked & needs refining" );
645 : : Assert( !(remote_lock_case > unlocked && remote_needs_refining),
646 : : "Invalid remote edge: locked & needs refining" );
647 : :
648 : : // compute lock from local and remote locks as most restrictive
649 : 42614 : local_lock_case = std::max( local_lock_case, remote_lock_case );
650 : :
651 [ + + ]: 42614 : if (local_lock_case > unlocked)
652 : 72 : local_needs_refining = 0;
653 : :
654 [ + + ][ + + ]: 42614 : if (local_lock_case == unlocked && remote_needs_refining)
655 : 28494 : local_needs_refining = 1;
656 : :
657 : : // if the remote sent us data that makes us change our local state,
658 : : // e.g., local{0,0} + remote(1,0} -> local{1,0}, i.e., I changed my
659 : : // state I need to tell the world about it
660 [ + - ]: 42614 : if ( (local_lock_case != local_lock_case_orig ||
661 [ + + ]: 42614 : local_needs_refining != local_needs_refining_orig) ||
662 : : // or if the remote data is inconsistent with what I think, e.g.,
663 : : // local{1,0} + remote(0,0} -> local{1,0}, i.e., the remote does not
664 : : // yet agree, I need to tell the world about it
665 [ + + ]: 42607 : (local_lock_case != remote_lock_case ||
666 [ + + ]: 42595 : local_needs_refining != remote_needs_refining) )
667 : : {
668 : 66 : auto l1 = tk::cref_find( m_lid, edge[0] );
669 : 33 : auto l2 = tk::cref_find( m_lid, edge[1] );
670 : : Assert( l1 != l2, "Edge end-points local ids are the same" );
671 [ + + ]: 33 : auto r1 = m_rid[ l1 ];
672 [ + + ]: 33 : auto r2 = m_rid[ l2 ];
673 : : Assert( r1 != r2, "Edge end-points refiner ids are the same" );
674 [ + + ][ + - ]: 59 : extra[ {{ std::min(r1,r2), std::max(r1,r2) }} ] =
675 : : { local_needs_refining, local_lock_case };
676 : : }
677 : : }
678 : : }
679 : : }
680 : :
681 : : m_remoteEdgeData.clear();
682 : 387 : m_extra = extra.size();
683 : :
684 [ + + ]: 387 : if (!extra.empty()) {
685 : : // Do refinement including edges that need to be corrected
686 [ + - ]: 3 : m_refiner.mark_error_refinement_corr( extra );
687 : : // Update our extra-edge store based on refiner
688 [ + - ]: 3 : updateEdgeData();
689 : : }
690 : :
691 : : // Aggregate number of extra edges that still need correction and some
692 : : // refinement/derefinement statistics
693 : : const auto& tet_store = m_refiner.tet_store;
694 : : std::vector< std::size_t >
695 : 387 : m{ m_meshid,
696 : 387 : m_extra,
697 : : tet_store.marked_refinements.size(),
698 : : tet_store.marked_derefinements.size(),
699 [ + - ]: 387 : static_cast< std::underlying_type_t< RefMode > >( m_mode ) };
700 [ + - ]: 387 : contribute( m, CkReduction::sum_ulong, m_cbr.get< tag::matched >() );
701 : 387 : }
702 : :
703 : : void
704 : 1109 : Refiner::updateEdgeData()
705 : : // *****************************************************************************
706 : : // Query AMR lib and update our local store of edge data
707 : : // *****************************************************************************
708 : : {
709 : : m_localEdgeData.clear();
710 : : m_intermediates.clear();
711 : :
712 : : // This currently takes ALL edges from the AMR lib, i.e., on the whole
713 : : // domain. We should eventually only collect edges here that are shared with
714 : : // other chares.
715 : : const auto& ref_edges = m_refiner.tet_store.edge_store.edges;
716 : 1109 : const auto& refinpoel = m_refiner.tet_store.get_active_inpoel();
717 : :
718 [ + + ]: 1082686 : for (std::size_t e=0; e<refinpoel.size()/4; ++e) {
719 : 1081577 : auto A = refinpoel[e*4+0];
720 : 1081577 : auto B = refinpoel[e*4+1];
721 : 1081577 : auto C = refinpoel[e*4+2];
722 : 1081577 : auto D = refinpoel[e*4+3];
723 : : std::array<Edge,6> edges{{ {{A,B}}, {{B,C}}, {{A,C}},
724 : 1081577 : {{A,D}}, {{B,D}}, {{C,D}} }};
725 [ + + ]: 7571039 : for (const auto& ed : edges) {
726 [ + + ]: 10579240 : auto ae = AMR::edge_t{{{ std::min(ed[0],ed[1]), std::max(ed[0],ed[1]) }}};
727 : 6489462 : auto r = tk::cref_find( ref_edges, ae );
728 : 12978924 : const auto ged = Edge{{ m_gid[ tk::cref_find( m_lref, ed[0] ) ],
729 : 12978924 : m_gid[ tk::cref_find( m_lref, ed[1] ) ] }};
730 : : m_localEdgeData[ ged ] = { r.needs_refining, r.lock_case };
731 : : }
732 : : }
733 : :
734 : : // Build intermediates to send. This currently takes ALL intermediates from
735 : : // the AMR lib, i.e., on the whole domain. We should eventually only collect
736 : : // edges here that are shared with other chares.
737 [ + + ]: 2497 : for (const auto& r : m_refiner.tet_store.intermediate_list) {
738 : 2776 : m_intermediates.insert( m_gid[ tk::cref_find( m_lref, r ) ] );
739 : : }
740 : 1109 : }
741 : :
742 : : std::tuple< std::vector< std::string >,
743 : : std::vector< std::vector< tk::real > >,
744 : : std::vector< std::string >,
745 : : std::vector< std::vector< tk::real > > >
746 : 2420 : Refiner::refinementFields() const
747 : : // *****************************************************************************
748 : : // Collect mesh output fields from refiner lib
749 : : //! \return Names and fields of mesh refinement data in mesh cells and nodes
750 : : // *****************************************************************************
751 : : {
752 : : // Find number of nodes in current mesh
753 : 2420 : auto npoin = tk::npoin_in_graph( m_inpoel );
754 : : // Generate edges surrounding points in current mesh
755 : 4840 : auto esup = tk::genEsup( m_inpoel, 4 );
756 : :
757 : : // Update solution on current mesh
758 [ + - ]: 2420 : auto u = solution( npoin, esup );
759 : : Assert( u.nunk() == npoin, "Solution uninitialized or wrong size" );
760 : :
761 : : // Compute error in edges on current mesh
762 [ + - ]: 2420 : auto edgeError = errorsInEdges( npoin, esup, u );
763 : :
764 : : // Transfer error from edges to cells for field output
765 [ + - ]: 2420 : std::vector< tk::real > error( m_inpoel.size()/4, 0.0 );
766 [ + + ]: 2363884 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
767 : 2361464 : auto A = m_inpoel[e*4+0];
768 : 2361464 : auto B = m_inpoel[e*4+1];
769 : 2361464 : auto C = m_inpoel[e*4+2];
770 : 2361464 : auto D = m_inpoel[e*4+3];
771 : : std::array<Edge,6> edges{{ {{A,B}}, {{B,C}}, {{A,C}},
772 : 2361464 : {{A,D}}, {{B,D}}, {{C,D}} }};
773 : : // sum error from edges to elements
774 [ + + ]: 16530248 : for (const auto& ed : edges) error[e] += tk::cref_find( edgeError, ed );
775 : 2361464 : error[e] /= 6.0; // assign edge-average error to element
776 : : }
777 : :
778 : : // Prepare element fields with mesh refinement data
779 : : std::vector< std::string >
780 [ + - ][ + - ]: 4840 : elemfieldnames{ "refinement level", "cell type", "error" };
[ + - ][ - + ]
[ - - ]
781 : 2420 : auto& tet_store = m_refiner.tet_store;
782 : : std::vector< std::vector< tk::real > > elemfields{
783 : : tet_store.get_refinement_level_list(),
784 : : tet_store.get_cell_type_list(),
785 [ + - ][ + - ]: 4840 : error };
[ + - ][ - + ]
[ + - ][ - - ]
786 : :
787 : : using tuple_t = std::tuple< std::vector< std::string >,
788 : : std::vector< std::vector< tk::real > >,
789 : : std::vector< std::string >,
790 : : std::vector< std::vector< tk::real > > >;
791 : 4840 : return tuple_t{ elemfieldnames, elemfields, {}, {} };
792 : : }
793 : :
794 : : void
795 : 187 : Refiner::writeMesh( const std::string& basefilename,
796 : : uint64_t itr,
797 : : tk::real t,
798 : : CkCallback c ) const
799 : : // *****************************************************************************
800 : : // Output mesh to file(s)
801 : : //! \param[in] basefilename File name to append to
802 : : //! \param[in] itr Iteration count since a new mesh
803 : : //! \param[in] t "Physical time" to write to file. "Time" here is used to
804 : : //! designate a new time step at which the mesh is saved.
805 : : //! \param[in] c Function to continue with after the write
806 : : // *****************************************************************************
807 : : {
808 : 187 : auto r = refinementFields();
809 : : auto& elemfieldnames = std::get< 0 >( r );
810 : : auto& elemfields = std::get< 1 >( r );
811 : : auto& nodefieldnames = std::get< 2 >( r );
812 : : auto& nodefields = std::get< 3 >( r );
813 : :
814 : : // Prepare solution field names: depvar + component id for all eqs
815 : 187 : auto nprop = g_inputdeck.get< tag::component >().nprop();
816 [ + - ]: 374 : auto solfieldnames = g_inputdeck.get<tag::component>().depvar( g_inputdeck );
817 : : Assert( solfieldnames.size() == nprop, "Size mismatch" );
818 : :
819 : 187 : const auto scheme = g_inputdeck.get< tag::discr, tag::scheme >();
820 [ + - ][ + - ]: 187 : const auto centering = ctr::Scheme().centering( scheme );
821 [ + + ]: 187 : auto t0 = g_inputdeck.get< tag::discr, tag::t0 >();
822 : :
823 : : // list of nodes/elements at which box ICs are defined
824 : 187 : std::vector< std::unordered_set< std::size_t > > inbox;
825 : : tk::real V = 1.0;
826 : :
827 : : // Prepare node or element fields for output to file
828 [ + + ]: 187 : if (centering == tk::Centering::NODE) {
829 : :
830 : : // Augment element field names with solution variable names + field ids
831 : : nodefieldnames.insert( end(nodefieldnames),
832 [ + - ]: 121 : begin(solfieldnames), end(solfieldnames) );
833 : :
834 : : // Evaluate initial conditions on current mesh at t0
835 [ + - ]: 121 : tk::Fields u( m_coord[0].size(), nprop );
836 [ + + ][ + - ]: 242 : for (auto& eq : g_cgpde) eq.initialize( m_coord, u, t0, V, inbox );
837 : :
838 : : // Extract all scalar components from solution for output to file
839 [ + + ]: 274 : for (std::size_t i=0; i<nprop; ++i)
840 [ + - ][ - - ]: 306 : nodefields.push_back( u.extract( i, 0 ) );
841 : :
842 [ + - ]: 66 : } else if (centering == tk::Centering::ELEM) {
843 : :
844 : : // Augment element field names with solution variable names + field ids
845 : : elemfieldnames.insert( end(elemfieldnames),
846 [ + - ]: 66 : begin(solfieldnames), end(solfieldnames) );
847 : :
848 : 66 : auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
849 [ + - ]: 66 : tk::Fields lhs( m_inpoel.size()/4, ndof*nprop );
850 : :
851 : : // Generate left hand side for DG initialize
852 [ + - ]: 66 : auto geoElem = tk::genGeoElemTet( m_inpoel, m_coord );
853 [ + + ]: 132 : for (const auto& eq : g_dgpde) eq.lhs( geoElem, lhs );
854 : :
855 : : // Evaluate initial conditions on current mesh at t0
856 : : auto u = lhs;
857 [ + + ]: 132 : for (const auto& eq : g_dgpde)
858 [ + - ]: 66 : eq.initialize( lhs, m_inpoel, m_coord, inbox, u, t0, m_inpoel.size()/4 );
859 : :
860 : : // Extract all scalar components from solution for output to file
861 [ + + ]: 148 : for (std::size_t i=0; i<nprop; ++i)
862 [ + - ][ - - ]: 164 : elemfields.push_back( u.extract( i, 0 ) );
863 : : }
864 : :
865 : : // Output mesh
866 : : m_meshwriter[ CkNodeFirst( CkMyNode() ) ].
867 [ + - ][ - + ]: 561 : write( m_meshid, /*meshoutput = */ true, /*fieldoutput = */ true, itr, 1, t,
[ - - ]
868 [ + - ]: 187 : thisIndex, basefilename, m_inpoel, m_coord, m_bface,
869 [ + - ][ + - ]: 448 : tk::remap(m_bnode,m_lid), tk::remap(m_triinpoel,m_lid),
[ + + ][ - - ]
870 : : elemfieldnames, nodefieldnames, {}, elemfields, nodefields, {},
871 : : {}, c );
872 : 187 : }
873 : :
874 : : void
875 : 383 : Refiner::perform()
876 : : // *****************************************************************************
877 : : // Perform mesh refinement and decide how to continue
878 : : //! \details First the mesh refiner object is called to perform a single step
879 : : //! of mesh refinement. Then, if this function is called during a step
880 : : //! (potentially multiple levels of) initial AMR, it evaluates whether to do
881 : : //! another one. If it is called during time stepping, this concludes the
882 : : //! single mesh refinement step and the new mesh is sent to the PDE worker
883 : : //! (Discretization).
884 : : // *****************************************************************************
885 : : {
886 : : // Save old tets and their ids before performing refinement. Outref is always
887 : : // followed by outderef, so to the outside world, the mesh is uchanged, thus
888 : : // no update.
889 [ + + ]: 383 : if (m_mode != RefMode::OUTREF && m_mode != RefMode::OUTDEREF) {
890 : 317 : m_oldntets = m_oldTets.size();
891 : : m_oldTets.clear();
892 [ + + ]: 278509 : for (const auto& [ id, tet ] : m_refiner.tet_store.tets)
893 : : m_oldTets.insert( tet );
894 : : }
895 : :
896 [ + + ]: 383 : if (m_mode == RefMode::T0REF) {
897 : :
898 : : // Refine mesh based on next initial refinement type
899 [ + - ]: 133 : if (!m_initref.empty()) {
900 : 133 : auto r = m_initref.back(); // consume (reversed) list from its back
901 [ + + ]: 133 : if (r == ctr::AMRInitialType::UNIFORM_DEREFINE)
902 : 39 : m_refiner.perform_derefinement();
903 : : else
904 : 94 : m_refiner.perform_refinement();
905 : : }
906 : :
907 : : } else {
908 : :
909 : : // TODO: does not work yet, fix as above
910 : 250 : m_refiner.perform_refinement();
911 : 250 : m_refiner.perform_derefinement();
912 : : }
913 : :
914 : : //auto& tet_store = m_refiner.tet_store;
915 : : //std::cout << "before ref: " << tet_store.marked_refinements.size() << ", " << tet_store.marked_derefinements.size() << ", " << tet_store.size() << ", " << tet_store.get_active_inpoel().size() << '\n';
916 : : //std::cout << "after ref: " << tet_store.marked_refinements.size() << ", " << tet_store.marked_derefinements.size() << ", " << tet_store.size() << ", " << tet_store.get_active_inpoel().size() << '\n';
917 : : //std::cout << "after deref: " << tet_store.marked_refinements.size() << ", " << tet_store.marked_derefinements.size() << ", " << tet_store.size() << ", " << tet_store.get_active_inpoel().size() << '\n';
918 : :
919 : : // Update volume and boundary mesh
920 : 383 : updateMesh();
921 : :
922 : : // Save mesh at every initial refinement step (mainly for debugging). Will
923 : : // replace with just a 'next()' in production.
924 [ + + ]: 383 : if (m_mode == RefMode::T0REF) {
925 : :
926 : 133 : auto l = m_ninitref - m_initref.size() + 1; // num initref steps completed
927 : 133 : auto t0 = g_inputdeck.get< tag::discr, tag::t0 >();
928 : : // Generate times equally subdividing t0-1...t0 to ninitref steps
929 : 133 : auto t =
930 : 133 : t0 - 1.0 + static_cast<tk::real>(l)/static_cast<tk::real>(m_ninitref);
931 : : auto itr = l;
932 : : // Output mesh after refinement step
933 [ + - ][ + - ]: 266 : writeMesh( "t0ref", itr, t,
934 [ + - ][ - + ]: 399 : CkCallback( CkIndex_Refiner::next(), thisProxy[thisIndex] ) );
[ - - ]
935 : :
936 : : } else {
937 : :
938 : 250 : next();
939 : :
940 : : }
941 : 383 : }
942 : :
943 : : void
944 : 383 : Refiner::next()
945 : : // *****************************************************************************
946 : : // Continue after finishing a refinement step
947 : : // *****************************************************************************
948 : : {
949 [ + + ]: 383 : if (m_mode == RefMode::T0REF) {
950 : :
951 : : // Remove initial mesh refinement step from list
952 [ + - ]: 133 : if (!m_initref.empty()) m_initref.pop_back();
953 : : // Continue to next initial AMR step or finish
954 [ + + ]: 133 : if (!m_initref.empty()) t0ref(); else endt0ref();
955 : :
956 [ + + ]: 250 : } else if (m_mode == RefMode::DTREF) {
957 : :
958 : : // Augment node communication map with newly added nodes on chare-boundary
959 [ + + ]: 1332 : for (const auto& [ neighborchare, edges ] : m_remoteEdges) {
960 : : auto& nodes = tk::ref_find( m_nodeCommMap, neighborchare );
961 [ + + ]: 236449 : for (const auto& e : edges) {
962 : : // If parent nodes were part of the node communication map for chare
963 [ + + ][ + + ]: 275995 : if (nodes.find(e[0]) != end(nodes) && nodes.find(e[1]) != end(nodes)) {
964 : : // Add new node if local id was generated for it
965 : 23147 : auto n = Hash<2>()( e );
966 [ + + ]: 46294 : if (m_lid.find(n) != end(m_lid)) nodes.insert( n );
967 : : }
968 : : }
969 : : }
970 : :
971 : : // Send new mesh, solution, and communication data back to PDE worker
972 : 368 : m_scheme[m_meshid].ckLocal< Scheme::resizePostAMR >( thisIndex, m_ginpoel,
973 : 184 : m_el, m_coord, m_addedNodes, m_addedTets, m_removedNodes, m_nodeCommMap,
974 : 184 : m_bface, m_bnode, m_triinpoel );
975 : :
976 [ + + ]: 66 : } else if (m_mode == RefMode::OUTREF) {
977 : :
978 : : // Augment node communication map with newly added nodes on chare-boundary
979 [ + + ]: 105 : for (const auto& [ neighborchare, edges ] : m_remoteEdges) {
980 : : auto& nodes = tk::ref_find( m_nodeCommMap, neighborchare );
981 [ + + ]: 24192 : for (const auto& e : edges) {
982 : : // If parent nodes were part of the node communication map for chare
983 [ + + ][ + + ]: 28392 : if (nodes.find(e[0]) != end(nodes) && nodes.find(e[1]) != end(nodes)) {
984 : : // Add new node if local id was generated for it
985 : 2418 : auto n = Hash<2>()( e );
986 [ + + ]: 4836 : if (m_lid.find(n) != end(m_lid)) nodes.insert( n );
987 : : }
988 : : }
989 : : }
990 : :
991 : : // Store field output mesh
992 : 33 : m_outref_ginpoel = m_ginpoel;
993 : : m_outref_el = m_el;
994 : : m_outref_coord = m_coord;
995 : : m_outref_addedNodes = m_addedNodes;
996 : : m_outref_addedTets = m_addedTets;
997 : : m_outref_nodeCommMap = m_nodeCommMap;
998 : : m_outref_bface = m_bface;
999 : : m_outref_bnode = m_bnode;
1000 : 33 : m_outref_triinpoel = m_triinpoel;
1001 : :
1002 : : // Derefine mesh to the state previous to field output
1003 [ + - ]: 66 : outref( m_outref_bface, m_outref_bnode, m_outref_triinpoel, m_writeCallback,
1004 : : RefMode::OUTDEREF );
1005 : :
1006 [ + - ]: 33 : } else if (m_mode == RefMode::OUTDEREF) {
1007 : :
1008 : : // Send field output mesh to PDE worker
1009 : 66 : m_scheme[m_meshid].ckLocal< Scheme::extractFieldOutput >( thisIndex,
1010 : 33 : m_outref_ginpoel, m_outref_el, m_outref_coord, m_outref_addedNodes,
1011 : 33 : m_outref_addedTets, m_outref_nodeCommMap, m_outref_bface, m_outref_bnode,
1012 : 33 : m_outref_triinpoel, m_writeCallback );
1013 : :
1014 [ - - ][ - - ]: 0 : } else Throw( "RefMode not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
1015 : 383 : }
1016 : :
1017 : : void
1018 : 2106 : Refiner::endt0ref()
1019 : : // *****************************************************************************
1020 : : // Finish initial mesh refinement
1021 : : //! \details This function is called as after initial mesh refinement has
1022 : : //! finished. If initial mesh reifnement was not configured by the user, this
1023 : : //! is the point where we continue after the constructor, by computing the
1024 : : //! total number of elements across the whole problem.
1025 : : // *****************************************************************************
1026 : : {
1027 : : // create sorter Charm++ chare array elements using dynamic insertion
1028 [ + - ]: 4212 : m_sorter[ thisIndex ].insert( m_meshid, m_host, m_meshwriter, m_cbs, m_scheme,
1029 [ + - ][ - + ]: 4212 : CkCallback(CkIndex_Refiner::reorder(), thisProxy[thisIndex]), m_ginpoel,
[ - + ][ - - ]
[ - - ]
1030 [ + - ]: 2106 : m_coordmap, m_el, m_bface, m_triinpoel, m_bnode, m_nchare );
1031 : :
1032 : : // Compute final number of cells across whole problem
1033 : : std::vector< std::size_t >
1034 : 2106 : meshdata{ m_meshid, m_ginpoel.size()/4, m_coord[0].size() };
1035 [ + - ]: 2106 : contribute( meshdata, CkReduction::sum_ulong, m_cbr.get< tag::refined >() );
1036 : :
1037 : : // // Free up memory if no dtref
1038 : : // if (!g_inputdeck.get< tag::amr, tag::dtref >()) {
1039 : : // tk::destroy( m_ginpoel );
1040 : : // tk::destroy( m_el );
1041 : : // tk::destroy( m_coordmap );
1042 : : // tk::destroy( m_coord );
1043 : : // tk::destroy( m_bface );
1044 : : // tk::destroy( m_bnode );
1045 : : // tk::destroy( m_triinpoel );
1046 : : // tk::destroy( m_initref );
1047 : : // tk::destroy( m_ch );
1048 : : // tk::destroy( m_edgech );
1049 : : // tk::destroy( m_chedge );
1050 : : // tk::destroy( m_localEdgeData );
1051 : : // tk::destroy( m_remoteEdgeData );
1052 : : // tk::destroy( m_remoteEdges );
1053 : : // tk::destroy( m_intermediates );
1054 : : // tk::destroy( m_nodeCommMap );
1055 : : // tk::destroy( m_oldTets );
1056 : : // tk::destroy( m_addedNodes );
1057 : : // tk::destroy( m_addedTets );
1058 : : // tk::destroy( m_coarseBndFaces );
1059 : : // tk::destroy( m_coarseBndNodes );
1060 : : // tk::destroy( m_rid );
1061 : : // tk::destroy( m_oldrid );
1062 : : // tk::destroy( m_lref );
1063 : : // tk::destroy( m_parent );
1064 : : // }
1065 : 2106 : }
1066 : :
1067 : : void
1068 : 287 : Refiner::uniformRefine()
1069 : : // *****************************************************************************
1070 : : // Do uniform mesh refinement
1071 : : // *****************************************************************************
1072 : : {
1073 : : // Do uniform refinement
1074 : 287 : m_refiner.mark_uniform_refinement();
1075 : :
1076 : : // Update our extra-edge store based on refiner
1077 : 287 : updateEdgeData();
1078 : :
1079 : : // Set number of extra edges to be zero, skipping correction (if this is the
1080 : : // only step in this refinement step)
1081 : 287 : m_extra = 0;
1082 : 287 : }
1083 : :
1084 : : void
1085 : 72 : Refiner::uniformDeRefine()
1086 : : // *****************************************************************************
1087 : : // Do uniform mesh derefinement
1088 : : // *****************************************************************************
1089 : : {
1090 : : // Do uniform derefinement
1091 : 72 : m_refiner.mark_uniform_derefinement();
1092 : :
1093 : : // Update our extra-edge store based on refiner
1094 : 72 : updateEdgeData();
1095 : :
1096 : : // Set number of extra edges to be zero, skipping correction (if this is the
1097 : : // only step in this refinement step)
1098 : 72 : m_extra = 0;
1099 : 72 : }
1100 : :
1101 : : Refiner::EdgeError
1102 : 2442 : Refiner::errorsInEdges(
1103 : : std::size_t npoin,
1104 : : const std::pair< std::vector<std::size_t>, std::vector<std::size_t> >& esup,
1105 : : const tk::Fields& u ) const
1106 : : // *****************************************************************************
1107 : : // Compute errors in edges
1108 : : //! \param[in] npoin Number nodes in current mesh (partition)
1109 : : //! \param[in] esup Elements surrounding points linked vectors
1110 : : //! \param[in] u Solution evaluated at mesh nodes for all scalar components
1111 : : //! \return A map associating errors (real values between 0.0 and 1.0 incusive)
1112 : : //! to edges (2 local node IDs)
1113 : : // *****************************************************************************
1114 : : {
1115 : : // Get the indices (in the system of systems) of refinement variables and the
1116 : : // error indicator configured
1117 : : const auto& refidx = g_inputdeck.get< tag::amr, tag::id >();
1118 : 2442 : auto errtype = g_inputdeck.get< tag::amr, tag::error >();
1119 : :
1120 : : // Compute points surrounding points
1121 : 4884 : auto psup = tk::genPsup( m_inpoel, 4, esup );
1122 : :
1123 : : // Compute errors in ICs and define refinement criteria for edges
1124 : : AMR::Error error;
1125 : : EdgeError edgeError;
1126 : :
1127 [ + + ]: 641479 : for (std::size_t p=0; p<npoin; ++p) { // for all mesh nodes on this chare
1128 [ + + ][ + + ]: 7533227 : for (auto q : tk::Around(psup,p)) { // for all nodes surrounding p
1129 : : tk::real cmax = 0.0;
1130 : 6894190 : AMR::edge_t e(p,q);
1131 [ + + ]: 8996028 : for (auto i : refidx) { // for all refinement variables
1132 [ + - ]: 2101838 : auto c = error.scalar( u, e, i, m_coord, m_inpoel, esup, errtype );
1133 [ + + ]: 2101838 : if (c > cmax) cmax = c; // find max error at edge
1134 : : }
1135 [ + - ]: 6894190 : edgeError[ {{p,q}} ] = cmax; // associate error to edge
1136 : : }
1137 : : }
1138 : :
1139 : 2442 : return edgeError;
1140 : : }
1141 : :
1142 : : tk::Fields
1143 [ + + ]: 2442 : Refiner::solution( std::size_t npoin,
1144 : : const std::pair< std::vector< std::size_t >,
1145 : : std::vector< std::size_t > >& esup ) const
1146 : : // *****************************************************************************
1147 : : // Update (or evaluate) solution on current mesh
1148 : : //! \param[in] npoin Number nodes in current mesh (partition)
1149 : : //! \param[in] esup Elements surrounding points linked vectors
1150 : : //! \return Solution updated/evaluated for all scalar components
1151 : : // *****************************************************************************
1152 : : {
1153 : : // Get solution whose error to evaluate
1154 : : tk::Fields u;
1155 : :
1156 [ + + ]: 2442 : if (m_mode == RefMode::T0REF) {
1157 : :
1158 : : // Evaluate initial conditions at mesh nodes
1159 [ + - ]: 4428 : u = nodeinit( npoin, esup );
1160 : :
1161 [ + - ]: 228 : } else if (m_mode == RefMode::DTREF) {
1162 : :
1163 : : // Query current solution
1164 [ + - ]: 456 : u = m_scheme[m_meshid].ckLocal< Scheme::solution >( thisIndex );
1165 : :
1166 : 228 : const auto scheme = g_inputdeck.get< tag::discr, tag::scheme >();
1167 [ + - ][ + - ]: 456 : const auto centering = ctr::Scheme().centering( scheme );
[ - - ]
1168 : : if (centering == tk::Centering::ELEM) {
1169 : :
1170 : : // ...
1171 : :
1172 : : }
1173 : :
1174 [ - - ]: 0 : } else if (m_mode == RefMode::OUTREF) {
1175 : :
1176 : :
1177 : :
1178 [ - - ]: 0 : } else if (m_mode == RefMode::OUTDEREF) {
1179 : :
1180 : :
1181 : :
1182 [ - - ][ - - ]: 0 : } else Throw( "RefMode not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
1183 : :
1184 : 2442 : return u;
1185 : : }
1186 : :
1187 : : void
1188 : 22 : Refiner::errorRefine()
1189 : : // *****************************************************************************
1190 : : // Do error-based mesh refinement and derefinement
1191 : : // *****************************************************************************
1192 : : {
1193 : : // Find number of nodes in old mesh
1194 : 22 : auto npoin = tk::npoin_in_graph( m_inpoel );
1195 : : // Generate edges surrounding points in old mesh
1196 : 44 : auto esup = tk::genEsup( m_inpoel, 4 );
1197 : :
1198 : : // Update solution on current mesh
1199 [ + - ]: 22 : auto u = solution( npoin, esup );
1200 : : Assert( u.nunk() == npoin, "Solution uninitialized or wrong size" );
1201 : :
1202 : : using AMR::edge_t;
1203 : : using AMR::edge_tag;
1204 : :
1205 : : // Compute error in edges. Tag edge for refinement if error exceeds
1206 : : // refinement tolerance, tag edge for derefinement if error is below
1207 : : // derefinement tolerance.
1208 : 22 : auto tolref = g_inputdeck.get< tag::amr, tag::tolref >();
1209 [ + - ]: 22 : auto tolderef = g_inputdeck.get< tag::amr, tag::tolderef >();
1210 : : std::vector< std::pair< edge_t, edge_tag > > tagged_edges;
1211 [ + - ][ + + ]: 7333 : for (const auto& e : errorsInEdges(npoin,esup,u)) {
1212 [ + + ]: 7289 : if (e.second > tolref) {
1213 [ + + ][ + - ]: 8092 : tagged_edges.push_back( { edge_t( m_rid[e.first[0]], m_rid[e.first[1]] ),
1214 : : edge_tag::REFINE } );
1215 [ + + ]: 3243 : } else if (e.second < tolderef) {
1216 [ + + ][ + - ]: 5976 : tagged_edges.push_back( { edge_t( m_rid[e.first[0]], m_rid[e.first[1]] ),
1217 : : edge_tag::DEREFINE } );
1218 : : }
1219 : : }
1220 : :
1221 : : // Do error-based refinement
1222 [ + - ]: 22 : m_refiner.mark_error_refinement( tagged_edges );
1223 : :
1224 : : // Update our extra-edge store based on refiner
1225 [ + - ]: 22 : updateEdgeData();
1226 : :
1227 : : // Set number of extra edges to a nonzero number, triggering correction
1228 [ + - ]: 22 : m_extra = 1;
1229 : 22 : }
1230 : :
1231 : : void
1232 [ - - ]: 0 : Refiner::edgelistRefine()
1233 : : // *****************************************************************************
1234 : : // Do mesh refinement based on user explicitly tagging edges
1235 : : // *****************************************************************************
1236 : : {
1237 : : // Get user-defined node-pairs (edges) to tag for refinement
1238 : : const auto& edgenodelist = g_inputdeck.get< tag::amr, tag::edge >();
1239 : :
1240 [ - - ]: 0 : if (!edgenodelist.empty()) { // if user explicitly tagged edges
1241 : : // Find number of nodes in old mesh
1242 : 0 : auto npoin = tk::npoin_in_graph( m_inpoel );
1243 : : // Generate edges surrounding points in old mesh
1244 : 0 : auto esup = tk::genEsup( m_inpoel, 4 );
1245 [ - - ]: 0 : auto psup = tk::genPsup( m_inpoel, 4, esup );
1246 : :
1247 : : EdgeSet useredges;
1248 [ - - ]: 0 : for (std::size_t i=0; i<edgenodelist.size()/2; ++i)
1249 [ - - ]: 0 : useredges.insert( {{ {edgenodelist[i*2+0], edgenodelist[i*2+1]} }} );
1250 : :
1251 : : using AMR::edge_t;
1252 : : using AMR::edge_tag;
1253 : :
1254 : : // Tag edges the user configured
1255 : : std::vector< std::pair< edge_t, edge_tag > > tagged_edges;
1256 [ - - ]: 0 : for (std::size_t p=0; p<npoin; ++p) // for all mesh nodes on this chare
1257 [ - - ]: 0 : for (auto q : tk::Around(psup,p)) { // for all nodes surrounding p
1258 : 0 : Edge e{{ m_gid[p], m_gid[q] }};
1259 : : auto f = useredges.find(e);
1260 [ - - ]: 0 : if (f != end(useredges)) { // tag edge if on user's list
1261 [ - - ][ - - ]: 0 : tagged_edges.push_back( { edge_t( m_rid[p], m_rid[q] ),
[ - - ]
1262 : : edge_tag::REFINE } );
1263 : : useredges.erase( f );
1264 : : }
1265 : : }
1266 : :
1267 : : // Do error-based refinement
1268 [ - - ]: 0 : m_refiner.mark_error_refinement( tagged_edges );
1269 : :
1270 : : // Update our extra-edge store based on refiner
1271 [ - - ]: 0 : updateEdgeData();
1272 : :
1273 : : // Set number of extra edges to a nonzero number, triggering correction
1274 [ - - ]: 0 : m_extra = 1;
1275 : : }
1276 : 0 : }
1277 : :
1278 : : void
1279 : 2 : Refiner::coordRefine()
1280 : : // *****************************************************************************
1281 : : // Do mesh refinement based on tagging edges based on end-point coordinates
1282 : : // *****************************************************************************
1283 : : {
1284 : : // Get user-defined half-world coordinates
1285 : 2 : auto xminus = g_inputdeck.get< tag::amr, tag::xminus >();
1286 : 2 : auto xplus = g_inputdeck.get< tag::amr, tag::xplus >();
1287 : 2 : auto yminus = g_inputdeck.get< tag::amr, tag::yminus >();
1288 : 2 : auto yplus = g_inputdeck.get< tag::amr, tag::yplus >();
1289 : 2 : auto zminus = g_inputdeck.get< tag::amr, tag::zminus >();
1290 : 2 : auto zplus = g_inputdeck.get< tag::amr, tag::zplus >();
1291 : :
1292 : : // The default is the largest representable double
1293 : : auto eps =
1294 : : std::numeric_limits< kw::amr_xminus::info::expect::type >::epsilon();
1295 : 2 : auto xminus_default = g_inputdeck_defaults.get< tag::amr, tag::xminus >();
1296 : 2 : auto xplus_default = g_inputdeck_defaults.get< tag::amr, tag::xplus >();
1297 : 2 : auto yminus_default = g_inputdeck_defaults.get< tag::amr, tag::yminus >();
1298 : 2 : auto yplus_default = g_inputdeck_defaults.get< tag::amr, tag::yplus >();
1299 : 2 : auto zminus_default = g_inputdeck_defaults.get< tag::amr, tag::zminus >();
1300 : 2 : auto zplus_default = g_inputdeck_defaults.get< tag::amr, tag::zplus >();
1301 : :
1302 : : // Decide if user has configured the half-world
1303 [ - + ]: 2 : bool xm = std::abs(xminus - xminus_default) > eps ? true : false;
1304 : 2 : bool xp = std::abs(xplus - xplus_default) > eps ? true : false;
1305 : 2 : bool ym = std::abs(yminus - yminus_default) > eps ? true : false;
1306 : 2 : bool yp = std::abs(yplus - yplus_default) > eps ? true : false;
1307 : 2 : bool zm = std::abs(zminus - zminus_default) > eps ? true : false;
1308 : 2 : bool zp = std::abs(zplus - zplus_default) > eps ? true : false;
1309 : :
1310 : : using AMR::edge_t;
1311 : : using AMR::edge_tag;
1312 : :
1313 [ - + ][ - - ]: 2 : if (xm || xp || ym || yp || zm || zp) { // if any half-world configured
[ - - ]
1314 : : // Find number of nodes in old mesh
1315 : 2 : auto npoin = tk::npoin_in_graph( m_inpoel );
1316 : : // Generate edges surrounding points in old mesh
1317 : 4 : auto esup = tk::genEsup( m_inpoel, 4 );
1318 [ + - ]: 4 : auto psup = tk::genPsup( m_inpoel, 4, esup );
1319 : : // Get access to node coordinates
1320 : : const auto& x = m_coord[0];
1321 : : const auto& y = m_coord[1];
1322 : : const auto& z = m_coord[2];
1323 : : // Compute edges to be tagged for refinement
1324 : : std::vector< std::pair< edge_t, edge_tag > > tagged_edges;
1325 [ + + ]: 621 : for (std::size_t p=0; p<npoin; ++p) // for all mesh nodes on this chare
1326 [ + + ]: 7181 : for (auto q : tk::Around(psup,p)) { // for all nodes surrounding p
1327 : : Edge e{{p,q}};
1328 : :
1329 : : bool t = true;
1330 [ + - ][ + + ]: 6562 : if (xm) { if (x[p]>xminus && x[q]>xminus) t = false; }
[ + + ]
1331 [ - + ][ - - ]: 6562 : if (xp) { if (x[p]<xplus && x[q]<xplus) t = false; }
[ - - ]
1332 [ - + ][ - - ]: 6562 : if (ym) { if (y[p]>yminus && y[q]>yminus) t = false; }
[ - - ]
1333 [ - + ][ - - ]: 6562 : if (yp) { if (y[p]<yplus && y[q]<yplus) t = false; }
[ - - ]
1334 [ - + ][ - - ]: 6562 : if (zm) { if (z[p]>zminus && z[q]>zminus) t = false; }
[ - - ]
1335 [ - + ][ - - ]: 6562 : if (zp) { if (z[p]<zplus && z[q]<zplus) t = false; }
[ - - ]
1336 : :
1337 [ + + ]: 6562 : if (t) {
1338 [ + + ][ + - ]: 9596 : tagged_edges.push_back( { edge_t( m_rid[e[0]], m_rid[e[1]] ),
[ - - ]
1339 : : edge_tag::REFINE } );
1340 : : }
1341 : : }
1342 : :
1343 : : // Do error-based refinement
1344 [ + - ]: 2 : m_refiner.mark_error_refinement( tagged_edges );
1345 : :
1346 : : // Update our extra-edge store based on refiner
1347 [ + - ]: 2 : updateEdgeData();
1348 : :
1349 : : // Set number of extra edges to a nonzero number, triggering correction
1350 [ + - ]: 2 : m_extra = 1;
1351 : : }
1352 : 2 : }
1353 : :
1354 : : tk::Fields
1355 : 2214 : Refiner::nodeinit( std::size_t npoin,
1356 : : const std::pair< std::vector< std::size_t >,
1357 : : std::vector< std::size_t > >& esup ) const
1358 : : // *****************************************************************************
1359 : : // Evaluate initial conditions (IC) at mesh nodes
1360 : : //! \param[in] npoin Number points in mesh (on this chare)
1361 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
1362 : : //! \return Initial conditions (evaluated at t0) at nodes
1363 : : // *****************************************************************************
1364 : : {
1365 : 2214 : auto t0 = g_inputdeck.get< tag::discr, tag::t0 >();
1366 : 2214 : auto nprop = g_inputdeck.get< tag::component >().nprop();
1367 : :
1368 : : // Will store nodal ICs
1369 : 2214 : tk::Fields u( m_coord[0].size(), nprop );
1370 : :
1371 : : // Evaluate ICs differently depending on nodal or cell-centered discretization
1372 : 2214 : const auto scheme = g_inputdeck.get< tag::discr, tag::scheme >();
1373 [ + - ][ + - ]: 4428 : const auto centering = ctr::Scheme().centering( scheme );
[ + + ][ - - ]
1374 : : // list of nodes/elements at which box ICs are defined
1375 : 2214 : std::vector< std::unordered_set< std::size_t > > inbox;
1376 : : tk::real V = 1.0;
1377 : :
1378 [ + + ]: 2214 : if (centering == tk::Centering::NODE) {
1379 : :
1380 : : // Evaluate ICs for all scalar components integrated
1381 [ + + ][ + - ]: 4296 : for (auto& eq : g_cgpde) eq.initialize( m_coord, u, t0, V, inbox );
1382 : :
1383 [ + - ]: 66 : } else if (centering == tk::Centering::ELEM) {
1384 : :
1385 [ + - ]: 66 : auto esuel = tk::genEsuelTet( m_inpoel, esup ); // elems surrounding elements
1386 : : // Initialize cell-based unknowns
1387 [ + - ]: 66 : tk::Fields ue( m_inpoel.size()/4, nprop );
1388 : : auto lhs = ue;
1389 [ + - ]: 66 : auto geoElem = tk::genGeoElemTet( m_inpoel, m_coord );
1390 [ + + ]: 132 : for (const auto& eq : g_dgpde)
1391 : : eq.lhs( geoElem, lhs );
1392 [ + + ]: 132 : for (const auto& eq : g_dgpde)
1393 [ + - ]: 66 : eq.initialize( lhs, m_inpoel, m_coord, inbox, ue, t0, esuel.size()/4 );
1394 : :
1395 : : // Transfer initial conditions from cells to nodes
1396 [ + + ]: 34172 : for (std::size_t p=0; p<npoin; ++p) { // for all mesh nodes on this chare
1397 [ + - ][ - - ]: 34106 : std::vector< tk::real > up( nprop, 0.0 );
1398 : : tk::real vol = 0.0;
1399 [ + + ]: 533606 : for (auto e : tk::Around(esup,p)) { // for all cells around node p
1400 : : // compute nodal volume: every element contributes their volume / 4
1401 : 499500 : vol += geoElem(e,0,0) / 4.0;
1402 : : // sum cell value to node weighed by cell volume / 4
1403 [ + + ]: 1209240 : for (std::size_t c=0; c<nprop; ++c)
1404 [ + - ]: 1419480 : up[c] += ue[e][c] * geoElem(e,0,0) / 4.0;
1405 : : }
1406 : : // store nodal value
1407 [ + + ]: 81412 : for (std::size_t c=0; c<nprop; ++c) u(p,c,0) = up[c] / vol;
1408 : : }
1409 : :
1410 [ - - ][ - - ]: 0 : } else Throw( "Scheme centring not handled for nodal initialization" );
[ - - ][ - - ]
[ - - ][ - - ]
1411 : :
1412 : : Assert( u.nunk() == m_coord[0].size(), "Size mismatch" );
1413 : : Assert( u.nprop() == nprop, "Size mismatch" );
1414 : :
1415 : 2214 : return u;
1416 : : }
1417 : :
1418 : : void
1419 : 383 : Refiner::updateMesh()
1420 : : // *****************************************************************************
1421 : : // Update old mesh after refinement
1422 : : // *****************************************************************************
1423 : : {
1424 : : // Get refined mesh connectivity
1425 : 383 : const auto& refinpoel = m_refiner.tet_store.get_active_inpoel();
1426 : : Assert( refinpoel.size()%4 == 0, "Inconsistent refined mesh connectivity" );
1427 : :
1428 : : // Generate unique node lists of old and refined mesh using local ids
1429 : 383 : auto rinpoel = m_inpoel;
1430 [ + - ]: 383 : tk::remap( rinpoel, m_rid );
1431 [ + - ][ - - ]: 383 : std::unordered_set< std::size_t > old( begin(rinpoel), end(rinpoel) );
1432 : 383 : std::unordered_set< std::size_t > ref( begin(refinpoel), end(refinpoel) );
1433 : :
1434 : : // Augment refiner id -> local node id map with newly added nodes
1435 : : std::size_t l = m_lref.size();
1436 [ + + ][ + + ]: 424109 : for (auto r : ref) if (old.find(r) == end(old)) m_lref[r] = l++;
[ + - ]
1437 : :
1438 : : // Get nodal communication map from Discretization worker
1439 : 383 : if ( m_mode == RefMode::DTREF ||
1440 [ + + ]: 383 : m_mode == RefMode::OUTREF ||
1441 : : m_mode == RefMode::OUTDEREF ) {
1442 : : m_nodeCommMap =
1443 [ + - ]: 500 : m_scheme[m_meshid].disc()[thisIndex].ckLocal()->NodeCommMap();
1444 : : }
1445 : :
1446 : : // Update mesh and solution after refinement
1447 [ + - ]: 383 : newVolMesh( old, ref );
1448 [ + - ]: 383 : newBndMesh( ref );
1449 : :
1450 : : // Update mesh connectivity from refiner lib, remapping refiner to local ids
1451 [ + - ][ + - ]: 383 : m_inpoel = m_refiner.tet_store.get_active_inpoel();
1452 [ + - ]: 383 : tk::remap( m_inpoel, m_lref );
1453 : :
1454 : : // Update mesh connectivity with new global node ids
1455 [ + - ]: 383 : m_ginpoel = m_inpoel;
1456 : : Assert( tk::uniquecopy(m_ginpoel).size() == m_coord[0].size(),
1457 : : "Size mismatch" );
1458 [ + - ]: 383 : tk::remap( m_ginpoel, m_gid );
1459 : :
1460 : : // Ensure valid mesh after refinement
1461 : : Assert( tk::positiveJacobians( m_inpoel, m_coord ),
1462 : : "Refined mesh cell Jacobian non-positive" );
1463 : :
1464 : : Assert( tk::conforming( m_inpoel, m_coord, true, m_rid ),
1465 : : "Mesh not conforming after updating mesh after mesh refinement" );
1466 : :
1467 : : // Perform leak test on new mesh
1468 : : Assert( !tk::leakyPartition(
1469 : : tk::genEsuelTet( m_inpoel, tk::genEsup(m_inpoel,4) ),
1470 : : m_inpoel, m_coord ),
1471 : : "Refined mesh partition leaky" );
1472 : 383 : }
1473 : :
1474 : : void
1475 : 383 : Refiner::newVolMesh( const std::unordered_set< std::size_t >& old,
1476 : : const std::unordered_set< std::size_t >& ref )
1477 : : // *****************************************************************************
1478 : : // Compute new volume mesh after mesh refinement
1479 : : //! \param[in] old Unique nodes of the old (unrefined) mesh using local ids
1480 : : //! \param[in] ref Unique nodes of the refined mesh using local ids
1481 : : // *****************************************************************************
1482 : : {
1483 : : const auto& x = m_coord[0];
1484 : : const auto& y = m_coord[1];
1485 : : const auto& z = m_coord[2];
1486 : :
1487 : : // Generate coordinates and ids to newly added nodes after refinement
1488 : : std::unordered_map< std::size_t, std::size_t > gid_add;
1489 : 383 : tk::destroy( m_addedNodes );
1490 : 383 : tk::destroy( m_removedNodes );
1491 [ + + ]: 212246 : for (auto r : ref) { // for all unique nodes of the refined mesh
1492 [ + + ]: 211863 : if (old.find(r) == end(old)) { // if node is newly added
1493 : : // get (local) parent ids of newly added node
1494 : 163383 : auto p = m_refiner.node_connectivity.get( r );
1495 : : Assert(p[0] != p[1], "Node without parent edge in newVolMesh");
1496 : : Assert( old.find(p[0]) != end(old) && old.find(p[1]) != end(old),
1497 : : "Parent(s) not in old mesh" );
1498 : : Assert( r >= old.size(), "Attempting to overwrite node with added one" );
1499 : : // local parent ids
1500 : 163383 : decltype(p) lp{{tk::cref_find(m_lref,p[0]), tk::cref_find(m_lref,p[1])}};
1501 : : // global parent ids
1502 : 163383 : decltype(p) gp{{m_gid[lp[0]], m_gid[lp[1]]}};
1503 : : // generate new global ID for newly added node
1504 : 163383 : auto g = Hash<2>()( gp );
1505 : :
1506 : : // if node added by AMR lib has not yet been added to Refiner's new mesh
1507 [ + - ]: 163383 : if (m_coordmap.find(g) == end(m_coordmap)) {
1508 : : Assert( g >= old.size(), "Hashed id overwriting old id" );
1509 : : Assert( m_lid.find(g) == end(m_lid),
1510 : : "Overwriting entry global->local node ID map" );
1511 [ + - ]: 163383 : auto l = tk::cref_find( m_lref, r );
1512 : : // store newly added node id and their parent ids (local ids)
1513 : 163383 : m_addedNodes[r] = lp; // key = r for later update to local
1514 : : // assign new node to refiner->global map
1515 [ + - ]: 163383 : gid_add[r] = g; // key = r for later search
1516 : : // assign new node to global->local map
1517 [ + - ]: 163383 : m_lid[g] = l;
1518 : : // generate and store coordinates for newly added node
1519 [ + - ]: 163383 : m_coordmap.insert( {g, {{ (x[lp[0]] + x[lp[1]])/2.0,
1520 : 163383 : (y[lp[0]] + y[lp[1]])/2.0,
1521 [ + - ]: 163383 : (z[lp[0]] + z[lp[1]])/2.0 }} } );
1522 : : }
1523 : : }
1524 : : }
1525 : 383 : tk::destroy( m_coord );
1526 : :
1527 : : // Remove coordinates and ids of removed nodes due to derefinement
1528 : : std::unordered_map< std::size_t, std::size_t > gid_rem;
1529 [ + + ]: 113625 : for (auto o : old) { // for all unique nodes of the old mesh
1530 [ + + ]: 113242 : if (ref.find(o) == end(ref)) { // if node is no longer in new mesh
1531 : 64762 : auto l = tk::cref_find( m_lref, o );
1532 [ + - ]: 64762 : auto g = m_gid[l];
1533 : : // store global-id to local-id map of removed nodes
1534 : : m_removedNodes.insert(l);
1535 [ + - ]: 64762 : gid_rem[l] = g;
1536 : 64762 : m_lid.erase( g );
1537 : : m_coordmap.erase( g );
1538 : : }
1539 : : }
1540 : :
1541 : : // Save previous states of refiner-local node id maps before update
1542 [ + - ]: 383 : m_oldrid = m_rid;
1543 : : //m_oldlref = m_lref;
1544 : :
1545 : : // Generate new node id maps for nodes kept
1546 : 383 : tk::destroy( m_lref );
1547 [ + - ][ + - ]: 383 : std::vector< std::size_t > rid( ref.size() );
1548 [ + - ][ - - ]: 383 : std::vector< std::size_t > gid( ref.size() );
1549 : 383 : std::size_t l = 0; // will generate new local node id
1550 [ + + ]: 113625 : for (std::size_t i=0; i<m_gid.size(); ++i) {
1551 [ + + ]: 113242 : if (gid_rem.find(i) == end(gid_rem)) {
1552 : 48480 : gid[l] = m_gid[i];
1553 : 48480 : rid[l] = m_rid[i];
1554 [ + - ]: 48480 : m_lref[ m_rid[i] ] = l;
1555 : 48480 : ++l;
1556 : : }
1557 : : }
1558 : : // Add newly added nodes due to refinement to node id maps
1559 [ - - ]: 383 : decltype(m_addedNodes) addedNodes( m_addedNodes.size() );
1560 [ + + ]: 163766 : for (const auto& n : gid_add) {
1561 : 163383 : auto r = n.first;
1562 : 163383 : auto g = n.second;
1563 [ + - ]: 163383 : gid[l] = g;
1564 : 163383 : rid[l] = r;
1565 : : Assert(m_lref.find(r) == m_lref.end(), "Overwriting lref");
1566 [ + - ]: 163383 : m_lref[r] = l;
1567 : : auto it = m_addedNodes.find( r );
1568 : : Assert( it != end(m_addedNodes), "Cannot find added node" );
1569 : 163383 : addedNodes[l] = std::move(it->second);
1570 : 163383 : ++l;
1571 : : }
1572 : : Assert( m_lref.size() == ref.size(), "Size mismatch" );
1573 : : m_rid = std::move( rid );
1574 : : Assert( m_rid.size() == ref.size(), "Size mismatch" );
1575 : : m_addedNodes = std::move( addedNodes );
1576 : :
1577 : : // Update node coordinates, ids, and id maps
1578 : : auto& rx = m_coord[0];
1579 : : auto& ry = m_coord[1];
1580 : : auto& rz = m_coord[2];
1581 [ + - ]: 383 : rx.resize( ref.size() );
1582 [ + - ]: 383 : ry.resize( ref.size() );
1583 [ + - ]: 383 : rz.resize( ref.size() );
1584 [ + + ]: 212246 : for (std::size_t i=0; i<gid.size(); ++i) {
1585 : 423726 : tk::ref_find( m_lid, gid[i] ) = i;
1586 : : const auto& c = tk::cref_find( m_coordmap, gid[i] );
1587 : 211863 : rx[i] = c[0];
1588 : 211863 : ry[i] = c[1];
1589 : 211863 : rz[i] = c[2];
1590 : : }
1591 [ + - ]: 383 : m_gid = std::move( gid );
1592 : : Assert( m_gid.size() == m_lid.size() && m_gid.size() == ref.size(),
1593 : : "Size mismatch" );
1594 : 383 : }
1595 : :
1596 : : std::unordered_set< std::size_t >
1597 : 16235820 : Refiner::ancestors( std::size_t n )
1598 : : // *****************************************************************************
1599 : : // Find the oldest parents of a mesh node in the AMR hierarchy
1600 : : //! \param[in] n Local node id whose ancestors to search
1601 : : //! \return Parents of local node id from the coarsest (original) mesh
1602 : : // *****************************************************************************
1603 : : {
1604 : 16235820 : auto d = m_refiner.node_connectivity.get( m_rid[n] );
1605 : 16235820 : decltype(d) p{{ tk::cref_find( m_lref, d[0] ),
1606 [ + + ]: 32471640 : tk::cref_find( m_lref, d[1] ) }};
1607 : :
1608 : : std::unordered_set< std::size_t > s;
1609 : :
1610 [ + + ]: 16235820 : if (p != AMR::node_pair_t{{n,n}}) {
1611 [ + - ]: 5402339 : auto q = ancestors( p[0] );
1612 : : s.insert( begin(q), end(q) );
1613 [ + - ]: 5402339 : auto r = ancestors( p[1] );
1614 : : s.insert( begin(r), end(r) );
1615 : : } else {
1616 : : s.insert( begin(p), end(p) );
1617 : : }
1618 : :
1619 : 16235820 : return s;
1620 : : }
1621 : :
1622 : : Refiner::BndFaceData
1623 : 383 : Refiner::boundary()
1624 : : // *****************************************************************************
1625 : : // Generate boundary data structures used to update refined/derefined boundary
1626 : : // faces and nodes of side sets
1627 : : //! \return A tuple of boundary face data
1628 : : //! \details The output of this function is used to regenerate physical boundary
1629 : : //! face and node data structures after refinement, see updateBndFaces() and
1630 : : //! updateBndNodes().
1631 : : // *****************************************************************************
1632 : : {
1633 : : // Generate the inverse of AMR's tet store
1634 : : std::unordered_map< Tet, std::size_t, Hash<4>, Eq<4> > invtets;
1635 [ + + ]: 943933 : for (const auto& [key, tet] : m_refiner.tet_store.tets)
1636 [ + - ]: 943550 : invtets[ tet ] = key;
1637 : :
1638 : : //std::cout << thisIndex << " invt: " << invtets.size() << '\n';
1639 : : //std::cout << thisIndex << " active inpoel size: " << m_refiner.tet_store.get_active_inpoel().size() << '\n';
1640 : : //std::cout << thisIndex << " marked derefinement size: " << m_refiner.tet_store.marked_derefinements.size() << '\n';
1641 : :
1642 : : // Generate data structure pcReFaceTets, that associates the id of a tet
1643 : : // adjacent to a refined boundary triangle face for all (physical and chare)
1644 : : // boundary faces in the old mesh (i.e., before the current
1645 : : // refinement/derefinement step). Also generate data structure pcDeFaceTets,
1646 : : // that associates the parent tetrahedron (given by four nodes) adjacent to a
1647 : : // derefined boundary face for all (physical and chare) boundary faces in the
1648 : : // old mesh (i.e., before the current refinement/derefinement step).
1649 : : std::unordered_map< Face, std::size_t, Hash<3>, Eq<3> > pcReFaceTets;
1650 : : std::unordered_map< Face, Tet, Hash<3>, Eq<3> > pcDeFaceTets;
1651 [ + - ][ + - ]: 766 : auto oldesuel = tk::genEsuelTet( m_inpoel, tk::genEsup(m_inpoel,4) );
1652 [ + + ]: 426073 : for (std::size_t e=0; e<oldesuel.size()/4; ++e) {
1653 : 425690 : auto m = e*4;
1654 [ + + ]: 2128450 : for (std::size_t f=0; f<4; ++f) {
1655 [ + + ]: 1702760 : if (oldesuel[m+f] == -1) { // if a face does not have an adjacent tet
1656 : 151770 : Face b{{ m_ginpoel[ m+tk::lpofa[f][0] ],
1657 : 151770 : m_ginpoel[ m+tk::lpofa[f][1] ],
1658 : 151770 : m_ginpoel[ m+tk::lpofa[f][2] ] }};
1659 : : Assert( m_inpoel[m+0] < m_oldrid.size() &&
1660 : : m_inpoel[m+1] < m_oldrid.size() &&
1661 : : m_inpoel[m+2] < m_oldrid.size() &&
1662 : : m_inpoel[m+3] < m_oldrid.size(), "Indexing out of rid" );
1663 : 151770 : Tet t{{ m_oldrid[ m_inpoel[m+0] ], m_oldrid[ m_inpoel[m+1] ],
1664 : 151770 : m_oldrid[ m_inpoel[m+2] ], m_oldrid[ m_inpoel[m+3] ] }};
1665 : : //Tet t{{ m_inpoel[m+0], m_inpoel[m+1],
1666 : : // m_inpoel[m+2], m_inpoel[m+3] }};
1667 : : // associate tet id to adjacent (physical or chare) boundary face
1668 : : auto i = invtets.find( t );
1669 [ + + ]: 151770 : if (i != end(invtets)) {
1670 [ + - ]: 61676 : pcReFaceTets[ b ] = i->second;
1671 : : } else {
1672 : : // find parent tet
1673 : 90094 : auto p = tk::cref_find( m_parent, t );
1674 : : // form all 4 faces of parent
1675 : : //auto A = p[0];
1676 : : //auto B = p[1];
1677 : : //auto C = p[2];
1678 : : //auto D = p[3];
1679 : 90094 : auto A = tk::cref_find( m_lref, p[0] );
1680 : 90094 : auto B = tk::cref_find( m_lref, p[1] );
1681 : 90094 : auto C = tk::cref_find( m_lref, p[2] );
1682 [ + - ]: 90094 : auto D = tk::cref_find( m_lref, p[3] );
1683 : : // assign parent tet faces to derefined child's face
1684 : 90094 : pcDeFaceTets[ b ] = {{ A, B, C, D }};
1685 : : }
1686 : : }
1687 : : }
1688 : : }
1689 : :
1690 : : // Generate child->parent tet and id maps after refinement/derefinement step
1691 : 383 : tk::destroy( m_parent );
1692 : : m_addedTets.clear();
1693 : : std::size_t p = 0;
1694 : : std::size_t c = 0;
1695 : : const auto& tet_store = m_refiner.tet_store;
1696 [ + + ]: 943933 : for (const auto& t : tet_store.tets) {
1697 : : // query number of children of tet
1698 [ + - ]: 943550 : auto nc = tet_store.data( t.first ).children.size();
1699 [ + + ]: 1805074 : for (decltype(nc) i=0; i<nc; ++i ) { // for all child tets
1700 : : // get child tet id
1701 [ + - ]: 861524 : auto childtet = tet_store.get_child_id( t.first, i );
1702 : : auto ct = tet_store.tets.find( childtet );
1703 : : //auto cA = tk::cref_find( m_lref, ct->second[0] );
1704 : : //auto cB = tk::cref_find( m_lref, ct->second[1] );
1705 : : //auto cC = tk::cref_find( m_lref, ct->second[2] );
1706 : : //auto cD = tk::cref_find( m_lref, ct->second[3] );
1707 : : // get nodes of parent tet
1708 : : //auto pA = tk::cref_find( m_lref, t.second[0] );
1709 : : //auto pB = tk::cref_find( m_lref, t.second[1] );
1710 : : //auto pC = tk::cref_find( m_lref, t.second[2] );
1711 : : //auto pD = tk::cref_find( m_lref, t.second[3] );
1712 : : // assign parent tet to child tet
1713 : : //m_parent[ {{cA,cB,cC,cD}} ] = {{pA,pB,pC,pD}};
1714 [ + - ]: 861524 : m_parent[ ct->second ] = t.second; //{{pA,pB,pC,pD}};
1715 [ + + ]: 861524 : if (m_oldTets.find(ct->second) == end(m_oldTets)) {
1716 [ + - ][ - - ]: 792886 : m_addedTets[ c++ ] = p - m_oldntets;
1717 : : }
1718 : : }
1719 : 943550 : ++p;
1720 : : }
1721 : :
1722 : : //std::cout << thisIndex << " added: " << m_addedTets.size() << '\n';
1723 : : //std::cout << thisIndex << " parent: " << m_parent.size() << '\n';
1724 : : //std::cout << thisIndex << " pcret: " << pcReFaceTets.size() << '\n';
1725 : : //std::cout << thisIndex << " pcdet: " << pcDeFaceTets.size() << '\n';
1726 : :
1727 : : // Generate unique set of faces for each side set
1728 : : std::unordered_map< int, FaceSet > bndFaces;
1729 [ + + ]: 1053 : for (const auto& [ setid, faceids ] : m_bface) {
1730 : : auto& faces = bndFaces[ setid ];
1731 [ + + ]: 89490 : for (auto f : faceids) {
1732 [ + - ]: 88820 : faces.insert(
1733 [ + - ]: 88820 : {{{ m_triinpoel[f*3+0], m_triinpoel[f*3+1], m_triinpoel[f*3+2] }}} );
1734 : : }
1735 : : }
1736 : :
1737 : 383 : return BndFaceData{ pcReFaceTets, pcDeFaceTets, bndFaces };
1738 : : }
1739 : :
1740 : : void
1741 : 383 : Refiner::newBndMesh( const std::unordered_set< std::size_t >& ref )
1742 : : // *****************************************************************************
1743 : : // Update boundary data structures after mesh refinement
1744 : : //! \param[in] ref Unique nodes of the refined mesh using local ids
1745 : : // *****************************************************************************
1746 : : {
1747 : : // Generate boundary face data structures used to regenerate boundary face
1748 : : // and node data after mesh refinement
1749 : 383 : auto bnd = boundary();
1750 : :
1751 : : // Regerate boundary faces and nodes after mesh refinement
1752 [ + - ]: 383 : updateBndFaces( ref, bnd );
1753 [ + - ]: 383 : updateBndNodes( ref, bnd );
1754 : 383 : }
1755 : :
1756 : : void
1757 : 383 : Refiner::updateBndFaces(
1758 : : [[maybe_unused]] const std::unordered_set< std::size_t >& ref,
1759 : : const BndFaceData& bnd )
1760 : : // *****************************************************************************
1761 : : // Regenerate boundary faces after mesh refinement step
1762 : : //! \param[in] ref Unique nodes of the refined mesh using local ids
1763 : : //! \param[in] bnd Boundary face data bundle
1764 : : // *****************************************************************************
1765 : : {
1766 : : // storage for boundary faces associated to side-set IDs of the refined mesh
1767 : 383 : tk::destroy( m_bface );
1768 : : // storage for boundary faces-node connectivity of the refined mesh
1769 : : tk::destroy( m_triinpoel );
1770 : :
1771 : : // face id counter
1772 : 383 : std::size_t facecnt = 0;
1773 : : // will collect unique faces added for each side set
1774 : : std::unordered_map< int, FaceSet > bf;
1775 : :
1776 : : // Lambda to associate a boundary face and connectivity to a side set.
1777 : : // Argument 's' is the list of faces (ids) to add the new face to. Argument
1778 : : // 'ss' is the side set id to which the face is added. Argument 'f' is the
1779 : : // triangle face connectivity to add.
1780 : 163992 : auto addBndFace = [&]( std::vector< std::size_t >& s, int ss, const Face& f )
1781 : : {
1782 : : // only add face if it has not yet been aded to this side set
1783 [ + + ]: 327984 : if (bf[ ss ].insert( f ).second) {
1784 : 114710 : s.push_back( facecnt++ );
1785 : 114710 : m_triinpoel.insert( end(m_triinpoel), begin(f), end(f) );
1786 : : }
1787 : 383 : };
1788 : :
1789 : : // Regenerate boundary faces for refined tets along side sets. For all
1790 : : // refined faces associated to side sets, we find the ancestors (parents of
1791 : : // nodes in the original/coarsest mesh) of the nodes comprising the faces of
1792 : : // the tetrahedron adjacent to the refined face.
1793 : : const auto& pcReFaceTets = std::get< 0 >( bnd );
1794 : : const auto& bndFaces = std::get< 2 >( bnd );
1795 : : const auto& tet_store = m_refiner.tet_store;
1796 [ + + ]: 62059 : for (const auto& [ face, tetid ] : pcReFaceTets) {
1797 : : // for all side sets of the face, match children's faces to side sets
1798 [ + - ][ + + ]: 148020 : for (const auto& ss : keys(bndFaces,face)) {
1799 : : // will associate to side set id of old (unrefined) mesh boundary face
1800 [ + - ]: 24668 : auto& faces = m_bface[ ss ];
1801 : : const auto& coarsefaces = tk::cref_find( m_coarseBndFaces, ss );
1802 : : // query number of children of boundary tet adjacent to boundary face
1803 [ + - ][ - + ]: 24668 : auto nc = tet_store.data( tetid ).children.size();
1804 [ - + ]: 24668 : if (nc == 0) { // if boundary tet is not refined, add its boundary face
1805 [ - - ]: 0 : addBndFace( faces, ss, face );
1806 : : } else { // if boundary tet is refined
1807 : : const auto& tets = tet_store.tets;
1808 [ + + ]: 222012 : for (decltype(nc) i=0; i<nc; ++i ) { // for all child tets
1809 : : // get child tet id
1810 [ + - ]: 197344 : auto childtet = tet_store.get_child_id( tetid, i );
1811 : : auto ct = tets.find( childtet );
1812 : : Assert( ct != end(tets), "Child tet not found" );
1813 : : // ensure all nodes of child tet are in refined mesh
1814 : : Assert( std::all_of( begin(ct->second), end(ct->second),
1815 : : [&]( std::size_t n ){ return ref.find(n) != end(ref); } ),
1816 : : "Boundary child tet node id not found in refined mesh" );
1817 : : // get nodes of child tet
1818 : 197344 : auto A = tk::cref_find( m_lref, ct->second[0] );
1819 : 197344 : auto B = tk::cref_find( m_lref, ct->second[1] );
1820 : 197344 : auto C = tk::cref_find( m_lref, ct->second[2] );
1821 : 197344 : auto D = tk::cref_find( m_lref, ct->second[3] );
1822 : : // form all 4 faces of child tet
1823 : 197344 : std::array<Face,4> g{{{{A,C,B}}, {{A,B,D}}, {{A,D,C}}, {{B,C,D}}}};
1824 : : // search all faces of child tet and match them to side sets of the
1825 : : // original (coarsest) mesh
1826 [ + + ]: 986720 : for (const auto& rf : g) { // for all faces of child tet
1827 [ + - ]: 789376 : auto a = ancestors( rf[0] );
1828 [ + - ]: 789376 : auto b = ancestors( rf[1] );
1829 [ + - ]: 789376 : auto c = ancestors( rf[2] );
1830 : : a.insert( begin(b), end(b) );
1831 : : a.insert( begin(c), end(c) );
1832 [ + + ]: 789376 : if (a.size() == 3) {
1833 [ + - ]: 381248 : std::vector< std::size_t > p( begin(a), end(a) );
1834 : 381248 : Face par{{ m_gid[p[0]], m_gid[p[1]], m_gid[p[2]] }};
1835 : : auto it = coarsefaces.find( par );
1836 [ + + ]: 381248 : if (it != end(coarsefaces))
1837 [ + - ][ - - ]: 99696 : addBndFace(faces,ss,{{m_gid[rf[0]],m_gid[rf[1]],m_gid[rf[2]]}});
1838 : : }
1839 : : }
1840 : : }
1841 : : }
1842 : : }
1843 : : }
1844 : :
1845 : : // Regenerate boundary faces for derefined tets along side sets. For all
1846 : : // derefined faces associated to side sets, we find the ancestors (parents of
1847 : : // nodes in the original/coarsest mesh) of the nodes comprising the faces of
1848 : : // the parent tetrahedron (previously) associated to the derefined face.
1849 : : const auto& pcDeFaceTets = std::get< 1 >( bnd );
1850 [ + + ]: 90477 : for (const auto& f : pcDeFaceTets) {
1851 [ + - ][ + + ]: 244340 : for (const auto& ss : keys(bndFaces,f.first)) {
1852 : : // will associate to side set id of old (refined) mesh boundary face
1853 [ + - ]: 64152 : auto& faces = m_bface[ ss ];
1854 : : const auto& coarsefaces = tk::cref_find( m_coarseBndFaces, ss );
1855 : : // form all 4 faces of parent tet
1856 : 64152 : auto A = f.second[0];
1857 : 64152 : auto B = f.second[1];
1858 : 64152 : auto C = f.second[2];
1859 : 64152 : auto D = f.second[3];
1860 : 64152 : std::array<Face,4> parf{{ {{A,C,B}}, {{A,B,D}}, {{A,D,C}}, {{B,C,D}} }};
1861 [ + + ]: 320760 : for (const auto& pf : parf) {
1862 [ + - ]: 256608 : auto a = ancestors( pf[0] );
1863 [ + - ]: 256608 : auto b = ancestors( pf[1] );
1864 [ + - ]: 256608 : auto c = ancestors( pf[2] );
1865 : : a.insert( begin(b), end(b) );
1866 : : a.insert( begin(c), end(c) );
1867 [ + - ]: 256608 : if (a.size() == 3) {
1868 [ + - ]: 256608 : std::vector< std::size_t > p( begin(a), end(a) );
1869 : 256608 : Face par{{ m_gid[p[0]], m_gid[p[1]], m_gid[p[2]] }};
1870 : : auto it = coarsefaces.find( par );
1871 [ + + ]: 256608 : if (it != end(coarsefaces))
1872 [ + - ][ - - ]: 64296 : addBndFace(faces,ss,{{m_gid[pf[0]],m_gid[pf[1]],m_gid[pf[2]]}});
1873 : : }
1874 : : }
1875 : : }
1876 : : }
1877 : :
1878 : : //std::cout << thisIndex << " bf: " << tk::sumvalsize( m_bface ) << '\n';
1879 : :
1880 : : // Perform leak-test on boundary face data just updated (only in DEBUG)
1881 : : Assert( bndIntegral(), "Partial boundary integral" );
1882 : 383 : }
1883 : :
1884 : : bool
1885 : 0 : Refiner::bndIntegral()
1886 : : // *****************************************************************************
1887 : : // Compute partial boundary surface integral and sum across all chares
1888 : : //! \return true so we don't trigger assert in client code
1889 : : //! \details This function computes a partial surface integral over the boundary
1890 : : //! of the faces of this mesh partition then sends its contribution to perform
1891 : : //! the integral acorss the total problem boundary. After the global sum a
1892 : : //! non-zero vector result indicates a leak, e.g., a hole in the boundary
1893 : : //! which indicates an error in the boundary face data structures used to
1894 : : //! compute the partial surface integrals.
1895 : : // *****************************************************************************
1896 : : {
1897 : : const auto& x = m_coord[0];
1898 : : const auto& y = m_coord[1];
1899 : : const auto& z = m_coord[2];
1900 : :
1901 : 0 : std::vector< tk::real > s{{ 0.0, 0.0, 0.0 }};
1902 : :
1903 [ - - ]: 0 : for (const auto& [ setid, faceids ] : m_bface) {
1904 [ - - ]: 0 : for (auto f : faceids) {
1905 : 0 : auto A = tk::cref_find( m_lid, m_triinpoel[f*3+0] );
1906 : 0 : auto B = tk::cref_find( m_lid, m_triinpoel[f*3+1] );
1907 : 0 : auto C = tk::cref_find( m_lid, m_triinpoel[f*3+2] );
1908 : : // Compute geometry data for face
1909 : 0 : auto geoface = tk::geoFaceTri( {{x[A], x[B], x[C]}},
1910 [ - - ]: 0 : {{y[A], y[B], y[C]}},
1911 [ - - ]: 0 : {{z[A], z[B], z[C]}} );
1912 : : // Sum up face area * face unit-normal
1913 : 0 : s[0] += geoface(0,0,0) * geoface(0,1,0);
1914 : 0 : s[1] += geoface(0,0,0) * geoface(0,2,0);
1915 : 0 : s[2] += geoface(0,0,0) * geoface(0,3,0);
1916 : : }
1917 : : }
1918 : :
1919 [ - - ]: 0 : s.push_back( -1.0 ); // negative: no call-back after reduction
1920 [ - - ][ - - ]: 0 : s.push_back( static_cast< tk::real >( m_meshid ) );
1921 : :
1922 : : // Send contribution to host summing partial surface integrals
1923 [ - - ]: 0 : contribute( s, CkReduction::sum_double, m_cbr.get< tag::bndint >() );
1924 : :
1925 : 0 : return true; // don't trigger the assert in client code
1926 : : }
1927 : :
1928 : : void
1929 : 383 : Refiner::updateBndNodes(
1930 : : [[maybe_unused]] const std::unordered_set< std::size_t >& ref,
1931 : : const BndFaceData& bnd )
1932 : : // *****************************************************************************
1933 : : // Update boundary nodes after mesh refinement
1934 : : //! \param[in] ref Unique nodes of the refined mesh using local ids
1935 : : //! \param[in] bnd Boundary face data bundle
1936 : : // *****************************************************************************
1937 : : {
1938 : : // storage for boundary nodes associated to side-set IDs of the refined mesh
1939 : 383 : tk::destroy( m_bnode );
1940 : :
1941 : : // Lambda to search the parents in the coarsest mesh of a mesh node and if
1942 : : // found, add its global id to boundary node lists associated to the side
1943 : : // set(s) of its parents. Argument 'n' is the local id of the mesh node id
1944 : : // whose parents to search.
1945 : 2293190 : auto search = [&]( std::size_t n ){
1946 : 2293190 : auto a = ancestors( n ); // find parents of n in coarse mesh
1947 [ + + ]: 2293190 : if (a.size() == 1) {
1948 : : // node was part of the coarse mesh
1949 [ + - ]: 473476 : auto ss = keys( m_coarseBndNodes, m_gid[*a.cbegin()] );
1950 [ + + ]: 492007 : for (auto s : ss)
1951 [ + - ][ + - ]: 18531 : m_bnode[ s ].push_back( m_gid[n] );
1952 [ + + ]: 1819714 : } else if (a.size() == 2) {
1953 : : // node was added to an edge of a coarse face
1954 [ + - ]: 1417566 : std::vector< std::size_t > p( begin(a), end(a) );
1955 [ + - ]: 1417566 : auto ss1 = keys( m_coarseBndNodes, m_gid[p[0]] );
1956 [ + - ]: 1417566 : auto ss2 = keys( m_coarseBndNodes, m_gid[p[1]] );
1957 [ + + ]: 1594268 : for (auto s : ss1) {
1958 [ + + ]: 176702 : if (ss2.find(s) != end(ss2)) {
1959 [ + - ][ + - ]: 130026 : m_bnode[ s ].push_back( m_gid[n] );
1960 : : }
1961 : : }
1962 [ + + ]: 402148 : } else if (a.size() == 3) {
1963 : : // node was added inside of a coarse face
1964 [ + - ]: 369738 : std::vector< std::size_t > p( begin(a), end(a) );
1965 [ + - ]: 369738 : auto ss1 = keys( m_coarseBndNodes, m_gid[p[0]] );
1966 [ + - ]: 369738 : auto ss2 = keys( m_coarseBndNodes, m_gid[p[1]] );
1967 [ + - ]: 369738 : auto ss3 = keys( m_coarseBndNodes, m_gid[p[2]] );
1968 [ + + ]: 528942 : for (auto s : ss1) {
1969 [ + + ][ + + ]: 274958 : if (ss2.find(s) != end(ss2) && ss3.find(s) != end(ss3)) {
1970 [ + - ][ + - ]: 53480 : m_bnode[ s ].push_back( m_gid[n] );
1971 : : }
1972 : : }
1973 : : }
1974 : 2293190 : };
1975 : :
1976 : : const auto& pcReFaceTets = std::get< 0 >( bnd );
1977 : :
1978 : : // Regenerate boundary node lists for refined tets along side sets. For all
1979 : : // refined faces associated to side sets, we find the ancestors (parents of
1980 : : // nodes in the original/coarsest mesh) of the nodes comprising the nodes of
1981 : : // the tetrahedron adjacent to the refined face.
1982 : : const auto& tet_store = m_refiner.tet_store;
1983 [ + + ]: 62059 : for (const auto& [ face, tetid ] : pcReFaceTets) {
1984 : : // query number of children of boundary tet adjacent to boundary face
1985 [ + + ]: 61676 : auto nc = tet_store.data( tetid ).children.size();
1986 [ + + ]: 61676 : if (nc == 0) {
1987 : : // if boundary tet is not refined, add its boundary nodes to the side
1988 : : // set(s) of their parent (in coarse mesh) nodes
1989 : 1066 : auto t = face;
1990 [ + + ]: 4264 : for (auto& n : t) n = tk::cref_find( m_lid, n );
1991 : 1066 : addBndNodes( t, search );
1992 : : } else { // if boundary tet is refined
1993 : : const auto& tets = tet_store.tets;
1994 [ + + ]: 543014 : for (decltype(nc) i=0; i<nc; ++i ) { // for all child tets
1995 : : // get child tet id
1996 : 482404 : auto childtet = tet_store.get_child_id( tetid, i );
1997 : : auto ct = tets.find( childtet );
1998 : : Assert( ct != end(tets), "Child tet not found" );
1999 : : // ensure all nodes of child tet are in refined mesh
2000 : : Assert( std::all_of( begin(ct->second), end(ct->second),
2001 : : [&]( std::size_t n ){ return ref.find(n) != end(ref); } ),
2002 : : "Boundary child tet node id not found in refined mesh" );
2003 : : // search each child tet of refined boundary tet and add their boundary
2004 : : // nodes to the side set(s) of their parent (in coarse mesh) nodes
2005 : 482404 : auto A = tk::cref_find( m_lref, ct->second[0] );
2006 : 482404 : auto B = tk::cref_find( m_lref, ct->second[1] );
2007 : 482404 : auto C = tk::cref_find( m_lref, ct->second[2] );
2008 : 482404 : auto D = tk::cref_find( m_lref, ct->second[3] );
2009 : 964808 : addBndNodes( Tet{{A,B,C,D}}, search );
2010 : : }
2011 : : }
2012 : : }
2013 : :
2014 : : const auto& pcDeFaceTets = std::get< 1 >( bnd );
2015 : :
2016 : : // Regenerate boundary node lists for derefined tets along side sets. For all
2017 : : // refined faces associated to side sets, we find the ancestors (parents of
2018 : : // nodes in the original/coarsest mesh) of the nodes comprising the nodes of
2019 : : // the parent tetrahedron (previously) associated to the derefined face.
2020 [ + + ]: 90477 : for (const auto& f : pcDeFaceTets) addBndNodes( f.second, search );
2021 : :
2022 : : // Make boundary node IDs unique for each physical boundary (side set)
2023 [ + + ]: 527 : for (auto& s : m_bnode) tk::unique( s.second );
2024 : :
2025 : : //std::cout << thisIndex << " bn: " << tk::sumvalsize( m_bnode ) << '\n';
2026 : 383 : }
2027 : :
2028 : : #include "NoWarning/refiner.def.h"
|
