refineable_mesh.cc

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================


#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif

#include <cstdlib>
#include <stdlib.h>
#include <limits>
 
#include "refineable_mesh.h"
//Include to fill in additional_synchronise_hanging_nodes() function
#include "refineable_mesh.template.cc"

namespace oomph
{

//========================================================================
/// Get refinement pattern of mesh: Consider the hypothetical mesh
/// obtained by truncating the refinement of the current mesh to a given level 
/// (where \c level=0 is the un-refined base mesh). To advance
/// to the next refinement level, we need to refine (split) the
/// \c to_be_refined[level].size() elements identified by the
/// element numbers contained in \c vector to_be_refined[level][...]
//========================================================================
void TreeBasedRefineableMeshBase::get_refinement_pattern(
 Vector<Vector<unsigned> >& to_be_refined)
{
 // Extract *all* elements from current (fully refined) mesh.
 Vector<Tree*> all_tree_nodes_pt;
 forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
 
 // Find out maximum refinement level
 unsigned max_level=0;
 unsigned nnodes=all_tree_nodes_pt.size();
 for (unsigned e=0;e<nnodes;e++)
  {
   unsigned level=all_tree_nodes_pt[e]->level();
   if (level>max_level) max_level=level;
  }

 // Assign storage for refinement pattern
 to_be_refined.clear();
 to_be_refined.resize(max_level);
 Vector<unsigned> el_count(max_level);
 
 // Initialise count of elements that exist in mesh when refinement
 // has proceeded to this level
 for (unsigned l=0;l<max_level;l++)
  {
   el_count[l]=0;
  }
   
 // Loop over all levels and extract all elements that exist
 // in reference mesh when refinement has proceeded to this level
 for (unsigned l=0;l<max_level;l++)
  {
   // Loop over all elements (tree nodes)
   for (unsigned e=0;e<nnodes;e++)
    {
     // What level does this element exist on?
     unsigned level=all_tree_nodes_pt[e]->level();
   
     // Element is part of the mesh at this refinement level
     // if it exists at this level OR if it exists at a lower level
     // and is a leaf
     if ((level==l)||((level<l)&&(all_tree_nodes_pt[e]->is_leaf())))
      {
       // If element exsts at this level and is not a leaf it will
       // be refined when we move to the next level:
       if ((level==l)&&(!all_tree_nodes_pt[e]->is_leaf()))
        {
         // Add element number (in mesh at current refinement level)
         // to the list of elements that need to be refined
         to_be_refined[l].push_back(el_count[l]);
        }
       // Element exists in this mesh: Add to counter
       el_count[l]++;
      }
    }
  }
}

//========================================================================
/// \short Extract the elements at a particular refinement level in
/// the refinement pattern (used in Mesh::redistribute or whatever it's
/// going to be called (RefineableMeshBase::reduce_halo_layers or something)
//========================================================================
void TreeBasedRefineableMeshBase::get_elements_at_refinement_level(
 unsigned& refinement_level,
 Vector<RefineableElement*>& level_elements)
{
 // Extract *all* elements from current (fully refined) mesh.
 Vector<Tree*> all_tree_nodes_pt;
 forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);

 // Add the element to the vector if its level matches refinement_level
 unsigned nnodes=all_tree_nodes_pt.size();
 for (unsigned e=0;e<nnodes;e++)
  {
   unsigned level=all_tree_nodes_pt[e]->level();
   if (level==refinement_level)
    {
     level_elements.push_back(dynamic_cast<RefineableElement*>
                              (all_tree_nodes_pt[e]->object_pt()));
    }
  }

}

//========================================================================
/// Refine original, unrefined mesh according to specified refinement 
/// pattern (relative to original, unrefined mesh).
//========================================================================
void TreeBasedRefineableMeshBase::refine_base_mesh(Vector<Vector<unsigned> >& 
                                                   to_be_refined)
{
 // Get mesh back to unrefined state
 unsigned my_max,my_min;
 get_refinement_levels(my_min,my_max);

 // Max refinement level:
 unsigned my_max_level=to_be_refined.size();

 unsigned global_max=0;
 unsigned global_max_level=0;
 Vector<unsigned> data(2,0);
 data[0]=my_max;
 data[1]=my_max_level;
 Vector<unsigned> global_data(2,0);
#ifdef OOMPH_HAS_MPI
 if (this->is_mesh_distributed())
  {
   MPI_Allreduce(&data[0],&global_data[0],2,MPI_UNSIGNED,MPI_MAX,
                 Comm_pt->mpi_comm());
   global_max=global_data[0];
   global_max_level=global_data[1];

  }
 else
#endif
  {
   global_max=my_max;
   global_max_level=my_max_level;
  }


 for (unsigned i=0;i<global_max;i++)
  {   
   unrefine_uniformly(); 
  }
 
 // Do refinement steps in current mesh
 for (unsigned l=0;l<global_max_level;l++)
  {
   // Loop over elements that need to be refined at this level
   unsigned n_to_be_refined=0;
   if (l<my_max_level) n_to_be_refined=to_be_refined[l].size();

   // Select relevant elements to be refined
   for (unsigned i=0;i<n_to_be_refined;i++)
    {
     dynamic_cast<RefineableElement*>(
      this->element_pt(to_be_refined[l][i]))->select_for_refinement();
    }
   
   // Now do the actual mesh refinement
   adapt_mesh();
  }

}


//========================================================================
/// Refine base mesh according to refinement pattern in restart file
//========================================================================
void TreeBasedRefineableMeshBase::refine(std::ifstream& restart_file)
{
 // Assign storage for refinement pattern
 Vector<Vector<unsigned> > to_be_refined;

 // Read refinement pattern
 read_refinement(restart_file,to_be_refined);

 // Refine
 refine_base_mesh(to_be_refined);
}


//========================================================================
/// Dump refinement pattern to allow for rebuild
///
//========================================================================
void TreeBasedRefineableMeshBase::dump_refinement(std::ostream &outfile)
{
 // Assign storage for refinement pattern
 Vector<Vector<unsigned> > to_be_refined;

 // Get refinement pattern of reference mesh:
 get_refinement_pattern(to_be_refined);

 // Dump max refinement level:
 unsigned max_level=to_be_refined.size();
 outfile << max_level << " # max. refinement level " << std::endl;

 // Doc the numbers of the elements that need to be refined at this level
 for (unsigned l=0;l<max_level;l++)
  {
   // Loop over elements that need to be refined at this level
   unsigned n_to_be_refined=to_be_refined[l].size();
   outfile << n_to_be_refined << " # number of elements to be refined. " 
           << "What follows are the numbers of the elements. " << std::endl;

   // Select relevant elements to be refined
   for (unsigned i=0;i<n_to_be_refined;i++)
    {
     outfile << to_be_refined[l][i] << std::endl;
    }
  }   
}


//========================================================================
/// Read refinement pattern to allow for rebuild
///
//========================================================================
void TreeBasedRefineableMeshBase::read_refinement(
 std::ifstream& restart_file, Vector<Vector<unsigned> >& to_be_refined)
{

 std::string input_string;

 // Read max refinement level:

 // Read line up to termination sign
 getline(restart_file,input_string,'#');

 // Ignore rest of line
 restart_file.ignore(80,'\n');

 // Convert
 unsigned max_level=std::atoi(input_string.c_str());

 // Assign storage for refinement pattern
 to_be_refined.resize(max_level);

 // Read the number of the elements that need to be refined at different levels
 for (unsigned l=0;l<max_level;l++)
  {

   // Read line up to termination sign
   getline(restart_file,input_string,'#');

   // Ignore rest of line
   restart_file.ignore(80,'\n');

   // Convert
   unsigned n_to_be_refined=atoi(input_string.c_str());;

   // Assign storage
   to_be_refined[l].resize(n_to_be_refined);

   // Read numbers of the elements that need to be refined
   for (unsigned i=0;i<n_to_be_refined;i++)
    {
     restart_file >> to_be_refined[l][i];
    }
  }   
}





//========================================================================
/// Do adaptive refinement for mesh.
/// - Pass Vector of error estimates for all elements.
/// - Refine those whose errors exceeds the threshold
/// - (Try to) unrefine those whose errors is less than
///   threshold (only possible if the three brothers also want to be 
///   unrefined, of course.)
/// - Update the nodal positions in the whole lot
/// - Store # of refined/unrefined elements.
/// - Doc refinement process (if required)
//========================================================================
void TreeBasedRefineableMeshBase::adapt(const Vector<double>& elemental_error)
 {
  //Set the refinement tolerance to be the max permissible error
  double refine_tol=max_permitted_error();

  //Set the unrefinement tolerance to be the min permissible error 
  double unrefine_tol=min_permitted_error();

  // Setup doc info
  DocInfo local_doc_info;
  if (doc_info_pt()==0) {local_doc_info.disable_doc();}
  else {local_doc_info=doc_info();}
                              

  // Check that the errors make sense
  if (refine_tol<=unrefine_tol)
   {
    std::ostringstream error_stream;
    error_stream << "Refinement tolerance <= Unrefinement tolerance" 
                 << refine_tol << " " << unrefine_tol << std::endl
                 << "doesn't make sense and will almost certainly crash" 
                 << std::endl
                 << "this beautiful code!" << std::endl;
    
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }


  // Select elements for refinement and unrefinement
  //================================================
  // Reset counter for number of elements that would like to be
  // refined further but can't
  nrefinement_overruled()=0;

  // Note: Yes, this needs to be a map because we'll have to check
  // the refinement wishes of brothers (who we only access via pointers)
  std::map<RefineableElement*,bool> wants_to_be_unrefined;

  // Initialise a variable to store the number of elements for refinment
  unsigned n_refine=0;
  
  // Loop over all elements and mark them according to the error criterion
  unsigned long Nelement=this->nelement();
  for (unsigned long e=0;e<Nelement;e++)
  {
   // (Cast) pointer to the element
   RefineableElement* el_pt =
    dynamic_cast<RefineableElement*>(this->element_pt(e));
 
   // Initially element is not to be refined
   el_pt->deselect_for_refinement();
    
   // If the element error exceeds the threshold ...
   if(elemental_error[e] > refine_tol)
   {
    // ... and its refinement level is less than the maximum desired level
    // mark is to be refined
    if ((el_pt->refinement_is_enabled())&&
        (el_pt->refinement_level() < max_refinement_level()))
    {
     el_pt->select_for_refinement();
     n_refine++;
    }
    // ... otherwise mark it as having been over-ruled
    else
    {
     nrefinement_overruled()+=1;
    }
   }
    
   // Now worry about unrefinement (first pass):

   // Is my error too small AND do I have a father?
   if ((elemental_error[e]<unrefine_tol)&&
       (el_pt->tree_pt()->father_pt()!=0))
   {
    // Flag to indicate whether to unrefine
    bool unrefine=true;
    unsigned n_sons = el_pt->tree_pt()->father_pt()->nsons();
    
    // Are all brothers leaf nodes?
    for (unsigned ison=0;ison<n_sons;ison++)
    {
     // (At least) one brother is not a leaf: end of story; we're not doing
     // it (= the unrefinement)
     if (!(el_pt->tree_pt()->father_pt()->son_pt(ison)->is_leaf()))
     {unrefine=false;}
    }

    // Don't allow unrefinement of elements that would become larger
    // than the minimum legal refinement level
    if (el_pt->refinement_level()-1<min_refinement_level())
    {unrefine=false;}

    // So, all things considered, is the element eligbible for refinement?
    if(unrefine) {wants_to_be_unrefined[el_pt]=true;}
    else {wants_to_be_unrefined[el_pt]=false;}
   }
  }
  
  oomph_info << " \n Number of elements to be refined: "
             << n_refine << std::endl;
  oomph_info << " \n Number of elements whose refinement was overruled: " 
             << nrefinement_overruled() << std::endl;

  // Second pass for unrefinement --- an element cannot be unrefined unless
  // all brothers want to be unrefined.
  // Loop over all elements again and let the first set of sons check if their
  // brothers also want to be unrefined
  unsigned n_unrefine=0;
  for (unsigned long e=0;e<Nelement;e++)
   {
    //(Cast) pointer to the element
    RefineableElement* el_pt = 
     dynamic_cast<RefineableElement*>(this->element_pt(e));
    
    // hierher: This is a bit naughty... We want to put the
    // first son in charge -- the statement below assumes (correctly) that the
    // enumeration of all (!) trees starts with son types.
    // This is correct for oc and quadtrees but will bite us if we
    // ever introduce other trees if/when we accidentally break this 
    // tacit assumption. Not sure what to do about it for
    // now other than leaving it hierher-ed...
    if (el_pt->tree_pt()->son_type()==OcTreeNames::LDB)
     {
      // Do all sons want to be unrefined?
      bool unrefine=true;
      unsigned n_sons = el_pt->tree_pt()->father_pt()->nsons();
      for (unsigned ison=0;ison<n_sons;ison++)
       {
        if (!(wants_to_be_unrefined[
               dynamic_cast<RefineableElement*>
               (el_pt->tree_pt()->father_pt()->
                son_pt(ison)->object_pt())]))
         {
          // One guy isn't cooperating and spoils the party.
          unrefine=false;
         }
       }

      // Tell father that his sons need to be merged
      if(unrefine)
       {
        el_pt->tree_pt()->father_pt()->object_pt()->
         select_sons_for_unrefinement();
        n_unrefine += n_sons; 
       }
      //Otherwise mark the sons as not to be touched
      else
       {
        el_pt->tree_pt()->father_pt()->object_pt()->
         deselect_sons_for_unrefinement();
       }
     }
   }
  oomph_info << " \n Number of elements to be merged : " 
             << n_unrefine << std::endl << std::endl;


  //Now do the actual mesh adaptation
  //---------------------------------

  // Check whether its worth our while
  // Either some elements want to be refined, 
  // or the number that want to be unrefined are greater than the
  // specified tolerance

  // In a parallel job, it is possible that one process may not have
  // any elements to refine, BUT a neighbouring process may refine an
  // element which changes the hanging status of a node that is on
  // both processes (i.e. a halo(ed) node).  To get around this issue, 
  // ALL processes need to call adapt_mesh if ANY refinement is to 
  // take place anywhere.

  unsigned total_n_refine=0;
#ifdef OOMPH_HAS_MPI
  // Sum n_refine across all processors
  if (this->is_mesh_distributed())
   {
    MPI_Allreduce(&n_refine,&total_n_refine,1,MPI_UNSIGNED,MPI_SUM,
                  Comm_pt->mpi_comm());
   }
  else 
   {
    total_n_refine=n_refine;
   }
#else
  total_n_refine=n_refine;
#endif

  // There may be some issues with unrefinement too, but I have not
  // been able to come up with an example (either in my head or in a
  // particular problem) where anything has arisen. I can see that
  // there may be an issue if n_unrefine differs across processes so
  // that (total_n_unrefine > max_keep_unrefined()) on some but not 
  // all processes. I haven't seen any examples of this yet so the 
  // following code may or may not work!  (Andy, 06/03/08)

  unsigned total_n_unrefine=0;
#ifdef OOMPH_HAS_MPI
  // Sum n_unrefine across all processors
  if (this->is_mesh_distributed())
   {
    MPI_Allreduce(&n_unrefine,&total_n_unrefine,1,MPI_UNSIGNED,MPI_SUM,
                  Comm_pt->mpi_comm());
   }
  else
   {
    total_n_unrefine=n_unrefine;
   }
#else
  total_n_unrefine=n_unrefine;
#endif

  oomph_info << "---> " << total_n_refine << " elements to be refined, and "
             << total_n_unrefine << " to be unrefined, in total.\n" << std::endl;

  if ((total_n_refine > 0) || (total_n_unrefine > max_keep_unrefined()))
   {

#ifdef PARANOID
#ifdef OOMPH_HAS_MPI
    

    // Sanity check: Each processor checks if the enforced unrefinement of
    // its haloed element is matched by enforced unrefinement of the
    // corresponding halo elements on the other processors.
    if (this->is_mesh_distributed())
     {
      // Store number of processors and current process
      MPI_Status status;
      int n_proc=Comm_pt->nproc();
      int my_rank=Comm_pt->my_rank();
    
      // Loop over processes: Each processor sends unrefinement pattern
      // for halo elements with processor d to processor d where it's
      // compared against the unrefinement pattern for the corresponding
      // haloed elements
      for (int d=0; d<n_proc; d++)
       {
        
        // No halo with self: Send unrefinement info to proc d
        if (d!=my_rank) 
         {
          // Get the vector of halo elements whose non-halo counterpart
          // are on processor d
          Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
          
          // Create vector containing (0)1 to indicate that
          // halo element is (not) to be unrefined
          unsigned nhalo=halo_elem_pt.size();
          Vector<int> halo_to_be_unrefined(nhalo,0);
          for (unsigned e=0;e<nhalo;e++)
           {
            if (dynamic_cast<RefineableElement*>(halo_elem_pt[e])
                ->sons_to_be_unrefined())
             {
              halo_to_be_unrefined[e]=1;
             }
           }
          
          //Trap the case when there are no halo elements
          //so that we don't get a segfault in the MPI send
          if(nhalo > 0)
           {
            // Send it across to proc d
            MPI_Send(&halo_to_be_unrefined[0],nhalo,MPI_INT,
                     d,0,Comm_pt->mpi_comm());
           }
         }
        // else (d=my_rank): Receive unrefinement pattern from all
        // other processors (dd)
        else
         {
          // Loop over other processors
          for (int dd=0; dd<n_proc; dd++)
           {
            // No halo with yourself
            if (dd!=d)
             {
              // Get the vector of haloed elements on current processor
              // with processor dd
              Vector<GeneralisedElement*> 
               haloed_elem_pt(this->haloed_element_pt(dd));
              
              // Ask processor dd to send vector containing (0)1 for 
              // halo element with current processor to be (not)unrefined
              unsigned nhaloed=haloed_elem_pt.size();
              Vector<int> halo_to_be_unrefined(nhaloed);
              //Trap to catch the case that there are no haloed elements
              if(nhaloed > 0)
               {
                // Receive unrefinement pattern of haloes from proc dd
                MPI_Recv(&halo_to_be_unrefined[0],nhaloed,MPI_INT,dd,0,
                         Comm_pt->mpi_comm(),&status);
               }
              
              // Check it
              for (unsigned e=0;e<nhaloed;e++)
               {
                if ( ( (halo_to_be_unrefined[e]==0)&&
                       (dynamic_cast<RefineableElement*>(haloed_elem_pt[e])->
                        sons_to_be_unrefined()) ) ||
                     ( (halo_to_be_unrefined[e]==1)&&
                       (!dynamic_cast<RefineableElement*>(haloed_elem_pt[e])->
                        sons_to_be_unrefined()) ) )
                 {
                  std::ostringstream error_message;
                  error_message 
                   << "Error in unrefinement: \n"
                   << "Haloed element: " << e << " on proc " << my_rank<< " \n"
                   << "wants to be unrefined whereas its halo counterpart on\n"
                   << "proc " << dd << " doesn't (or vice versa)...\n"
                   << "This is most likely because the error estimator\n"
                   << "has not assigned the same errors to halo and haloed\n"
                   << "elements -- it ought to!\n";
                  throw OomphLibError(
                   error_message.str(),
                   OOMPH_CURRENT_FUNCTION,
                   OOMPH_EXCEPTION_LOCATION);
                 }
               }
             }
           }       
         }
       }



      // Loop over processes: Each processor sends refinement pattern
      // for halo elements with processor d to processor d where it's
      // compared against the refinement pattern for the corresponding
      // haloed elements
      for (int d=0; d<n_proc; d++)
       {
        
        // No halo with self: Send refinement info to proc d
        if (d!=my_rank) 
         {
          // Get the vector of halo elements whose non-halo counterpart
          // are on processor d
          Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
          
          // Create vector containing (0)1 to indicate that
          // halo element is (not) to be refined
          unsigned nhalo=halo_elem_pt.size();
          Vector<int> halo_to_be_refined(nhalo,0);
          for (unsigned e=0;e<nhalo;e++)
           {
            if (dynamic_cast<RefineableElement*>(halo_elem_pt[e])
                ->to_be_refined())
             {
              halo_to_be_refined[e]=1;
             }
           }
          
          //Trap the case when there are no halo elements
          //so that we don't get a segfault in the MPI send
          if(nhalo > 0)
           {
            // Send it across to proc d
            MPI_Send(&halo_to_be_refined[0],nhalo,MPI_INT,
                     d,0,Comm_pt->mpi_comm());
           }
         }
        // else (d=my_rank): Receive refinement pattern from all
        // other processors (dd)
        else
         {
          // Loop over other processors
          for (int dd=0; dd<n_proc; dd++)
           {
            // No halo with yourself
            if (dd!=d)
             {
              // Get the vector of haloed elements on current processor
              // with processor dd
              Vector<GeneralisedElement*> 
               haloed_elem_pt(this->haloed_element_pt(dd));
              
              // Ask processor dd to send vector containing (0)1 for 
              // halo element with current processor to be (not)refined
              unsigned nhaloed=haloed_elem_pt.size();
              Vector<int> halo_to_be_refined(nhaloed);
              //Trap to catch the case that there are no haloed elements
              if(nhaloed > 0)
               {
                // Receive unrefinement pattern of haloes from proc dd
                MPI_Recv(&halo_to_be_refined[0],nhaloed,MPI_INT,dd,0,
                         Comm_pt->mpi_comm(),&status);
               }
              
              // Check it
              for (unsigned e=0;e<nhaloed;e++)
               {
                if ( ( (halo_to_be_refined[e]==0)&&
                       (dynamic_cast<RefineableElement*>(haloed_elem_pt[e])->
                        to_be_refined()) ) ||
                     ( (halo_to_be_refined[e]==1)&&
                       (!dynamic_cast<RefineableElement*>(haloed_elem_pt[e])->
                        to_be_refined()) ) )
                 {
                  std::ostringstream error_message;
                  error_message 
                   << "Error in refinement: \n"
                   << "Haloed element: " << e << " on proc " << my_rank<< " \n"
                   << "wants to be refined whereas its halo counterpart on\n"
                   << "proc " << dd << " doesn't (or vice versa)...\n"
                   << "This is most likely because the error estimator\n"
                   << "has not assigned the same errors to halo and haloed\n"
                   << "elements -- it ought to!\n";
                  throw OomphLibError(
                   error_message.str(),
                   OOMPH_CURRENT_FUNCTION,
                   OOMPH_EXCEPTION_LOCATION);
                 } 
               } 
             }
           }       
         }
       }
     }

#endif
#endif

    // Perform the actual adaptation
    adapt_mesh(local_doc_info);

    // The number of refineable elements is still local to each process
    Nunrefined=n_unrefine;
    Nrefined=n_refine;
   }
  // If not worthwhile, say so but still reorder nodes and kill 
  // external storage for consistency in parallel computations
  else
   {

#ifdef OOMPH_HAS_MPI
    // Delete any external element storage - any interaction will still
    // be set up on the fly again, so we need to get rid of old information.
    // This particularly causes problems in multi-domain examples where
    // we decide not to refine one of the meshes
    this->delete_all_external_storage();
#endif

    // Reorder the nodes within the mesh's node vector
    // to establish a standard ordering regardless of the sequence
    // of mesh refinements -- this is required to allow dump/restart
    // on refined meshes
    this->reorder_nodes();

#ifdef OOMPH_HAS_MPI

    // Now (re-)classify halo and haloed nodes and synchronise hanging
    // nodes
    // This is required in cases where delete_all_external_storage()
    // made slave nodes to external halo nodes nonhanging.
    if (this->is_mesh_distributed())
     {
      DocInfo doc_info;
      doc_info.disable_doc();
      classify_halo_and_haloed_nodes(doc_info,doc_info.is_doc_enabled());
     } 

#endif

    if (n_refine==0)
     {
      oomph_info 
       << " Not enough benefit in adapting mesh."
       << std::endl << std::endl;
     }
    Nunrefined=0;
    Nrefined=0;
   }


 }

//========================================================================
/// Get max/min refinement level
//========================================================================
void TreeBasedRefineableMeshBase:: 
get_refinement_levels(unsigned& min_refinement_level,
                      unsigned& max_refinement_level)
{ 
 // Initialise
 min_refinement_level=UINT_MAX;
 max_refinement_level=0;
 
 // Loop over all elements
 unsigned long n_element=this->nelement();
 if (n_element==0)
  {
   min_refinement_level=0;
   max_refinement_level=0;
  }
 else
  {
   for (unsigned long e=0;e<n_element;e++)
    {
     //Get the refinement level of the element
     unsigned level =
      dynamic_cast<RefineableElement*>(this->element_pt(e))->refinement_level();
     
     if (level>max_refinement_level) max_refinement_level=level;
     if (level<min_refinement_level) min_refinement_level=level;
    }
  }

}


//================================================================
/// Adapt mesh, which exists in two representations,
/// namely as:
///  - a FE mesh 
///  - a forest of Oc or QuadTrees 
///
/// Refinement/derefinement process is documented (in tecplot-able form)
/// if requested. 
///
/// Procedure:
/// - Loop over all elements and do the refinement for those who want to
///   be refined. Note: Refinement/splitting only allocates elements but 
///   doesn't build them.
/// - Build the new elements (i.e. give them nodes (create new ones where
///   necessary), assign boundary conditions, and add nodes to mesh
///   and mesh boundaries.
/// - For all nodes that were hanging on the previous mesh (and are still
///   marked as such), fill in their nodal values (consistent
///   with the current hanging node scheme) to make sure they are fully 
///   functional, should they have become non-hanging during the 
///   mesh-adaptation. Then mark the nodes as non-hanging.
/// - Unrefine selected elements (which may cause nodes to be re-built).
/// - Add the new elements to the mesh (by completely overwriting 
///   the old Vector of elements).
/// - Delete any nodes that have become obsolete.
/// - Mark up hanging nodes and setup hanging node scheme (incl.
///   recursive cleanup for hanging nodes that depend on other
///   hanging nodes).
/// - Adjust position of hanging nodes to make sure their position
///   is consistent with the FE-based represenetation of their larger 
///   neighbours.
/// - run a quick self-test on the neighbour finding scheme and
///   check the integrity of the elements (if PARANOID)
/// - doc hanging node status, boundary conditions, neighbour
///   scheme if requested.
///
///
/// After adaptation, all nodes (whether new or old) have up-to-date
/// current and previous values. 
///
/// If refinement process is being documented, the following information
/// is documented:
/// - The files
///   - "neighbours.dat" 
///   - "all_nodes.dat" 
///   - "new_nodes.dat" 
///   - "hang_nodes_*.dat" 
///     where the * denotes a direction (n,s,e,w) in 2D 
///     or (r,l,u,d,f,b) in 3D
///. 
///   can be viewed with 
///   - QHangingNodes.mcr
///   .
/// - The file
///    - "hangnodes_withmasters.dat"
///    .
///    can be viewed with 
///    - QHangingNodesWithMasters.mcr
///    .
///    to check the hanging node status.
/// - The neighbour status of the elements is documented in
///   - "neighbours.dat"
///   .
///   and can be viewed with 
///   - QuadTreeNeighbours.mcr
///   .
//=================================================================
void TreeBasedRefineableMeshBase::adapt_mesh(DocInfo& doc_info)
{

#ifdef OOMPH_HAS_MPI
 // Delete any external element storage before performing the adaptation
 // (in particular, external halo nodes that are on mesh boundaries)
 this->delete_all_external_storage();
#endif

 //Only perform the adapt step if the mesh has any elements.  This is relevant
 //in a distributed problem with multiple meshes, where a particular
 // process may not have any elements on a particular submesh.
 if (this->nelement()>0)
  {
   // Pointer to mesh needs to be passed to some functions
   Mesh* mesh_pt=this;
 
   double t_start = 0.0;
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_start=TimingHelpers::timer();
    }

   // Do refinement(=splitting) of elements that have been selected
   // This function encapsulates the template parameter
   this->split_elements_if_required();


   if (Global_timings::Doc_comprehensive_timings)
    {
     double t_end = TimingHelpers::timer();
     oomph_info << "Time for split_elements_if_required: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   // Now elements have been created -- build all the leaves
   //-------------------------------------------------------
   //Firstly put all the elements into a vector
   Vector<Tree*> leaf_nodes_pt;
   Forest_pt->stick_leaves_into_vector(leaf_nodes_pt);


   if (Global_timings::Doc_comprehensive_timings)
    {
     double t_end = TimingHelpers::timer();
     oomph_info << "Time for stick_leaves_into_vector: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   //If we are documenting the output, create the filename
   std::ostringstream fullname;
   std::ofstream new_nodes_file;
   if(doc_info.is_doc_enabled())
    {
     fullname << doc_info.directory() << "/new_nodes" 
              << doc_info.number() << ".dat";
     new_nodes_file.open(fullname.str().c_str());
    }


  
   // Build all elements and store vector of pointers to new nodes
   // (Note: build() checks if the element has been built 
   // already, i.e. if it's not a new element).
   Vector<Node*> new_node_pt;
   bool was_already_built;
   unsigned long num_tree_nodes=leaf_nodes_pt.size();
   for (unsigned long e=0;e<num_tree_nodes;e++)
    {
     // Pre-build must be performed before any elements are built
     leaf_nodes_pt[e]->object_pt()->pre_build(mesh_pt,new_node_pt);
    }
   for (unsigned long e=0;e<num_tree_nodes;e++)
    {
     // Now do the actual build of the new elements
     leaf_nodes_pt[e]->object_pt()
      ->build(mesh_pt,new_node_pt,was_already_built,new_nodes_file);
    }


   double t_end=0.0;
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for building " << num_tree_nodes << " new elements: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
  
   //Close the new nodes files, if it was opened
   if(doc_info.is_doc_enabled()) {new_nodes_file.close();}

   // Loop over all nodes in mesh and free the dofs of those that were
   //-----------------------------------------------------------------
   // pinned only because they were hanging nodes. Also update their
   //-----------------------------------------------------------------
   // nodal values so that they contain data that is consistent
   //----------------------------------------------------------
   // with the hanging node representation
   //-------------------------------------
   // (Even if the nodal data isn't actually accessed because the node 
   // is still hanging -- we don't know this yet, and this step makes
   // sure that all nodes are fully functional and up-to-date, should
   // they become non-hanging below).
   //
   //
   // However, if we have a fixed mesh and hanging nodes on the boundary 
   // become non-hanging they will not necessarily respect the curvilinear
   // boundaries. This can only happen in 3D of course because it is not
   // possible to have hanging nodes on boundaries in 2D.
   // The solution is to store those nodes on the boundaries that are
   // currently hanging and then check to see whether they have changed
   // status at the end of the refinement procedure.
   // If it has changed, then we need to adjust their positions.
   const unsigned n_boundary = this->nboundary();
   const unsigned mesh_dim = this->finite_element_pt(0)->dim();
   Vector<std::set<Node*> > hanging_nodes_on_boundary_pt(n_boundary);

   unsigned long n_node=this->nnode();
   for (unsigned long n=0;n<n_node;n++)
    {
     //Get the pointer to the node
     Node* nod_pt=this->node_pt(n);
  
     //Get the number of values in the node
     unsigned n_value=nod_pt->nvalue();
   
     //We need to find if any of the values are hanging
     bool is_hanging = nod_pt->is_hanging();
     //Loop over the values and find out whether any are hanging
     for(unsigned n=0;n<n_value;n++)
      {is_hanging |= nod_pt->is_hanging(n);}

     //If the node is hanging then ...
     if(is_hanging)
      {       
       // Unless they are turned into hanging nodes again below
       // (this might or might not happen), fill in all the necessary
       // data to make them 'proper' nodes again.
     
       // Reconstruct the nodal values/position from the node's 
       // hanging node representation
       unsigned nt=nod_pt->ntstorage();
       Vector<double> values(n_value);
       unsigned n_dim=nod_pt->ndim();
       Vector<double> position(n_dim);
       // Loop over all history values
       for(unsigned t=0;t<nt;t++)
        {
         nod_pt->value(t,values);
         for(unsigned i=0;i<n_value;i++) 
          {
           nod_pt->set_value(t,i,values[i]);
          }
         nod_pt->position(t,position);
         for(unsigned i=0;i<n_dim;i++) 
          {
           nod_pt->x(t,i)=position[i];
          }
        }

       // If it's an algebraic node: Update its previous nodal positions too
       AlgebraicNode* alg_node_pt=dynamic_cast<AlgebraicNode*>(nod_pt);
       if (alg_node_pt!=0)
        {
         bool update_all_time_levels=true;
         alg_node_pt->node_update(update_all_time_levels);
        }
     
     
       //If it's a Solid node, update Lagrangian coordinates
       // from its hanging node representation
       SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
       if(solid_node_pt!=0)
        {
         unsigned n_lagrangian = solid_node_pt->nlagrangian();
         for(unsigned i=0;i<n_lagrangian;i++)
          {
           solid_node_pt->xi(i) = solid_node_pt->lagrangian_position(i);
          }
        }

       //Now store geometrically hanging nodes on boundaries that
       //may need updating after refinement.
       //There will only be a problem if we have 3 spatial dimensions
       if((mesh_dim > 2) && (nod_pt->is_hanging()))
        {
         //If the node is on a boundary then add a pointer to the node
         //to our lookup scheme
         if(nod_pt->is_on_boundary())
          {
           //Storage for the boundaries on which the Node is located
           std::set<unsigned>* boundaries_pt;
           nod_pt->get_boundaries_pt(boundaries_pt);
           if(boundaries_pt!=0)
            {
             //Loop over the boundaries and add a pointer to the node
             //to the appropriate storage scheme
             for(std::set<unsigned>::iterator it=boundaries_pt->begin();
                 it!=boundaries_pt->end();++it)
              {
               hanging_nodes_on_boundary_pt[*it].insert(nod_pt);
              }
            }
          }
        }

      } //End of is_hanging

     // Initially mark all nodes as 'non-hanging' and `obsolete' 
     nod_pt->set_nonhanging();
     nod_pt->set_obsolete();
    }

   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for sorting out initial hanging status: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   // Unrefine all the selected elements: This needs to be
   //-----------------------------------------------------
   // all elements, because the father elements are not actually leaves.
   //-------------------------------------------------------------------

   // Unrefine
   for (unsigned long e=0;e<Forest_pt->ntree();e++)
    {
     Forest_pt->tree_pt(e)->traverse_all(&Tree::merge_sons_if_required,
                                         mesh_pt);
    }

   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for unrefinement: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   // Add the newly created elements to mesh
   //---------------------------------------

   // Stick all elements into a new vector
   //(note the leaves may have changed, so this is not duplicated work)
   Vector<Tree*> tree_nodes_pt;
   Forest_pt->stick_leaves_into_vector(tree_nodes_pt);

   //Copy the elements into the mesh Vector
   num_tree_nodes=tree_nodes_pt.size();
   Element_pt.resize(num_tree_nodes);
   for (unsigned long e=0;e<num_tree_nodes;e++)
    {
     Element_pt[e]=tree_nodes_pt[e]->object_pt();

     // Now loop over all nodes in element and mark them as non-obsolete
     // Logic: Initially all nodes in the unrefined mesh were labeled
     // as deleteable. Then we create new elements (whose newly created
     // nodes are obviously non-obsolete), and killed some other elements (by
     // by deleting them and marking the nodes that were not shared by
     // their father as obsolete. Now we loop over all the remaining
     // elements and (re-)label all their nodes as non-obsolete. This
     // saves some nodes that were regarded as obsolete by deleted
     // elements but are still required in some surviving ones
     // from a tragic early death...
     FiniteElement* this_el_pt=this->finite_element_pt(e);
     unsigned n_node=this_el_pt->nnode(); // caching pre-loop
     for (unsigned n=0;n<n_node;n++)
      {
       this_el_pt->node_pt(n)->set_non_obsolete();
      }
    }


   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for adding elements to mesh: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
  
   // Cannot delete nodes that are still marked as obsolete
   // because they may still be required to assemble the hanging schemes
   //-------------------------------------------------------------------

   // Mark up hanging nodes
   //----------------------
 
   //Output streams for the hanging nodes
   Vector<std::ofstream*> hanging_output_files;
   //Setup the output files for hanging nodes, this must be called
   //precisely once for the forest. Note that the files will only
   //actually be opened if doc_info.is_doc_enabled() is true
   Forest_pt->open_hanging_node_files(doc_info,hanging_output_files);

   for(unsigned long e=0;e<num_tree_nodes;e++)
    {
     //Generic setup
     tree_nodes_pt[e]->object_pt()->setup_hanging_nodes(hanging_output_files);
     //Element specific setup
     tree_nodes_pt[e]->object_pt()->further_setup_hanging_nodes();
    }

   //Close the hanging node files and delete the memory allocated 
   //for the streams
   Forest_pt->close_hanging_node_files(doc_info,hanging_output_files);


   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info 
     <<"Time for setup_hanging_nodes() and further_setup_hanging_nodes() for " 
     << num_tree_nodes << " elements: "
     << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   // Read out the number of continously interpolated values
   // from one of the elements (assuming it's the same in all elements)
   unsigned ncont_interpolated_values=
    tree_nodes_pt[0]->object_pt()->ncont_interpolated_values();

   // Complete the hanging nodes schemes by dealing with the
   // recursively hanging nodes
   complete_hanging_nodes(ncont_interpolated_values);
  
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info 
      <<"Time for complete_hanging_nodes: "
      << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   /// Update the boundary element info -- this can be a costly procedure
   /// and for this reason the mesh writer might have decided not to 
   /// set up this scheme. If so, we won't change this and suppress 
   /// its creation...
   if (Lookup_for_elements_next_boundary_is_setup)
    {
     this->setup_boundary_element_info(); 
    }

   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info
      <<"Time for boundary element info: " 
      << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
  
#ifdef PARANOID
 
   // Doc/check the neighbours
   //-------------------------
   Vector<Tree*> all_tree_nodes_pt;
   Forest_pt->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);

   //Check the neighbours
   Forest_pt->check_all_neighbours(doc_info);
 
   // Check the integrity of the elements
   // -----------------------------------
 
   // Loop over elements and get the elemental integrity
   double max_error=0.0;
   for (unsigned long e=0;e<num_tree_nodes;e++)
    {
     double max_el_error;
     tree_nodes_pt[e]->object_pt()->check_integrity(max_el_error);
     //If the elemental error is greater than our maximum error
     //reset the maximum
     if(max_el_error > max_error) {max_error=max_el_error;}
    }

   if (max_error>RefineableElement::max_integrity_tolerance())
    {
     std::ostringstream error_stream;
     error_stream << "Mesh refined: Max. error in integrity check: " 
                  << max_error << " is too big\n";
     error_stream
      << "i.e. bigger than RefineableElement::max_integrity_tolerance()="
      << RefineableElement::max_integrity_tolerance() << "\n" << std::endl;

     std::ofstream some_file;
     some_file.open("ProblemMesh.dat");
     for (unsigned long n=0;n<n_node;n++)
      {
       //Get the pointer to the node
       Node* nod_pt = this->node_pt(n);
       //Get the dimension
       unsigned n_dim = nod_pt->ndim();
       //Output the coordinates
       for(unsigned i=0;i<n_dim;i++)
        {
         some_file << this->node_pt(n)->x(i) << " ";
        }
       some_file << std::endl;
      }
     some_file.close();

     error_stream << "Doced problem mesh in ProblemMesh.dat" << std::endl;

     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   else
    {
     oomph_info << "Mesh refined: Max. error in integrity check: " 
                << max_error << " is OK" << std::endl;
     oomph_info 
      << "i.e. less than RefineableElement::max_integrity_tolerance()="
      << RefineableElement::max_integrity_tolerance() << "\n" << std::endl;
    }
 

   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info
      <<"Time for (paranoid only) checking of integrity: " 
      << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

#endif
  
   //Loop over all elements other than the final level and deactivate the
   //objects, essentially set the pointer that point to nodes that are
   //about to be deleted to NULL. This must take place here because nodes
   //addressed by elements that are dead but still living in the tree might
   //have been made obsolete in the last round of refinement
   for (unsigned long e=0;e<Forest_pt->ntree();e++)
    {
     Forest_pt->tree_pt(e)->
      traverse_all_but_leaves(&Tree::deactivate_object);
    }
 
   //Now we can prune the dead nodes from the mesh.
   Vector<Node*> deleted_node_pt=this->prune_dead_nodes();

   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for deactivating objects and pruning nodes: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   // Finally: Reorder the nodes within the mesh's node vector
   // to establish a standard ordering regardless of the sequence
   // of mesh refinements -- this is required to allow dump/restart
   // on refined meshes
   this->reorder_nodes();
  
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for reordering " << nnode() << " nodes: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
  
   //Now we can correct the nodes on boundaries that were hanging that
   //are no longer hanging
   //Only bother if we have more than two dimensions
   if(mesh_dim > 2)
    {
     //Loop over the boundaries
     for(unsigned b=0;b<n_boundary;b++)
      {

       // Remove deleted nodes from the set
       unsigned n_del=deleted_node_pt.size();
       for (unsigned j=0;j<n_del;j++)
        {
         hanging_nodes_on_boundary_pt[b].erase(deleted_node_pt[j]);
        }
       
       //If the nodes that were hanging are still hanging then remove them
       //from the set (note increment is not in for command for efficiencty)
       for(std::set<Node*>::iterator 
            it=hanging_nodes_on_boundary_pt[b].begin();
            it!=hanging_nodes_on_boundary_pt[b].end();)
        {
         if((*it)->is_hanging()) 
          {hanging_nodes_on_boundary_pt[b].erase(it++);}
         else {++it;}
        }

       //Are there any nodes that have changed geometric hanging status
       //on the boundary
       //The slightly painful part is that we must adjust the position
       //via the macro-elements which are only available through the
       //elements and not the nodes.
       if(hanging_nodes_on_boundary_pt[b].size()>0)
        {
         //If so we loop over all elements adjacent to the boundary
         unsigned n_boundary_element = this->nboundary_element(b);
         for(unsigned e=0;e<n_boundary_element;++e)
          {
           //Get a pointer to the element
           FiniteElement* el_pt = this->boundary_element_pt(b,e);

           //Do we have a solid element
           SolidFiniteElement* solid_el_pt = 
            dynamic_cast<SolidFiniteElement*>(el_pt);
           
           //Determine whether there is a macro element
           bool macro_present = (el_pt->macro_elem_pt()!=0);
           //Or a solid macro element
           if(solid_el_pt!=0)
            {
             macro_present |= (solid_el_pt->undeformed_macro_elem_pt()!=0);
            }
           
           //Only bother to do anything if there is a macro element
           //or undeformed macro element in a SolidElement
           if(macro_present)
            {
             //Loop over the nodes 
             //ALH: (could optimise to only loop over
             //node associated with the boundary with more effort)
             unsigned n_el_node = el_pt->nnode();
             for(unsigned n=0;n<n_el_node;n++)
              {
               //Cache pointer to the node
               Node* nod_pt = el_pt->node_pt(n);
               if(nod_pt->is_on_boundary(b))
                {
                 //Is the Node in our set
                 std::set<Node*>::iterator it = 
                  hanging_nodes_on_boundary_pt[b].find(nod_pt);

                 //If we have found the Node then update the position
                 //to be consistent with the macro-element representation
                 if(it!=hanging_nodes_on_boundary_pt[b].end())
                  {
                   //Specialise local and global coordinates to 3D
                   //because there is only a problem in 3D.
                   Vector<double> s(3), x(3);

                   //Find the local coordinate of the ndoe
                   el_pt->local_coordinate_of_node(n,s);

                   //Find the number of time history values
                   const unsigned ntstorage = nod_pt->ntstorage();
                   
                   //Do we have a solid node
                   SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
                   if(solid_node_pt)
                    {
                     // Assign Lagrangian coordinates from undeformed
                     // macro element (if it has one -- get_x_and_xi()
                     // does "the right thing" anyway. Leave actual
                     // nodal positions alone -- we're doing a solid
                     // mechanics problem and once we're going
                     // the nodal positions are always computed, never
                     // (automatically) reset to macro-element based
                     // positions; not even on pinned boundaries
                     // because the user may have other ideas about where
                     // these should go -- pinning means "don't touch the
                     // value", not "leave where the macro-element thinks
                     // it should be"
                     Vector<double> x_fe(3), xi(3), xi_fe(3);
                     solid_el_pt->get_x_and_xi(s,x_fe,x,xi_fe,xi);
                     for(unsigned i=0;i<3;i++) {solid_node_pt->xi(i) = xi[i];}
                    }
                   else
                    {
                     //Set position and history values from the macro-element
                     //representation
                     for(unsigned t=0;t<ntstorage;t++)
                      {
                       //Get the history value from the macro element
                       el_pt->get_x(t,s,x);

                       //Set the coordinate to that of the macroelement
                       //representation
                       for(unsigned i=0;i<3;i++) {nod_pt->x(t,i) = x[i];}
                      }
                    } //End of non-solid node case
                   
                   //Now remove the node from the list
                   hanging_nodes_on_boundary_pt[b].erase(it);

                   //If there are no Nodes left then exit the loops
                   if(hanging_nodes_on_boundary_pt[b].size()==0)
                    {e=n_boundary_element; break;}
                  }
                }
              }
            } //End of macro element case
          }
        }
      }
    } //End of case when we have fixed nodal positions

  
   // Final doc
   //-----------
   if (doc_info.is_doc_enabled())
    {
     // Doc the boundary conditions ('0' for non-existent, '1' for free, 
     //----------------------------------------------------------------
     // '2' for pinned -- ideal for tecplot scatter sizing.
     //----------------------------------------------------
     //num_tree_nodes=tree_nodes_pt.size();
   
     // Determine maximum number of values at any node in this type of element
     RefineableElement* el_pt = tree_nodes_pt[0]->object_pt();
     //Initalise max_nval
     unsigned max_nval=0;
     for (unsigned n=0;n<el_pt->nnode();n++)
      {
       if (el_pt->node_pt(n)->nvalue()>max_nval)
        {max_nval=el_pt->node_pt(n)->nvalue();}
      }
   
     //Open the output file
     std::ofstream bcs_file;
     fullname.str("");
     fullname << doc_info.directory() << "/bcs" << doc_info.number()
              << ".dat";
     bcs_file.open(fullname.str().c_str());  
   
     // Loop over elements
     for(unsigned long e=0;e<num_tree_nodes;e++)
      {
       el_pt = tree_nodes_pt[e]->object_pt();
       // Loop over nodes in element
       unsigned n_nod=el_pt->nnode();
       for(unsigned n=0;n<n_nod;n++)
        {
         //Get pointer to the node
         Node* nod_pt=el_pt->node_pt(n);
         //Find the dimension of the node
         unsigned n_dim = nod_pt->ndim();
         //Write the nodal coordinates to the file
         for(unsigned i=0;i<n_dim;i++)
          {bcs_file << nod_pt->x(i) << " ";}
       
         // Loop over all values in this element
         for(unsigned i=0;i<max_nval;i++)
          {
           // Value exists at this node:
           if (i<nod_pt->nvalue())
            {
             bcs_file << " " << 1+nod_pt->is_pinned(i);
            }
           // ...if not just dump out a zero
           else
            {
             bcs_file << " 0 ";
            }
          }
         bcs_file << std::endl;
        }  
      }
     bcs_file.close();
   
     // Doc all nodes
     //---------------
     std::ofstream all_nodes_file;
     fullname.str("");
     fullname << doc_info.directory() << "/all_nodes"
              << doc_info.number()  << ".dat";
     all_nodes_file.open(fullname.str().c_str());  
   
     all_nodes_file << "ZONE \n"; 
   
     // Need to recompute the number of nodes since it may have
     // changed during mesh refinement/unrefinement
     n_node = this->nnode();
     for(unsigned long n=0;n<n_node;n++)
      {
       Node* nod_pt = this->node_pt(n);
       unsigned n_dim = nod_pt->ndim();
       for(unsigned i=0;i<n_dim;i++)
        {
         all_nodes_file << this->node_pt(n)->x(i) << " ";
        }
       all_nodes_file << std::endl;
      }
   
     all_nodes_file.close();
   
   
     // Doc all hanging nodes:
     //-----------------------
     std::ofstream some_file;
     fullname.str("");
     fullname << doc_info.directory() << "/all_hangnodes"
              << doc_info.number() << ".dat";
     some_file.open(fullname.str().c_str());
     for(unsigned long n=0;n<n_node;n++)
      {
       Node* nod_pt=this->node_pt(n);
     
       if (nod_pt->is_hanging())
        {
         unsigned n_dim = nod_pt->ndim();       
         for(unsigned i=0;i<n_dim;i++)
          {
           some_file << nod_pt->x(i) << " ";
          }
       
         //ALH: Added this to stop Solid problems seg-faulting
         if(this->node_pt(n)->nvalue() > 0)
          {
           some_file  << " " << nod_pt->raw_value(0);
          }
         some_file << std::endl;
        }
      }
     some_file.close();
   
     // Doc all hanging nodes and their masters 
     // View with QHangingNodesWithMasters.mcr
     fullname.str("");
     fullname << doc_info.directory() 
              << "/geometric_hangnodes_withmasters" 
              << doc_info.number() << ".dat";
     some_file.open(fullname.str().c_str());
     for(unsigned long n=0;n<n_node;n++)
      {
       Node* nod_pt=this->node_pt(n);
       if (nod_pt->is_hanging())
        {
         unsigned n_dim = nod_pt->ndim();
         unsigned nmaster=nod_pt->hanging_pt()->nmaster();
         some_file << "ZONE I="<<nmaster+1 << std::endl;
         for(unsigned i=0;i<n_dim;i++)
          {
           some_file << nod_pt->x(i) << " ";
          }
         some_file << " 2 " <<  std::endl;
       
         for (unsigned imaster=0;imaster<nmaster;imaster++)
          {
           Node* master_nod_pt = 
            nod_pt->hanging_pt()->master_node_pt(imaster);
           unsigned n_dim = master_nod_pt->ndim();
           for(unsigned i=0;i<n_dim;i++)
            {
             some_file << master_nod_pt->x(i) << " ";
            }
           some_file << " 1 " << std::endl;
          }
        }
      }
     some_file.close();
   
     // Doc all hanging nodes and their masters 
     // View with QHangingNodesWithMasters.mcr
     for(unsigned i=0;i<ncont_interpolated_values;i++)
      {
       fullname.str("");
       fullname << doc_info.directory()
                <<"/nonstandard_hangnodes_withmasters" << i << "_" 
                << doc_info.number() << ".dat";
       some_file.open(fullname.str().c_str());
       unsigned n_nod=this->nnode();
       for(unsigned long n=0;n<n_nod;n++)
        {
         Node* nod_pt=this->node_pt(n);
         if (nod_pt->is_hanging(i))
          {
           if (nod_pt->hanging_pt(i)!=nod_pt->hanging_pt())
            {
             unsigned nmaster=nod_pt->hanging_pt(i)->nmaster();
             some_file << "ZONE I="<<nmaster+1 << std::endl;
             unsigned n_dim = nod_pt->ndim();
             for(unsigned j=0;j<n_dim;j++)
              {
               some_file << nod_pt->x(j) << " ";
              }
             some_file << " 2 " << std::endl;
             for (unsigned imaster=0;imaster<nmaster;imaster++)
              {
               Node* master_nod_pt = 
                nod_pt->hanging_pt(i)->master_node_pt(imaster);
               unsigned n_dim = master_nod_pt->ndim();
               for(unsigned j=0;j<n_dim;j++)
                {
//               some_file << master_nod_pt->x(i) << " ";
                }
               some_file << " 1 " << std::endl;
              }
            }
          }
        }
       some_file.close();
      }
    } //End of documentation
  } // End if (this->nelement()>0)


#ifdef OOMPH_HAS_MPI

 // Now (re-)classify halo and haloed nodes and synchronise hanging
 // nodes
 if (this->is_mesh_distributed())
  {
   classify_halo_and_haloed_nodes(doc_info,doc_info.is_doc_enabled());
  }

#endif

}


//========================================================================
/// Refine mesh uniformly
//========================================================================
void TreeBasedRefineableMeshBase::refine_uniformly(DocInfo& doc_info)
{ 
 //Select all elements for refinement
 unsigned long Nelement=this->nelement();
 for (unsigned long e=0;e<Nelement;e++)
  {
   dynamic_cast<RefineableElement*>
    (this->element_pt(e))->select_for_refinement();
  }
 
 // Do the actual mesh adaptation
 adapt_mesh(doc_info);  

 }


//========================================================================
/// p-refine mesh uniformly
//========================================================================
void TreeBasedRefineableMeshBase::p_refine_uniformly(DocInfo& doc_info)
{ 
 //Select all elements for refinement
 unsigned long Nelement=this->nelement();
 for (unsigned long e=0;e<Nelement;e++)
  {
   //Get pointer to p-refineable element
   PRefineableElement* el_pt
    = dynamic_cast<PRefineableElement*>(this->element_pt(e));
   //Mark for p-refinement if possible. If not then p_adapt_mesh() will
   //report the error.
   if(el_pt!=0)
    {
     el_pt->select_for_p_refinement();
    }
  }
 
 // Do the actual mesh adaptation
 p_adapt_mesh(doc_info);  
}


//========================================================================
/// Refine mesh by splitting the elements identified
/// by their numbers.
//========================================================================
void TreeBasedRefineableMeshBase::refine_selected_elements(
 const Vector<unsigned>& elements_to_be_refined)
{ 
 
#ifdef OOMPH_HAS_MPI
 if(this->is_mesh_distributed())
  {
   std::ostringstream warn_stream;
   warn_stream << "You are attempting to refine selected elements of a "
               << std::endl
               << "distributed mesh. This may have undesired effects."
               << std::endl;
   
   OomphLibWarning(warn_stream.str(),
                   "TreeBasedRefineableMeshBase::refine_selected_elements()",
                   OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 //Select elements for refinement
 unsigned long nref=elements_to_be_refined.size();
 for (unsigned long e=0;e<nref;e++)
  {
   dynamic_cast<RefineableElement*>
    (this->element_pt(elements_to_be_refined[e]))->select_for_refinement();
  }

 // Do the actual mesh adaptation
 adapt_mesh(); 

}



//========================================================================
/// Refine mesh by splitting the elements identified
/// by their pointers
//========================================================================
void TreeBasedRefineableMeshBase::refine_selected_elements(
 const Vector<RefineableElement*>& elements_to_be_refined_pt)
{ 
 
#ifdef OOMPH_HAS_MPI
 if(this->is_mesh_distributed())
  {
   std::ostringstream warn_stream;
   warn_stream << "You are attempting to refine selected elements of a "
               << std::endl
               << "distributed mesh. This may have undesired effects."
               << std::endl;
   
   OomphLibWarning(warn_stream.str(),
                   "TreeBasedRefineableMeshBase::refine_selected_elements()",
                   OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 //Select elements for refinement
 unsigned long nref=elements_to_be_refined_pt.size();
 for (unsigned long e=0;e<nref;e++)
  {
   elements_to_be_refined_pt[e]->select_for_refinement();
  }
 
 // Do the actual mesh adaptation
 adapt_mesh();  
}



//========================================================================
/// \short Refine to same degree as the reference mesh.
//========================================================================
void TreeBasedRefineableMeshBase::refine_base_mesh_as_in_reference_mesh(
 TreeBasedRefineableMeshBase* const &ref_mesh_pt)
{
 // Assign storage for refinement pattern
 Vector<Vector<unsigned> > to_be_refined;

 // Get refinement pattern of reference mesh:
 ref_mesh_pt->get_refinement_pattern(to_be_refined);

 // Refine mesh according to given refinement pattern
 refine_base_mesh(to_be_refined);
}

 
//========================================================================
/// \short Refine to same degree as the reference mesh minus one. Useful
/// function for multigrid solvers; allows the easy copy of a mesh
/// to the level of refinement just below the current one. Returns
/// a boolean variable which indicates if the reference mesh has not
/// been refined at all
//========================================================================
 bool TreeBasedRefineableMeshBase::
 refine_base_mesh_as_in_reference_mesh_minus_one(
  TreeBasedRefineableMeshBase* const &ref_mesh_pt)
 {
  // Assign storage for refinement pattern
  Vector<Vector<unsigned> > to_be_refined;

  // Get refinement pattern of reference mesh:
  ref_mesh_pt->get_refinement_pattern(to_be_refined);

  // Find the length of the vector
  unsigned nrefinement_levels=to_be_refined.size();

  // If the reference mesh has not been refined a single time then
  // we cannot create an unrefined copy so stop here
  if (nrefinement_levels==0)
  {
   return false;
  }
  // If the reference mesh has been refined at least once
  else
  {  
   // Remove the last entry of the vector to make sure we refine to
   // the same level minus one
   to_be_refined.resize(nrefinement_levels-1);  

   // Refine mesh according to given refinement pattern
   refine_base_mesh(to_be_refined);

   // Indicate that it was possible to create an unrefined copy
   return true;
  }
 }

 
//========================================================================
/// Refine mesh once so that its topology etc becomes that of the 
/// (finer!) reference mesh -- if possible! Useful for meshes in multigrid 
/// hierarchies. If the meshes are too different and the conversion
/// cannot be performed, the code dies (provided PARANOID is enabled).
//========================================================================
void TreeBasedRefineableMeshBase::refine_as_in_reference_mesh(
 TreeBasedRefineableMeshBase* const &ref_mesh_pt)
{
 oomph_info << "WARNING : This has not been checked comprehensively yet"
           << std::endl
           << "Check it and remove this break " << std::endl;
 pause("Yes really pause");

#ifdef PARANOID
 // The max. refinement levels of the two meshes need to differ 
 // by one, otherwise what we're doing here doesn't make sense.
 unsigned my_min,my_max;
 get_refinement_levels(my_min,my_max);

 unsigned ref_min,ref_max;
 ref_mesh_pt->get_refinement_levels(ref_min,ref_max);

 if (ref_max!=my_max+1)
  {
   std::ostringstream error_stream;
   error_stream
    << "Meshes definitely don't differ by one refinement level \n"
    << "max. refinement levels: "<< ref_max << " " << my_max << std::endl;

   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Vector storing the elements of the uniformly unrefined mesh
 Vector<Tree*> coarse_elements_pt;

 // Map storing which father elements have already been added to coarse mesh
 // (Default return is 0).
 std::map<Tree*,unsigned> father_element_included;

 // Extract active elements (=leaf nodes in the quadtree) from reference mesh.
 Vector<Tree*> leaf_nodes_pt;
 ref_mesh_pt->forest_pt()->
  stick_leaves_into_vector(leaf_nodes_pt);

 // Loop over all elements (in their quadtree impersonation) and
 // check if their fathers's sons are all leaves too:
 unsigned nelem=leaf_nodes_pt.size();
 for (unsigned e=0;e<nelem;e++)
  {
   //Pointer to leaf node
   Tree* leaf_pt=leaf_nodes_pt[e];

   // Get pointer to father:
   Tree* father_pt=leaf_pt->father_pt();
   
   // If we don't have a father we're at the root level in which
   // case this element can't be unrefined.
   if (0==father_pt)
    {
     coarse_elements_pt.push_back(leaf_pt);
    }
   else
    {
     // Loop over the father's sons to check if they're
     // all non-leafs, i.e. if they can be unrefined
     bool can_unrefine=true;
     unsigned n_sons = father_pt->nsons();
     for (unsigned i=0;i<n_sons;i++)
      {
       // If (at least) one of the sons is not a leaf, we can't unrefine
       if (!father_pt->son_pt(i)->is_leaf()) can_unrefine=false;
      }

     // If we can unrefine, the father element will be 
     // an element in the coarse mesh, the sons won't
     if (can_unrefine)
      {
       if (father_element_included[father_pt]==0)
        {
         coarse_elements_pt.push_back(father_pt);
         father_element_included[father_pt]=1;
        }
      }
     // Son will still be there on the coarse mesh
     else
      {
       coarse_elements_pt.push_back(leaf_pt);
      }
    }
  }

 // Number of elements in ref mesh if it was unrefined uniformly:
 unsigned nel_coarse=coarse_elements_pt.size();


#ifdef PARANOID
 bool stop_it=false;
 // The numbers had better match otherwise we might as well stop now...
 if (nel_coarse!=this->nelement())
  {
   oomph_info << "Number of elements in uniformly unrefined reference mesh: " 
        << nel_coarse<< std::endl;
   oomph_info << "Number of elements in 'this' mesh: " 
        << nel_coarse<< std::endl;
   oomph_info << "don't match" << std::endl;
   stop_it=true;
  }
#endif

 // Now loop over all elements in uniformly coarsened reference mesh
 // and check if add the number of any element that was created
 // by having had its sons merged to the vector of elements that
 // need to get refined if we go the other way
 Vector<unsigned> elements_to_be_refined;
 for (unsigned i=0;i<nel_coarse;i++)
  {
   if (father_element_included[coarse_elements_pt[i]]==1)
    {
     elements_to_be_refined.push_back(i);
    }
  }

 
#ifdef PARANOID
  // Doc troublesome meshes:
 if (stop_it)
  {
   std::ofstream some_file;
   some_file.open("orig_mesh.dat");
   this->output(some_file);
   some_file.close();
   oomph_info << "Documented original ('this')mesh in orig_mesh.dat" 
             << std::endl;
  }
#endif


 // Now refine precisely these elements in "this" mesh.
 refine_selected_elements(elements_to_be_refined);


#ifdef PARANOID

   // Check if the nodal positions of all element's nodes agree
   // in the two fine meshes:
   double tol=1.0e-5;
   for (unsigned e=0;e<nelem;e++)
    {
     // Get elements
     FiniteElement* ref_el_pt=ref_mesh_pt->finite_element_pt(e);
     FiniteElement* el_pt=this->finite_element_pt(e);

     // Loop over nodes
     unsigned nnod=ref_el_pt->nnode();
     for (unsigned j=0;j<nnod;j++)
      {
       // Get nodes
       Node* ref_node_pt=ref_el_pt->node_pt(j);
       Node* node_pt=el_pt->node_pt(j);

       // Check error in position
       double error=0.0;
       unsigned ndim=node_pt->ndim();
       for (unsigned i=0;i<ndim;i++)
        {
         error+=pow(node_pt->x(i)-ref_node_pt->x(i),2);
        }
       error=sqrt(error);

       if (error>tol)
        {
         oomph_info << "Error in nodal position of node " << j << ": " 
              << error << "           [tol=" << tol<< "]" << std::endl; 
         stop_it=true;
        }
      }
    }

   // Do we have a death wish?
   if (stop_it)
    {
     // Doc troublesome meshes:
     std::ofstream some_file;
     some_file.open("refined_mesh.dat");
     this->output(some_file);
     some_file.close();
     
     some_file.open("finer_mesh.dat");
     ref_mesh_pt->output(some_file);
     some_file.close();
 
     throw OomphLibError(
      "Bailing out. Doced refined_mesh.dat finer_mesh.dat\n",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }

#endif
   
}


//========================================================================
/// Unrefine mesh uniformly. Return 0 for success,
/// 1 for failure (if unrefinement has reached the coarsest permitted
/// level)
//========================================================================
unsigned TreeBasedRefineableMeshBase::unrefine_uniformly()
{ 

 // We can't just select all elements for unrefinement
 // because they need to merge with their brothers.
 // --> Rather than repeating the convoluted logic of
 // RefineableQuadMesh<ELEMENT>::adapt(Vector<double>& elemental_error) 
 // here (code duplication!) hack it by filling the error
 // vector with values that ensure unrefinement for all
 // elements where this is possible

 // Create dummy vector for elemental errors
 unsigned long Nelement=this->nelement();
 Vector<double> elemental_error(Nelement);

 //In order to force unrefinement, set the min permitted error to 
 //be the default and then set the actual error to be below this.
 //This avoids problems when the actual min error is zero (or small)
 //For sanity's sake, also set the max permitted error back to the default
 //so that we have a max error bigger than a min error
 const double current_min_error = this->min_permitted_error();
 const double current_max_error = this->max_permitted_error();

 this->min_permitted_error() = 1.0e-5;
 this->max_permitted_error() = 1.0e-3;

 double error=min_permitted_error()/100.0;
 for (unsigned long e=0;e<Nelement;e++)
  {
   elemental_error[e]=error;
  }

 // Temporarily lift any restrictions on the minimum number of
 // elements that need to be unrefined to make it worthwhile
 unsigned backup=max_keep_unrefined();
 max_keep_unrefined()=0;

  // Do the actual mesh adaptation with fake error vector
 adapt(elemental_error); 

 // Reset the minimum number of elements that need to be unrefined 
 // to make it worthwhile
 max_keep_unrefined()=backup;

 //Now restore the error tolerances
 this->min_permitted_error() = current_min_error;
 this->max_permitted_error() = current_max_error;

 // Has the unrefinement actually changed anything?
 if (Nelement==this->nelement())
  {
   return 1;
  }
 else
  {
   return 0;
  }

}

//==================================================================
/// Given a node, return a vector of pointers to master nodes and a
/// vector of the associated weights.
/// This is done recursively, so if a node is not hanging, 
/// the node is regarded as its own master node which has weight 1.0.
//==================================================================
void TreeBasedRefineableMeshBase::
complete_hanging_nodes_recursively(Node*& nod_pt,
                                   Vector<Node*>& master_nodes,
                                   Vector<double>& hang_weights,
                                   const int& i)
{
 // Is the node hanging in the variable i
 if(nod_pt->is_hanging(i))
  {
   // Loop over all master nodes
   HangInfo* const hang_pt = nod_pt->hanging_pt(i);
   unsigned nmaster=hang_pt->nmaster();

   for(unsigned m=0;m<nmaster;m++)
    {
     // Get the master node
     Node* master_nod_pt=hang_pt->master_node_pt(m);

     // Keep in memory the size of the list before adding the nodes this 
     // master node depends on. This is required so that the recursion is
     // only performed on these particular master nodes. A master node
     // could contain contributions from two separate pseudo-masters.
     // These contributions must be summed, not multiplied. 
     int first_new_node=master_nodes.size();

     // Now check which master nodes this master node depends on
     complete_hanging_nodes_recursively(master_nod_pt,
                                        master_nodes,hang_weights,i);
     
     // Multiply old weight by new weight for all the nodes this master 
     // node depends on
     unsigned n_new_master_node = master_nodes.size();

     double mtr_weight=hang_pt->master_weight(m);

     for(unsigned k=first_new_node;k<n_new_master_node;k++)
      {
       hang_weights[k]=mtr_weight*hang_weights[k];
      }
    }
  }
 else
  // Node isn't hanging so it enters itself with the full weight
  {
   master_nodes.push_back(nod_pt);
   hang_weights.push_back(1.0);
  }
 
}


//==================================================================
/// Complete the hanging node scheme recursively. 
/// After the initial markup scheme, hanging nodes
/// can depend on other hanging nodes ---> AAAAAAAAARGH!
/// Need to translate this into a scheme where all
/// hanging  nodes only depend on non-hanging nodes...
//==================================================================
void TreeBasedRefineableMeshBase::
complete_hanging_nodes(const int& ncont_interpolated_values)
 {
  //Number of nodes in mesh
  unsigned long n_node=this->nnode();
  double min_weight=1.0e-8; //RefineableBrickElement::min_weight_value();

  //Loop over the nodes in the mesh
  for (unsigned long n=0;n<n_node;n++)
   {
    //Assign a local pointer to the node
    Node* nod_pt=this->node_pt(n);

    //Loop over the values,
    //N.B. geometric hanging data is stored at the index -1
    for(int i=-1;i<ncont_interpolated_values;i++)
     {
      // Is the node hanging?
      if (nod_pt->is_hanging(i))
       {
        //If it is geometric OR has hanging node data that differs
        //from the geometric data, we must do some work
        if((i==-1) || (nod_pt->hanging_pt(i)!=nod_pt->hanging_pt()))
         {
          // Find out the ultimate map of dependencies: Master nodes
          // and associated weights
          Vector<Node*> master_nodes;
          Vector<double> hanging_weights;
          complete_hanging_nodes_recursively(nod_pt,
                                             master_nodes,hanging_weights,i);
          
          // put them into a map to merge all the occurences of the same node 
          // (add the weights)
          std::map<Node*,double> hang_weights;   
          unsigned n_master=master_nodes.size();
          for(unsigned k=0;k<n_master;k++)
           {
            if(std::fabs(hanging_weights[k])>min_weight)
             hang_weights[master_nodes[k]]+=hanging_weights[k];
           }
          
          //Create new hanging data (we know how many data there are)
          HangInfo* hang_pt = new HangInfo(hang_weights.size());
          
          unsigned hang_weights_index=0;
          //Copy the map into the HangInfo object
          typedef std::map<Node*,double>::iterator IT;
          for (IT it=hang_weights.begin();it!=hang_weights.end();++it)
           {
            hang_pt->set_master_node_pt(hang_weights_index,it->first,
                                        it->second);
            ++hang_weights_index;
           }
          
          //Assign the new hanging pointer to the appropriate value
          nod_pt->set_hanging_pt(hang_pt,i);
         }
       }
     }
   }

#ifdef PARANOID

  // Check hanging node scheme: The weights need to add up to one
  //-------------------------------------------------------------
  //Loop over all values indices
  for (int i=-1;i<ncont_interpolated_values;i++)
   {
    //Loop over all nodes in mesh
    for (unsigned long n=0;n<n_node;n++)
     {
      //Set a local pointer to the node
      Node* nod_pt=this->node_pt(n);
      
      // Is it hanging?
      if (nod_pt->is_hanging(i))
       {
        unsigned nmaster= nod_pt->hanging_pt(i)->nmaster();
        double sum=0.0;
        for (unsigned imaster=0;imaster<nmaster;imaster++)
         {
          sum+=nod_pt->hanging_pt(i)->master_weight(imaster);
         }
        if (std::fabs(sum-1.0)>1.0e-7)
         {
          oomph_info << "WARNING: Sum of master node weights fabs(sum-1.0) " 
               << std::fabs(sum-1.0) << " for node number " << n 
               << " at value " << i << std::endl;
         }
       }
     }
   }
#endif
  
 }

// Sorting out the cases of nodes that are hanging on at least one but not
// all of the processors for which they are part of the local mesh.

#ifdef OOMPH_HAS_MPI

//========================================================================
/// Deal with nodes that are hanging on one process but not another
/// (i.e. the hanging status of the haloed and halo layers disagrees)
//========================================================================
void TreeBasedRefineableMeshBase::synchronise_hanging_nodes
(const unsigned& ncont_interpolated_values)
{
 // Store number of processors and current process
 MPI_Status status;
 int n_proc=Comm_pt->nproc();
 int my_rank=Comm_pt->my_rank();

 double t_start = 0.0;
 double t_end = 0.0;
 
 // Storage for the hanging status of halo/haloed nodes on elements
 Vector<Vector<int> > haloed_hanging(n_proc);
 Vector<Vector<int> > halo_hanging(n_proc);

 // Counter for the number of nodes which require additional synchronisation
 unsigned nnode_still_requiring_synchronisation=0;

 if (Global_timings::Doc_comprehensive_timings)
  {
   t_start = TimingHelpers::timer();
  }

 // Store number of continuosly interpolated values as int
 int ncont_inter_values=ncont_interpolated_values;
 
 // Loop over processes: Each processor checks that is haloed nodes
 // with proc d have consistent hanging stats with halo counterparts.
 for (int d=0; d<n_proc; d++)
  {

   // No halo with self: Setup hang info for my haloed nodes with proc d 
   // then get ready to receive halo info from processor d.
   if (d!=my_rank) 
    {
     
     // Loop over haloed nodes
     unsigned nh=nhaloed_node(d);     
     for (unsigned j=0;j<nh;j++)
      {
       // Get node
       Node* nod_pt=haloed_node_pt(d,j);
       
       // Loop over the hanging status for each interpolated variable
       // (and the geometry)
       for (int icont=-1; icont<ncont_inter_values; icont++)
        { 
         // Store the hanging status of this haloed node
         if (nod_pt->is_hanging(icont))
          {
           unsigned n_master=nod_pt->hanging_pt(icont)->nmaster();
           haloed_hanging[d].push_back(n_master);
          }
         else
          {
           haloed_hanging[d].push_back(0);
          }
        }
      }
     
     // Receive the hanging status information from the corresponding process
     unsigned count_haloed=haloed_hanging[d].size();
     
#ifdef PARANOID
     // Check that number of halo and haloed data match
     unsigned tmp=0;     
     MPI_Recv(&tmp,1,MPI_UNSIGNED,d,0,Comm_pt->mpi_comm(),&status);
     if (tmp!=count_haloed)
      {
       std::ostringstream error_stream;
       error_stream  << "Number of halo data, " << tmp  
                     << ", does not match number of haloed data, " 
                     << count_haloed << std::endl;
       throw OomphLibError(
        error_stream.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
#endif
     
     // Get the data (if any)
     if (count_haloed!=0)
      {
       halo_hanging[d].resize(count_haloed);
       MPI_Recv(&halo_hanging[d][0],count_haloed,MPI_INT,d,0,
                Comm_pt->mpi_comm(),&status);       
      }
    }
   else // d==my_rank, i.e. current process: Send halo hanging status 
        // to process dd where it's received (see above) and compared
        // and compared against the hang status of the haloed nodes
    {
     for (int dd=0; dd<n_proc; dd++)
      {
       // No halo with yourself
       if (dd!=d)
        {
       
         // Storage for halo hanging status and counter
         Vector<int> local_halo_hanging;

         // Loop over halo nodes
         unsigned nh=nhalo_node(dd);
         for (unsigned j=0;j<nh;j++)
          {
           // Get node
           Node* nod_pt=halo_node_pt(dd,j);
           
           // Loop over the hanging status for each interpolated variable
           // (and the geometry)
           for (int icont=-1; icont<ncont_inter_values; icont++)
            { 
             // Store hanging status of halo node
             if (nod_pt->is_hanging(icont))
              {
               unsigned n_master=nod_pt->hanging_pt(icont)->nmaster();
               local_halo_hanging.push_back(n_master);
              }
             else
              {
               local_halo_hanging.push_back(0);
              }
            }
          }
         
         
         // Send the information to the relevant process
         unsigned count_halo=local_halo_hanging.size();
         
#ifdef PARANOID
         // Check that number of halo and haloed data match
         MPI_Send(&count_halo,1,MPI_UNSIGNED,dd,0,Comm_pt->mpi_comm());
#endif
         
         // Send data (if any)
         if (count_halo!=0)
          {
           MPI_Send(&local_halo_hanging[0],count_halo,MPI_INT,
                    dd,0,Comm_pt->mpi_comm());
          }
        }
      }
    }
  }

 if (Global_timings::Doc_comprehensive_timings)
  {
   t_end = TimingHelpers::timer();
   oomph_info << "Time for first all-to-all in synchronise_hanging_nodes(): " 
              << t_end-t_start << std::endl;
   t_start = TimingHelpers::timer();
  }
 
  
 // Now compare equivalent halo and haloed vectors to find discrepancies.
 // It is possible that a master node may not be on either process involved
 // in the halo-haloed scheme; to work round this, we use the shared_node
 // storage scheme, which stores all nodes that are on each pair of processors
 // in the same order on each of the two processors

 // Vector to store info that will help us locate master nodes on "other" 
 // processor
 Vector<HangHelperStruct> hang_info;

 // Copy vector-based representation of shared nodes into
 // map for faster search
 Vector<std::map<Node*,unsigned> > shared_node_map(n_proc);
 for (int d=0;d<n_proc;d++)
  {
   unsigned n=Shared_node_pt[d].size();
   for (unsigned jj=0;jj<n;jj++)
    {
     shared_node_map[d][Shared_node_pt[d][jj]]=jj;
    }
  }


 // Loop over domains: Each processor checks consistency of hang status
 // of its haloed nodes with proc d against the halo counterpart. Haloed
 // wins if there are any discrepancies.
 for (int d=0; d<n_proc; d++)
  {
   // No halo with yourself
   if (d!=my_rank)      
    {
 
     unsigned discrepancy_count=0;
     unsigned discrepancy_count_buff=0;
     
     // Storage for hanging information that needs to be sent to the
     // relevant process if there is a discrepancy in the hanging status
     Vector<int> send_data;
     Vector<double> send_double_data;
     //Buffer storage for data to be sent
     //(We need this because we cannot tell until the end of the loop over
     //a node's master nodes whether we can reconcile its hanging status)
     Vector<int> send_data_buff(0);
     Vector<double> send_double_data_buff(0);
     
     // Counter for traversing haloed data
     unsigned count=0;
     
     // Indicate presence of discrepancy. Default: there's none
     unsigned discrepancy=0;
     unsigned discrepancy_buff=0;
     
     // Loop over haloed nodes
     unsigned nh=nhaloed_node(d);
     for (unsigned j=0;j<nh;j++)
      {
       // Get node
       Node* nod_pt=haloed_node_pt(d,j);
       
       // Loop over the hanging status for each interpolated variable
       // (and the geometry)
       for (int icont=-1; icont<ncont_inter_values; icont++)
        { 
         // Compare hanging status of halo/haloed counterpart structure
         
         // Haloed is is hanging and haloed has different number
         // of master nodes (which includes none in which case it isn't
         // hanging)
         if ((haloed_hanging[d][count]>0)&&
             (haloed_hanging[d][count]!=halo_hanging[d][count]))
          {
           discrepancy_buff=1;
           discrepancy_count_buff++;

           //Flag to check if all masters of this node have been found
           bool found_all_masters=true;

           // Find master nodes of haloed node
           HangInfo* hang_pt=nod_pt->hanging_pt(icont);
           unsigned nhd_master=hang_pt->nmaster();
           
           // Add the number of master nodes to the hanging_nodes vector
           send_data_buff.push_back(nhd_master);
           
           // Loop over master nodes required for HangInfo
           for (unsigned m=0; m<nhd_master; m++)
            {
             // Get mth master node
             Node* master_nod_pt=hang_pt->master_node_pt(m);
             


//              //------------------------------------------
//              // Old direct search is much slower!
//              // Keep this code alive to allow comparisons
//
//              double t_start_search=TimingHelpers::timer();
//
//              // This node will be shared: find it!
//              bool found=false;
//
//              // Look in shared node storage with domain d first
//              unsigned nnod_shared=nshared_node(d);
//              for (unsigned k=0; k<nnod_shared; k++)
//               {
//                if (master_nod_pt==shared_node_pt(d,k))
//                 {
//                  // Found a master: Put its number in the shared 
//                  // scheme into send data
//                  send_data.push_back(k);
//
//                  // Add the weight
//                  send_double_data.push_back(hang_pt->master_weight(m));
//                 
//                  // Done
//                  found=true;
//                  break;
//                 }
//               }
//
//              double t_end_search=TimingHelpers::timer();
//              t_start_search_total+=(t_end_search-t_start_search);
//
//              // end direct search demo
//              //-----------------------

             
             // This node will be shared: find it!
             bool found=false;

             // Which processor holds the non-halo counterpart of this
             // node?
             int non_halo_proc_id=master_nod_pt->non_halo_proc_ID();

             // Try to find node in map with proc d and get iterator to entry
             std::map<Node*,unsigned>::iterator it=
              shared_node_map[d].find(master_nod_pt); 
             
             // If it's not in there iterator points to end of
             // set
             if (it!=shared_node_map[d].end())
              {
               // Found a master: When looking up the node in the shared
               // node scheme on processor d, which processor do I work
               // with? The current one
               send_data_buff.push_back(my_rank);

               // Found a master: Send its index in the shared 
               // node scheme
               send_data_buff.push_back((*it).second);
               
               // Add the weight
               send_double_data_buff.push_back(hang_pt->master_weight(m));
               
               // Done
               found=true;
              }

             // If we haven't found it in the shared node scheme with proc d
             // find it in the shared node scheme with the processor that holds
             // the non-halo version
             if (!found)
              {

               // This is odd -- can't currently handle the case where
               // node is owned by current processor (indicated by 
               // non_halo_proc_id being negative
               if (non_halo_proc_id<0)
                {
                 // This case is now handled by the function
                 // additional_synchronise_hanging_nodes()
                 // called (if necessary) at the end
                 //OomphLibWarning(
                 // "Odd: missing master node is owned by current proc. Will crash below.",
                 // "TreeBasedRefineableMeshBase::synchronise_hanging_nodes(...)",
                 // OOMPH_EXCEPTION_LOCATION);
                }
               else // i.e. (non_halo_proc_id>=0)
                {
                 if (shared_node_map[non_halo_proc_id].size()>0)
                  {
                   std::map<Node*,unsigned>::iterator it=
                    shared_node_map[non_halo_proc_id].find(master_nod_pt); 
                   
                   // If it's not in there iterator points to end of
                   // set
                   if (it!=shared_node_map[non_halo_proc_id].end())
                    {
                     // Found a master: Send ID of processor that holds 
                     // non-halo (the fact that this is different from 
                     // my_rank (the processor that sends this) will alert 
                     // the other processor to the fact that it needs to 
                     send_data_buff.push_back(non_halo_proc_id);
                     
                     // Found a master: Send its index in the shared 
                     // node scheme
                     send_data_buff.push_back((*it).second);
                     
                     // Add the weight
                     send_double_data_buff.push_back(hang_pt->master_weight(m));
                     
                     // Done
                     found=true;

                     // It is possible that the master node found in the shared
                     // storage with processor <non_halo_proc_id> does not
                     // actually exist on processor d. If it does turn out to
                     // exist, we are ok, but if not then we must remember to
                     // create it later (in the
                     // additional_synchronise_hanging_nodes() function)
                    }
                  }                 
                }
              }

             /*
             // Don't throw error, we will construct missing master nodes in
             // additional_synchronise_hanging_nodes() below
             
             // Paranoid check: if we haven't found the master node
             // then throw an error
             if (!found)
              {
               char filename[100];
               std::ofstream some_file;
               sprintf(filename,"sync_hanging_node_crash_mesh_proc%i.dat",
                       my_rank);
               some_file.open(filename);
               this->output(some_file);
               some_file.close();
               
               sprintf(filename,
                       "sync_hanging_node_crash_mesh_with_haloes_proc%i.dat",
                       my_rank);
               some_file.open(filename);
               this->enable_output_of_halo_elements();
               this->output(some_file);
               this->disable_output_of_halo_elements();
               some_file.close();
               

               std::set<unsigned> other_proc_id;
               other_proc_id.insert(d);
               other_proc_id.insert(non_halo_proc_id);
               for (std::set<unsigned>::iterator it=other_proc_id.begin();
                    it!=other_proc_id.end();it++)
                {
                 unsigned d_doc=(*it);
                 
                 sprintf(
                  filename,
                  "sync_hanging_node_crash_halo_elements_with_proc%i_proc%i.dat",
                  d_doc,my_rank);
                 some_file.open(filename);
                 Vector<GeneralisedElement*> 
                  halo_elem_pt(this->halo_element_pt(d_doc));
                 unsigned nelem=halo_elem_pt.size();
                 for (unsigned e=0;e<nelem;e++)
                  {
                   FiniteElement* el_pt=
                    dynamic_cast<FiniteElement*>(halo_elem_pt[e]);
                   el_pt->output(some_file);
                  }
                 some_file.close();
                 
                 sprintf(
                  filename,
                  "sync_hanging_node_crash_haloed_elements_with_proc%i_proc%i.dat",
                  d_doc,my_rank);
                 some_file.open(filename);
                 Vector<GeneralisedElement*> 
                  haloed_elem_pt(this->haloed_element_pt(d_doc));
                 nelem=haloed_elem_pt.size();
                 for (unsigned e=0;e<nelem;e++)
                  {
                   FiniteElement* el_pt=
                    dynamic_cast<FiniteElement*>(haloed_elem_pt[e]);
                   el_pt->output(some_file);
                  }
                 some_file.close();
                 
                 
                 sprintf(
                  filename,
                  "sync_hanging_node_crash_shared_nodes_with_proc%i_proc%i.dat",
                  d_doc,my_rank);
                 some_file.open(filename);
                 unsigned n=nshared_node(d_doc);
                 for (unsigned j=0;j<n;j++)
                  {
                   Node* nod_pt=shared_node_pt(d_doc,j);
                   unsigned nd=nod_pt->ndim();
                   for (unsigned i=0;i<nd;i++)
                    {
                     some_file << nod_pt->x(i) << " ";
                    }
                   some_file << "\n";
                  }
                 some_file.close();
                 
                 
                 sprintf(
                  filename,
                  "sync_hanging_node_crash_halo_nodes_with_proc%i_proc%i.dat",
                  d_doc,my_rank);
                 some_file.open(filename);
                 n=nhalo_node(d_doc);
                 for (unsigned j=0;j<n;j++)
                  {
                   Node* nod_pt=halo_node_pt(d_doc,j);
                   unsigned nd=nod_pt->ndim();
                   for (unsigned i=0;i<nd;i++)
                    {
                     some_file << nod_pt->x(i) << " ";
                    }
                   some_file << "\n";
                  }
                 some_file.close();
                 
                 
                 sprintf(
                  filename,
                  "sync_hanging_node_crash_haloed_nodes_with_proc%i_proc%i.dat",
                  d_doc,my_rank);
                 some_file.open(filename);
                 n=nhaloed_node(d_doc);
                 for (unsigned j=0;j<n;j++)
                  {
                   Node* nod_pt=haloed_node_pt(d_doc,j);
                   unsigned nd=nod_pt->ndim();
                   for (unsigned i=0;i<nd;i++)
                    {
                     some_file << nod_pt->x(i) << " ";
                    }
                   some_file << "\n";
                  }
                 some_file.close();
                 
                } // end of loop over all inter-processor lookup schemes
               
               std::ostringstream error_stream;
               unsigned n=master_nod_pt->ndim();
               error_stream  << "Error: Master node at:\n\n";
               for (unsigned i=0;i<n;i++)
                {
                 error_stream <<  master_nod_pt->x(i) << " ";
                }
               error_stream   << "\n\nnot found for icont="
                              << icont << "in  " 
                              << "shared node storage with proc " << d << "\n"
                              << "or in shared node storage with proc " 
                              << non_halo_proc_id 
                              << " which is where its non-halo counterpart lives.\n"
                              << "Relevant files: sync_hanging_node_crash*.dat\n\n" 
                              << "Hanging node itself: \n\n";
               n=nod_pt->ndim();
               for (unsigned i=0;i<n;i++)
                {
                 error_stream << nod_pt->x(i) << " ";
                }
               error_stream << nod_pt->non_halo_proc_ID();
               error_stream << "\n\nMaster nodes:\n\n";
               for (unsigned m=0; m<nhd_master; m++)
                {
                 Node* master_nod_pt=hang_pt->master_node_pt(m);
                 n=master_nod_pt->ndim();
                 for (unsigned i=0;i<n;i++)
                  {
                   error_stream << master_nod_pt->x(i) << " ";
                  }
                 error_stream << master_nod_pt->non_halo_proc_ID();
                 error_stream << "\n";
                }
               
               // try to find it somewhere else -- sub-optimal search but
               // we're about to die on our arses (yes, plural -- this is 
               // a parallel run!) anyway...
               for (int dddd=0;dddd<n_proc;dddd++)
                {
                 bool loc_found=false;
                 unsigned nnnod_shared=nshared_node(dddd);
                 for (unsigned k=0; k<nnnod_shared; k++)
                  {
                   if (master_nod_pt==shared_node_pt(dddd,k))
                    {
                     loc_found=true;
                     error_stream 
                      << "Found that master node as " << k 
                      << "-th entry in shared node storage with proc " 
                      << dddd << "\n";
                    }
                  }
                 if (!loc_found)
                  {
                   error_stream 
                    << "Did not find that master node in shared node storage with proc " 
                    << dddd << "\n";
                  }
                }
               error_stream << "\n\n";

               throw OomphLibError(
                error_stream.str(),
                OOMPH_CURRENT_FUNCTION,
                OOMPH_EXCEPTION_LOCATION);
              }
             */

             //Check if the master has been found
             if (!found)
              {
               //If this master hasn't been found then set the flag
               found_all_masters=false;
               //No need to continue searching for masters
               break;
              }
             
             
            } // loop over master nodes
           

           // Check if we need to send the data
           if(found_all_masters)
            {
             // All masters were found, so populate send data from buffer
             discrepancy = discrepancy_buff;
             discrepancy_count += discrepancy_count_buff;
             for(unsigned i=0; i<send_data_buff.size(); i++)
              {
               send_data.push_back(send_data_buff[i]);
              }
             for(unsigned i=0; i<send_double_data_buff.size(); i++)
              {
               send_double_data.push_back(send_double_data_buff[i]);
              }

             // Clear buffers and reset
             discrepancy_buff = 0;
             discrepancy_count_buff = 0;
             send_data_buff.clear();
             send_double_data_buff.clear();
            }
           else
            {
             // At least one master node was not found, so we can't
             // reconcile the hanging status of this node yet. We tell
             // the other processor to do nothing for now.
             send_data.push_back(0);

             // Clear buffers and reset
             discrepancy_buff = 0;
             discrepancy_count_buff = 0;
             send_data_buff.clear();
             send_double_data_buff.clear();

             // Set flag to trigger another round of synchronisation
             nnode_still_requiring_synchronisation++;
            }
             
          }
         // Haloed node isn't hanging but halo is: the latter
         // shouldn't so send a -1 to indicate that it's to be made
         // non-hanging
         else if ((haloed_hanging[d][count]==0) && 
                  (halo_hanging[d][count]>0))
          {
           discrepancy=1;
           discrepancy_count++;
           send_data.push_back(-1);
          }
         // Both halo and haloed node have the same number of masters
         // we're happy!
         else if (haloed_hanging[d][count]==
                  halo_hanging[d][count]) 
          {
           send_data.push_back(0);
          }
         else
          {
           std::ostringstream error_stream;
           error_stream  << "Never get here!\n " 
                         << "haloed_hanging[d][count]=" << haloed_hanging[d][count]
                         << "; halo_hanging[d][count]=" <<   halo_hanging[d][count]
                         << std::endl;
           throw OomphLibError(
            error_stream.str(),
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
          }
         // Increment counter for number of haloed data
         count++;
        } // end of loop over icont
      } // end of loop over haloed nodes
     
     // Now send all the required info to the equivalent halo layer -
     // If there are no discrepancies, no need to send anything
     Vector<unsigned> n_all_send(2,0);
     if (discrepancy==0)
      {
       MPI_Send(&n_all_send[0],2,MPI_UNSIGNED,d,0,Comm_pt->mpi_comm());
      }
     else
      {       
       // How much data is there to be sent?
       n_all_send[0]=send_data.size();         
       n_all_send[1]=send_double_data.size();         
       MPI_Send(&n_all_send[0],2,MPI_UNSIGNED,d,0,Comm_pt->mpi_comm());
       
       // Send flat-packed ints
       if (n_all_send[0]!=0)
        {
         MPI_Send(&send_data[0],n_all_send[0],MPI_INT,d,1,
                  Comm_pt->mpi_comm());
        }
       
       // Send flat-packed double data
       if (n_all_send[1]!=0)
        {
         MPI_Send(&send_double_data[0],n_all_send[1],MPI_DOUBLE,d,1,
                  Comm_pt->mpi_comm());
        }
      }
    }
   else // (d==my_rank), current process
    {
     // Receive the master nodes and weights in order to modify the
     // hanging status of nodes in the halo layer
     for (int dd=0; dd<n_proc; dd++)
      {
       if (dd!=d) // don't talk to yourself
        {
         // Are we going to receive anything? This is zero
         // either if there are no discrepancies or there's zero
         // data to be sent.
         Vector<unsigned> n_all_recv(2,0);
         MPI_Recv(&n_all_recv[0],2,MPI_UNSIGNED,dd,0,
                  Comm_pt->mpi_comm(),&status);

         // Storage for received information
         Vector<int> receive_data;
         Vector<double> receive_double_data;
         
         // Receive unsigneds (if any)
         if (n_all_recv[0]!=0)
          {
           // Receive the data
           receive_data.resize(n_all_recv[0]);
           MPI_Recv(&receive_data[0],n_all_recv[0],MPI_INT,dd,1,
                    Comm_pt->mpi_comm(),&status);
          }

         // Receive doubles (if any)
         if (n_all_recv[1]!=0)
          {
           // Receive the data
           receive_double_data.resize(n_all_recv[1]);
           MPI_Recv(&receive_double_data[0],n_all_recv[1],MPI_DOUBLE,dd,1,
                    Comm_pt->mpi_comm(),&status);
          }
         
         // If no information, no need to do anything else  
         if (n_all_recv[0]!=0)
          {
           // Counters for traversing received data
           unsigned count=0;
           unsigned count_double=0;
           
           // Loop over halo nodes
           unsigned nh=nhalo_node(dd);
           for (unsigned j=0;j<nh;j++)
            {
             // Get node
             Node* nod_pt=halo_node_pt(dd,j);           
             
             // Loop over the hanging status for each interpolated variable
             // (and the geometry)
             for (int icont=-1; icont<ncont_inter_values; icont++)
              { 
               
               // Read next entry
               int next_entry=receive_data[count++];
               
               // If it's positive, then the number tells us how 
               // many master nodes we have
               if (next_entry>0) 
                {
                 unsigned nhd_master=unsigned(next_entry);
                 
                 // Set up a new HangInfo for this node
                 HangInfo* hang_pt = new HangInfo(nhd_master);
                 
                 // Now set up the master nodes and weights
                 for (unsigned m=0; m<nhd_master; m++)
                  {
                   // Get the sent master node (a shared node) and 
                   // the weight

                   // ID of proc in whose shared node lookup scheme
                   // the sending processor found the node
                   unsigned shared_node_proc=unsigned(receive_data[count++]);

                   // Index of node in the shared node lookup scheme 
                   unsigned shared_node_id=unsigned(receive_data[count++]);
                   
                   // Get weight
                   double mtr_weight=receive_double_data[count_double++];

                   // If the shared node processor is the same as the
                   // the sending processor we can processor everything here
                   if (shared_node_proc==unsigned(dd))
                    {
                     // Get node
                     Node* master_nod_pt=shared_node_pt(dd,shared_node_id);
                     
                     // Set as a master node (with corresponding weight)
                     hang_pt->set_master_node_pt(m,master_nod_pt,mtr_weight);
                    }
                   // ...otherwise we have do another communication with
                   // intermediate processor that holds the non-halo
                   // version of the master node -- only that processor can
                   // translate the index of the node the share node 
                   // lookup scheme with the sending processor to the
                   // index in the shared node lookup scheme with this
                   // processor
                   else
                    {
                     // Store
                     HangHelperStruct tmp;
                     tmp.Sending_processor=dd;
                     tmp.Shared_node_id_on_sending_processor=shared_node_id;
                     tmp.Shared_node_proc=shared_node_proc;
                     tmp.Weight=mtr_weight;
                     tmp.Hang_pt=hang_pt;
                     tmp.Master_node_index=m;
                     tmp.Node_pt=nod_pt;
                     tmp.icont=icont;
                     hang_info.push_back(tmp);
                    }
                  }

                 // Set the hanging pointer for the current halo node
                 // (does delete any previous hang data)
                 nod_pt->set_hanging_pt(hang_pt,icont);
                }
               // Negative entry: the hanging node already exists, 
               // but it shouldn't, so set it to nonhanging
               else if (next_entry<0)
                {
                 nod_pt->set_hanging_pt(0,icont);
                }
               
              } // end of loop over icont
            } // end of loop over nodes
          } // end of anything to receive
        }
      }
    }
  }// end loop over all processors
 


 if (Global_timings::Doc_comprehensive_timings)
  {
   t_end = TimingHelpers::timer();
   oomph_info << "Time for second all-to-all in synchronise_hanging_nodes() " 
              << t_end-t_start << std::endl;
   t_start = TimingHelpers::timer();
  }  
 

 // Now identify master nodes by translating index in shared 
 // node lookup scheme from the lookup scheme with the sending
 // processor to that with the current processor
 unsigned n=hang_info.size();

 // Is there anything to do be done?
 unsigned global_n=0;
 MPI_Allreduce(&n,&global_n,1,MPI_UNSIGNED,MPI_MAX,
               Comm_pt->mpi_comm());
 if (global_n==0)
  {
   oomph_info 
    << "No need for reconciliation of wrongly synchronised hang nodes\n";
  }
 // Reconcilation required
 else  
 {

   oomph_info 
    << "Need to reconcile of wrongly syncronised hang nodes\n";

   // Storage for the translated entries (in order) for/from
   // other processors
   // Must be ints so that an entry of -1 tells the other processor
   // that the node could not be found. See comment above for why
   // this may be necessary.
   Vector<Vector<int> > translated_entry(n_proc);
   
   // Storage for how-many-th entry in this processor's 
   // hang_info vector will be completed by processor rank.
   Vector<Vector<unsigned> > hang_info_index_for_proc(n_proc);
   
   // Send information to intermediate processor that holds
   // non-halo version of missing master node
   {
    // Storage for number of data to be sent to each processor
    Vector<int> send_n(n_proc,0);
    
    // Storage for all values to be sent to all processors
    Vector<int> send_data;
    
    // Start location within send_data for data to be sent to each processor 
    Vector<int> send_displacement(n_proc,0);
    
    // Loop over all processors
    for(int rank=0;rank<n_proc;rank++)
     {  
      //Set the offset for the current processor
      send_displacement[rank] = send_data.size();
      
      //Don't bother to do anything if the processor in the loop is the 
      //current processor
      if(rank!=my_rank)
       {
        
        // Search through the (typically few) entries
        // in hang info vector to find the ones that
        // must be sent to proc rank (the proc that holds
        // the non-halo version of the "missing" master node
        for (unsigned i=0;i<n;i++)
         {
          HangHelperStruct tmp=hang_info[i];
          if (tmp.Shared_node_proc==unsigned(rank))
           {
            // Add the sending processor
            send_data.push_back(tmp.Sending_processor);
            
            // Add the index of the missing master node
            // in the shared node lookup scheme between
            // sending processor and processor rank
            send_data.push_back(tmp.Shared_node_id_on_sending_processor);
            
            // Record the how-many-th entry in this processor's 
            // hang_info vector will be completed by processor rank.
            hang_info_index_for_proc[rank].push_back(i);
           }
         }
       }
      
      //Find the number of data added to the vector
      send_n[rank] = send_data.size() - send_displacement[rank];
     }
    
    
    //Storage for the number of data to be received from each processor
    Vector<int> receive_n(n_proc,0);
    
    //Now send numbers of data to be sent between all processors
    MPI_Alltoall(&send_n[0],1,MPI_INT,&receive_n[0],1,MPI_INT,
                 Comm_pt->mpi_comm());
    
    //We now prepare the data to be received
    //by working out the displacements from the received data
    Vector<int> receive_displacement(n_proc,0);
    int receive_data_count=0;
    for(int rank=0;rank<n_proc;++rank)
     {
      //Displacement is number of data received so far
      receive_displacement[rank] = receive_data_count;
      receive_data_count += receive_n[rank];
     }
    
    //Now resize the receive buffer for all data from all processors
    //Make sure that it has a size of at least one
    if(receive_data_count==0) {++receive_data_count;}
    Vector<unsigned> receive_data(receive_data_count);
    
    //Make sure that the send buffer has size at least one
    //so that we don't get a segmentation fault
    if(send_data.size()==0) {send_data.resize(1);}
    
    //Now send the data between all the processors
    MPI_Alltoallv(&send_data[0],&send_n[0],&send_displacement[0],
                  MPI_INT,
                  &receive_data[0],&receive_n[0],
                  &receive_displacement[0],
                  MPI_INT,
                  Comm_pt->mpi_comm());
    
    //Now use the received data to update the halo nodes
    for (int send_rank=0;send_rank<n_proc;send_rank++)
     {
      //Don't bother to do anything for the processor corresponding to the
      //current processor or if no data were received from this processor
      if((send_rank != my_rank) && (receive_n[send_rank] != 0))
       {
        //Counter for the data within the large array
        unsigned count=receive_displacement[send_rank];
        
        // We're reading two numbers per missing halo node
        unsigned n_rec=unsigned(receive_n[send_rank]);
        for (unsigned i=0;i<n_rec/2;i++)
         {
          
          // Receive orig sending proc
          unsigned orig_sending_proc=receive_data[count];
          count++;
          
          // Receive the index of the missing master node
          // in the shared node lookup scheme between
          // orig sending processor and current processor
          unsigned shared_node_id_on_orig_sending_proc=
           receive_data[count];
          count++;
          
          // Extract node from shared node lookup scheme
          Node* master_nod_pt=
           shared_node_pt(orig_sending_proc,
                          shared_node_id_on_orig_sending_proc);
          
          // Now find it in shared halo scheme with the processor
          // that's sent the request
          
          std::map<Node*,unsigned>::iterator it=
           shared_node_map[send_rank].find(master_nod_pt); 
          
          // If it's not in there iterator points to end of
          // set
          if (it!=shared_node_map[send_rank].end())
           {          
            // Store it so we can send it back
            translated_entry[send_rank].push_back((*it).second);
           }
          else
           {
            // This node has not been found in the shared scheme, so
            // the translation query has failed. We send a -1 to tell
            // the other processor the bad news.
            translated_entry[send_rank].push_back(-1);

            /*
            // We don't need to crash anymore because the function
            // additional_synchronise_hanging_nodes() will magically
            // sort out all the problems!
            std::ostringstream error_stream;
            error_stream
             << "Received translation query for shared node"
             << " entry " << shared_node_id_on_orig_sending_proc 
             << " with processor " << orig_sending_proc 
             << " from proc " << send_rank << std::endl
             << "but did not find node in shared node scheme with proc "
             << send_rank << std::endl;
            throw OomphLibError(
             error_stream.str(),
             OOMPH_CURRENT_FUNCTION,
             OOMPH_EXCEPTION_LOCATION);
            */
           }
          
         }
        
       }
     } //End of data is received
    
   }   // end of sending stuff to intermediate processor that holds
       // non halo version of missing master node
   
   
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for third all-to-all in synchronise_hanging_nodes() " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }  
   
   
   
   // Send information back to processor that needs to identify
   // missing master node via shared node lookup scheme with
   // this processor
   {
    // Storage for number of data to be sent to each processor
    Vector<int> send_n(n_proc,0);
    
    // Storage for all values to be sent to all processors
    Vector<int> send_data;
    
    // Start location within send_data for data to be sent to each processor 
    Vector<int> send_displacement(n_proc,0);
    
    // Loop over all processors
    for(int rank=0;rank<n_proc;rank++)
     {  
      //Set the offset for the current processor
      send_displacement[rank] = send_data.size();
      
      //Don't bother to do anything if the processor in the loop is the 
      //current processor
      if(rank!=my_rank)
       {
        // Put the translated entries for processor rank into
        // send data
        unsigned n=translated_entry[rank].size();
        for (unsigned j=0;j<n;j++)
         {
          send_data.push_back(translated_entry[rank][j]);
         }
       }
      
      //Find the number of data added to the vector
      send_n[rank] = send_data.size() - send_displacement[rank];
     }
    
    
    //Storage for the number of data to be received from each processor
    Vector<int> receive_n(n_proc,0);
    
    //Now send numbers of data to be sent between all processors
    MPI_Alltoall(&send_n[0],1,MPI_INT,&receive_n[0],1,MPI_INT,
                 Comm_pt->mpi_comm());
    
    //We now prepare the data to be received
    //by working out the displacements from the received data
    Vector<int> receive_displacement(n_proc,0);
    int receive_data_count=0;
    for(int rank=0;rank<n_proc;++rank)
     {
      //Displacement is number of data received so far
      receive_displacement[rank] = receive_data_count;
      receive_data_count += receive_n[rank];
     }
    
    //Now resize the receive buffer for all data from all processors
    //Make sure that it has a size of at least one
    if(receive_data_count==0) {++receive_data_count;}
    Vector<unsigned> receive_data(receive_data_count);
    
    //Make sure that the send buffer has size at least one
    //so that we don't get a segmentation fault
    if(send_data.size()==0) {send_data.resize(1);}
    
    //Now send the data between all the processors
    MPI_Alltoallv(&send_data[0],&send_n[0],&send_displacement[0],
                  MPI_INT,
                  &receive_data[0],&receive_n[0],
                  &receive_displacement[0],
                  MPI_INT,
                  Comm_pt->mpi_comm());
    
    //Now use the received data to update the halo nodes
    for (int send_rank=0;send_rank<n_proc;send_rank++)
     {
      //Don't bother to do anything for the processor corresponding to the
      //current processor or if no data were received from this processor
      if((send_rank != my_rank) && (receive_n[send_rank] != 0))
       {
        //Counter for the data within the large array
        unsigned count=receive_displacement[send_rank];
        
        // We're reading one number per missing halo node
        unsigned n_rec=unsigned(receive_n[send_rank]);
        for (unsigned i=0;i<n_rec;i++)
         {
          // Index of missing master node in shared node lookup scheme
          // with processor send_rank:
          // Must be an int because failure returns -1
          int index=receive_data[count];
          count++;

          // Translation query has been successful if index >= 0
          if(index >= 0)
           {
            // Recall information associated with that missing master
            unsigned hang_info_index=hang_info_index_for_proc[send_rank][i];
            HangHelperStruct tmp=hang_info[hang_info_index];

            // Extract node from shared node lookup scheme
            Node* master_nod_pt=
             shared_node_pt(send_rank,index);

            // Set as a master node (with corresponding weight)
            tmp.Hang_pt->set_master_node_pt(tmp.Master_node_index,
                                            master_nod_pt,tmp.Weight);
           }
          else
           {
            // Translation query has failed. This is the processor
            // on which the node was a halo, so we must delete the
            // partial hang info.

            // Recall information associated with that missing master
            unsigned hang_info_index=hang_info_index_for_proc[send_rank][i];
            HangHelperStruct tmp=hang_info[hang_info_index];

            // Delete partial hanging information
            tmp.Node_pt->set_hanging_pt(0,tmp.icont);
            
            // Set flag to trigger another round of synchronisation
            // This works even though we don't own the node that
            // still requires synchrionisation because this variable
            // is reduced over all processors at the end
            nnode_still_requiring_synchronisation++;
           }
         }
        
       }
     } //End of data is received
    
   } // end of completing hang info for missing master nodes
   
   
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info 
      << "Time for fourth all-to-all in synchronise_hanging_nodes() " 
      << t_end-t_start << std::endl;
    }  
  } // end  of reconciliation required


 //Get global number of nodes still requiring synchronisation due to
 //missing master nodes
 //This will only be necessary for meshes involving elements
 //with nonuniformly spaced nodes. All other cases will continue to
 //work as before because all nodes will have been successfully
 //synchronised by now
 unsigned global_nnode_still_requiring_synchronisation=0;
 MPI_Allreduce(&nnode_still_requiring_synchronisation,
               &global_nnode_still_requiring_synchronisation,
               1,MPI_UNSIGNED,MPI_MAX,Comm_pt->mpi_comm());
 if (global_nnode_still_requiring_synchronisation>0)
  {

   double tt_start = 0.0;
   if (Global_timings::Doc_comprehensive_timings)
    {
     tt_start=TimingHelpers::timer();
    }

   oomph_info << "Need to do additional synchronisation of hanging nodes"
              << std::endl;

   // Do additional synchronisation
   additional_synchronise_hanging_nodes(ncont_interpolated_values);

   double tt_end=0.0;
   if (Global_timings::Doc_comprehensive_timings)
    {
     tt_end=TimingHelpers::timer();
     oomph_info 
      << "Time for RefineableMesh::additional_synchronise_hanging_nodes() "
      << "in TreeBasedRefineableMeshBase::synchronise_hanging_nodes(): " 
      << tt_end-tt_start << std::endl;
     tt_start = TimingHelpers::timer();
    }

  }
 else
  {
   oomph_info << "No need to do additional synchronisation of hanging nodes"
              << std::endl;
  }
   
}

//========================================================================
/// Synchronise the positions of non-hanging nodes that depend on
/// non-existent neighbours (e.g. h-refinement of neighbouring elements
/// with different p-orders where the shared edge is on the outer edge of
/// the halo layer)
//========================================================================
void TreeBasedRefineableMeshBase::synchronise_nonhanging_nodes()
{
 // Store number of processors and current process
 MPI_Status status;
 int n_proc=Comm_pt->nproc();
 int my_rank=Comm_pt->my_rank();

 double t_start = 0.0;
 double t_end = 0.0;
 
 // Storage for the hanging status of halo/haloed nodes on elements
 Vector<Vector<unsigned> > recv_unsigneds(n_proc);
 Vector<Vector<double> > recv_doubles(n_proc);

 if (Global_timings::Doc_comprehensive_timings)
  {
   t_start = TimingHelpers::timer();
  }

 // Loop over processes: Each processor checks if its nonhaning nodes in
 // haloed elements with proc d require additional information to determine
 // their positions.
 for (int d=0; d<n_proc; d++)
  {

   // No halo with self: Setup hang info for my haloed nodes with proc d 
   // then get ready to receive halo info from processor d.
   if (d!=my_rank) 
    {
     
     // Receive the position information from the corresponding process
     unsigned recv_unsigneds_count=0;
     MPI_Recv(&recv_unsigneds_count,1,MPI_UNSIGNED,d,0,Comm_pt->mpi_comm(),&status);
     unsigned recv_doubles_count=0;
     MPI_Recv(&recv_doubles_count,1,MPI_UNSIGNED,d,1,Comm_pt->mpi_comm(),&status); 

     // Get the data (if any)
     if (recv_unsigneds_count!=0)
      {
       recv_unsigneds[d].resize(recv_unsigneds_count);
       MPI_Recv(&recv_unsigneds[d][0],recv_unsigneds_count,MPI_UNSIGNED,d,0,
                Comm_pt->mpi_comm(),&status);
      }
     if (recv_doubles_count!=0)
      {
       recv_doubles[d].resize(recv_doubles_count);
       MPI_Recv(&recv_doubles[d][0],recv_doubles_count,MPI_DOUBLE,d,1,
                Comm_pt->mpi_comm(),&status);      
      }

     // Counters for received data
     unsigned recv_unsigneds_index = 0;
     double recv_doubles_index = 0;

     // Get halo elements with processor d
     Vector<GeneralisedElement*> halo_element_pt(this->halo_element_pt(d));

     // Loop over recieved indices
     while(recv_unsigneds_index<recv_unsigneds_count)
      {
       // Get (finite) element
       FiniteElement* el_pt=
        dynamic_cast<FiniteElement*>(
         halo_element_pt[recv_unsigneds[d][recv_unsigneds_index++]]);

       // If we have a finite element...
       if(el_pt!=0)
        {
         // Get dimension
         unsigned n_dim = el_pt->dim();

         // Get node
         Node* nod_pt=el_pt->node_pt(recv_unsigneds[d][recv_unsigneds_index++]);
         
         // Get current position
         Vector<double> x_cur(n_dim);
         for(unsigned dir=0; dir<n_dim; dir++)
          {
           x_cur[dir] = nod_pt->x(dir);
          }
         
         // Get recieved position
         Vector<double> x_rec(n_dim);
         for(unsigned dir=0; dir<n_dim; dir++)
          {
           x_rec[dir] = recv_doubles[d][recv_doubles_index+dir];
          }

         // Compare actual and expected positions
         bool node_pos_differs=false;
         for(unsigned dir=0; dir<n_dim; dir++)
          {
           node_pos_differs = node_pos_differs
            || (std::fabs(x_cur[dir]-x_rec[dir])>1.0e-14);
          }
         
         // Set the actual position
         Vector<double> x_act(n_dim);
         for(unsigned dir=0; dir<n_dim; dir++)
          {
           nod_pt->x(dir) = recv_doubles[d][recv_doubles_index++];
          }
        }
      }

     if(recv_unsigneds_count!=recv_unsigneds_index)
      {
       std::ostringstream error_stream;
       error_stream << "recv_unsigneds_count != recv_unsigneds_index ( "
                    << recv_unsigneds_count << " != "
                    << recv_unsigneds_index << ")" << std::endl;
       throw OomphLibError(
              error_stream.str(),
              "TreeBasedRefineableMeshBase::synchronise_nonhanging_nodes()",
              OOMPH_EXCEPTION_LOCATION);
      }

    }
   else // d==my_rank, i.e. current process: Send halo hanging status 
        // to process dd where it's received (see above) and compared
        // and compared against the hang status of the haloed nodes
    {
     for (int dd=0; dd<n_proc; dd++)
      {
       // No halo with yourself
       if (dd!=d)
        {
       
         // Storage for halo hanging status and counter
         Vector<int> send_unsigneds;
         Vector<double> send_doubles;

         // Set to store nodes whose position requires adjustment
         std::set<Node*> nodes_requiring_adjustment;

         // Get haloed elements with processor dd
         Vector<GeneralisedElement*> haloed_element_pt(this->haloed_element_pt(dd));
     
         // Loop over haloed elements with processor dd
         unsigned nh=haloed_element_pt.size();
         for (unsigned e=0;e<nh;e++)
          {
           // Get (finite) element
           FiniteElement* el_pt=dynamic_cast<FiniteElement*>(haloed_element_pt[e]);

           // If we have a finite element...
           if(el_pt!=0)
            {
             // Get dimension
             unsigned n_dim = el_pt->dim();

             // Loop over element nodes
             unsigned n_node = el_pt->nnode();
             for (unsigned j=0;j<n_node;j++)
              {
               // Get node
               Node* nod_pt=el_pt->node_pt(j);
         
               // Only do non-hanging nodes
               if (!nod_pt->is_hanging())
                {
                 // Check if node's position is the same as that interpolated using its
                 // local coordinate in the haloed element

                 // Loop over all history values
                 unsigned nt=nod_pt->ntstorage();
                 for(unsigned t=0;t<nt;t++)
                  {
                   // Get expected position
                   Vector<double> s(n_dim), x_exp(n_dim);
                   el_pt->local_coordinate_of_node(j,s);
                   el_pt->get_x(t,s,x_exp);
               
                   // Get actual position
                   Vector<double> x_act(n_dim);
                   for(unsigned dir=0; dir<n_dim; dir++)
                    {
                     x_act[dir] = nod_pt->x(dir);
                    }

                   // Compare actual and expected positions
                   bool node_pos_differs=false;
                   for(unsigned dir=0; dir<n_dim; dir++)
                    {
                     node_pos_differs = node_pos_differs
                      || (std::fabs(x_act[dir]-x_exp[dir])>1.0e-14);
                    }

                   // If the node's actual position differs from its
                   // expected position we need to communicate this
                   // information to processors on which this is a halo node
                   if(node_pos_differs)
                    {
                     // Check that node has not been done already
                     if(nodes_requiring_adjustment.insert(nod_pt).second)
                      {
                       // Send index of haloed element
                       send_unsigneds.push_back(e);
                       // Send index of node in the element
                       send_unsigneds.push_back(j);
                       // Send actual position of node
                       for(unsigned dir=0; dir<n_dim; dir++)
                        {
                         send_doubles.push_back(x_act[dir]);
                        }
                      }
                    }
                  }
                }
              }
            }
          }
         
         // Send the information to the relevant process
         unsigned send_unsigneds_count=send_unsigneds.size();
         unsigned send_doubles_count=send_doubles.size();

         if(send_unsigneds_count>0)
          {
           //exit(1);
          }

         // Tell processor dd how much data to receive
         MPI_Send(&send_unsigneds_count,1,MPI_UNSIGNED,dd,0,Comm_pt->mpi_comm());
         MPI_Send(&send_doubles_count,1,MPI_UNSIGNED,dd,1,Comm_pt->mpi_comm());
         
         // Send data (if any)
         if (send_unsigneds_count!=0)
          {
           MPI_Send(&send_unsigneds[0],send_unsigneds_count,MPI_UNSIGNED,
                    dd,0,Comm_pt->mpi_comm());
          }
         if (send_doubles_count!=0)
          {
           MPI_Send(&send_doubles[0],send_doubles_count,MPI_DOUBLE,
                    dd,1,Comm_pt->mpi_comm());
          }

        }
      }
    }
  }

 if (Global_timings::Doc_comprehensive_timings)
  {
   t_end = TimingHelpers::timer();
   oomph_info << "Time for synchronise_nonhanging_nodes(): " 
              << t_end-t_start << std::endl;
   t_start = TimingHelpers::timer();
  }
 
}

#endif


////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////




//========================================================================
/// Do adaptive p-refinement for mesh.
/// - Pass Vector of error estimates for all elements.
/// - p-refine those whose errors exceeds the threshold
/// - p-unrefine those whose errors is less than
///   threshold.
//========================================================================
void TreeBasedRefineableMeshBase::p_adapt(const Vector<double>& elemental_error)
{
 //Set the refinement tolerance to be the max permissible error
 double refine_tol=this->max_permitted_error();
 
 //Set the unrefinement tolerance to be the min permissible error 
 double unrefine_tol=this->min_permitted_error();
 
 // Setup doc info
 DocInfo local_doc_info;
 if (doc_info_pt()==0) {local_doc_info.disable_doc();}
 else {local_doc_info=this->doc_info();}
 
 
 // Check that the errors make sense
 if (refine_tol<=unrefine_tol)
  {
   std::ostringstream error_stream;
   error_stream << "Refinement tolerance <= Unrefinement tolerance" 
                << refine_tol << " " << unrefine_tol << std::endl
                << "doesn't make sense and will almost certainly crash" 
                << std::endl
                << "this beautiful code!" << std::endl;
   
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
 
 
 //Select elements for refinement and unrefinement
 //==============================================
 // Reset counter for number of elements that would like to be
 // refined further but can't
 this->nrefinement_overruled()=0;

 unsigned n_refine=0;
 unsigned n_unrefine=0;
 // Loop over all elements and mark them according to the error criterion
 unsigned long Nelement=this->nelement();
 for (unsigned long e=0;e<Nelement;e++)
  {
   //(Cast) pointer to the element
   PRefineableElement* el_pt =
    dynamic_cast<PRefineableElement*>(this->element_pt(e));

   //Check that we can p-refine the element
   if(el_pt!=0)
    {
     //Initially element is not to be refined
     el_pt->deselect_for_p_refinement();
     el_pt->deselect_for_p_unrefinement();
      
     //If the element error exceeds the threshold ...
     if(elemental_error[e] > refine_tol)
      {
       // ... and its refinement level is less than the maximum desired level
       //mark is to be refined
       if((el_pt->p_refinement_is_enabled())&&
          (el_pt->p_order() < this->max_p_refinement_level()))
        {
         el_pt->select_for_p_refinement();
         n_refine++;
        }
       // ... otherwise mark it as having been over-ruled
       else
        {
         this->nrefinement_overruled()+=1;
        }
      }
     if(elemental_error[e] < unrefine_tol)
      {
       // ... and its refinement level is more than the minimum desired level
       //mark is to be refined
       //(Also check that we don't unrefine past the initial refinement level)
       if((el_pt->p_refinement_is_enabled())&&
          (el_pt->p_order() > this->min_p_refinement_level())&&
          (el_pt->p_order() > el_pt->initial_p_order()))
        {
         el_pt->select_for_p_unrefinement();
         n_unrefine++;
        }
       // Don't mark as overruled - it's misleading
       //// ... otherwise mark it as having been over-ruled
       //else
       // {
       //  this->nrefinement_overruled()+=1;
       // }
      }
    }//End of check for p-refineability of element
   else
    {
     oomph_info << "p-refinement is not possible for these elements"
                << std::endl;
     //Don't try to p-refine any more elements
     break;
    }
  }
  
 oomph_info 
  << " \n Number of elements to be refined: " << n_refine << std::endl;
 oomph_info << " \n Number of elements whose refinement was overruled: " 
            << this->nrefinement_overruled() << std::endl;
  
 oomph_info << " \n Number of elements to be unrefined : " 
            << n_unrefine << std::endl << std::endl;



 //Now do the actual mesh adaptation
 //---------------------------------

 // Check whether its worth our while
 // Either some elements want to be refined, 
 // or the number that want to be unrefined are greater than the
 // specified tolerance

 // In a parallel job, it is possible that one process may not have
 // any elements to refine, BUT a neighbouring process may refine an
 // element which changes the hanging status of a node that is on
 // both processes (i.e. a halo(ed) node).  To get around this issue, 
 // ALL processes need to call adapt_mesh if ANY refinement is to 
 // take place anywhere.

 unsigned total_n_refine=0;
#ifdef OOMPH_HAS_MPI
 // Sum n_refine across all processors
 if (this->is_mesh_distributed())
  {
   MPI_Allreduce(&n_refine,&total_n_refine,1,MPI_INT,MPI_SUM,
                 Comm_pt->mpi_comm());
  }
 else 
  {
   total_n_refine=n_refine;
  }
#else
 total_n_refine=n_refine;
#endif

 // There may be some issues with unrefinement too, but I have not
 // been able to come up with an example (either in my head or in a
 // particular problem) where anything has arisen.  I can see that
 // there may be an issue if n_unrefine differs across processes so
 // that (total_n_unrefine > max_keep_unrefined()) on some but not 
 // all processes. I haven't seen any examples of this yet so the 
 // following code may or may not work!  (Andy, 06/03/08)

 unsigned total_n_unrefine=0;
#ifdef OOMPH_HAS_MPI
 // Sum n_unrefine across all processors
 if (this->is_mesh_distributed())
  {
   MPI_Allreduce(&n_unrefine,&total_n_unrefine,1,MPI_INT,MPI_SUM,
                 Comm_pt->mpi_comm());
  }
 else
  {
   total_n_unrefine=n_unrefine;
  }
#else
 total_n_unrefine=n_unrefine;
#endif

 oomph_info << "---> " << total_n_refine << " elements to be refined, and "
            << total_n_unrefine << " to be unrefined, in total." << std::endl;

 if ((total_n_refine > 0) || (total_n_unrefine > this->max_keep_unrefined()))
  {

#ifdef PARANOID
#ifdef OOMPH_HAS_MPI
   
   // Sanity check: Each processor checks if the enforced unrefinement of
   // its haloed element is matched by enforced unrefinement of the
   // corresponding halo elements on the other processors.
   if (this->is_mesh_distributed())
    {
     // Store number of processors and current process
     MPI_Status status;
     int n_proc=Comm_pt->nproc();
     int my_rank=Comm_pt->my_rank();
     
     // Loop over all other domains/processors
     for (int d=0;d<n_proc;d++)
      {
       // Don't talk to yourself
       if (d!=my_rank)
        {

         {
          // Get the vector of halo elements whose non-halo counterpart
          // are on processor d
          Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
          
          // Create vector containing (0)1 to indicate that
          // halo element is (not) to be unrefined
          unsigned nhalo=halo_elem_pt.size();
          Vector<int> halo_to_be_unrefined(nhalo,0);
          for (unsigned e=0;e<nhalo;e++)
           {
            if (dynamic_cast<PRefineableElement*>(halo_elem_pt[e])
                ->to_be_p_unrefined())
             {
              halo_to_be_unrefined[e]=1;
             }
           }

          //Trap the case when there are no halo elements
          //so that we don't get a segfault in the MPI send
          if(nhalo > 0)
           {
            // Send it across
            MPI_Send(&halo_to_be_unrefined[0],nhalo,MPI_INT,
                     d,0,Comm_pt->mpi_comm());
           }
         }

         {
          
          // Get the vector of haloed elements on current processor
          Vector<GeneralisedElement*> 
           haloed_elem_pt(this->haloed_element_pt(d));
          
          // Ask processor d to send vector containing (0)1 for 
          // halo element with current processor to be (not)unrefined
          unsigned nhaloed=haloed_elem_pt.size();
          Vector<int> halo_to_be_unrefined(nhaloed);
          //Trap to catch the case that there are no haloed elements
          if(nhaloed > 0)
           {
            MPI_Recv(&halo_to_be_unrefined[0],nhaloed,MPI_INT,d,0,
                     Comm_pt->mpi_comm(),&status);
           }

          // Check it
          for (unsigned e=0;e<nhaloed;e++)
           {
            if ( ( (halo_to_be_unrefined[e]==0)&&
                   (dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])->
                    to_be_p_unrefined()) ) ||
                 ( (halo_to_be_unrefined[e]==1)&&
                   (!dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])->
                    to_be_p_unrefined()) ) )
             {
              std::ostringstream error_message;
              error_message 
               << "Error in refinement: \n"
               << "Haloed element: " << e << " on proc " << my_rank << " \n"
               << "wants to be unrefined whereas its halo counterpart on\n"
               << "proc " << d << " doesn't (or vice versa)...\n"
               << "This is most likely because the error estimator\n"
               << "has not assigned the same errors to halo and haloed\n"
               << "elements -- it ought to!\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
             }
           }
          }        
        }
       
      }



     // Loop over all other domains/processors
     for (int d=0;d<n_proc;d++)
      {
       // Don't talk to yourself
       if (d!=my_rank)
        {

         {
          // Get the vector of halo elements whose non-halo counterpart
          // are on processor d
          Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
          
          // Create vector containing (0)1 to indicate that
          // halo element is (not) to be refined
          unsigned nhalo=halo_elem_pt.size();
          Vector<int> halo_to_be_refined(nhalo,0);
          for (unsigned e=0;e<nhalo;e++)
           {
            if (dynamic_cast<PRefineableElement*>(halo_elem_pt[e])
                ->to_be_p_refined())
             {
              halo_to_be_refined[e]=1;
             }
           }
          
          // Send it across
          if(nhalo > 0)
           {
            MPI_Send(&halo_to_be_refined[0],nhalo,MPI_INT,
                     d,0,Comm_pt->mpi_comm());
           }
         }

         {
          
          // Get the vector of haloed elements on current processor
          Vector<GeneralisedElement*> 
           haloed_elem_pt(this->haloed_element_pt(d));
          
          // Ask processor d to send vector containing (0)1 for 
          // halo element with current processor to be (not)refined
          unsigned nhaloed=haloed_elem_pt.size();
          Vector<int> halo_to_be_refined(nhaloed);
          if(nhaloed > 0)
           {
            MPI_Recv(&halo_to_be_refined[0],nhaloed,MPI_INT,d,0,
                     Comm_pt->mpi_comm(),&status);
           }

          // Check it
          for (unsigned e=0;e<nhaloed;e++)
           {
            if ( ( (halo_to_be_refined[e]==0)&&
                   (dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])->
                    to_be_p_refined()) ) ||
                 ( (halo_to_be_refined[e]==1)&&
                   (!dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])->
                    to_be_p_refined()) ) )
             {
              std::ostringstream error_message;
              error_message 
               << "Error in refinement: \n"
               << "Haloed element: " << e << " on proc " << my_rank << " \n"
               << "wants to be refined whereas its halo counterpart on\n"
               << "proc " << d << " doesn't (or vice versa)...\n"
               << "This is most likely because the error estimator\n"
               << "has not assigned the same errors to halo and haloed\n"
               << "elements -- it ought to!\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
             }
           }
          }        
        }
       
      }
    }   
#endif
#endif
   
   //Perform the actual adaptation
   p_adapt_mesh(local_doc_info);
   
   //The number of refineable elements is still local to each process
   this->Nunrefined=n_unrefine;
   this->Nrefined=n_refine;
  }
 // If not worthwhile, say so but still reorder nodes and kill external storage
 // for consistency in parallel computations
 else
  {

#ifdef OOMPH_HAS_MPI
   // Delete any external element storage - any interaction will still
   // be set up on the fly again, so we need to get rid of old information.
   // This particularly causes problems in multi-domain examples where
   // we decide not to refine one of the meshes
   this->delete_all_external_storage();
#endif

   // Reorder the nodes within the mesh's node vector
   // to establish a standard ordering regardless of the sequence
   // of mesh refinements -- this is required to allow dump/restart
   // on refined meshes
   this->reorder_nodes();

#ifdef OOMPH_HAS_MPI

    // Now (re-)classify halo and haloed nodes and synchronise hanging
    // nodes
    // This is required in cases where delete_all_external_storage()
    // made slave nodes to external halo nodes nonhanging.
    if (this->is_mesh_distributed())
     {
      DocInfo doc_info;
      doc_info.disable_doc();
      classify_halo_and_haloed_nodes(doc_info,doc_info.is_doc_enabled());
     } 

#endif
   
   if (n_refine==0)
    {
     oomph_info 
      << "\n Not enough benefit in adapting mesh. " 
      << std::endl << std::endl;
    }
   this->Nunrefined=0;
   this->Nrefined=0;
  }

}


//================================================================
/// p-adapt mesh, which exists in two representations,
/// namely as:
///  - a FE mesh 
///  - a forest of Oc or QuadTrees 
///
/// p-refinement/derefinement process is documented (in tecplot-able form)
/// if requested. 
///
/// Procedure:
/// - Loop over all elements and do the p-refinement/unrefinement for
///   those who want to be refined. Note: p-refinement builds fully-
///   functional elements.
/// - For all nodes that were hanging on the previous mesh (and are still
///   marked as such), fill in their nodal values (consistent
///   with the current hanging node scheme) to make sure they are fully 
///   functional, should they have become non-hanging during the 
///   mesh-adaptation. Then mark the nodes as non-hanging.
/// - Delete any nodes that have become obsolete.
/// - Mark up hanging nodes and setup hanging node scheme (incl.
///   recursive cleanup for hanging nodes that depend on other
///   hanging nodes).
/// - run a quick self-test on the neighbour finding scheme and
///   check the integrity of the elements (if PARANOID)
/// - doc hanging node status, boundary conditions, neighbour
///   scheme if requested.
///
///
/// After adaptation, all nodes (whether new or old) have up-to-date
/// current and previous values. 
///
/// If refinement process is being documented, the following information
/// is documented:
/// - The files
///   - "neighbours.dat" 
///   - "all_nodes.dat" 
///   - "new_nodes.dat" 
///   - "hang_nodes_*.dat" 
///     where the * denotes a direction (n,s,e,w) in 2D 
///     or (r,l,u,d,f,b) in 3D
///. 
///   can be viewed with 
///   - QHangingNodes.mcr
///   .
/// - The file
///    - "hangnodes_withmasters.dat"
///    .
///    can be viewed with 
///    - QHangingNodesWithMasters.mcr
///    .
///    to check the hanging node status.
/// - The neighbour status of the elements is documented in
///   - "neighbours.dat"
///   .
///   and can be viewed with 
///   - QuadTreeNeighbours.mcr
///   .
//=================================================================
void TreeBasedRefineableMeshBase::p_adapt_mesh(DocInfo& doc_info)
{
   
#ifdef OOMPH_HAS_MPI
 // Delete any external element storage before performing the adaptation
 // (in particular, external halo nodes that are on mesh boundaries)
 this->delete_all_external_storage();
#endif
  
 //Only perform the adapt step if the mesh has any elements.  This is relevant
 //in a distributed problem with multiple meshes, where a particular
 // process may not have any elements on a particular submesh.
 if (this->nelement()>0)
  {
   double t_start = 0.0;
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_start=TimingHelpers::timer();
    }
    
   // Do refinement/unrefinement if required
   this->p_refine_elements_if_required();

   if (Global_timings::Doc_comprehensive_timings)
    {
     double t_end = TimingHelpers::timer();
     oomph_info << "Time for p-refinement/unrefinement: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   // Loop over all nodes in mesh and free the dofs of those that were
   //-----------------------------------------------------------------
   // pinned only because they were hanging nodes. Also update their
   //-----------------------------------------------------------------
   // nodal values so that they contain data that is consistent
   //----------------------------------------------------------
   // with the hanging node representation
   //-------------------------------------
   // (Even if the nodal data isn't actually accessed because the node 
   // is still hanging -- we don't know this yet, and this step makes
   // sure that all nodes are fully functional and up-to-date, should
   // they become non-hanging below).
   //
   //
   // However, if we have a fixed mesh and hanging nodes on the boundary 
   // become non-hanging they will not necessarily respect the curvilinear
   // boundaries. This can only happen in 3D of course because it is not
   // possible to have hanging nodes on boundaries in 2D.
   // The solution is to store those nodes on the boundaries that are
   // currently hanging and then check to see whether they have changed
   // status at the end of the refinement procedure.
   // If it has changed, then we need to adjust their positions.
   const unsigned n_boundary = this->nboundary();
   const unsigned mesh_dim = this->finite_element_pt(0)->dim();
   Vector<std::set<Node*> > hanging_nodes_on_boundary_pt(n_boundary);

   unsigned long n_node=this->nnode();
   for (unsigned long n=0;n<n_node;n++)
    {
     //Get the pointer to the node
     Node* nod_pt=this->node_pt(n);
      
     //Get the number of values in the node
     unsigned n_value=nod_pt->nvalue();
      
     //We need to find if any of the values are hanging
     bool is_hanging = nod_pt->is_hanging();
     //Loop over the values and find out whether any are hanging
     for(unsigned n=0;n<n_value;n++)
      {is_hanging |= nod_pt->is_hanging(n);}
      
     //If the node is hanging then ...
     if(is_hanging)
      {
       // Unless they are turned into hanging nodes again below
       // (this might or might not happen), fill in all the necessary
       // data to make them 'proper' nodes again.
        
       // Reconstruct the nodal values/position from the node's 
       // hanging node representation
       unsigned nt=nod_pt->ntstorage();
       Vector<double> values(n_value);
       unsigned n_dim=nod_pt->ndim();
       Vector<double> position(n_dim);
       // Loop over all history values
       for(unsigned t=0;t<nt;t++)
        {
         nod_pt->value(t,values);
         for(unsigned i=0;i<n_value;i++) {nod_pt->set_value(t,i,values[i]);}
         nod_pt->position(t,position);
         for(unsigned i=0;i<n_dim;i++) {nod_pt->x(t,i)=position[i];}
        }
      
       // If it's an algebraic node: Update its previous nodal positions too
       AlgebraicNode* alg_node_pt=dynamic_cast<AlgebraicNode*>(nod_pt);
       if (alg_node_pt!=0)
        {
         bool update_all_time_levels=true;
         alg_node_pt->node_update(update_all_time_levels);
        }
        
        
       //If it's a Solid node, update Lagrangian coordinates
       // from its hanging node representation
       SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
       if(solid_node_pt!=0)
        {
         unsigned n_lagrangian = solid_node_pt->nlagrangian();
         for(unsigned i=0;i<n_lagrangian;i++)
          {
           solid_node_pt->xi(i) = solid_node_pt->lagrangian_position(i);
          }
        }

       //Now store geometrically hanging nodes on boundaries that
       //may need updating after refinement.
       //There will only be a problem if we have 3 spatial dimensions
       if((mesh_dim > 2) && (nod_pt->is_hanging()))
        {
         //If the node is on a boundary then add a pointer to the node
         //to our lookup scheme
         if(nod_pt->is_on_boundary())
          {
           //Storage for the boundaries on which the Node is located
           std::set<unsigned>* boundaries_pt;
           nod_pt->get_boundaries_pt(boundaries_pt);
           if(boundaries_pt!=0)
            {
             //Loop over the boundaries and add a pointer to the node
             //to the appropriate storage scheme
             for(std::set<unsigned>::iterator it=boundaries_pt->begin();
                 it!=boundaries_pt->end();++it)
              {
               hanging_nodes_on_boundary_pt[*it].insert(nod_pt);
              }
            }
          }
        }

      } //End of is_hanging
      
     // Initially mark all nodes as 'non-hanging' and `obsolete' 
     nod_pt->set_nonhanging();
     nod_pt->set_obsolete();
    }
   
   double t_end=0.0;
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for sorting out initial hanging status: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
      
   // Stick all elements into a vector
   Vector<Tree*> tree_nodes_pt;
   this->forest_pt()->stick_leaves_into_vector(tree_nodes_pt);
     
   //Copy the elements into the mesh Vector
   unsigned long num_tree_nodes=tree_nodes_pt.size();
   this->element_pt().resize(num_tree_nodes);
   for (unsigned long e=0;e<num_tree_nodes;e++)
    {
     this->element_pt(e)=tree_nodes_pt[e]->object_pt();
     
     // Now loop over all nodes in element and mark them as non-obsolete
     FiniteElement* this_el_pt=this->finite_element_pt(e);
     unsigned n_node=this_el_pt->nnode(); // caching pre-loop
     for (unsigned n=0;n<n_node;n++)
      {
       this_el_pt->node_pt(n)->set_non_obsolete();
      }
       
     // Mark up so that repeated refinements do not occur
     // (Required because refined element is the same element as the original)
     PRefineableElement* cast_el_pt =
      dynamic_cast<PRefineableElement*>(this->element_pt(e));
     cast_el_pt->deselect_for_p_refinement();
     cast_el_pt->deselect_for_p_unrefinement();
    }
    
   // Cannot delete nodes that are still marked as obsolete
   // because they may still be required to assemble the hanging schemes
   //-------------------------------------------------------------------
 
   // Mark up hanging nodes
   //----------------------
  
   //Output streams for the hanging nodes
   Vector<std::ofstream*> hanging_output_files;
   //Setup the output files for hanging nodes, this must be called
   //precisely once for the forest. Note that the files will only
   //actually be opened if doc_info.Doc_flag is true
   this->forest_pt()->open_hanging_node_files(doc_info,hanging_output_files);
 
   for(unsigned long e=0;e<num_tree_nodes;e++)
    {
     //Generic setup
     tree_nodes_pt[e]->object_pt()->setup_hanging_nodes(hanging_output_files);
     //Element specific setup
     tree_nodes_pt[e]->object_pt()->further_setup_hanging_nodes();
    }
 
   //Close the hanging node files and delete the memory allocated 
   //for the streams
   this->forest_pt()->close_hanging_node_files(doc_info,hanging_output_files);


   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info 
     <<"Time for setup_hanging_nodes() and further_setup_hanging_nodes() for " 
     << num_tree_nodes << " elements: "
     << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
 
   // Read out the number of continously interpolated values
   // from one of the elements (assuming it's the same in all elements)
   unsigned ncont_interpolated_values=
    tree_nodes_pt[0]->object_pt()->ncont_interpolated_values();
 
   // Complete the hanging nodes schemes by dealing with the
   // recursively hanging nodes
   this->complete_hanging_nodes(ncont_interpolated_values);


   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info 
      <<"Time for complete_hanging_nodes: "
      << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
 
   /// Update the boundary element info -- this can be a costly procedure
   /// and for this reason the mesh writer might have decided not to set up this
   /// scheme. If so, we won't change this and suppress its creation...
   if (Lookup_for_elements_next_boundary_is_setup)
    {
     this->setup_boundary_element_info(); 
    }
   
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info
      <<"Time for boundary element info: " 
      << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   //BENFLAG: Reset all the node update elements.
   //         This is necessary to prevent the following case: A node N is shared between two elements,
   //         A and B. The update element for the node is set to A, say. Element A is p-refined and now
   //         nolonger has N as a node. However the node update element for N is still A but the node
   //         doesn't exist in A.
   MacroElementNodeUpdateElementBase* first_macro_el_pt = dynamic_cast<MacroElementNodeUpdateElementBase*>(this->element_pt(0));
   if(first_macro_el_pt!=0)
    {
     // Now set the node update info elementwise
     for(unsigned e=0; e<this->nelement(); e++)
      {
       // Cast to macro element
       MacroElementNodeUpdateElementBase* macro_el_pt
        = dynamic_cast<MacroElementNodeUpdateElementBase*>(this->element_pt(e));
       if(macro_el_pt!=0)
        {
         // Get vector of geometric objects from element (construct vector
         // via copy operation)
         Vector<GeomObject*> geom_object_pt(macro_el_pt->geom_object_pt());

         // (Re)set node update info for all the nodes in the element
         macro_el_pt->set_node_update_info(geom_object_pt);
        }
      }
    }
 
#ifdef PARANOID
  
   // Doc/check the neighbours
   //-------------------------
   Vector<Tree*> all_tree_nodes_pt;
   this->forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
 
   //Check the neighbours
   this->forest_pt()->check_all_neighbours(doc_info);
  
   // Check the integrity of the elements
   // -----------------------------------
  
   // Loop over elements and get the elemental integrity
   double max_error=0.0;
   for (unsigned long e=0;e<num_tree_nodes;e++)
    {
     double max_el_error;
     tree_nodes_pt[e]->object_pt()->check_integrity(max_el_error);
     //If the elemental error is greater than our maximum error
     //reset the maximum
     if(max_el_error > max_error) {max_error=max_el_error;}
    }
 
   if (max_error>RefineableElement::max_integrity_tolerance())
    {
     std::ostringstream error_stream;
     error_stream << "Mesh refined: Max. error in integrity check: " 
                  << max_error << " is too big"
                  << "\ni.e. bigger than RefineableElement::"
                  << "max_integrity_tolerance()="
                  << RefineableElement::max_integrity_tolerance()
                  << std::endl;
 
     std::ofstream some_file;
     some_file.open("ProblemMesh.dat");
     for (unsigned long n=0;n<n_node;n++)
      {
       //Get the pointer to the node
       Node* nod_pt = this->node_pt(n);
       //Get the dimension
       unsigned n_dim = nod_pt->ndim();
       //Output the coordinates
       for(unsigned i=0;i<n_dim;i++)
        {
         some_file << this->node_pt(n)->x(i) << " ";
        }
       some_file << std::endl;
      }
     some_file.close();
 
     error_stream << "Documented problem mesh in ProblemMesh.dat" << std::endl;
 
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   else
    {
     oomph_info << "Mesh refined: Max. error in integrity check: " 
                << max_error << " is OK" << std::endl;
     oomph_info << "i.e. less than RefineableElement::max_integrity_tolerance()="
                << RefineableElement::max_integrity_tolerance()
                << "\n" << std::endl;
    }
 

   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info
      << "Time for (paranoid only) checking of integrity: " 
      << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
  
#endif
    
   //Loop over all elements other than the final level and deactivate the
   //objects, essentially set the pointer that point to nodes that are
   //about to be deleted to NULL. This must take place here because nodes
   //addressed by elements that are dead but still living in the tree might
   //have been made obsolete in the last round of refinement
   //(Not strictly required, as tree structure has not changed, but does no harm)
   for (unsigned long e=0;e<this->forest_pt()->ntree();e++)
    {
     this->forest_pt()->tree_pt(e)->
      traverse_all(&Tree::deactivate_object);
    }
 
   //Now we can prune the dead nodes from the mesh.
   Vector<Node*> deleted_node_pt=this->prune_dead_nodes();

   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for deactivating objects and pruning nodes: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }

   // Finally: Reorder the nodes within the mesh's node vector
   // to establish a standard ordering regardless of the sequence
   // of mesh refinements -- this is required to allow dump/restart
   // on refined meshes
   this->reorder_nodes();
  
   if (Global_timings::Doc_comprehensive_timings)
    {
     t_end = TimingHelpers::timer();
     oomph_info << "Time for reordering " << nnode() << " nodes: " 
                << t_end-t_start << std::endl;
     t_start = TimingHelpers::timer();
    }
   
   //Now we can correct the nodes on boundaries that were hanging that
   //are no longer hanging
   //Only bother if we have more than two dimensions
   if(mesh_dim > 2)
    {
     //Loop over the boundaries
     for(unsigned b=0;b<n_boundary;b++)
      {

       // Remove deleted nodes from the set
       unsigned n_del=deleted_node_pt.size();
       for (unsigned j=0;j<n_del;j++)
        {
         hanging_nodes_on_boundary_pt[b].erase(deleted_node_pt[j]);
        }
       
       //If the nodes that were hanging are still hanging then remove them
       //from the set (note increment is not in for command for efficiencty)
       for(std::set<Node*>::iterator 
            it=hanging_nodes_on_boundary_pt[b].begin();
            it!=hanging_nodes_on_boundary_pt[b].end();)
        {
         if((*it)->is_hanging()) 
          {hanging_nodes_on_boundary_pt[b].erase(it++);}
         else {++it;}
        }

       //Are there any nodes that have changed geometric hanging status
       //on the boundary
       //The slightly painful part is that we must adjust the position
       //via the macro-elements which are only available through the
       //elements and not the nodes.
       if(hanging_nodes_on_boundary_pt[b].size()>0)
        {
         //If so we loop over all elements adjacent to the boundary
         unsigned n_boundary_element = this->nboundary_element(b);
         for(unsigned e=0;e<n_boundary_element;++e)
          {
           //Get a pointer to the element
           FiniteElement* el_pt = this->boundary_element_pt(b,e);

           //Do we have a solid element
           SolidFiniteElement* solid_el_pt = 
            dynamic_cast<SolidFiniteElement*>(el_pt);
           
           //Determine whether there is a macro element
           bool macro_present = (el_pt->macro_elem_pt()!=0);
           //Or a solid macro element
           if(solid_el_pt!=0)
            {
             macro_present |= (solid_el_pt->undeformed_macro_elem_pt()!=0);
            }
           
           //Only bother to do anything if there is a macro element
           //or undeformed macro element in a SolidElement
           if(macro_present)
            {
             //Loop over the nodes 
             //ALH: (could optimise to only loop over
             //node associated with the boundary with more effort)
             unsigned n_el_node = el_pt->nnode();
             for(unsigned n=0;n<n_el_node;n++)
              {
               //Cache pointer to the node
               Node* nod_pt = el_pt->node_pt(n);
               if(nod_pt->is_on_boundary(b))
                {
                 //Is the Node in our set
                 std::set<Node*>::iterator it = 
                  hanging_nodes_on_boundary_pt[b].find(nod_pt);
                 //If we have found the Node then update the position
                 //to be consistent with the macro-element representation
                 if(it!=hanging_nodes_on_boundary_pt[b].end())
                  {
                   //Specialise local and global coordinates to 3D
                   //because there is only a problem in 3D.
                   Vector<double> s(3), x(3);
                   //Find the local coordinate of the ndoe
                   el_pt->local_coordinate_of_node(n,s);
                   //Find the number of time history values
                   const unsigned ntstorage = nod_pt->ntstorage();
                   
                   //Do we have a solid node
                   SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
                   if(solid_node_pt)
                    {

                     // Assign Lagrangian coordinates from undeformed
                     // macro element (if it has one -- get_x_and_xi()
                     // does "the right thing" anyway. Leave actual
                     // nodal positions alone -- we're doing a solid
                     // mechanics problem and once we're going
                     // the nodal positions are always computed, never
                     // (automatically) reset to macro-element based
                     // positions; not even on pinned boundaries
                     // because the user may have other ideas about where
                     // these should go -- pinning means "don't touch the
                     // value", not "leave where the macro-element thinks
                     // it should be"
                     Vector<double> x_fe(3), xi(3), xi_fe(3);
                     solid_el_pt->get_x_and_xi(s,x_fe,x,xi_fe,xi);
                     for(unsigned i=0;i<3;i++) {solid_node_pt->xi(i) = xi[i];}
                    }
                   else
                    {
                     //Set position and history values from the macro-element
                     //representation
                     for(unsigned t=0;t<ntstorage;t++)
                      {
                       //Get the history value from the macro element
                       el_pt->get_x(t,s,x);

                       //Set the coordinate to that of the macroelement
                       //representation
                       for(unsigned i=0;i<3;i++) {nod_pt->x(t,i) = x[i];}
                      }
                    } //End of non-solid node case
                   
                   //Now remove the node from the list
                   hanging_nodes_on_boundary_pt[b].erase(it);
                   //If there are no Nodes left then exit the loops
                   if(hanging_nodes_on_boundary_pt[b].size()==0)
                    {e=n_boundary_element; break;}
                  }
                }
              }
            } //End of macro element case
          }
        }
      }
    } //End of case when we have fixed nodal positions
    
   // Final doc
   //-----------
   if (doc_info.is_doc_enabled())
    {
     // Doc the boundary conditions ('0' for non-existent, '1' for free, 
     //----------------------------------------------------------------
     // '2' for pinned -- ideal for tecplot scatter sizing.
     //----------------------------------------------------
     //num_tree_nodes=tree_nodes_pt.size();
 
     // Determine maximum number of values at any node in this type of element
     RefineableElement* el_pt = tree_nodes_pt[0]->object_pt();
     //Initalise max_nval
     unsigned max_nval=0;
     for (unsigned n=0;n<el_pt->nnode();n++)
      {
       if (el_pt->node_pt(n)->nvalue()>max_nval)
        {max_nval=el_pt->node_pt(n)->nvalue();}
      }
 
     //Open the output file
     std::ostringstream fullname;
     std::ofstream bcs_file;
     fullname.str("");
     fullname << doc_info.directory() << "/bcs" << doc_info.number()
              << ".dat";
     bcs_file.open(fullname.str().c_str());  
    
     // Loop over elements
     for(unsigned long e=0;e<num_tree_nodes;e++)
      {
       el_pt = tree_nodes_pt[e]->object_pt();
       // Loop over nodes in element
       unsigned n_nod=el_pt->nnode();
       for(unsigned n=0;n<n_nod;n++)
        {
         //Get pointer to the node
         Node* nod_pt=el_pt->node_pt(n);
         //Find the dimension of the node
         unsigned n_dim = nod_pt->ndim();
         //Write the nodal coordinates to the file
         for(unsigned i=0;i<n_dim;i++)
          {bcs_file << nod_pt->x(i) << " ";}
        
         // Loop over all values in this element
         for(unsigned i=0;i<max_nval;i++)
          {
           // Value exists at this node:
           if (i<nod_pt->nvalue())
            {
             bcs_file << " " << 1+nod_pt->is_pinned(i);
            }
           // ...if not just dump out a zero
           else
            {
             bcs_file << " 0 ";
            }
          }
         bcs_file << std::endl;
        }  
      }
     bcs_file.close();
    
     // Doc all nodes
     //---------------
     std::ofstream all_nodes_file;
     fullname.str("");
     fullname << doc_info.directory() << "/all_nodes"
              << doc_info.number()  << ".dat";
     all_nodes_file.open(fullname.str().c_str());  
   
     all_nodes_file << "ZONE \n"; 
      
     // Need to recompute the number of nodes since it may have
     // changed during mesh refinement/unrefinement
     n_node = this->nnode();
     for(unsigned long n=0;n<n_node;n++)
      {
       Node* nod_pt = this->node_pt(n);
       unsigned n_dim = nod_pt->ndim();
       for(unsigned i=0;i<n_dim;i++)
        {
         all_nodes_file << this->node_pt(n)->x(i) << " ";
        }
       all_nodes_file << std::endl;
      }
      
     all_nodes_file.close();
 
 
     // Doc all hanging nodes:
     //-----------------------
     std::ofstream some_file;
     fullname.str("");
     fullname << doc_info.directory() << "/all_hangnodes"
              << doc_info.number() << ".dat";
     some_file.open(fullname.str().c_str());
     for(unsigned long n=0;n<n_node;n++)
      {
       Node* nod_pt=this->node_pt(n);
 
       if (nod_pt->is_hanging())
        {
         unsigned n_dim = nod_pt->ndim();       
         for(unsigned i=0;i<n_dim;i++)
          {
           some_file << nod_pt->x(i) << " ";
          }
        
         //ALH: Added this to stop Solid problems seg-faulting
         if(this->node_pt(n)->nvalue() > 0)
          {
           some_file  << " " << nod_pt->raw_value(0);
          }
         some_file << std::endl;
        }
      }
     some_file.close();
 
     // Doc all hanging nodes and their masters 
     // View with QHangingNodesWithMasters.mcr
     fullname.str("");
     fullname << doc_info.directory() 
              << "/geometric_hangnodes_withmasters" 
              << doc_info.number() << ".dat";
     some_file.open(fullname.str().c_str());
     for(unsigned long n=0;n<n_node;n++)
      {
       Node* nod_pt=this->node_pt(n);
       if (nod_pt->is_hanging())
        {
         unsigned n_dim = nod_pt->ndim();
         unsigned nmaster=nod_pt->hanging_pt()->nmaster();
         some_file << "ZONE I="<<nmaster+1 << std::endl;
         for(unsigned i=0;i<n_dim;i++)
          {
           some_file << nod_pt->x(i) << " ";
          }
         some_file << " 2 " <<  std::endl;
        
         for (unsigned imaster=0;imaster<nmaster;imaster++)
          {
           Node* master_nod_pt = 
            nod_pt->hanging_pt()->master_node_pt(imaster);
           unsigned n_dim = master_nod_pt->ndim();
           for(unsigned i=0;i<n_dim;i++)
            {
             some_file << master_nod_pt->x(i) << " ";
            }
           some_file << " 1 " << std::endl;
          }
        }
      }
     some_file.close();
    
     // Doc all hanging nodes and their masters 
     // View with QHangingNodesWithMasters.mcr
     for(unsigned i=0;i<ncont_interpolated_values;i++)
      {
       fullname.str("");
       fullname << doc_info.directory()
                <<"/nonstandard_hangnodes_withmasters" << i << "_" 
                << doc_info.number() << ".dat";
       some_file.open(fullname.str().c_str());
       unsigned n_nod=this->nnode();
       for(unsigned long n=0;n<n_nod;n++)
        {
         Node* nod_pt=this->node_pt(n);
         if (nod_pt->is_hanging(i))
          {
           if (nod_pt->hanging_pt(i)!=nod_pt->hanging_pt())
            {
             unsigned nmaster=nod_pt->hanging_pt(i)->nmaster();
             some_file << "ZONE I="<<nmaster+1 << std::endl;
             unsigned n_dim = nod_pt->ndim();
             for(unsigned j=0;j<n_dim;j++)
              {
               some_file << nod_pt->x(j) << " ";
              }
             some_file << " 2 " << std::endl;
             for (unsigned imaster=0;imaster<nmaster;imaster++)
              {
               Node* master_nod_pt = 
                nod_pt->hanging_pt(i)->master_node_pt(imaster);
               unsigned n_dim = master_nod_pt->ndim();
               for(unsigned j=0;j<n_dim;j++)
                {
//               some_file << master_nod_pt->x(i) << " ";
                }
               some_file << " 1 " << std::endl;
              }
            }
          }
        }
       some_file.close();
      }
    
    } //End of documentation
  } // End if (this->nelement()>0)

 ////BENFLAG: Check that all the nodes belong to their update elements
 //std::cout << "p_adapt_mesh(): Checking stuff works!" << std::endl;
 //for(unsigned j=0; j<this->nnode(); j++)
 // {
 //  MacroElementNodeUpdateNode* macro_nod_pt = dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j));
 //  if(macro_nod_pt!=0)
 //   {
 //    bool big_problem = true;
 //    std::cout << "Node " << macro_nod_pt << " at [ " << macro_nod_pt->x(0) << ", " << macro_nod_pt->x(1) << " ]" << std::endl;
 //    FiniteElement* up_el_pt = dynamic_cast<FiniteElement*>(macro_nod_pt->node_update_element_pt());
 //    for(unsigned l=0; l<up_el_pt->nnode(); l++)
 //     {
 //      if(up_el_pt->node_pt(l)==macro_nod_pt)
 //       {
 //        big_problem = false;
 //        break;
 //       }
 //     }
 //    if(big_problem)
 //     {
 //      std::cout << "  This node doesn't exist in it's update element!" << std::endl;
 //     }
 //   }
 // }

#ifdef OOMPH_HAS_MPI

 // Now (re-)classify halo and haloed nodes and synchronise hanging
 // nodes
 if (this->is_mesh_distributed())
  {
   classify_halo_and_haloed_nodes(doc_info,doc_info.is_doc_enabled());
  } 

#endif

}

//========================================================================
/// p-unrefine mesh uniformly
/// Unlike in h-refinement, we can simply p-unrefine each element in the mesh
//========================================================================
void TreeBasedRefineableMeshBase::p_unrefine_uniformly(DocInfo& doc_info)
{ 
 //Select all elements for unrefinement
 unsigned long Nelement=this->nelement();
 for (unsigned long e=0;e<Nelement;e++)
  {
   //Get pointer to p-refineable element
   PRefineableElement* el_pt
    = dynamic_cast<PRefineableElement*>(this->element_pt(e));
   //Mark for p-refinement if possible. If not then p_adapt_mesh() will
   //report the error.
   if(el_pt!=0)
    {
     el_pt->select_for_p_unrefinement();
    }
  }
 
 // Do the actual mesh adaptation
 p_adapt_mesh(doc_info);  
}

//========================================================================
/// p-refine mesh by refining the elements identified by their numbers.
//========================================================================
void TreeBasedRefineableMeshBase::p_refine_selected_elements(
 const Vector<unsigned>& elements_to_be_refined)
{
 
#ifdef OOMPH_HAS_MPI
 if(this->is_mesh_distributed())
  {
   std::ostringstream warn_stream;
   warn_stream << "You are attempting to refine selected elements of a "
               << std::endl
               << "distributed mesh. This may have undesired effects."
               << std::endl;
   
   OomphLibWarning(warn_stream.str(),
                   "TreeBasedRefineableMeshBase::refine_selected_elements()",
                   OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 //Select elements for refinement
 unsigned long nref=elements_to_be_refined.size();
 for (unsigned long e=0;e<nref;e++)
  {
   //Get pointer to p-refineable element
   PRefineableElement* el_pt
    = dynamic_cast<PRefineableElement*>
       (this->element_pt(elements_to_be_refined[e]));
   //Mark for p-refinement if possible. If not then p_adapt_mesh() will
   //report the error.
   if(el_pt!=0)
    {
     el_pt->select_for_p_refinement();
    }
  }
 
 // Do the actual mesh adaptation
 p_adapt_mesh(); 
}

//========================================================================
/// p-refine mesh by refining the elements identified by their pointers.
//========================================================================
void TreeBasedRefineableMeshBase::p_refine_selected_elements(
 const Vector<PRefineableElement*>& elements_to_be_refined_pt)
{
 
#ifdef OOMPH_HAS_MPI
 if(this->is_mesh_distributed())
  {
   std::ostringstream warn_stream;
   warn_stream << "You are attempting to refine selected elements of a "
               << std::endl
               << "distributed mesh. This may have undesired effects."
               << std::endl;
   
   OomphLibWarning(warn_stream.str(),
                   "TreeBasedRefineableMeshBase::refine_selected_elements()",
                   OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 //Select elements for refinement
 unsigned long nref=elements_to_be_refined_pt.size();
 for (unsigned long e=0;e<nref;e++)
  {
   elements_to_be_refined_pt[e]->select_for_p_refinement();
  }
 
 // Do the actual mesh adaptation
 p_adapt_mesh(); 
}

}