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-2022 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 element objects
 
// Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_GENERIC_ELEMENTS_HEADER
#define OOMPH_GENERIC_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include <map>
#include <deque>
#include <string>
#include <list>
 
// oomph-lib includes
#include "integral.h"
#include "nodes.h"
#include "geom_objects.h"
 
namespace oomph
{
  // Forward definition for time
  class Time;
 
 
  //========================================================================
  //========================================================================
  class GeneralisedElement
  {
  private:
    unsigned long* Eqn_number;
 
    double** Dof_pt;
 
    Data** Data_pt;
 
    int** Data_local_eqn;
 
    unsigned Ndof;
 
    unsigned Ninternal_data;
 
    unsigned Nexternal_data;
 
    std::vector<bool> Data_fd;
 
  protected:
#ifdef OOMPH_HAS_MPI
 
    int Non_halo_proc_ID;
 
    bool Must_be_kept_as_halo;
 
#endif
 
    unsigned add_internal_data(Data* const& data_pt, const bool& fd = true);
 
 
    inline bool internal_data_fd(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      if (i >= Ninternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Internal data " << i
                      << " is not in the range (0," << Ninternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return the value
      return Data_fd[i];
    }
 
 
    inline void exclude_internal_data_fd(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      if (i >= Ninternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Internal data " << i
                      << " is not in the range (0," << Ninternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Set the value
      Data_fd[i] = false;
    }
 
    inline void include_internal_data_fd(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      if (i >= Ninternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Internal data " << i
                      << " is not in the range (0," << Ninternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Set the value
      Data_fd[i] = true;
    }
 
    void clear_global_eqn_numbers()
    {
      delete[] Eqn_number;
      Eqn_number = 0;
      delete[] Dof_pt;
      Dof_pt = 0;
      Ndof = 0;
    }
 
    void add_global_eqn_numbers(
      std::deque<unsigned long> const& global_eqn_numbers,
      std::deque<double*> const& global_dof_pt);
 
    static DenseMatrix<double> Dummy_matrix;
 
    static std::deque<double*> Dof_pt_deque;
 
    virtual void assign_internal_and_external_local_eqn_numbers(
      const bool& store_local_dof_pt);
 
    virtual inline void assign_all_generic_local_eqn_numbers(
      const bool& store_local_dof_pt)
    {
      this->assign_internal_and_external_local_eqn_numbers(store_local_dof_pt);
    }
 
    virtual void assign_additional_local_eqn_numbers() {}
 
    inline int internal_local_eqn(const unsigned& i, const unsigned& j) const
    {
#ifdef RANGE_CHECKING
      if (i >= Ninternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Internal data " << i
                      << " is not in the range (0," << Ninternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else
      {
        unsigned n_value = internal_data_pt(i)->nvalue();
        if (j >= n_value)
        {
          std::ostringstream error_message;
          error_message << "Range Error: value " << j << " at internal data "
                        << i << " is not in the range (0," << n_value - 1
                        << ")";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
#endif
      // Check that the data has been allocated
#ifdef PARANOID
      if (Data_local_eqn == 0)
      {
        throw OomphLibError(
          "Internal local equation numbers have not been allocated",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return the local equation number as stored in the Data_local_eqn array
      return Data_local_eqn[i][j];
    }
 
    inline int external_local_eqn(const unsigned& i, const unsigned& j)
    {
      // Check that the data has been allocated
#ifdef RANGE_CHECKING
      if (i >= Nexternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: External data " << i
                      << " is not in the range (0," << Nexternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else
      {
        unsigned n_value = external_data_pt(i)->nvalue();
        if (j >= n_value)
        {
          std::ostringstream error_message;
          error_message << "Range Error: value " << j << " at internal data "
                        << i << " is not in the range (0," << n_value - 1
                        << ")";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
#endif
#ifdef PARANOID
      if (Data_local_eqn == 0)
      {
        throw OomphLibError(
          "External local equation numbers have not been allocated",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return the  local equation number stored as the j-th value of the
      // i-th external data object.
      return Data_local_eqn[Ninternal_data + i][j];
    }
 
    virtual void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      std::string error_message =
        "Empty fill_in_contribution_to_residuals() has been called.\n";
      error_message +=
        "This function is called from the default implementations of\n";
      error_message += "get_residuals() and get_jacobian();\n";
      error_message +=
        "and must calculate the residuals vector without initialising any of ";
      error_message += "its entries.\n\n";
 
      error_message +=
        "If you do not wish to use these defaults, you must overload both\n";
      error_message += "get_residuals() and get_jacobian(), which must "
                       "initialise the entries\n";
      error_message += "of the residuals vector and jacobian matrix to zero.\n";
 
      error_message += "N.B. the default get_jacobian() function employs "
                       "finite differencing\n";
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    void fill_in_jacobian_from_internal_by_fd(Vector<double>& residuals,
                                              DenseMatrix<double>& jacobian,
                                              const bool& fd_all_data = false);
 
    void fill_in_jacobian_from_internal_by_fd(DenseMatrix<double>& jacobian,
                                              const bool& fd_all_data = false)
    {
      // Allocate storage for the residuals
      Vector<double> residuals(Ndof, 0.0);
      // Get the residuals for the entire element
      get_residuals(residuals);
      // Call the jacobian calculation
      fill_in_jacobian_from_internal_by_fd(residuals, jacobian, fd_all_data);
    }
 
 
    void fill_in_jacobian_from_external_by_fd(Vector<double>& residuals,
                                              DenseMatrix<double>& jacobian,
                                              const bool& fd_all_data = false);
 
 
    void fill_in_jacobian_from_external_by_fd(DenseMatrix<double>& jacobian,
                                              const bool& fd_all_data = false)
    {
      // Allocate storage for a residuals vector
      Vector<double> residuals(Ndof);
      // Get the residuals for the entire element
      get_residuals(residuals);
      // Call the jacobian calculation
      fill_in_jacobian_from_external_by_fd(residuals, jacobian, fd_all_data);
    }
 
    virtual inline void update_before_internal_fd() {}
 
    virtual inline void reset_after_internal_fd() {}
 
    virtual inline void update_in_internal_fd(const unsigned& i) {}
 
    virtual inline void reset_in_internal_fd(const unsigned& i)
    {
      update_in_internal_fd(i);
    }
 
 
    virtual inline void update_before_external_fd() {}
 
    virtual inline void reset_after_external_fd() {}
 
    virtual inline void update_in_external_fd(const unsigned& i) {}
 
    virtual inline void reset_in_external_fd(const unsigned& i)
    {
      update_in_external_fd(i);
    }
 
    virtual void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                                  DenseMatrix<double>& jacobian)
    {
      // Add the contribution to the residuals
      fill_in_contribution_to_residuals(residuals);
      // Allocate storage for the full residuals (residuals of entire element)
      Vector<double> full_residuals(Ndof);
      // Get the residuals for the entire element
      get_residuals(full_residuals);
      // Calculate the contributions from the internal dofs
      //(finite-difference the lot by default)
      fill_in_jacobian_from_internal_by_fd(full_residuals, jacobian, true);
      // Calculate the contributions from the external dofs
      //(finite-difference the lot by default)
      fill_in_jacobian_from_external_by_fd(full_residuals, jacobian, true);
    }
 
    virtual void fill_in_contribution_to_mass_matrix(
      Vector<double>& residuals, DenseMatrix<double>& mass_matrix);
 
    virtual void fill_in_contribution_to_jacobian_and_mass_matrix(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix);
 
    virtual void fill_in_contribution_to_dresiduals_dparameter(
      double* const& parameter_pt, Vector<double>& dres_dparam);
 
 
    virtual void fill_in_contribution_to_djacobian_dparameter(
      double* const& parameter_pt,
      Vector<double>& dres_dparam,
      DenseMatrix<double>& djac_dparam);
 
    virtual void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(
      double* const& parameter_pt,
      Vector<double>& dres_dparam,
      DenseMatrix<double>& djac_dparam,
      DenseMatrix<double>& dmass_matrix_dparam);
 
 
    virtual void fill_in_contribution_to_hessian_vector_products(
      Vector<double> const& Y,
      DenseMatrix<double> const& C,
      DenseMatrix<double>& product);
 
    virtual void fill_in_contribution_to_inner_products(
      Vector<std::pair<unsigned, unsigned>> const& history_index,
      Vector<double>& inner_product);
 
    virtual void fill_in_contribution_to_inner_product_vectors(
      Vector<unsigned> const& history_index,
      Vector<Vector<double>>& inner_product_vector);
 
  public:
    GeneralisedElement()
      : Eqn_number(0),
        Dof_pt(0),
        Data_pt(0),
        Data_local_eqn(0),
        Ndof(0),
        Ninternal_data(0),
        Nexternal_data(0)
#ifdef OOMPH_HAS_MPI
        ,
        Non_halo_proc_ID(-1),
        Must_be_kept_as_halo(false)
#endif
    {
    }
 
    virtual ~GeneralisedElement();
 
    GeneralisedElement(const GeneralisedElement&) = delete;
 
    void operator=(const GeneralisedElement&) = delete;
 
    Data*& internal_data_pt(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      if (i >= Ninternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Internal data " << i
                      << " is not in the range (0," << Ninternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Data_pt[i];
    }
 
    Data* const& internal_data_pt(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      if (i >= Ninternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Internal data " << i
                      << " is not in the range (0," << Ninternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Data_pt[i];
    }
 
 
    Data*& external_data_pt(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      if (i >= Nexternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: External data " << i
                      << " is not in the range (0," << Nexternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Data_pt[Ninternal_data + i];
    }
 
    Data* const& external_data_pt(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      if (i >= Nexternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: External data " << i
                      << " is not in the range (0," << Nexternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Data_pt[Ninternal_data + i];
    }
 
    static bool Suppress_warning_about_repeated_internal_data;
 
    static bool Suppress_warning_about_repeated_external_data;
 
    unsigned long eqn_number(const unsigned& ieqn_local) const
    {
#ifdef RANGE_CHECKING
      if (ieqn_local >= Ndof)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Equation number " << ieqn_local
                      << " is not in the range (0," << Ndof - 1 << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Eqn_number[ieqn_local];
    }
 
 
    int local_eqn_number(const unsigned long& ieqn_global) const
    {
      // Cache the number of degrees of freedom in the element
      const unsigned n_dof = this->Ndof;
      // Loop over the local equation numbers
      for (unsigned n = 0; n < n_dof; n++)
      {
        // If the global equation numbers match return
        if (ieqn_global == Eqn_number[n])
        {
          return n;
        }
      }
 
      // If we've got all the way to the end the number has not been found
      // return minus one.
      return -1;
    }
 
 
    unsigned add_external_data(Data* const& data_pt, const bool& fd = true);
 
    inline bool external_data_fd(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      if (i >= Nexternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Internal data " << i
                      << " is not in the range (0," << Nexternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return the value
      return Data_fd[Ninternal_data + i];
    }
 
 
    inline void exclude_external_data_fd(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      if (i >= Nexternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: External data " << i
                      << " is not in the range (0," << Nexternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Set the value
      Data_fd[Ninternal_data + i] = false;
    }
 
    inline void include_external_data_fd(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      if (i >= Nexternal_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error: External data " << i
                      << " is not in the range (0," << Nexternal_data - 1
                      << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Set the value
      Data_fd[Ninternal_data + i] = true;
    }
 
    void flush_external_data();
 
    void flush_external_data(Data* const& data_pt);
 
    unsigned ninternal_data() const
    {
      return Ninternal_data;
    }
 
    unsigned nexternal_data() const
    {
      return Nexternal_data;
    }
 
    unsigned ndof() const
    {
      return Ndof;
    }
 
    void dof_vector(const unsigned& t, Vector<double>& dof)
    {
      // Check that the internal storage has been set up
#ifdef PARANOID
      if (Dof_pt == 0)
      {
        std::stringstream error_stream;
        error_stream << "Internal dof array not set up in element.\n"
                     << "In order to set it up you must call\n"
                     << "   Problem::enable_store_local_dof_in_elements()\n"
                     << "before the call to Problem::assign_eqn_numbers()\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Resize the vector
      dof.resize(this->Ndof);
      // Loop over the vector and fill in the desired values
      for (unsigned i = 0; i < this->Ndof; ++i)
      {
        dof[i] = Dof_pt[i][t];
      }
    }
 
    void dof_pt_vector(Vector<double*>& dof_pt)
    {
      // Check that the internal storage has been set up
#ifdef PARANOID
      if (Dof_pt == 0)
      {
        std::stringstream error_stream;
        error_stream << "Internal dof array not set up in element.\n"
                     << "In order to set it up you must call\n"
                     << "   Problem::enable_store_local_dof_in_elements()\n"
                     << "before the call to Problem::assign_eqn_numbers()\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Resize the vector
      dof_pt.resize(this->Ndof);
      // Loop over the vector and fill in the desired values
      for (unsigned i = 0; i < this->Ndof; ++i)
      {
        dof_pt[i] = Dof_pt[i];
      }
    }
 
 
    void set_internal_data_time_stepper(const unsigned& i,
                                        TimeStepper* const& time_stepper_pt,
                                        const bool& preserve_existing_data)
    {
      this->internal_data_pt(i)->set_time_stepper(time_stepper_pt,
                                                  preserve_existing_data);
    }
 
    void assign_internal_eqn_numbers(unsigned long& global_number,
                                     Vector<double*>& Dof_pt);
 
    void describe_dofs(std::ostream& out,
                       const std::string& current_string) const;
 
    virtual void describe_local_dofs(std::ostream& out,
                                     const std::string& current_string) const;
 
    void add_internal_value_pt_to_map(
      std::map<unsigned, double*>& map_of_value_pt);
 
#ifdef OOMPH_HAS_MPI
    void add_internal_data_values_to_vector(Vector<double>& vector_of_values);
 
    void read_internal_data_values_from_vector(
      const Vector<double>& vector_of_values, unsigned& index);
 
    void add_internal_eqn_numbers_to_vector(
      Vector<long>& vector_of_eqn_numbers);
 
    void read_internal_eqn_numbers_from_vector(
      const Vector<long>& vector_of_eqn_numbers, unsigned& index);
 
#endif
 
    virtual void assign_local_eqn_numbers(const bool& store_local_dof_pt);
 
    virtual void complete_setup_of_dependencies() {}
 
    virtual void get_residuals(Vector<double>& residuals)
    {
      // Zero the residuals vector
      residuals.initialise(0.0);
      // Add the elemental contribution to the residuals vector
      fill_in_contribution_to_residuals(residuals);
    }
 
    virtual void get_jacobian(Vector<double>& residuals,
                              DenseMatrix<double>& jacobian)
    {
      // Zero the residuals vector
      residuals.initialise(0.0);
      // Zero the jacobian matrix
      jacobian.initialise(0.0);
      // Add the elemental contribution to the residuals vector
      fill_in_contribution_to_jacobian(residuals, jacobian);
    }
 
    virtual void get_mass_matrix(Vector<double>& residuals,
                                 DenseMatrix<double>& mass_matrix)
    {
      // Zero the residuals vector
      residuals.initialise(0.0);
      // Zero the mass matrix
      mass_matrix.initialise(0.0);
      // Add the elemental contribution to the vector and matrix
      fill_in_contribution_to_mass_matrix(residuals, mass_matrix);
    }
 
    virtual void get_jacobian_and_mass_matrix(Vector<double>& residuals,
                                              DenseMatrix<double>& jacobian,
                                              DenseMatrix<double>& mass_matrix)
    {
      // Zero the residuals vector
      residuals.initialise(0.0);
      // Zero the jacobian matrix
      jacobian.initialise(0.0);
      // Zero the mass matrix
      mass_matrix.initialise(0.0);
      // Add the elemental contribution to the vector and matrix
      fill_in_contribution_to_jacobian_and_mass_matrix(
        residuals, jacobian, mass_matrix);
    }
 
 
    virtual void get_dresiduals_dparameter(double* const& parameter_pt,
                                           Vector<double>& dres_dparam)
    {
      // Zero the dres_dparam vector
      dres_dparam.initialise(0.0);
      // Add the elemental contribution to the residuals vector
      this->fill_in_contribution_to_dresiduals_dparameter(parameter_pt,
                                                          dres_dparam);
    }
 
    virtual void get_djacobian_dparameter(double* const& parameter_pt,
                                          Vector<double>& dres_dparam,
                                          DenseMatrix<double>& djac_dparam)
    {
      // Zero the residuals vector
      dres_dparam.initialise(0.0);
      // Zero the jacobian matrix
      djac_dparam.initialise(0.0);
      // Add the elemental contribution to the residuals vector
      this->fill_in_contribution_to_djacobian_dparameter(
        parameter_pt, dres_dparam, djac_dparam);
    }
 
    virtual void get_djacobian_and_dmass_matrix_dparameter(
      double* const& parameter_pt,
      Vector<double>& dres_dparam,
      DenseMatrix<double>& djac_dparam,
      DenseMatrix<double>& dmass_matrix_dparam)
    {
      // Zero the residuals derivative vector
      dres_dparam.initialise(0.0);
      // Zero the jacobian derivative  matrix
      djac_dparam.initialise(0.0);
      // Zero the mass matrix derivative
      dmass_matrix_dparam.initialise(0.0);
      // Add the elemental contribution to the residuals vector and matrices
      this->fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(
        parameter_pt, dres_dparam, djac_dparam, dmass_matrix_dparam);
    }
 
 
    virtual void get_hessian_vector_products(Vector<double> const& Y,
                                             DenseMatrix<double> const& C,
                                             DenseMatrix<double>& product)
    {
      // Initialise the product to zero
      product.initialise(0.0);
      this->fill_in_contribution_to_hessian_vector_products(Y, C, product);
    }
 
    virtual void get_inner_products(
      Vector<std::pair<unsigned, unsigned>> const& history_index,
      Vector<double>& inner_product)
    {
      inner_product.initialise(0.0);
      this->fill_in_contribution_to_inner_products(history_index,
                                                   inner_product);
    }
 
    virtual void get_inner_product_vectors(
      Vector<unsigned> const& history_index,
      Vector<Vector<double>>& inner_product_vector)
    {
      const unsigned n_inner_product = inner_product_vector.size();
      for (unsigned i = 0; i < n_inner_product; ++i)
      {
        inner_product_vector[i].initialise(0.0);
      }
      this->fill_in_contribution_to_inner_product_vectors(history_index,
                                                          inner_product_vector);
    }
 
 
    virtual unsigned self_test();
 
 
    virtual void compute_norm(Vector<double>& norm)
    {
      std::string error_message =
        "compute_norm(...) undefined for this element\n";
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void compute_norm(double& norm)
    {
      std::string error_message = "compute_norm undefined for this element \n";
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
#ifdef OOMPH_HAS_MPI
 
    void set_halo(const unsigned& non_halo_proc_ID)
    {
      Non_halo_proc_ID = non_halo_proc_ID;
    }
 
    void set_nonhalo()
    {
      Non_halo_proc_ID = -1;
    }
 
    bool is_halo() const
    {
      return (Non_halo_proc_ID != -1);
    }
 
    int non_halo_proc_ID()
    {
      return Non_halo_proc_ID;
    }
 
    void set_must_be_kept_as_halo()
    {
      Must_be_kept_as_halo = true;
    }
 
    void unset_must_be_kept_as_halo()
    {
      Must_be_kept_as_halo = false;
    }
 
    bool must_be_kept_as_halo() const
    {
      return Must_be_kept_as_halo;
    }
 
#endif
 
    static double Default_fd_jacobian_step;
 
    virtual unsigned ndof_types() const
    {
      // error message stream
      std::ostringstream error_message;
      error_message << "ndof_types() const has not been implemented for this \n"
                    << "element\n"
                    << std::endl;
      // throw error
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
    {
      // error message stream
      std::ostringstream error_message;
      error_message << "get_dof_numbers_for_unknowns() const has not been \n"
                    << " implemented for this element\n"
                    << std::endl;
      // throw error
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  };
 
  namespace ElementGeometry
  {
    enum ElementGeometry
    {
      Q,
      T
    };
  }
 
  // Forward class definitions that are used in FiniteElements
  class FaceElement;
  class MacroElement;
  class TimeStepper;
  class Shape;
  class DShape;
  class Integral;
 
  //===========================================================================
  //===========================================================================
  namespace Locate_zeta_helpers
  {
    extern double Newton_tolerance;
 
    extern unsigned Max_newton_iterations;
 
    extern double Radius_multiplier_for_fast_exit_from_locate_zeta;
 
    extern unsigned N_local_points;
 
  } // namespace Locate_zeta_helpers
 
 
  typedef void (*CoordinateMappingFctPt)(const Vector<double>& s,
                                         Vector<double>& s_bulk);
 
  typedef void (*BulkCoordinateDerivativesFctPt)(
    const Vector<double>& s,
    DenseMatrix<double>& ds_bulk_dsface,
    unsigned& interior_direction);
 
 
  //========================================================================
  //========================================================================
  class FiniteElement : public virtual GeneralisedElement, public GeomObject
  {
  private:
    Integral* Integral_pt;
 
    Node** Node_pt;
 
    int** Nodal_local_eqn;
 
    unsigned Nnode;
 
    unsigned Elemental_dimension;
 
    unsigned Nodal_dimension;
 
    unsigned Nnodal_position_type;
 
  protected:
    virtual void assemble_local_to_eulerian_jacobian(
      const DShape& dpsids, DenseMatrix<double>& jacobian) const;
 
    virtual void assemble_local_to_eulerian_jacobian2(
      const DShape& d2psids, DenseMatrix<double>& jacobian2) const;
 
    virtual void assemble_eulerian_base_vectors(
      const DShape& dpsids, DenseMatrix<double>& interpolated_G) const;
 
    static const unsigned Default_Initial_Nvalue;
 
    static const double Node_location_tolerance;
 
  public:
    void set_dimension(const unsigned& dim)
    {
      Elemental_dimension = dim;
      Nodal_dimension = dim;
    }
 
    void set_nodal_dimension(const unsigned& nodal_dim)
    {
      Nodal_dimension = nodal_dim;
    }
 
    void set_nnodal_position_type(const unsigned& nposition_type)
    {
      Nnodal_position_type = nposition_type;
    }
 
 
    void set_n_node(const unsigned& n)
    {
      // This should only be done once, in a Node constructor
      //#ifdef PARANOID
      // if(Node_pt)
      // {
      //  OomphLibWarning(
      //   "You are resetting the number of nodes in an element -- Be careful",
      //   "FiniteElement::set_n_node()",
      //   OOMPH_EXCEPTION_LOCATION);
      // }
      //#endif
      // Delete any previous storage to avoid memory leaks
      // This will only happen in very special cases
      delete[] Node_pt;
      // Set the number of nodes
      Nnode = n;
      // Allocate the storage
      Node_pt = new Node*[n];
      // Initialise all the pointers to NULL
      for (unsigned i = 0; i < n; i++)
      {
        Node_pt[i] = 0;
      }
    }
 
    inline int nodal_local_eqn(const unsigned& n, const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      if (n >= Nnode)
      {
        std::ostringstream error_message;
        error_message << "Range Error: Node number " << n
                      << " is not in the range (0," << Nnode - 1 << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else
      {
        unsigned n_value = node_pt(n)->nvalue();
        if (i >= n_value)
        {
          std::ostringstream error_message;
          error_message << "Range Error: value " << i << " at node " << n
                        << " is not in the range (0," << n_value - 1 << ")";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
#endif
#ifdef PARANOID
      // Check that the equations have been allocated
      if (Nodal_local_eqn == 0)
      {
        throw OomphLibError(
          "Nodal local equation numbers have not been allocated",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Nodal_local_eqn[n][i];
    }
 
 
    double dJ_eulerian_at_knot(const unsigned& ipt,
                               Shape& psi,
                               DenseMatrix<double>& djacobian_dX) const;
 
  protected:
    static const unsigned N2deriv[];
 
    template<unsigned DIM>
    double invert_jacobian(const DenseMatrix<double>& jacobian,
                           DenseMatrix<double>& inverse_jacobian) const;
 
 
    virtual double invert_jacobian_mapping(
      const DenseMatrix<double>& jacobian,
      DenseMatrix<double>& inverse_jacobian) const;
 
 
    virtual double local_to_eulerian_mapping(
      const DShape& dpsids,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& inverse_jacobian) const
    {
      // Assemble the jacobian
      assemble_local_to_eulerian_jacobian(dpsids, jacobian);
      // Invert the jacobian (use the template-free interface)
      return invert_jacobian_mapping(jacobian, inverse_jacobian);
    }
 
 
    double local_to_eulerian_mapping(
      const DShape& dpsids, DenseMatrix<double>& inverse_jacobian) const
    {
      // Find the dimension of the element
      const unsigned el_dim = dim();
      // Assign memory for the jacobian
      DenseMatrix<double> jacobian(el_dim);
      // Calculate the jacobian and inverse
      return local_to_eulerian_mapping(dpsids, jacobian, inverse_jacobian);
    }
 
    virtual double local_to_eulerian_mapping_diagonal(
      const DShape& dpsids,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& inverse_jacobian) const;
 
    virtual void dJ_eulerian_dnodal_coordinates(
      const DenseMatrix<double>& jacobian,
      const DShape& dpsids,
      DenseMatrix<double>& djacobian_dX) const;
 
    template<unsigned DIM>
    void dJ_eulerian_dnodal_coordinates_templated_helper(
      const DenseMatrix<double>& jacobian,
      const DShape& dpsids,
      DenseMatrix<double>& djacobian_dX) const;
 
    virtual void d_dshape_eulerian_dnodal_coordinates(
      const double& det_jacobian,
      const DenseMatrix<double>& jacobian,
      const DenseMatrix<double>& djacobian_dX,
      const DenseMatrix<double>& inverse_jacobian,
      const DShape& dpsids,
      RankFourTensor<double>& d_dpsidx_dX) const;
 
    template<unsigned DIM>
    void d_dshape_eulerian_dnodal_coordinates_templated_helper(
      const double& det_jacobian,
      const DenseMatrix<double>& jacobian,
      const DenseMatrix<double>& djacobian_dX,
      const DenseMatrix<double>& inverse_jacobian,
      const DShape& dpsids,
      RankFourTensor<double>& d_dpsidx_dX) const;
 
    virtual void transform_derivatives(
      const DenseMatrix<double>& inverse_jacobian, DShape& dbasis) const;
 
    void transform_derivatives_diagonal(
      const DenseMatrix<double>& inverse_jacobian, DShape& dbasis) const;
 
    virtual void transform_second_derivatives(
      const DenseMatrix<double>& jacobian,
      const DenseMatrix<double>& inverse_jacobian,
      const DenseMatrix<double>& jacobian2,
      DShape& dbasis,
      DShape& d2basis) const;
 
    template<unsigned DIM>
    void transform_second_derivatives_template(
      const DenseMatrix<double>& jacobian,
      const DenseMatrix<double>& inverse_jacobian,
      const DenseMatrix<double>& jacobian2,
      DShape& dbasis,
      DShape& d2basis) const;
 
    template<unsigned DIM>
    void transform_second_derivatives_diagonal(
      const DenseMatrix<double>& jacobian,
      const DenseMatrix<double>& inverse_jacobian,
      const DenseMatrix<double>& jacobian2,
      DShape& dbasis,
      DShape& d2basis) const;
 
    MacroElement* Macro_elem_pt;
 
    virtual void fill_in_jacobian_from_nodal_by_fd(
      Vector<double>& residuals, DenseMatrix<double>& jacobian);
 
    void fill_in_jacobian_from_nodal_by_fd(DenseMatrix<double>& jacobian)
    {
      // Allocate storage for a residuals vector and initialise to zero
      unsigned n_dof = ndof();
      Vector<double> residuals(n_dof, 0.0);
      // Get the residuals for the entire element
      get_residuals(residuals);
      // Call the jacobian calculation
      fill_in_jacobian_from_nodal_by_fd(residuals, jacobian);
    }
 
    virtual inline void update_before_nodal_fd() {}
 
    virtual inline void reset_after_nodal_fd() {}
 
    virtual inline void update_in_nodal_fd(const unsigned& i) {}
 
    virtual inline void reset_in_nodal_fd(const unsigned& i)
    {
      update_in_nodal_fd(i);
    }
 
 
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Add the contribution to the residuals
      fill_in_contribution_to_residuals(residuals);
      // Allocate storage for the full residuals (residuals of entire element)
      unsigned n_dof = ndof();
      Vector<double> full_residuals(n_dof);
      // Get the residuals for the entire element
      get_residuals(full_residuals);
      // Calculate the contributions from the internal dofs
      //(finite-difference the lot by default)
      fill_in_jacobian_from_internal_by_fd(full_residuals, jacobian, true);
      // Calculate the contributions from the external dofs
      //(finite-difference the lot by default)
      fill_in_jacobian_from_external_by_fd(full_residuals, jacobian, true);
      // Calculate the contributions from the nodal dofs
      fill_in_jacobian_from_nodal_by_fd(full_residuals, jacobian);
    }
 
  public:
    typedef void (*SteadyExactSolutionFctPt)(const Vector<double>&,
                                             Vector<double>&);
 
    typedef void (*UnsteadyExactSolutionFctPt)(const double&,
                                               const Vector<double>&,
                                               Vector<double>&);
 
    static double Tolerance_for_singular_jacobian;
 
    static bool Accept_negative_jacobian;
 
    static bool Suppress_output_while_checking_for_inverted_elements;
 
    FiniteElement()
      : GeneralisedElement(),
        Integral_pt(0),
        Node_pt(0),
        Nodal_local_eqn(0),
        Nnode(0),
        Elemental_dimension(0),
        Nodal_dimension(0),
        Nnodal_position_type(1),
        Macro_elem_pt(0)
    {
    }
 
    virtual ~FiniteElement();
 
    FiniteElement(const FiniteElement&) = delete;
 
    // Commented out broken assignment operator because this can lead to a
    // conflict warning when used in the virtual inheritence hierarchy.
    // Essentially the compiler doesn't realise that two separate
    // implementations of the broken function are the same and so, quite
    // rightly, it shouts.
    /*void operator=(const FiniteElement&) = delete;*/
 
    virtual bool local_coord_is_valid(const Vector<double>& s)
    {
      throw OomphLibError(
        "local_coord_is_valid is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      return true;
    }
 
    virtual void move_local_coord_back_into_element(Vector<double>& s) const
    {
      throw OomphLibError("move_local_coords_back_into_element() is not "
                          "implemented for this element\n",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void get_centre_of_gravity_and_max_radius_in_terms_of_zeta(
      Vector<double>& cog, double& max_radius) const;
 
    virtual void local_coordinate_of_node(const unsigned& j,
                                          Vector<double>& s) const
    {
      throw OomphLibError(
        "local_coordinate_of_node(...) is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void local_fraction_of_node(const unsigned& j,
                                        Vector<double>& s_fraction);
 
    virtual double local_one_d_fraction_of_node(const unsigned& n1d,
                                                const unsigned& i)
    {
      std::string error_message =
        "local one_d_fraction_of_node is not implemented for this element\n";
      error_message +=
        "It only makes sense for elements that use tensor-product expansions\n";
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void set_macro_elem_pt(MacroElement* macro_elem_pt)
    {
      Macro_elem_pt = macro_elem_pt;
    }
 
    MacroElement* macro_elem_pt()
    {
      return Macro_elem_pt;
    }
 
    void get_x(const Vector<double>& s, Vector<double>& x) const
    {
      // If there is no macro element then return interpolated x
      if (Macro_elem_pt == 0)
      {
        interpolated_x(s, x);
      }
      // Otherwise call the macro element representation
      else
      {
        get_x_from_macro_element(s, x);
      }
    }
 
 
    void get_x(const unsigned& t, const Vector<double>& s, Vector<double>& x)
    {
      // Get timestepper from first node
      TimeStepper* time_stepper_pt = node_pt(0)->time_stepper_pt();
 
      // Number of previous values
      const unsigned nprev = time_stepper_pt->nprev_values();
 
      // If t > nprev_values(), we're not dealing with a previous value
      // but a generalised history value -- this cannot be recovered from
      // macro element but must be determined by finite element interpolation
 
      // If there is no macro element, or we're dealing with a generalised
      // history value then use the FE representation
      if ((Macro_elem_pt == 0) || (t > nprev))
      {
        interpolated_x(t, s, x);
      }
      // Otherwise use the macro element representation
      else
      {
        get_x_from_macro_element(t, s, x);
      }
    }
 
 
    virtual void get_x_from_macro_element(const Vector<double>& s,
                                          Vector<double>& x) const
    {
      throw OomphLibError(
        "get_x_from_macro_element(...) is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void get_x_from_macro_element(const unsigned& t,
                                          const Vector<double>& s,
                                          Vector<double>& x)
    {
      throw OomphLibError(
        "get_x_from_macro_element(...) is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void set_integration_scheme(Integral* const& integral_pt);
 
    Integral* const& integral_pt() const
    {
      return Integral_pt;
    }
 
    virtual void shape(const Vector<double>& s, Shape& psi) const = 0;
 
    virtual void shape_at_knot(const unsigned& ipt, Shape& psi) const;
 
    virtual void dshape_local(const Vector<double>& s,
                              Shape& psi,
                              DShape& dpsids) const
    {
      throw OomphLibError(
        "dshape_local() is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void dshape_local_at_knot(const unsigned& ipt,
                                      Shape& psi,
                                      DShape& dpsids) const;
 
    virtual void d2shape_local(const Vector<double>& s,
                               Shape& psi,
                               DShape& dpsids,
                               DShape& d2psids) const
    {
      throw OomphLibError(
        "d2shape_local() is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void d2shape_local_at_knot(const unsigned& ipt,
                                       Shape& psi,
                                       DShape& dpsids,
                                       DShape& d2psids) const;
 
    virtual double J_eulerian(const Vector<double>& s) const;
 
    virtual double J_eulerian_at_knot(const unsigned& ipt) const;
 
    void check_J_eulerian_at_knots(bool& passed) const;
 
    void check_jacobian(const double& jacobian) const;
 
    double dshape_eulerian(const Vector<double>& s,
                           Shape& psi,
                           DShape& dpsidx) const;
 
    virtual double dshape_eulerian_at_knot(const unsigned& ipt,
                                           Shape& psi,
                                           DShape& dpsidx) const;
 
    virtual double dshape_eulerian_at_knot(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsi,
      DenseMatrix<double>& djacobian_dX,
      RankFourTensor<double>& d_dpsidx_dX) const;
 
    double d2shape_eulerian(const Vector<double>& s,
                            Shape& psi,
                            DShape& dpsidx,
                            DShape& d2psidx) const;
 
 
    virtual double d2shape_eulerian_at_knot(const unsigned& ipt,
                                            Shape& psi,
                                            DShape& dpsidx,
                                            DShape& d2psidx) const;
 
    virtual void assign_nodal_local_eqn_numbers(const bool& store_local_dof_pt);
 
    virtual void describe_local_dofs(std::ostream& out,
                                     const std::string& current_string) const;
 
    virtual void describe_nodal_local_dofs(
      std::ostream& out, const std::string& current_string) const;
 
    virtual inline void assign_all_generic_local_eqn_numbers(
      const bool& store_local_dof_pt)
    {
      // GeneralisedElement's version assigns internal and external data
      GeneralisedElement::assign_all_generic_local_eqn_numbers(
        store_local_dof_pt);
      // Need simply to add the nodal numbering scheme
      assign_nodal_local_eqn_numbers(store_local_dof_pt);
    }
 
    Node*& node_pt(const unsigned& n)
    {
#ifdef RANGE_CHECKING
      if (n >= Nnode)
      {
        std::ostringstream error_message;
        error_message << "Range Error: " << n << " is not in the range (0,"
                      << Nnode - 1 << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Node_pt[n];
    }
 
    Node* const& node_pt(const unsigned& n) const
    {
#ifdef RANGE_CHECKING
      if (n >= Nnode)
      {
        std::ostringstream error_message;
        error_message << "Range Error: " << n << " is not in the range (0,"
                      << Nnode - 1 << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      return Node_pt[n];
    }
 
    unsigned nnode() const
    {
      return Nnode;
    }
 
    virtual unsigned nnode_1d() const
    {
      return 0;
    }
 
    double raw_nodal_position(const unsigned& n, const unsigned& i) const;
    /* { */
    /*  /\* oomph_info << "n: "<< n << std::endl; *\/ */
    /*  /\* oomph_info << "i: "<< i << std::endl; *\/ */
    /*  /\* oomph_info << "node_pt(n): "<< node_pt(n) << std::endl; *\/ */
    /*  /\* oomph_info << "node_pt(n)->x(i): "<< node_pt(n)->x(i) << std::endl;
     * *\/ */
    /*  double tmp=node_pt(n)->x(i); */
    /*  //oomph_info << "tmp: "<< tmp << std::endl; */
    /*  return tmp; // node_pt(n)->x(i); */
    /* } */
 
    double raw_nodal_position(const unsigned& t,
                              const unsigned& n,
                              const unsigned& i) const
    {
      return node_pt(n)->x(t, i);
    }
 
    double raw_dnodal_position_dt(const unsigned& n, const unsigned& i) const
    {
      return node_pt(n)->dx_dt(i);
    }
 
    double raw_dnodal_position_dt(const unsigned& n,
                                  const unsigned& j,
                                  const unsigned& i) const
    {
      return node_pt(n)->dx_dt(j, i);
    }
 
    double raw_nodal_position_gen(const unsigned& n,
                                  const unsigned& k,
                                  const unsigned& i) const
    {
      return node_pt(n)->x_gen(k, i);
    }
 
    double raw_nodal_position_gen(const unsigned& t,
                                  const unsigned& n,
                                  const unsigned& k,
                                  const unsigned& i) const
    {
      return node_pt(n)->x_gen(t, k, i);
    }
 
    double raw_dnodal_position_gen_dt(const unsigned& n,
                                      const unsigned& k,
                                      const unsigned& i) const
    {
      return node_pt(n)->dx_gen_dt(k, i);
    }
 
    double raw_dnodal_position_gen_dt(const unsigned& j,
                                      const unsigned& n,
                                      const unsigned& k,
                                      const unsigned& i) const
    {
      return node_pt(n)->dx_gen_dt(j, k, i);
    }
 
 
    double nodal_position(const unsigned& n, const unsigned& i) const
    {
      return node_pt(n)->position(i);
    }
 
    double nodal_position(const unsigned& t,
                          const unsigned& n,
                          const unsigned& i) const
    {
      return node_pt(n)->position(t, i);
    }
 
    double dnodal_position_dt(const unsigned& n, const unsigned& i) const
    {
      return node_pt(n)->dposition_dt(i);
    }
 
    double dnodal_position_dt(const unsigned& n,
                              const unsigned& j,
                              const unsigned& i) const
    {
      return node_pt(n)->dposition_dt(j, i);
    }
 
    double nodal_position_gen(const unsigned& n,
                              const unsigned& k,
                              const unsigned& i) const
    {
      return node_pt(n)->position_gen(k, i);
    }
 
    double nodal_position_gen(const unsigned& t,
                              const unsigned& n,
                              const unsigned& k,
                              const unsigned& i) const
    {
      return node_pt(n)->position_gen(t, k, i);
    }
 
    double dnodal_position_gen_dt(const unsigned& n,
                                  const unsigned& k,
                                  const unsigned& i) const
    {
      return node_pt(n)->dposition_gen_dt(k, i);
    }
 
    double dnodal_position_gen_dt(const unsigned& j,
                                  const unsigned& n,
                                  const unsigned& k,
                                  const unsigned& i) const
    {
      return node_pt(n)->dposition_gen_dt(j, k, i);
    }
 
 
    virtual void get_dresidual_dnodal_coordinates(
      RankThreeTensor<double>& dresidual_dnodal_coordinates);
 
    virtual void disable_ALE()
    {
      std::ostringstream warn_message;
      warn_message
        << "Warning: You have just called the default (empty) function \n\n"
        << "    FiniteElement::disable_ALE() \n\n"
        << "This suggests that you've either tried to call it for an element\n"
        << "that \n"
        << "(1) does not involve time-derivatives (e.g. a Poisson element) \n"
        << "(2) an element for which the time-derivatives aren't implemented \n"
        << "    in ALE form \n"
        << "(3) an element for which the ALE form of the equations can't be \n"
        << "    be disabled (yet).\n";
      OomphLibWarning(
        warn_message.str(), "Problem::disable_ALE()", OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void enable_ALE()
    {
      std::ostringstream warn_message;
      warn_message
        << "Warning: You have just called the default (empty) function \n\n"
        << "    FiniteElement::enable_ALE() \n\n"
        << "This suggests that you've either tried to call it for an element\n"
        << "that \n"
        << "(1) does not involve time-derivatives (e.g. a Poisson element) \n"
        << "(2) an element for which the time-derivatives aren't implemented \n"
        << "    in ALE form \n"
        << "(3) an element for which the ALE form of the equations can't be \n"
        << "    be disabled (yet)\n"
        << "(4) an element for which this function has not been (properly) \n "
        << "    implemented. This is likely to be a bug!\n ";
      OomphLibWarning(
        warn_message.str(), "Problem::enable_ALE()", OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual unsigned required_nvalue(const unsigned& n) const
    {
      return Default_Initial_Nvalue;
    }
 
    unsigned nnodal_position_type() const
    {
      return Nnodal_position_type;
    }
 
    bool has_hanging_nodes() const
    {
      unsigned nnod = nnode();
      for (unsigned j = 0; j < nnod; j++)
      {
        if (node_pt(j)->is_hanging())
        {
          return true;
        }
      }
      return false;
    }
 
    unsigned nodal_dimension() const
    {
      return Nodal_dimension;
    }
 
    virtual unsigned nvertex_node() const
    {
      std::string error_msg = "Not implemented for FiniteElement.";
      throw OomphLibError(
        error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual Node* vertex_node_pt(const unsigned& j) const
    {
      std::string error_msg = "Not implemented for FiniteElement.";
      throw OomphLibError(
        error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual Node* construct_node(const unsigned& n)
    {
      // Create a node and assign it to the local node pointer
      // The dimension and number of values are taken from internal element data
      node_pt(n) =
        new Node(Nodal_dimension, Nnodal_position_type, required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can use it
      return node_pt(n);
    }
 
    virtual Node* construct_node(const unsigned& n,
                                 TimeStepper* const& time_stepper_pt)
    {
      // Create a node and assign it to the local node pointer.
      // The dimension and number of values are taken from internal element data
      node_pt(n) = new Node(time_stepper_pt,
                            Nodal_dimension,
                            Nnodal_position_type,
                            required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return node_pt(n);
    }
 
    virtual Node* construct_boundary_node(const unsigned& n)
    {
      // Create a node and assign it to the local node pointer
      // The dimension and number of values are taken from internal element data
      node_pt(n) = new BoundaryNode<Node>(
        Nodal_dimension, Nnodal_position_type, required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can use it
      return node_pt(n);
    }
 
    virtual Node* construct_boundary_node(const unsigned& n,
                                          TimeStepper* const& time_stepper_pt)
    {
      // Create a node and assign it to the local node pointer.
      // The dimension and number of values are taken from internal element data
      node_pt(n) = new BoundaryNode<Node>(time_stepper_pt,
                                          Nodal_dimension,
                                          Nnodal_position_type,
                                          required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return node_pt(n);
    }
 
    int get_node_number(Node* const& node_pt) const;
 
 
    virtual Node* get_node_at_local_coordinate(const Vector<double>& s) const;
 
    double raw_nodal_value(const unsigned& n, const unsigned& i) const
    {
      return node_pt(n)->raw_value(i);
    }
 
    double raw_nodal_value(const unsigned& t,
                           const unsigned& n,
                           const unsigned& i) const
    {
      return node_pt(n)->raw_value(t, i);
    }
 
    double nodal_value(const unsigned& n, const unsigned& i) const
    {
      return node_pt(n)->value(i);
    }
 
    double nodal_value(const unsigned& t,
                       const unsigned& n,
                       const unsigned& i) const
    {
      return node_pt(n)->value(t, i);
    }
 
    unsigned dim() const
    {
      return Elemental_dimension;
    }
 
    virtual ElementGeometry::ElementGeometry element_geometry() const
    {
      std::string err = "Broken virtual function.";
      throw OomphLibError(
        err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual double interpolated_x(const Vector<double>& s,
                                  const unsigned& i) const;
 
    virtual double interpolated_x(const unsigned& t,
                                  const Vector<double>& s,
                                  const unsigned& i) const;
 
    virtual void interpolated_x(const Vector<double>& s,
                                Vector<double>& x) const;
 
    virtual void interpolated_x(const unsigned& t,
                                const Vector<double>& s,
                                Vector<double>& x) const;
 
    virtual double interpolated_dxdt(const Vector<double>& s,
                                     const unsigned& i,
                                     const unsigned& t);
 
    virtual void interpolated_dxdt(const Vector<double>& s,
                                   const unsigned& t,
                                   Vector<double>& dxdt);
 
    inline unsigned ngeom_data() const
    {
      return 0;
    }
 
    inline Data* geom_data_pt(const unsigned& j)
    {
      return 0;
    }
 
    void position(const Vector<double>& zeta, Vector<double>& r) const
    {
      this->interpolated_x(zeta, r);
    }
 
    void position(const unsigned& t,
                  const Vector<double>& zeta,
                  Vector<double>& r) const
    {
      this->interpolated_x(t, zeta, r);
    }
 
    void dposition_dt(const Vector<double>& zeta,
                      const unsigned& t,
                      Vector<double>& drdt)
    {
      this->interpolated_dxdt(zeta, t, drdt);
    }
 
    virtual double zeta_nodal(const unsigned& n,
                              const unsigned& k,
                              const unsigned& i) const
    {
      // By default return the value for nodal_position_gen
      return nodal_position_gen(n, k, i);
    }
 
 
    void interpolated_zeta(const Vector<double>& s, Vector<double>& zeta) const;
 
    void locate_zeta(const Vector<double>& zeta,
                     GeomObject*& geom_object_pt,
                     Vector<double>& s,
                     const bool& use_coordinate_as_initial_guess = false);
 
 
    virtual void node_update();
 
    virtual void identify_field_data_for_interactions(
      std::set<std::pair<Data*, unsigned>>& paired_field_data);
 
 
    virtual void identify_geometric_data(std::set<Data*>& geometric_data_pt) {}
 
 
    virtual double s_min() const
    {
      throw OomphLibError("s_min() isn't implemented for this element\n",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
      // Dummy return
      return 0.0;
    }
 
    virtual double s_max() const
    {
      throw OomphLibError("s_max() isn't implemented for this element\n",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
      // Dummy return
      return 0.0;
    }
 
    double size() const;
 
 
    virtual double compute_physical_size() const
    {
      throw OomphLibError(
        "compute_physical_size() isn't implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      // Dummy return
      return 0.0;
    }
 
    virtual void point_output_data(const Vector<double>& s,
                                   Vector<double>& data)
    {
    }
 
    void point_output(std::ostream& outfile, const Vector<double>& s)
    {
      // Get point data
      Vector<double> data;
      this->point_output_data(s, data);
 
      // Output
      unsigned n = data.size();
      for (unsigned i = 0; i < n; i++)
      {
        outfile << data[i] << " ";
      }
    }
 
    virtual unsigned nplot_points_paraview(const unsigned& nplot) const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
 
      // Dummy unsigned
      return 0;
    }
 
    virtual unsigned nsub_elements_paraview(const unsigned& nplot) const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
 
      // Dummy unsigned
      return 0;
    }
 
    void output_paraview(std::ofstream& file_out, const unsigned& nplot) const
    {
      // Decide the dimensions of the nodes
      unsigned nnod = nnode();
      if (nnod == 0) return;
      unsigned n = node_pt(0)->ndim();
 
      // Vector for local coordinates
      Vector<double> s(n, 0.0);
 
      // Vector for cartesian coordinates
      Vector<double> x(n, 0.0);
 
      // Determine the total number of plotpoints
      unsigned plot = nplot_points_paraview(nplot);
 
      // Loop over the local points
      for (unsigned j = 0; j < plot; j++)
      {
        // Determine where in the local element the point is
        this->get_s_plot(j, nplot, s);
 
        // Update the cartesian coordinates vector
        this->interpolated_x(s, x);
 
        // Print the global coordinates. Note: no whitespace after last
        // coordinate or paraview is very unhappy.
        for (unsigned i = 0; i < n - 1; i++)
        {
          file_out << x[i] << " ";
        }
        file_out << x[n - 1];
 
        // Since unstructured grid always needs 3 components for each
        // point, output 0's by default
        switch (n)
        {
          case 1:
            file_out << " 0"
                     << " 0" << std::endl;
            break;
 
          case 2:
            file_out << " 0" << std::endl;
            break;
 
          case 3:
            file_out << std::endl;
            break;
 
          // Paraview can't handle more than 3 dimensions, output error
          default:
            throw OomphLibError(
              "Printing PlotPoint to .vtu failed; it has >3 dimensions.",
              OOMPH_CURRENT_FUNCTION,
              OOMPH_EXCEPTION_LOCATION);
        }
      }
    }
 
    virtual void write_paraview_output_offset_information(
      std::ofstream& file_out, const unsigned& nplot, unsigned& counter) const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void write_paraview_type(std::ofstream& file_out,
                                     const unsigned& nplot) const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void write_paraview_offsets(std::ofstream& file_out,
                                        const unsigned& nplot,
                                        unsigned& offset_sum) const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual unsigned nscalar_paraview() const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
 
      // Dummy unsigned
      return 0;
    }
 
    virtual void scalar_value_paraview(std::ofstream& file_out,
                                       const unsigned& i,
                                       const unsigned& nplot) const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void scalar_value_fct_paraview(
      std::ofstream& file_out,
      const unsigned& i,
      const unsigned& nplot,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void scalar_value_fct_paraview(
      std::ofstream& file_out,
      const unsigned& i,
      const unsigned& nplot,
      const double& time,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
    {
      throw OomphLibError(
        "This function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual std::string scalar_name_paraview(const unsigned& i) const
    {
      return "V" + StringConversion::to_string(i);
    }
 
    virtual void output(std::ostream& outfile)
    {
      throw OomphLibError(
        "Output function function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void output(std::ostream& outfile, const unsigned& n_plot)
    {
      throw OomphLibError(
        "Output function function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void output(const unsigned& t,
                        std::ostream& outfile,
                        const unsigned& n_plot) const
    {
      throw OomphLibError(
        "Output function function hasn't been implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void output(FILE* file_pt)
    {
      throw OomphLibError("C-style otput function function hasn't been "
                          "implemented for this element",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void output(FILE* file_pt, const unsigned& n_plot)
    {
      throw OomphLibError("C-style output function function hasn't been "
                          "implemented for this element",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void output_fct(
      std::ostream& outfile,
      const unsigned& n_plot,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
    {
      throw OomphLibError(
        "Output function function hasn't been implemented for exact solution",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void output_fct(
      std::ostream& outfile,
      const unsigned& n_plot,
      const double& time,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
    {
      throw OomphLibError(
        "Output function function hasn't been implemented for exact solution",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void output_fct(std::ostream& outfile,
                            const unsigned& n_plot,
                            const double& time,
                            const SolutionFunctorBase& exact_soln) const
    {
      throw OomphLibError(
        "Output function function hasn't been implemented for exact solution",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void get_s_plot(const unsigned& i,
                            const unsigned& nplot,
                            Vector<double>& s,
                            const bool& shifted_to_interior = false) const
    {
      throw OomphLibError(
        "get_s_plot(...) is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    };
 
    virtual std::string tecplot_zone_string(const unsigned& nplot) const
    {
      throw OomphLibError(
        "tecplot_zone_string(...) is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      return "dummy return";
    }
 
    virtual void write_tecplot_zone_footer(std::ostream& outfile,
                                           const unsigned& nplot) const {};
 
    virtual void write_tecplot_zone_footer(FILE* file_pt,
                                           const unsigned& nplot) const {};
 
    virtual unsigned nplot_points(const unsigned& nplot) const
    {
      throw OomphLibError(
        "nplot_points(...) is not implemented for this element",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      // Dummy return
      return 0;
    }
 
 
    virtual void compute_error(
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      double& error,
      double& norm)
    {
      std::string error_message = "compute_error undefined for this element \n";
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void compute_error(
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
      const double& time,
      double& error,
      double& norm)
    {
      std::string error_message = "time-dependent compute_error ";
      error_message += "undefined for this element \n";
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void compute_error(
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      Vector<double>& error,
      Vector<double>& norm)
    {
      std::string error_message = "compute_error undefined for this element \n";
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void compute_error(
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
      const double& time,
      Vector<double>& error,
      Vector<double>& norm)
    {
      std::string error_message = "time-dependent compute_error ";
      error_message += "undefined for this element \n";
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void compute_error(
      std::ostream& outfile,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      double& error,
      double& norm)
    {
      std::string error_message = "compute_error undefined for this element \n";
      outfile << error_message;
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void compute_error(
      std::ostream& outfile,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
      const double& time,
      double& error,
      double& norm)
    {
      std::string error_message =
        "time-dependent compute_error undefined for this element \n";
      outfile << error_message;
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void compute_error(
      std::ostream& outfile,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      Vector<double>& error,
      Vector<double>& norm)
    {
      std::string error_message = "compute_error undefined for this element \n";
      outfile << error_message;
 
      throw OomphLibError(error_message,
                          "FiniteElement::compute_error()",
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void compute_error(
      std::ostream& outfile,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
      const double& time,
      Vector<double>& error,
      Vector<double>& norm)
    {
      std::string error_message =
        "time-dependent compute_error undefined for this element \n";
      outfile << error_message;
 
      throw OomphLibError(error_message,
                          "FiniteElement::compute_error()",
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void compute_abs_error(
      std::ostream& outfile,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      double& error)
    {
      std::string error_message =
        "compute_abs_error undefined for this element \n";
      outfile << error_message;
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    void integrate_fct(FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt,
                       Vector<double>& integral);
 
 
    void integrate_fct(
      FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt,
      const double& time,
      Vector<double>& integral);
 
    virtual void build_face_element(const int& face_index,
                                    FaceElement* face_element_pt);
 
 
    virtual unsigned self_test();
 
    virtual unsigned get_bulk_node_number(const int& face_index,
                                          const unsigned& i) const
    {
      std::string err = "Not implemented for this element.";
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
 
    virtual int face_outer_unit_normal_sign(const int& face_index) const
    {
      std::string err = "Not implemented for this element.";
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
 
    virtual unsigned nnode_on_face() const
    {
      std::string err = "Not implemented for this element.";
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
 
    void face_node_number_error_check(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      if (i > nnode_on_face())
      {
        std::string err = "Face node index i out of range on face.";
        throw OomphLibError(
          err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
      }
#endif
    }
 
    virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt(
      const int& face_index) const
    {
      std::string err = "Not implemented for this element.";
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
 
    virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt(
      const int& face_index) const
    {
      std::string err = "Not implemented for this element.";
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  class PointElement : public virtual FiniteElement
  {
  public:
    PointElement()
    {
      this->set_n_node(1);
 
      // Set the integration scheme
      this->set_integration_scheme(&Default_integration_scheme);
    }
 
    PointElement(const PointElement&) = delete;
 
    /*void operator=(const PointElement&) = delete;*/
 
    void shape(const Vector<double>& s, Shape& psi) const;
 
    void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
    {
      s.resize(0);
    }
 
    // unsigned dim() const {return 0;}
 
  private:
    static PointIntegral Default_integration_scheme;
  };
 
 
 
 
  class GeomObject;
 
  //=======================================================================
  //=======================================================================
  class SolidInitialCondition
  {
  public:
    SolidInitialCondition(GeomObject* geom_object_pt)
      : Geom_object_pt(geom_object_pt), IC_time_deriv(0){};
 
 
    SolidInitialCondition(const SolidInitialCondition&) = delete;
 
    void operator=(const SolidInitialCondition&) = delete;
 
    GeomObject*& geom_object_pt()
    {
      return Geom_object_pt;
    }
 
    unsigned& ic_time_deriv()
    {
      return IC_time_deriv;
    }
 
 
  private:
    GeomObject* Geom_object_pt;
 
    unsigned IC_time_deriv;
  };
 
 
 
 
  //============================================================================
  //============================================================================
  class SolidFiniteElement : public virtual FiniteElement
  {
  public:
    void set_lagrangian_dimension(const unsigned& lagrangian_dimension)
    {
      Lagrangian_dimension = lagrangian_dimension;
    }
 
    virtual bool has_internal_solid_data()
    {
      return false;
    }
 
    typedef double (*MultiplierFctPt)(const Vector<double>& xi);
 
    SolidFiniteElement()
      : FiniteElement(),
        Undeformed_macro_elem_pt(0),
        Solid_ic_pt(0),
        Multiplier_fct_pt(0),
        Position_local_eqn(0),
        Lagrangian_dimension(0),
        Nnodal_lagrangian_type(1),
        Solve_for_consistent_newmark_accel_flag(false)
    {
    }
 
    virtual ~SolidFiniteElement();
 
    SolidFiniteElement(const SolidFiniteElement&) = delete;
 
    /*void operator=(const SolidFiniteElement&) = delete;*/
 
    inline unsigned ngeom_data() const
    {
      return nnode();
    }
 
    inline Data* geom_data_pt(const unsigned& j)
    {
      return static_cast<SolidNode*>(node_pt(j))->variable_position_pt();
    }
 
    void identify_geometric_data(std::set<Data*>& geometric_data_pt)
    {
      // Loop over the node update data and add to the set
      const unsigned n_node = this->nnode();
      for (unsigned n = 0; n < n_node; n++)
      {
        geometric_data_pt.insert(
          dynamic_cast<SolidNode*>(this->node_pt(n))->variable_position_pt());
      }
    }
 
 
    inline double zeta_nodal(const unsigned& n,
                             const unsigned& k,
                             const unsigned& i) const
    {
      return lagrangian_position_gen(n, k, i);
    }
 
    virtual void get_x_and_xi(const Vector<double>& s,
                              Vector<double>& x_fe,
                              Vector<double>& x,
                              Vector<double>& xi_fe,
                              Vector<double>& xi) const
    {
      throw OomphLibError(
        "get_x_and_xi(...) is not implemented for this element\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void set_macro_elem_pt(MacroElement* macro_elem_pt)
    {
      Macro_elem_pt = macro_elem_pt;
      Undeformed_macro_elem_pt = macro_elem_pt;
    }
 
    virtual void set_macro_elem_pt(MacroElement* macro_elem_pt,
                                   MacroElement* undeformed_macro_elem_pt)
    {
      Macro_elem_pt = macro_elem_pt;
      Undeformed_macro_elem_pt = undeformed_macro_elem_pt;
    }
 
    void set_undeformed_macro_elem_pt(MacroElement* undeformed_macro_elem_pt)
    {
      Undeformed_macro_elem_pt = undeformed_macro_elem_pt;
    }
 
 
    MacroElement* undeformed_macro_elem_pt()
    {
      return Undeformed_macro_elem_pt;
    }
 
    double dshape_lagrangian(const Vector<double>& s,
                             Shape& psi,
                             DShape& dpsidxi) const;
 
    virtual double dshape_lagrangian_at_knot(const unsigned& ipt,
                                             Shape& psi,
                                             DShape& dpsidxi) const;
 
    double d2shape_lagrangian(const Vector<double>& s,
                              Shape& psi,
                              DShape& dpsidxi,
                              DShape& d2psidxi) const;
 
    virtual double d2shape_lagrangian_at_knot(const unsigned& ipt,
                                              Shape& psi,
                                              DShape& dpsidxi,
                                              DShape& d2psidxi) const;
 
    unsigned lagrangian_dimension() const
    {
      return Lagrangian_dimension;
    }
 
    unsigned nnodal_lagrangian_type() const
    {
      return Nnodal_lagrangian_type;
    }
 
    Node* construct_node(const unsigned& n)
    {
      // Construct a solid node and assign it to the local node pointer vector.
      // The dimension and number of values are taken from internal element data
      // The number of solid pressure dofs are also taken from internal data
      // The number of timesteps to be stored comes from the problem!
      node_pt(n) = new SolidNode(lagrangian_dimension(),
                                 nnodal_lagrangian_type(),
                                 nodal_dimension(),
                                 nnodal_position_type(),
                                 required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return node_pt(n);
    }
 
    Node* construct_node(const unsigned& n, TimeStepper* const& time_stepper_pt)
    {
      // Construct a solid node and assign it to the local node pointer vector
      // The dimension and number of values are taken from internal element data
      // The number of solid pressure dofs are also taken from internal data
      node_pt(n) = new SolidNode(time_stepper_pt,
                                 lagrangian_dimension(),
                                 nnodal_lagrangian_type(),
                                 nodal_dimension(),
                                 nnodal_position_type(),
                                 required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return node_pt(n);
    }
 
    Node* construct_boundary_node(const unsigned& n)
    {
      // Construct a solid node and assign it to the local node pointer vector.
      // The dimension and number of values are taken from internal element data
      // The number of solid pressure dofs are also taken from internal data
      // The number of timesteps to be stored comes from the problem!
      node_pt(n) = new BoundaryNode<SolidNode>(lagrangian_dimension(),
                                               nnodal_lagrangian_type(),
                                               nodal_dimension(),
                                               nnodal_position_type(),
                                               required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return node_pt(n);
    }
 
    Node* construct_boundary_node(const unsigned& n,
                                  TimeStepper* const& time_stepper_pt)
    {
      // Construct a solid node and assign it to the local node pointer vector
      // The dimension and number of values are taken from internal element data
      // The number of solid pressure dofs are also taken from internal data
      node_pt(n) = new BoundaryNode<SolidNode>(time_stepper_pt,
                                               lagrangian_dimension(),
                                               nnodal_lagrangian_type(),
                                               nodal_dimension(),
                                               nnodal_position_type(),
                                               required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return node_pt(n);
    }
 
    virtual inline void assign_all_generic_local_eqn_numbers(
      const bool& store_local_dof_pt)
    {
      // Call the standard finite element equation numbering
      //(internal, external and nodal data).
      FiniteElement::assign_all_generic_local_eqn_numbers(store_local_dof_pt);
      // Assign the numbering for the solid dofs
      assign_solid_local_eqn_numbers(store_local_dof_pt);
    }
 
    void describe_local_dofs(std::ostream& out,
                             const std::string& current_string) const;
 
    double raw_lagrangian_position(const unsigned& n, const unsigned& i) const
    {
      return static_cast<SolidNode*>(node_pt(n))->xi(i);
    }
 
    double raw_lagrangian_position_gen(const unsigned& n,
                                       const unsigned& k,
                                       const unsigned& i) const
    {
      return static_cast<SolidNode*>(node_pt(n))->xi_gen(k, i);
    }
 
    double lagrangian_position(const unsigned& n, const unsigned& i) const
    {
      return static_cast<SolidNode*>(node_pt(n))->lagrangian_position(i);
    }
 
    double lagrangian_position_gen(const unsigned& n,
                                   const unsigned& k,
                                   const unsigned& i) const
    {
      return static_cast<SolidNode*>(node_pt(n))->lagrangian_position_gen(k, i);
    }
 
    virtual double interpolated_xi(const Vector<double>& s,
                                   const unsigned& i) const;
 
    virtual void interpolated_xi(const Vector<double>& s,
                                 Vector<double>& xi) const;
 
    virtual void interpolated_dxids(const Vector<double>& s,
                                    DenseMatrix<double>& dxids) const;
 
    virtual void J_lagrangian(const Vector<double>& s) const
    {
      // Must be implemented and overloaded in FaceElements
      throw OomphLibError("Function not implemented yet",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual double J_lagrangian_at_knot(const unsigned& ipt) const
    {
      // Must be implemented and overloaded in FaceElements
      throw OomphLibError("Function not implemented yet",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    SolidInitialCondition*& solid_ic_pt()
    {
      return Solid_ic_pt;
    }
 
    void enable_solve_for_consistent_newmark_accel()
    {
      Solve_for_consistent_newmark_accel_flag = true;
    }
 
    void disable_solve_for_consistent_newmark_accel()
    {
      Solve_for_consistent_newmark_accel_flag = false;
    }
 
    MultiplierFctPt& multiplier_fct_pt()
    {
      return Multiplier_fct_pt;
    }
 
 
    MultiplierFctPt multiplier_fct_pt() const
    {
      return Multiplier_fct_pt;
    }
 
 
    virtual void get_residuals_for_solid_ic(Vector<double>& residuals)
    {
      residuals.initialise(0.0);
      fill_in_residuals_for_solid_ic(residuals);
    }
 
    void fill_in_residuals_for_solid_ic(Vector<double>& residuals)
    {
      // Call the generic residuals function with flag set to 0
      // using a dummy matrix argument
      fill_in_generic_jacobian_for_solid_ic(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    void fill_in_jacobian_for_solid_ic(Vector<double>& residuals,
                                       DenseMatrix<double>& jacobian)
    {
      // Call the generic routine with the flag set to 1
      fill_in_generic_jacobian_for_solid_ic(residuals, jacobian, 1);
    }
 
 
    void fill_in_jacobian_for_newmark_accel(DenseMatrix<double>& jacobian);
 
 
    void compute_norm(double& el_norm);
 
  protected:
    void fill_in_generic_jacobian_for_solid_ic(Vector<double>& residuals,
                                               DenseMatrix<double>& jacobian,
                                               const unsigned& flag);
 
    void set_nnodal_lagrangian_type(const unsigned& nlagrangian_type)
    {
      Nnodal_lagrangian_type = nlagrangian_type;
    }
 
    MacroElement* Undeformed_macro_elem_pt;
 
    virtual double local_to_lagrangian_mapping(
      const DShape& dpsids,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& inverse_jacobian) const
    {
      // Assemble the jacobian
      assemble_local_to_lagrangian_jacobian(dpsids, jacobian);
      // Invert the jacobian
      return invert_jacobian_mapping(jacobian, inverse_jacobian);
    }
 
    double local_to_lagrangian_mapping(
      const DShape& dpsids, DenseMatrix<double>& inverse_jacobian) const
    {
      // Find the dimension of the element
      unsigned el_dim = dim();
      // Assign memory for the jacobian
      DenseMatrix<double> jacobian(el_dim);
      // Calculate the jacobian and inverse
      return local_to_lagrangian_mapping(dpsids, jacobian, inverse_jacobian);
    }
 
    virtual double local_to_lagrangian_mapping_diagonal(
      const DShape& dpsids,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& inverse_jacobian) const;
 
 
    virtual void assign_solid_local_eqn_numbers(const bool& store_local_dof);
 
    void describe_solid_local_dofs(std::ostream& out,
                                   const std::string& current_string) const;
 
    SolidInitialCondition* Solid_ic_pt;
 
  public:
    inline int position_local_eqn(const unsigned& n,
                                  const unsigned& k,
                                  const unsigned& j) const
    {
#ifdef RANGE_CHECKING
      std::ostringstream error_message;
      bool error = false;
      if (n >= nnode())
      {
        error = true;
        error_message << "Range Error: Nodal number " << n
                      << " is not in the range (0," << nnode() << ")";
      }
 
      if (k >= nnodal_position_type())
      {
        error = true;
        error_message << "Range Error: Position type " << k
                      << " is not in the range (0," << nnodal_position_type()
                      << ")";
      }
 
      if (j >= nodal_dimension())
      {
        error = true;
        error_message << "Range Error: Nodal coordinate " << j
                      << " is not in the range (0," << nodal_dimension() << ")";
      }
 
      if (error)
      {
        // Throw the error
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Return the value
      return Position_local_eqn[(n * nnodal_position_type() + k) *
                                  nodal_dimension() +
                                j];
    }
 
  protected:
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Add the contribution to the residuals
      fill_in_contribution_to_residuals(residuals);
 
      // Solve for the consistent acceleration in Newmark scheme?
      if (Solve_for_consistent_newmark_accel_flag)
      {
        fill_in_jacobian_for_newmark_accel(jacobian);
        return;
      }
 
      // Allocate storage for the full residuals (residuals of entire element)
      unsigned n_dof = ndof();
      Vector<double> full_residuals(n_dof);
      // Get the residuals for the entire element
      get_residuals(full_residuals);
      // Get the solid entries in the jacobian using finite differences
      fill_in_jacobian_from_solid_position_by_fd(full_residuals, jacobian);
      // There could be internal data
      //(finite-difference the lot by default)
      fill_in_jacobian_from_internal_by_fd(full_residuals, jacobian, true);
      // There could also be external data
      //(finite-difference the lot by default)
      fill_in_jacobian_from_external_by_fd(full_residuals, jacobian, true);
      // There could also be nodal data
      fill_in_jacobian_from_nodal_by_fd(full_residuals, jacobian);
    }
 
 
    virtual void fill_in_jacobian_from_solid_position_by_fd(
      Vector<double>& residuals, DenseMatrix<double>& jacobian);
 
 
    void fill_in_jacobian_from_solid_position_by_fd(
      DenseMatrix<double>& jacobian)
    {
      // Allocate storage for a residuals vector and initialise to zero
      unsigned n_dof = ndof();
      Vector<double> residuals(n_dof, 0.0);
      // Get the residuals for the entire element
      get_residuals(residuals);
      // Call the jacobian calculation
      fill_in_jacobian_from_solid_position_by_fd(residuals, jacobian);
    }
 
    virtual inline void update_before_solid_position_fd() {}
 
    virtual inline void reset_after_solid_position_fd() {}
 
    virtual inline void update_in_solid_position_fd(const unsigned& i) {}
 
    virtual inline void reset_in_solid_position_fd(const unsigned& i)
    {
      update_in_solid_position_fd(i);
    }
 
 
  private:
    virtual void assemble_local_to_lagrangian_jacobian(
      const DShape& dpsids, DenseMatrix<double>& jacobian) const;
 
 
    virtual void assemble_local_to_lagrangian_jacobian2(
      const DShape& d2psids, DenseMatrix<double>& jacobian2) const;
 
    MultiplierFctPt Multiplier_fct_pt;
 
    // ALH: (This is here so that the numbering can be done automatically)
    int* Position_local_eqn;
 
    unsigned Lagrangian_dimension;
 
    unsigned Nnodal_lagrangian_type;
 
  protected:
    bool Solve_for_consistent_newmark_accel_flag;
 
  private:
    inline double multiplier(const Vector<double>& xi)
    {
      // If no function has been set, return 1
      if (Multiplier_fct_pt == 0)
      {
        return 1.0;
      }
      else
      {
        // Evaluate function pointer
        return (*Multiplier_fct_pt)(xi);
      }
    }
  };
 
 
  //========================================================================
  //========================================================================
  class FaceElement : public virtual FiniteElement
  {
    typedef void (*CoordinateMappingFctPt)(const Vector<double>& s,
                                           Vector<double>& s_bulk);
 
    typedef void (*BulkCoordinateDerivativesFctPt)(
      const Vector<double>& s,
      DenseMatrix<double>& ds_bulk_dsface,
      unsigned& interior_direction);
 
  private:
    CoordinateMappingFctPt Face_to_bulk_coordinate_fct_pt;
 
    BulkCoordinateDerivativesFctPt Bulk_coordinate_derivatives_fct_pt;
 
    Vector<unsigned> Bulk_position_type;
 
    int Normal_sign;
 
    int Face_index;
 
    static bool Ignore_discontinuous_tangent_warning;
 
  protected:
    unsigned Boundary_number_in_bulk_mesh;
 
#ifdef PARANOID
 
    bool Boundary_number_in_bulk_mesh_has_been_set;
 
#endif
 
    FiniteElement* Bulk_element_pt;
 
    Vector<unsigned> Bulk_node_number;
 
    // NOTE: This breaks if the nodes have already been resized
    // before calling the face element, i.e. two separate sets of equations on
    // the face!
    Vector<unsigned> Nbulk_value;
 
    Vector<double>* Tangent_direction_pt;
 
    void add_additional_values(const Vector<unsigned>& nadditional_values,
                               const unsigned& id)
    {
      // How many nodes?
      const unsigned n_node = nnode();
 
      // loop over the nodes
      for (unsigned n = 0; n < n_node; n++)
      {
        // Assign the required number of additional nodes to the node
        dynamic_cast<BoundaryNodeBase*>(this->node_pt(n))
          ->assign_additional_values_with_face_id(nadditional_values[n], id);
      }
    }
 
 
  public:
    FaceElement()
      : Face_to_bulk_coordinate_fct_pt(0),
        Bulk_coordinate_derivatives_fct_pt(0),
        Normal_sign(0),
        Face_index(0),
        Boundary_number_in_bulk_mesh(0),
        Bulk_element_pt(0),
        Tangent_direction_pt(0)
    {
      // Check whether things have been set
#ifdef PARANOID
      Boundary_number_in_bulk_mesh_has_been_set = false;
#endif
 
      // Bulk_position_type[0] is always 0 (the position)
      Bulk_position_type.push_back(0);
    }
 
    virtual ~FaceElement() {}
 
 
    FaceElement(const FaceElement&) = delete;
 
    /*void operator=(const FaceElement&) = delete;*/
 
    inline const unsigned& boundary_number_in_bulk_mesh() const
    {
      return Boundary_number_in_bulk_mesh;
    }
 
 
    inline void set_boundary_number_in_bulk_mesh(const unsigned& b)
    {
      Boundary_number_in_bulk_mesh = b;
#ifdef PARANOID
      Boundary_number_in_bulk_mesh_has_been_set = true;
#endif
    }
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      // Vector in which to hold the intrinsic coordinate
      Vector<double> zeta(this->dim());
 
      // Get the k-th generalised boundary coordinate at node n
      this->node_pt(n)->get_coordinates_on_boundary(
        Boundary_number_in_bulk_mesh, k, zeta);
 
      // Return the individual coordinate
      return zeta[i];
    }
 
    double J_eulerian(const Vector<double>& s) const;
 
    double J_eulerian_at_knot(const unsigned& ipt) const;
 
    void check_J_eulerian_at_knots(bool& passed) const;
 
    double interpolated_x(const Vector<double>& s, const unsigned& i) const
    {
      // Local coordinates in bulk element
      Vector<double> s_bulk(dim() + 1);
      s_bulk = local_coordinate_in_bulk(s);
 
      // Return Eulerian coordinate as computed by bulk
      return bulk_element_pt()->interpolated_x(s_bulk, i);
    }
 
    double interpolated_x(const unsigned& t,
                          const Vector<double>& s,
                          const unsigned& i) const
    {
      // Local coordinates in bulk element
      Vector<double> s_bulk(dim() + 1);
      s_bulk = local_coordinate_in_bulk(s);
 
      // Return Eulerian coordinate as computed by bulk
      return bulk_element_pt()->interpolated_x(t, s_bulk, i);
    }
 
    void interpolated_x(const Vector<double>& s, Vector<double>& x) const
    {
      // Local coordinates in bulk element
      Vector<double> s_bulk(dim() + 1);
      s_bulk = local_coordinate_in_bulk(s);
 
      // Get Eulerian position vector
      bulk_element_pt()->interpolated_x(s_bulk, x);
    }
 
    void interpolated_x(const unsigned& t,
                        const Vector<double>& s,
                        Vector<double>& x) const
    {
      // Local coordinates in bulk element
      Vector<double> s_bulk(dim() + 1);
      s_bulk = local_coordinate_in_bulk(s);
 
      // Get Eulerian position vector
      bulk_element_pt()->interpolated_x(t, s_bulk, x);
    }
 
    double interpolated_dxdt(const Vector<double>& s,
                             const unsigned& i,
                             const unsigned& t)
    {
      // Local coordinates in bulk element
      Vector<double> s_bulk(dim() + 1);
      s_bulk = local_coordinate_in_bulk(s);
 
      // Return Eulerian coordinate as computed by bulk
      return bulk_element_pt()->interpolated_dxdt(s_bulk, i, t);
    }
 
    void interpolated_dxdt(const Vector<double>& s,
                           const unsigned& t,
                           Vector<double>& dxdt)
    {
      // Local coordinates in bulk element
      Vector<double> s_bulk(dim() + 1);
      s_bulk = local_coordinate_in_bulk(s);
 
      // Get Eulerian position vector
      bulk_element_pt()->interpolated_dxdt(s_bulk, t, dxdt);
    }
 
    int& normal_sign()
    {
      return Normal_sign;
    }
 
    int normal_sign() const
    {
      return Normal_sign;
    }
 
    int& face_index()
    {
      return Face_index;
    }
 
    int face_index() const
    {
      return Face_index;
    }
 
    const Vector<double>* tangent_direction_pt() const
    {
      return Tangent_direction_pt;
    }
 
    void set_tangent_direction(Vector<double>* tangent_direction_pt)
    {
#ifdef PARANOID
      // Check that tangent_direction_pt is not null.
      if (tangent_direction_pt == 0)
      {
        std::ostringstream error_message;
        error_message << "The pointer tangent_direction_pt is null.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that the vector is the correct size.
      // The size of the tangent vector.
      unsigned tang_dir_size = tangent_direction_pt->size();
      unsigned spatial_dimension = this->nodal_dimension();
      if (tang_dir_size != spatial_dimension)
      {
        std::ostringstream error_message;
        error_message << "The tangent direction vector has size "
                      << tang_dir_size << "\n"
                      << "but this element has spatial dimension "
                      << spatial_dimension << ".\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      if (tang_dir_size == 2)
      {
        std::ostringstream warning_message;
        warning_message
          << "The spatial dimension is " << spatial_dimension << ".\n"
          << "I do not need a tangent direction vector to create \n"
          << "continuous tangent vectors in two spatial dimensions.";
        OomphLibWarning(warning_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Set the direction vector for the tangent.
      Tangent_direction_pt = tangent_direction_pt;
    }
 
    void turn_on_warning_for_discontinuous_tangent()
    {
      FaceElement::Ignore_discontinuous_tangent_warning = true;
    }
 
    void turn_off_warning_for_discontinuous_tangent()
    {
      FaceElement::Ignore_discontinuous_tangent_warning = false;
    }
 
    void continuous_tangent_and_outer_unit_normal(
      const Vector<double>& s,
      Vector<Vector<double>>& tang_vec,
      Vector<double>& unit_normal) const;
 
    void continuous_tangent_and_outer_unit_normal(
      const unsigned& ipt,
      Vector<Vector<double>>& tang_vec,
      Vector<double>& unit_normal) const;
 
    void outer_unit_normal(const Vector<double>& s,
                           Vector<double>& unit_normal) const;
 
    void outer_unit_normal(const unsigned& ipt,
                           Vector<double>& unit_normal) const;
 
    FiniteElement*& bulk_element_pt()
    {
      return Bulk_element_pt;
    }
 
 
    FiniteElement* bulk_element_pt() const
    {
      return Bulk_element_pt;
    }
 
    // Clang specific pragma's
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Woverloaded-virtual"
#endif
 
    CoordinateMappingFctPt& face_to_bulk_coordinate_fct_pt()
    {
      return Face_to_bulk_coordinate_fct_pt;
    }
 
    CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt() const
    {
      return Face_to_bulk_coordinate_fct_pt;
    }
 
 
    BulkCoordinateDerivativesFctPt& bulk_coordinate_derivatives_fct_pt()
    {
      return Bulk_coordinate_derivatives_fct_pt;
    }
 
    BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt() const
    {
      return Bulk_coordinate_derivatives_fct_pt;
    }
 
#ifdef __clang__
#pragma clang diagnostic pop
#endif
 
    Vector<double> local_coordinate_in_bulk(const Vector<double>& s) const;
 
    void get_local_coordinate_in_bulk(const Vector<double>& s,
                                      Vector<double>& s_bulk) const;
 
    void get_ds_bulk_ds_face(const Vector<double>& s,
                             DenseMatrix<double>& dsbulk_dsface,
                             unsigned& interior_direction) const;
 
    unsigned& bulk_position_type(const unsigned& i)
    {
      return Bulk_position_type[i];
    }
 
    const unsigned& bulk_position_type(const unsigned& i) const
    {
      return Bulk_position_type[i];
    }
 
    void bulk_node_number_resize(const unsigned& i)
    {
      Bulk_node_number.resize(i);
    }
 
    unsigned& bulk_node_number(const unsigned& n)
    {
      return Bulk_node_number[n];
    }
 
    const unsigned& bulk_node_number(const unsigned& n) const
    {
      return Bulk_node_number[n];
    }
 
    void bulk_position_type_resize(const unsigned& i)
    {
      Bulk_position_type.resize(i);
    }
 
    unsigned& nbulk_value(const unsigned& n)
    {
      return Nbulk_value[n];
    }
 
    unsigned nbulk_value(const unsigned& n) const
    {
      return Nbulk_value[n];
    }
 
    void nbulk_value_resize(const unsigned& i)
    {
      Nbulk_value.resize(i);
    }
 
    void resize_nodes(Vector<unsigned>& nadditional_data_values)
    {
      // Locally cache the number of node
      unsigned n_node = nnode();
      // Resize the storage for values at the nodes:
      for (unsigned l = 0; l < n_node; l++)
      {
        // Find number of values stored at the node
        unsigned Initial_Nvalue = node_pt(l)->nvalue();
        // Read out the number of additional values
        unsigned Nadditional = nadditional_data_values[l];
        // If the node has not already been resized, resize it
        if ((Initial_Nvalue == Nbulk_value[l]) && (Nadditional > 0))
        {
          // Resize the node according to the number of additional values
          node_pt(l)->resize(Nbulk_value[l] + Nadditional);
        }
      } // End of loop over nodes
    }
 
    void output_zeta(std::ostream& outfile, const unsigned& nplot);
  };
 
 
  //========================================================================
  //========================================================================
  class SolidFaceElement : public virtual FaceElement,
                           public virtual SolidFiniteElement
  {
  public:
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
 
 
    double interpolated_xi(const Vector<double>& s, const unsigned& i) const
    {
      // Local coordinates in bulk element
      Vector<double> s_bulk(dim() + 1);
      s_bulk = local_coordinate_in_bulk(s);
 
      // Return Lagrangian coordinate as computed by bulk
      return dynamic_cast<SolidFiniteElement*>(bulk_element_pt())
        ->interpolated_xi(s_bulk, i);
    }
 
 
    void interpolated_xi(const Vector<double>& s, Vector<double>& xi) const
    {
      // Local coordinates in bulk element
      Vector<double> s_bulk(dim() + 1);
      s_bulk = local_coordinate_in_bulk(s);
 
      // Get Lagrangian position vector
      dynamic_cast<SolidFiniteElement*>(bulk_element_pt())
        ->interpolated_xi(s_bulk, xi);
    }
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  class SolidPointElement : public virtual SolidFiniteElement,
                            public virtual PointElement
  {
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry
  {
  };
 
 
 
 
  //======================================================================
  //======================================================================
  template<class ELEMENT>
  class DummyFaceElement : public virtual FaceGeometry<ELEMENT>,
                           public virtual FaceElement
  {
  public:
    DummyFaceElement(FiniteElement* const& element_pt, const int& face_index)
      : FaceGeometry<ELEMENT>(), FaceElement()
    {
      // Attach the geometrical information to the element. N.B. This function
      // also assigns nbulk_value from the required_nvalue of the bulk element
      element_pt->build_face_element(face_index, this);
    }
 
 
    DummyFaceElement() : FaceGeometry<ELEMENT>(), FaceElement() {}
 
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
 
    void output(std::ostream& outfile)
    {
      outfile << "ZONE" << std::endl;
      unsigned nnod = nnode();
      for (unsigned j = 0; j < nnod; j++)
      {
        Node* nod_pt = node_pt(j);
        unsigned dim = nod_pt->ndim();
        for (unsigned i = 0; i < dim; i++)
        {
          outfile << nod_pt->x(i) << " ";
        }
        outfile << std::endl;
      }
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      FiniteElement::output(outfile, n_plot);
    }
 
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      FiniteElement::output(file_pt, n_plot);
    }
  };
 
 
 
  //=========================================================================
  //=========================================================================
  class ElementWithDragFunction
  {
  public:
    ElementWithDragFunction() {}
 
    virtual ~ElementWithDragFunction() {}
 
    virtual void get_drag_and_torque(Vector<double>& drag_force,
                                     Vector<double>& drag_torque) = 0;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  template<class ELEMENT>
  class FreeStandingFaceElement : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FreeStandingFaceElement()
      : FaceGeometry<ELEMENT>(), Boundary_number_in_bulk_mesh(0)
    {
      // Check whether things have been set
#ifdef PARANOID
      Boundary_number_in_bulk_mesh_has_been_set = false;
#endif
    }
 
    inline const unsigned& boundary_number_in_bulk_mesh() const
    {
      return Boundary_number_in_bulk_mesh;
    }
 
 
    inline void set_boundary_number_in_bulk_mesh(const unsigned& b)
    {
      Boundary_number_in_bulk_mesh = b;
#ifdef PARANOID
      Boundary_number_in_bulk_mesh_has_been_set = true;
#endif
    }
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      // Vector in which to hold the intrinsic coordinate
      Vector<double> zeta(this->dim());
 
      // Get the k-th generalised boundary coordinate at node n
      this->node_pt(n)->get_coordinates_on_boundary(
        Boundary_number_in_bulk_mesh, k, zeta);
 
      // Return the individual coordinate
      return zeta[i];
    }
 
  protected:
    unsigned Boundary_number_in_bulk_mesh;
 
#ifdef PARANOID
 
    bool Boundary_number_in_bulk_mesh_has_been_set;
 
#endif
  };
 
 
 
 
  //======================================================================
  //======================================================================
  class SolidElementWithDiagonalMassMatrix
  {
  public:
    SolidElementWithDiagonalMassMatrix() {}
 
    virtual ~SolidElementWithDiagonalMassMatrix() {}
 
    SolidElementWithDiagonalMassMatrix(
      const SolidElementWithDiagonalMassMatrix&) = delete;
 
    void operator=(const SolidElementWithDiagonalMassMatrix&) = delete;
 
    virtual void get_mass_matrix_diagonal(Vector<double>& mass_diag) = 0;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  class NavierStokesElementWithDiagonalMassMatrices
  {
  public:
    NavierStokesElementWithDiagonalMassMatrices() {}
 
    virtual ~NavierStokesElementWithDiagonalMassMatrices() {}
 
    NavierStokesElementWithDiagonalMassMatrices(
      const NavierStokesElementWithDiagonalMassMatrices&) = delete;
 
    void operator=(const NavierStokesElementWithDiagonalMassMatrices&) = delete;
 
    virtual void get_pressure_and_velocity_mass_matrix_diagonal(
      Vector<double>& press_mass_diag,
      Vector<double>& veloc_mass_diag,
      const unsigned& which_one = 0) = 0;
  };
 
 
 
} // namespace oomph
 
#endif