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