refineable_mesh.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26
27
28#ifdef OOMPH_HAS_MPI
29#include "mpi.h"
30#endif
31
32#include <cstdlib>
33#include <stdlib.h>
34#include <limits>
35
36#include "refineable_mesh.h"
37// Include to fill in additional_synchronise_hanging_nodes() function
39
40namespace oomph
41{
42 //========================================================================
43 /// Get refinement pattern of mesh: Consider the hypothetical mesh
44 /// obtained by truncating the refinement of the current mesh to a given level
45 /// (where \c level=0 is the un-refined base mesh). To advance
46 /// to the next refinement level, we need to refine (split) the
47 /// \c to_be_refined[level].size() elements identified by the
48 /// element numbers contained in \c vector to_be_refined[level][...]
49 //========================================================================
51 Vector<Vector<unsigned>>& to_be_refined)
52 {
53 // Extract *all* elements from current (fully refined) mesh.
54 Vector<Tree*> all_tree_nodes_pt;
55 forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
56
57 // Find out maximum refinement level
58 unsigned max_level = 0;
59 unsigned nnodes = all_tree_nodes_pt.size();
60 for (unsigned e = 0; e < nnodes; e++)
61 {
62 unsigned level = all_tree_nodes_pt[e]->level();
63 if (level > max_level) max_level = level;
64 }
65
66 // Assign storage for refinement pattern
67 to_be_refined.clear();
68 to_be_refined.resize(max_level);
69 Vector<unsigned> el_count(max_level);
70
71 // Initialise count of elements that exist in mesh when refinement
72 // has proceeded to this level
73 for (unsigned l = 0; l < max_level; l++)
74 {
75 el_count[l] = 0;
76 }
77
78 // Loop over all levels and extract all elements that exist
79 // in reference mesh when refinement has proceeded to this level
80 for (unsigned l = 0; l < max_level; l++)
81 {
82 // Loop over all elements (tree nodes)
83 for (unsigned e = 0; e < nnodes; e++)
84 {
85 // What level does this element exist on?
86 unsigned level = all_tree_nodes_pt[e]->level();
87
88 // Element is part of the mesh at this refinement level
89 // if it exists at this level OR if it exists at a lower level
90 // and is a leaf
91 if ((level == l) || ((level < l) && (all_tree_nodes_pt[e]->is_leaf())))
92 {
93 // If element exsts at this level and is not a leaf it will
94 // be refined when we move to the next level:
95 if ((level == l) && (!all_tree_nodes_pt[e]->is_leaf()))
96 {
97 // Add element number (in mesh at current refinement level)
98 // to the list of elements that need to be refined
99 to_be_refined[l].push_back(el_count[l]);
100 }
101 // Element exists in this mesh: Add to counter
102 el_count[l]++;
103 }
104 }
105 }
106 }
107
108 //========================================================================
109 /// Extract the elements at a particular refinement level in
110 /// the refinement pattern (used in Mesh::redistribute or whatever it's
111 /// going to be called (RefineableMeshBase::reduce_halo_layers or something)
112 //========================================================================
114 unsigned& refinement_level, Vector<RefineableElement*>& level_elements)
115 {
116 // Extract *all* elements from current (fully refined) mesh.
117 Vector<Tree*> all_tree_nodes_pt;
118 forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
119
120 // Add the element to the vector if its level matches refinement_level
121 unsigned nnodes = all_tree_nodes_pt.size();
122 for (unsigned e = 0; e < nnodes; e++)
123 {
124 unsigned level = all_tree_nodes_pt[e]->level();
125 if (level == refinement_level)
126 {
127 level_elements.push_back(
128 dynamic_cast<RefineableElement*>(all_tree_nodes_pt[e]->object_pt()));
129 }
130 }
131 }
132
133 //========================================================================
134 /// Refine original, unrefined mesh according to specified refinement
135 /// pattern (relative to original, unrefined mesh).
136 //========================================================================
138 Vector<Vector<unsigned>>& to_be_refined)
139 {
140 // Get mesh back to unrefined state
141 unsigned my_max, my_min;
142 get_refinement_levels(my_min, my_max);
143
144 // Max refinement level:
145 unsigned my_max_level = to_be_refined.size();
146
147 unsigned global_max = 0;
148 unsigned global_max_level = 0;
149 Vector<unsigned> data(2, 0);
150 data[0] = my_max;
151 data[1] = my_max_level;
152 Vector<unsigned> global_data(2, 0);
153#ifdef OOMPH_HAS_MPI
154 if (this->is_mesh_distributed())
155 {
156 MPI_Allreduce(&data[0],
157 &global_data[0],
158 2,
159 MPI_UNSIGNED,
160 MPI_MAX,
161 Comm_pt->mpi_comm());
162 global_max = global_data[0];
163 global_max_level = global_data[1];
164 }
165 else
166#endif
167 {
168 global_max = my_max;
169 global_max_level = my_max_level;
170 }
171
172
173 for (unsigned i = 0; i < global_max; i++)
174 {
176 }
177
178 // Do refinement steps in current mesh
179 for (unsigned l = 0; l < global_max_level; l++)
180 {
181 // Loop over elements that need to be refined at this level
182 unsigned n_to_be_refined = 0;
183 if (l < my_max_level) n_to_be_refined = to_be_refined[l].size();
184
185 // Select relevant elements to be refined
186 for (unsigned i = 0; i < n_to_be_refined; i++)
187 {
188 dynamic_cast<RefineableElement*>(this->element_pt(to_be_refined[l][i]))
189 ->select_for_refinement();
190 }
191
192 // Now do the actual mesh refinement
193 adapt_mesh();
194 }
195 }
196
197
198 //========================================================================
199 /// Refine base mesh according to refinement pattern in restart file
200 //========================================================================
201 void TreeBasedRefineableMeshBase::refine(std::ifstream& restart_file)
202 {
203 // Assign storage for refinement pattern
204 Vector<Vector<unsigned>> to_be_refined;
205
206 // Read refinement pattern
207 read_refinement(restart_file, to_be_refined);
208
209 // Refine
210 refine_base_mesh(to_be_refined);
211 }
212
213
214 //========================================================================
215 /// Dump refinement pattern to allow for rebuild
216 ///
217 //========================================================================
219 {
220 // Assign storage for refinement pattern
221 Vector<Vector<unsigned>> to_be_refined;
222
223 // Get refinement pattern of reference mesh:
224 get_refinement_pattern(to_be_refined);
225
226 // Dump max refinement level:
227 unsigned max_level = to_be_refined.size();
228 outfile << max_level << " # max. refinement level " << std::endl;
229
230 // Doc the numbers of the elements that need to be refined at this level
231 for (unsigned l = 0; l < max_level; l++)
232 {
233 // Loop over elements that need to be refined at this level
234 unsigned n_to_be_refined = to_be_refined[l].size();
235 outfile << n_to_be_refined << " # number of elements to be refined. "
236 << "What follows are the numbers of the elements. " << std::endl;
237
238 // Select relevant elements to be refined
239 for (unsigned i = 0; i < n_to_be_refined; i++)
240 {
241 outfile << to_be_refined[l][i] << std::endl;
242 }
243 }
244 }
245
246
247 //========================================================================
248 /// Read refinement pattern to allow for rebuild
249 ///
250 //========================================================================
252 std::ifstream& restart_file, Vector<Vector<unsigned>>& to_be_refined)
253 {
254 std::string input_string;
255
256 // Read max refinement level:
257
258 // Read line up to termination sign
259 getline(restart_file, input_string, '#');
260
261 // Ignore rest of line
262 restart_file.ignore(80, '\n');
263
264 // Convert
265 unsigned max_level = std::atoi(input_string.c_str());
266
267 // Assign storage for refinement pattern
268 to_be_refined.resize(max_level);
269
270 // Read the number of the elements that need to be refined at different
271 // levels
272 for (unsigned l = 0; l < max_level; l++)
273 {
274 // Read line up to termination sign
275 getline(restart_file, input_string, '#');
276
277 // Ignore rest of line
278 restart_file.ignore(80, '\n');
279
280 // Convert
281 unsigned n_to_be_refined = atoi(input_string.c_str());
282 ;
283
284 // Assign storage
285 to_be_refined[l].resize(n_to_be_refined);
286
287 // Read numbers of the elements that need to be refined
288 for (unsigned i = 0; i < n_to_be_refined; i++)
289 {
290 restart_file >> to_be_refined[l][i];
291 }
292 }
293 }
294
295
296 //========================================================================
297 /// Do adaptive refinement for mesh.
298 /// - Pass Vector of error estimates for all elements.
299 /// - Refine those whose errors exceeds the threshold
300 /// - (Try to) unrefine those whose errors is less than
301 /// threshold (only possible if the three brothers also want to be
302 /// unrefined, of course.)
303 /// - Update the nodal positions in the whole lot
304 /// - Store # of refined/unrefined elements.
305 /// - Doc refinement process (if required)
306 //========================================================================
308 {
309 // Set the refinement tolerance to be the max permissible error
310 double refine_tol = max_permitted_error();
311
312 // Set the unrefinement tolerance to be the min permissible error
313 double unrefine_tol = min_permitted_error();
314
315 // Setup doc info
316 DocInfo local_doc_info;
317 if (doc_info_pt() == 0)
318 {
319 local_doc_info.disable_doc();
320 }
321 else
322 {
323 local_doc_info = doc_info();
324 }
325
326
327 // Check that the errors make sense
328 if (refine_tol <= unrefine_tol)
329 {
330 std::ostringstream error_stream;
331 error_stream << "Refinement tolerance <= Unrefinement tolerance"
332 << refine_tol << " " << unrefine_tol << std::endl
333 << "doesn't make sense and will almost certainly crash"
334 << std::endl
335 << "this beautiful code!" << std::endl;
336
337 throw OomphLibError(
338 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
339 }
340
341
342 // Select elements for refinement and unrefinement
343 //================================================
344 // Reset counter for number of elements that would like to be
345 // refined further but can't
347
348 // Note: Yes, this needs to be a map because we'll have to check
349 // the refinement wishes of brothers (who we only access via pointers)
350 std::map<RefineableElement*, bool> wants_to_be_unrefined;
351
352 // Initialise a variable to store the number of elements for refinment
353 unsigned n_refine = 0;
354
355 // Loop over all elements and mark them according to the error criterion
356 unsigned long Nelement = this->nelement();
357 for (unsigned long e = 0; e < Nelement; e++)
358 {
359 // (Cast) pointer to the element
360 RefineableElement* el_pt =
361 dynamic_cast<RefineableElement*>(this->element_pt(e));
362
363 // Initially element is not to be refined
365
366 // If the element error exceeds the threshold ...
367 if (elemental_error[e] > refine_tol)
368 {
369 // ... and its refinement level is less than the maximum desired level
370 // mark is to be refined
371 if ((el_pt->refinement_is_enabled()) &&
373 {
374 el_pt->select_for_refinement();
375 n_refine++;
376 }
377 // ... otherwise mark it as having been over-ruled
378 else
379 {
381 }
382 }
383
384 // Now worry about unrefinement (first pass):
385
386 // Is my error too small AND do I have a father?
387 if ((elemental_error[e] < unrefine_tol) &&
388 (el_pt->tree_pt()->father_pt() != 0))
389 {
390 // Flag to indicate whether to unrefine
391 bool unrefine = true;
392 unsigned n_sons = el_pt->tree_pt()->father_pt()->nsons();
393
394 // Are all brothers leaf nodes?
395 for (unsigned ison = 0; ison < n_sons; ison++)
396 {
397 // (At least) one brother is not a leaf: end of story; we're not doing
398 // it (= the unrefinement)
399 if (!(el_pt->tree_pt()->father_pt()->son_pt(ison)->is_leaf()))
400 {
401 unrefine = false;
402 }
403 }
404
405 // Don't allow unrefinement of elements that would become larger
406 // than the minimum legal refinement level
407 if (el_pt->refinement_level() - 1 < min_refinement_level())
408 {
409 unrefine = false;
410 }
411
412 // So, all things considered, is the element eligbible for refinement?
413 if (unrefine)
414 {
415 wants_to_be_unrefined[el_pt] = true;
416 }
417 else
418 {
419 wants_to_be_unrefined[el_pt] = false;
420 }
421 }
422 }
423
424 oomph_info << " \n Number of elements to be refined: " << n_refine
425 << std::endl;
426 oomph_info << " \n Number of elements whose refinement was overruled: "
427 << nrefinement_overruled() << std::endl;
428
429 // Second pass for unrefinement --- an element cannot be unrefined unless
430 // all brothers want to be unrefined.
431 // Loop over all elements again and let the first set of sons check if their
432 // brothers also want to be unrefined
433 unsigned n_unrefine = 0;
434 for (unsigned long e = 0; e < Nelement; e++)
435 {
436 //(Cast) pointer to the element
437 RefineableElement* el_pt =
438 dynamic_cast<RefineableElement*>(this->element_pt(e));
439
440 // hierher: This is a bit naughty... We want to put the
441 // first son in charge -- the statement below assumes (correctly) that the
442 // enumeration of all (!) trees starts with son types.
443 // This is correct for oc and quadtrees but will bite us if we
444 // ever introduce other trees if/when we accidentally break this
445 // tacit assumption. Not sure what to do about it for
446 // now other than leaving it hierher-ed...
447 if (el_pt->tree_pt()->son_type() == OcTreeNames::LDB)
448 {
449 // Do all sons want to be unrefined?
450 bool unrefine = true;
451 unsigned n_sons = el_pt->tree_pt()->father_pt()->nsons();
452 for (unsigned ison = 0; ison < n_sons; ison++)
453 {
454 if (!(wants_to_be_unrefined[dynamic_cast<RefineableElement*>(
455 el_pt->tree_pt()->father_pt()->son_pt(ison)->object_pt())]))
456 {
457 // One guy isn't cooperating and spoils the party.
458 unrefine = false;
459 }
460 }
461
462 // Tell father that his sons need to be merged
463 if (unrefine)
464 {
465 el_pt->tree_pt()
466 ->father_pt()
467 ->object_pt()
469 n_unrefine += n_sons;
470 }
471 // Otherwise mark the sons as not to be touched
472 else
473 {
474 el_pt->tree_pt()
475 ->father_pt()
476 ->object_pt()
478 }
479 }
480 }
481 oomph_info << " \n Number of elements to be merged : " << n_unrefine
482 << std::endl
483 << std::endl;
484
485
486 // Now do the actual mesh adaptation
487 //---------------------------------
488
489 // Check whether its worth our while
490 // Either some elements want to be refined,
491 // or the number that want to be unrefined are greater than the
492 // specified tolerance
493
494 // In a parallel job, it is possible that one process may not have
495 // any elements to refine, BUT a neighbouring process may refine an
496 // element which changes the hanging status of a node that is on
497 // both processes (i.e. a halo(ed) node). To get around this issue,
498 // ALL processes need to call adapt_mesh if ANY refinement is to
499 // take place anywhere.
500
501 unsigned total_n_refine = 0;
502#ifdef OOMPH_HAS_MPI
503 // Sum n_refine across all processors
504 if (this->is_mesh_distributed())
505 {
506 MPI_Allreduce(&n_refine,
507 &total_n_refine,
508 1,
509 MPI_UNSIGNED,
510 MPI_SUM,
511 Comm_pt->mpi_comm());
512 }
513 else
514 {
515 total_n_refine = n_refine;
516 }
517#else
518 total_n_refine = n_refine;
519#endif
520
521 // There may be some issues with unrefinement too, but I have not
522 // been able to come up with an example (either in my head or in a
523 // particular problem) where anything has arisen. I can see that
524 // there may be an issue if n_unrefine differs across processes so
525 // that (total_n_unrefine > max_keep_unrefined()) on some but not
526 // all processes. I haven't seen any examples of this yet so the
527 // following code may or may not work! (Andy, 06/03/08)
528
529 unsigned total_n_unrefine = 0;
530#ifdef OOMPH_HAS_MPI
531 // Sum n_unrefine across all processors
532 if (this->is_mesh_distributed())
533 {
534 MPI_Allreduce(&n_unrefine,
535 &total_n_unrefine,
536 1,
537 MPI_UNSIGNED,
538 MPI_SUM,
539 Comm_pt->mpi_comm());
540 }
541 else
542 {
543 total_n_unrefine = n_unrefine;
544 }
545#else
546 total_n_unrefine = n_unrefine;
547#endif
548
549 oomph_info << "---> " << total_n_refine << " elements to be refined, and "
550 << total_n_unrefine << " to be unrefined, in total.\n"
551 << std::endl;
552
553 if ((total_n_refine > 0) || (total_n_unrefine > max_keep_unrefined()))
554 {
555#ifdef PARANOID
556#ifdef OOMPH_HAS_MPI
557
558
559 // Sanity check: Each processor checks if the enforced unrefinement of
560 // its haloed element is matched by enforced unrefinement of the
561 // corresponding halo elements on the other processors.
562 if (this->is_mesh_distributed())
563 {
564 // Store number of processors and current process
565 MPI_Status status;
566 int n_proc = Comm_pt->nproc();
567 int my_rank = Comm_pt->my_rank();
568
569 // Loop over processes: Each processor sends unrefinement pattern
570 // for halo elements with processor d to processor d where it's
571 // compared against the unrefinement pattern for the corresponding
572 // haloed elements
573 for (int d = 0; d < n_proc; d++)
574 {
575 // No halo with self: Send unrefinement info to proc d
576 if (d != my_rank)
577 {
578 // Get the vector of halo elements whose non-halo counterpart
579 // are on processor d
580 Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
581
582 // Create vector containing (0)1 to indicate that
583 // halo element is (not) to be unrefined
584 unsigned nhalo = halo_elem_pt.size();
585 Vector<int> halo_to_be_unrefined(nhalo, 0);
586 for (unsigned e = 0; e < nhalo; e++)
587 {
588 if (dynamic_cast<RefineableElement*>(halo_elem_pt[e])
589 ->sons_to_be_unrefined())
590 {
591 halo_to_be_unrefined[e] = 1;
592 }
593 }
594
595 // Trap the case when there are no halo elements
596 // so that we don't get a segfault in the MPI send
597 if (nhalo > 0)
598 {
599 // Send it across to proc d
600 MPI_Send(&halo_to_be_unrefined[0],
601 nhalo,
602 MPI_INT,
603 d,
604 0,
605 Comm_pt->mpi_comm());
606 }
607 }
608 // else (d=my_rank): Receive unrefinement pattern from all
609 // other processors (dd)
610 else
611 {
612 // Loop over other processors
613 for (int dd = 0; dd < n_proc; dd++)
614 {
615 // No halo with yourself
616 if (dd != d)
617 {
618 // Get the vector of haloed elements on current processor
619 // with processor dd
620 Vector<GeneralisedElement*> haloed_elem_pt(
621 this->haloed_element_pt(dd));
622
623 // Ask processor dd to send vector containing (0)1 for
624 // halo element with current processor to be (not)unrefined
625 unsigned nhaloed = haloed_elem_pt.size();
626 Vector<int> halo_to_be_unrefined(nhaloed);
627 // Trap to catch the case that there are no haloed elements
628 if (nhaloed > 0)
629 {
630 // Receive unrefinement pattern of haloes from proc dd
631 MPI_Recv(&halo_to_be_unrefined[0],
632 nhaloed,
633 MPI_INT,
634 dd,
635 0,
636 Comm_pt->mpi_comm(),
637 &status);
638 }
639
640 // Check it
641 for (unsigned e = 0; e < nhaloed; e++)
642 {
643 if (((halo_to_be_unrefined[e] == 0) &&
644 (dynamic_cast<RefineableElement*>(haloed_elem_pt[e])
645 ->sons_to_be_unrefined())) ||
646 ((halo_to_be_unrefined[e] == 1) &&
647 (!dynamic_cast<RefineableElement*>(haloed_elem_pt[e])
648 ->sons_to_be_unrefined())))
649 {
650 std::ostringstream error_message;
651 error_message
652 << "Error in unrefinement: \n"
653 << "Haloed element: " << e << " on proc " << my_rank
654 << " \n"
655 << "wants to be unrefined whereas its halo counterpart "
656 "on\n"
657 << "proc " << dd << " doesn't (or vice versa)...\n"
658 << "This is most likely because the error estimator\n"
659 << "has not assigned the same errors to halo and haloed\n"
660 << "elements -- it ought to!\n";
661 throw OomphLibError(error_message.str(),
662 OOMPH_CURRENT_FUNCTION,
663 OOMPH_EXCEPTION_LOCATION);
664 }
665 }
666 }
667 }
668 }
669 }
670
671
672 // Loop over processes: Each processor sends refinement pattern
673 // for halo elements with processor d to processor d where it's
674 // compared against the refinement pattern for the corresponding
675 // haloed elements
676 for (int d = 0; d < n_proc; d++)
677 {
678 // No halo with self: Send refinement info to proc d
679 if (d != my_rank)
680 {
681 // Get the vector of halo elements whose non-halo counterpart
682 // are on processor d
683 Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
684
685 // Create vector containing (0)1 to indicate that
686 // halo element is (not) to be refined
687 unsigned nhalo = halo_elem_pt.size();
688 Vector<int> halo_to_be_refined(nhalo, 0);
689 for (unsigned e = 0; e < nhalo; e++)
690 {
691 if (dynamic_cast<RefineableElement*>(halo_elem_pt[e])
692 ->to_be_refined())
693 {
694 halo_to_be_refined[e] = 1;
695 }
696 }
697
698 // Trap the case when there are no halo elements
699 // so that we don't get a segfault in the MPI send
700 if (nhalo > 0)
701 {
702 // Send it across to proc d
703 MPI_Send(&halo_to_be_refined[0],
704 nhalo,
705 MPI_INT,
706 d,
707 0,
708 Comm_pt->mpi_comm());
709 }
710 }
711 // else (d=my_rank): Receive refinement pattern from all
712 // other processors (dd)
713 else
714 {
715 // Loop over other processors
716 for (int dd = 0; dd < n_proc; dd++)
717 {
718 // No halo with yourself
719 if (dd != d)
720 {
721 // Get the vector of haloed elements on current processor
722 // with processor dd
723 Vector<GeneralisedElement*> haloed_elem_pt(
724 this->haloed_element_pt(dd));
725
726 // Ask processor dd to send vector containing (0)1 for
727 // halo element with current processor to be (not)refined
728 unsigned nhaloed = haloed_elem_pt.size();
729 Vector<int> halo_to_be_refined(nhaloed);
730 // Trap to catch the case that there are no haloed elements
731 if (nhaloed > 0)
732 {
733 // Receive unrefinement pattern of haloes from proc dd
734 MPI_Recv(&halo_to_be_refined[0],
735 nhaloed,
736 MPI_INT,
737 dd,
738 0,
739 Comm_pt->mpi_comm(),
740 &status);
741 }
742
743 // Check it
744 for (unsigned e = 0; e < nhaloed; e++)
745 {
746 if (((halo_to_be_refined[e] == 0) &&
747 (dynamic_cast<RefineableElement*>(haloed_elem_pt[e])
748 ->to_be_refined())) ||
749 ((halo_to_be_refined[e] == 1) &&
750 (!dynamic_cast<RefineableElement*>(haloed_elem_pt[e])
751 ->to_be_refined())))
752 {
753 std::ostringstream error_message;
754 error_message
755 << "Error in refinement: \n"
756 << "Haloed element: " << e << " on proc " << my_rank
757 << " \n"
758 << "wants to be refined whereas its halo counterpart on\n"
759 << "proc " << dd << " doesn't (or vice versa)...\n"
760 << "This is most likely because the error estimator\n"
761 << "has not assigned the same errors to halo and haloed\n"
762 << "elements -- it ought to!\n";
763 throw OomphLibError(error_message.str(),
764 OOMPH_CURRENT_FUNCTION,
765 OOMPH_EXCEPTION_LOCATION);
766 }
767 }
768 }
769 }
770 }
771 }
772 }
773
774#endif
775#endif
776
777 // Perform the actual adaptation
778 adapt_mesh(local_doc_info);
779
780 // The number of refineable elements is still local to each process
781 Nunrefined = n_unrefine;
782 Nrefined = n_refine;
783 }
784 // If not worthwhile, say so but still reorder nodes and kill
785 // external storage for consistency in parallel computations
786 else
787 {
788#ifdef OOMPH_HAS_MPI
789 // Delete any external element storage - any interaction will still
790 // be set up on the fly again, so we need to get rid of old information.
791 // This particularly causes problems in multi-domain examples where
792 // we decide not to refine one of the meshes
794#endif
795
796 // Reorder the nodes within the mesh's node vector
797 // to establish a standard ordering regardless of the sequence
798 // of mesh refinements -- this is required to allow dump/restart
799 // on refined meshes
800 this->reorder_nodes();
801
802#ifdef OOMPH_HAS_MPI
803
804 // Now (re-)classify halo and haloed nodes and synchronise hanging
805 // nodes
806 // This is required in cases where delete_all_external_storage()
807 // made dependent nodes to external halo nodes nonhanging.
808 if (this->is_mesh_distributed())
809 {
813 }
814
815#endif
816
817 if (n_refine == 0)
818 {
819 oomph_info << " Not enough benefit in adapting mesh." << std::endl
820 << std::endl;
821 }
822 Nunrefined = 0;
823 Nrefined = 0;
824 }
825 }
826
827 //========================================================================
828 /// Get max/min refinement level
829 //========================================================================
831 unsigned& min_refinement_level, unsigned& max_refinement_level)
832 {
833 // Initialise
834 min_refinement_level = UINT_MAX;
836
837 // Loop over all elements
838 unsigned long n_element = this->nelement();
839 if (n_element == 0)
840 {
843 }
844 else
845 {
846 for (unsigned long e = 0; e < n_element; e++)
847 {
848 // Get the refinement level of the element
849 unsigned level = dynamic_cast<RefineableElement*>(this->element_pt(e))
850 ->refinement_level();
851
852 if (level > max_refinement_level) max_refinement_level = level;
853 if (level < min_refinement_level) min_refinement_level = level;
854 }
855 }
856 }
857
858
859 //================================================================
860 /// Adapt mesh, which exists in two representations,
861 /// namely as:
862 /// - a FE mesh
863 /// - a forest of Oc or QuadTrees
864 ///
865 /// Refinement/derefinement process is documented (in tecplot-able form)
866 /// if requested.
867 ///
868 /// Procedure:
869 /// - Loop over all elements and do the refinement for those who want to
870 /// be refined. Note: Refinement/splitting only allocates elements but
871 /// doesn't build them.
872 /// - Build the new elements (i.e. give them nodes (create new ones where
873 /// necessary), assign boundary conditions, and add nodes to mesh
874 /// and mesh boundaries.
875 /// - For all nodes that were hanging on the previous mesh (and are still
876 /// marked as such), fill in their nodal values (consistent
877 /// with the current hanging node scheme) to make sure they are fully
878 /// functional, should they have become non-hanging during the
879 /// mesh-adaptation. Then mark the nodes as non-hanging.
880 /// - Unrefine selected elements (which may cause nodes to be re-built).
881 /// - Add the new elements to the mesh (by completely overwriting
882 /// the old Vector of elements).
883 /// - Delete any nodes that have become obsolete.
884 /// - Mark up hanging nodes and setup hanging node scheme (incl.
885 /// recursive cleanup for hanging nodes that depend on other
886 /// hanging nodes).
887 /// - Adjust position of hanging nodes to make sure their position
888 /// is consistent with the FE-based represenetation of their larger
889 /// neighbours.
890 /// - run a quick self-test on the neighbour finding scheme and
891 /// check the integrity of the elements (if PARANOID)
892 /// - doc hanging node status, boundary conditions, neighbour
893 /// scheme if requested.
894 ///
895 ///
896 /// After adaptation, all nodes (whether new or old) have up-to-date
897 /// current and previous values.
898 ///
899 /// If refinement process is being documented, the following information
900 /// is documented:
901 /// - The files
902 /// - "neighbours.dat"
903 /// - "all_nodes.dat"
904 /// - "new_nodes.dat"
905 /// - "hang_nodes_*.dat"
906 /// where the * denotes a direction (n,s,e,w) in 2D
907 /// or (r,l,u,d,f,b) in 3D
908 /// .
909 /// can be viewed with
910 /// - QHangingNodes.mcr
911 /// .
912 /// - The file
913 /// - "hangnodes_withmasters.dat"
914 /// .
915 /// can be viewed with
916 /// - QHangingNodesWithMasters.mcr
917 /// .
918 /// to check the hanging node status.
919 /// - The neighbour status of the elements is documented in
920 /// - "neighbours.dat"
921 /// .
922 /// and can be viewed with
923 /// - QuadTreeNeighbours.mcr
924 /// .
925 //=================================================================
927 {
928#ifdef OOMPH_HAS_MPI
929 // Delete any external element storage before performing the adaptation
930 // (in particular, external halo nodes that are on mesh boundaries)
932#endif
933
934 // Only perform the adapt step if the mesh has any elements. This is
935 // relevant in a distributed problem with multiple meshes, where a
936 // particular process may not have any elements on a particular submesh.
937 if (this->nelement() > 0)
938 {
939 // Pointer to mesh needs to be passed to some functions
940 Mesh* mesh_pt = this;
941
942 double t_start = 0.0;
944 {
945 t_start = TimingHelpers::timer();
946 }
947
948 // Do refinement(=splitting) of elements that have been selected
949 // This function encapsulates the template parameter
951
952
954 {
955 double t_end = TimingHelpers::timer();
956 oomph_info << "Time for split_elements_if_required: " << t_end - t_start
957 << std::endl;
958 t_start = TimingHelpers::timer();
959 }
960
961 // Now elements have been created -- build all the leaves
962 //-------------------------------------------------------
963 // Firstly put all the elements into a vector
964 Vector<Tree*> leaf_nodes_pt;
965 Forest_pt->stick_leaves_into_vector(leaf_nodes_pt);
966
967
969 {
970 double t_end = TimingHelpers::timer();
971 oomph_info << "Time for stick_leaves_into_vector: " << t_end - t_start
972 << std::endl;
973 t_start = TimingHelpers::timer();
974 }
975
976 // If we are documenting the output, create the filename
977 std::ostringstream fullname;
978 std::ofstream new_nodes_file;
980 {
981 fullname << doc_info.directory() << "/new_nodes" << doc_info.number()
982 << ".dat";
983 new_nodes_file.open(fullname.str().c_str());
984 }
985
986
987 // Build all elements and store vector of pointers to new nodes
988 // (Note: build() checks if the element has been built
989 // already, i.e. if it's not a new element).
990 Vector<Node*> new_node_pt;
991 bool was_already_built;
992 unsigned long num_tree_nodes = leaf_nodes_pt.size();
993 for (unsigned long e = 0; e < num_tree_nodes; e++)
994 {
995 // Pre-build must be performed before any elements are built
996 leaf_nodes_pt[e]->object_pt()->pre_build(mesh_pt, new_node_pt);
997 }
998 for (unsigned long e = 0; e < num_tree_nodes; e++)
999 {
1000 // Now do the actual build of the new elements
1001 leaf_nodes_pt[e]->object_pt()->build(
1002 mesh_pt, new_node_pt, was_already_built, new_nodes_file);
1003 }
1004
1005
1006 double t_end = 0.0;
1008 {
1009 t_end = TimingHelpers::timer();
1010 oomph_info << "Time for building " << num_tree_nodes
1011 << " new elements: " << t_end - t_start << std::endl;
1012 t_start = TimingHelpers::timer();
1013 }
1014
1015 // Close the new nodes files, if it was opened
1017 {
1018 new_nodes_file.close();
1019 }
1020
1021 // Loop over all nodes in mesh and free the dofs of those that were
1022 //-----------------------------------------------------------------
1023 // pinned only because they were hanging nodes. Also update their
1024 //-----------------------------------------------------------------
1025 // nodal values so that they contain data that is consistent
1026 //----------------------------------------------------------
1027 // with the hanging node representation
1028 //-------------------------------------
1029 // (Even if the nodal data isn't actually accessed because the node
1030 // is still hanging -- we don't know this yet, and this step makes
1031 // sure that all nodes are fully functional and up-to-date, should
1032 // they become non-hanging below).
1033 //
1034 //
1035 // However, if we have a fixed mesh and hanging nodes on the boundary
1036 // become non-hanging they will not necessarily respect the curvilinear
1037 // boundaries. This can only happen in 3D of course because it is not
1038 // possible to have hanging nodes on boundaries in 2D.
1039 // The solution is to store those nodes on the boundaries that are
1040 // currently hanging and then check to see whether they have changed
1041 // status at the end of the refinement procedure.
1042 // If it has changed, then we need to adjust their positions.
1043 const unsigned n_boundary = this->nboundary();
1044 const unsigned mesh_dim = this->finite_element_pt(0)->dim();
1045 Vector<std::set<Node*>> hanging_nodes_on_boundary_pt(n_boundary);
1046
1047 unsigned long n_node = this->nnode();
1048 for (unsigned long n = 0; n < n_node; n++)
1049 {
1050 // Get the pointer to the node
1051 Node* nod_pt = this->node_pt(n);
1052
1053 // Get the number of values in the node
1054 unsigned n_value = nod_pt->nvalue();
1055
1056 // We need to find if any of the values are hanging
1057 bool is_hanging = nod_pt->is_hanging();
1058 // Loop over the values and find out whether any are hanging
1059 for (unsigned n = 0; n < n_value; n++)
1060 {
1061 is_hanging |= nod_pt->is_hanging(n);
1062 }
1063
1064 // If the node is hanging then ...
1065 if (is_hanging)
1066 {
1067 // Unless they are turned into hanging nodes again below
1068 // (this might or might not happen), fill in all the necessary
1069 // data to make them 'proper' nodes again.
1070
1071 // Reconstruct the nodal values/position from the node's
1072 // hanging node representation
1073 unsigned nt = nod_pt->ntstorage();
1074 Vector<double> values(n_value);
1075 unsigned n_dim = nod_pt->ndim();
1076 Vector<double> position(n_dim);
1077 // Loop over all history values
1078 for (unsigned t = 0; t < nt; t++)
1079 {
1080 nod_pt->value(t, values);
1081 for (unsigned i = 0; i < n_value; i++)
1082 {
1083 nod_pt->set_value(t, i, values[i]);
1084 }
1085 nod_pt->position(t, position);
1086 for (unsigned i = 0; i < n_dim; i++)
1087 {
1088 nod_pt->x(t, i) = position[i];
1089 }
1090 }
1091
1092 // If it's an algebraic node: Update its previous nodal positions too
1093 AlgebraicNode* alg_node_pt = dynamic_cast<AlgebraicNode*>(nod_pt);
1094 if (alg_node_pt != 0)
1095 {
1096 bool update_all_time_levels = true;
1097 alg_node_pt->node_update(update_all_time_levels);
1098 }
1099
1100
1101 // If it's a Solid node, update Lagrangian coordinates
1102 // from its hanging node representation
1103 SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
1104 if (solid_node_pt != 0)
1105 {
1106 unsigned n_lagrangian = solid_node_pt->nlagrangian();
1107 for (unsigned i = 0; i < n_lagrangian; i++)
1108 {
1109 solid_node_pt->xi(i) = solid_node_pt->lagrangian_position(i);
1110 }
1111 }
1112
1113 // Now store geometrically hanging nodes on boundaries that
1114 // may need updating after refinement.
1115 // There will only be a problem if we have 3 spatial dimensions
1116 if ((mesh_dim > 2) && (nod_pt->is_hanging()))
1117 {
1118 // If the node is on a boundary then add a pointer to the node
1119 // to our lookup scheme
1120 if (nod_pt->is_on_boundary())
1121 {
1122 // Storage for the boundaries on which the Node is located
1123 std::set<unsigned>* boundaries_pt;
1124 nod_pt->get_boundaries_pt(boundaries_pt);
1125 if (boundaries_pt != 0)
1126 {
1127 // Loop over the boundaries and add a pointer to the node
1128 // to the appropriate storage scheme
1129 for (std::set<unsigned>::iterator it = boundaries_pt->begin();
1130 it != boundaries_pt->end();
1131 ++it)
1132 {
1133 hanging_nodes_on_boundary_pt[*it].insert(nod_pt);
1134 }
1135 }
1136 }
1137 }
1138
1139 } // End of is_hanging
1140
1141 // Initially mark all nodes as 'non-hanging' and `obsolete'
1142 nod_pt->set_nonhanging();
1143 nod_pt->set_obsolete();
1144 }
1145
1147 {
1148 t_end = TimingHelpers::timer();
1149 oomph_info << "Time for sorting out initial hanging status: "
1150 << t_end - t_start << std::endl;
1151 t_start = TimingHelpers::timer();
1152 }
1153
1154 // Unrefine all the selected elements: This needs to be
1155 //-----------------------------------------------------
1156 // all elements, because the father elements are not actually leaves.
1157 //-------------------------------------------------------------------
1158
1159 // Unrefine
1160 for (unsigned long e = 0; e < Forest_pt->ntree(); e++)
1161 {
1163 mesh_pt);
1164 }
1165
1167 {
1168 t_end = TimingHelpers::timer();
1169 oomph_info << "Time for unrefinement: " << t_end - t_start << std::endl;
1170 t_start = TimingHelpers::timer();
1171 }
1172
1173 // Add the newly created elements to mesh
1174 //---------------------------------------
1175
1176 // Stick all elements into a new vector
1177 //(note the leaves may have changed, so this is not duplicated work)
1178 Vector<Tree*> tree_nodes_pt;
1179 Forest_pt->stick_leaves_into_vector(tree_nodes_pt);
1180
1181 // Copy the elements into the mesh Vector
1182 num_tree_nodes = tree_nodes_pt.size();
1183 Element_pt.resize(num_tree_nodes);
1184 for (unsigned long e = 0; e < num_tree_nodes; e++)
1185 {
1186 Element_pt[e] = tree_nodes_pt[e]->object_pt();
1187
1188 // Now loop over all nodes in element and mark them as non-obsolete
1189 // Logic: Initially all nodes in the unrefined mesh were labeled
1190 // as deleteable. Then we create new elements (whose newly created
1191 // nodes are obviously non-obsolete), and killed some other elements (by
1192 // by deleting them and marking the nodes that were not shared by
1193 // their father as obsolete. Now we loop over all the remaining
1194 // elements and (re-)label all their nodes as non-obsolete. This
1195 // saves some nodes that were regarded as obsolete by deleted
1196 // elements but are still required in some surviving ones
1197 // from a tragic early death...
1198 FiniteElement* this_el_pt = this->finite_element_pt(e);
1199 unsigned n_node = this_el_pt->nnode(); // caching pre-loop
1200 for (unsigned n = 0; n < n_node; n++)
1201 {
1202 this_el_pt->node_pt(n)->set_non_obsolete();
1203 }
1204 }
1205
1206
1208 {
1209 t_end = TimingHelpers::timer();
1210 oomph_info << "Time for adding elements to mesh: " << t_end - t_start
1211 << std::endl;
1212 t_start = TimingHelpers::timer();
1213 }
1214
1215 // Cannot delete nodes that are still marked as obsolete
1216 // because they may still be required to assemble the hanging schemes
1217 //-------------------------------------------------------------------
1218
1219 // Mark up hanging nodes
1220 //----------------------
1221
1222 // Output streams for the hanging nodes
1223 Vector<std::ofstream*> hanging_output_files;
1224 // Setup the output files for hanging nodes, this must be called
1225 // precisely once for the forest. Note that the files will only
1226 // actually be opened if doc_info.is_doc_enabled() is true
1227 Forest_pt->open_hanging_node_files(doc_info, hanging_output_files);
1228
1229 for (unsigned long e = 0; e < num_tree_nodes; e++)
1230 {
1231 // Generic setup
1232 tree_nodes_pt[e]->object_pt()->setup_hanging_nodes(
1233 hanging_output_files);
1234 // Element specific setup
1235 tree_nodes_pt[e]->object_pt()->further_setup_hanging_nodes();
1236 }
1237
1238 // Close the hanging node files and delete the memory allocated
1239 // for the streams
1240 Forest_pt->close_hanging_node_files(doc_info, hanging_output_files);
1241
1242
1244 {
1245 t_end = TimingHelpers::timer();
1246 oomph_info << "Time for setup_hanging_nodes() and "
1247 "further_setup_hanging_nodes() for "
1248 << num_tree_nodes << " elements: " << t_end - t_start
1249 << std::endl;
1250 t_start = TimingHelpers::timer();
1251 }
1252
1253 // Read out the number of continously interpolated values
1254 // from one of the elements (assuming it's the same in all elements)
1255 unsigned ncont_interpolated_values =
1256 tree_nodes_pt[0]->object_pt()->ncont_interpolated_values();
1257
1258 // Complete the hanging nodes schemes by dealing with the
1259 // recursively hanging nodes
1260 complete_hanging_nodes(ncont_interpolated_values);
1261
1263 {
1264 t_end = TimingHelpers::timer();
1265 oomph_info << "Time for complete_hanging_nodes: " << t_end - t_start
1266 << std::endl;
1267 t_start = TimingHelpers::timer();
1268 }
1269
1270 /// Update the boundary element info -- this can be a costly procedure
1271 /// and for this reason the mesh writer might have decided not to
1272 /// set up this scheme. If so, we won't change this and suppress
1273 /// its creation...
1275 {
1277 }
1278
1280 {
1281 t_end = TimingHelpers::timer();
1282 oomph_info << "Time for boundary element info: " << t_end - t_start
1283 << std::endl;
1284 t_start = TimingHelpers::timer();
1285 }
1286
1287#ifdef PARANOID
1288
1289 // Doc/check the neighbours
1290 //-------------------------
1291 Vector<Tree*> all_tree_nodes_pt;
1292 Forest_pt->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
1293
1294 // Check the neighbours
1296
1297 // Check the integrity of the elements
1298 // -----------------------------------
1299
1300 // Loop over elements and get the elemental integrity
1301 double max_error = 0.0;
1302 for (unsigned long e = 0; e < num_tree_nodes; e++)
1303 {
1304 double max_el_error;
1305 tree_nodes_pt[e]->object_pt()->check_integrity(max_el_error);
1306 // If the elemental error is greater than our maximum error
1307 // reset the maximum
1308 if (max_el_error > max_error)
1309 {
1310 max_error = max_el_error;
1311 }
1312 }
1313
1315 {
1316 std::ostringstream error_stream;
1317 error_stream << "Mesh refined: Max. error in integrity check: "
1318 << max_error << " is too big\n";
1319 error_stream
1320 << "i.e. bigger than RefineableElement::max_integrity_tolerance()="
1322 << std::endl;
1323
1324 std::ofstream some_file;
1325 some_file.open("ProblemMesh.dat");
1326 for (unsigned long n = 0; n < n_node; n++)
1327 {
1328 // Get the pointer to the node
1329 Node* nod_pt = this->node_pt(n);
1330 // Get the dimension
1331 unsigned n_dim = nod_pt->ndim();
1332 // Output the coordinates
1333 for (unsigned i = 0; i < n_dim; i++)
1334 {
1335 some_file << this->node_pt(n)->x(i) << " ";
1336 }
1337 some_file << std::endl;
1338 }
1339 some_file.close();
1340
1341 error_stream << "Doced problem mesh in ProblemMesh.dat" << std::endl;
1342
1343 throw OomphLibError(
1344 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1345 }
1346 else
1347 {
1348 oomph_info << "Mesh refined: Max. error in integrity check: "
1349 << max_error << " is OK" << std::endl;
1351 << "i.e. less than RefineableElement::max_integrity_tolerance()="
1353 << std::endl;
1354 }
1355
1356
1358 {
1359 t_end = TimingHelpers::timer();
1360 oomph_info << "Time for (paranoid only) checking of integrity: "
1361 << t_end - t_start << std::endl;
1362 t_start = TimingHelpers::timer();
1363 }
1364
1365#endif
1366
1367 // Loop over all elements other than the final level and deactivate the
1368 // objects, essentially set the pointer that point to nodes that are
1369 // about to be deleted to NULL. This must take place here because nodes
1370 // addressed by elements that are dead but still living in the tree might
1371 // have been made obsolete in the last round of refinement
1372 for (unsigned long e = 0; e < Forest_pt->ntree(); e++)
1373 {
1376 }
1377
1378 // Now we can prune the dead nodes from the mesh.
1379 Vector<Node*> deleted_node_pt = this->prune_dead_nodes();
1380
1382 {
1383 t_end = TimingHelpers::timer();
1384 oomph_info << "Time for deactivating objects and pruning nodes: "
1385 << t_end - t_start << std::endl;
1386 t_start = TimingHelpers::timer();
1387 }
1388
1389 // Finally: Reorder the nodes within the mesh's node vector
1390 // to establish a standard ordering regardless of the sequence
1391 // of mesh refinements -- this is required to allow dump/restart
1392 // on refined meshes
1393 this->reorder_nodes();
1394
1396 {
1397 t_end = TimingHelpers::timer();
1398 oomph_info << "Time for reordering " << nnode()
1399 << " nodes: " << t_end - t_start << std::endl;
1400 t_start = TimingHelpers::timer();
1401 }
1402
1403 // Now we can correct the nodes on boundaries that were hanging that
1404 // are no longer hanging
1405 // Only bother if we have more than two dimensions
1406 if (mesh_dim > 2)
1407 {
1408 // Loop over the boundaries
1409 for (unsigned b = 0; b < n_boundary; b++)
1410 {
1411 // Remove deleted nodes from the set
1412 unsigned n_del = deleted_node_pt.size();
1413 for (unsigned j = 0; j < n_del; j++)
1414 {
1415 hanging_nodes_on_boundary_pt[b].erase(deleted_node_pt[j]);
1416 }
1417
1418 // If the nodes that were hanging are still hanging then remove them
1419 // from the set (note increment is not in for command for efficiencty)
1420 for (std::set<Node*>::iterator it =
1421 hanging_nodes_on_boundary_pt[b].begin();
1422 it != hanging_nodes_on_boundary_pt[b].end();)
1423 {
1424 if ((*it)->is_hanging())
1425 {
1426 hanging_nodes_on_boundary_pt[b].erase(it++);
1427 }
1428 else
1429 {
1430 ++it;
1431 }
1432 }
1433
1434 // Are there any nodes that have changed geometric hanging status
1435 // on the boundary
1436 // The slightly painful part is that we must adjust the position
1437 // via the macro-elements which are only available through the
1438 // elements and not the nodes.
1439 if (hanging_nodes_on_boundary_pt[b].size() > 0)
1440 {
1441 // If so we loop over all elements adjacent to the boundary
1442 unsigned n_boundary_element = this->nboundary_element(b);
1443 for (unsigned e = 0; e < n_boundary_element; ++e)
1444 {
1445 // Get a pointer to the element
1446 FiniteElement* el_pt = this->boundary_element_pt(b, e);
1447
1448 // Do we have a solid element
1449 SolidFiniteElement* solid_el_pt =
1450 dynamic_cast<SolidFiniteElement*>(el_pt);
1451
1452 // Determine whether there is a macro element
1453 bool macro_present = (el_pt->macro_elem_pt() != 0);
1454 // Or a solid macro element
1455 if (solid_el_pt != 0)
1456 {
1457 macro_present |= (solid_el_pt->undeformed_macro_elem_pt() != 0);
1458 }
1459
1460 // Only bother to do anything if there is a macro element
1461 // or undeformed macro element in a SolidElement
1462 if (macro_present)
1463 {
1464 // Loop over the nodes
1465 // ALH: (could optimise to only loop over
1466 // node associated with the boundary with more effort)
1467 unsigned n_el_node = el_pt->nnode();
1468 for (unsigned n = 0; n < n_el_node; n++)
1469 {
1470 // Cache pointer to the node
1471 Node* nod_pt = el_pt->node_pt(n);
1472 if (nod_pt->is_on_boundary(b))
1473 {
1474 // Is the Node in our set
1475 std::set<Node*>::iterator it =
1476 hanging_nodes_on_boundary_pt[b].find(nod_pt);
1477
1478 // If we have found the Node then update the position
1479 // to be consistent with the macro-element representation
1480 if (it != hanging_nodes_on_boundary_pt[b].end())
1481 {
1482 // Specialise local and global coordinates to 3D
1483 // because there is only a problem in 3D.
1484 Vector<double> s(3), x(3);
1485
1486 // Find the local coordinate of the ndoe
1487 el_pt->local_coordinate_of_node(n, s);
1488
1489 // Find the number of time history values
1490 const unsigned ntstorage = nod_pt->ntstorage();
1491
1492 // Do we have a solid node
1493 SolidNode* solid_node_pt =
1494 dynamic_cast<SolidNode*>(nod_pt);
1495 if (solid_node_pt)
1496 {
1497 // Assign Lagrangian coordinates from undeformed
1498 // macro element (if it has one -- get_x_and_xi()
1499 // does "the right thing" anyway. Leave actual
1500 // nodal positions alone -- we're doing a solid
1501 // mechanics problem and once we're going
1502 // the nodal positions are always computed, never
1503 // (automatically) reset to macro-element based
1504 // positions; not even on pinned boundaries
1505 // because the user may have other ideas about where
1506 // these should go -- pinning means "don't touch the
1507 // value", not "leave where the macro-element thinks
1508 // it should be"
1509 Vector<double> x_fe(3), xi(3), xi_fe(3);
1510 solid_el_pt->get_x_and_xi(s, x_fe, x, xi_fe, xi);
1511 for (unsigned i = 0; i < 3; i++)
1512 {
1513 solid_node_pt->xi(i) = xi[i];
1514 }
1515 }
1516 else
1517 {
1518 // Set position and history values from the
1519 // macro-element representation
1520 for (unsigned t = 0; t < ntstorage; t++)
1521 {
1522 // Get the history value from the macro element
1523 el_pt->get_x(t, s, x);
1524
1525 // Set the coordinate to that of the macroelement
1526 // representation
1527 for (unsigned i = 0; i < 3; i++)
1528 {
1529 nod_pt->x(t, i) = x[i];
1530 }
1531 }
1532 } // End of non-solid node case
1533
1534 // Now remove the node from the list
1535 hanging_nodes_on_boundary_pt[b].erase(it);
1536
1537 // If there are no Nodes left then exit the loops
1538 if (hanging_nodes_on_boundary_pt[b].size() == 0)
1539 {
1540 e = n_boundary_element;
1541 break;
1542 }
1543 }
1544 }
1545 }
1546 } // End of macro element case
1547 }
1548 }
1549 }
1550 } // End of case when we have fixed nodal positions
1551
1552
1553 // Final doc
1554 //-----------
1556 {
1557 // Doc the boundary conditions ('0' for non-existent, '1' for free,
1558 //----------------------------------------------------------------
1559 // '2' for pinned -- ideal for tecplot scatter sizing.
1560 //----------------------------------------------------
1561 // num_tree_nodes=tree_nodes_pt.size();
1562
1563 // Determine maximum number of values at any node in this type of
1564 // element
1565 RefineableElement* el_pt = tree_nodes_pt[0]->object_pt();
1566 // Initalise max_nval
1567 unsigned max_nval = 0;
1568 for (unsigned n = 0; n < el_pt->nnode(); n++)
1569 {
1570 if (el_pt->node_pt(n)->nvalue() > max_nval)
1571 {
1572 max_nval = el_pt->node_pt(n)->nvalue();
1573 }
1574 }
1575
1576 // Open the output file
1577 std::ofstream bcs_file;
1578 fullname.str("");
1579 fullname << doc_info.directory() << "/bcs" << doc_info.number()
1580 << ".dat";
1581 bcs_file.open(fullname.str().c_str());
1582
1583 // Loop over elements
1584 for (unsigned long e = 0; e < num_tree_nodes; e++)
1585 {
1586 el_pt = tree_nodes_pt[e]->object_pt();
1587 // Loop over nodes in element
1588 unsigned n_nod = el_pt->nnode();
1589 for (unsigned n = 0; n < n_nod; n++)
1590 {
1591 // Get pointer to the node
1592 Node* nod_pt = el_pt->node_pt(n);
1593 // Find the dimension of the node
1594 unsigned n_dim = nod_pt->ndim();
1595 // Write the nodal coordinates to the file
1596 for (unsigned i = 0; i < n_dim; i++)
1597 {
1598 bcs_file << nod_pt->x(i) << " ";
1599 }
1600
1601 // Loop over all values in this element
1602 for (unsigned i = 0; i < max_nval; i++)
1603 {
1604 // Value exists at this node:
1605 if (i < nod_pt->nvalue())
1606 {
1607 bcs_file << " " << 1 + nod_pt->is_pinned(i);
1608 }
1609 // ...if not just dump out a zero
1610 else
1611 {
1612 bcs_file << " 0 ";
1613 }
1614 }
1615 bcs_file << std::endl;
1616 }
1617 }
1618 bcs_file.close();
1619
1620 // Doc all nodes
1621 //---------------
1622 std::ofstream all_nodes_file;
1623 fullname.str("");
1624 fullname << doc_info.directory() << "/all_nodes" << doc_info.number()
1625 << ".dat";
1626 all_nodes_file.open(fullname.str().c_str());
1627
1628 all_nodes_file << "ZONE \n";
1629
1630 // Need to recompute the number of nodes since it may have
1631 // changed during mesh refinement/unrefinement
1632 n_node = this->nnode();
1633 for (unsigned long n = 0; n < n_node; n++)
1634 {
1635 Node* nod_pt = this->node_pt(n);
1636 unsigned n_dim = nod_pt->ndim();
1637 for (unsigned i = 0; i < n_dim; i++)
1638 {
1639 all_nodes_file << this->node_pt(n)->x(i) << " ";
1640 }
1641 all_nodes_file << std::endl;
1642 }
1643
1644 all_nodes_file.close();
1645
1646
1647 // Doc all hanging nodes:
1648 //-----------------------
1649 std::ofstream some_file;
1650 fullname.str("");
1651 fullname << doc_info.directory() << "/all_hangnodes"
1652 << doc_info.number() << ".dat";
1653 some_file.open(fullname.str().c_str());
1654 for (unsigned long n = 0; n < n_node; n++)
1655 {
1656 Node* nod_pt = this->node_pt(n);
1657
1658 if (nod_pt->is_hanging())
1659 {
1660 unsigned n_dim = nod_pt->ndim();
1661 for (unsigned i = 0; i < n_dim; i++)
1662 {
1663 some_file << nod_pt->x(i) << " ";
1664 }
1665
1666 // ALH: Added this to stop Solid problems seg-faulting
1667 if (this->node_pt(n)->nvalue() > 0)
1668 {
1669 some_file << " " << nod_pt->raw_value(0);
1670 }
1671 some_file << std::endl;
1672 }
1673 }
1674 some_file.close();
1675
1676 // Doc all hanging nodes and their masters
1677 // View with QHangingNodesWithMasters.mcr
1678 fullname.str("");
1679 fullname << doc_info.directory() << "/geometric_hangnodes_withmasters"
1680 << doc_info.number() << ".dat";
1681 some_file.open(fullname.str().c_str());
1682 for (unsigned long n = 0; n < n_node; n++)
1683 {
1684 Node* nod_pt = this->node_pt(n);
1685 if (nod_pt->is_hanging())
1686 {
1687 unsigned n_dim = nod_pt->ndim();
1688 unsigned nmaster = nod_pt->hanging_pt()->nmaster();
1689 some_file << "ZONE I=" << nmaster + 1 << std::endl;
1690 for (unsigned i = 0; i < n_dim; i++)
1691 {
1692 some_file << nod_pt->x(i) << " ";
1693 }
1694 some_file << " 2 " << std::endl;
1695
1696 for (unsigned imaster = 0; imaster < nmaster; imaster++)
1697 {
1698 Node* master_nod_pt =
1699 nod_pt->hanging_pt()->master_node_pt(imaster);
1700 unsigned n_dim = master_nod_pt->ndim();
1701 for (unsigned i = 0; i < n_dim; i++)
1702 {
1703 some_file << master_nod_pt->x(i) << " ";
1704 }
1705 some_file << " 1 " << std::endl;
1706 }
1707 }
1708 }
1709 some_file.close();
1710
1711 // Doc all hanging nodes and their masters
1712 // View with QHangingNodesWithMasters.mcr
1713 for (unsigned i = 0; i < ncont_interpolated_values; i++)
1714 {
1715 fullname.str("");
1716 fullname << doc_info.directory()
1717 << "/nonstandard_hangnodes_withmasters" << i << "_"
1718 << doc_info.number() << ".dat";
1719 some_file.open(fullname.str().c_str());
1720 unsigned n_nod = this->nnode();
1721 for (unsigned long n = 0; n < n_nod; n++)
1722 {
1723 Node* nod_pt = this->node_pt(n);
1724 if (nod_pt->is_hanging(i))
1725 {
1726 if (nod_pt->hanging_pt(i) != nod_pt->hanging_pt())
1727 {
1728 unsigned nmaster = nod_pt->hanging_pt(i)->nmaster();
1729 some_file << "ZONE I=" << nmaster + 1 << std::endl;
1730 unsigned n_dim = nod_pt->ndim();
1731 for (unsigned j = 0; j < n_dim; j++)
1732 {
1733 some_file << nod_pt->x(j) << " ";
1734 }
1735 some_file << " 2 " << std::endl;
1736 for (unsigned imaster = 0; imaster < nmaster; imaster++)
1737 {
1738 Node* master_nod_pt =
1739 nod_pt->hanging_pt(i)->master_node_pt(imaster);
1740 unsigned n_dim = master_nod_pt->ndim();
1741 for (unsigned j = 0; j < n_dim; j++)
1742 {
1743 // some_file << master_nod_pt->x(i) << " ";
1744 }
1745 some_file << " 1 " << std::endl;
1746 }
1747 }
1748 }
1749 }
1750 some_file.close();
1751 }
1752 } // End of documentation
1753 } // End if (this->nelement()>0)
1754
1755
1756#ifdef OOMPH_HAS_MPI
1757
1758 // Now (re-)classify halo and haloed nodes and synchronise hanging
1759 // nodes
1760 if (this->is_mesh_distributed())
1761 {
1763 }
1764
1765#endif
1766 }
1767
1768
1769 //========================================================================
1770 /// Refine mesh uniformly
1771 //========================================================================
1773 {
1774 // Select all elements for refinement
1775 unsigned long Nelement = this->nelement();
1776 for (unsigned long e = 0; e < Nelement; e++)
1777 {
1778 dynamic_cast<RefineableElement*>(this->element_pt(e))
1779 ->select_for_refinement();
1780 }
1781
1782 // Do the actual mesh adaptation
1784 }
1785
1786
1787 //========================================================================
1788 /// p-refine mesh uniformly
1789 //========================================================================
1791 {
1792 // Select all elements for refinement
1793 unsigned long Nelement = this->nelement();
1794 for (unsigned long e = 0; e < Nelement; e++)
1795 {
1796 // Get pointer to p-refineable element
1797 PRefineableElement* el_pt =
1798 dynamic_cast<PRefineableElement*>(this->element_pt(e));
1799 // Mark for p-refinement if possible. If not then p_adapt_mesh() will
1800 // report the error.
1801 if (el_pt != 0)
1802 {
1803 el_pt->select_for_p_refinement();
1804 }
1805 }
1806
1807 // Do the actual mesh adaptation
1809 }
1810
1811
1812 //========================================================================
1813 /// Refine mesh by splitting the elements identified
1814 /// by their numbers.
1815 //========================================================================
1817 const Vector<unsigned>& elements_to_be_refined)
1818 {
1819#ifdef OOMPH_HAS_MPI
1820 if (this->is_mesh_distributed())
1821 {
1822 std::ostringstream warn_stream;
1823 warn_stream << "You are attempting to refine selected elements of a "
1824 << std::endl
1825 << "distributed mesh. This may have undesired effects."
1826 << std::endl;
1827
1828 OomphLibWarning(warn_stream.str(),
1829 "TreeBasedRefineableMeshBase::refine_selected_elements()",
1830 OOMPH_EXCEPTION_LOCATION);
1831 }
1832#endif
1833
1834 // Select elements for refinement
1835 unsigned long nref = elements_to_be_refined.size();
1836 for (unsigned long e = 0; e < nref; e++)
1837 {
1838 dynamic_cast<RefineableElement*>(
1839 this->element_pt(elements_to_be_refined[e]))
1840 ->select_for_refinement();
1841 }
1842
1843 // Do the actual mesh adaptation
1844 adapt_mesh();
1845 }
1846
1847
1848 //========================================================================
1849 /// Refine mesh by splitting the elements identified
1850 /// by their pointers
1851 //========================================================================
1853 const Vector<RefineableElement*>& elements_to_be_refined_pt)
1854 {
1855#ifdef OOMPH_HAS_MPI
1856 if (this->is_mesh_distributed())
1857 {
1858 std::ostringstream warn_stream;
1859 warn_stream << "You are attempting to refine selected elements of a "
1860 << std::endl
1861 << "distributed mesh. This may have undesired effects."
1862 << std::endl;
1863
1864 OomphLibWarning(warn_stream.str(),
1865 "TreeBasedRefineableMeshBase::refine_selected_elements()",
1866 OOMPH_EXCEPTION_LOCATION);
1867 }
1868#endif
1869
1870 // Select elements for refinement
1871 unsigned long nref = elements_to_be_refined_pt.size();
1872 for (unsigned long e = 0; e < nref; e++)
1873 {
1874 elements_to_be_refined_pt[e]->select_for_refinement();
1875 }
1876
1877 // Do the actual mesh adaptation
1878 adapt_mesh();
1879 }
1880
1881
1882 //========================================================================
1883 /// Refine to same degree as the reference mesh.
1884 //========================================================================
1886 TreeBasedRefineableMeshBase* const& ref_mesh_pt)
1887 {
1888 // Assign storage for refinement pattern
1889 Vector<Vector<unsigned>> to_be_refined;
1890
1891 // Get refinement pattern of reference mesh:
1892 ref_mesh_pt->get_refinement_pattern(to_be_refined);
1893
1894 // Refine mesh according to given refinement pattern
1895 refine_base_mesh(to_be_refined);
1896 }
1897
1898
1899 //========================================================================
1900 /// Refine to same degree as the reference mesh minus one. Useful
1901 /// function for multigrid solvers; allows the easy copy of a mesh
1902 /// to the level of refinement just below the current one. Returns
1903 /// a boolean variable which indicates if the reference mesh has not
1904 /// been refined at all
1905 //========================================================================
1908 TreeBasedRefineableMeshBase* const& ref_mesh_pt)
1909 {
1910 // Assign storage for refinement pattern
1911 Vector<Vector<unsigned>> to_be_refined;
1912
1913 // Get refinement pattern of reference mesh:
1914 ref_mesh_pt->get_refinement_pattern(to_be_refined);
1915
1916 // Find the length of the vector
1917 unsigned nrefinement_levels = to_be_refined.size();
1918
1919 // If the reference mesh has not been refined a single time then
1920 // we cannot create an unrefined copy so stop here
1921 if (nrefinement_levels == 0)
1922 {
1923 return false;
1924 }
1925 // If the reference mesh has been refined at least once
1926 else
1927 {
1928 // Remove the last entry of the vector to make sure we refine to
1929 // the same level minus one
1930 to_be_refined.resize(nrefinement_levels - 1);
1931
1932 // Refine mesh according to given refinement pattern
1933 refine_base_mesh(to_be_refined);
1934
1935 // Indicate that it was possible to create an unrefined copy
1936 return true;
1937 }
1938 }
1939
1940
1941 //========================================================================
1942 /// Refine mesh once so that its topology etc becomes that of the
1943 /// (finer!) reference mesh -- if possible! Useful for meshes in multigrid
1944 /// hierarchies. If the meshes are too different and the conversion
1945 /// cannot be performed, the code dies (provided PARANOID is enabled).
1946 //========================================================================
1948 TreeBasedRefineableMeshBase* const& ref_mesh_pt)
1949 {
1950 oomph_info << "WARNING : This has not been checked comprehensively yet"
1951 << std::endl
1952 << "Check it and remove this break " << std::endl;
1953 pause("Yes really pause");
1954
1955#ifdef PARANOID
1956 // The max. refinement levels of the two meshes need to differ
1957 // by one, otherwise what we're doing here doesn't make sense.
1958 unsigned my_min, my_max;
1959 get_refinement_levels(my_min, my_max);
1960
1961 unsigned ref_min, ref_max;
1962 ref_mesh_pt->get_refinement_levels(ref_min, ref_max);
1963
1964 if (ref_max != my_max + 1)
1965 {
1966 std::ostringstream error_stream;
1967 error_stream
1968 << "Meshes definitely don't differ by one refinement level \n"
1969 << "max. refinement levels: " << ref_max << " " << my_max << std::endl;
1970
1971 throw OomphLibError(
1972 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1973 }
1974#endif
1975
1976 // Vector storing the elements of the uniformly unrefined mesh
1977 Vector<Tree*> coarse_elements_pt;
1978
1979 // Map storing which father elements have already been added to coarse mesh
1980 // (Default return is 0).
1981 std::map<Tree*, unsigned> father_element_included;
1982
1983 // Extract active elements (=leaf nodes in the quadtree) from reference
1984 // mesh.
1985 Vector<Tree*> leaf_nodes_pt;
1986 ref_mesh_pt->forest_pt()->stick_leaves_into_vector(leaf_nodes_pt);
1987
1988 // Loop over all elements (in their quadtree impersonation) and
1989 // check if their fathers's sons are all leaves too:
1990 unsigned nelem = leaf_nodes_pt.size();
1991 for (unsigned e = 0; e < nelem; e++)
1992 {
1993 // Pointer to leaf node
1994 Tree* leaf_pt = leaf_nodes_pt[e];
1995
1996 // Get pointer to father:
1997 Tree* father_pt = leaf_pt->father_pt();
1998
1999 // If we don't have a father we're at the root level in which
2000 // case this element can't be unrefined.
2001 if (0 == father_pt)
2002 {
2003 coarse_elements_pt.push_back(leaf_pt);
2004 }
2005 else
2006 {
2007 // Loop over the father's sons to check if they're
2008 // all non-leafs, i.e. if they can be unrefined
2009 bool can_unrefine = true;
2010 unsigned n_sons = father_pt->nsons();
2011 for (unsigned i = 0; i < n_sons; i++)
2012 {
2013 // If (at least) one of the sons is not a leaf, we can't unrefine
2014 if (!father_pt->son_pt(i)->is_leaf()) can_unrefine = false;
2015 }
2016
2017 // If we can unrefine, the father element will be
2018 // an element in the coarse mesh, the sons won't
2019 if (can_unrefine)
2020 {
2021 if (father_element_included[father_pt] == 0)
2022 {
2023 coarse_elements_pt.push_back(father_pt);
2024 father_element_included[father_pt] = 1;
2025 }
2026 }
2027 // Son will still be there on the coarse mesh
2028 else
2029 {
2030 coarse_elements_pt.push_back(leaf_pt);
2031 }
2032 }
2033 }
2034
2035 // Number of elements in ref mesh if it was unrefined uniformly:
2036 unsigned nel_coarse = coarse_elements_pt.size();
2037
2038
2039#ifdef PARANOID
2040 bool stop_it = false;
2041 // The numbers had better match otherwise we might as well stop now...
2042 if (nel_coarse != this->nelement())
2043 {
2044 oomph_info << "Number of elements in uniformly unrefined reference mesh: "
2045 << nel_coarse << std::endl;
2046 oomph_info << "Number of elements in 'this' mesh: " << nel_coarse
2047 << std::endl;
2048 oomph_info << "don't match" << std::endl;
2049 stop_it = true;
2050 }
2051#endif
2052
2053 // Now loop over all elements in uniformly coarsened reference mesh
2054 // and check if add the number of any element that was created
2055 // by having had its sons merged to the vector of elements that
2056 // need to get refined if we go the other way
2057 Vector<unsigned> elements_to_be_refined;
2058 for (unsigned i = 0; i < nel_coarse; i++)
2059 {
2060 if (father_element_included[coarse_elements_pt[i]] == 1)
2061 {
2062 elements_to_be_refined.push_back(i);
2063 }
2064 }
2065
2066
2067#ifdef PARANOID
2068 // Doc troublesome meshes:
2069 if (stop_it)
2070 {
2071 std::ofstream some_file;
2072 some_file.open("orig_mesh.dat");
2073 this->output(some_file);
2074 some_file.close();
2075 oomph_info << "Documented original ('this')mesh in orig_mesh.dat"
2076 << std::endl;
2077 }
2078#endif
2079
2080
2081 // Now refine precisely these elements in "this" mesh.
2082 refine_selected_elements(elements_to_be_refined);
2083
2084
2085#ifdef PARANOID
2086
2087 // Check if the nodal positions of all element's nodes agree
2088 // in the two fine meshes:
2089 double tol = 1.0e-5;
2090 for (unsigned e = 0; e < nelem; e++)
2091 {
2092 // Get elements
2093 FiniteElement* ref_el_pt = ref_mesh_pt->finite_element_pt(e);
2094 FiniteElement* el_pt = this->finite_element_pt(e);
2095
2096 // Loop over nodes
2097 unsigned nnod = ref_el_pt->nnode();
2098 for (unsigned j = 0; j < nnod; j++)
2099 {
2100 // Get nodes
2101 Node* ref_node_pt = ref_el_pt->node_pt(j);
2102 Node* node_pt = el_pt->node_pt(j);
2103
2104 // Check error in position
2105 double error = 0.0;
2106 unsigned ndim = node_pt->ndim();
2107 for (unsigned i = 0; i < ndim; i++)
2108 {
2109 error += pow(node_pt->x(i) - ref_node_pt->x(i), 2);
2110 }
2111 error = sqrt(error);
2112
2113 if (error > tol)
2114 {
2115 oomph_info << "Error in nodal position of node " << j << ": " << error
2116 << " [tol=" << tol << "]" << std::endl;
2117 stop_it = true;
2118 }
2119 }
2120 }
2121
2122 // Do we have a death wish?
2123 if (stop_it)
2124 {
2125 // Doc troublesome meshes:
2126 std::ofstream some_file;
2127 some_file.open("refined_mesh.dat");
2128 this->output(some_file);
2129 some_file.close();
2130
2131 some_file.open("finer_mesh.dat");
2132 ref_mesh_pt->output(some_file);
2133 some_file.close();
2134
2135 throw OomphLibError(
2136 "Bailing out. Doced refined_mesh.dat finer_mesh.dat\n",
2137 OOMPH_CURRENT_FUNCTION,
2138 OOMPH_EXCEPTION_LOCATION);
2139 }
2140
2141#endif
2142 }
2143
2144
2145 //========================================================================
2146 /// Unrefine mesh uniformly. Return 0 for success,
2147 /// 1 for failure (if unrefinement has reached the coarsest permitted
2148 /// level)
2149 //========================================================================
2151 {
2152 // We can't just select all elements for unrefinement
2153 // because they need to merge with their brothers.
2154 // --> Rather than repeating the convoluted logic of
2155 // RefineableQuadMesh<ELEMENT>::adapt(Vector<double>& elemental_error)
2156 // here (code duplication!) hack it by filling the error
2157 // vector with values that ensure unrefinement for all
2158 // elements where this is possible
2159
2160 // Create dummy vector for elemental errors
2161 unsigned long Nelement = this->nelement();
2162 Vector<double> elemental_error(Nelement);
2163
2164 // In order to force unrefinement, set the min permitted error to
2165 // be the default and then set the actual error to be below this.
2166 // This avoids problems when the actual min error is zero (or small)
2167 // For sanity's sake, also set the max permitted error back to the default
2168 // so that we have a max error bigger than a min error
2169 const double current_min_error = this->min_permitted_error();
2170 const double current_max_error = this->max_permitted_error();
2171
2172 this->min_permitted_error() = 1.0e-5;
2173 this->max_permitted_error() = 1.0e-3;
2174
2175 double error = min_permitted_error() / 100.0;
2176 for (unsigned long e = 0; e < Nelement; e++)
2177 {
2178 elemental_error[e] = error;
2179 }
2180
2181 // Temporarily lift any restrictions on the minimum number of
2182 // elements that need to be unrefined to make it worthwhile
2183 unsigned backup = max_keep_unrefined();
2184 max_keep_unrefined() = 0;
2185
2186 // Do the actual mesh adaptation with fake error vector
2187 adapt(elemental_error);
2188
2189 // Reset the minimum number of elements that need to be unrefined
2190 // to make it worthwhile
2191 max_keep_unrefined() = backup;
2192
2193 // Now restore the error tolerances
2194 this->min_permitted_error() = current_min_error;
2195 this->max_permitted_error() = current_max_error;
2196
2197 // Has the unrefinement actually changed anything?
2198 if (Nelement == this->nelement())
2199 {
2200 return 1;
2201 }
2202 else
2203 {
2204 return 0;
2205 }
2206 }
2207
2208 //==================================================================
2209 /// Given a node, return a vector of pointers to master nodes and a
2210 /// vector of the associated weights.
2211 /// This is done recursively, so if a node is not hanging,
2212 /// the node is regarded as its own master node which has weight 1.0.
2213 //==================================================================
2215 Node*& nod_pt,
2216 Vector<Node*>& master_nodes,
2217 Vector<double>& hang_weights,
2218 const int& i)
2219 {
2220 // Is the node hanging in the variable i
2221 if (nod_pt->is_hanging(i))
2222 {
2223 // Loop over all master nodes
2224 HangInfo* const hang_pt = nod_pt->hanging_pt(i);
2225 unsigned nmaster = hang_pt->nmaster();
2226
2227 for (unsigned m = 0; m < nmaster; m++)
2228 {
2229 // Get the master node
2230 Node* master_nod_pt = hang_pt->master_node_pt(m);
2231
2232 // Keep in memory the size of the list before adding the nodes this
2233 // master node depends on. This is required so that the recursion is
2234 // only performed on these particular master nodes. A master node
2235 // could contain contributions from two separate pseudo-masters.
2236 // These contributions must be summed, not multiplied.
2237 int first_new_node = master_nodes.size();
2238
2239 // Now check which master nodes this master node depends on
2241 master_nod_pt, master_nodes, hang_weights, i);
2242
2243 // Multiply old weight by new weight for all the nodes this master
2244 // node depends on
2245 unsigned n_new_master_node = master_nodes.size();
2246
2247 double mtr_weight = hang_pt->master_weight(m);
2248
2249 for (unsigned k = first_new_node; k < n_new_master_node; k++)
2250 {
2251 hang_weights[k] = mtr_weight * hang_weights[k];
2252 }
2253 }
2254 }
2255 else
2256 // Node isn't hanging so it enters itself with the full weight
2257 {
2258 master_nodes.push_back(nod_pt);
2259 hang_weights.push_back(1.0);
2260 }
2261 }
2262
2263
2264 //==================================================================
2265 /// Complete the hanging node scheme recursively.
2266 /// After the initial markup scheme, hanging nodes
2267 /// can depend on other hanging nodes ---> AAAAAAAAARGH!
2268 /// Need to translate this into a scheme where all
2269 /// hanging nodes only depend on non-hanging nodes...
2270 //==================================================================
2272 const int& ncont_interpolated_values)
2273 {
2274 // Number of nodes in mesh
2275 unsigned long n_node = this->nnode();
2276 double min_weight = 1.0e-8; // RefineableBrickElement::min_weight_value();
2277
2278 // Loop over the nodes in the mesh
2279 for (unsigned long n = 0; n < n_node; n++)
2280 {
2281 // Assign a local pointer to the node
2282 Node* nod_pt = this->node_pt(n);
2283
2284 // Loop over the values,
2285 // N.B. geometric hanging data is stored at the index -1
2286 for (int i = -1; i < ncont_interpolated_values; i++)
2287 {
2288 // Is the node hanging?
2289 if (nod_pt->is_hanging(i))
2290 {
2291 // If it is geometric OR has hanging node data that differs
2292 // from the geometric data, we must do some work
2293 if ((i == -1) || (nod_pt->hanging_pt(i) != nod_pt->hanging_pt()))
2294 {
2295 // Find out the ultimate map of dependencies: Master nodes
2296 // and associated weights
2297 Vector<Node*> master_nodes;
2298 Vector<double> hanging_weights;
2300 nod_pt, master_nodes, hanging_weights, i);
2301
2302 // put them into a map to merge all the occurences of the same node
2303 // (add the weights)
2304 std::map<Node*, double> hang_weights;
2305 unsigned n_master = master_nodes.size();
2306 for (unsigned k = 0; k < n_master; k++)
2307 {
2308 if (std::fabs(hanging_weights[k]) > min_weight)
2309 hang_weights[master_nodes[k]] += hanging_weights[k];
2310 }
2311
2312 // Create new hanging data (we know how many data there are)
2313 HangInfo* hang_pt = new HangInfo(hang_weights.size());
2314
2315 unsigned hang_weights_index = 0;
2316 // Copy the map into the HangInfo object
2317 typedef std::map<Node*, double>::iterator IT;
2318 for (IT it = hang_weights.begin(); it != hang_weights.end(); ++it)
2319 {
2320 hang_pt->set_master_node_pt(
2321 hang_weights_index, it->first, it->second);
2322 ++hang_weights_index;
2323 }
2324
2325 // Assign the new hanging pointer to the appropriate value
2326 nod_pt->set_hanging_pt(hang_pt, i);
2327 }
2328 }
2329 }
2330 }
2331
2332#ifdef PARANOID
2333
2334 // Check hanging node scheme: The weights need to add up to one
2335 //-------------------------------------------------------------
2336 // Loop over all values indices
2337 for (int i = -1; i < ncont_interpolated_values; i++)
2338 {
2339 // Loop over all nodes in mesh
2340 for (unsigned long n = 0; n < n_node; n++)
2341 {
2342 // Set a local pointer to the node
2343 Node* nod_pt = this->node_pt(n);
2344
2345 // Is it hanging?
2346 if (nod_pt->is_hanging(i))
2347 {
2348 unsigned nmaster = nod_pt->hanging_pt(i)->nmaster();
2349 double sum = 0.0;
2350 for (unsigned imaster = 0; imaster < nmaster; imaster++)
2351 {
2352 sum += nod_pt->hanging_pt(i)->master_weight(imaster);
2353 }
2354 if (std::fabs(sum - 1.0) > 1.0e-7)
2355 {
2356 oomph_info << "WARNING: Sum of master node weights fabs(sum-1.0) "
2357 << std::fabs(sum - 1.0) << " for node number " << n
2358 << " at value " << i << std::endl;
2359 }
2360 }
2361 }
2362 }
2363#endif
2364 }
2365
2366 // Sorting out the cases of nodes that are hanging on at least one but not
2367 // all of the processors for which they are part of the local mesh.
2368
2369#ifdef OOMPH_HAS_MPI
2370
2371 //========================================================================
2372 /// Deal with nodes that are hanging on one process but not another
2373 /// (i.e. the hanging status of the haloed and halo layers disagrees)
2374 //========================================================================
2376 const unsigned& ncont_interpolated_values)
2377 {
2378 // Store number of processors and current process
2379 MPI_Status status;
2380 int n_proc = Comm_pt->nproc();
2381 int my_rank = Comm_pt->my_rank();
2382
2383 double t_start = 0.0;
2384 double t_end = 0.0;
2385
2386 // Storage for the hanging status of halo/haloed nodes on elements
2387 Vector<Vector<int>> haloed_hanging(n_proc);
2388 Vector<Vector<int>> halo_hanging(n_proc);
2389
2390 // Counter for the number of nodes which require additional synchronisation
2391 unsigned nnode_still_requiring_synchronisation = 0;
2392
2394 {
2395 t_start = TimingHelpers::timer();
2396 }
2397
2398 // Store number of continuosly interpolated values as int
2399 int ncont_inter_values = ncont_interpolated_values;
2400
2401 // Loop over processes: Each processor checks that is haloed nodes
2402 // with proc d have consistent hanging stats with halo counterparts.
2403 for (int d = 0; d < n_proc; d++)
2404 {
2405 // No halo with self: Setup hang info for my haloed nodes with proc d
2406 // then get ready to receive halo info from processor d.
2407 if (d != my_rank)
2408 {
2409 // Loop over haloed nodes
2410 unsigned nh = nhaloed_node(d);
2411 for (unsigned j = 0; j < nh; j++)
2412 {
2413 // Get node
2414 Node* nod_pt = haloed_node_pt(d, j);
2415
2416 // Loop over the hanging status for each interpolated variable
2417 // (and the geometry)
2418 for (int icont = -1; icont < ncont_inter_values; icont++)
2419 {
2420 // Store the hanging status of this haloed node
2421 if (nod_pt->is_hanging(icont))
2422 {
2423 unsigned n_master = nod_pt->hanging_pt(icont)->nmaster();
2424 haloed_hanging[d].push_back(n_master);
2425 }
2426 else
2427 {
2428 haloed_hanging[d].push_back(0);
2429 }
2430 }
2431 }
2432
2433 // Receive the hanging status information from the corresponding process
2434 unsigned count_haloed = haloed_hanging[d].size();
2435
2436#ifdef PARANOID
2437 // Check that number of halo and haloed data match
2438 unsigned tmp = 0;
2439 MPI_Recv(&tmp, 1, MPI_UNSIGNED, d, 0, Comm_pt->mpi_comm(), &status);
2440 if (tmp != count_haloed)
2441 {
2442 std::ostringstream error_stream;
2443 error_stream << "Number of halo data, " << tmp
2444 << ", does not match number of haloed data, "
2445 << count_haloed << std::endl;
2446 throw OomphLibError(error_stream.str(),
2447 OOMPH_CURRENT_FUNCTION,
2448 OOMPH_EXCEPTION_LOCATION);
2449 }
2450#endif
2451
2452 // Get the data (if any)
2453 if (count_haloed != 0)
2454 {
2455 halo_hanging[d].resize(count_haloed);
2456 MPI_Recv(&halo_hanging[d][0],
2457 count_haloed,
2458 MPI_INT,
2459 d,
2460 0,
2461 Comm_pt->mpi_comm(),
2462 &status);
2463 }
2464 }
2465 else // d==my_rank, i.e. current process: Send halo hanging status
2466 // to process dd where it's received (see above) and compared
2467 // and compared against the hang status of the haloed nodes
2468 {
2469 for (int dd = 0; dd < n_proc; dd++)
2470 {
2471 // No halo with yourself
2472 if (dd != d)
2473 {
2474 // Storage for halo hanging status and counter
2475 Vector<int> local_halo_hanging;
2476
2477 // Loop over halo nodes
2478 unsigned nh = nhalo_node(dd);
2479 for (unsigned j = 0; j < nh; j++)
2480 {
2481 // Get node
2482 Node* nod_pt = halo_node_pt(dd, j);
2483
2484 // Loop over the hanging status for each interpolated variable
2485 // (and the geometry)
2486 for (int icont = -1; icont < ncont_inter_values; icont++)
2487 {
2488 // Store hanging status of halo node
2489 if (nod_pt->is_hanging(icont))
2490 {
2491 unsigned n_master = nod_pt->hanging_pt(icont)->nmaster();
2492 local_halo_hanging.push_back(n_master);
2493 }
2494 else
2495 {
2496 local_halo_hanging.push_back(0);
2497 }
2498 }
2499 }
2500
2501
2502 // Send the information to the relevant process
2503 unsigned count_halo = local_halo_hanging.size();
2504
2505#ifdef PARANOID
2506 // Check that number of halo and haloed data match
2507 MPI_Send(&count_halo, 1, MPI_UNSIGNED, dd, 0, Comm_pt->mpi_comm());
2508#endif
2509
2510 // Send data (if any)
2511 if (count_halo != 0)
2512 {
2513 MPI_Send(&local_halo_hanging[0],
2514 count_halo,
2515 MPI_INT,
2516 dd,
2517 0,
2518 Comm_pt->mpi_comm());
2519 }
2520 }
2521 }
2522 }
2523 }
2524
2526 {
2527 t_end = TimingHelpers::timer();
2528 oomph_info << "Time for first all-to-all in synchronise_hanging_nodes(): "
2529 << t_end - t_start << std::endl;
2530 t_start = TimingHelpers::timer();
2531 }
2532
2533
2534 // Now compare equivalent halo and haloed vectors to find discrepancies.
2535 // It is possible that a master node may not be on either process involved
2536 // in the halo-haloed scheme; to work round this, we use the shared_node
2537 // storage scheme, which stores all nodes that are on each pair of
2538 // processors in the same order on each of the two processors
2539
2540 // Vector to store info that will help us locate master nodes on "other"
2541 // processor
2542 Vector<HangHelperStruct> hang_info;
2543
2544 // Copy vector-based representation of shared nodes into
2545 // map for faster search
2546 Vector<std::map<Node*, unsigned>> shared_node_map(n_proc);
2547 for (int d = 0; d < n_proc; d++)
2548 {
2549 unsigned n = Shared_node_pt[d].size();
2550 for (unsigned jj = 0; jj < n; jj++)
2551 {
2552 shared_node_map[d][Shared_node_pt[d][jj]] = jj;
2553 }
2554 }
2555
2556
2557 // Loop over domains: Each processor checks consistency of hang status
2558 // of its haloed nodes with proc d against the halo counterpart. Haloed
2559 // wins if there are any discrepancies.
2560 for (int d = 0; d < n_proc; d++)
2561 {
2562 // No halo with yourself
2563 if (d != my_rank)
2564 {
2565 unsigned discrepancy_count = 0;
2566 unsigned discrepancy_count_buff = 0;
2567
2568 // Storage for hanging information that needs to be sent to the
2569 // relevant process if there is a discrepancy in the hanging status
2570 Vector<int> send_data;
2571 Vector<double> send_double_data;
2572 // Buffer storage for data to be sent
2573 //(We need this because we cannot tell until the end of the loop over
2574 // a node's master nodes whether we can reconcile its hanging status)
2575 Vector<int> send_data_buff(0);
2576 Vector<double> send_double_data_buff(0);
2577
2578 // Counter for traversing haloed data
2579 unsigned count = 0;
2580
2581 // Indicate presence of discrepancy. Default: there's none
2582 unsigned discrepancy = 0;
2583 unsigned discrepancy_buff = 0;
2584
2585 // Loop over haloed nodes
2586 unsigned nh = nhaloed_node(d);
2587 for (unsigned j = 0; j < nh; j++)
2588 {
2589 // Get node
2590 Node* nod_pt = haloed_node_pt(d, j);
2591
2592 // Loop over the hanging status for each interpolated variable
2593 // (and the geometry)
2594 for (int icont = -1; icont < ncont_inter_values; icont++)
2595 {
2596 // Compare hanging status of halo/haloed counterpart structure
2597
2598 // Haloed is is hanging and haloed has different number
2599 // of master nodes (which includes none in which case it isn't
2600 // hanging)
2601 if ((haloed_hanging[d][count] > 0) &&
2602 (haloed_hanging[d][count] != halo_hanging[d][count]))
2603 {
2604 discrepancy_buff = 1;
2605 discrepancy_count_buff++;
2606
2607 // Flag to check if all masters of this node have been found
2608 bool found_all_masters = true;
2609
2610 // Find master nodes of haloed node
2611 HangInfo* hang_pt = nod_pt->hanging_pt(icont);
2612 unsigned nhd_master = hang_pt->nmaster();
2613
2614 // Add the number of master nodes to the hanging_nodes vector
2615 send_data_buff.push_back(nhd_master);
2616
2617 // Loop over master nodes required for HangInfo
2618 for (unsigned m = 0; m < nhd_master; m++)
2619 {
2620 // Get mth master node
2621 Node* master_nod_pt = hang_pt->master_node_pt(m);
2622
2623
2624 // //------------------------------------------
2625 // // Old direct search is much slower!
2626 // // Keep this code alive to allow comparisons
2627 //
2628 // double t_start_search=TimingHelpers::timer();
2629 //
2630 // // This node will be shared: find it!
2631 // bool found=false;
2632 //
2633 // // Look in shared node storage with domain d
2634 // first unsigned nnod_shared=nshared_node(d); for
2635 // (unsigned k=0; k<nnod_shared; k++)
2636 // {
2637 // if (master_nod_pt==shared_node_pt(d,k))
2638 // {
2639 // // Found a master: Put its number in the
2640 // shared
2641 // // scheme into send data
2642 // send_data.push_back(k);
2643 //
2644 // // Add the weight
2645 // send_double_data.push_back(hang_pt->master_weight(m));
2646 //
2647 // // Done
2648 // found=true;
2649 // break;
2650 // }
2651 // }
2652 //
2653 // double t_end_search=TimingHelpers::timer();
2654 // t_start_search_total+=(t_end_search-t_start_search);
2655 //
2656 // // end direct search demo
2657 // //-----------------------
2658
2659
2660 // This node will be shared: find it!
2661 bool found = false;
2662
2663 // Which processor holds the non-halo counterpart of this
2664 // node?
2665 int non_halo_proc_id = master_nod_pt->non_halo_proc_ID();
2666
2667 // Try to find node in map with proc d and get iterator to entry
2668 std::map<Node*, unsigned>::iterator it =
2669 shared_node_map[d].find(master_nod_pt);
2670
2671 // If it's not in there iterator points to end of
2672 // set
2673 if (it != shared_node_map[d].end())
2674 {
2675 // Found a master: When looking up the node in the shared
2676 // node scheme on processor d, which processor do I work
2677 // with? The current one
2678 send_data_buff.push_back(my_rank);
2679
2680 // Found a master: Send its index in the shared
2681 // node scheme
2682 send_data_buff.push_back((*it).second);
2683
2684 // Add the weight
2685 send_double_data_buff.push_back(hang_pt->master_weight(m));
2686
2687 // Done
2688 found = true;
2689 }
2690
2691 // If we haven't found it in the shared node scheme with proc d
2692 // find it in the shared node scheme with the processor that
2693 // holds the non-halo version
2694 if (!found)
2695 {
2696 // This is odd -- can't currently handle the case where
2697 // node is owned by current processor (indicated by
2698 // non_halo_proc_id being negative
2699 if (non_halo_proc_id < 0)
2700 {
2701 // This case is now handled by the function
2702 // additional_synchronise_hanging_nodes()
2703 // called (if necessary) at the end
2704 // OomphLibWarning(
2705 // "Odd: missing master node is owned by current proc. Will
2706 // crash below.",
2707 // "TreeBasedRefineableMeshBase::synchronise_hanging_nodes(...)",
2708 // OOMPH_EXCEPTION_LOCATION);
2709 }
2710 else // i.e. (non_halo_proc_id>=0)
2711 {
2712 if (shared_node_map[non_halo_proc_id].size() > 0)
2713 {
2714 std::map<Node*, unsigned>::iterator it =
2715 shared_node_map[non_halo_proc_id].find(master_nod_pt);
2716
2717 // If it's not in there iterator points to end of
2718 // set
2719 if (it != shared_node_map[non_halo_proc_id].end())
2720 {
2721 // Found a master: Send ID of processor that holds
2722 // non-halo (the fact that this is different from
2723 // my_rank (the processor that sends this) will alert
2724 // the other processor to the fact that it needs to
2725 send_data_buff.push_back(non_halo_proc_id);
2726
2727 // Found a master: Send its index in the shared
2728 // node scheme
2729 send_data_buff.push_back((*it).second);
2730
2731 // Add the weight
2732 send_double_data_buff.push_back(
2733 hang_pt->master_weight(m));
2734
2735 // Done
2736 found = true;
2737
2738 // It is possible that the master node found in the
2739 // shared storage with processor <non_halo_proc_id> does
2740 // not actually exist on processor d. If it does turn
2741 // out to exist, we are ok, but if not then we must
2742 // remember to create it later (in the
2743 // additional_synchronise_hanging_nodes() function)
2744 }
2745 }
2746 }
2747 }
2748
2749 /*
2750 // Don't throw error, we will construct missing master nodes in
2751 // additional_synchronise_hanging_nodes() below
2752
2753 // Paranoid check: if we haven't found the master node
2754 // then throw an error
2755 if (!found)
2756 {
2757 char filename[100];
2758 std::ofstream some_file;
2759 sprintf(filename,"sync_hanging_node_crash_mesh_proc%i.dat",
2760 my_rank);
2761 some_file.open(filename);
2762 this->output(some_file);
2763 some_file.close();
2764
2765 sprintf(filename,
2766 "sync_hanging_node_crash_mesh_with_haloes_proc%i.dat",
2767 my_rank);
2768 some_file.open(filename);
2769 this->enable_output_of_halo_elements();
2770 this->output(some_file);
2771 this->disable_output_of_halo_elements();
2772 some_file.close();
2773
2774
2775 std::set<unsigned> other_proc_id;
2776 other_proc_id.insert(d);
2777 other_proc_id.insert(non_halo_proc_id);
2778 for (std::set<unsigned>::iterator it=other_proc_id.begin();
2779 it!=other_proc_id.end();it++)
2780 {
2781 unsigned d_doc=(*it);
2782
2783 sprintf(
2784 filename,
2785 "sync_hanging_node_crash_halo_elements_with_proc%i_proc%i.dat",
2786 d_doc,my_rank);
2787 some_file.open(filename);
2788 Vector<GeneralisedElement*>
2789 halo_elem_pt(this->halo_element_pt(d_doc));
2790 unsigned nelem=halo_elem_pt.size();
2791 for (unsigned e=0;e<nelem;e++)
2792 {
2793 FiniteElement* el_pt=
2794 dynamic_cast<FiniteElement*>(halo_elem_pt[e]);
2795 el_pt->output(some_file);
2796 }
2797 some_file.close();
2798
2799 sprintf(
2800 filename,
2801 "sync_hanging_node_crash_haloed_elements_with_proc%i_proc%i.dat",
2802 d_doc,my_rank);
2803 some_file.open(filename);
2804 Vector<GeneralisedElement*>
2805 haloed_elem_pt(this->haloed_element_pt(d_doc));
2806 nelem=haloed_elem_pt.size();
2807 for (unsigned e=0;e<nelem;e++)
2808 {
2809 FiniteElement* el_pt=
2810 dynamic_cast<FiniteElement*>(haloed_elem_pt[e]);
2811 el_pt->output(some_file);
2812 }
2813 some_file.close();
2814
2815
2816 sprintf(
2817 filename,
2818 "sync_hanging_node_crash_shared_nodes_with_proc%i_proc%i.dat",
2819 d_doc,my_rank);
2820 some_file.open(filename);
2821 unsigned n=nshared_node(d_doc);
2822 for (unsigned j=0;j<n;j++)
2823 {
2824 Node* nod_pt=shared_node_pt(d_doc,j);
2825 unsigned nd=nod_pt->ndim();
2826 for (unsigned i=0;i<nd;i++)
2827 {
2828 some_file << nod_pt->x(i) << " ";
2829 }
2830 some_file << "\n";
2831 }
2832 some_file.close();
2833
2834
2835 sprintf(
2836 filename,
2837 "sync_hanging_node_crash_halo_nodes_with_proc%i_proc%i.dat",
2838 d_doc,my_rank);
2839 some_file.open(filename);
2840 n=nhalo_node(d_doc);
2841 for (unsigned j=0;j<n;j++)
2842 {
2843 Node* nod_pt=halo_node_pt(d_doc,j);
2844 unsigned nd=nod_pt->ndim();
2845 for (unsigned i=0;i<nd;i++)
2846 {
2847 some_file << nod_pt->x(i) << " ";
2848 }
2849 some_file << "\n";
2850 }
2851 some_file.close();
2852
2853
2854 sprintf(
2855 filename,
2856 "sync_hanging_node_crash_haloed_nodes_with_proc%i_proc%i.dat",
2857 d_doc,my_rank);
2858 some_file.open(filename);
2859 n=nhaloed_node(d_doc);
2860 for (unsigned j=0;j<n;j++)
2861 {
2862 Node* nod_pt=haloed_node_pt(d_doc,j);
2863 unsigned nd=nod_pt->ndim();
2864 for (unsigned i=0;i<nd;i++)
2865 {
2866 some_file << nod_pt->x(i) << " ";
2867 }
2868 some_file << "\n";
2869 }
2870 some_file.close();
2871
2872 } // end of loop over all inter-processor lookup schemes
2873
2874 std::ostringstream error_stream;
2875 unsigned n=master_nod_pt->ndim();
2876 error_stream << "Error: Master node at:\n\n";
2877 for (unsigned i=0;i<n;i++)
2878 {
2879 error_stream << master_nod_pt->x(i) << " ";
2880 }
2881 error_stream << "\n\nnot found for icont="
2882 << icont << "in "
2883 << "shared node storage with proc " << d <<
2884 "\n"
2885 << "or in shared node storage with proc "
2886 << non_halo_proc_id
2887 << " which is where its non-halo counterpart
2888 lives.\n"
2889 << "Relevant files:
2890 sync_hanging_node_crash*.dat\n\n"
2891 << "Hanging node itself: \n\n";
2892 n=nod_pt->ndim();
2893 for (unsigned i=0;i<n;i++)
2894 {
2895 error_stream << nod_pt->x(i) << " ";
2896 }
2897 error_stream << nod_pt->non_halo_proc_ID();
2898 error_stream << "\n\nMaster nodes:\n\n";
2899 for (unsigned m=0; m<nhd_master; m++)
2900 {
2901 Node* master_nod_pt=hang_pt->master_node_pt(m);
2902 n=master_nod_pt->ndim();
2903 for (unsigned i=0;i<n;i++)
2904 {
2905 error_stream << master_nod_pt->x(i) << " ";
2906 }
2907 error_stream << master_nod_pt->non_halo_proc_ID();
2908 error_stream << "\n";
2909 }
2910
2911 // try to find it somewhere else -- sub-optimal search but
2912 // we're about to die on our arses (yes, plural -- this is
2913 // a parallel run!) anyway...
2914 for (int dddd=0;dddd<n_proc;dddd++)
2915 {
2916 bool loc_found=false;
2917 unsigned nnnod_shared=nshared_node(dddd);
2918 for (unsigned k=0; k<nnnod_shared; k++)
2919 {
2920 if (master_nod_pt==shared_node_pt(dddd,k))
2921 {
2922 loc_found=true;
2923 error_stream
2924 << "Found that master node as " << k
2925 << "-th entry in shared node storage with proc "
2926 << dddd << "\n";
2927 }
2928 }
2929 if (!loc_found)
2930 {
2931 error_stream
2932 << "Did not find that master node in shared node storage
2933 with proc "
2934 << dddd << "\n";
2935 }
2936 }
2937 error_stream << "\n\n";
2938
2939 throw OomphLibError(
2940 error_stream.str(),
2941 OOMPH_CURRENT_FUNCTION,
2942 OOMPH_EXCEPTION_LOCATION);
2943 }
2944 */
2945
2946 // Check if the master has been found
2947 if (!found)
2948 {
2949 // If this master hasn't been found then set the flag
2950 found_all_masters = false;
2951 // No need to continue searching for masters
2952 break;
2953 }
2954
2955
2956 } // loop over master nodes
2957
2958
2959 // Check if we need to send the data
2960 if (found_all_masters)
2961 {
2962 // All masters were found, so populate send data from buffer
2963 discrepancy = discrepancy_buff;
2964 discrepancy_count += discrepancy_count_buff;
2965 for (unsigned i = 0; i < send_data_buff.size(); i++)
2966 {
2967 send_data.push_back(send_data_buff[i]);
2968 }
2969 for (unsigned i = 0; i < send_double_data_buff.size(); i++)
2970 {
2971 send_double_data.push_back(send_double_data_buff[i]);
2972 }
2973
2974 // Clear buffers and reset
2975 discrepancy_buff = 0;
2976 discrepancy_count_buff = 0;
2977 send_data_buff.clear();
2978 send_double_data_buff.clear();
2979 }
2980 else
2981 {
2982 // At least one master node was not found, so we can't
2983 // reconcile the hanging status of this node yet. We tell
2984 // the other processor to do nothing for now.
2985 send_data.push_back(0);
2986
2987 // Clear buffers and reset
2988 discrepancy_buff = 0;
2989 discrepancy_count_buff = 0;
2990 send_data_buff.clear();
2991 send_double_data_buff.clear();
2992
2993 // Set flag to trigger another round of synchronisation
2994 nnode_still_requiring_synchronisation++;
2995 }
2996 }
2997 // Haloed node isn't hanging but halo is: the latter
2998 // shouldn't so send a -1 to indicate that it's to be made
2999 // non-hanging
3000 else if ((haloed_hanging[d][count] == 0) &&
3001 (halo_hanging[d][count] > 0))
3002 {
3003 discrepancy = 1;
3004 discrepancy_count++;
3005 send_data.push_back(-1);
3006 }
3007 // Both halo and haloed node have the same number of masters
3008 // we're happy!
3009 else if (haloed_hanging[d][count] == halo_hanging[d][count])
3010 {
3011 send_data.push_back(0);
3012 }
3013 else
3014 {
3015 std::ostringstream error_stream;
3016 error_stream << "Never get here!\n "
3017 << "haloed_hanging[d][count]="
3018 << haloed_hanging[d][count]
3019 << "; halo_hanging[d][count]="
3020 << halo_hanging[d][count] << std::endl;
3021 throw OomphLibError(error_stream.str(),
3022 OOMPH_CURRENT_FUNCTION,
3023 OOMPH_EXCEPTION_LOCATION);
3024 }
3025 // Increment counter for number of haloed data
3026 count++;
3027 } // end of loop over icont
3028 } // end of loop over haloed nodes
3029
3030 // Now send all the required info to the equivalent halo layer -
3031 // If there are no discrepancies, no need to send anything
3032 Vector<unsigned> n_all_send(2, 0);
3033 if (discrepancy == 0)
3034 {
3035 MPI_Send(&n_all_send[0], 2, MPI_UNSIGNED, d, 0, Comm_pt->mpi_comm());
3036 }
3037 else
3038 {
3039 // How much data is there to be sent?
3040 n_all_send[0] = send_data.size();
3041 n_all_send[1] = send_double_data.size();
3042 MPI_Send(&n_all_send[0], 2, MPI_UNSIGNED, d, 0, Comm_pt->mpi_comm());
3043
3044 // Send flat-packed ints
3045 if (n_all_send[0] != 0)
3046 {
3047 MPI_Send(
3048 &send_data[0], n_all_send[0], MPI_INT, d, 1, Comm_pt->mpi_comm());
3049 }
3050
3051 // Send flat-packed double data
3052 if (n_all_send[1] != 0)
3053 {
3054 MPI_Send(&send_double_data[0],
3055 n_all_send[1],
3056 MPI_DOUBLE,
3057 d,
3058 1,
3059 Comm_pt->mpi_comm());
3060 }
3061 }
3062 }
3063 else // (d==my_rank), current process
3064 {
3065 // Receive the master nodes and weights in order to modify the
3066 // hanging status of nodes in the halo layer
3067 for (int dd = 0; dd < n_proc; dd++)
3068 {
3069 if (dd != d) // don't talk to yourself
3070 {
3071 // Are we going to receive anything? This is zero
3072 // either if there are no discrepancies or there's zero
3073 // data to be sent.
3074 Vector<unsigned> n_all_recv(2, 0);
3075 MPI_Recv(&n_all_recv[0],
3076 2,
3077 MPI_UNSIGNED,
3078 dd,
3079 0,
3080 Comm_pt->mpi_comm(),
3081 &status);
3082
3083 // Storage for received information
3084 Vector<int> receive_data;
3085 Vector<double> receive_double_data;
3086
3087 // Receive unsigneds (if any)
3088 if (n_all_recv[0] != 0)
3089 {
3090 // Receive the data
3091 receive_data.resize(n_all_recv[0]);
3092 MPI_Recv(&receive_data[0],
3093 n_all_recv[0],
3094 MPI_INT,
3095 dd,
3096 1,
3097 Comm_pt->mpi_comm(),
3098 &status);
3099 }
3100
3101 // Receive doubles (if any)
3102 if (n_all_recv[1] != 0)
3103 {
3104 // Receive the data
3105 receive_double_data.resize(n_all_recv[1]);
3106 MPI_Recv(&receive_double_data[0],
3107 n_all_recv[1],
3108 MPI_DOUBLE,
3109 dd,
3110 1,
3111 Comm_pt->mpi_comm(),
3112 &status);
3113 }
3114
3115 // If no information, no need to do anything else
3116 if (n_all_recv[0] != 0)
3117 {
3118 // Counters for traversing received data
3119 unsigned count = 0;
3120 unsigned count_double = 0;
3121
3122 // Loop over halo nodes
3123 unsigned nh = nhalo_node(dd);
3124 for (unsigned j = 0; j < nh; j++)
3125 {
3126 // Get node
3127 Node* nod_pt = halo_node_pt(dd, j);
3128
3129 // Loop over the hanging status for each interpolated variable
3130 // (and the geometry)
3131 for (int icont = -1; icont < ncont_inter_values; icont++)
3132 {
3133 // Read next entry
3134 int next_entry = receive_data[count++];
3135
3136 // If it's positive, then the number tells us how
3137 // many master nodes we have
3138 if (next_entry > 0)
3139 {
3140 unsigned nhd_master = unsigned(next_entry);
3141
3142 // Set up a new HangInfo for this node
3143 HangInfo* hang_pt = new HangInfo(nhd_master);
3144
3145 // Now set up the master nodes and weights
3146 for (unsigned m = 0; m < nhd_master; m++)
3147 {
3148 // Get the sent master node (a shared node) and
3149 // the weight
3150
3151 // ID of proc in whose shared node lookup scheme
3152 // the sending processor found the node
3153 unsigned shared_node_proc =
3154 unsigned(receive_data[count++]);
3155
3156 // Index of node in the shared node lookup scheme
3157 unsigned shared_node_id = unsigned(receive_data[count++]);
3158
3159 // Get weight
3160 double mtr_weight = receive_double_data[count_double++];
3161
3162 // If the shared node processor is the same as the
3163 // the sending processor we can processor everything here
3164 if (shared_node_proc == unsigned(dd))
3165 {
3166 // Get node
3167 Node* master_nod_pt =
3168 shared_node_pt(dd, shared_node_id);
3169
3170 // Set as a master node (with corresponding weight)
3171 hang_pt->set_master_node_pt(
3172 m, master_nod_pt, mtr_weight);
3173 }
3174 // ...otherwise we have do another communication with
3175 // intermediate processor that holds the non-halo
3176 // version of the master node -- only that processor can
3177 // translate the index of the node the share node
3178 // lookup scheme with the sending processor to the
3179 // index in the shared node lookup scheme with this
3180 // processor
3181 else
3182 {
3183 // Store
3184 HangHelperStruct tmp;
3185 tmp.Sending_processor = dd;
3187 shared_node_id;
3188 tmp.Shared_node_proc = shared_node_proc;
3189 tmp.Weight = mtr_weight;
3190 tmp.Hang_pt = hang_pt;
3191 tmp.Master_node_index = m;
3192 tmp.Node_pt = nod_pt;
3193 tmp.icont = icont;
3194 hang_info.push_back(tmp);
3195 }
3196 }
3197
3198 // Set the hanging pointer for the current halo node
3199 // (does delete any previous hang data)
3200 nod_pt->set_hanging_pt(hang_pt, icont);
3201 }
3202 // Negative entry: the hanging node already exists,
3203 // but it shouldn't, so set it to nonhanging
3204 else if (next_entry < 0)
3205 {
3206 nod_pt->set_hanging_pt(0, icont);
3207 }
3208
3209 } // end of loop over icont
3210 } // end of loop over nodes
3211 } // end of anything to receive
3212 }
3213 }
3214 }
3215 } // end loop over all processors
3216
3217
3219 {
3220 t_end = TimingHelpers::timer();
3221 oomph_info << "Time for second all-to-all in synchronise_hanging_nodes() "
3222 << t_end - t_start << std::endl;
3223 t_start = TimingHelpers::timer();
3224 }
3225
3226
3227 // Now identify master nodes by translating index in shared
3228 // node lookup scheme from the lookup scheme with the sending
3229 // processor to that with the current processor
3230 unsigned n = hang_info.size();
3231
3232 // Is there anything to do be done?
3233 unsigned global_n = 0;
3234 MPI_Allreduce(&n, &global_n, 1, MPI_UNSIGNED, MPI_MAX, Comm_pt->mpi_comm());
3235 if (global_n == 0)
3236 {
3238 << "No need for reconciliation of wrongly synchronised hang nodes\n";
3239 }
3240 // Reconcilation required
3241 else
3242 {
3243 oomph_info << "Need to reconcile of wrongly syncronised hang nodes\n";
3244
3245 // Storage for the translated entries (in order) for/from
3246 // other processors
3247 // Must be ints so that an entry of -1 tells the other processor
3248 // that the node could not be found. See comment above for why
3249 // this may be necessary.
3250 Vector<Vector<int>> translated_entry(n_proc);
3251
3252 // Storage for how-many-th entry in this processor's
3253 // hang_info vector will be completed by processor rank.
3254 Vector<Vector<unsigned>> hang_info_index_for_proc(n_proc);
3255
3256 // Send information to intermediate processor that holds
3257 // non-halo version of missing master node
3258 {
3259 // Storage for number of data to be sent to each processor
3260 Vector<int> send_n(n_proc, 0);
3261
3262 // Storage for all values to be sent to all processors
3263 Vector<int> send_data;
3264
3265 // Start location within send_data for data to be sent to each processor
3266 Vector<int> send_displacement(n_proc, 0);
3267
3268 // Loop over all processors
3269 for (int rank = 0; rank < n_proc; rank++)
3270 {
3271 // Set the offset for the current processor
3272 send_displacement[rank] = send_data.size();
3273
3274 // Don't bother to do anything if the processor in the loop is the
3275 // current processor
3276 if (rank != my_rank)
3277 {
3278 // Search through the (typically few) entries
3279 // in hang info vector to find the ones that
3280 // must be sent to proc rank (the proc that holds
3281 // the non-halo version of the "missing" master node
3282 for (unsigned i = 0; i < n; i++)
3283 {
3284 HangHelperStruct tmp = hang_info[i];
3285 if (tmp.Shared_node_proc == unsigned(rank))
3286 {
3287 // Add the sending processor
3288 send_data.push_back(tmp.Sending_processor);
3289
3290 // Add the index of the missing master node
3291 // in the shared node lookup scheme between
3292 // sending processor and processor rank
3293 send_data.push_back(tmp.Shared_node_id_on_sending_processor);
3294
3295 // Record the how-many-th entry in this processor's
3296 // hang_info vector will be completed by processor rank.
3297 hang_info_index_for_proc[rank].push_back(i);
3298 }
3299 }
3300 }
3301
3302 // Find the number of data added to the vector
3303 send_n[rank] = send_data.size() - send_displacement[rank];
3304 }
3305
3306
3307 // Storage for the number of data to be received from each processor
3308 Vector<int> receive_n(n_proc, 0);
3309
3310 // Now send numbers of data to be sent between all processors
3311 MPI_Alltoall(&send_n[0],
3312 1,
3313 MPI_INT,
3314 &receive_n[0],
3315 1,
3316 MPI_INT,
3317 Comm_pt->mpi_comm());
3318
3319 // We now prepare the data to be received
3320 // by working out the displacements from the received data
3321 Vector<int> receive_displacement(n_proc, 0);
3322 int receive_data_count = 0;
3323 for (int rank = 0; rank < n_proc; ++rank)
3324 {
3325 // Displacement is number of data received so far
3326 receive_displacement[rank] = receive_data_count;
3327 receive_data_count += receive_n[rank];
3328 }
3329
3330 // Now resize the receive buffer for all data from all processors
3331 // Make sure that it has a size of at least one
3332 if (receive_data_count == 0)
3333 {
3334 ++receive_data_count;
3335 }
3336 Vector<unsigned> receive_data(receive_data_count);
3337
3338 // Make sure that the send buffer has size at least one
3339 // so that we don't get a segmentation fault
3340 if (send_data.size() == 0)
3341 {
3342 send_data.resize(1);
3343 }
3344
3345 // Now send the data between all the processors
3346 MPI_Alltoallv(&send_data[0],
3347 &send_n[0],
3348 &send_displacement[0],
3349 MPI_INT,
3350 &receive_data[0],
3351 &receive_n[0],
3352 &receive_displacement[0],
3353 MPI_INT,
3354 Comm_pt->mpi_comm());
3355
3356 // Now use the received data to update the halo nodes
3357 for (int send_rank = 0; send_rank < n_proc; send_rank++)
3358 {
3359 // Don't bother to do anything for the processor corresponding to the
3360 // current processor or if no data were received from this processor
3361 if ((send_rank != my_rank) && (receive_n[send_rank] != 0))
3362 {
3363 // Counter for the data within the large array
3364 unsigned count = receive_displacement[send_rank];
3365
3366 // We're reading two numbers per missing halo node
3367 unsigned n_rec = unsigned(receive_n[send_rank]);
3368 for (unsigned i = 0; i < n_rec / 2; i++)
3369 {
3370 // Receive orig sending proc
3371 unsigned orig_sending_proc = receive_data[count];
3372 count++;
3373
3374 // Receive the index of the missing master node
3375 // in the shared node lookup scheme between
3376 // orig sending processor and current processor
3377 unsigned shared_node_id_on_orig_sending_proc =
3378 receive_data[count];
3379 count++;
3380
3381 // Extract node from shared node lookup scheme
3382 Node* master_nod_pt = shared_node_pt(
3383 orig_sending_proc, shared_node_id_on_orig_sending_proc);
3384
3385 // Now find it in shared halo scheme with the processor
3386 // that's sent the request
3387
3388 std::map<Node*, unsigned>::iterator it =
3389 shared_node_map[send_rank].find(master_nod_pt);
3390
3391 // If it's not in there iterator points to end of
3392 // set
3393 if (it != shared_node_map[send_rank].end())
3394 {
3395 // Store it so we can send it back
3396 translated_entry[send_rank].push_back((*it).second);
3397 }
3398 else
3399 {
3400 // This node has not been found in the shared scheme, so
3401 // the translation query has failed. We send a -1 to tell
3402 // the other processor the bad news.
3403 translated_entry[send_rank].push_back(-1);
3404
3405 /*
3406 // We don't need to crash anymore because the function
3407 // additional_synchronise_hanging_nodes() will magically
3408 // sort out all the problems!
3409 std::ostringstream error_stream;
3410 error_stream
3411 << "Received translation query for shared node"
3412 << " entry " << shared_node_id_on_orig_sending_proc
3413 << " with processor " << orig_sending_proc
3414 << " from proc " << send_rank << std::endl
3415 << "but did not find node in shared node scheme with proc "
3416 << send_rank << std::endl;
3417 throw OomphLibError(
3418 error_stream.str(),
3419 OOMPH_CURRENT_FUNCTION,
3420 OOMPH_EXCEPTION_LOCATION);
3421 */
3422 }
3423 }
3424 }
3425 } // End of data is received
3426
3427 } // end of sending stuff to intermediate processor that holds
3428 // non halo version of missing master node
3429
3430
3432 {
3433 t_end = TimingHelpers::timer();
3435 << "Time for third all-to-all in synchronise_hanging_nodes() "
3436 << t_end - t_start << std::endl;
3437 t_start = TimingHelpers::timer();
3438 }
3439
3440
3441 // Send information back to processor that needs to identify
3442 // missing master node via shared node lookup scheme with
3443 // this processor
3444 {
3445 // Storage for number of data to be sent to each processor
3446 Vector<int> send_n(n_proc, 0);
3447
3448 // Storage for all values to be sent to all processors
3449 Vector<int> send_data;
3450
3451 // Start location within send_data for data to be sent to each processor
3452 Vector<int> send_displacement(n_proc, 0);
3453
3454 // Loop over all processors
3455 for (int rank = 0; rank < n_proc; rank++)
3456 {
3457 // Set the offset for the current processor
3458 send_displacement[rank] = send_data.size();
3459
3460 // Don't bother to do anything if the processor in the loop is the
3461 // current processor
3462 if (rank != my_rank)
3463 {
3464 // Put the translated entries for processor rank into
3465 // send data
3466 unsigned n = translated_entry[rank].size();
3467 for (unsigned j = 0; j < n; j++)
3468 {
3469 send_data.push_back(translated_entry[rank][j]);
3470 }
3471 }
3472
3473 // Find the number of data added to the vector
3474 send_n[rank] = send_data.size() - send_displacement[rank];
3475 }
3476
3477
3478 // Storage for the number of data to be received from each processor
3479 Vector<int> receive_n(n_proc, 0);
3480
3481 // Now send numbers of data to be sent between all processors
3482 MPI_Alltoall(&send_n[0],
3483 1,
3484 MPI_INT,
3485 &receive_n[0],
3486 1,
3487 MPI_INT,
3488 Comm_pt->mpi_comm());
3489
3490 // We now prepare the data to be received
3491 // by working out the displacements from the received data
3492 Vector<int> receive_displacement(n_proc, 0);
3493 int receive_data_count = 0;
3494 for (int rank = 0; rank < n_proc; ++rank)
3495 {
3496 // Displacement is number of data received so far
3497 receive_displacement[rank] = receive_data_count;
3498 receive_data_count += receive_n[rank];
3499 }
3500
3501 // Now resize the receive buffer for all data from all processors
3502 // Make sure that it has a size of at least one
3503 if (receive_data_count == 0)
3504 {
3505 ++receive_data_count;
3506 }
3507 Vector<unsigned> receive_data(receive_data_count);
3508
3509 // Make sure that the send buffer has size at least one
3510 // so that we don't get a segmentation fault
3511 if (send_data.size() == 0)
3512 {
3513 send_data.resize(1);
3514 }
3515
3516 // Now send the data between all the processors
3517 MPI_Alltoallv(&send_data[0],
3518 &send_n[0],
3519 &send_displacement[0],
3520 MPI_INT,
3521 &receive_data[0],
3522 &receive_n[0],
3523 &receive_displacement[0],
3524 MPI_INT,
3525 Comm_pt->mpi_comm());
3526
3527 // Now use the received data to update the halo nodes
3528 for (int send_rank = 0; send_rank < n_proc; send_rank++)
3529 {
3530 // Don't bother to do anything for the processor corresponding to the
3531 // current processor or if no data were received from this processor
3532 if ((send_rank != my_rank) && (receive_n[send_rank] != 0))
3533 {
3534 // Counter for the data within the large array
3535 unsigned count = receive_displacement[send_rank];
3536
3537 // We're reading one number per missing halo node
3538 unsigned n_rec = unsigned(receive_n[send_rank]);
3539 for (unsigned i = 0; i < n_rec; i++)
3540 {
3541 // Index of missing master node in shared node lookup scheme
3542 // with processor send_rank:
3543 // Must be an int because failure returns -1
3544 int index = receive_data[count];
3545 count++;
3546
3547 // Translation query has been successful if index >= 0
3548 if (index >= 0)
3549 {
3550 // Recall information associated with that missing master
3551 unsigned hang_info_index =
3552 hang_info_index_for_proc[send_rank][i];
3553 HangHelperStruct tmp = hang_info[hang_info_index];
3554
3555 // Extract node from shared node lookup scheme
3556 Node* master_nod_pt = shared_node_pt(send_rank, index);
3557
3558 // Set as a master node (with corresponding weight)
3560 tmp.Master_node_index, master_nod_pt, tmp.Weight);
3561 }
3562 else
3563 {
3564 // Translation query has failed. This is the processor
3565 // on which the node was a halo, so we must delete the
3566 // partial hang info.
3567
3568 // Recall information associated with that missing master
3569 unsigned hang_info_index =
3570 hang_info_index_for_proc[send_rank][i];
3571 HangHelperStruct tmp = hang_info[hang_info_index];
3572
3573 // Delete partial hanging information
3574 tmp.Node_pt->set_hanging_pt(0, tmp.icont);
3575
3576 // Set flag to trigger another round of synchronisation
3577 // This works even though we don't own the node that
3578 // still requires synchrionisation because this variable
3579 // is reduced over all processors at the end
3580 nnode_still_requiring_synchronisation++;
3581 }
3582 }
3583 }
3584 } // End of data is received
3585
3586 } // end of completing hang info for missing master nodes
3587
3588
3590 {
3591 t_end = TimingHelpers::timer();
3593 << "Time for fourth all-to-all in synchronise_hanging_nodes() "
3594 << t_end - t_start << std::endl;
3595 }
3596 } // end of reconciliation required
3597
3598
3599 // Get global number of nodes still requiring synchronisation due to
3600 // missing master nodes
3601 // This will only be necessary for meshes involving elements
3602 // with nonuniformly spaced nodes. All other cases will continue to
3603 // work as before because all nodes will have been successfully
3604 // synchronised by now
3605 unsigned global_nnode_still_requiring_synchronisation = 0;
3606 MPI_Allreduce(&nnode_still_requiring_synchronisation,
3607 &global_nnode_still_requiring_synchronisation,
3608 1,
3609 MPI_UNSIGNED,
3610 MPI_MAX,
3611 Comm_pt->mpi_comm());
3612 if (global_nnode_still_requiring_synchronisation > 0)
3613 {
3614 double tt_start = 0.0;
3616 {
3617 tt_start = TimingHelpers::timer();
3618 }
3619
3620 oomph_info << "Need to do additional synchronisation of hanging nodes"
3621 << std::endl;
3622
3623 // Do additional synchronisation
3624 additional_synchronise_hanging_nodes(ncont_interpolated_values);
3625
3626 double tt_end = 0.0;
3628 {
3629 tt_end = TimingHelpers::timer();
3631 << "Time for RefineableMesh::additional_synchronise_hanging_nodes() "
3632 << "in TreeBasedRefineableMeshBase::synchronise_hanging_nodes(): "
3633 << tt_end - tt_start << std::endl;
3634 tt_start = TimingHelpers::timer();
3635 }
3636 }
3637 else
3638 {
3639 oomph_info << "No need to do additional synchronisation of hanging nodes"
3640 << std::endl;
3641 }
3642 }
3643
3644 //========================================================================
3645 /// Synchronise the positions of non-hanging nodes that depend on
3646 /// non-existent neighbours (e.g. h-refinement of neighbouring elements
3647 /// with different p-orders where the shared edge is on the outer edge of
3648 /// the halo layer)
3649 //========================================================================
3651 {
3652 // Store number of processors and current process
3653 MPI_Status status;
3654 int n_proc = Comm_pt->nproc();
3655 int my_rank = Comm_pt->my_rank();
3656
3657 double t_start = 0.0;
3658 double t_end = 0.0;
3659
3660 // Storage for the hanging status of halo/haloed nodes on elements
3661 Vector<Vector<unsigned>> recv_unsigneds(n_proc);
3662 Vector<Vector<double>> recv_doubles(n_proc);
3663
3665 {
3666 t_start = TimingHelpers::timer();
3667 }
3668
3669 // Loop over processes: Each processor checks if its nonhaning nodes in
3670 // haloed elements with proc d require additional information to determine
3671 // their positions.
3672 for (int d = 0; d < n_proc; d++)
3673 {
3674 // No halo with self: Setup hang info for my haloed nodes with proc d
3675 // then get ready to receive halo info from processor d.
3676 if (d != my_rank)
3677 {
3678 // Receive the position information from the corresponding process
3679 unsigned recv_unsigneds_count = 0;
3680 MPI_Recv(&recv_unsigneds_count,
3681 1,
3682 MPI_UNSIGNED,
3683 d,
3684 0,
3685 Comm_pt->mpi_comm(),
3686 &status);
3687 unsigned recv_doubles_count = 0;
3688 MPI_Recv(&recv_doubles_count,
3689 1,
3690 MPI_UNSIGNED,
3691 d,
3692 1,
3693 Comm_pt->mpi_comm(),
3694 &status);
3695
3696 // Get the data (if any)
3697 if (recv_unsigneds_count != 0)
3698 {
3699 recv_unsigneds[d].resize(recv_unsigneds_count);
3700 MPI_Recv(&recv_unsigneds[d][0],
3701 recv_unsigneds_count,
3702 MPI_UNSIGNED,
3703 d,
3704 0,
3705 Comm_pt->mpi_comm(),
3706 &status);
3707 }
3708 if (recv_doubles_count != 0)
3709 {
3710 recv_doubles[d].resize(recv_doubles_count);
3711 MPI_Recv(&recv_doubles[d][0],
3712 recv_doubles_count,
3713 MPI_DOUBLE,
3714 d,
3715 1,
3716 Comm_pt->mpi_comm(),
3717 &status);
3718 }
3719
3720 // Counters for received data
3721 unsigned recv_unsigneds_index = 0;
3722 double recv_doubles_index = 0;
3723
3724 // Get halo elements with processor d
3726
3727 // Loop over recieved indices
3728 while (recv_unsigneds_index < recv_unsigneds_count)
3729 {
3730 // Get (finite) element
3731 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(
3732 halo_element_pt[recv_unsigneds[d][recv_unsigneds_index++]]);
3733
3734 // If we have a finite element...
3735 if (el_pt != 0)
3736 {
3737 // Get dimension
3738 unsigned n_dim = el_pt->dim();
3739
3740 // Get node
3741 Node* nod_pt =
3742 el_pt->node_pt(recv_unsigneds[d][recv_unsigneds_index++]);
3743
3744 // Get current position
3745 Vector<double> x_cur(n_dim);
3746 for (unsigned dir = 0; dir < n_dim; dir++)
3747 {
3748 x_cur[dir] = nod_pt->x(dir);
3749 }
3750
3751 // Get recieved position
3752 Vector<double> x_rec(n_dim);
3753 for (unsigned dir = 0; dir < n_dim; dir++)
3754 {
3755 x_rec[dir] = recv_doubles[d][recv_doubles_index + dir];
3756 }
3757
3758 // Compare actual and expected positions
3759 bool node_pos_differs = false;
3760 for (unsigned dir = 0; dir < n_dim; dir++)
3761 {
3762 node_pos_differs = node_pos_differs ||
3763 (std::fabs(x_cur[dir] - x_rec[dir]) > 1.0e-14);
3764 }
3765
3766 // Set the actual position
3767 Vector<double> x_act(n_dim);
3768 for (unsigned dir = 0; dir < n_dim; dir++)
3769 {
3770 nod_pt->x(dir) = recv_doubles[d][recv_doubles_index++];
3771 }
3772 }
3773 }
3774
3775 if (recv_unsigneds_count != recv_unsigneds_index)
3776 {
3777 std::ostringstream error_stream;
3778 error_stream << "recv_unsigneds_count != recv_unsigneds_index ( "
3779 << recv_unsigneds_count << " != " << recv_unsigneds_index
3780 << ")" << std::endl;
3781 throw OomphLibError(
3782 error_stream.str(),
3783 "TreeBasedRefineableMeshBase::synchronise_nonhanging_nodes()",
3784 OOMPH_EXCEPTION_LOCATION);
3785 }
3786 }
3787 else // d==my_rank, i.e. current process: Send halo hanging status
3788 // to process dd where it's received (see above) and compared
3789 // and compared against the hang status of the haloed nodes
3790 {
3791 for (int dd = 0; dd < n_proc; dd++)
3792 {
3793 // No halo with yourself
3794 if (dd != d)
3795 {
3796 // Storage for halo hanging status and counter
3797 Vector<int> send_unsigneds;
3798 Vector<double> send_doubles;
3799
3800 // Set to store nodes whose position requires adjustment
3801 std::set<Node*> nodes_requiring_adjustment;
3802
3803 // Get haloed elements with processor dd
3805 this->haloed_element_pt(dd));
3806
3807 // Loop over haloed elements with processor dd
3808 unsigned nh = haloed_element_pt.size();
3809 for (unsigned e = 0; e < nh; e++)
3810 {
3811 // Get (finite) element
3812 FiniteElement* el_pt =
3813 dynamic_cast<FiniteElement*>(haloed_element_pt[e]);
3814
3815 // If we have a finite element...
3816 if (el_pt != 0)
3817 {
3818 // Get dimension
3819 unsigned n_dim = el_pt->dim();
3820
3821 // Loop over element nodes
3822 unsigned n_node = el_pt->nnode();
3823 for (unsigned j = 0; j < n_node; j++)
3824 {
3825 // Get node
3826 Node* nod_pt = el_pt->node_pt(j);
3827
3828 // Only do non-hanging nodes
3829 if (!nod_pt->is_hanging())
3830 {
3831 // Check if node's position is the same as that interpolated
3832 // using its local coordinate in the haloed element
3833
3834 // Loop over all history values
3835 unsigned nt = nod_pt->ntstorage();
3836 for (unsigned t = 0; t < nt; t++)
3837 {
3838 // Get expected position
3839 Vector<double> s(n_dim), x_exp(n_dim);
3840 el_pt->local_coordinate_of_node(j, s);
3841 el_pt->get_x(t, s, x_exp);
3842
3843 // Get actual position
3844 Vector<double> x_act(n_dim);
3845 for (unsigned dir = 0; dir < n_dim; dir++)
3846 {
3847 x_act[dir] = nod_pt->x(dir);
3848 }
3849
3850 // Compare actual and expected positions
3851 bool node_pos_differs = false;
3852 for (unsigned dir = 0; dir < n_dim; dir++)
3853 {
3854 node_pos_differs =
3855 node_pos_differs ||
3856 (std::fabs(x_act[dir] - x_exp[dir]) > 1.0e-14);
3857 }
3858
3859 // If the node's actual position differs from its
3860 // expected position we need to communicate this
3861 // information to processors on which this is a halo node
3862 if (node_pos_differs)
3863 {
3864 // Check that node has not been done already
3865 if (nodes_requiring_adjustment.insert(nod_pt).second)
3866 {
3867 // Send index of haloed element
3868 send_unsigneds.push_back(e);
3869 // Send index of node in the element
3870 send_unsigneds.push_back(j);
3871 // Send actual position of node
3872 for (unsigned dir = 0; dir < n_dim; dir++)
3873 {
3874 send_doubles.push_back(x_act[dir]);
3875 }
3876 }
3877 }
3878 }
3879 }
3880 }
3881 }
3882 }
3883
3884 // Send the information to the relevant process
3885 unsigned send_unsigneds_count = send_unsigneds.size();
3886 unsigned send_doubles_count = send_doubles.size();
3887
3888 if (send_unsigneds_count > 0)
3889 {
3890 // exit(1);
3891 }
3892
3893 // Tell processor dd how much data to receive
3894 MPI_Send(&send_unsigneds_count,
3895 1,
3896 MPI_UNSIGNED,
3897 dd,
3898 0,
3899 Comm_pt->mpi_comm());
3900 MPI_Send(
3901 &send_doubles_count, 1, MPI_UNSIGNED, dd, 1, Comm_pt->mpi_comm());
3902
3903 // Send data (if any)
3904 if (send_unsigneds_count != 0)
3905 {
3906 MPI_Send(&send_unsigneds[0],
3907 send_unsigneds_count,
3908 MPI_UNSIGNED,
3909 dd,
3910 0,
3911 Comm_pt->mpi_comm());
3912 }
3913 if (send_doubles_count != 0)
3914 {
3915 MPI_Send(&send_doubles[0],
3916 send_doubles_count,
3917 MPI_DOUBLE,
3918 dd,
3919 1,
3920 Comm_pt->mpi_comm());
3921 }
3922 }
3923 }
3924 }
3925 }
3926
3928 {
3929 t_end = TimingHelpers::timer();
3930 oomph_info << "Time for synchronise_nonhanging_nodes(): "
3931 << t_end - t_start << std::endl;
3932 t_start = TimingHelpers::timer();
3933 }
3934 }
3935
3936#endif
3937
3938
3939 /// /////////////////////////////////////////////////////////////////
3940 /// /////////////////////////////////////////////////////////////////
3941 /// /////////////////////////////////////////////////////////////////
3942
3943
3944 //========================================================================
3945 /// Do adaptive p-refinement for mesh.
3946 /// - Pass Vector of error estimates for all elements.
3947 /// - p-refine those whose errors exceeds the threshold
3948 /// - p-unrefine those whose errors is less than
3949 /// threshold.
3950 //========================================================================
3952 const Vector<double>& elemental_error)
3953 {
3954 // Set the refinement tolerance to be the max permissible error
3955 double refine_tol = this->max_permitted_error();
3956
3957 // Set the unrefinement tolerance to be the min permissible error
3958 double unrefine_tol = this->min_permitted_error();
3959
3960 // Setup doc info
3961 DocInfo local_doc_info;
3962 if (doc_info_pt() == 0)
3963 {
3964 local_doc_info.disable_doc();
3965 }
3966 else
3967 {
3968 local_doc_info = this->doc_info();
3969 }
3970
3971
3972 // Check that the errors make sense
3973 if (refine_tol <= unrefine_tol)
3974 {
3975 std::ostringstream error_stream;
3976 error_stream << "Refinement tolerance <= Unrefinement tolerance"
3977 << refine_tol << " " << unrefine_tol << std::endl
3978 << "doesn't make sense and will almost certainly crash"
3979 << std::endl
3980 << "this beautiful code!" << std::endl;
3981
3982 throw OomphLibError(
3983 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3984 }
3985
3986
3987 // Select elements for refinement and unrefinement
3988 //==============================================
3989 // Reset counter for number of elements that would like to be
3990 // refined further but can't
3991 this->nrefinement_overruled() = 0;
3992
3993 unsigned n_refine = 0;
3994 unsigned n_unrefine = 0;
3995 // Loop over all elements and mark them according to the error criterion
3996 unsigned long Nelement = this->nelement();
3997 for (unsigned long e = 0; e < Nelement; e++)
3998 {
3999 //(Cast) pointer to the element
4000 PRefineableElement* el_pt =
4001 dynamic_cast<PRefineableElement*>(this->element_pt(e));
4002
4003 // Check that we can p-refine the element
4004 if (el_pt != 0)
4005 {
4006 // Initially element is not to be refined
4009
4010 // If the element error exceeds the threshold ...
4011 if (elemental_error[e] > refine_tol)
4012 {
4013 // ... and its refinement level is less than the maximum desired level
4014 // mark is to be refined
4015 if ((el_pt->p_refinement_is_enabled()) &&
4016 (el_pt->p_order() < this->max_p_refinement_level()))
4017 {
4018 el_pt->select_for_p_refinement();
4019 n_refine++;
4020 }
4021 // ... otherwise mark it as having been over-ruled
4022 else
4023 {
4024 this->nrefinement_overruled() += 1;
4025 }
4026 }
4027 if (elemental_error[e] < unrefine_tol)
4028 {
4029 // ... and its refinement level is more than the minimum desired level
4030 // mark is to be refined
4031 //(Also check that we don't unrefine past the initial refinement
4032 // level)
4033 if ((el_pt->p_refinement_is_enabled()) &&
4034 (el_pt->p_order() > this->min_p_refinement_level()) &&
4035 (el_pt->p_order() > el_pt->initial_p_order()))
4036 {
4038 n_unrefine++;
4039 }
4040 // Don't mark as overruled - it's misleading
4041 /// / ... otherwise mark it as having been over-ruled
4042 // else
4043 // {
4044 // this->nrefinement_overruled()+=1;
4045 // }
4046 }
4047 } // End of check for p-refineability of element
4048 else
4049 {
4050 oomph_info << "p-refinement is not possible for these elements"
4051 << std::endl;
4052 // Don't try to p-refine any more elements
4053 break;
4054 }
4055 }
4056
4057 oomph_info << " \n Number of elements to be refined: " << n_refine
4058 << std::endl;
4059 oomph_info << " \n Number of elements whose refinement was overruled: "
4060 << this->nrefinement_overruled() << std::endl;
4061
4062 oomph_info << " \n Number of elements to be unrefined : " << n_unrefine
4063 << std::endl
4064 << std::endl;
4065
4066
4067 // Now do the actual mesh adaptation
4068 //---------------------------------
4069
4070 // Check whether its worth our while
4071 // Either some elements want to be refined,
4072 // or the number that want to be unrefined are greater than the
4073 // specified tolerance
4074
4075 // In a parallel job, it is possible that one process may not have
4076 // any elements to refine, BUT a neighbouring process may refine an
4077 // element which changes the hanging status of a node that is on
4078 // both processes (i.e. a halo(ed) node). To get around this issue,
4079 // ALL processes need to call adapt_mesh if ANY refinement is to
4080 // take place anywhere.
4081
4082 unsigned total_n_refine = 0;
4083#ifdef OOMPH_HAS_MPI
4084 // Sum n_refine across all processors
4085 if (this->is_mesh_distributed())
4086 {
4087 MPI_Allreduce(
4088 &n_refine, &total_n_refine, 1, MPI_INT, MPI_SUM, Comm_pt->mpi_comm());
4089 }
4090 else
4091 {
4092 total_n_refine = n_refine;
4093 }
4094#else
4095 total_n_refine = n_refine;
4096#endif
4097
4098 // There may be some issues with unrefinement too, but I have not
4099 // been able to come up with an example (either in my head or in a
4100 // particular problem) where anything has arisen. I can see that
4101 // there may be an issue if n_unrefine differs across processes so
4102 // that (total_n_unrefine > max_keep_unrefined()) on some but not
4103 // all processes. I haven't seen any examples of this yet so the
4104 // following code may or may not work! (Andy, 06/03/08)
4105
4106 unsigned total_n_unrefine = 0;
4107#ifdef OOMPH_HAS_MPI
4108 // Sum n_unrefine across all processors
4109 if (this->is_mesh_distributed())
4110 {
4111 MPI_Allreduce(&n_unrefine,
4112 &total_n_unrefine,
4113 1,
4114 MPI_INT,
4115 MPI_SUM,
4116 Comm_pt->mpi_comm());
4117 }
4118 else
4119 {
4120 total_n_unrefine = n_unrefine;
4121 }
4122#else
4123 total_n_unrefine = n_unrefine;
4124#endif
4125
4126 oomph_info << "---> " << total_n_refine << " elements to be refined, and "
4127 << total_n_unrefine << " to be unrefined, in total."
4128 << std::endl;
4129
4130 if ((total_n_refine > 0) || (total_n_unrefine > this->max_keep_unrefined()))
4131 {
4132#ifdef PARANOID
4133#ifdef OOMPH_HAS_MPI
4134
4135 // Sanity check: Each processor checks if the enforced unrefinement of
4136 // its haloed element is matched by enforced unrefinement of the
4137 // corresponding halo elements on the other processors.
4138 if (this->is_mesh_distributed())
4139 {
4140 // Store number of processors and current process
4141 MPI_Status status;
4142 int n_proc = Comm_pt->nproc();
4143 int my_rank = Comm_pt->my_rank();
4144
4145 // Loop over all other domains/processors
4146 for (int d = 0; d < n_proc; d++)
4147 {
4148 // Don't talk to yourself
4149 if (d != my_rank)
4150 {
4151 {
4152 // Get the vector of halo elements whose non-halo counterpart
4153 // are on processor d
4154 Vector<GeneralisedElement*> halo_elem_pt(
4155 this->halo_element_pt(d));
4156
4157 // Create vector containing (0)1 to indicate that
4158 // halo element is (not) to be unrefined
4159 unsigned nhalo = halo_elem_pt.size();
4160 Vector<int> halo_to_be_unrefined(nhalo, 0);
4161 for (unsigned e = 0; e < nhalo; e++)
4162 {
4163 if (dynamic_cast<PRefineableElement*>(halo_elem_pt[e])
4164 ->to_be_p_unrefined())
4165 {
4166 halo_to_be_unrefined[e] = 1;
4167 }
4168 }
4169
4170 // Trap the case when there are no halo elements
4171 // so that we don't get a segfault in the MPI send
4172 if (nhalo > 0)
4173 {
4174 // Send it across
4175 MPI_Send(&halo_to_be_unrefined[0],
4176 nhalo,
4177 MPI_INT,
4178 d,
4179 0,
4180 Comm_pt->mpi_comm());
4181 }
4182 }
4183
4184 {
4185 // Get the vector of haloed elements on current processor
4186 Vector<GeneralisedElement*> haloed_elem_pt(
4187 this->haloed_element_pt(d));
4188
4189 // Ask processor d to send vector containing (0)1 for
4190 // halo element with current processor to be (not)unrefined
4191 unsigned nhaloed = haloed_elem_pt.size();
4192 Vector<int> halo_to_be_unrefined(nhaloed);
4193 // Trap to catch the case that there are no haloed elements
4194 if (nhaloed > 0)
4195 {
4196 MPI_Recv(&halo_to_be_unrefined[0],
4197 nhaloed,
4198 MPI_INT,
4199 d,
4200 0,
4201 Comm_pt->mpi_comm(),
4202 &status);
4203 }
4204
4205 // Check it
4206 for (unsigned e = 0; e < nhaloed; e++)
4207 {
4208 if (((halo_to_be_unrefined[e] == 0) &&
4209 (dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])
4210 ->to_be_p_unrefined())) ||
4211 ((halo_to_be_unrefined[e] == 1) &&
4212 (!dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])
4213 ->to_be_p_unrefined())))
4214 {
4215 std::ostringstream error_message;
4216 error_message
4217 << "Error in refinement: \n"
4218 << "Haloed element: " << e << " on proc " << my_rank
4219 << " \n"
4220 << "wants to be unrefined whereas its halo counterpart on\n"
4221 << "proc " << d << " doesn't (or vice versa)...\n"
4222 << "This is most likely because the error estimator\n"
4223 << "has not assigned the same errors to halo and haloed\n"
4224 << "elements -- it ought to!\n";
4225 throw OomphLibError(error_message.str(),
4226 OOMPH_CURRENT_FUNCTION,
4227 OOMPH_EXCEPTION_LOCATION);
4228 }
4229 }
4230 }
4231 }
4232 }
4233
4234
4235 // Loop over all other domains/processors
4236 for (int d = 0; d < n_proc; d++)
4237 {
4238 // Don't talk to yourself
4239 if (d != my_rank)
4240 {
4241 {
4242 // Get the vector of halo elements whose non-halo counterpart
4243 // are on processor d
4244 Vector<GeneralisedElement*> halo_elem_pt(
4245 this->halo_element_pt(d));
4246
4247 // Create vector containing (0)1 to indicate that
4248 // halo element is (not) to be refined
4249 unsigned nhalo = halo_elem_pt.size();
4250 Vector<int> halo_to_be_refined(nhalo, 0);
4251 for (unsigned e = 0; e < nhalo; e++)
4252 {
4253 if (dynamic_cast<PRefineableElement*>(halo_elem_pt[e])
4254 ->to_be_p_refined())
4255 {
4256 halo_to_be_refined[e] = 1;
4257 }
4258 }
4259
4260 // Send it across
4261 if (nhalo > 0)
4262 {
4263 MPI_Send(&halo_to_be_refined[0],
4264 nhalo,
4265 MPI_INT,
4266 d,
4267 0,
4268 Comm_pt->mpi_comm());
4269 }
4270 }
4271
4272 {
4273 // Get the vector of haloed elements on current processor
4274 Vector<GeneralisedElement*> haloed_elem_pt(
4275 this->haloed_element_pt(d));
4276
4277 // Ask processor d to send vector containing (0)1 for
4278 // halo element with current processor to be (not)refined
4279 unsigned nhaloed = haloed_elem_pt.size();
4280 Vector<int> halo_to_be_refined(nhaloed);
4281 if (nhaloed > 0)
4282 {
4283 MPI_Recv(&halo_to_be_refined[0],
4284 nhaloed,
4285 MPI_INT,
4286 d,
4287 0,
4288 Comm_pt->mpi_comm(),
4289 &status);
4290 }
4291
4292 // Check it
4293 for (unsigned e = 0; e < nhaloed; e++)
4294 {
4295 if (((halo_to_be_refined[e] == 0) &&
4296 (dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])
4297 ->to_be_p_refined())) ||
4298 ((halo_to_be_refined[e] == 1) &&
4299 (!dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])
4300 ->to_be_p_refined())))
4301 {
4302 std::ostringstream error_message;
4303 error_message
4304 << "Error in refinement: \n"
4305 << "Haloed element: " << e << " on proc " << my_rank
4306 << " \n"
4307 << "wants to be refined whereas its halo counterpart on\n"
4308 << "proc " << d << " doesn't (or vice versa)...\n"
4309 << "This is most likely because the error estimator\n"
4310 << "has not assigned the same errors to halo and haloed\n"
4311 << "elements -- it ought to!\n";
4312 throw OomphLibError(error_message.str(),
4313 OOMPH_CURRENT_FUNCTION,
4314 OOMPH_EXCEPTION_LOCATION);
4315 }
4316 }
4317 }
4318 }
4319 }
4320 }
4321#endif
4322#endif
4323
4324 // Perform the actual adaptation
4325 p_adapt_mesh(local_doc_info);
4326
4327 // The number of refineable elements is still local to each process
4328 this->Nunrefined = n_unrefine;
4329 this->Nrefined = n_refine;
4330 }
4331 // If not worthwhile, say so but still reorder nodes and kill external
4332 // storage for consistency in parallel computations
4333 else
4334 {
4335#ifdef OOMPH_HAS_MPI
4336 // Delete any external element storage - any interaction will still
4337 // be set up on the fly again, so we need to get rid of old information.
4338 // This particularly causes problems in multi-domain examples where
4339 // we decide not to refine one of the meshes
4341#endif
4342
4343 // Reorder the nodes within the mesh's node vector
4344 // to establish a standard ordering regardless of the sequence
4345 // of mesh refinements -- this is required to allow dump/restart
4346 // on refined meshes
4347 this->reorder_nodes();
4348
4349#ifdef OOMPH_HAS_MPI
4350
4351 // Now (re-)classify halo and haloed nodes and synchronise hanging
4352 // nodes
4353 // This is required in cases where delete_all_external_storage()
4354 // made dependent nodes to external halo nodes nonhanging.
4355 if (this->is_mesh_distributed())
4356 {
4360 }
4361
4362#endif
4363
4364 if (n_refine == 0)
4365 {
4366 oomph_info << "\n Not enough benefit in adapting mesh. " << std::endl
4367 << std::endl;
4368 }
4369 this->Nunrefined = 0;
4370 this->Nrefined = 0;
4371 }
4372 }
4373
4374
4375 //================================================================
4376 /// p-adapt mesh, which exists in two representations,
4377 /// namely as:
4378 /// - a FE mesh
4379 /// - a forest of Oc or QuadTrees
4380 ///
4381 /// p-refinement/derefinement process is documented (in tecplot-able form)
4382 /// if requested.
4383 ///
4384 /// Procedure:
4385 /// - Loop over all elements and do the p-refinement/unrefinement for
4386 /// those who want to be refined. Note: p-refinement builds fully-
4387 /// functional elements.
4388 /// - For all nodes that were hanging on the previous mesh (and are still
4389 /// marked as such), fill in their nodal values (consistent
4390 /// with the current hanging node scheme) to make sure they are fully
4391 /// functional, should they have become non-hanging during the
4392 /// mesh-adaptation. Then mark the nodes as non-hanging.
4393 /// - Delete any nodes that have become obsolete.
4394 /// - Mark up hanging nodes and setup hanging node scheme (incl.
4395 /// recursive cleanup for hanging nodes that depend on other
4396 /// hanging nodes).
4397 /// - run a quick self-test on the neighbour finding scheme and
4398 /// check the integrity of the elements (if PARANOID)
4399 /// - doc hanging node status, boundary conditions, neighbour
4400 /// scheme if requested.
4401 ///
4402 ///
4403 /// After adaptation, all nodes (whether new or old) have up-to-date
4404 /// current and previous values.
4405 ///
4406 /// If refinement process is being documented, the following information
4407 /// is documented:
4408 /// - The files
4409 /// - "neighbours.dat"
4410 /// - "all_nodes.dat"
4411 /// - "new_nodes.dat"
4412 /// - "hang_nodes_*.dat"
4413 /// where the * denotes a direction (n,s,e,w) in 2D
4414 /// or (r,l,u,d,f,b) in 3D
4415 /// .
4416 /// can be viewed with
4417 /// - QHangingNodes.mcr
4418 /// .
4419 /// - The file
4420 /// - "hangnodes_withmasters.dat"
4421 /// .
4422 /// can be viewed with
4423 /// - QHangingNodesWithMasters.mcr
4424 /// .
4425 /// to check the hanging node status.
4426 /// - The neighbour status of the elements is documented in
4427 /// - "neighbours.dat"
4428 /// .
4429 /// and can be viewed with
4430 /// - QuadTreeNeighbours.mcr
4431 /// .
4432 //=================================================================
4434 {
4435#ifdef OOMPH_HAS_MPI
4436 // Delete any external element storage before performing the adaptation
4437 // (in particular, external halo nodes that are on mesh boundaries)
4439#endif
4440
4441 // Only perform the adapt step if the mesh has any elements. This is
4442 // relevant in a distributed problem with multiple meshes, where a
4443 // particular process may not have any elements on a particular submesh.
4444 if (this->nelement() > 0)
4445 {
4446 double t_start = 0.0;
4448 {
4449 t_start = TimingHelpers::timer();
4450 }
4451
4452 // Do refinement/unrefinement if required
4454
4456 {
4457 double t_end = TimingHelpers::timer();
4458 oomph_info << "Time for p-refinement/unrefinement: " << t_end - t_start
4459 << std::endl;
4460 t_start = TimingHelpers::timer();
4461 }
4462
4463 // Loop over all nodes in mesh and free the dofs of those that were
4464 //-----------------------------------------------------------------
4465 // pinned only because they were hanging nodes. Also update their
4466 //-----------------------------------------------------------------
4467 // nodal values so that they contain data that is consistent
4468 //----------------------------------------------------------
4469 // with the hanging node representation
4470 //-------------------------------------
4471 // (Even if the nodal data isn't actually accessed because the node
4472 // is still hanging -- we don't know this yet, and this step makes
4473 // sure that all nodes are fully functional and up-to-date, should
4474 // they become non-hanging below).
4475 //
4476 //
4477 // However, if we have a fixed mesh and hanging nodes on the boundary
4478 // become non-hanging they will not necessarily respect the curvilinear
4479 // boundaries. This can only happen in 3D of course because it is not
4480 // possible to have hanging nodes on boundaries in 2D.
4481 // The solution is to store those nodes on the boundaries that are
4482 // currently hanging and then check to see whether they have changed
4483 // status at the end of the refinement procedure.
4484 // If it has changed, then we need to adjust their positions.
4485 const unsigned n_boundary = this->nboundary();
4486 const unsigned mesh_dim = this->finite_element_pt(0)->dim();
4487 Vector<std::set<Node*>> hanging_nodes_on_boundary_pt(n_boundary);
4488
4489 unsigned long n_node = this->nnode();
4490 for (unsigned long n = 0; n < n_node; n++)
4491 {
4492 // Get the pointer to the node
4493 Node* nod_pt = this->node_pt(n);
4494
4495 // Get the number of values in the node
4496 unsigned n_value = nod_pt->nvalue();
4497
4498 // We need to find if any of the values are hanging
4499 bool is_hanging = nod_pt->is_hanging();
4500 // Loop over the values and find out whether any are hanging
4501 for (unsigned n = 0; n < n_value; n++)
4502 {
4503 is_hanging |= nod_pt->is_hanging(n);
4504 }
4505
4506 // If the node is hanging then ...
4507 if (is_hanging)
4508 {
4509 // Unless they are turned into hanging nodes again below
4510 // (this might or might not happen), fill in all the necessary
4511 // data to make them 'proper' nodes again.
4512
4513 // Reconstruct the nodal values/position from the node's
4514 // hanging node representation
4515 unsigned nt = nod_pt->ntstorage();
4516 Vector<double> values(n_value);
4517 unsigned n_dim = nod_pt->ndim();
4518 Vector<double> position(n_dim);
4519 // Loop over all history values
4520 for (unsigned t = 0; t < nt; t++)
4521 {
4522 nod_pt->value(t, values);
4523 for (unsigned i = 0; i < n_value; i++)
4524 {
4525 nod_pt->set_value(t, i, values[i]);
4526 }
4527 nod_pt->position(t, position);
4528 for (unsigned i = 0; i < n_dim; i++)
4529 {
4530 nod_pt->x(t, i) = position[i];
4531 }
4532 }
4533
4534 // If it's an algebraic node: Update its previous nodal positions too
4535 AlgebraicNode* alg_node_pt = dynamic_cast<AlgebraicNode*>(nod_pt);
4536 if (alg_node_pt != 0)
4537 {
4538 bool update_all_time_levels = true;
4539 alg_node_pt->node_update(update_all_time_levels);
4540 }
4541
4542
4543 // If it's a Solid node, update Lagrangian coordinates
4544 // from its hanging node representation
4545 SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
4546 if (solid_node_pt != 0)
4547 {
4548 unsigned n_lagrangian = solid_node_pt->nlagrangian();
4549 for (unsigned i = 0; i < n_lagrangian; i++)
4550 {
4551 solid_node_pt->xi(i) = solid_node_pt->lagrangian_position(i);
4552 }
4553 }
4554
4555 // Now store geometrically hanging nodes on boundaries that
4556 // may need updating after refinement.
4557 // There will only be a problem if we have 3 spatial dimensions
4558 if ((mesh_dim > 2) && (nod_pt->is_hanging()))
4559 {
4560 // If the node is on a boundary then add a pointer to the node
4561 // to our lookup scheme
4562 if (nod_pt->is_on_boundary())
4563 {
4564 // Storage for the boundaries on which the Node is located
4565 std::set<unsigned>* boundaries_pt;
4566 nod_pt->get_boundaries_pt(boundaries_pt);
4567 if (boundaries_pt != 0)
4568 {
4569 // Loop over the boundaries and add a pointer to the node
4570 // to the appropriate storage scheme
4571 for (std::set<unsigned>::iterator it = boundaries_pt->begin();
4572 it != boundaries_pt->end();
4573 ++it)
4574 {
4575 hanging_nodes_on_boundary_pt[*it].insert(nod_pt);
4576 }
4577 }
4578 }
4579 }
4580
4581 } // End of is_hanging
4582
4583 // Initially mark all nodes as 'non-hanging' and `obsolete'
4584 nod_pt->set_nonhanging();
4585 nod_pt->set_obsolete();
4586 }
4587
4588 double t_end = 0.0;
4590 {
4591 t_end = TimingHelpers::timer();
4592 oomph_info << "Time for sorting out initial hanging status: "
4593 << t_end - t_start << std::endl;
4594 t_start = TimingHelpers::timer();
4595 }
4596
4597 // Stick all elements into a vector
4598 Vector<Tree*> tree_nodes_pt;
4599 this->forest_pt()->stick_leaves_into_vector(tree_nodes_pt);
4600
4601 // Copy the elements into the mesh Vector
4602 unsigned long num_tree_nodes = tree_nodes_pt.size();
4603 this->element_pt().resize(num_tree_nodes);
4604 for (unsigned long e = 0; e < num_tree_nodes; e++)
4605 {
4606 this->element_pt(e) = tree_nodes_pt[e]->object_pt();
4607
4608 // Now loop over all nodes in element and mark them as non-obsolete
4609 FiniteElement* this_el_pt = this->finite_element_pt(e);
4610 unsigned n_node = this_el_pt->nnode(); // caching pre-loop
4611 for (unsigned n = 0; n < n_node; n++)
4612 {
4613 this_el_pt->node_pt(n)->set_non_obsolete();
4614 }
4615
4616 // Mark up so that repeated refinements do not occur
4617 // (Required because refined element is the same element as the
4618 // original)
4619 PRefineableElement* cast_el_pt =
4620 dynamic_cast<PRefineableElement*>(this->element_pt(e));
4621 cast_el_pt->deselect_for_p_refinement();
4622 cast_el_pt->deselect_for_p_unrefinement();
4623 }
4624
4625 // Cannot delete nodes that are still marked as obsolete
4626 // because they may still be required to assemble the hanging schemes
4627 //-------------------------------------------------------------------
4628
4629 // Mark up hanging nodes
4630 //----------------------
4631
4632 // Output streams for the hanging nodes
4633 Vector<std::ofstream*> hanging_output_files;
4634 // Setup the output files for hanging nodes, this must be called
4635 // precisely once for the forest. Note that the files will only
4636 // actually be opened if doc_info.Doc_flag is true
4637 this->forest_pt()->open_hanging_node_files(doc_info,
4638 hanging_output_files);
4639
4640 for (unsigned long e = 0; e < num_tree_nodes; e++)
4641 {
4642 // Generic setup
4643 tree_nodes_pt[e]->object_pt()->setup_hanging_nodes(
4644 hanging_output_files);
4645 // Element specific setup
4646 tree_nodes_pt[e]->object_pt()->further_setup_hanging_nodes();
4647 }
4648
4649 // Close the hanging node files and delete the memory allocated
4650 // for the streams
4651 this->forest_pt()->close_hanging_node_files(doc_info,
4652 hanging_output_files);
4653
4654
4656 {
4657 t_end = TimingHelpers::timer();
4658 oomph_info << "Time for setup_hanging_nodes() and "
4659 "further_setup_hanging_nodes() for "
4660 << num_tree_nodes << " elements: " << t_end - t_start
4661 << std::endl;
4662 t_start = TimingHelpers::timer();
4663 }
4664
4665 // Read out the number of continously interpolated values
4666 // from one of the elements (assuming it's the same in all elements)
4667 unsigned ncont_interpolated_values =
4668 tree_nodes_pt[0]->object_pt()->ncont_interpolated_values();
4669
4670 // Complete the hanging nodes schemes by dealing with the
4671 // recursively hanging nodes
4672 this->complete_hanging_nodes(ncont_interpolated_values);
4673
4674
4676 {
4677 t_end = TimingHelpers::timer();
4678 oomph_info << "Time for complete_hanging_nodes: " << t_end - t_start
4679 << std::endl;
4680 t_start = TimingHelpers::timer();
4681 }
4682
4683 /// Update the boundary element info -- this can be a costly procedure
4684 /// and for this reason the mesh writer might have decided not to set up
4685 /// this scheme. If so, we won't change this and suppress its creation...
4687 {
4689 }
4690
4692 {
4693 t_end = TimingHelpers::timer();
4694 oomph_info << "Time for boundary element info: " << t_end - t_start
4695 << std::endl;
4696 t_start = TimingHelpers::timer();
4697 }
4698
4699 // BENFLAG: Reset all the node update elements.
4700 // This is necessary to prevent the following case: A node N is
4701 // shared between two elements, A and B. The update element for
4702 // the node is set to A, say. Element A is p-refined and now
4703 // nolonger has N as a node. However the node update element for N
4704 // is still A but the node doesn't exist in A.
4705 MacroElementNodeUpdateElementBase* first_macro_el_pt =
4706 dynamic_cast<MacroElementNodeUpdateElementBase*>(this->element_pt(0));
4707 if (first_macro_el_pt != 0)
4708 {
4709 // Now set the node update info elementwise
4710 for (unsigned e = 0; e < this->nelement(); e++)
4711 {
4712 // Cast to macro element
4714 dynamic_cast<MacroElementNodeUpdateElementBase*>(
4715 this->element_pt(e));
4716 if (macro_el_pt != 0)
4717 {
4718 // Get vector of geometric objects from element (construct vector
4719 // via copy operation)
4720 Vector<GeomObject*> geom_object_pt(macro_el_pt->geom_object_pt());
4721
4722 // (Re)set node update info for all the nodes in the element
4723 macro_el_pt->set_node_update_info(geom_object_pt);
4724 }
4725 }
4726 }
4727
4728#ifdef PARANOID
4729
4730 // Doc/check the neighbours
4731 //-------------------------
4732 Vector<Tree*> all_tree_nodes_pt;
4733 this->forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
4734
4735 // Check the neighbours
4736 this->forest_pt()->check_all_neighbours(doc_info);
4737
4738 // Check the integrity of the elements
4739 // -----------------------------------
4740
4741 // Loop over elements and get the elemental integrity
4742 double max_error = 0.0;
4743 for (unsigned long e = 0; e < num_tree_nodes; e++)
4744 {
4745 double max_el_error;
4746 tree_nodes_pt[e]->object_pt()->check_integrity(max_el_error);
4747 // If the elemental error is greater than our maximum error
4748 // reset the maximum
4749 if (max_el_error > max_error)
4750 {
4751 max_error = max_el_error;
4752 }
4753 }
4754
4756 {
4757 std::ostringstream error_stream;
4758 error_stream << "Mesh refined: Max. error in integrity check: "
4759 << max_error << " is too big"
4760 << "\ni.e. bigger than RefineableElement::"
4761 << "max_integrity_tolerance()="
4763 << std::endl;
4764
4765 std::ofstream some_file;
4766 some_file.open("ProblemMesh.dat");
4767 for (unsigned long n = 0; n < n_node; n++)
4768 {
4769 // Get the pointer to the node
4770 Node* nod_pt = this->node_pt(n);
4771 // Get the dimension
4772 unsigned n_dim = nod_pt->ndim();
4773 // Output the coordinates
4774 for (unsigned i = 0; i < n_dim; i++)
4775 {
4776 some_file << this->node_pt(n)->x(i) << " ";
4777 }
4778 some_file << std::endl;
4779 }
4780 some_file.close();
4781
4782 error_stream << "Documented problem mesh in ProblemMesh.dat"
4783 << std::endl;
4784
4785 throw OomphLibError(
4786 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4787 }
4788 else
4789 {
4790 oomph_info << "Mesh refined: Max. error in integrity check: "
4791 << max_error << " is OK" << std::endl;
4793 << "i.e. less than RefineableElement::max_integrity_tolerance()="
4795 << std::endl;
4796 }
4797
4798
4800 {
4801 t_end = TimingHelpers::timer();
4802 oomph_info << "Time for (paranoid only) checking of integrity: "
4803 << t_end - t_start << std::endl;
4804 t_start = TimingHelpers::timer();
4805 }
4806
4807#endif
4808
4809 // Loop over all elements other than the final level and deactivate the
4810 // objects, essentially set the pointer that point to nodes that are
4811 // about to be deleted to NULL. This must take place here because nodes
4812 // addressed by elements that are dead but still living in the tree might
4813 // have been made obsolete in the last round of refinement
4814 //(Not strictly required, as tree structure has not changed, but does no
4815 // harm)
4816 for (unsigned long e = 0; e < this->forest_pt()->ntree(); e++)
4817 {
4819 }
4820
4821 // Now we can prune the dead nodes from the mesh.
4822 Vector<Node*> deleted_node_pt = this->prune_dead_nodes();
4823
4825 {
4826 t_end = TimingHelpers::timer();
4827 oomph_info << "Time for deactivating objects and pruning nodes: "
4828 << t_end - t_start << std::endl;
4829 t_start = TimingHelpers::timer();
4830 }
4831
4832 // Finally: Reorder the nodes within the mesh's node vector
4833 // to establish a standard ordering regardless of the sequence
4834 // of mesh refinements -- this is required to allow dump/restart
4835 // on refined meshes
4836 this->reorder_nodes();
4837
4839 {
4840 t_end = TimingHelpers::timer();
4841 oomph_info << "Time for reordering " << nnode()
4842 << " nodes: " << t_end - t_start << std::endl;
4843 t_start = TimingHelpers::timer();
4844 }
4845
4846 // Now we can correct the nodes on boundaries that were hanging that
4847 // are no longer hanging
4848 // Only bother if we have more than two dimensions
4849 if (mesh_dim > 2)
4850 {
4851 // Loop over the boundaries
4852 for (unsigned b = 0; b < n_boundary; b++)
4853 {
4854 // Remove deleted nodes from the set
4855 unsigned n_del = deleted_node_pt.size();
4856 for (unsigned j = 0; j < n_del; j++)
4857 {
4858 hanging_nodes_on_boundary_pt[b].erase(deleted_node_pt[j]);
4859 }
4860
4861 // If the nodes that were hanging are still hanging then remove them
4862 // from the set (note increment is not in for command for efficiencty)
4863 for (std::set<Node*>::iterator it =
4864 hanging_nodes_on_boundary_pt[b].begin();
4865 it != hanging_nodes_on_boundary_pt[b].end();)
4866 {
4867 if ((*it)->is_hanging())
4868 {
4869 hanging_nodes_on_boundary_pt[b].erase(it++);
4870 }
4871 else
4872 {
4873 ++it;
4874 }
4875 }
4876
4877 // Are there any nodes that have changed geometric hanging status
4878 // on the boundary
4879 // The slightly painful part is that we must adjust the position
4880 // via the macro-elements which are only available through the
4881 // elements and not the nodes.
4882 if (hanging_nodes_on_boundary_pt[b].size() > 0)
4883 {
4884 // If so we loop over all elements adjacent to the boundary
4885 unsigned n_boundary_element = this->nboundary_element(b);
4886 for (unsigned e = 0; e < n_boundary_element; ++e)
4887 {
4888 // Get a pointer to the element
4889 FiniteElement* el_pt = this->boundary_element_pt(b, e);
4890
4891 // Do we have a solid element
4892 SolidFiniteElement* solid_el_pt =
4893 dynamic_cast<SolidFiniteElement*>(el_pt);
4894
4895 // Determine whether there is a macro element
4896 bool macro_present = (el_pt->macro_elem_pt() != 0);
4897 // Or a solid macro element
4898 if (solid_el_pt != 0)
4899 {
4900 macro_present |= (solid_el_pt->undeformed_macro_elem_pt() != 0);
4901 }
4902
4903 // Only bother to do anything if there is a macro element
4904 // or undeformed macro element in a SolidElement
4905 if (macro_present)
4906 {
4907 // Loop over the nodes
4908 // ALH: (could optimise to only loop over
4909 // node associated with the boundary with more effort)
4910 unsigned n_el_node = el_pt->nnode();
4911 for (unsigned n = 0; n < n_el_node; n++)
4912 {
4913 // Cache pointer to the node
4914 Node* nod_pt = el_pt->node_pt(n);
4915 if (nod_pt->is_on_boundary(b))
4916 {
4917 // Is the Node in our set
4918 std::set<Node*>::iterator it =
4919 hanging_nodes_on_boundary_pt[b].find(nod_pt);
4920 // If we have found the Node then update the position
4921 // to be consistent with the macro-element representation
4922 if (it != hanging_nodes_on_boundary_pt[b].end())
4923 {
4924 // Specialise local and global coordinates to 3D
4925 // because there is only a problem in 3D.
4926 Vector<double> s(3), x(3);
4927 // Find the local coordinate of the ndoe
4928 el_pt->local_coordinate_of_node(n, s);
4929 // Find the number of time history values
4930 const unsigned ntstorage = nod_pt->ntstorage();
4931
4932 // Do we have a solid node
4933 SolidNode* solid_node_pt =
4934 dynamic_cast<SolidNode*>(nod_pt);
4935 if (solid_node_pt)
4936 {
4937 // Assign Lagrangian coordinates from undeformed
4938 // macro element (if it has one -- get_x_and_xi()
4939 // does "the right thing" anyway. Leave actual
4940 // nodal positions alone -- we're doing a solid
4941 // mechanics problem and once we're going
4942 // the nodal positions are always computed, never
4943 // (automatically) reset to macro-element based
4944 // positions; not even on pinned boundaries
4945 // because the user may have other ideas about where
4946 // these should go -- pinning means "don't touch the
4947 // value", not "leave where the macro-element thinks
4948 // it should be"
4949 Vector<double> x_fe(3), xi(3), xi_fe(3);
4950 solid_el_pt->get_x_and_xi(s, x_fe, x, xi_fe, xi);
4951 for (unsigned i = 0; i < 3; i++)
4952 {
4953 solid_node_pt->xi(i) = xi[i];
4954 }
4955 }
4956 else
4957 {
4958 // Set position and history values from the
4959 // macro-element representation
4960 for (unsigned t = 0; t < ntstorage; t++)
4961 {
4962 // Get the history value from the macro element
4963 el_pt->get_x(t, s, x);
4964
4965 // Set the coordinate to that of the macroelement
4966 // representation
4967 for (unsigned i = 0; i < 3; i++)
4968 {
4969 nod_pt->x(t, i) = x[i];
4970 }
4971 }
4972 } // End of non-solid node case
4973
4974 // Now remove the node from the list
4975 hanging_nodes_on_boundary_pt[b].erase(it);
4976 // If there are no Nodes left then exit the loops
4977 if (hanging_nodes_on_boundary_pt[b].size() == 0)
4978 {
4979 e = n_boundary_element;
4980 break;
4981 }
4982 }
4983 }
4984 }
4985 } // End of macro element case
4986 }
4987 }
4988 }
4989 } // End of case when we have fixed nodal positions
4990
4991 // Final doc
4992 //-----------
4994 {
4995 // Doc the boundary conditions ('0' for non-existent, '1' for free,
4996 //----------------------------------------------------------------
4997 // '2' for pinned -- ideal for tecplot scatter sizing.
4998 //----------------------------------------------------
4999 // num_tree_nodes=tree_nodes_pt.size();
5000
5001 // Determine maximum number of values at any node in this type of
5002 // element
5003 RefineableElement* el_pt = tree_nodes_pt[0]->object_pt();
5004 // Initalise max_nval
5005 unsigned max_nval = 0;
5006 for (unsigned n = 0; n < el_pt->nnode(); n++)
5007 {
5008 if (el_pt->node_pt(n)->nvalue() > max_nval)
5009 {
5010 max_nval = el_pt->node_pt(n)->nvalue();
5011 }
5012 }
5013
5014 // Open the output file
5015 std::ostringstream fullname;
5016 std::ofstream bcs_file;
5017 fullname.str("");
5018 fullname << doc_info.directory() << "/bcs" << doc_info.number()
5019 << ".dat";
5020 bcs_file.open(fullname.str().c_str());
5021
5022 // Loop over elements
5023 for (unsigned long e = 0; e < num_tree_nodes; e++)
5024 {
5025 el_pt = tree_nodes_pt[e]->object_pt();
5026 // Loop over nodes in element
5027 unsigned n_nod = el_pt->nnode();
5028 for (unsigned n = 0; n < n_nod; n++)
5029 {
5030 // Get pointer to the node
5031 Node* nod_pt = el_pt->node_pt(n);
5032 // Find the dimension of the node
5033 unsigned n_dim = nod_pt->ndim();
5034 // Write the nodal coordinates to the file
5035 for (unsigned i = 0; i < n_dim; i++)
5036 {
5037 bcs_file << nod_pt->x(i) << " ";
5038 }
5039
5040 // Loop over all values in this element
5041 for (unsigned i = 0; i < max_nval; i++)
5042 {
5043 // Value exists at this node:
5044 if (i < nod_pt->nvalue())
5045 {
5046 bcs_file << " " << 1 + nod_pt->is_pinned(i);
5047 }
5048 // ...if not just dump out a zero
5049 else
5050 {
5051 bcs_file << " 0 ";
5052 }
5053 }
5054 bcs_file << std::endl;
5055 }
5056 }
5057 bcs_file.close();
5058
5059 // Doc all nodes
5060 //---------------
5061 std::ofstream all_nodes_file;
5062 fullname.str("");
5063 fullname << doc_info.directory() << "/all_nodes" << doc_info.number()
5064 << ".dat";
5065 all_nodes_file.open(fullname.str().c_str());
5066
5067 all_nodes_file << "ZONE \n";
5068
5069 // Need to recompute the number of nodes since it may have
5070 // changed during mesh refinement/unrefinement
5071 n_node = this->nnode();
5072 for (unsigned long n = 0; n < n_node; n++)
5073 {
5074 Node* nod_pt = this->node_pt(n);
5075 unsigned n_dim = nod_pt->ndim();
5076 for (unsigned i = 0; i < n_dim; i++)
5077 {
5078 all_nodes_file << this->node_pt(n)->x(i) << " ";
5079 }
5080 all_nodes_file << std::endl;
5081 }
5082
5083 all_nodes_file.close();
5084
5085
5086 // Doc all hanging nodes:
5087 //-----------------------
5088 std::ofstream some_file;
5089 fullname.str("");
5090 fullname << doc_info.directory() << "/all_hangnodes"
5091 << doc_info.number() << ".dat";
5092 some_file.open(fullname.str().c_str());
5093 for (unsigned long n = 0; n < n_node; n++)
5094 {
5095 Node* nod_pt = this->node_pt(n);
5096
5097 if (nod_pt->is_hanging())
5098 {
5099 unsigned n_dim = nod_pt->ndim();
5100 for (unsigned i = 0; i < n_dim; i++)
5101 {
5102 some_file << nod_pt->x(i) << " ";
5103 }
5104
5105 // ALH: Added this to stop Solid problems seg-faulting
5106 if (this->node_pt(n)->nvalue() > 0)
5107 {
5108 some_file << " " << nod_pt->raw_value(0);
5109 }
5110 some_file << std::endl;
5111 }
5112 }
5113 some_file.close();
5114
5115 // Doc all hanging nodes and their masters
5116 // View with QHangingNodesWithMasters.mcr
5117 fullname.str("");
5118 fullname << doc_info.directory() << "/geometric_hangnodes_withmasters"
5119 << doc_info.number() << ".dat";
5120 some_file.open(fullname.str().c_str());
5121 for (unsigned long n = 0; n < n_node; n++)
5122 {
5123 Node* nod_pt = this->node_pt(n);
5124 if (nod_pt->is_hanging())
5125 {
5126 unsigned n_dim = nod_pt->ndim();
5127 unsigned nmaster = nod_pt->hanging_pt()->nmaster();
5128 some_file << "ZONE I=" << nmaster + 1 << std::endl;
5129 for (unsigned i = 0; i < n_dim; i++)
5130 {
5131 some_file << nod_pt->x(i) << " ";
5132 }
5133 some_file << " 2 " << std::endl;
5134
5135 for (unsigned imaster = 0; imaster < nmaster; imaster++)
5136 {
5137 Node* master_nod_pt =
5138 nod_pt->hanging_pt()->master_node_pt(imaster);
5139 unsigned n_dim = master_nod_pt->ndim();
5140 for (unsigned i = 0; i < n_dim; i++)
5141 {
5142 some_file << master_nod_pt->x(i) << " ";
5143 }
5144 some_file << " 1 " << std::endl;
5145 }
5146 }
5147 }
5148 some_file.close();
5149
5150 // Doc all hanging nodes and their masters
5151 // View with QHangingNodesWithMasters.mcr
5152 for (unsigned i = 0; i < ncont_interpolated_values; i++)
5153 {
5154 fullname.str("");
5155 fullname << doc_info.directory()
5156 << "/nonstandard_hangnodes_withmasters" << i << "_"
5157 << doc_info.number() << ".dat";
5158 some_file.open(fullname.str().c_str());
5159 unsigned n_nod = this->nnode();
5160 for (unsigned long n = 0; n < n_nod; n++)
5161 {
5162 Node* nod_pt = this->node_pt(n);
5163 if (nod_pt->is_hanging(i))
5164 {
5165 if (nod_pt->hanging_pt(i) != nod_pt->hanging_pt())
5166 {
5167 unsigned nmaster = nod_pt->hanging_pt(i)->nmaster();
5168 some_file << "ZONE I=" << nmaster + 1 << std::endl;
5169 unsigned n_dim = nod_pt->ndim();
5170 for (unsigned j = 0; j < n_dim; j++)
5171 {
5172 some_file << nod_pt->x(j) << " ";
5173 }
5174 some_file << " 2 " << std::endl;
5175 for (unsigned imaster = 0; imaster < nmaster; imaster++)
5176 {
5177 Node* master_nod_pt =
5178 nod_pt->hanging_pt(i)->master_node_pt(imaster);
5179 unsigned n_dim = master_nod_pt->ndim();
5180 for (unsigned j = 0; j < n_dim; j++)
5181 {
5182 // some_file << master_nod_pt->x(i) << " ";
5183 }
5184 some_file << " 1 " << std::endl;
5185 }
5186 }
5187 }
5188 }
5189 some_file.close();
5190 }
5191
5192 } // End of documentation
5193 } // End if (this->nelement()>0)
5194
5195 /// /BENFLAG: Check that all the nodes belong to their update elements
5196 // std::cout << "p_adapt_mesh(): Checking stuff works!" << std::endl;
5197 // for(unsigned j=0; j<this->nnode(); j++)
5198 // {
5199 // MacroElementNodeUpdateNode* macro_nod_pt =
5200 // dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j));
5201 // if(macro_nod_pt!=0)
5202 // {
5203 // bool big_problem = true;
5204 // std::cout << "Node " << macro_nod_pt << " at [ " << macro_nod_pt->x(0)
5205 // << ", " << macro_nod_pt->x(1) << " ]" << std::endl; FiniteElement*
5206 // up_el_pt =
5207 // dynamic_cast<FiniteElement*>(macro_nod_pt->node_update_element_pt());
5208 // for(unsigned l=0; l<up_el_pt->nnode(); l++)
5209 // {
5210 // if(up_el_pt->node_pt(l)==macro_nod_pt)
5211 // {
5212 // big_problem = false;
5213 // break;
5214 // }
5215 // }
5216 // if(big_problem)
5217 // {
5218 // std::cout << " This node doesn't exist in it's update element!" <<
5219 // std::endl;
5220 // }
5221 // }
5222 // }
5223
5224#ifdef OOMPH_HAS_MPI
5225
5226 // Now (re-)classify halo and haloed nodes and synchronise hanging
5227 // nodes
5228 if (this->is_mesh_distributed())
5229 {
5231 }
5232
5233#endif
5234 }
5235
5236 //========================================================================
5237 /// p-unrefine mesh uniformly
5238 /// Unlike in h-refinement, we can simply p-unrefine each element in the mesh
5239 //========================================================================
5241 {
5242 // Select all elements for unrefinement
5243 unsigned long Nelement = this->nelement();
5244 for (unsigned long e = 0; e < Nelement; e++)
5245 {
5246 // Get pointer to p-refineable element
5247 PRefineableElement* el_pt =
5248 dynamic_cast<PRefineableElement*>(this->element_pt(e));
5249 // Mark for p-refinement if possible. If not then p_adapt_mesh() will
5250 // report the error.
5251 if (el_pt != 0)
5252 {
5254 }
5255 }
5256
5257 // Do the actual mesh adaptation
5259 }
5260
5261 //========================================================================
5262 /// p-refine mesh by refining the elements identified by their numbers.
5263 //========================================================================
5265 const Vector<unsigned>& elements_to_be_refined)
5266 {
5267#ifdef OOMPH_HAS_MPI
5268 if (this->is_mesh_distributed())
5269 {
5270 std::ostringstream warn_stream;
5271 warn_stream << "You are attempting to refine selected elements of a "
5272 << std::endl
5273 << "distributed mesh. This may have undesired effects."
5274 << std::endl;
5275
5276 OomphLibWarning(warn_stream.str(),
5277 "TreeBasedRefineableMeshBase::refine_selected_elements()",
5278 OOMPH_EXCEPTION_LOCATION);
5279 }
5280#endif
5281
5282 // Select elements for refinement
5283 unsigned long nref = elements_to_be_refined.size();
5284 for (unsigned long e = 0; e < nref; e++)
5285 {
5286 // Get pointer to p-refineable element
5287 PRefineableElement* el_pt = dynamic_cast<PRefineableElement*>(
5288 this->element_pt(elements_to_be_refined[e]));
5289 // Mark for p-refinement if possible. If not then p_adapt_mesh() will
5290 // report the error.
5291 if (el_pt != 0)
5292 {
5293 el_pt->select_for_p_refinement();
5294 }
5295 }
5296
5297 // Do the actual mesh adaptation
5298 p_adapt_mesh();
5299 }
5300
5301 //========================================================================
5302 /// p-refine mesh by refining the elements identified by their pointers.
5303 //========================================================================
5305 const Vector<PRefineableElement*>& elements_to_be_refined_pt)
5306 {
5307#ifdef OOMPH_HAS_MPI
5308 if (this->is_mesh_distributed())
5309 {
5310 std::ostringstream warn_stream;
5311 warn_stream << "You are attempting to refine selected elements of a "
5312 << std::endl
5313 << "distributed mesh. This may have undesired effects."
5314 << std::endl;
5315
5316 OomphLibWarning(warn_stream.str(),
5317 "TreeBasedRefineableMeshBase::refine_selected_elements()",
5318 OOMPH_EXCEPTION_LOCATION);
5319 }
5320#endif
5321
5322 // Select elements for refinement
5323 unsigned long nref = elements_to_be_refined_pt.size();
5324 for (unsigned long e = 0; e < nref; e++)
5325 {
5326 elements_to_be_refined_pt[e]->select_for_p_refinement();
5327 }
5328
5329 // Do the actual mesh adaptation
5330 p_adapt_mesh();
5331 }
5332
5333} // namespace oomph
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
////////////////////////////////////////////////////////////////////
void node_update(const bool &update_all_time_levels_for_new_node=false)
Broken assignment operator.
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:879
int non_halo_proc_ID()
ID of processor ID that holds non-halo counterpart of halo node; negative if not a halo.
Definition: nodes.h:539
bool is_pinned(const unsigned &i) const
Test whether the i-th variable is pinned (1: true; 0: false).
Definition: nodes.h:417
Information for documentation of results: Directory and file number to enable output in the form RESL...
bool is_doc_enabled() const
Are we documenting?
void disable_doc()
Disable documentation.
std::string directory() const
Output directory.
unsigned & number()
Number used (e.g.) for labeling output files.
A general Finite Element class.
Definition: elements.h:1313
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
Definition: elements.h:1842
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void get_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
Definition: elements.h:1885
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Definition: elements.h:1878
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
void set_master_node_pt(const unsigned &i, Node *const &master_node_pt, const double &weight)
Set the pointer to the i-th master node and its weight.
Definition: nodes.cc:1474
////////////////////////////////////////////////////////////////////
virtual void set_node_update_info(const Vector< GeomObject * > &geom_object_pt)=0
Set node update information: Pass the vector of (pointers to) the geometric objects that affect the n...
Vector< GeomObject * > & geom_object_pt()
Vector of (pointers to) geometric objects involved in node update function.
A general mesh class.
Definition: mesh.h:67
std::map< unsigned, Vector< Node * > > Shared_node_pt
Map of vectors holding the pointers to the shared nodes. These are all the nodes that are on two "nei...
Definition: mesh.h:116
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
bool Lookup_for_elements_next_boundary_is_setup
Flag to indicate that the lookup schemes for elements that are adjacent to the boundaries has been se...
Definition: mesh.h:87
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
virtual void setup_boundary_element_info()
Interface for function that is used to setup the boundary information (Empty virtual function – imple...
Definition: mesh.h:275
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
OomphCommunicator * Comm_pt
Pointer to communicator – set to NULL if mesh is not distributed.
Definition: mesh.h:119
virtual void reorder_nodes(const bool &use_old_ordering=true)
Re-order nodes in the order in which they appear in elements – can be overloaded for more efficient r...
Definition: mesh.cc:508
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
Vector< GeneralisedElement * > halo_element_pt(const unsigned &p)
Return vector of halo elements in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:1740
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
void delete_all_external_storage()
Wipe the storage for all externally-based elements.
Definition: mesh.cc:9190
unsigned nhalo_node()
Total number of halo nodes in this Mesh.
Definition: mesh.h:1878
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
Node * haloed_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th haloed node in this Mesh whose halo counterpart is held on processor p.
Definition: mesh.h:2014
Vector< Node * > prune_dead_nodes()
Prune nodes. Nodes that have been marked as obsolete are removed from the mesh (and its boundary-node...
Definition: mesh.cc:966
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
unsigned nhaloed_node()
Total number of haloed nodes in this Mesh.
Definition: mesh.h:1985
Node * halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th halo node in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:1914
Vector< GeneralisedElement * > haloed_element_pt(const unsigned &p)
Return vector of haloed elements in this Mesh whose haloing counterpart is held on processor p.
Definition: mesh.h:1779
Node * shared_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th shared node in this Mesh who has a counterpart on processor p.
Definition: mesh.h:2110
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
void set_nonhanging()
Label node as non-hanging node by removing all hanging node data.
Definition: nodes.cc:2281
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition: nodes.h:1373
void set_obsolete()
Mark node as obsolete.
Definition: nodes.h:1436
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1365
void position(Vector< double > &pos) const
Compute Vector of nodal positions either directly or via hanging node representation.
Definition: nodes.cc:2499
void set_non_obsolete()
Mark node as non-obsolete.
Definition: nodes.h:1442
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
void set_hanging_pt(HangInfo *const &hang_pt, const int &i)
Set the hanging data for the i-th value. (hang_pt=0 to make non-hanging)
Definition: nodes.cc:2068
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
double raw_value(const unsigned &i) const
Return the i-th value stored at the Node. This interface does NOT take the hanging status of the Node...
Definition: nodes.h:1455
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
Definition: nodes.cc:2408
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
bool p_refinement_is_enabled()
Flag to indicate suppression of any refinement.
virtual unsigned initial_p_order() const =0
Get the initial P_order This is required so that elements which are constructed with a higher p-order...
unsigned & p_order()
Access function to P_order.
void select_for_p_refinement()
Select the element for p-refinement.
void deselect_for_p_unrefinement()
Deselect the element for p-unrefinement.
void select_for_p_unrefinement()
Select the element for p-unrefinement.
void deselect_for_p_refinement()
Deselect the element for p-refinement.
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
void select_sons_for_unrefinement()
Unrefinement will be performed by merging the four sons of this element.
void select_for_refinement()
Select the element for refinement.
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
static double & max_integrity_tolerance()
Max. allowed discrepancy in element integrity check.
unsigned refinement_level() const
Return the Refinement level.
bool refinement_is_enabled()
Flag to indicate suppression of any refinement.
void deselect_for_refinement()
Deselect the element for refinement.
void deselect_sons_for_unrefinement()
No unrefinement will be performed by merging the four sons of this element.
unsigned & max_keep_unrefined()
Max. number of elements that we allow to remain unrefined if no other mesh adaptation is required (to...
unsigned Nrefined
Stats: Number of elements that were refined.
double & min_permitted_error()
Access fct for min. error (i.e. (try to) merge elements if their error is smaller)
DocInfo *& doc_info_pt()
Access fct for pointer to DocInfo.
DocInfo doc_info()
Access fct for DocInfo.
unsigned & nrefinement_overruled()
Number of elements that would have liked to be refined further but can't because they've reached the ...
unsigned Nunrefined
Stats: Number of elements that were unrefined.
double & max_error()
Access fct for max. actual error in present solution (i.e. before re-solve on adapted mesh)
double & max_permitted_error()
Access fct for max. error (i.e. split elements if their error is larger)
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: elements.h:3561
MacroElement * undeformed_macro_elem_pt()
Access function to pointer to "undeformed" macro element.
Definition: elements.h:3706
virtual void get_x_and_xi(const Vector< double > &s, Vector< double > &x_fe, Vector< double > &x, Vector< double > &xi_fe, Vector< double > &xi) const
Eulerian and Lagrangian coordinates as function of the local coordinates: The Eulerian position is re...
Definition: elements.h:3659
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
double lagrangian_position(const unsigned &i) const
Return lagrangian coordinate either directly or via hanging node representation.
Definition: nodes.cc:3586
double & xi(const unsigned &i)
Reference to i-th Lagrangian position.
Definition: nodes.h:1883
unsigned nlagrangian() const
Return number of lagrangian coordinates.
Definition: nodes.h:1870
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
TreeForest * forest_pt()
Return pointer to the Forest represenation of the mesh.
void synchronise_hanging_nodes(const unsigned &ncont_interpolated_values)
Synchronise the hanging nodes if the mesh is distributed.
virtual void read_refinement(std::ifstream &restart_file, Vector< Vector< unsigned > > &to_be_refined)
Read refinement pattern to allow for rebuild.
virtual void get_elements_at_refinement_level(unsigned &refinement_level, Vector< RefineableElement * > &level_elements)
Extract the elements at a particular refinement level in the refinement pattern - used in Mesh::redis...
virtual void refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
Refine mesh by splitting the elements identified by their numbers.
unsigned & max_p_refinement_level()
Access fct for max. permissible p-refinement level (relative to base mesh)
virtual void dump_refinement(std::ostream &outfile)
Dump refinement pattern to allow for rebuild.
unsigned unrefine_uniformly()
Unrefine mesh uniformly: Return 0 for success, 1 for failure (if unrefinement has reached the coarses...
virtual void refine(std::ifstream &restart_file)
Refine mesh according to refinement pattern in restart file.
virtual void additional_synchronise_hanging_nodes(const unsigned &ncont_interpolated_values)=0
Additional synchronisation of hanging nodes Required for reconcilliation of hanging nodes on the oute...
unsigned & min_refinement_level()
Access fct for min. permissible refinement level (relative to base mesh)
virtual void adapt_mesh()
Perform the actual tree-based mesh adaptation. A simple wrapper to call the function without document...
virtual void refine_base_mesh_as_in_reference_mesh(TreeBasedRefineableMeshBase *const &ref_mesh_pt)
Refine base mesh to same degree as reference mesh (relative to original unrefined mesh).
virtual void get_refinement_pattern(Vector< Vector< unsigned > > &to_be_refined)
Extract refinement pattern: Consider the hypothetical mesh obtained by truncating the refinement of t...
void complete_hanging_nodes_recursively(Node *&nod_pt, Vector< Node * > &master_nodes, Vector< double > &hang_weights, const int &ival)
Auxiliary routine for recursive hanging node completion.
unsigned & max_refinement_level()
Access fct for max. permissible refinement level (relative to base mesh)
virtual void refine_as_in_reference_mesh(TreeBasedRefineableMeshBase *const &ref_mesh_pt)
Refine mesh once so that its topology etc becomes that of the (finer!) reference mesh – if possible!...
void refine_base_mesh(Vector< Vector< unsigned > > &to_be_refined)
Refine base mesh according to specified refinement pattern.
void p_refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
p-refine mesh by refining the elements identified by their numbers.
virtual void split_elements_if_required()=0
Split all the elements in the mesh if required. This template free interface will be overloaded in Re...
void refine_uniformly()
Refine mesh uniformly.
void p_unrefine_uniformly(DocInfo &doc_info)
p-unrefine mesh uniformly
virtual void get_refinement_levels(unsigned &min_refinement_level, unsigned &max_refinement_level)
Get max/min refinement levels in mesh.
virtual void p_refine_elements_if_required()=0
p-refine all the elements in the mesh if required. This template free interface will be overloaded in...
void complete_hanging_nodes(const int &ncont_interpolated_values)
Complete the hanging node scheme recursively.
TreeForest * Forest_pt
Forest representation of the mesh.
void synchronise_nonhanging_nodes()
Synchronise the positions of non-hanging nodes that depend on non-existent neighbours (e....
void adapt(const Vector< double > &elemental_error)
Adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose error...
void p_adapt(const Vector< double > &elemental_error)
p-adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose err...
void p_adapt_mesh()
Perform the actual tree-based mesh p-adaptation. A simple wrapper to call the function without docume...
virtual bool refine_base_mesh_as_in_reference_mesh_minus_one(TreeBasedRefineableMeshBase *const &ref_mesh_pt)
Refine base mesh to same degree as reference mesh minus one level of refinement (relative to original...
void p_refine_uniformly()
p-refine mesh uniformly
void classify_halo_and_haloed_nodes(DocInfo &doc_info, const bool &report_stats)
Classify all halo and haloed information in the mesh (overloaded version from Mesh base class....
virtual void check_all_neighbours(DocInfo &doc_info)=0
Document/check the neighbours of all the nodes in the forest. This must be overloaded for different t...
void stick_all_tree_nodes_into_vector(Vector< Tree * > &all_forest_nodes)
Traverse forest and stick pointers to all "nodes" into Vector.
Definition: tree.cc:405
virtual void open_hanging_node_files(DocInfo &doc_info, Vector< std::ofstream * > &output_stream)=0
Open output files that will store any hanging nodes in the forest. Return a vector of the output stre...
unsigned ntree()
Number of trees in forest.
Definition: tree.h:460
void stick_leaves_into_vector(Vector< Tree * > &forest_nodes)
Traverse forst and stick pointers to leaf "nodes" into Vector.
Definition: tree.cc:379
TreeRoot * tree_pt(const unsigned &i) const
Return pointer to i-th tree in forest.
Definition: tree.h:466
void close_hanging_node_files(DocInfo &doc_info, Vector< std::ofstream * > &output_stream)
Close output files that will store any hanging nodes in the forest and delete any associated storage....
Definition: tree.cc:432
A generalised tree base class that abstracts the common functionality between the quad- and octrees u...
Definition: tree.h:74
unsigned nsons() const
Return number of sons (zero if it's a leaf node)
Definition: tree.h:129
void traverse_all(Tree::VoidMemberFctPt member_function)
Traverse the tree and execute void Tree member function member_function() at all its "nodes".
Definition: tree.cc:145
bool is_leaf()
Return true if the tree is a leaf node.
Definition: tree.h:220
int son_type() const
Return son type.
Definition: tree.h:214
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:88
Tree * son_pt(const int &son_index) const
Return pointer to the son for a given index. Note that to aid code readability specific enums have be...
Definition: tree.h:103
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:235
void deactivate_object()
Call the RefineableElement's deactivate_element() function.
Definition: tree.cc:341
void traverse_all_but_leaves(Tree::VoidMemberFctPt member_function)
Traverse the tree and execute void Tree member function member_function() at all its "nodes" aparat f...
Definition: tree.cc:184
void merge_sons_if_required(Mesh *&mesh_pt)
If required, merge the four sons for unrefinement – criterion: bool object_pt()-> sons_to_be_unrefine...
Definition: tree.cc:302
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
bool Doc_comprehensive_timings
Global boolean to switch on comprehensive timing – can probably be declared const false when developm...
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void pause(std::string message)
Pause and display message.
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...
Helper struct to collate data required during TreeBasedRefineableMeshBase::synchronise_hanging_nodes.