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