refineable_elements.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//====================================================================
//Header file for classes that define refineable element objects

//Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_REFINEABLE_ELEMENTS_HEADER
#define OOMPH_REFINEABLE_ELEMENTS_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

#include "elements.h"
#include "tree.h"

namespace oomph
{

class Mesh;

//=======================================================================
/// RefineableElements are FiniteElements that may be subdivided into 
/// children to provide a better local approximation to the solution. 
/// After non-uniform refinement adjacent elements need not necessarily have
/// nodes in common. A node that does not have a counterpart in its
/// neighbouring element is known as a hanging node and its position and
/// any data that it stores must be constrained to ensure inter-element
/// continuity. 
///
/// Generic data and function interfaces associated with refinement 
/// are defined in this class.
/// 
/// Additional data includes:
/// - a pointer to a general Tree object that is used to track the
///   refinement history,
/// - a refinement level (not necessarily the same as the level in the tree!),
/// - a flag indicating whether the element should be refined,
/// - a flag indicating whether the element should be de-refined,
/// - a global element number for plotting/validation  purposes,
/// - storage for local equation numbers associated with hanging nodes.
///
/// Additional functions perform the following generic tasks:
/// - provide access to  additional data,
/// - setup local equation numbering for data associated with hanging nodes,
/// - generic finite-difference calculation of contributions to the 
///   elemental jacobian from nodal data to include hanging nodes,
/// - split of the element into its sons,
/// - select and deselect the element for refinement,
/// - select and deselect the sons of the element for de-refinement (merging),
///
/// In addition, there are a number of interfaces that specify element-specific
/// tasks. These should be overloaded in RefineableElements of particular
/// geometric types and perform the following tasks:
/// - return a pointer to the root and father elements in the tree structure,
/// - define the number of sons into which the element is divided,
/// - build the element: construct nodes, assign their positions, 
///   values and any boundary conditions,
/// - recreate the element from its sons if they are merged,
/// - deactivate the element, perform any operations that are 
///   required when the element is still in the tree, but no longer active 
/// - set the number and provide access to the values interpolated
///   by the nodes,
/// - setup the hanging nodes
///  
/// In mixed element different sets of nodes are used to interpolate different
/// unknowns. Interfaces are provided for functions that can be used to find
/// the position of the nodes that interpolate the different unknowns. These
/// functions are used to setup hanging node information automatically in
/// particular elements, e.g. Taylor Hood Navier--Stokes. 
/// The default implementation assumes that the elements are isoparameteric.
///
//======================================================================
class RefineableElement : public virtual FiniteElement
{
  protected:

 /// A pointer to a general tree object
 Tree* Tree_pt;

 /// Refinement level
 unsigned Refine_level;

 /// Flag for refinement
 bool To_be_refined;
 
 /// Flag to indicate suppression of any refinement
 bool Refinement_is_enabled;
 
 /// Flag for unrefinement
 bool Sons_to_be_unrefined;

 /// Global element number -- for plotting/validation purposes
 long Number;
  
 /// \short Max. allowed discrepancy in element integrity check
 static double Max_integrity_tolerance;

 /// \short Static helper function that is used to check that the value_id
 /// is in range
 static void check_value_id(const int &n_continuously_interpolated_values,
                            const int &value_id);


  /// \short Assemble the jacobian matrix for the mapping from local
 /// to Eulerian coordinates, given the derivatives of the shape function
 /// w.r.t the local coordinates.
 /// Overload the standard version to use the hanging information for
 /// the Eulerian coordinates.
 void assemble_local_to_eulerian_jacobian(
  const DShape &dpsids, DenseMatrix<double> &jacobian) const;
 
 /// \short Assemble the the "jacobian" matrix of second derivatives of the
 /// mapping from local to Eulerian coordinates, given
 /// the second derivatives of the shape functions w.r.t. local coordinates.
 /// Overload the standard version to use the hanging information for
 /// the Eulerian coordinates.
 void assemble_local_to_eulerian_jacobian2(
  const DShape &d2psids, DenseMatrix<double> &jacobian2) const;
 
 /// \short Assemble the covariant Eulerian base vectors, assuming that
 /// the derivatives of the shape functions with respect to the local 
 /// coordinates have already been constructed.
 /// Overload the standard version to account for hanging nodes.
 void assemble_eulerian_base_vectors(
  const DShape &dpsids, DenseMatrix<double> &interpolated_G) const;
 
 /// \short Calculate the mapping from local to Eulerian coordinates given 
 /// the derivatives of the shape functions w.r.t the local coordinates.
 /// assuming that the coordinates are aligned in the direction of the local 
 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
 /// This funciton returns the determinant of the jacobian, the jacobian 
 /// and the inverse jacobian. Overload the standard version to take
 /// hanging info into account.
 double local_to_eulerian_mapping_diagonal(
  const DShape &dpsids,DenseMatrix<double> &jacobian,
  DenseMatrix<double> &inverse_jacobian) const;

  private:

 /// \short Storage for local equation numbers of hanging node variables
 /// (values stored at master nodes). It is
 /// essential that these are indexed by a Node pointer because the Node 
 /// may be internal or external to the element.
 /// local equation number = Local_hang_eqn(master_node_pt,ival)
  std::map<Node*,int> *Local_hang_eqn;

 /// \short Lookup scheme for unique number associated with any of the nodes
 /// that actively control the shape of the element (i.e. they are either
 /// non-hanging nodes of this element or master nodes of hanging nodes.
 std::map<Node*,unsigned> Shape_controlling_node_lookup;

  protected:
 

 /// \short Assign the local equation numbers for hanging node variables
 void assign_hanging_local_eqn_numbers(const bool &store_local_dof_pt);

 /// \short Calculate the contributions to the jacobian from the nodal
 /// degrees of freedom using finite differences.
 /// This version is overloaded to take hanging node information into
 /// account
 virtual void fill_in_jacobian_from_nodal_by_fd(Vector<double> &residuals,
                                            DenseMatrix<double> &jacobian);

  public:

 /// \short Constructor, calls the FiniteElement constructor and initialises
 /// the member data
 RefineableElement() : FiniteElement(), Tree_pt(0), Refine_level(0),
  To_be_refined(false), Refinement_is_enabled(true),
  Sons_to_be_unrefined(false), Number(-1),
  Local_hang_eqn(0) {}

 /// Destructor, delete the allocated storage for the hanging equations
 // (The body is now in the cc file to keep the xlC compiler happy under AIX)
 virtual ~RefineableElement();

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

  /// Access function: Pointer to quadtree representation of this element
 Tree* tree_pt() {return Tree_pt;}

 /// Set pointer to quadtree representation of this element
 void set_tree_pt(Tree* my_tree_pt) {Tree_pt=my_tree_pt;}

 /// \short Set the number of sons that can be constructed by the element
 /// The default is none
 virtual unsigned required_nsons() const {return 0;}


 /// Flag to indicate suppression of any refinement
 bool refinement_is_enabled(){return Refinement_is_enabled;}

 /// Suppress of any refinement for this element 
 void disable_refinement(){Refinement_is_enabled=false;}

 /// Emnable refinement for this element 
 void enable_refinement(){Refinement_is_enabled=true;}

 /// \short Split the element into the  number of sons to be
 /// constructed and return a 
 /// vector of pointers to the sons. Elements are allocated, but they are
 /// not given any properties. The refinement level of the sons is one
 /// higher than that of the father elemern.
 template<class ELEMENT>
  void split(Vector<ELEMENT*>& son_pt) const
  {
   // Increase refinement level
   int son_refine_level=Refine_level+1;
   
   //How many sons are to be constructed
   unsigned n_sons = required_nsons();
   //Resize the son pointer
   son_pt.resize(n_sons);

   //Loop over the sons and construct
   for(unsigned i=0;i<n_sons;i++)
    {
     son_pt[i]=new ELEMENT;
     //Set the refinement level of the newly constructed son.
     son_pt[i]->set_refinement_level(son_refine_level);
    }
  }


 /// \short Access function that returns the local equation number for the
 /// hanging node variables (values stored at master nodes). The local
 /// equation number corresponds to the i-th unknown stored at the node
 /// addressed by node_pt
 inline int local_hang_eqn(Node* const &node_pt, const unsigned &i)
  {
#ifdef RANGE_CHECKING
   if(i > ncont_interpolated_values())
    {
     std::ostringstream error_message;
     error_message << "Range Error: Value " << i
                   << " is not in the range (0,"
                   << ncont_interpolated_values()-1 << ")";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   return Local_hang_eqn[i][node_pt];
  }

 /// \short Interface to function that builds the element: i.e.  construct
 /// the nodes, assign their positions, apply boundary conditions, etc. The
 /// required procedures depend on the geometrical type of the element and
 /// must be implemented in specific refineable elements. Any new nodes
 /// created during the build process are returned in the vector
 /// new_node_pt.
 virtual void build(Mesh* &mesh_pt, Vector<Node*> &new_node_pt,
                    bool &was_already_built, std::ofstream &new_nodes_file)=0;

 /// Set the refinement level
 void set_refinement_level(const int& refine_level) 
  {Refine_level=refine_level;}
 
 /// Return the Refinement level
 unsigned refinement_level() const {return Refine_level;}
 
 /// Select the element for refinement
 void select_for_refinement() {To_be_refined=true;}

 /// Deselect the element for refinement
 void deselect_for_refinement() {To_be_refined=false;}
 
 /// Unrefinement will be performed by merging the four sons of this element
 void select_sons_for_unrefinement() {Sons_to_be_unrefined=true;}

 /// No unrefinement will be performed by merging the four sons of this element
 void deselect_sons_for_unrefinement() {Sons_to_be_unrefined=false;}

 /// Has the element been selected for refinement?
 bool to_be_refined() {return To_be_refined;}

 /// Has the element been selected for unrefinement?
 bool sons_to_be_unrefined() {return Sons_to_be_unrefined;}

 /// \short Rebuild the element, e.g. set internal values in line with
 /// those of the sons that have now merged
 virtual void rebuild_from_sons(Mesh* &mesh_pt)=0;

 /// \short Unbuild the element, i.e. mark the nodes that were created
 /// during its creation for possible deletion
 virtual void unbuild()
  {
   //Get pointer to father element
   RefineableElement* father_pt = father_element_pt();
   //If there is no father, nothing to do
   if(father_pt==0) {return;}

   //Loop over all the nodes
   unsigned n_node = this->nnode();
   for(unsigned n=0;n<n_node;n++)
    {
     //If any node in this element is in the father, it can't be deleted
     if(father_pt->get_node_number(this->node_pt(n)) >= 0)
      {
       node_pt(n)->set_non_obsolete();
      }
    }
  }

 /// \short Final operations that must be performed when the element is no
 /// longer active in the mesh, but still resident in the QuadTree.
 virtual void deactivate_element();
 
 /// Return true if all the nodes have been built, false if not
 virtual bool nodes_built() {return node_pt(0)!=0;}

 /// Element number (for debugging/plotting)
 long number() const {return Number;}
 
 /// Set element number (for debugging/plotting)
 void set_number(const long& mynumber) {Number=mynumber;}

 /// \short Number of continuously interpolated values. Note: We assume
 /// that they are located at the beginning of the value_pt Vector!
 /// (Used for interpolation to son elements, for integrity check
 /// and post-processing -- we can only expect
 /// the continously interpolated values to be continous across
 /// element boundaries).
 virtual unsigned ncont_interpolated_values() const=0;

 /// \short Get all continously interpolated function values in this 
 /// element as a Vector. Note: Vector sets is own size to ensure that
 /// that this function can be used in black-box fashion.
 virtual void get_interpolated_values(const Vector<double>&s, 
                                      Vector<double>& values)
 {
  get_interpolated_values(0, s, values);
 }
 
 /// \short Get all continously interpolated function values at previous
 /// timestep in this element as a Vector. (t=0: present; t>0:
 /// prev. timestep) Note: Vector sets is own size to ensure that that this
 /// function can be used in black-box fashion.
 virtual void get_interpolated_values(const unsigned& t,
                                      const Vector<double>&s, 
                                      Vector<double>& values)=0;

 /// \short In mixed elements, different sets of nodes are used to interpolate
 /// different unknowns. This function returns the n-th node that interpolates
 /// the value_id-th unknown. Default implementation is that all
 /// variables use the positional nodes, i.e. isoparametric elements. Note
 /// that any overloaded versions of this function MUST provide a set
 /// of nodes for the position, which always has the value_id -1.
 virtual Node* interpolating_node_pt(const unsigned &n,
                                     const int &value_id)
  
  {return node_pt(n);}

 /// \short Return the local one dimensional fraction of the n1d-th node
 /// in the direction of the local coordinate s[i] that is used to interpolate
 /// the value_id-th continuously interpolated unknown. Default assumes 
 /// isoparametric interpolation for all unknowns
 virtual double local_one_d_fraction_of_interpolating_node(
  const unsigned &n1d, const unsigned &i, const int &value_id)
  {return local_one_d_fraction_of_node(n1d,i);}

 /// \short Return a pointer to the node that interpolates the value-id-th
 /// unknown at local coordinate s. If there is not a node at that position,
 /// then return 0.
 virtual Node* 
  get_interpolating_node_at_local_coordinate(const Vector<double> &s,
                                             const int &value_id)
                                             
  {return get_node_at_local_coordinate(s);}


 /// \short Return the number of nodes that are used to interpolate the
 /// value_id-th unknown. Default is to assume isoparametric elements.
 virtual unsigned ninterpolating_node(const int &value_id) {return nnode();}

 /// \short Return the number of nodes in a one_d direction that are 
 /// used to interpolate the value_id-th unknown. Default is to assume
 /// an isoparametric mapping.
 virtual unsigned ninterpolating_node_1d(const int &value_id) 
  {return nnode_1d();}

 /// \short Return the basis functions that are used to interpolate
 /// the value_id-th unknown. By default assume isoparameteric interpolation
 virtual void interpolating_basis(const Vector<double> &s,
                                  Shape &psi,
                                  const int &value_id) const
  {shape(s,psi);}

 /// \short Check the integrity of the element: Continuity of positions
 /// values, etc. Essentially, check that the approximation of the functions
 /// is consistent when viewed from both sides of the element boundaries
 /// Must be overloaded for each different geometric element
 virtual void check_integrity(double &max_error)=0;

 /// \short Max. allowed discrepancy in element integrity check
 static double& max_integrity_tolerance()
  { return Max_integrity_tolerance;}

  /// \short The purpose of this function is to identify all possible
 /// Data that can affect the fields interpolated by the FiniteElement.
 /// This must be overloaded to include data from any hanging nodes 
 /// correctly
 void identify_field_data_for_interactions(
  std::set<std::pair<Data*,unsigned> > &paired_field_data);


 /// \short Overload the function that assigns local equation numbers
 /// for the Data stored at the nodes so that hanging data is taken 
 /// into account
 inline void assign_nodal_local_eqn_numbers(const bool &store_local_dof_pt)
  {
   FiniteElement::assign_nodal_local_eqn_numbers(store_local_dof_pt);
   assign_hanging_local_eqn_numbers(store_local_dof_pt);
  }

 /// \short Pointer to the root element in refinement hierarchy (must be 
 /// implemented in specific elements that do refinement via
 /// tree-like refinement structure. Here we provide a default
 /// implementation that is appropriate for cases where tree-like
 /// refinement doesn't exist or if the element doesn't have 
 /// root in that tree (i.e. if it's a root itself): We return
 /// "this". 
 virtual RefineableElement* root_element_pt()
  {
   //If there is no tree -- the element is its own root
   if(Tree_pt==0) {return this;}
   //Otherwise it's the tree's root object
   else {return Tree_pt->root_pt()->object_pt();}
  }

 /// Return a pointer to the father element.
 virtual RefineableElement* father_element_pt() const 
  {
   //If we have no tree, we have no father
   if(Tree_pt==0) {return 0;}
   else
    {
     //Otherwise get the father of the tree
     Tree* father_pt = Tree_pt->father_pt();
     //If the tree has no father then return null, no father
     if(father_pt==0) {return 0;}
     else {return father_pt->object_pt();}
    }
  }

 /// Return a pointer to the "father" element at the specified refinement level
 void get_father_at_refinement_level(unsigned& refinement_level, 
                                     RefineableElement*& father_at_reflevel_pt)
  {
   // Get the father in the tree (it shouldn't try to get a null Tree...)
   Tree* father_pt = Tree_pt->father_pt();
   // Get the refineable element associated with this father
   RefineableElement* father_el_pt=dynamic_cast<RefineableElement*>
    (father_pt->object_pt());
   // Get the refinement level
   unsigned level=father_el_pt->refinement_level();
   // If the level matches the required one then return, if not call again
   if (level==refinement_level) 
    {
     father_at_reflevel_pt=father_el_pt;
    }
   else
    {
     // Recursive call
     father_el_pt->get_father_at_refinement_level(refinement_level,
                                                  father_at_reflevel_pt);
    }
  }

 /// \short Initial setup of the element: e.g. set the appropriate internal
 /// p-order. If an adopted father is specified, information from this is
 /// used instead of using the father found from the tree.
 virtual void initial_setup(Tree* const &adopted_father_pt=0, const unsigned &initial_p_order=0) {}

 /// \short Pre-build the element
 virtual void pre_build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt) {}

 /// \short Further build: e.g. deal with interpolation of internal values
 virtual void further_build() {}
 
 /// \short Mark up any hanging nodes that arise as a result of non-uniform
 /// refinement. Any hanging nodes will be documented in files addressed by
 /// the streams in the vector output_stream, if the streams are open.
 virtual void setup_hanging_nodes(Vector<std::ofstream*> &output_stream) { }

 /// \short Perform additional hanging node procedures for variables
 /// that are not interpolated by all nodes (e.g. lower order interpolations
 /// for the pressure in Taylor Hood). 
 virtual void further_setup_hanging_nodes() { }

 /// \short Compute derivatives of elemental residual vector with respect
 /// to nodal coordinates. Default implementation by FD can be overwritten
 /// for specific elements. 
 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
 /// This version is overloaded from the version in FiniteElement
 /// and takes hanging nodes into account -- j in the above loop 
 /// loops over all the nodes that actively control the
 /// shape of the element (i.e. they are non-hanging or master nodes of 
 /// hanging nodes in this element).
 void get_dresidual_dnodal_coordinates(RankThreeTensor<double>&
                                       dresidual_dnodal_coordinates);
 

 /// \short Number of shape-controlling nodes = the number
 /// of non-hanging nodes plus the number of master nodes associated
 /// with hanging nodes.
 unsigned nshape_controlling_nodes()
  {
   return Shape_controlling_node_lookup.size();
  }

 /// \short Return lookup scheme for unique number associated
 /// with any of the nodes that actively control the shape of the
 /// element (i.e. they are either non-hanging nodes of this element
 /// or master nodes of hanging nodes.
 std::map<Node*,unsigned> shape_controlling_node_lookup()
  {
   return Shape_controlling_node_lookup;
  }


};
 

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




//======================================================================
/// p-refineable version of RefineableElement
//======================================================================
class PRefineableElement : public virtual RefineableElement
{
protected:
 
 /// The polynomial expansion order of the elemental basis functions
 unsigned P_order;

 /// Flag for p-refinement
 bool To_be_p_refined;
 
 /// Flag to indicate suppression of any refinement
 bool P_refinement_is_enabled;
 
 /// Flag for unrefinement
 bool To_be_p_unrefined;
 
public:

 /// \short Constructor, calls the RefineableElement constructor
 PRefineableElement() : RefineableElement(), P_order(2),
  To_be_p_refined(false), P_refinement_is_enabled(true),
  To_be_p_unrefined(false) {}

 /// Destructor, empty
 virtual ~PRefineableElement() {}

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


 /// Flag to indicate suppression of any refinement
 bool p_refinement_is_enabled(){return P_refinement_is_enabled;}

 /// Suppress of any refinement for this element 
 void disable_p_refinement(){P_refinement_is_enabled=false;}

 /// Emnable refinement for this element 
 void enable_p_refinement(){P_refinement_is_enabled=true;}
 
 /// Access function to P_order
 unsigned &p_order() {return P_order;}
 
 /// Access function to P_order (const version)
 unsigned p_order() const {return P_order;}
 
 /// Get the initial P_order
 /// This is required so that elements which are constructed with a
 /// higher p-order initially are not un-refined past this level (e.g.
 /// in fluid problems where elements initially use quadratic velocity
 /// and linear pressure)
 /// Virtual because this needs to be set in templated derived classes
 virtual unsigned initial_p_order() const=0;
 
 /// Select the element for p-refinement
 void select_for_p_refinement() {To_be_p_refined=true;}
 
 /// Deselect the element for p-refinement
 void deselect_for_p_refinement() {To_be_p_refined=false;}
 
 /// Select the element for p-unrefinement
 void select_for_p_unrefinement() {To_be_p_unrefined=true;}
 
 /// Deselect the element for p-unrefinement
 void deselect_for_p_unrefinement() {To_be_p_unrefined=false;}

 /// Has the element been selected for refinement?
 bool to_be_p_refined() {return To_be_p_refined;}
 
 /// Has the element been selected for p-unrefinement?
 bool to_be_p_unrefined() {return To_be_p_unrefined;}
 
 /// p-refine the element
 virtual void p_refine(const int &inc,
                       Mesh* const &mesh_pt,
                       GeneralisedElement* const &clone_pt)=0;
 
 // Overload the nodes_built function to check every node
 bool nodes_built()
  {
   // Must check that EVERY node exists
   unsigned n_node = this->nnode();
   for(unsigned n=0; n<n_node; n++)
    {
     if(this->node_pt(n)==0) {return false;}
    }
   // If we get here then all the nodes are built
   return true;
  }

};

 

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





//=======================================================================
/// A base class for elements that can have hanging nodes
/// but are not refineable as such. This class is usually used as a
/// base class for FaceElements that are attached to refineable
/// bulk elements (and stripped out before adapting the bulk
/// mesh, so they don't participate in the refimenent process
/// itself). We therefore simply break the pure virtual functions
/// that don't make any sense for such elements
//======================================================================
 class NonRefineableElementWithHangingNodes : public virtual RefineableElement
  {
   
    public:
   
   
   /// Broken build function -- shouldn't really be needed
   void build(Mesh* &mesh_pt, Vector<Node*> &new_node_pt,
              bool &was_already_built, std::ofstream &new_nodes_file)
    {
     std::ostringstream error_message;
     error_message << "This function is broken as it's only needed/used \n"
                   << "during \"proper\" refinement\n";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   
   
   
   /// \short Broken function -- this shouldn't really be needed.
   void get_interpolated_values(const Vector<double>&s, 
                                Vector<double>& values)
    {
     std::ostringstream error_message;
     error_message << "This function is broken as it's only needed/used \n"
                   << "during \"proper\" refinement\n";
     throw OomphLibError(
      error_message.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
   
   /// \short Broken function -- this shouldn't really be needed.
   virtual void get_interpolated_values(const unsigned& t,
                                        const Vector<double>&s, 
                                        Vector<double>& values)
    {
     std::ostringstream error_message;
     error_message << "This function is broken as it's only needed/used \n"
                   << "during \"proper\" refinement\n";
     throw OomphLibError(
      error_message.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
   
   
   /// Broken function -- this shouldn't really be needed.
   void check_integrity(double &max_error)
    {    
     std::ostringstream error_message;
     error_message << "This function is broken as it's only needed/used \n"
                   << "during \"proper\" refinement\n";
     throw OomphLibError(
      error_message.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
   
   /// Broken function -- this shouldn't really be needed.
   void rebuild_from_sons(Mesh* &mesh_pt)
    {
     std::ostringstream error_message;
     error_message << "This function is broken as it's only needed/used \n"
                   << "during \"proper\" refinement\n";
     throw OomphLibError(
      error_message.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
   
   
   
  };


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




//=======================================================================
/// RefineableSolidElements are SolidFiniteElements that may 
/// be subdivided into children to provide a better local approximation 
/// to the solution. The distinction is required to keep a clean 
/// separation between problems that alter nodal positions and others.
/// A number of procedures are generic and are included in this class.
//======================================================================
class RefineableSolidElement : public virtual RefineableElement,
                               public virtual SolidFiniteElement
{

  private:

 /// \short Storage for local equation numbers of 
 /// hanging node variables associated with nodal positions.
 /// local position equation number = 
 /// Local_position_hang_eqn(master_node_pt,ival)
 std::map<Node*,DenseMatrix<int> > Local_position_hang_eqn; 


 /// \short Assign local equation numbers to the hanging values associated
 /// with positions or additional solid values.
 void assign_solid_hanging_local_eqn_numbers(const bool &store_local_dof_pt);


 protected:

 /// \short Flag deciding if the Lagrangian coordinates of newly-created interior 
 /// SolidNodes are to be determined by the father element's undeformed 
 /// MacroElement representation (if it has one). Default: False as it means that, 
 /// following a refinement an element is no longer in equilbrium (as the Lagrangian
 /// coordinate is determined differently from the Eulerian one). However, basing
 /// the Lagrangian coordinates on the undeformed MacroElement can be
 /// more stable numerically and for steady problems it's not a big deal
 /// either way as the difference between the two formulations only matters
 /// at finite resolution so we have no right to say that one is "more correct"
 /// than the other...
 bool Use_undeformed_macro_element_for_new_lagrangian_coords;
 
 /// \short Assemble the jacobian matrix for the mapping from local
 /// to lagrangian coordinates, given the derivatives of the shape function
 /// Overload the standard version to use the hanging information for
 /// the lagrangian coordinates.
void assemble_local_to_lagrangian_jacobian(
  const DShape &dpsids, DenseMatrix<double> &jacobian) const;

 /// \short Assemble the the "jacobian" matrix of second derivatives, given
 /// the second derivatives of the shape functions w.r.t. local coordinates
 /// Overload the standard version to use the hanging information for
 /// the lagrangian coordinates.
 void assemble_local_to_lagrangian_jacobian2(
  const DShape &d2psids, DenseMatrix<double> &jacobian2) const;
 
 /// \short Calculate the mapping from local to Lagrangian coordinates given 
 /// the derivatives of the shape functions w.r.t the local coorindates.
 /// assuming that the coordinates are aligned in the direction of the local 
 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
 /// This function returns the determinant of the jacobian, the jacobian 
 /// and the inverse jacobian.
 double local_to_lagrangian_mapping_diagonal(
  const DShape &dpsids,DenseMatrix<double> &jacobian,
  DenseMatrix<double> &inverse_jacobian) const;

  public:
 
 /// Constructor
 RefineableSolidElement() : RefineableElement(), SolidFiniteElement(),
  Use_undeformed_macro_element_for_new_lagrangian_coords(false){}

 /// Virtual Destructor, delete any allocated storage
 virtual ~RefineableSolidElement() { }
 
 /// \short Overload the local equation numbers for Data stored as part
 /// of solid nodes to include hanging node data
 void assign_solid_local_eqn_numbers(const bool &store_local_dof_pt)
  {
   SolidFiniteElement::assign_solid_local_eqn_numbers(store_local_dof_pt);
   assign_solid_hanging_local_eqn_numbers(store_local_dof_pt);
  }

 ///\short The number of geometric data affecting a 
 /// RefineableSolidFiniteElement is the positional Data of all
 /// non-hanging nodes plus the geometric Data of all distinct 
 /// master nodes. Recomputed on the fly.
 unsigned ngeom_data() const;

 /// \short Return pointer to the j-th Data item that the object's 
 /// shape depends on: Positional data of non-hanging nodes and
 /// positional data of master nodes. Recomputed on the fly.
  Data* geom_data_pt(const unsigned& j);

 /// \short Specify Data that affects the geometry of the element
 /// by adding the position Data to the set that's passed in.
 /// (This functionality is required in FSI problems; set is used to
 /// avoid double counting). Refineable version includes hanging nodes
  void identify_geometric_data(std::set<Data*> &geometric_data_pt);

 /// \short Compute element residual Vector and element Jacobian matrix
 /// corresponding to the solid positions. Overloaded version to take
 /// the hanging nodes into account
 void fill_in_jacobian_from_solid_position_by_fd(Vector<double> & residuals,
                                             DenseMatrix<double> &jacobian);

 /// \short Return the flag deciding if the Lagrangian coordinates of 
 /// newly-created interior SolidNodes are to be determined by the father 
 /// element's undeformed MacroElement representation (if it has one). 
 /// Default: False as it means that, following a refinement an element 
 /// is no longer in equilbrium (as the Lagrangian coordinate is 
 /// determined differently from the Eulerian one). However, basing
 /// the Lagrangian coordinates on the undeformed MacroElement can be
 /// more stable numerically and for steady problems it's not a big deal
 /// either way as the difference between the two formulations only matters
 /// at finite resolution so we have no right to say that one is "more correct"
 /// than the other...
 bool is_undeformed_macro_element_used_for_new_lagrangian_coords() const
 {return  Use_undeformed_macro_element_for_new_lagrangian_coords;}

 /// \short Set the flag deciding if the Lagrangian coordinates of 
 /// newly-created interior SolidNodes are to be determined by the father 
 /// element's undeformed MacroElement representation (if it has one). 
 void enable_use_of_undeformed_macro_element_for_new_lagrangian_coords()
  {Use_undeformed_macro_element_for_new_lagrangian_coords=true;}

 /// \short Unset the flag deciding if the Lagrangian coordinates of 
 /// newly-created interior SolidNodes are to be determined by the father 
 /// element's undeformed MacroElement representation (if it has one). 
 void disable_use_of_undeformed_macro_element_for_new_lagrangian_coords()
  {Use_undeformed_macro_element_for_new_lagrangian_coords=false;}

 /// \short Access the local equation number of of hanging node variables
 /// associated with nodal positions. The function returns a dense
 /// matrix that contains all the local equation numbers corresponding to
 /// the positional degrees of freedom.
 DenseMatrix<int> &local_position_hang_eqn(Node* const &node_pt)
  {return Local_position_hang_eqn[node_pt];}

 /// \short Further build: Pass the father's 
 /// Use_undeformed_macro_element_for_new_lagrangian_coords
 /// flag down, then call the underlying RefineableElement's
 /// version.
 virtual void further_build()
  {
   Use_undeformed_macro_element_for_new_lagrangian_coords=
    dynamic_cast<RefineableSolidElement*>(father_element_pt())
    ->is_undeformed_macro_element_used_for_new_lagrangian_coords();
   
   RefineableElement::further_build();
  }

};


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





//=======================================================================
/// A base class for SolidElements that can have hanging nodes
/// but are not refineable as such. This class is usually used as a
/// base class for FaceElements that are attached to refineable
/// bulk elements (and stripped out before adapting the bulk
/// mesh, so they don't participate in the refimenent process
/// itself). We therefore simply break the pure virtual functions
/// that don't make any sense for such elements
//======================================================================
class NonRefineableSolidElementWithHangingNodes : 
public virtual NonRefineableElementWithHangingNodes,
public virtual RefineableSolidElement
{
  
};


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


}

#endif