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