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// Copyright (C) 2006-2023 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;
 
    /// Max. allowed discrepancy in element integrity check
    static double Max_integrity_tolerance;
 
    /// 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);
 
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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:
    /// 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;
 
    /// 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:
    /// Assign the local equation numbers for hanging node variables
    void assign_hanging_local_eqn_numbers(const bool& store_local_dof_pt);
 
    /// 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:
    /// 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&) = delete;
 
    /// Broken assignment operator
    void operator=(const RefineableElement&) = delete;
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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);
      }
    }
 
 
    /// 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];
    }
 
    /// 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;
    }
 
    /// 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;
 
    /// 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();
        }
      }
    }
 
    /// 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;
    }
 
    /// 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;
 
    /// 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);
    }
 
    /// 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;
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
 
    /// 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();
    }
 
    /// 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();
    }
 
    /// 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);
    }
 
    /// 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;
 
    /// Max. allowed discrepancy in element integrity check
    static double& max_integrity_tolerance()
    {
      return Max_integrity_tolerance;
    }
 
    /// 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);
 
 
    /// 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);
    }
 
    /// 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);
      }
    }
 
    /// 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)
    {
    }
 
    /// Pre-build the element
    virtual void pre_build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt) {}
 
    /// Further build: e.g. deal with interpolation of internal values
    virtual void further_build() {}
 
    /// 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) {}
 
    /// 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() {}
 
    /// 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);
 
 
    /// 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();
    }
 
    /// 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:
    /// 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&) = delete;
 
    /// Broken assignment operator
    void operator=(const PRefineableElement&) = delete;
 
    /// 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);
    }
 
 
    /// 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);
    }
 
    /// 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:
    /// 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;
 
 
    /// 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:
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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() {}
 
    /// 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);
    }
 
    /// 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;
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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];
    }
 
    /// 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
  {
  };
 
 
  /// ///////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////
 
 
} // namespace oomph
 
#endif