refineable_mesh.h

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//====================================================================
#ifndef OOMPH_REFINEABLE_MESH_HEADER
#define OOMPH_REFINEABLE_MESH_HEADER

#include "mesh.h"
#include "refineable_elements.h"
//Include the tree template to fill in the C++ split function
//Must be called after refineable_element.h
#include "tree.template.cc"
#include "error_estimator.h" 

namespace oomph
{


  



//=======================================================================
/// Base class for refineable meshes. Provides standardised interfaces
/// for the following standard mesh adaptation routines:
/// - Uniform mesh refinement.
/// - Adaptation, based on the elemental error estimates obtained
///   from an error estimator function which is accessed via a function 
///   pointer. 
///
//=======================================================================
class RefineableMeshBase : public virtual Mesh
{
public:

 bool adapt_flag() {return Adapt_flag;}

 /// Constructor sets default values for refinement targets etc.
 /// and initialises pointer to spatial error estimator to NULL.
 RefineableMeshBase()
  {
   // Initialise pointer to spatial error estimator
   Spatial_error_estimator_pt=0;
 
   // Targets
   Min_permitted_error=1.0e-5;
   Max_permitted_error=1.0e-3;

   //Actual errors
   Min_error = 0.0;
   Max_error = 0.0;
    
   // Do I adapt or not?
   Adapt_flag=true;
    
   // Do I p-adapt or not?
   P_adapt_flag=true;
    
   // Do I disable additional synchronisation of hanging nodes?
   Additional_synchronisation_of_hanging_nodes_not_required=false;
 
   // Where do I write the documatation of the refinement process?
   Doc_info_pt=0;

   // If only a few elements are scheduled for unrefinement, don't bother
   // By default unrefine all
   Max_keep_unrefined=0;

   // Number of elements where the refinement is over-ruled
   Nrefinement_overruled = 0;
  };


 /// Broken copy constructor
 RefineableMeshBase(const RefineableMeshBase& dummy) 
  { 
   BrokenCopy::broken_copy("RefineableMeshBase");
  } 
 
 /// Broken assignment operator
 void operator=(const RefineableMeshBase&) 
  {
   BrokenCopy::broken_assign("RefineableMeshBase");
  }

 /// Empty Destructor: 
 virtual ~RefineableMeshBase()
  {
  }

 ///  Access fct for number of elements that were refined
 unsigned nrefined() {return Nrefined;}
 
 /// Access fct for  number of elements that were unrefined
 unsigned nunrefined() {return Nunrefined;}

 /// \short Number of elements that would have liked to be refined further 
 /// but can't because they've reached the max. refinement level
 unsigned& nrefinement_overruled(){return Nrefinement_overruled;}

 /// \short Max. number of elements that we allow to remain unrefined
 /// if no other mesh adaptation is required (to avoid 
 /// mesh-adaptations that would only unrefine a few elements
 /// and then force a new solve -- this can't be worth our while!)
 unsigned& max_keep_unrefined(){return Max_keep_unrefined;}
 
 /// Doc the targets for mesh adaptation
 virtual void doc_adaptivity_targets(std::ostream &outfile)
  {
   outfile << std::endl;
   outfile << "Targets for mesh adaptation: " << std::endl;
   outfile << "---------------------------- " << std::endl;
   outfile << "Target for max. error: " << Max_permitted_error << std::endl;
   outfile << "Target for min. error: " << Min_permitted_error << std::endl;
   outfile << "Don't unrefine if less than " << Max_keep_unrefined 
           << " elements need unrefinement." << std::endl;
   outfile << std::endl;
  }


 /// Access to spatial error estimator
 ErrorEstimator*& spatial_error_estimator_pt()
  {
   return Spatial_error_estimator_pt;
  }

 /// Access to spatial error estimator (const version
 ErrorEstimator* spatial_error_estimator_pt() const
  {
   return Spatial_error_estimator_pt;
  }

 /// \short Access fct for min. error (i.e. (try to) merge elements if 
 /// their error is smaller)
 double& min_permitted_error() {return Min_permitted_error;}

 /// \short Access fct for max. error (i.e. split elements if their 
 /// error is larger)
 double& max_permitted_error() {return Max_permitted_error;}

 /// \short Access fct for min. actual error in present solution (i.e. before
 /// re-solve on adapted mesh)
 double& min_error() {return Min_error;}

 /// \short Access fct for max. actual error in present solution (i.e. before
 /// re-solve on adapted mesh)
 double& max_error() {return Max_error;}
 
 ///  Access fct for pointer to DocInfo
 DocInfo*& doc_info_pt() {return Doc_info_pt;}
 
 /// Enable adaptation
 void enable_adaptation() {Adapt_flag=true;}

 /// Disable adaptation
 void disable_adaptation() {Adapt_flag=false;}
 
 /// Enable adaptation
 void enable_p_adaptation() {P_adapt_flag=true;}

 /// Disable adaptation
 void disable_p_adaptation() {P_adapt_flag=false;}
 
 /// Enable additional synchronisation of hanging nodes
 void enable_additional_synchronisation_of_hanging_nodes()
  {Additional_synchronisation_of_hanging_nodes_not_required=false;}

 /// Disable additional synchronisation of hanging nodes
 void disable_additional_synchronisation_of_hanging_nodes()
  {Additional_synchronisation_of_hanging_nodes_not_required=true;}

 /// Return whether the mesh is to be adapted
 bool is_adaptation_enabled() const {return Adapt_flag;}

 /// Return whether the mesh is to be adapted
 bool is_p_adaptation_enabled() const {return P_adapt_flag;}

 /// Return whether additional synchronisation is enabled
 bool is_additional_synchronisation_of_hanging_nodes_disabled() const
  {return Additional_synchronisation_of_hanging_nodes_not_required;}
 
 ///  Access fct for DocInfo
 DocInfo doc_info()
  {
   return *Doc_info_pt;
  }
 
 /// \short Adapt mesh: Refine elements whose error is lager than err_max
 /// and (try to) unrefine those whose error is smaller than err_min
 virtual void adapt(const Vector<double>& elemental_error)=0;
 
 /// \short p-adapt mesh: Refine elements whose error is lager than err_max
 /// and (try to) unrefine those whose error is smaller than err_min
 virtual void p_adapt(const Vector<double>& elemental_error)
  {
   //Derived classes must implement this as required. Default throws an error.
   std::ostringstream err_stream;
   err_stream << "p_adapt() called in base class RefineableMeshBase."
              << std::endl
              << "This needs to be implemented in the derived class."
              << std::endl;
   throw OomphLibError(err_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

 /// Refine mesh uniformly and doc process
 virtual void refine_uniformly(DocInfo& doc_info)=0;
 
 /// Refine mesh uniformly
 virtual void refine_uniformly()
  {
   DocInfo doc_info;
   doc_info.directory()="";
   doc_info.disable_doc();
   refine_uniformly(doc_info);
  }

 /// p-refine mesh uniformly and doc process
 virtual void p_refine_uniformly(DocInfo& doc_info)
  {
   //Derived classes must implement this as required. Default throws an error.
   std::ostringstream err_stream;
   err_stream << "p_refine_uniformly() called in base class RefineableMeshBase."
              << std::endl
              << "This needs to be implemented in the derived class."
              << std::endl;
   throw OomphLibError(err_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
 
 /// p-refine mesh uniformly
 virtual void p_refine_uniformly()
  {
   DocInfo doc_info;
   doc_info.directory()="";
   doc_info.disable_doc();
   p_refine_uniformly(doc_info);
  }

 /// \short Unrefine mesh uniformly: Return 0 for success,
 /// 1 for failure (if unrefinement has reached the coarsest permitted
 /// level)
 virtual unsigned unrefine_uniformly()=0;

 /// \short p-unrefine mesh uniformly
 void p_unrefine_uniformly(DocInfo& doc_info)
  {
   //Derived classes must implement this as required. Default throws an error.
   std::ostringstream err_stream;
   err_stream << "p_unrefine_uniformly() called in base class RefineableMeshBase."
              << std::endl
              << "This needs to be implemented in the derived class."
              << std::endl;
   throw OomphLibError(err_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

protected:
 
 /// \short Pointer to spatial error estimator
 ErrorEstimator* Spatial_error_estimator_pt;

 /// Max. error (i.e. split elements if their error is larger)
 double Max_permitted_error;
 
 /// Min. error (i.e. (try to) merge elements if their error is smaller)
 double Min_permitted_error;
 
 /// Min.actual error
 double Min_error;
 
 /// Max. actual error 
 double Max_error;
 
 /// Stats: Number of elements that were refined
 unsigned Nrefined;

 /// Stats: Number of elements that were unrefined
 unsigned Nunrefined;

  /// Flag that requests adaptation
 bool Adapt_flag;

  /// Flag that requests p-adaptation
 bool P_adapt_flag;

  /// Flag that disables additional synchronisation of hanging nodes
 bool Additional_synchronisation_of_hanging_nodes_not_required;

 /// Pointer to DocInfo
 DocInfo* Doc_info_pt;

 /// \short Max. number of elements that can remain unrefined
 /// if no other mesh adaptation is required (to avoid 
 /// mesh-adaptations that would only unrefine a few elements
 /// and then force a new solve -- this can't be worth our while!)
 unsigned Max_keep_unrefined;
 
 /// \short Number of elements that would like to be refined further but can't
 /// because they've reached the max. refinement level
 unsigned Nrefinement_overruled;

};




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



//=======================================================================
/// Base class for tree-based refineable meshes. 
//=======================================================================
class TreeBasedRefineableMeshBase : public virtual RefineableMeshBase
{
public:


 /// Constructor 
 TreeBasedRefineableMeshBase()
  {
   // Max/min refinement levels
   Max_refinement_level=5;
   Min_refinement_level=0;
   
   // Max/min p-refinement levels
   Max_p_refinement_level=7;
   Min_p_refinement_level=2;
 
   // Stats
   Nrefined=0;
   Nunrefined=0;

   // Where do I write the documatation of the refinement process?
   Doc_info_pt=0;
   
   // Initialise the forest pointer to NULL
   Forest_pt = 0;

   // Mesh hasn't been pruned yet
   Uniform_refinement_level_when_pruned=0;
  }


 /// Broken copy constructor
 TreeBasedRefineableMeshBase(const TreeBasedRefineableMeshBase& dummy) 
  { 
   BrokenCopy::broken_copy("TreeBasedRefineableMeshBase");
  } 
 
 /// Broken assignment operator
 void operator=(const TreeBasedRefineableMeshBase&) 
  {
   BrokenCopy::broken_assign("TreeBasedRefineableMeshBase");
  }

 /// Empty Destructor: 
 virtual ~TreeBasedRefineableMeshBase()
  {
   //Kill the forest if there is one
   if(Forest_pt!=0)
    {
     delete Forest_pt;
     Forest_pt=0;
    }
  }

 /// \short Adapt mesh: Refine elements whose error is lager than err_max
 /// and (try to) unrefine those whose error is smaller than err_min
 void adapt(const Vector<double>& elemental_error);

 /// \short p-adapt mesh: Refine elements whose error is lager than err_max
 /// and (try to) unrefine those whose error is smaller than err_min
 void p_adapt(const Vector<double>& elemental_error);

 /// Refine mesh uniformly and doc process
 void refine_uniformly(DocInfo& doc_info);
 
 /// Refine mesh uniformly
 void refine_uniformly()
 {
  RefineableMeshBase::refine_uniformly();
 }
 
 /// p-refine mesh uniformly and doc process
 void p_refine_uniformly(DocInfo& doc_info);
 
 /// p-refine mesh uniformly
 void p_refine_uniformly()
  {
   RefineableMeshBase::p_refine_uniformly();
  }
 
 /// \short Unrefine mesh uniformly: Return 0 for success,
 /// 1 for failure (if unrefinement has reached the coarsest permitted
 /// level)
 unsigned unrefine_uniformly();
 
 /// \short p-unrefine mesh uniformly
 void p_unrefine_uniformly(DocInfo& doc_info);

 /// Set up the tree forest associated with the Mesh (if any)
 virtual void setup_tree_forest()=0;

 /// Return pointer to the Forest represenation of the mesh
 TreeForest* forest_pt(){return Forest_pt;}


 /// Doc the targets for mesh adaptation
 void doc_adaptivity_targets(std::ostream &outfile)
  {
   outfile << std::endl;
   outfile << "Targets for mesh adaptation: " << std::endl;
   outfile << "---------------------------- " << std::endl;
   outfile << "Target for max. error: " << Max_permitted_error << std::endl;
   outfile << "Target for min. error: " << Min_permitted_error << std::endl;
   outfile << "Min. refinement level: " << Min_refinement_level << std::endl;
   outfile << "Max. refinement level: " << Max_refinement_level << std::endl;
   outfile << "Min. p-refinement level: " << Min_p_refinement_level
           << std::endl;
   outfile << "Max. p-refinement level: " << Max_p_refinement_level
           << std::endl;
   outfile << "Don't unrefine if less than " << Max_keep_unrefined 
           << " elements need unrefinement." << std::endl;
   outfile << std::endl;
  }
 
 /// Access fct for max. permissible refinement level (relative to base mesh)
 unsigned& max_refinement_level() {return Max_refinement_level;}
 
 /// Access fct for min. permissible refinement level (relative to base mesh)
 unsigned& min_refinement_level() {return Min_refinement_level;}
 
 /// Access fct for max. permissible p-refinement level (relative to base mesh)
 unsigned& max_p_refinement_level() {return Max_p_refinement_level;}
 
 /// Access fct for min. permissible p-refinement level (relative to base mesh)
 unsigned& min_p_refinement_level() {return Min_p_refinement_level;}
 
 /// \short Perform the actual tree-based mesh adaptation, 
 /// documenting the progress in the directory specified in DocInfo object.
 virtual void adapt_mesh(DocInfo& doc_info);

 /// \short Perform the actual tree-based mesh adaptation. A simple wrapper to 
 /// call the function without documentation.
 virtual void adapt_mesh()
  {
   //Create a dummy doc_info object
   DocInfo doc_info;
   doc_info.directory()="";
   doc_info.disable_doc();
   //Call the other adapt mesh
   adapt_mesh(doc_info);
  }
 
 /// \short Perform the actual tree-based mesh p-adaptation, 
 /// documenting the progress in the directory specified in DocInfo object.
 void p_adapt_mesh(DocInfo& doc_info);

 /// \short Perform the actual tree-based mesh p-adaptation. A simple wrapper
 /// to call the function without documentation.
 void p_adapt_mesh()
  {
   //Create a dummy doc_info object
   DocInfo doc_info;
   doc_info.directory()="";
   doc_info.disable_doc();
   //Call the other adapt mesh
   p_adapt_mesh(doc_info);
  }

 /// \short Refine mesh by splitting the elements identified
 /// by their numbers.
 virtual void refine_selected_elements(const Vector<unsigned>& 
                                       elements_to_be_refined);

 /// \short Refine mesh by splitting the elements identified
 /// by their pointers.
 virtual void refine_selected_elements(const Vector<RefineableElement*>& 
                                       elements_to_be_refined);

 /// \short p-refine mesh by refining the elements identified
 /// by their numbers.
 void p_refine_selected_elements(const Vector<unsigned>& 
                                 elements_to_be_refined);

 /// \short p-refine mesh by refining the elements identified
 /// by their pointers.
 void p_refine_selected_elements(const Vector<PRefineableElement*>& 
                                 elements_to_be_refined_pt);

 
 /// \short Refine base mesh to same degree as reference mesh (relative
 /// to original unrefined mesh).
 virtual void refine_base_mesh_as_in_reference_mesh(
  TreeBasedRefineableMeshBase* const &ref_mesh_pt);
 
 /// \short Refine base mesh to same degree as reference mesh minus one
 /// level of refinement (relative to original unrefined mesh). Useful
 /// function for multigrid solvers; allows the easy copy of a mesh
 /// to the level of refinement just below the current one
 virtual bool refine_base_mesh_as_in_reference_mesh_minus_one(
  TreeBasedRefineableMeshBase* const &ref_mesh_pt);

 /// \short 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).
 virtual void refine_as_in_reference_mesh(TreeBasedRefineableMeshBase* 
                                          const &ref_mesh_pt);

 /// Get max/min refinement levels in mesh
 virtual void get_refinement_levels(unsigned& min_refinement_level,
                                    unsigned& max_refinement_level);

 /// \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)
 virtual void get_elements_at_refinement_level(
  unsigned& refinement_level, Vector<RefineableElement*>& level_elements);

 /// \short Extract refinement pattern: 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][...]
 virtual void get_refinement_pattern(Vector<Vector<unsigned> >& 
                                     to_be_refined);

 /// Refine base mesh according to specified refinement pattern
 void refine_base_mesh(Vector<Vector<unsigned> >& to_be_refined);

 /// Refine mesh according to refinement pattern in restart file
 virtual void refine(std::ifstream& restart_file);

 /// Dump refinement pattern to allow for rebuild
 virtual void dump_refinement(std::ostream &outfile);

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


 /// \short Level to which the mesh was uniformly refined when it was pruned
 /// (const version)
 unsigned uniform_refinement_level_when_pruned() const
 {
  return Uniform_refinement_level_when_pruned;
 }


 /// \short Level to which the mesh was uniformly refined when it was pruned
 unsigned& uniform_refinement_level_when_pruned()
 {
  return Uniform_refinement_level_when_pruned;
 }
 
#ifdef OOMPH_HAS_MPI

 /// Classify all halo and haloed information in the mesh (overloaded
 /// version from Mesh base class. Calls that one first, then synchronises
 /// hanging nodes)
 void classify_halo_and_haloed_nodes(DocInfo& doc_info, 
                                     const bool& report_stats)
 {
  // Call version in base class but don't bother to call
  // resize_halo_nodes() -- we'll do it ourselves below
  bool backup=Resize_halo_nodes_not_required;
  Resize_halo_nodes_not_required=false;
  Mesh::classify_halo_and_haloed_nodes(doc_info,
                                       report_stats);
  Resize_halo_nodes_not_required=backup;

  double t_start=0.0;
  if (Global_timings::Doc_comprehensive_timings)
   {
    t_start = TimingHelpers::timer();
   }
  
  // Communicate 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)
  synchronise_nonhanging_nodes();

  // Get number of continously interpolated variables
  unsigned local_ncont_interpolated_values=0;
  
  // Bypass if we have no elements and then get overall
  // value by reduction
  if (nelement()>0)
   {         
    local_ncont_interpolated_values=dynamic_cast<RefineableElement*>
     (element_pt(0))->ncont_interpolated_values();
   }
  unsigned ncont_interpolated_values=0;
  MPI_Allreduce(&local_ncont_interpolated_values,
                &ncont_interpolated_values,1,
                MPI_UNSIGNED,MPI_MAX,
                Comm_pt->mpi_comm());
  
  // Synchronise the hanging nodes 
  synchronise_hanging_nodes(ncont_interpolated_values);

  if (Global_timings::Doc_comprehensive_timings)
   {
    double t_end = TimingHelpers::timer();
    oomph_info 
        << "Time for "
        << "TreeBasedRefineableMeshBase::synchronise_hanging_nodes() "
        << "incl. time for initial allreduce in " 
        << "TreeBasedRefineableMeshBase::classify_halo_and_haloed_nodes(): " 
        << t_end-t_start << std::endl;
   }

  // Now do resize halo nodes after all (unless it's been 
  // declared to be unnecessary by somebody...)
  if (!Resize_halo_nodes_not_required)
   {
    resize_halo_nodes();
   }
 }


 /// Classify the halo and haloed nodes in the mesh. 
 void classify_halo_and_haloed_nodes(const bool& report_stats=false)
 {
  DocInfo doc_info;
  doc_info.disable_doc();
  classify_halo_and_haloed_nodes(doc_info,report_stats);
 }

#endif

protected:
 
#ifdef OOMPH_HAS_MPI

 /// Synchronise the hanging nodes if the mesh is distributed
 void synchronise_hanging_nodes(const unsigned& ncont_interpolated_values);

 /// Additional synchronisation of hanging nodes
 /// Required for reconcilliation of hanging nodes on the outer edge of the
 /// halo layer when using elements with nonuniformly spaced nodes
 /// Must be implemented in templated derived class because ELEMENT template
 /// parameter is required for the node-creation functions in the
 /// Missing_masters_functions namespace
 virtual void additional_synchronise_hanging_nodes(
               const unsigned& ncont_interpolated_values)=0;

 /// 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 synchronise_nonhanging_nodes();
 
#endif

 /// \short Split all the elements in the mesh if required. This template free
 /// interface will be overloaded in RefineableMesh<ELEMENT> so that
 /// any new elements that are created will be of the correct type.
 virtual void split_elements_if_required()=0;

 /// \short p-refine all the elements in the mesh if required. This template
 /// free interface will be overloaded in RefineableMesh<ELEMENT> so that
 /// any temporary copies of the element that are created will be of the
 /// correct type.
 virtual void p_refine_elements_if_required()=0;

 /// \short Complete the hanging node scheme recursively
 void complete_hanging_nodes(const int& ncont_interpolated_values);

 
 /// Auxiliary routine for recursive hanging node completion
 void complete_hanging_nodes_recursively(Node*& nod_pt,
                                         Vector<Node*>& master_nodes,
                                         Vector<double>& hang_weights,
                                         const int &ival);
 
 /// \short Level to which the mesh was uniformly refined when it was pruned
 unsigned Uniform_refinement_level_when_pruned;
 
 /// Max. permissible refinement level (relative to base mesh)
 unsigned Max_refinement_level;

 /// Min. permissible refinement level (relative to base mesh)
 unsigned Min_refinement_level;
 
 /// Max. permissible p-refinement level (relative to base mesh)
 unsigned Max_p_refinement_level;

 /// Min. permissible p-refinement level (relative to base mesh)
 unsigned Min_p_refinement_level;

 /// Forest representation of the mesh
 TreeForest* Forest_pt;

  private:

#ifdef OOMPH_HAS_MPI
 
 /// \short Helper struct to collate data required during 
 /// TreeBasedRefineableMeshBase::synchronise_hanging_nodes
 struct HangHelperStruct
 {
  unsigned Sending_processor;
  unsigned Shared_node_id_on_sending_processor;
  unsigned Shared_node_proc;
  double Weight;
  HangInfo* Hang_pt;
  unsigned Master_node_index;
  // Store pointer to node and index of value in case translation attempt fails
  // (see synchronise_hanging_nodes())
  Node* Node_pt;
  int icont;
 };
 
#endif

};


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


//=======================================================================
/// Templated base class for refineable meshes. The use of the template
/// parameter is required only for creating new elements during mesh
/// adaptation. This class overloaded the template-free inteface to
/// the function split_elements_if_required() to make use of the template
/// parameter.
/// All refineable meshes should inherit directly from 
/// TreeBasedRefineableMesh<ELEMENT>
//=======================================================================
template <class ELEMENT>
class TreeBasedRefineableMesh : public virtual TreeBasedRefineableMeshBase
 {
   private:
   
  /// \short Split all the elements if required. Overload the template-free
  /// interface so that any new elements that are created 
  /// will be of the correct type.
  void split_elements_if_required()
   {
    //Find the number of trees in the forest
    unsigned n_tree = this->Forest_pt->ntree();
    //Loop over all "active" elements in the forest and split them
    //if required
    for (unsigned long e=0;e<n_tree;e++)
     {
      this->Forest_pt->tree_pt(e)->
       traverse_leaves(&Tree::split_if_required<ELEMENT>);
     }
   }
   
  /// \short p-refine all the elements if required. Overload the template-free
  /// interface so that any temporary copies of the element that are created
  /// will be of the correct type.
  void p_refine_elements_if_required()
   {
    //BENFLAG: Make a non const pointer to the mesh so it can be passed (HACK)
    Mesh* mesh_pt=this;
    //Find the number of trees in the forest
    unsigned n_tree = this->Forest_pt->ntree();
    //Loop over all "active" elements in the forest and p-refine them
    //if required
    for (unsigned long e=0;e<n_tree;e++)
     {
      this->Forest_pt->tree_pt(e)->
       traverse_leaves(&Tree::p_refine_if_required<ELEMENT>,mesh_pt);
     }
   }

   protected:

#ifdef OOMPH_HAS_MPI
  
  /// Additional setup of shared node scheme
  /// This is Required for reconcilliation of hanging nodes acrross processor
  /// boundaries when using elements with nonuniformly spaced nodes.
  /// ELEMENT template parameter is required so that
  /// MacroElementNodeUpdateNodes which are added as external halo master nodes
  /// can be made fully functional
  void additional_synchronise_hanging_nodes(
        const unsigned& ncont_interpolated_values);

#endif

   public:

  ///Constructor, call the constructor of the base class
  TreeBasedRefineableMesh() : TreeBasedRefineableMeshBase() {}
  
  /// Broken copy constructor
  TreeBasedRefineableMesh(const TreeBasedRefineableMesh& dummy) 
   { 
    BrokenCopy::broken_copy("TreeBasedRefineableMesh");
   } 
  
  ///Empty virtual destructor
  virtual ~TreeBasedRefineableMesh() { }

 };





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


//=======================================================================
/// Base class for refineable tet meshes
//=======================================================================
class RefineableTetMeshBase : public virtual RefineableMeshBase
 {


   public:

   /// Max element size allowed during adaptation
   double& max_element_size(){return Max_element_size;}
   
   /// Min element size allowed during adaptation
   double& min_element_size(){return Min_element_size;}
   
   /// Min edge ratio before remesh gets triggered
   double& max_permitted_edge_ratio(){return Max_permitted_edge_ratio;}
   

   /// Doc the targets for mesh adaptation
   void doc_adaptivity_targets(std::ostream &outfile)
   {
    outfile << std::endl;
    outfile << "Targets for mesh adaptation: " << std::endl;
    outfile << "---------------------------- " << std::endl;
    outfile << "Target for max. error: " << Max_permitted_error << std::endl;
    outfile << "Target for min. error: " << Min_permitted_error << std::endl;
    outfile << "Target max edge ratio: " 
            << Max_permitted_edge_ratio << std::endl;
    outfile << "Min. allowed element size: " << Min_element_size << std::endl;
    outfile << "Max. allowed element size: " << Max_element_size << std::endl;
    outfile << "Don't unrefine if less than " << Max_keep_unrefined 
            << " elements need unrefinement." << std::endl;
    outfile << std::endl;
   }
   


   /// \short Compute target volume based on the elements' error and the
   /// error target; return max edge ratio
   double compute_volume_target(const Vector<double> &elem_error,
                                Vector<double> &target_volume)
   {
     double max_edge_ratio=0.0;
     unsigned count_unrefined=0;
     unsigned count_refined=0;
     this->Nrefinement_overruled=0;
     
     unsigned nel=this->nelement();
     for (unsigned e=0;e<nel;e++)
      {
       // Get element
       FiniteElement* el_pt=this->finite_element_pt(e);
       
       // Calculate the volume of the element
       double volume=el_pt->size();
       
       //Find the vertex coordinates
       // (vertices are enumerated first)
       double vertex[4][3];
       for(unsigned n=0;n<4;++n)
        {
         for(unsigned i=0;i<3;++i)
          {
           vertex[n][i] = el_pt->node_pt(n)->x(i);
          }
        }
       
       //Compute the radius of the circumsphere of the tetrahedron
       //Algorithm stolen from tetgen for consistency
       DenseDoubleMatrix A(3);
       for(unsigned i=0;i<3;++i)
        {
         A(0,i) = vertex[1][i] - vertex[0][i];
         A(1,i) = vertex[2][i] - vertex[0][i];
         A(2,i) = vertex[3][i] - vertex[0][i];
        }
       
       Vector<double> rhs(3);
       // Compute the right hand side vector b (3x1).
       for(unsigned i=0;i<3;++i)
        {
         rhs[i] = 0.0;
         for(unsigned k=0;k<3;++k)
          {
           rhs[i] += A(i,k)*A(i,k);
          }
         rhs[i] *= 0.5;
        }
       
       //Solve the linear system, in which the rhs is over-written with 
       //the solution
       A.solve(rhs);

       //Calculate the circum-radius
       double circum_radius = 
        sqrt(rhs[0]*rhs[0] + rhs[1]*rhs[1] + rhs[2]*rhs[2]);
       
       //Now find the shortest edge length
       Vector<double> edge(3);
       double min_length = DBL_MAX;
       for(unsigned start=0;start<4;++start)
        {
         for(unsigned end=start+1;end<4;++end)
          {
           for(unsigned i=0;i<3;++i)
            {
             edge[i] = vertex[start][i] - vertex[end][i];
            }
           double length = 
            sqrt(edge[0]*edge[0] + edge[1]*edge[1] + edge[2]*edge[2]);
           if(length < min_length) {min_length = length;}
          }
        }
       
       //Now calculate the minimum edge ratio for this element
       double local_max_edge_ratio = circum_radius/min_length;
       if(local_max_edge_ratio > max_edge_ratio)
        {
         max_edge_ratio = local_max_edge_ratio;
        }

       // Mimick refinement in tree-based procedure: Target volumes
       // for elements that exceed permitted error is 1/4 of their
       // current volume, corresponding to a uniform sub-division.
       if (elem_error[e] > this->max_permitted_error())
        {
         target_volume[e]=std::max(volume/4.0,Min_element_size);
         if (target_volume[e]!=Min_element_size)
          {
           count_refined++;
          }
         else
          {
           this->Nrefinement_overruled++;
          }
        }
       else if (elem_error[e] < this->min_permitted_error())
        {
         target_volume[e]=std::min(4.0*volume,Max_element_size);
         if (target_volume[e]!=Max_element_size)
          {
           count_unrefined++;
          }
        }
       else
        {
         // Leave it alone
         target_volume[e] = std::max(volume,Min_element_size); 
        }
       
      } //End of loop over elements
     
     // Tell everybody
     this->Nrefined=count_refined;
     this->Nunrefined=count_unrefined;
     
     return max_edge_ratio;
   }

   /// Max permitted element size
   double Max_element_size;
   
   /// Min permitted element size
   double Min_element_size;
   
   /// Max edge ratio before remesh gets triggered
   double Max_permitted_edge_ratio;

 };




}

#endif