elements.h

Go to the documentation of this file.

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


//========================================================================
/// \short A Generalised Element class.
///
/// The main components of a GeneralisedElement are:
/// - pointers to its internal and external Data, which are stored together
///   in a single array so that we need only store one pointer. 
///   The internal Data are placed at the beginning
///   of the array and the external Data are stored at the end.
/// - a pointer to a global Time
/// - a lookup table that establishes the relation between local
///   and global equation numbers.
/// 
/// We also provide interfaces for functions that compute the
/// element's Jacobian matrix and/or the Vector of residuals.
/// In addition, an interface that returns a mass matrix --- the matrix
/// of terms that multiply any time derivatives in the problem --- is
/// also provided to permit explicit time-stepping and the solution
/// of the generalised eigenproblems. 
//========================================================================
class GeneralisedElement
{
  private:
  
 /// \short Storage for the global equation numbers represented by the 
 /// degrees of freedom in the element.
 unsigned long *Eqn_number;

 /// \short Storage for array of pointers to degrees of freedom ordered
 /// by local equation number. This information is not needed, except in
 /// continuation, bifurcation tracking and periodic orbit calculations.
 /// It is not set up unless the control flag 
 /// Problem::Store_local_dof_pts_in_elements = true
 double **Dof_pt;

 /// \short Storage for pointers to internal and external data.
 /// The data is of the same type and so can be addressed by 
 /// a single pointer. The idea is that the array is of 
 /// total size Ninternal_data + Nexternal_data.
 /// The internal data are stored at the beginning of the array
 /// and the external data are stored at the end of the array.
 Data** Data_pt;

 /// \short Pointer to array storage for local equation numbers associated
 /// with internal and external data. Again, we save storage by using 
 /// a single pointer to access this information. The first index of the
 /// array is of dimension Nineternal_data + Nexternal_data and the second
 /// index varies with the number of values stored at the data object.
 /// In the most general case, however, the scheme is as memory efficient
 /// as possible.
 int** Data_local_eqn;

 /// Number of degrees of freedom 
 unsigned Ndof;

 /// Number of internal data
 unsigned Ninternal_data;

 /// Number of external data
 unsigned Nexternal_data;

 /// \short Storage for boolean flags of size Ninternal_data + Nexternal_data
 /// that correspond to the data used in the element. The flags
 /// indicate whether the particular
 /// internal or external data should be included in the general 
 /// finite-difference loop in fill_in_jacobian_from_internal_by_fd() or
 /// fill_in_jacobian_from_external_by_fd(). The default is that all
 /// data WILL be included in the finite-difference loop, but in many 
 /// circumstances it is possible to treat certain (external) data 
 /// analytically. The use of an STL vector is optimal for memory
 /// use because the booleans are represented as single-bits.
 std::vector<bool> Data_fd;

  protected:

#ifdef OOMPH_HAS_MPI

 /// \short Non-halo processor ID for Data; -1 if it's not a halo.
 int Non_halo_proc_ID;

 /// Does this element need to be kept as a halo element during a distribute?
 bool Must_be_kept_as_halo;

#endif

 /// \short Add a (pointer to an) internal data object to the element and
 /// return the index required to obtain it from the access
 /// function \c internal_data_pt(). The boolean indicates whether
 /// the datum should be included in the general finite-difference loop
 /// when calculating the jacobian. The default value is true, i.e.
 /// the data will be included in the finite differencing.
 unsigned add_internal_data(Data* const &data_pt, const bool &fd=true);
 
 
 /// \short Return the status of the boolean flag indicating whether
 /// the internal data is included in the finite difference loop
 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];
  }


 /// \short Set the boolean flag to exclude the internal datum from
 /// the finite difference loop when computing the jacobian matrix
 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;
  }

 /// \short Set the boolean flag to include the internal datum in
 /// the finite difference loop when computing the jacobian matrix
 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;
  }

 /// \short Clear the storage for the global equation numbers
 /// and pointers to dofs (if stored)
 void clear_global_eqn_numbers() 
 {delete[] Eqn_number; Eqn_number=0; delete[] Dof_pt; Dof_pt=0; Ndof=0;}

 /// \short Add the contents of the queue global_eqn_numbers 
 /// to the local storage for the local-to-global translation scheme.
 /// It is essential that the entries in the queue are added IN ORDER
 /// i.e. from the front.
 void add_global_eqn_numbers(std::deque<unsigned long> 
                             const &global_eqn_numbers,
                             std::deque<double*>
                             const &global_dof_pt);
 
 /// \short Empty dense matrix used as a dummy argument to combined
 /// residual and jacobian functions in the case when only the residuals
 /// are being assembled
 static DenseMatrix<double> Dummy_matrix;
 
 /// \short Static storage for deque used to add_global_equation_numbers
 /// when pointers to the dofs in each element are not required
 static std::deque<double*> Dof_pt_deque;

 /// \short Assign the local equation numbers for the internal 
 /// and external Data
 /// This must be called after the global equation numbers have all been
 /// assigned. It is virtual so that it can be overloaded by
 /// ElementWithExternalElements so that any external data from the
 /// external elements in included in the numbering scheme.
 /// If the boolean argument is true then pointers to the dofs will be 
 /// stored in Dof_pt
 virtual void assign_internal_and_external_local_eqn_numbers(
  const bool &store_local_dof_pt);

 /// \short Assign all the local equation numbering schemes that can
 /// be applied generically for the element. In most cases, this is the
 /// function that will be overloaded by inherited classes. It is required
 /// to ensure that assign_additional_local_eqn_numbers() can always be 
 /// called after ALL other local equation numbering has been performed.
 /// The default for the GeneralisedElement is simply to call internal
 /// and external local equation numbering.
 /// If the boolean argument is true then pointers to the dofs will be stored
 /// in 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);
  }
 
 /// \short Setup any additional look-up schemes for local equation numbers.
 /// Examples of use include using local storage to refer to explicit
 /// degrees of freedom. The additional memory cost of such storage may
 /// or may not be offset by fast local access.
 virtual void assign_additional_local_eqn_numbers() {}

 /// \short Return the local equation number corresponding to the j-th
 /// value stored at the i-th internal data
 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];
  }

 /// \short Return the local equation number corresponding to the j-th
 /// value stored at the i-th external data
 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];
  }
 
 /// \short Add the elemental contribution to the residuals vector. Note that
 /// this function will NOT initialise the residuals vector. It must be
 /// called after the residuals vector has been initialised to zero.
 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);
  }
 
 /// \short Calculate the contributions to the jacobian from the internal
 /// degrees of freedom using finite differences.
 /// This version of the function assumes that the residuals vector has 
 /// already been calculated. If the boolean argument is true, the finite
 /// differencing will be performed for all internal data, irrespective of
 /// the information in Data_fd. The default value (false) 
 /// uses the information in Data_fd to selectively difference only
 /// certain data.
 void fill_in_jacobian_from_internal_by_fd(Vector<double> &residuals,
                                           DenseMatrix<double> &jacobian,
                                           const bool &fd_all_data=false);

 /// \short Calculate the contributions to the jacobian from the internal
 /// degrees of freedom using finite differences. This version computes
 /// the residuals vector before calculating the jacobian terms.
 /// If the boolean argument is true, the finite
 /// differencing will be performed for all internal data, irrespective of
 /// the information in Data_fd. The default value (false) 
 /// uses the information in Data_fd to selectively difference only
 /// certain data.
 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);
  }


 /// \short Calculate the contributions to the jacobian from the external
 /// degrees of freedom using finite differences.
 /// This version of the function assumes that the residuals vector has 
 /// already been calculated.
 /// If the boolean argument is true, the finite
 /// differencing will be performed for all external data, irrespective of
 /// the information in Data_fd. The default value (false) 
 /// uses the information in Data_fd to selectively difference only
 /// certain data.
 void fill_in_jacobian_from_external_by_fd(Vector<double> &residuals,
                                           DenseMatrix<double> &jacobian,
                                           const bool &fd_all_data=false);
 
  
 /// \short Calculate the contributions to the jacobian from the external
 /// degrees of freedom using finite differences. This version computes
 /// the residuals vector before calculating the jacobian terms.
 /// If the boolean argument is true, the finite
 /// differencing will be performed for all internal data, irrespective of
 /// the information in Data_fd. The default value (false) 
 /// uses the information in Data_fd to selectively difference only
 /// certain data.
 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);
  }

 /// \short Function that is called before the finite differencing of
 /// any internal data. This may be overloaded to update any slaved
 /// data before finite differencing takes place.
 virtual inline void update_before_internal_fd() { }

 /// \short Function that is call after the finite differencing of
 /// the internal data. This may be overloaded to reset any slaved
 /// variables that may have changed during the finite differencing.
 virtual inline void reset_after_internal_fd() { }

 /// \short Function called within the finite difference loop for 
 /// internal data after a change in any values in the i-th 
 /// internal data object.
 virtual inline void update_in_internal_fd(const unsigned &i) { }

 /// \short Function called within the finite difference loop for
 /// internal data after the values in the i-th external data object
 /// are reset. The default behaviour is to call the update function.
 virtual inline void reset_in_internal_fd(const unsigned &i)
  {update_in_internal_fd(i);}


 /// \short Function that is called before the finite differencing of
 /// any external data. This may be overloaded to update any slaved
 /// data before finite differencing takes place.
 virtual inline void update_before_external_fd() { }

 /// \short Function that is call after the finite differencing of
 /// the external data. This may be overloaded to reset any slaved
 /// variables that may have changed during the finite differencing.
 virtual inline void reset_after_external_fd() { }

 /// \short Function called within the finite difference loop for 
 /// external data after a change in any values in the i-th 
 /// external data object
 virtual inline void update_in_external_fd(const unsigned &i) { }

 /// \short Function called within the finite difference loop for
 /// external data after the values in the i-th external data object
 /// are reset. The default behaviour is to call the update function.
 virtual inline void reset_in_external_fd(const unsigned &i)
  {update_in_external_fd(i);}

 /// \short Add the elemental contribution to the jacobian matrix.
 /// and the residuals vector. Note that
 /// this function will NOT initialise the residuals vector or the jacobian
 /// matrix. It must be called after the residuals vector and 
 /// jacobian matrix have been initialised to zero. The default
 /// is to use finite differences to calculate the jacobian
 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);
  }

 /// \short Add the elemental contribution to the mass matrix matrix.
 /// and the residuals vector. Note that
 /// this function should NOT initialise the residuals vector or the mass
 /// matrix. It must be called after the residuals vector and 
 /// jacobian matrix have been initialised to zero. The default
 /// is deliberately broken
 virtual void fill_in_contribution_to_mass_matrix(
  Vector<double> &residuals, DenseMatrix<double> &mass_matrix);
 
 /// \short Add the elemental contribution to the jacobian matrix,
 /// mass matrix and the residuals vector. Note that
 /// this function should NOT initialise any entries.
 /// It must be called after the residuals vector and 
 /// matrices have been initialised to zero. 
 virtual void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals,
  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix);

 /// \short Add the elemental contribution to the derivatives of 
 /// the residuals with respect to a parameter. This function should
 /// NOT initialise any entries and must be called after the entries
 /// have been initialised to zero
 /// The default implementation is to use finite differences to calculate 
 /// the derivatives.
 virtual void fill_in_contribution_to_dresiduals_dparameter(
  double* const &parameter_pt,
  Vector<double> &dres_dparam);


 /// \short Add the elemental contribution to the derivatives of 
 /// the elemental Jacobian matrix 
 /// and residuals with respect to a parameter. This function should
 /// NOT initialise any entries and must be called after the entries
 /// have been initialised to zero
 /// The default implementation is to use finite differences to calculate 
 /// the derivatives.
 virtual void fill_in_contribution_to_djacobian_dparameter(
  double* const &parameter_pt,
  Vector<double> &dres_dparam, 
  DenseMatrix<double> &djac_dparam);

 /// \short Add the elemental contribution to the derivative of the
 /// jacobian matrix,
 /// mass matrix and the residuals vector with respect to the 
 /// passed parameter. Note that
 /// this function should NOT initialise any entries.
 /// It must be called after the residuals vector and 
 /// matrices have been initialised to zero. 
 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);


 /// \short Fill in contribution to the product of the Hessian 
 /// (derivative of Jacobian with
 /// respect to all variables) an eigenvector, Y, and 
 /// other specified vectors, C
 /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
 virtual void fill_in_contribution_to_hessian_vector_products(
  Vector<double> const &Y,
  DenseMatrix<double> const &C,
  DenseMatrix<double> &product);

 /// \short Fill in the contribution to the inner products between given
 /// pairs of history values
 virtual void fill_in_contribution_to_inner_products(
  Vector<std::pair<unsigned,unsigned> > const &history_index,
  Vector<double> &inner_product);

 /// \short Fill in the contributions to the vectors that when taken
 /// as dot product with other history values give the inner product
 /// over the element
 virtual void fill_in_contribution_to_inner_product_vectors(
  Vector<unsigned> const &history_index,
  Vector<Vector<double> >  &inner_product_vector);

public:

 /// \short Constructor: Initialise all pointers and all values to zero
  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 destructor to clean up any memory allocated by the object.
 virtual ~GeneralisedElement();

 /// Broken copy constructor
 GeneralisedElement(const GeneralisedElement&) 
  { 
   BrokenCopy::broken_copy("GeneralisedElement");
  } 
 
 /// Broken assignment operator
 void operator=(const GeneralisedElement&) 
  {
   BrokenCopy::broken_assign("GeneralisedElement");
  }
 
 /// Return a pointer to i-th internal data object.
 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];
  }

 /// Return a pointer to i-th internal data object (const version)
 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];
  }
 

 /// Return a pointer to i-th external data object.
 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];
  }

 /// Return a pointer to i-th external data object (const version)
 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];
  }

 /// \short Static boolean to suppress warnings about repeated internal
 /// data. Defaults to false.
 static bool Suppress_warning_about_repeated_internal_data;

 /// \short Static boolean to suppress warnings about repeated external
 /// data. Defaults to true.
 static bool Suppress_warning_about_repeated_external_data;

 /// \short Return the global equation number corresponding to the 
 /// ieqn_local-th local equation number
 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];
  }


 ///\short Return the local equation number corresponding to the ieqn_global-th
 ///global equation number. Returns minus one (-1) if there is no
 ///local degree of freedom corresponding to the chosen global equation
 ///number
 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;
  }


 /// Add a (pointer to an) external data object to the element and return its
 /// index (i.e. the index required to obtain it from 
 /// the access function \c external_data_pt(...). The optional boolean
 /// flag indicates whether the data should be included in the
 /// general finite-difference loop when calculating the jacobian. 
 /// The default value is true, i.e. the data will be included in the 
 /// finite-differencing
 unsigned add_external_data(Data* const &data_pt, const bool &fd=true); 
 
 /// \short Return the status of the boolean flag indicating whether
 /// the external data is included in the finite difference loop
 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];
  }



 /// \short Set the boolean flag to exclude the external datum from the 
 /// the finite difference loop when computing the jacobian matrix
 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;
  }

 /// \short Set the boolean flag to include the external datum in the 
 /// the finite difference loop when computing the jacobian matrix
 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;
  }

 /// Flush all external data
 void flush_external_data();
 
 /// Flush the object addressed by data_pt from the external data array
 void flush_external_data(Data* const &data_pt);

 /// Return the number of internal data objects.
 unsigned ninternal_data() const {return Ninternal_data;}

 /// Return the number of external data objects.
 unsigned nexternal_data() const {return Nexternal_data;}

 /// Return the number of equations/dofs in the element.
 unsigned ndof() const {return Ndof;}

 ///\short Return the vector of dof values at time level t
 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];
   }
 }

 ///\short Return the vector of pointers to dof values
 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];
   }
 }

 

 /// \short Set the timestepper associated with the i-th internal data
 /// object
 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);}

 /// \short Assign the global equation numbers to the internal Data. 
 /// The arguments are the current highest global equation number 
 /// (which will be incremented) and a Vector of pointers to the 
 /// global variables (to which any unpinned values in the internal Data are
 /// added). 
 void assign_internal_eqn_numbers(unsigned long &global_number, 
                                  Vector<double *> &Dof_pt);

 /// \short Function to describe the dofs of the element. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 void describe_dofs(std::ostream& out,const std::string& current_string) const;
 
 /// \short Function to describe the local dofs of the element. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 virtual void describe_local_dofs(std::ostream& out,
                                  const std::string& current_string) const;

 /// \short Add pointers to the internal data values to map indexed
 /// by the global equation number.
 void add_internal_value_pt_to_map(std::map<unsigned,double*> &
                                   map_of_value_pt);

#ifdef OOMPH_HAS_MPI
 /// \short Add all internal data and time history values to the vector in 
 /// the internal storage order
 void add_internal_data_values_to_vector(Vector<double> &vector_of_values);

 /// \short Read all internal data and time history values from the vector
 /// starting from index. On return the index will be
 /// set to the value at the end of the data that has been read in
 void read_internal_data_values_from_vector(
  const Vector<double> & vector_of_values, unsigned &index);

 /// \short Add all equation numbers associated with internal data 
 /// to the vector in the internal storage order
 void add_internal_eqn_numbers_to_vector(Vector<long> &vector_of_eqn_numbers);

 /// \short Read all equation numbers associated with internal data
 /// from the vector
 /// starting from index. On return the index will be
 /// set to the value at the end of the data that has been read in
 void read_internal_eqn_numbers_from_vector(
  const Vector<long> & vector_of_eqn_numbers, unsigned &index);
 
#endif

 /// \short Setup the arrays of local equation numbers for the element. 
 /// If the optional boolean argument is true, then pointers to the associated
 /// degrees of freedom are stored locally in the array Dof_pt
 virtual void assign_local_eqn_numbers(const bool &store_local_dof_pt);
 
 /// \short Complete the setup of any additional dependencies
 /// that the element may have. Empty virtual function that may be
 /// overloaded for specific derived elements. Used, e.g., for elements 
 /// with algebraic node update functions to determine the "geometric 
 /// Data", i.e. the Data that affects the element's shape.
 /// This function is called (for all elements) at the very beginning of the
 /// equation numbering procedure to ensure that all dependencies
 /// are accounted for.
 virtual void complete_setup_of_dependencies(){}

 /// \short Calculate the vector of residuals of the equations in the element.
 /// By default initialise the vector to zero and then call the 
 /// fill_in_contribution_to_residuals() function. Note that this entire
 /// function can be overloaded if desired.
 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);
  }

 /// \short Calculate the elemental Jacobian matrix "d equation / d variable".
 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);
 }

 /// \short Calculate the residuals and the elemental "mass" matrix, the
 /// matrix that multiplies the time derivative terms in a problem.
 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);
  }

 /// \short Calculate the residuals and jacobian and elemental "mass" matrix,
 /// the matrix that multiplies the time derivative terms.
 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);
  }


 /// \short Calculate the derivatives of the residuals with respect to
 /// a parameter
 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);
  }

 /// \short Calculate the derivatives of the elemental Jacobian matrix 
 /// and residuals with respect to a parameter
 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);
  }

 /// \short Calculate the derivatives of the elemental Jacobian matrix 
 /// mass matrix and residuals with respect to a parameter
 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);
  }


 /// \short Calculate the product of the Hessian (derivative of Jacobian with
 /// respect to all variables) an eigenvector, Y, and 
 /// other specified vectors, C
 /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
 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);
   ///Add the elemental contribution to the product
   this->fill_in_contribution_to_hessian_vector_products(Y,C,product);
  }
 
 /// \short Return the vector of inner product of the given pairs of 
 /// history values
 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);
 }

 /// \short  Compute the vectors that when taken as a dot product with
 /// other history values give the inner product over the element.
 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);
 }


 /// \short Self-test: Have all internal values been classified as 
 /// pinned/unpinned? Return 0 if OK. 
 virtual unsigned self_test();


 /// \short Compute norm of solution -- broken virtual can be overloaded
 /// by element writer to implement whatever norm is desired for
 /// the specific element
 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

 /// \short Label the element as halo and specify processor that holds
 /// non-halo counterpart
 void set_halo(const unsigned& non_halo_proc_ID) 
 {
  Non_halo_proc_ID=non_halo_proc_ID;
 }

 /// \short Label the element as not being a halo
 void set_nonhalo() {Non_halo_proc_ID=-1;}
 
 /// \short Is this element a halo?
 bool is_halo() const {return (Non_halo_proc_ID!=-1);} 
 
 /// \short ID of processor ID that holds non-halo counterpart
 /// of halo element; negative if not a halo.
 int non_halo_proc_ID() 
 {
  return Non_halo_proc_ID;
 }
 
 /// Insist that this element be kept as a halo element during a distribute?
 void set_must_be_kept_as_halo() {Must_be_kept_as_halo = true;}

 /// \short Do not insist that this element be kept as a halo element during
 /// distribution
 void unset_must_be_kept_as_halo() {Must_be_kept_as_halo = false;}

 /// \short Test whether the element must be kept as a halo element
 bool must_be_kept_as_halo() const {return Must_be_kept_as_halo;}

#endif

 /// \short Double used for the default finite difference step in elemental
 /// jacobian calculations
 static double Default_fd_jacobian_step;

 /// \short The number of types of degrees of freedom in this element
 /// are sub-divided into
 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);
  }

 /// \short Create a list of pairs for the unknowns that this element
 /// is "in charge of" -- ignore any unknowns associated with external
 /// \c Data. The first entry in each pair must contain the global equation 
 /// number of the unknown, while the second one contains the number
 /// of the DOF type that this unknown is associated with. 
 /// (The function can obviously only be called if the equation numbering
 /// scheme has been set up.)
 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);
  }

};

 /// Enumeration a finite element's geometry "type". Either "Q" (square,
 /// cubeoid like) or "T" (triangle, tetrahedron).
 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;

//===========================================================================
/// \short Helper namespace for tolerances and iterations within the Newton
/// method used in the locate_zeta function in FiniteElements.
//===========================================================================
namespace Locate_zeta_helpers
{
 /// Convergence tolerance for the Newton solver
 extern double Newton_tolerance;

 /// Maximum number of Newton iterations
 extern unsigned Max_newton_iterations;
 
 /// \short Multiplier for (zeta-based) outer radius of element used for
 /// deciding that point is outside element. Set to negative value
 /// to suppress test.
 extern double Radius_multiplier_for_fast_exit_from_locate_zeta;

 /// Number of points along one dimension of each element used
 /// to create coordinates within the element in order to see
 /// which has the smallest initial residual (and is therefore used
 /// as the initial guess in the Newton method for locate_zeta)
 extern unsigned N_local_points;

}


 /// \short Typedef for the function that translates the face coordinate
 /// to the coordinate in the bulk element
 typedef void (*CoordinateMappingFctPt)(const Vector<double> &s,
                                        Vector<double> &s_bulk);

 /// \short Typedef for the function that returns the partial derivative
 /// of the local coordinates in the bulk element
 /// with respect to the coordinates along the face.
 /// In addition this function returns an index of one of the 
 /// bulk local coordinates that varies away from the edge
 typedef void (*BulkCoordinateDerivativesFctPt)(const Vector<double> &s,
                                                DenseMatrix<double> &ds_bulk_dsface,
                                                unsigned &interior_direction);


//========================================================================
/// \short A general Finite Element class.
///
/// The main components of a FiniteElement are:
/// - pointers to its Nodes
/// - pointers to its internal Data (inherited from GeneralisedElement)
/// - pointers to its external Data (inherited from GeneralisedElement)
/// - a pointer to a spatial integration scheme
/// - a pointer to the global  Time object (inherited from GeneralisedElement)
/// - a lookup table which establishes the relation between local
///   and global equation numbers (inherited from GeneralisedElement)
///
/// We also provide interfaces for functions that compute the
/// element's Jacobian matrix and/or the Vector of residuals 
/// (inherited from GeneralisedElement) plus various output routines. 
//========================================================================
class FiniteElement : public virtual GeneralisedElement, public GeomObject
{
  private:

 /// Pointer to the spatial integration scheme
 Integral* Integral_pt;
 
 /// Storage for pointers to the nodes in the element
 Node** Node_pt;

 /// \short Storage for the local equation numbers associated with
 /// the values stored at the nodes
 int** Nodal_local_eqn;
   
 /// Number of nodes in the element
 unsigned Nnode;

 /// \short The spatial dimension of the element, i.e. the number
 /// of local coordinates used to parametrize it.
 unsigned Elemental_dimension;

 /// \short The spatial dimension of the nodes in the element.
 /// We assume that nodes have the same spatial dimension, because
 /// we cannot think of any "real" problems for which that would not
 /// be the case.
 unsigned Nodal_dimension;
 
 /// \short The number of coordinate types required to interpolate
 /// the element's geometry between the nodes. For Lagrange elements 
 /// it is 1 (the default). It must be over-ridden by using
 /// the set_nposition_type() function in the constructors of elements
 /// that use generalised coordinate, e.g. for 1D Hermite elements
 /// Nnodal_position_types =2.
 unsigned Nnodal_position_type;

  protected:

 /// \short Assemble the jacobian matrix for the mapping from local
 /// to Eulerian coordinates, given the derivatives of the shape function
 /// w.r.t the local coordinates.
 virtual void assemble_local_to_eulerian_jacobian(
  const DShape &dpsids, DenseMatrix<double> &jacobian) const;
 
 /// \short Assemble the the "jacobian" matrix of second derivatives of the
 /// mapping from local to Eulerian coordinates, given
 /// the second derivatives of the shape functions w.r.t. local coordinates.
 virtual void assemble_local_to_eulerian_jacobian2(
  const DShape &d2psids, DenseMatrix<double> &jacobian2) const;

 /// \short Assemble the covariant Eulerian base vectors, assuming that
 /// the derivatives of the shape functions with respect to the local 
 /// coordinates have already been constructed.
 virtual void assemble_eulerian_base_vectors(
  const DShape &dpsids, DenseMatrix<double> &interpolated_G) const;
                                   
 /// \short Default return value for required_nvalue(n) which gives the number 
 /// of "data" values required by the element at node n; for example, solving
 /// a Poisson equation would required only one "data" value at each node. The
 /// defaults is set to zero, because a general element is problem-less. 
 static const unsigned Default_Initial_Nvalue;

 /// \short Default value for the tolerance to be used when locating nodes
 /// via local coordinates
 static const double Node_location_tolerance;
 
public:

 /// \short Set the dimension of the element and initially set
 /// the dimension of the nodes to be the same as the dimension of the
 /// element.
 void set_dimension(const unsigned &dim)
  {Elemental_dimension = dim; Nodal_dimension = dim;}
 
 /// \short Set the dimension of the nodes in the element. This will
 /// typically only be required when constructing FaceElements or
 /// in beam and shell type elements where a lower dimensional surface
 /// is embedded in a higher dimensional space.
 void set_nodal_dimension(const unsigned &nodal_dim)
  {Nodal_dimension = nodal_dim;}
 
 /// \short Set the number of types required to interpolate the coordinate
 void set_nnodal_position_type(const unsigned &nposition_type)
  {Nnodal_position_type = nposition_type;}


 /// \short Set the number of nodes in the element to n, by resizing
 /// the storage for pointers to the Node objects.
 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;}
  }

 /// \short Return the local equation number corresponding to the i-th
 /// value at the n-th local node.
 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];
  }



 /// \short Compute the geometric shape functions (psi) at integration point
 /// ipt. Return the determinant of the jacobian of the mapping (detJ).
 /// Additionally calculate the derivatives of "detJ" w.r.t. the 
 /// nodal coordinates.
 double dJ_eulerian_at_knot(const unsigned &ipt,
                            Shape &psi, 
                            DenseMatrix<double> &djacobian_dX) const;
  protected: 

 /// \short Static array that holds the number of second derivatives
 /// as a function of the dimension of the element
 static const unsigned N2deriv[];

 /// \short Take the matrix passed as jacobian and return its inverse in
 /// inverse_jacobian. This function is templated by the dimension of the 
 /// element because matrix inversion cannot be written efficiently in a 
 /// generic manner.
 template<unsigned DIM>
 double invert_jacobian(const DenseMatrix<double> &jacobian,
                        DenseMatrix<double> &inverse_jacobian) const;
 

 /// \short A template-free interface that takes the matrix passed as jacobian
 /// and return its inverse in inverse_jacobian. By default the function will
 /// use the dimension of the element to call the correct invert_jacobian(..)
 /// function. This should be overloaded for efficiency (removal of a switch
 /// statement) in specific elements.
 virtual double invert_jacobian_mapping(const DenseMatrix<double> &jacobian,
                                        DenseMatrix<double> &inverse_jacobian) 
  const;

 
 /// \short Calculate the mapping from local to Eulerian coordinates,
 /// given the derivatives of the shape functions w.r.t. local coordinates.
 /// Returns the determinant of the jacobian, the jacobian and inverse jacobian
 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);
  }
 

 /// \short Calculate the mapping from local to Eulerian coordinates,
 /// given the derivatives of the shape functions w.r.t. local coordinates,
 /// Return only the determinant of the jacobian and the inverse of the
 /// mapping (ds/dx).
 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);
  }

 /// \short Calculate the mapping from local to Eulerian coordinates given 
 /// the derivatives of the shape functions w.r.t the local coordinates.
 /// assuming that the coordinates are aligned in the direction of the local 
 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
 /// This function returns the determinant of the jacobian, the jacobian 
 /// and the inverse jacobian.
 virtual double local_to_eulerian_mapping_diagonal(
  const DShape &dpsids,DenseMatrix<double> &jacobian,
  DenseMatrix<double> &inverse_jacobian) const;

 /// \short A template-free interface that calculates the derivative of the
 /// jacobian of a mapping with respect to the nodal coordinates X_ij.
 /// To do this it requires the jacobian matrix and the derivatives of the
 /// shape functions w.r.t. the local coordinates. By default the function
 /// will use the dimension of the element to call the correct
 /// dJ_eulerian_dnodal_coordinates_templated_helper(..) function. This
 /// should be overloaded for efficiency (removal of a switch statement)
 /// in specific elements.
 virtual void dJ_eulerian_dnodal_coordinates(
  const DenseMatrix<double> &jacobian,const DShape &dpsids,
  DenseMatrix<double> &djacobian_dX) const;
 
 /// \short Calculate the derivative of the jacobian of a mapping with
 /// respect to the nodal coordinates X_ij using the jacobian matrix
 /// and the derivatives of the shape functions w.r.t. the local
 /// coordinates. This function is templated by the dimension of the
 /// element.
 template<unsigned DIM>
  void dJ_eulerian_dnodal_coordinates_templated_helper(
   const DenseMatrix<double> &jacobian,const DShape &dpsids,
   DenseMatrix<double> &djacobian_dX) const;

 /// \short A template-free interface that calculates the derivative w.r.t.
 /// the nodal coordinates \f$ X_{pq} \f$ of the derivative of the shape
 /// functions \f$ \psi_j \f$ w.r.t. the global eulerian coordinates
 /// \f$ x_i \f$. I.e. this function calculates
 /// \f[
 /// \frac{\partial}{\partial X_{pq}}
 /// \left( \frac{\partial \psi_j}{\partial x_i} \right).
 /// \f]
 /// To do this it requires the determinant of the jacobian mapping, its
 /// derivative w.r.t. the nodal coordinates \f$ X_{pq} \f$, the inverse
 /// jacobian and the derivatives of the shape functions w.r.t. the local
 /// coordinates. The result is returned as a tensor of rank four.
 ///  Numbering: 
 /// d_dpsidx_dX(p,q,j,i) = \f$ \frac{\partial}{\partial X_{pq}}
 /// \left( \frac{\partial \psi_j}{\partial x_i} \right) \f$ 
 /// By default the function will use the dimension of the element to call
 /// the correct d_dshape_eulerian_dnodal_coordinates_templated_helper(..)
 /// function. This should be overloaded for efficiency (removal of a
 /// switch statement) in specific elements.
 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;
 
 /// \short Calculate the derivative w.r.t. the nodal coordinates
 /// \f$ X_{pq} \f$ of the derivative of the shape functions w.r.t. the
 /// global eulerian coordinates \f$ x_i \f$, using the determinant of
 /// the jacobian mapping, its derivative w.r.t. the nodal coordinates
 /// \f$ X_{pq} \f$, the inverse jacobian and the derivatives of the
 /// shape functions w.r.t. the local coordinates. The result is returned
 /// as a tensor of rank four.
 ///  Numbering: 
 /// d_dpsidx_dX(p,q,j,i) = \f$ \frac{\partial}{\partial X_{pq}}
 /// \left( \frac{\partial \psi_j}{\partial x_i} \right) \f$ 
 /// This function is templated by the dimension of the element.
 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;

 /// \short Convert derivative w.r.t.local coordinates to derivatives
 /// w.r.t the coordinates used to assemble the inverse_jacobian passed
 /// in the mapping. On entry, dbasis must contain the basis function 
 /// derivatives w.r.t. the local coordinates; it will contain the 
 /// derivatives w.r.t. the new coordinates on exit. 
 /// This is virtual so that it may be overloaded if desired 
 /// for efficiency reasons.
 virtual void transform_derivatives(const DenseMatrix<double>
                                    &inverse_jacobian, DShape &dbasis) const;
 
 /// \short Convert derivative w.r.t local coordinates to derivatives
 /// w.r.t the coordinates used to assemble the inverse jacobian passed
 /// in the mapping, assuming that the coordinates are aligned in the 
 /// direction of the local coordinates. On entry dbasis must contain 
 /// the derivatives of the basis functions w.r.t. the local coordinates;
 /// it will contain the derivatives w.r.t. the new coordinates.
 /// are converted into the new  using the mapping inverse_jacobian.
 void transform_derivatives_diagonal(const DenseMatrix<double>
                                     &inverse_jacobian, DShape &dbasis) const;

 /// \short Convert derivatives and second derivatives w.r.t. local coordiantes
 /// to derivatives and second derivatives w.r.t. the coordinates used to
 /// assemble the jacobian, inverse jacobian and jacobian2 passed to the
 /// function. By default this function will call 
 /// transform_second_derivatives_template<>(...) using the dimension of 
 /// the element as the template parameter. It is virtual so that it can
 /// be overloaded by a specific element to save using a switch statement.
 /// Optionally, the element writer may wish to use the 
 /// transform_second_derivatives_diagonal<>(...) function
 /// On entry dbasis and d2basis must contain the derivatives w.r.t. the
 /// local coordinates; on exit they will be the derivatives w.r.t. the
 /// transformed coordinates.
 virtual void transform_second_derivatives(const DenseMatrix<double>
                                           &jacobian,
                                           const DenseMatrix<double> 
                                           &inverse_jacobian,
                                           const DenseMatrix<double> 
                                           &jacobian2,
                                           DShape &dbasis,
                                           DShape &d2basis) const;

 /// \short Convert derivatives and second derivatives w.r.t. local coordinates
 /// to derivatives and second derivatives w.r.t. the coordinates used to
 /// asssmble the jacobian, inverse jacobian and jacobian2 passed in the
 /// mapping. This is templated by dimension because the method of calculation
 /// varies significantly with the dimension.
 /// On entry dbasis and d2basis must contain the derivatives w.r.t. the
 /// local coordinates; on exit they will be the derivatives w.r.t. the
 /// transformed coordinates.
 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;
 
 /// \short Convert derivatives and second derivatives w.r.t. local coordinates
 /// to derivatives and second derivatives w.r.t. the coordinates used to
 /// asssmble the jacobian, inverse jacobian and jacobian2 passed in the
 /// mapping. This version of the function assumes that the local coordinates
 /// are aligned with the global coordinates, i.e. the jacobians are diagonal
 /// On entry dbasis and d2basis must contain the derivatives w.r.t. the
 /// local coordinates; on exit they will be the derivatives w.r.t. the
 /// transformed coordinates.
 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;

 /// Pointer to the element's macro element (NULL by default)
 MacroElement* Macro_elem_pt;

 /// \short Calculate the contributions to the jacobian from the nodal
 /// degrees of freedom using finite differences.
 /// This version of the function assumes that the residuals vector has 
 /// already been calculated
 virtual void fill_in_jacobian_from_nodal_by_fd(Vector<double> &residuals,
                                                DenseMatrix<double> &jacobian);
 
 /// \short Calculate the contributions to the jacobian from the nodal
 /// degrees of freedom using finite differences. This version computes
 /// the residuals vector before calculating the jacobian terms.
 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);
  }

 /// \short Function that is called before the finite differencing of
 /// any nodal data. This may be overloaded to update any slaved
 /// data before finite differencing takes place.
 virtual inline void update_before_nodal_fd() { }

 /// \short Function that is call after the finite differencing of
 /// the nodal data. This may be overloaded to reset any slaved
 /// variables that may have changed during the finite differencing.
 virtual inline void reset_after_nodal_fd() { }

 /// \short Function called within the finite difference loop for 
 /// nodal data after a change in the i-th nodal value.
 virtual inline void update_in_nodal_fd(const unsigned &i) { }

 /// \short Function called within the finite difference loop for
 /// nodal data after the i-th nodal values is reset.
 /// The default behaviour is to call the update function.
 virtual inline void reset_in_nodal_fd(const unsigned &i)
  {update_in_nodal_fd(i);}


 /// \short Add the elemental contribution to the jacobian matrix.
 /// and the residuals vector. Note that
 /// this function will NOT initialise the residuals vector or the jacobian
 /// matrix. It must be called after the residuals vector and 
 /// jacobian matrix have been initialised to zero. The default
 /// is to use finite differences to calculate the jacobian
 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:


 /// \short Function pointer for function that computes vector-valued
 /// steady "exact solution" \f$ {\bf f}({\bf x}) \f$
 /// as \f$ \mbox{\tt fct}({\bf x}, {\bf f}) \f$.
 typedef void (*SteadyExactSolutionFctPt)(const Vector<double>&,
                                                Vector<double>&);

 /// \short Function pointer for function that computes Vector-valued
 /// time-dependent function \f$ {\bf f}(t,{\bf x}) \f$
 /// as \f$ \mbox{\tt fct}(t, {\bf x}, {\bf f}) \f$.
 typedef void (*UnsteadyExactSolutionFctPt)(const double&, 
                                            const Vector<double>&, 
                                            Vector<double>&);

 /// \short Tolerance below which the jacobian is considered singular
 static double Tolerance_for_singular_jacobian;

 /// \short Boolean that if set to true allows a negative jacobian 
 /// in the transform between global and local coordinates (negative surface
 /// area = left-handed coordinate system).
 static bool Accept_negative_jacobian;

 /// \short Static boolean to suppress output while checking 
 /// for inverted elements
 static bool Suppress_output_while_checking_for_inverted_elements;

 /// Constructor
 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) {}

 /// \short The destructor cleans up the static memory allocated
 /// for shape function storage. Internal and external data get
 /// wiped by the GeneralisedElement destructor; nodes get
 /// killed in mesh destructor. 
 virtual ~FiniteElement();

 /// Broken copy constructor
 FiniteElement(const FiniteElement&) 
  { 
   BrokenCopy::broken_copy("FiniteElement");
  } 
 
 /// Broken assignment operator
//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&) 
  {
   BrokenCopy::broken_assign("FiniteElement");
   }*/

 ///Check whether the local coordinate are valid or not
 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;
  }
 
 ///\short Adjust local coordinates so that they're located inside
 /// the element
 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);
  }


 /// \short Compute centre of gravity of all nodes and radius of node that
 /// is furthest from it. Used to assess approximately if a point
 /// is likely to be contained with an element in locate_zeta-like
 /// operations
 void get_centre_of_gravity_and_max_radius_in_terms_of_zeta(
  Vector<double>& cog, double& max_radius) const;

 /// \short Get local coordinates of node j in the element; vector 
 /// sets its own size (broken virtual)
 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);
  }

 /// \short Get the local fraction of the node j in the element
 /// A dumb, but correct default implementation is provided.
 virtual void local_fraction_of_node(const unsigned &j, 
                                     Vector<double> &s_fraction);

 /// \short Get the local fraction of any node in the n-th position 
 /// in a one dimensional expansion along the i-th local coordinate
 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);
  }

 /// \short Set pointer to macro element -- can be overloaded in derived
 /// elements to perform additional tasks
 virtual void set_macro_elem_pt(MacroElement* macro_elem_pt)
  {Macro_elem_pt=macro_elem_pt;}

 ///Access function to pointer to macro element
 MacroElement* macro_elem_pt() {return Macro_elem_pt;}

 /// \short Global coordinates as function of local coordinates.
 /// Either via FE representation or via macro-element (if Macro_elem_pt!=0)
 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);}
  }


 /// \short Global coordinates as function of local coordinates
 /// at previous time "level" t (t=0: present; t>0: previous).
 /// Either via FE representation of QElement or
 /// via macro-element (if Macro_elem_pt!=0).
 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);}
  }


 /// \short Global coordinates as function of local coordinates
 /// using macro element representation. 
 /// (Broken virtual --- this  must be overloaded in specific geometric
 /// element classes)
 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);
  }
 
 /// \short Global coordinates as function of local coordinates
 /// at previous time "level" t (t=0: present; t>0: previous).
 /// using macro element representation
 /// (Broken virtual -- overload in specific geometric element class
 /// if you want to use this functionality.)
 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);
  }


 ///Set the spatial integration scheme
 virtual void set_integration_scheme(Integral* const &integral_pt); 

 /// Return the pointer to the integration scheme (const version)
 Integral* const &integral_pt() const {return Integral_pt;}
 
 /// \short Calculate the geometric shape functions 
 /// at local coordinate s. This function must be overloaded for each specific 
 /// geometric element.
 virtual void shape(const Vector<double> &s, Shape &psi) const=0;

 /// \short Return the geometric shape function at the ipt-th integration point
 virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const;

 /// \short Function to compute the geometric shape functions and
 /// derivatives w.r.t. local coordinates at local coordinate s. 
 /// This function must be overloaded for each specific geometric element.
 /// (Broken virtual function --- specifies the interface)
 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);
  }

 /// \short Return the geometric shape function and its derivative w.r.t. 
 /// the local coordinates at the ipt-th integration point.
 virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi,
                                   DShape &dpsids) const;

 /// \short Function to compute the geometric shape functions and also 
 /// first and second derivatives w.r.t. 
 /// local coordinates at local coordinate s. This function must
 /// be overloaded for each specific geometric element (if required).
 /// (Broken virtual function --- specifies the interface). 
 ///  Numbering:
 ///  \b 1D: 
 /// d2psids(i,0) = \f$ d^2 \psi_j / ds^2 \f$
 ///  \b 2D: 
 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$ 
 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$ 
 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$ 
 ///  \b 3D: 
 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$ 
 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$ 
 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$ 
 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$ 
 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$ 
 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$ 
 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);
  }

 /// \short Return the geometric shape function and its first and
 /// second derivatives w.r.t. 
 /// the local coordinates at the ipt-th integration point.
 ///  Numbering:
 ///  \b 1D: 
 /// d2psids(i,0) = \f$ d^2 \psi_j / ds^2 \f$
 ///  \b 2D: 
 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$ 
 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$ 
 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$ 
 ///  \b 3D: 
 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$ 
 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$ 
 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$ 
 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$ 
 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$ 
 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$ 
 virtual void d2shape_local_at_knot(const unsigned &ipt, Shape &psi,
                                    DShape &dpsids, DShape &d2psids) const;

 /// \short Return the Jacobian of mapping from local to global 
 /// coordinates at local position s.
 virtual double J_eulerian(const Vector<double> &s) const;

 /// \short Return the Jacobian of the mapping from local to global
 /// coordinates at the ipt-th integration point
 virtual double J_eulerian_at_knot(const unsigned &ipt) const;
 
/// \short Check that Jacobian of mapping between local and Eulerian
/// coordinates at all integration points is positive.
 void check_J_eulerian_at_knots(bool& passed) const;

 /// \short Helper function used to check for singular or negative
 /// Jacobians in the transform from local to global or Lagrangian
 /// coordinates.
 void check_jacobian(const double &jacobian) const;

 /// \short Compute the geometric shape functions and also
 /// first derivatives w.r.t. global coordinates at local coordinate s;
 /// Returns Jacobian of mapping from global to local coordinates.
 double dshape_eulerian(const Vector<double> &s, Shape &psi,
                        DShape &dpsidx) const;

 /// \short Return the geometric shape functions and also first
 /// derivatives w.r.t. global coordinates at the ipt-th integration point.
 virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi,
                                        DShape &dpsidx) const;

 /// \short Compute the geometric shape functions (psi) and first derivatives
 /// w.r.t. global coordinates (dpsidx) at the ipt-th integration point.
 /// Return the determinant of the jacobian of the mapping (detJ).
 /// Additionally calculate the derivatives of both "detJ" and "dpsidx"
 /// w.r.t. the nodal coordinates.
 virtual double dshape_eulerian_at_knot(
  const unsigned &ipt,Shape &psi,DShape &dpsi,
  DenseMatrix<double> &djacobian_dX,RankFourTensor<double> &d_dpsidx_dX) const;

 /// \short Compute the geometric shape functions and also first
 /// and second derivatives w.r.t. global coordinates at local coordinate s;
 /// Returns Jacobian of mapping from global to local coordinates.
 ///  Numbering:
 ///  \b 1D: 
 /// d2psidx(i,0) = \f$ d^2 \psi_j / d x^2 \f$
 ///  \b 2D: 
 /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$ 
 /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$ 
 /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$ 
 ///  \b 3D: 
 /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$ 
 /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$ 
 /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_2^2 \f$ 
 /// d2psidx(i,3) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$ 
 /// d2psidx(i,4) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_2 \f$ 
 /// d2psidx(i,5) = \f$ \partial^2 \psi_j / \partial x_1 \partial x_2 \f$ 
 double d2shape_eulerian(const Vector<double> &s, Shape &psi, 
                         DShape &dpsidx, DShape &d2psidx) const;


 /// \short Return the geometric shape functions and also first
 /// and second derivatives w.r.t. global coordinates at ipt-th integration
 /// point. 
 ///  Numbering:
 ///  \b 1D: 
 /// d2psidx(i,0) = \f$ d^2 \psi_j / d s^2 \f$
 ///  \b 2D: 
 /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$ 
 /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$ 
 /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$ 
 ///  \b 3D: 
 /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$ 
 /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$ 
 /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_2^2 \f$ 
 /// d2psidx(i,3) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$ 
 /// d2psidx(i,4) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_2 \f$ 
 /// d2psidx(i,5) = \f$ \partial^2 \psi_j / \partial x_1 \partial x_2 \f$ 
 virtual double d2shape_eulerian_at_knot(const unsigned &ipt, Shape &psi, 
                                         DShape &dpsidx, 
                                         DShape &d2psidx) const;
 
 /// \short Assign the local equation numbers for Data stored at the nodes
 /// Virtual so that it can be overloaded by RefineableFiniteElements.
 /// If the boolean is true then the pointers to the degrees of freedom
 /// associated with each equation number are stored in Dof_pt
 virtual void assign_nodal_local_eqn_numbers(
  const bool &store_local_dof_pt);

 /// \short Function to describe the local dofs of the element[s]. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 virtual void describe_local_dofs(std::ostream& out,
                                  const std::string& current_string) const;

 /// \short Function to describe the local dofs of the element[s]. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 virtual void describe_nodal_local_dofs(
  std::ostream& out,
  const std::string& current_string) const;

 /// \short Overloaded version of the calculation of the local equation
 /// numbers. If the boolean argument is true then pointers to the degrees
 /// of freedom associated with each equation number are stored locally
 /// in the array Dof_pt.
 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);
  }

 /// Return a pointer to the local node n
 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];
  }

 /// Return a pointer to the local node n (const version)
 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];
  }

 /// Return the number of nodes
 unsigned nnode() const {return Nnode;}
 
 /// \short Return the number of nodes along one edge of the element
 /// Default is to return zero --- must be overloaded by geometric
 /// elements
 virtual unsigned nnode_1d() const {return 0;}

 /// \short Return the i-th coordinate at local node n. Do not use
 /// the hanging node representation.
 /// NOTE: Moved to cc file because of a possible compiler bug
 /// in gcc (yes, really!). The move to the cc file avoids inlining
 /// which appears to cause problems (only) when compiled with gcc
 /// and -O3; offensive "illegal read" is in optimised-out section
 /// of code and data that is allegedly illegal is readily readable
 /// (by other means) just before this function is called so I can't
 /// really see how we could possibly be responsible for this...
 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); */
  /* } */
 
 /// \short Return the i-th coordinate at local node n, at time level t
 /// (t=0: present; t>0: previous time level).
 /// Do not use the hanging node representation.
 double raw_nodal_position(const unsigned &t, const unsigned &n, 
              const unsigned &i) const
  {return node_pt(n)->x(t,i);}
 
 /// \short Return the i-th component of nodal velocity: dx/dt at local node n.
 /// Do not use the hanging node representation.
 double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
  {return node_pt(n)->dx_dt(i);}
 
 /// \short Return the i-th component of j-th derivative of nodal position: 
 /// d^jx/dt^j at node n. Do not use the hanging node representation.
 double raw_dnodal_position_dt(const unsigned &n, const unsigned &j,
                  const unsigned &i) const
  {return node_pt(n)->dx_dt(j,i);}
 
 /// \short Return the value of the k-th type of the i-th positional variable
 /// at the local node n. Do not use the hanging node representation.
 double raw_nodal_position_gen(const unsigned &n, const unsigned &k, 
                  const unsigned &i) const
  {return node_pt(n)->x_gen(k,i);}
 
 /// \short Return the generalised nodal position (type k, i-th variable) 
 /// at previous timesteps at local node n. Do not use the hanging node
 /// representation.
 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);}
 
 /// \short i-th component of time derivative (velocity) of the 
 /// generalised position, dx(k,i)/dt at local node n.
 /// `Type': k; Coordinate direction: i. Do not use the hanging node
 /// representation.
 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);}
 
 /// \short i-th component of j-th time derivative  of the 
 /// generalised position, dx(k,i)/dt at local node n.
 /// `Type': k; Coordinate direction: i. Do not use the hanging node
 /// representation.
 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);}


 /// \short Return the i-th coordinate at local node n. If the 
 /// node is hanging, the appropriate interpolation is handled by the
 /// position function in the Node class.
 double nodal_position(const unsigned &n, const unsigned &i) const
  {return node_pt(n)->position(i);}
 
 /// \short Return the i-th coordinate at local node n, at time level t
 /// (t=0: present; t>0: previous time level)
 /// Returns suitably interpolated version for hanging nodes.
 double nodal_position(const unsigned &t, const unsigned &n, 
          const unsigned &i) const
  {return node_pt(n)->position(t,i);}

 /// Return the i-th component of nodal velocity: dx/dt at local node n.
 double dnodal_position_dt(const unsigned &n, const unsigned &i) const
  {return node_pt(n)->dposition_dt(i);}

 /// \short Return the i-th component of j-th derivative of nodal position: 
 /// d^jx/dt^j at node n.
 double dnodal_position_dt(const unsigned &n, const unsigned &j,
              const unsigned &i) const
  {return node_pt(n)->dposition_dt(j,i);}

 /// \short Return the value of the k-th type of the i-th positional variable
 /// at the local node n.
 double nodal_position_gen(const unsigned &n, const unsigned &k, 
               const unsigned &i) const
  {return node_pt(n)->position_gen(k,i);}

 /// \short Return the generalised nodal position (type k, i-th variable) 
 /// at previous timesteps at local node n
 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);}

 /// \short i-th component of time derivative (velocity) of the 
 /// generalised position, dx(k,i)/dt at local node n.
 /// `Type': k; Coordinate direction: 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);}

 /// \short i-th component of j-th time derivative  of the 
 /// generalised position, dx(k,i)/dt at local node n.
 /// `Type': k; Coordinate direction: 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);}


 /// \short Compute derivatives of elemental residual vector with respect
 /// to nodal coordinates. Default implementation by FD can be overwritten
 /// for specific elements. 
 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
 virtual void get_dresidual_dnodal_coordinates(RankThreeTensor<double>&
                                               dresidual_dnodal_coordinates);

 /// \short This is an empty function that establishes a uniform
 /// interface for all (derived) elements that involve time-derivatives.
 /// Such elements are/should be implemented in ALE form to allow
 /// mesh motions. The additional expense associated with the
 /// computation of the mesh velocities is, of course, superfluous
 /// if the elements are used in problems in which the mesh is
 /// stationary. This function should therefore be overloaded
 /// in all derived elements that are formulated in ALE form
 /// to suppress the computation of the mesh velocities.
 /// The user disables the ALE functionality at his/her own risk!
 /// If the mesh does move after all, then the results will be wrong.
 /// Here we simply issue a warning message stating that the empty
 /// function has been called.
 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);
  }


 /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
 /// when evaluating the time-derivative. This function is empty
 /// and simply establishes a common interface for all derived
 /// elements that are formulated in ALE form.
 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);
  }

 /// \short Number of values that must be stored at local node n by
 /// the element. The default is 0, until over-ridden by a particular element.
 /// For example, a Poisson equation requires only one value to be stored
 /// at each node; 2D Navier--Stokes equations require two values (velocity
 /// components) to be stored at each Node 
 /// (provided that the pressure interpolation is discontinuous).
 virtual unsigned required_nvalue(const unsigned &n) const
  {return Default_Initial_Nvalue;} 

 /// \short Return the number of coordinate types that the element requires
 /// to interpolate the geometry between
 /// the nodes. For Lagrange elements it is 1.
 unsigned nnodal_position_type() const {return Nnodal_position_type;}

 /// \short Return boolean to indicate if any of the element's nodes
 /// are geometrically hanging.
 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;
 }

 /// Return the required Eulerian dimension of the nodes in this element
 unsigned nodal_dimension() const {return Nodal_dimension;}

 /// Return the number of vertex nodes in this element. Broken virtual
 /// function in "pure" finite elements.
 virtual unsigned nvertex_node() const
  {
     std::string error_msg = "Not implemented for FiniteElement.";
     throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
  }

 /// \short Pointer to the j-th vertex node in the element. Broken virtual
 /// function in "pure" finite elements.
 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);
 }

 /// \short Construct the local node n and return a pointer to the newly 
 /// created node object.
 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);
  }

 /// \short Construct the local node n, including
 /// storage for history values required by timestepper, and return a pointer
 /// to the newly created node object.
 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);
  }
 
 /// \short Construct the local node n as a boundary node;
 /// that is a node that MAY be placed on a mesh boundary
 /// and return a pointer to the newly created node object.
 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);
  }

 /// \short Construct the local node n, including
 /// storage for history values required by timestepper, 
 /// as a boundary node; that is a node that MAY be placed on a mesh
 /// boundary and return a pointer to the newly created node object.
 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);
  }
  
 /// \short Return the number of the node *node_pt if this node
 /// is in the element, else return -1; 
 int get_node_number(Node* const &node_pt) const;


 /// \short If there is a node at this local coordinate, return the pointer to
 /// the node
 virtual Node* get_node_at_local_coordinate(const Vector<double> &s) const;

 /// \short Return the i-th value stored at local node n
 /// but do NOT take hanging nodes into account
 double raw_nodal_value(const unsigned &n, const unsigned &i) const
  {return node_pt(n)->raw_value(i);}

 /// \short Return the i-th value stored at local node n, at time level t
 /// (t=0: present; t>0 previous timesteps), but do NOT take hanging nodes
 /// into account
 double raw_nodal_value(const unsigned &t, const unsigned &n, 
                        const unsigned &i) 
  const {return node_pt(n)->raw_value(t,i);}
 
 /// \short Return the i-th value stored at local node n.
 /// Produces suitably interpolated values for hanging nodes.
 double nodal_value(const unsigned &n, const unsigned &i) const
  {return node_pt(n)->value(i);}

 /// \short Return the i-th value stored at local node n, at time level t
 /// (t=0: present; t>0 previous timesteps). Produces suitably interpolated 
 /// values for hanging nodes.
 double nodal_value(const unsigned &t, const unsigned &n, const unsigned &i) 
  const {return node_pt(n)->value(t,i);}

 /// \short Return the spatial dimension of the element, i.e.
 /// the number of local coordinates required to parametrise its 
 /// geometry.
 unsigned dim() const {return Elemental_dimension;}

 /// Return the geometry type of the element (either Q or T usually).
 virtual ElementGeometry::ElementGeometry element_geometry() const
 {
  std::string err = "Broken virtual function.";
  throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
 }

 /// Return FE interpolated coordinate x[i] at local coordinate s
 virtual double interpolated_x(const Vector<double> &s, 
                               const unsigned &i) const;

 ///\short Return FE interpolated coordinate x[i] at local coordinate s
 ///at previous timestep t (t=0: present; t>0: previous timestep)
 virtual double interpolated_x(const unsigned& t, const Vector<double> &s,
                               const unsigned &i) const;
 
 /// Return FE interpolated position x[] at local coordinate s as Vector
 virtual void interpolated_x(const Vector<double> &s,Vector<double>& x) const;

 /// \short Return FE interpolated position x[] at local coordinate s
 /// at previous timestep t as Vector (t=0: present; t>0: previous timestep)
 virtual void interpolated_x(const unsigned& t, const Vector<double> &s,
                             Vector<double>& x) const;

 /// \short Return t-th time-derivative of the
 /// i-th FE-interpolated Eulerian coordinate at
 /// local coordinate s. 
 virtual double interpolated_dxdt(const Vector<double> &s,
                                  const unsigned &i,
                                  const unsigned &t);
 
 /// \short Compte t-th time-derivative of the
 /// FE-interpolated Eulerian coordinate vector at
 /// local coordinate s. 
 virtual void interpolated_dxdt(const Vector<double> &s,
                                const unsigned &t,
                                Vector<double>& dxdt);
 
 /// \short A standard FiniteElement is fixed, so there are no geometric
 /// data when viewed in its GeomObject incarnation
 inline unsigned ngeom_data() const {return 0;}

 /// \short A standard FiniteElement is fixed, so there are no geometric
 /// data when viewed in its GeomObject incarnation
 inline Data* geom_data_pt(const unsigned &j) {return 0;}

 /// \short Return the parametrised position of the FiniteElement in 
 /// its incarnation as a GeomObject, r(zeta). 
 /// The position is given by the Eulerian coordinate and the intrinsic
 /// coordinate (zeta) is the local coordinate of the element (s).
 void position(const Vector<double>& zeta, Vector<double>& r) 
  const {this->interpolated_x(zeta,r);}

 /// \short Return the parametrised position of the FiniteElement
 /// in its GeomObject incarnation: 
 /// r(zeta). The position is given by the Eulerian coordinate and the
 /// intrinsic coordinate (zeta) is the local coordinate of the element (s)
 /// This version of the function returns the position as a function
 /// of time t=0: current time; t>0: previous
 /// timestep. Works for t=0 but needs to be overloaded
 /// if genuine time-dependence is required.
 void position(const unsigned& t, const Vector<double>& zeta, 
               Vector<double>& r) const
  {this->interpolated_x(t,zeta,r);}

 /// \short Return the t-th time derivative of the 
 /// parametrised position of the FiniteElement
 /// in its GeomObject incarnation: 
 /// \f$ \frac{d^{t} dr(zeta)}{d t^{t}} \f$. 
 /// Call the t-th time derivative of the FE-interpolated Eulerian coordinate 
 void dposition_dt(const Vector<double> &zeta, const unsigned &t,
                  Vector<double>  &drdt)
  {this->interpolated_dxdt(zeta,t,drdt);}
   
 /// \short Specify the values of the "global" intrinsic coordinate, zeta, 
 /// of a compound geometric object (a mesh of elements) when 
 /// the element is viewied as a sub-geometric object. 
 /// The default assumption is that the element will be 
 /// treated as a sub-geometric object in a bulk Mesh of other elements
 /// (geometric objects). The "global" coordinate of the compound geometric
 /// object is simply the Eulerian coordinate, x.  
 /// The second default assumption is that the coordinate zeta will
 /// be stored at the nodes and interpolated using the shape functions
 /// of the element. This function returns the value of zeta stored at
 /// local node n, where k is the type of coordinate and i is the 
 /// coordinate direction. The function is virtual so that it can
 /// be overloaded by different types of element: FaceElements
 /// and SolidFiniteElements
 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);                             
  }     


 /// \short Calculate the interpolated value of zeta, the intrinsic coordinate
 /// of the element when viewed as a compound geometric object within a Mesh
 /// as a function of the local coordinate of the element, s. The default 
 /// assumption is the zeta is interpolated using the shape functions of
 /// the element with the values given by zeta_nodal().
  /// A MacroElement representation of the intrinsic coordinate parametrised
 /// by the local coordinate s is used if available. 
 /// Choosing the MacroElement representation of zeta (Eulerian x by default)
 /// allows a correspondence to be established between elements on different
 /// Meshes covering the same curvilinear domain in cases where one element
 /// is much coarser than the other.
 void interpolated_zeta(const Vector<double> &s,                              
                        Vector<double> &zeta) const;       

 /// \short For a given value of zeta, the "global" intrinsic coordinate of
 /// a mesh of FiniteElements represented as a compound geometric object,
 /// find the local coordinate in this element that corresponds to the 
 /// requested value of zeta. 
 /// If zeta cannot be located in this element, geom_object_pt is set
 /// to NULL. If zeta is located in this element, we return its "this"
 /// pointer.
 /// By default don't use any value passed in to the local coordinate s
 /// as the initial guess in the Newton method
 void locate_zeta(const Vector<double> &zeta,
                  GeomObject*& geom_object_pt, Vector<double> &s,
                  const bool& use_coordinate_as_initial_guess=false);


 /// \short Update the positions of all nodes in the element using
 /// each node update function. The default implementation may
 /// be overloaded so that more efficient versions can be written
 virtual void node_update();

 /// \short The purpose of this function is to identify all possible
 /// Data that can affect the fields interpolated by the FiniteElement.
 /// The information will typically be used in interaction problems in
 /// which the FiniteElement provides a forcing term for an 
 /// ElementWithExternalElement. The Data must be provided as 
 /// \c paired_load data containing
 ///  - the pointer to a Data object
 /// and
 /// - the index of the value in that Data object
 /// .
 /// The generic implementation (should be overloaded in more specific
 /// applications) is to include all nodal and internal Data stored in
 /// the FiniteElement. The geometric data, which includes the positions
 /// of SolidNodes, is treated separately by the function 
 /// \c identify_geometric_data()
 virtual void identify_field_data_for_interactions(
  std::set<std::pair<Data*,unsigned> > &paired_field_data);


 /// \short The purpose of this
 /// function is to identify all \c Data objects that affect the elements'
 /// geometry. This function is implemented as an empty virtual
 /// function since it can only be implemented in conjunction
 /// with a node-update strategy. A specific implementation
 /// is provided in the ElementWithMovingNodes class.
 virtual void identify_geometric_data(std::set<Data*> &geometric_data_pt){}


 /// Min value of local coordinate
 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;
  }

 /// Max. value of local coordinate
 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;
  }

  /// Calculate the size of the element (length, area, volume,...)
  /// in Eulerian computational 
  /// coordinates. Use suitably overloaded compute_physical_size()
  /// function to compute the actual size (taking into account
  /// factors such as 2pi or radii the integrand) -- such function
  /// can only be implemented on an equation-by-equation basis.
  double size() const; 
 

  /// \short Broken virtual 
  /// function to compute the actual size (taking into account
  /// factors such as 2pi or radii the integrand) -- such function
  /// can only be implemented on an equation-by-equation basis.
 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;
  }
 
 /// \short Virtual function to write the double precision numbers that
 /// appear in a single line of output into the data vector. Empty virtual,
 /// can be overloaded for specific elements; used e.g. by LineVisualiser.
 virtual void point_output_data(const Vector<double> &s, Vector<double>& data)
 {}
 
 /// \short Output solution (as defined by point_output_data()) 
 /// at local cordinates s
 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] << " ";
   }
 }
 
 /// \short Return the number of actual plot points for paraview
 /// plot with parameter nplot. Broken virtual; can be overloaded
 /// in specific elements.
 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;
  }
 
 /// \short Return the number of local sub-elements for paraview plot with 
 /// parameter nplot. Broken virtual; can be overloaded
 /// in specific elements.
 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;
  }

 /// \short Paraview output -- this outputs the coordinates at the plot
 /// points (for parameter nplot) to specified output file.
 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);
      }
    }
  }

 /// \short Fill in the offset information for paraview plot. Broken virtual.
 /// Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 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);
  }

 /// \short Return the paraview element type.  Broken virtual.
 /// Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 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);
  }
 
 /// \short Return the offsets for the paraview sub-elements. Broken
 /// virtual. Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 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);
  }
 
 /// \short Number of scalars/fields output by this element. Broken
 /// virtual. Needs to be implemented for each new specific element type.
 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;
  }

 /// \short Write values of the i-th scalar field at the plot points. Broken
 /// virtual. Needs to be implemented for each new specific element type.
 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);
  }
 
 /// \short Write values of the i-th scalar field at the plot points. Broken
 /// virtual. Needs to be implemented for each new specific element type.
 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);
  }

 /// \short Name of the i-th scalar field. Default implementation
 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
 /// overloaded with more meaningful names in specific elements.
 virtual std::string scalar_name_paraview(const unsigned& i) const
  {
   return "V"+StringConversion::to_string(i);
  }

 /// \short Output the element data --- typically the values at the 
 /// nodes in a format suitable for post-processing.
 virtual void output(std::ostream &outfile)
  {
    throw OomphLibError(
    "Output function function hasn't been implemented for this element",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// \short Output the element data --- pass (some measure of) 
 /// the number of plot points per element
 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);
  }

 /// \short Output the element data at time step t. This is const because
 /// it is newly added and so can be done easily. Really all the
 /// output(...) functions should be const!
 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);
 }

 /// \short Output the element data --- typically the values at the 
 /// nodes in a format suitable for post-processing.
 /// (C style output)
 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);
  }

 /// \short Output the element data --- pass (some measure of) 
 /// the number of plot points per element
 /// (C style output)
 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);
  }

 /// Output an exact solution over the element.
 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);
  }


 /// Output a time-dependent exact solution over the element.
 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);
  }

 /// Output a time-dependent exact solution over the element.
 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);
 }



 /// \short  Get cector of local coordinates of plot point i (when plotting 
 /// nplot points in each "coordinate direction").
 /// Generally these plot points will be uniformly spaced across the element.
 /// The optional final boolean flag (default: false) allows them to
 /// be shifted inwards to avoid duplication of plot point points between
 /// elements -- useful when they are used in locate_zeta, say.
 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);
  };
 
 /// \short Return string for tecplot zone header (when plotting 
 /// nplot points in each "coordinate direction")
 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";
  }

 /// \short Add tecplot zone "footer" to output stream (when plotting 
 /// nplot points in each "coordinate direction").
 /// Empty by default -- can be used, e.g., to add FE connectivity
 /// lists to elements that need it.
 virtual void write_tecplot_zone_footer(std::ostream& outfile, 
                                        const unsigned& nplot)
 const {};
 
 /// \short Add tecplot zone "footer" to C-style output. (when plotting 
 /// nplot points in each "coordinate direction").
 /// Empty by default -- can be used, e.g., to add FE connectivity
 /// lists to elements that need it.
 virtual void write_tecplot_zone_footer(FILE* file_pt,
                                        const unsigned& nplot)
 const {};
 
 /// \short Return total number of plot points (when plotting 
 /// nplot points in each "coordinate direction")
 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;
  }



 /// \short Plot the error when compared 
 /// against a given exact solution \f$ {\bf f}({\bf x}) \f$.
 /// Also calculates the norm of the
 /// error and that of the exact solution. 
 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);

  }

 /// \short Plot the error when compared 
 /// against a given time-dependent exact solution \f$ {\bf f}(t,{\bf x}) \f$.
 /// Also calculates the norm of the error and that of the exact solution.
 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);
  }

 /// \short Plot the error when compared 
 /// against a given exact solution \f$ {\bf f}({\bf x}) \f$.
 /// Also calculates the norm of the
 /// error and that of the exact solution. 
 /// Version with vectors of norms and errors so that different variables' norms
 /// and errors can be returned individually
 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);

  }

 /// \short Plot the error when compared 
 /// against a given time-dependent exact solution \f$ {\bf f}(t,{\bf x}) \f$.
 /// Also calculates the norm of the error and that of the exact solution.
 /// Version with vectors of norms and errors so that different variables' norms
 /// and errors can be returned individually
 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);

  }

 /// \short Plot the error when compared 
 /// against a given exact solution \f$ {\bf f}({\bf x}) \f$.
 /// Also calculates the maximum absolute error
 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);
  }

 /// \short Evaluate integral of a Vector-valued function
 /// \f$ {\bf f}({\bf x}) \f$ over the element.
 void integrate_fct(FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt,
                    Vector<double>& integral);


 /// \short Evaluate integral of a Vector-valued, time-dependent function 
 /// \f$ {\bf f}(t,{\bf x}) \f$ over the element.
 void integrate_fct(FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt,
                    const double& time,
                    Vector<double>& integral);

 /// \short Function for building a lower dimensional FaceElement
 /// on the specified face of the FiniteElement.
 /// The arguments are the index of the face, an integer whose value
 /// depends on the particular element type, and a pointer to the
 /// FaceElement.
 virtual void build_face_element(const int &face_index, 
                                 FaceElement* face_element_pt);


 /// \short Self-test: Check inversion of element & do self-test for
 /// GeneralisedElement. Return 0 if OK. 
 virtual unsigned self_test();

 /// Get the number of the ith node on face face_index (in the bulk node
 /// vector). 
 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);
  }

 /// Get the sign of the outer unit normal on the face given by face_index.
 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);
  }

 /// Range check for face node numbers
 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
 }

 /// Get a pointer to the function mapping face coordinates to bulk coordinates
 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);
  }

 /// Get a pointer to the derivative of the mapping from face to bulk
 /// coordinates.
 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);
  }

};


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



//=======================================================================
/// Point element has just a single node and a single shape function
/// which is identically equal to one.
//=======================================================================
class PointElement : public virtual FiniteElement
{

  public:
 
 /// Constructor
 PointElement()
  {
   /// Allocate storage for the pointers to the nodes,
   /// There is only one nodes
   this->set_n_node(1);

   // Set the integration scheme
   this->set_integration_scheme(&Default_integration_scheme);
  }

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

 /// Calculate the geometric shape functions at local coordinate s
 void shape(const Vector<double> &s, Shape &psi) const;

 /// Get local coordinates of node j in the element; vector sets its own size
 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
  {
   s.resize(0);
  }

 /// \short  Return spatial dimension of element (=number of local coordinates
 /// needed to parametrise the element)
 //unsigned dim() const {return 0;}

  private:
 
 /// \short Default integration scheme
 static PointIntegral Default_integration_scheme;


};



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


class GeomObject;

//=======================================================================
/// \short A class to specify the initial conditions for a solid body.
/// Solid bodies are often discretised with
/// Hermite-type elements, for which the assignment
/// of the generalised nodal values is nontrivial since they represent
/// derivatives  w.r.t. to the local coordinates. A SolidInitialCondition
/// object specifies initial position (i.e. shape), velocity and acceleration
/// of the structure with a geometric object.
/// An integer specifies which time-derivative derivative is currently
/// assigned. See example codes for a demonstration of its use.
//=======================================================================
class SolidInitialCondition
{

  public:

 /// Constructor: Pass geometric object; initialise time deriv to 0
 SolidInitialCondition(GeomObject* geom_object_pt) :
  Geom_object_pt(geom_object_pt), IC_time_deriv(0)
  {};


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


 /// (Reference to) pointer to geom object that specifies the initial condition
 GeomObject*& geom_object_pt()
  {
   return Geom_object_pt;
  }

 /// Which time derivative are we currently assigning?
 unsigned& ic_time_deriv()
  {
   return IC_time_deriv;
  }


  private:

 /// \short Pointer to the GeomObject that specifies the initial
 /// condition (shape, veloc and accel)
 GeomObject* Geom_object_pt;

 /// \short Which time derivative (0,1,2) are we currently assigning
 unsigned IC_time_deriv;

};


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


//============================================================================
/// \short SolidFiniteElement class. 
///  
/// Solid elements are elements whose nodal positions are unknowns
/// in the problem -- their nodes are SolidNodes. In such elements,
/// the nodes not only have a variable (Eulerian) but also a fixed
/// (Lagrangian) position. The positional variables have their own
/// local equation numbering scheme which is set up with
/// \code assign_solid_local_eqn_numbers () \endcode
/// The derivatives of the
/// `solid equations' (i.e. the equations that determine the
/// nodal positions) with respect to the nodal positions, required 
/// in the Jacobian matrix, are determined by finite differencing.
///  
/// In the present form, the SolidFiniteElement represents 
/// a fully functional base class for `proper' solid mechanics elements,
/// but it can also be combined 
/// (via inheritance) with elements that solve additional equations.
/// This is particularly useful in cases where the solid equations
/// are merely used to update the nodal positions in a moving mesh
/// problem, although this can prove costly in practice. 
//============================================================================
class SolidFiniteElement : public virtual FiniteElement
{

public:
 
 /// \short Set the lagrangian dimension of the element --- the number
 /// of lagrangian coordinates stored at the nodes in the element.
 void set_lagrangian_dimension(const unsigned &lagrangian_dimension)
  {Lagrangian_dimension = lagrangian_dimension;}

 /// \short Return whether there is internal solid data (e.g. discontinuous
 /// solid pressure). At present, this is used to report an error in
 /// setting initial conditions for ElasticProblems which can't
 /// handle such cases. The default is false.
 virtual bool has_internal_solid_data() {return false;}

 /// \short Pointer to function that computes the "multiplier" for the
 /// inertia terms in the consistent determination of the initial
 /// conditions for Newmark timestepping. 
 typedef double (*MultiplierFctPt)(const Vector<double>& xi);

 ///Constructor: Set defaults
 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)
  {}
  
 ///Destructor to clean up any allocated memory
 virtual ~SolidFiniteElement();

 /// Broken copy constructor
 SolidFiniteElement(const SolidFiniteElement&) 
  { 
   BrokenCopy::broken_copy("SolidFiniteElement");
  } 
 

 /// Broken assignment operator
 /*void operator=(const SolidFiniteElement&) 
  {
   BrokenCopy::broken_assign("SolidFiniteElement");
   }*/

 ///\short The number of geometric data affecting a SolidFiniteElemnet is
 ///the same as the number of nodes (one variable position data per node)
 inline unsigned ngeom_data() const {return nnode();}

 /// \short Return pointer to the j-th Data item that the object's 
 /// shape depends on. (Redirects to the nodes' positional Data).
 inline Data* geom_data_pt(const unsigned& j)
  {return static_cast<SolidNode*>(node_pt(j))->variable_position_pt();}

 /// \short Specify Data that affects the geometry of the element 
 /// by adding the position Data to the set that's passed in.
 /// (This functionality is required in FSI problems; set is used to
 /// avoid double counting). 
 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());
    }
  }


 /// \short In a SolidFiniteElement, the "global" intrinsic coordinate
 /// of the element when viewed as part of a compound geometric
 /// object (a Mesh) is, by default, the Lagrangian coordinate
 /// Note the assumption here is that we are always using isoparameteric
 /// elements in which the lagrangian coordinate is interpolated by the
 /// same shape functions as the eulerian coordinate.
 inline double zeta_nodal(const unsigned &n, const unsigned &k,           
                          const unsigned &i) const 
  {return lagrangian_position_gen(n,k,i);}     

 /// \short Eulerian and Lagrangian coordinates as function of the
 /// local coordinates: The Eulerian position is returned in 
 /// FE-interpolated form (\c x_fe) and then in the form obtained
 /// from the "current" MacroElement representation (if it exists -- if not,
 /// \c x is the same as \c x_fe). This allows the Domain/MacroElement-
 /// based representation to be used to apply displacement boundary
 /// conditions exactly. Ditto for the Lagrangian coordinates returned
 /// in xi_fe and xi.
 /// (Broken virtual -- overload in specific geometric element class
 /// if you want to use this functionality.)
 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);
  }

 /// \short Set pointer to MacroElement -- overloads generic version
 /// and uses the MacroElement 
 /// also as the default for the "undeformed" configuration.
 /// This assignment must be overwritten with 
 /// set_undeformed_macro_elem_pt(...) if the deformation of 
 /// the solid body is driven by a deformation of the
 /// "current" Domain/MacroElement representation of it's boundary.
 /// Can be overloaded in derived classes to perform additional 
 /// tasks
 virtual void set_macro_elem_pt(MacroElement* macro_elem_pt)
  {
   Macro_elem_pt=macro_elem_pt;
   Undeformed_macro_elem_pt=macro_elem_pt;
  }
 
 /// \short Set pointers to "current" and "undeformed" MacroElements.
 /// Can be overloaded in derived classes to perform additional 
 /// tasks
 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;
  }
 
 /// \short Set pointer to "undeformed" macro element.
 /// Can be overloaded in derived classes to perform additional 
 /// tasks.
 void set_undeformed_macro_elem_pt(MacroElement* 
                                   undeformed_macro_elem_pt)
  {Undeformed_macro_elem_pt=undeformed_macro_elem_pt;}


 /// Access function to pointer to "undeformed" macro element
 MacroElement* undeformed_macro_elem_pt()
  {return Undeformed_macro_elem_pt;}

 /// \short Calculate shape functions and derivatives w.r.t. Lagrangian 
 /// coordinates at local coordinate s. Returns the Jacobian of the mapping
 /// from Lagrangian to local coordinates.
 double dshape_lagrangian(const Vector<double> &s, Shape &psi,
                          DShape &dpsidxi) const;
 
 /// \short Return the geometric shape functions and also first
 /// derivatives w.r.t. Lagrangian coordinates at ipt-th integration point.
 virtual double dshape_lagrangian_at_knot(const unsigned &ipt,
                                          Shape &psi, DShape &dpsidxi) const;
  
 /// \short Compute the geometric shape functions and also first
 /// and second derivatives w.r.t. Lagrangian coordinates at 
 /// local coordinate s;
 /// Returns Jacobian of mapping from Lagrangian to local coordinates.
 ///  Numbering:
 ///  \b 1D: 
 /// d2pidxi(i,0) = \f$ d^2 \psi_j / d \xi^2 \f$
 ///  \b 2D: 
 /// d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial \xi_0^2 \f$ 
 /// d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial \xi_1^2 \f$ 
 /// d2psidxi(i,2) = \f$ \partial^2 \psi_j / \partial \xi_0 \partial \xi_1 \f$ 
 ///  \b 3D: 
 /// d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial \xi_0^2 \f$ 
 /// d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial \xi_1^2 \f$ 
 /// d2psidxi(i,2) = \f$ \partial^2 \psi_j / \partial \xi_2^2 \f$ 
 /// d2psidxi(i,3) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$ 
 /// d2psidxi(i,4) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_2 \f$ 
 /// d2psidxi(i,5) = \f$ \partial^2 \psi_j/\partial \xi_1 \partial \xi_2 \f$ 
 double d2shape_lagrangian(const Vector<double> &s, Shape &psi, 
                           DShape &dpsidxi, DShape &d2psidxi) const;

 /// \short Return  the geometric shape functions and also first
 /// and second derivatives w.r.t. Lagrangian coordinates at 
 /// the ipt-th integration point. 
 /// Returns Jacobian of mapping from Lagrangian to local coordinates.
 ///  Numbering:
 ///  \b 1D: 
 /// d2pidxi(i,0) = \f$ d^2 \psi_j / d^2 \xi^2 \f$
 ///  \b 2D: 
 /// d2psidxi(i,0) = \f$ \partial^2 \psi_j/\partial^2 \xi_0^2 \f$ 
 /// d2psidxi(i,1) = \f$ \partial^2 \psi_j/\partial^2 \xi_1^2 \f$ 
 /// d2psidxi(i,2) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$ 
 ///  \b 3D: 
 /// d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial^2 \xi_0^2 \f$ 
 /// d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial^2 \xi_1^2 \f$ 
 /// d2psidxi(i,2) = \f$ \partial^2 \psi_j / \partial^2 \xi_2^2 \f$ 
 /// d2psidxi(i,3) = \f$ \partial^2 \psi_j / \partial \xi_0 \partial \xi_1 \f$ 
 /// d2psidxi(i,4) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_2 \f$ 
 /// d2psidxi(i,5) = \f$ \partial^2 \psi_j/\partial \xi_1 \partial \xi_2 \f$ 
 virtual double d2shape_lagrangian_at_knot(const unsigned &ipt, 
                                           Shape &psi, 
                                           DShape &dpsidxi, 
                                           DShape &d2psidxi) const;

 /// \short Return the number of Lagrangian coordinates that the
 /// element requires at all nodes.
 /// This is by default the elemental dimension. If we ever need any
 /// other case, it can be implemented.
 unsigned lagrangian_dimension() const {return Lagrangian_dimension;}

 /// \short Return the number of types of (generalised)
 /// nodal Lagrangian coordinates required to
 /// interpolate the Lagrangian coordinates in the element. (E.g.
 /// 1 for Lagrange-type elements; 2 for Hermite beam elements;
 /// 4 for Hermite shell elements). Default value is 1. Needs to
 /// be overloaded for any other element.
 unsigned nnodal_lagrangian_type() const {return Nnodal_lagrangian_type;}

 /// Construct the local node n and return a pointer to it.
 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);
  }

 ///\short  Construct the local node n and return
 /// a pointer to it. Additionally, create storage for `history'
 /// values as required by timestepper
 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);
  }

 /// \short Construct the local node n and return a pointer to it.
 /// in the case when it is a boundary node; that is it MAY be
 /// located on a Mesh boundary
 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);
  }

 ///\short  Construct the local node n and return
 /// a pointer to it, in the case when the node MAY be located
 /// on a boundary. Additionally, create storage for `history'
 /// values as required by timestepper
 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);
  }
 
 /// \short Overload assign_all_generic_local_equation numbers to 
 /// include the data associated with solid dofs. 
 /// It remains virtual so that it can be overloaded
 /// by RefineableSolidElements. If the boolean argument is true
 /// then the degrees of freedom are stored in Dof_pt
 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);
  }

 /// \short Function to describe the local dofs of the element. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 void describe_local_dofs(std::ostream& out,
                          const std::string& current_string) const;
 
 /// \short Return i-th Lagrangian coordinate at local node n without using
 /// the hanging representation
 double raw_lagrangian_position(const unsigned &n, const unsigned &i) const
  {return static_cast<SolidNode*>(node_pt(n))->xi(i);}

 /// \short Return Generalised Lagrangian coordinate at local node n.
 /// `Direction' i, `Type' k. Does not use the hanging node representation 
 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);}

 /// Return i-th Lagrangian coordinate at local node n 
 double lagrangian_position(const unsigned &n, const unsigned &i) const
  {return static_cast<SolidNode*>(node_pt(n))->lagrangian_position(i);}

 /// \short Return Generalised Lagrangian coordinate at local node n.
 /// `Direction' i, `Type' k. 
 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);}
 
 /// \short Return i-th FE-interpolated Lagrangian coordinate xi[i] at
 /// local coordinate s.
 virtual double interpolated_xi(const Vector<double> &s, 
                                const unsigned &i) const;

 /// \short Compute FE interpolated Lagrangian coordinate vector xi[] at 
 /// local coordinate s as Vector
 virtual void interpolated_xi(const Vector<double> &s, 
                              Vector<double>& xi) const;

 /// \short Compute derivatives of FE-interpolated Lagrangian coordinates xi
 /// with respect to local coordinates: dxids[i][j]=dxi_i/ds_j.
 virtual void interpolated_dxids(const Vector<double> &s, 
                                 DenseMatrix<double> &dxids) const;
  
 /// \short Return the Jacobian of mapping from local to Lagrangian
 /// coordinates at local position s.  NOT YET IMPLEMENTED
 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);
  }
 
 /// \short Return the Jacobian of the mapping from local to Lagrangian
 /// coordinates at the ipt-th integration point. NOT YET IMPLEMENTED
 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);
  }
 
 /// \short Pointer to object that describes the initial condition.
 SolidInitialCondition*& solid_ic_pt()
  {
  return Solid_ic_pt;
  }

 /// \short Set to alter the problem being solved when 
 /// assigning the initial conditions for time-dependent problems:
 /// solve for the history value that
 /// corresponds to the acceleration in the Newmark scheme by demanding
 /// that the PDE is satisifed at the initial time.  In this case
 /// the Jacobian is replaced by the mass matrix.
 void enable_solve_for_consistent_newmark_accel()
 {Solve_for_consistent_newmark_accel_flag=true;}

 /// \short Set to reset the problem being solved to be the standard problem
 void disable_solve_for_consistent_newmark_accel()
 {Solve_for_consistent_newmark_accel_flag=false;}

 /// \short Access function: Pointer to multiplicator function for assignment
 /// of consistent assignement of initial conditions for Newmark scheme
 MultiplierFctPt& multiplier_fct_pt() 
  {
   return Multiplier_fct_pt;
  }


 /// \short Access function: Pointer to multiplicator function for assignment
 /// of consistent assignement of initial conditions for Newmark scheme
 /// (const version)
 MultiplierFctPt multiplier_fct_pt() const  
  {
   return Multiplier_fct_pt;
  }


 /// \short Compute the residuals for the setup of an initial condition.
 /// The global equations are:
 /// \f[
 /// 0 = \int \left( \sum_{j=1}^N \sum_{k=1}^K X_{ijk} \psi_{jk}(\xi_n) 
 /// - \frac{\partial^D R^{(IC)}_i(\xi_n)}{\partial t^D}
 /// \right) \psi_{lm}(\xi_n) \ dv
 /// \mbox{ \ \ \ \ for \ \ \ $l=1,...,N, \ \ m=1,...,K$}
 /// \f]
 /// where \f$ N \f$ is the number of nodes in the mesh and \f$ K \f$
 /// the number of generalised nodal coordinates. The initial shape
 /// of the solid body, \f$ {\bf R}^{(IC)},\f$ and its time-derivatives
 /// are specified via the \c GeomObject that is stored in the 
 /// \c SolidFiniteElement::SolidInitialCondition object. The latter also
 /// stores the order of the time-derivative \f$ D \f$ to be assigned.
 virtual void get_residuals_for_solid_ic(Vector<double>& residuals)
  {
   residuals.initialise(0.0);
   fill_in_residuals_for_solid_ic(residuals);
  }

 /// \short Fill in the residuals for the setup of an initial condition.
 /// The global equations are:
 /// \f[
 /// 0 = \int \left( \sum_{j=1}^N \sum_{k=1}^K X_{ijk} \psi_{jk}(\xi_n) 
 /// - \frac{\partial^D R^{(IC)}_i(\xi_n)}{\partial t^D}
 /// \right) \psi_{lm}(\xi_n) \ dv
 /// \mbox{ \ \ \ \ for \ \ \ $l=1,...,N, \ \ m=1,...,K$}
 /// \f]
 /// where \f$ N \f$ is the number of nodes in the mesh and \f$ K \f$
 /// the number of generalised nodal coordinates. The initial shape
 /// of the solid body, \f$ {\bf R}^{(IC)},\f$ and its time-derivatives
 /// are specified via the \c GeomObject that is stored in the 
 /// \c SolidFiniteElement::SolidInitialCondition object. The latter also
 /// stores the order of the time-derivative \f$ D \f$ to be assigned.
 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);
  }

 /// \short Fill in the residuals and Jacobian for the setup of an 
 /// initial condition. The global equations are:
 /// \f[
 /// 0 = \int \left( \sum_{j=1}^N \sum_{k=1}^K X_{ijk} \psi_{jk}(\xi_n) 
 /// - \frac{\partial^D R^{(IC)}_i(\xi_n)}{\partial t^D}
 /// \right) \psi_{lm}(\xi_n) \ dv
 /// \mbox{ \ \ \ \ for \ \ \ $l=1,...,N, \ \ m=1,...,K$}
 /// \f]
 /// where \f$ N \f$ is the number of nodes in the mesh and \f$ K \f$
 /// the number of generalised nodal coordinates. The initial shape
 /// of the solid body, \f$ {\bf R}^{(IC)},\f$ and its time-derivatives
 /// are specified via the \c GeomObject that is stored in the 
 /// \c SolidFiniteElement::SolidInitialCondition object. The latter also
 /// stores the order of the time-derivative \f$ D \f$ to be assigned.
 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);
  }


 /// \short Fill in the contributions of the Jacobian matrix 
 /// for the consistent assignment of the initial "accelerations" in 
 /// Newmark scheme. In this case the Jacobian is the mass matrix.
  void fill_in_jacobian_for_newmark_accel(DenseMatrix<double> &jacobian);
 


 /// \short Calculate the L2 norm of the displacement u=R-r to overload the
 /// compute_norm function in the GeneralisedElement base class
 void compute_norm(double& el_norm);

  protected:

 /// \short Helper function to fill in the residuals and (if flag==1) the 
 /// Jacobian for the setup of an initial condition. The global equations are:
 /// \f[
 /// 0 = \int \left( \sum_{j=1}^N \sum_{k=1}^K X_{ijk} \psi_{jk}(\xi_n) 
 /// - \frac{\partial^D R^{(IC)}_i(\xi_n)}{\partial t^D}
 /// \right) \psi_{lm}(\xi_n) \ dv
 /// \mbox{ \ \ \ \ for \ \ \ $l=1,...,N, \ \ m=1,...,K$}
 /// \f]
 /// where \f$ N \f$ is the number of nodes in the mesh and \f$ K \f$
 /// the number of generalised nodal coordinates. The initial shape
 /// of the solid body, \f$ {\bf R}^{(IC)},\f$ and its time-derivatives
 /// are specified via the \c GeomObject that is stored in the 
 /// \c SolidFiniteElement::SolidInitialCondition object. The latter also
 /// stores the order of the time-derivative \f$ D \f$ to be assigned.
 void fill_in_generic_jacobian_for_solid_ic(Vector<double> &residuals,
                                            DenseMatrix<double> &jacobian,
                                            const unsigned& flag); 

 /// \short Set the number of types required to interpolate the 
 /// Lagrangian coordinates
 void set_nnodal_lagrangian_type(const unsigned &nlagrangian_type)
  {Nnodal_lagrangian_type = nlagrangian_type;}
 
 /// Pointer to the element's "undeformed" macro element (NULL by default)
 MacroElement* Undeformed_macro_elem_pt;

 /// \short Calculate the mapping from local to lagrangian coordinates,
 /// given the derivatives of the shape functions w.r.t. local coorindates.
 /// Return the determinant of the jacobian, the jacobian  and inverse jacobian
 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);
  }
 
 /// \short Calculate the mapping from local to lagrangian coordinates,
 /// given the derivatives of the shape functions w.r.t. local coordinates,
 /// Return only the determinant of the jacobian and the inverse of the
 /// mapping (ds/dx)
 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);
  }
 
 /// \short Calculate the mapping from local to Lagrangian coordinates given 
 /// the derivatives of the shape functions w.r.t the local coorindates.
 /// assuming that the coordinates are aligned in the direction of the local 
 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
 /// This function returns the determinant of the jacobian, the jacobian 
 /// and the inverse jacobian.
 virtual double local_to_lagrangian_mapping_diagonal(
  const DShape &dpsids,DenseMatrix<double> &jacobian,
  DenseMatrix<double> &inverse_jacobian) const;
 

 /// \short Assigns local equation numbers for the generic solid 
 /// local equation numbering schemes.
 /// If the boolean flag is true the the degrees of freedom are stored in
 /// Dof_pt
 virtual void assign_solid_local_eqn_numbers(const bool &store_local_dof);

 /// \short Classifies dofs locally for solid specific aspects.
 void describe_solid_local_dofs(std::ostream& out,
                                const std::string& current_string) const;

 /// \short Pointer to object that specifies the initial condition
 SolidInitialCondition* Solid_ic_pt;

  public:

 /// \short Access function that returns the local equation number that
 /// corresponds to the j-th coordinate of the k-th position-type at the 
 /// n-th local node.
 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:


 /// \short Overload the fill_in_contribution_to_jacobian() function to
 /// use finite
 /// differences to calculate the solid residuals by default
 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);
  }



 /// \short Use finite differences to calculate the Jacobian entries
 /// corresponding to the solid positions. This version assumes
 /// that the residuals vector has already been computed
virtual void 
 fill_in_jacobian_from_solid_position_by_fd(Vector<double> &residuals,
                                         DenseMatrix<double> &jacobian);


 /// \short Use finite differences to calculate the Jacobian entries 
 /// corresponding to the solid positions.
 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);
  }

 /// \short Function that is called before the finite differencing of
 /// any solid position data. This may be overloaded to update any slaved
 /// data before finite differencing takes place.
 virtual inline void update_before_solid_position_fd() { }

 /// \short Function that is call after the finite differencing of
 /// the solid position data. This may be overloaded to reset any slaved
 /// variables that may have changed during the finite differencing.
 virtual inline void reset_after_solid_position_fd() { }

 /// \short Function called within the finite difference loop for 
 /// the solid position dat after a change in any values in the n-th 
 ///  node.
 virtual inline void update_in_solid_position_fd(const unsigned &i) { }

 /// \short Function called within the finite difference loop for
 /// solid position data after the values in the i-th node
 /// are reset. The default behaviour is to call the update function.
 virtual inline void reset_in_solid_position_fd(const unsigned &i)
  {update_in_solid_position_fd(i);}


private:
 
 /// \short Assemble the jacobian matrix for the mapping from local
 /// to lagrangian coordinates, given the derivatives of the shape function
 virtual void assemble_local_to_lagrangian_jacobian(
  const DShape &dpsids, DenseMatrix<double> &jacobian) const;


 /// \short Assemble the the "jacobian" matrix of second derivatives, given
 /// the second derivatives of the shape functions w.r.t. local coordinates
 virtual void assemble_local_to_lagrangian_jacobian2(
  const DShape &d2psids, DenseMatrix<double> &jacobian2) const;
 
 /// \short Pointer to function that computes the "multiplier" for the
 /// inertia terms in the consistent determination of the initial
 /// conditions for Newmark timestepping. 
 MultiplierFctPt Multiplier_fct_pt;
 
 /// \short Array to hold the local equation number information for the 
 /// solid equations, whatever they may be.
 // ALH: (This is here so that the numbering can be done automatically)
int *Position_local_eqn;

/// \short The Lagrangian dimension of the nodes stored in the element,
//// i.e. the number of Lagrangian coordinates
unsigned Lagrangian_dimension;

 /// \short The number of coordinate types requried to intepolate the
 /// Lagrangian coordinates in the element.  For Lagrange elements 
 /// it is 1 (the default). It must be over-ridden by using
 /// the set_nlagrangian_type() function in the constructors of elements
 /// that use generalised coordinate, e.g. for 1D Hermite elements
 /// Nnodal_position_types =2.
unsigned Nnodal_lagrangian_type;

protected:

/// \short Flag to indicate which system of equations to solve
/// when assigning initial conditions for time-dependent problems.
/// If true, solve for the history value that
/// corresponds to the acceleration in the Newmark scheme by demanding
/// that the PDE is satisifed at the initial time. In this case
/// the Jacobian is replaced by the mass matrix.
bool Solve_for_consistent_newmark_accel_flag;

private:

 /// \short Access to the "multiplier" for the
 /// inertia terms in the consistent determination of the initial
 /// conditions for Newmark timestepping. 
 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);
    }
  }


};




//========================================================================
/// FaceElements are elements that coincide with the faces of
/// higher-dimensional "bulk" elements. They are used on 
/// boundaries where additional non-trivial boundary conditions need to be
/// applied. Examples include free surfaces, and applied traction conditions.
/// In many cases, FaceElements need to evaluate to quantities
/// in the associated bulk elements. For instance, the evaluation
/// of a shear stresses on 2D FaceElement requires the evaluation
/// of velocity derivatives in the associated 3D volume element etc.
/// Therefore we store a pointer to the associated bulk element, 
/// and information about the relation between the local
/// coordinates in the face and bulk elements. 
//========================================================================
class FaceElement: public virtual FiniteElement
{
 /// \short Typedef for the function that translates the face coordinate
 /// to the coordinate in the bulk element
 typedef void (*CoordinateMappingFctPt)(const Vector<double> &s,
                                       Vector<double> &s_bulk);

 /// \short Typedef for the function that returns the partial derivative
 /// of the local coordinates in the bulk element
 /// with respect to the coordinates along the face.
 /// In addition this function returns an index of one of the 
 /// bulk local coordinates that varies away from the edge
 typedef void (*BulkCoordinateDerivativesFctPt)(
  const Vector<double> &s, DenseMatrix<double> &ds_bulk_dsface,
  unsigned &interior_direction);

  private:

 /// \short Pointer to a function that translates the face coordinate
 /// to the coordinate in the bulk element
 CoordinateMappingFctPt Face_to_bulk_coordinate_fct_pt;

 /// \short Pointer to a function that returns the derivatives of the local
 /// "bulk" coordinates with respect to the local face coordinates.
 BulkCoordinateDerivativesFctPt Bulk_coordinate_derivatives_fct_pt;

 /// \short Vector holding integers to translate additional position types
 /// to those of the "bulk" element; e.g. Bulk_position_type(1) is
 /// the position type of the "bulk" element that corresponds to face position
 /// type 1. This is required in QHermiteElements, where the slope in the 
 /// direction of the 1D element is face position type 1, but may be 
 /// position type 1 or 2 in the "bulk" element, depending upon which face,
 /// we are located.
 Vector<unsigned> Bulk_position_type;

 /// \short Sign of outer unit normal (relative to cross-products of tangent
 /// vectors in the corresponding "bulk" element.
 int Normal_sign;

 /// \short Index of the face
 int Face_index;
 
 /// \short Ignores the warning when the tangent vectors from
 /// continuous_tangent_and_outer_unit_normal may not be continuous
 /// as a result of the unit normal being aligned with the z axis.
 /// This can be avoided by supplying a general direction vector for
 /// the tangent vector via set_tangent_direction(...).
 static bool Ignore_discontinuous_tangent_warning;

  protected:
 
 /// The boundary number in the bulk mesh to which this element is attached
 unsigned Boundary_number_in_bulk_mesh;
 
#ifdef PARANOID
 
 /// \short Has the Boundary_number_in_bulk_mesh been set? Only included if
 /// compiled with PARANOID switched on.
 bool Boundary_number_in_bulk_mesh_has_been_set;

#endif

 /// \short Pointer to the associated higher-dimensional "bulk" element
 FiniteElement* Bulk_element_pt;
 
 /// \short List of indices of the local node numbers in the "bulk" element
 /// that correspond to the local node numbers in the FaceElement
 Vector<unsigned> Bulk_node_number;

 /// \short A vector that will hold the number of data values at the nodes that
 /// are associated with the "bulk" element. i.e. not including any additional
 /// degrees of freedom that might be required for extra equations
 /// that are being solved by
 /// the FaceElement. 
 // 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;

 /// \short A general direction pointer for the tangent vectors.
 /// This is used in the function continuous_tangent_and_outer_unit_normal()
 /// for creating continuous tangent vectors in spatial dimensions. 
 /// The general direction is projected on to the surface. This technique is
 /// not required in two spatial dimensions.
 Vector<double>* Tangent_direction_pt;

 /// \short Helper function adding additional values for the unknowns
 /// associated with the FaceElement. This function also sets the map 
 /// containing the position of the first entry of this face element's
 /// additional values.The inputs are the number of additional values 
 /// and the face element's ID. Note the same number of additonal values
 /// are allocated at ALL nodes.
 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:
 
 /// Constructor: Initialise all appropriate member data
 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);
  }

 /// Empty virtual destructor 
 virtual ~FaceElement() {}


 /// Broken copy constructor
 FaceElement(const FaceElement&) 
  { 
   BrokenCopy::broken_copy("FaceElement");
  } 
 
 /// Broken assignment operator
 /*void operator=(const FaceElement&) 
  {
   BrokenCopy::broken_assign("FaceElement");
   }*/
 
 /// Access function for the boundary number in bulk mesh
 inline const unsigned& boundary_number_in_bulk_mesh() const
  {return Boundary_number_in_bulk_mesh;}
 
 
 /// Set function for the 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
  }

 /// \short In a FaceElement, the "global" intrinsic coordinate
 /// of the element along the boundary, when viewed as part of
 /// a compound geometric object is specified using the
 /// boundary coordinate defined by the mesh. 
 /// Note: Boundary coordinates will have been set up when
 /// creating the underlying mesh, and their values will have 
 /// been stored at the nodes.
 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];
  }

 /// \short Return the Jacobian of mapping from local to global 
 /// coordinates at local position s.
 /// Overloaded from FiniteElement.
 double J_eulerian(const Vector<double> &s) const;

 /// \short Return the Jacobian of the mapping from local to global
 /// coordinates at the ipt-th integration point
 /// Overloaded from FiniteElement.
 double J_eulerian_at_knot(const unsigned &ipt) const;
 
 /// \short Check that Jacobian of mapping between local and Eulerian
/// coordinates at all integration points is positive.
 void check_J_eulerian_at_knots(bool& passed) const;
 
 /// \short Return FE interpolated coordinate x[i] at local coordinate s.
 /// Overloaded to get information from bulk.
 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);
  }
 
 ///\short Return FE interpolated coordinate x[i] at local coordinate s
 ///at previous timestep t (t=0: present; t>0: previous timestep).
 /// Overloaded to get information from bulk.
 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);
  }

 /// \short Return FE interpolated position x[] at local coordinate s as Vector
 /// Overloaded to get information from bulk.
 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);
  }

 /// \short Return FE interpolated position x[] at local coordinate s
 /// at previous timestep t as Vector (t=0: present; t>0: previous timestep).
 /// Overloaded to get information from bulk.
 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);
  }

 /// \short Return t-th time-derivative of the
 /// i-th FE-interpolated Eulerian coordinate at
 /// local coordinate s. Overloaded to get information from bulk.
 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);
  }

 /// \short Compte t-th time-derivative of the
 /// FE-interpolated Eulerian coordinate vector at
 /// local coordinate s.  Overloaded to get information from bulk.
 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);
  }

 /// \short Sign of outer unit normal (relative to cross-products of tangent
 /// vectors in the corresponding "bulk" element.
 int& normal_sign(){return Normal_sign;}

 /// \short Return sign of outer unit normal (relative to cross-products 
 /// of tangent vectors in the corresponding "bulk" element. (const version)
 int normal_sign() const {return Normal_sign;}
 
 /// \short Index of the face (a number that uniquely identifies the face
 /// in the element)
 int& face_index(){return Face_index;}

 /// \short Index of the face (a number that uniquely identifies the face
 /// in the element) (const version)
 int face_index() const {return Face_index;}

 /// \short Public access function for the tangent direction pointer.
 const Vector<double> * tangent_direction_pt() const
 {
  return Tangent_direction_pt;
 }

 /// \short Set the tangent direction vector.
 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;
 }

 /// \short Turn on warning for when there may be discontinuous tangent vectors
 /// from continuous_tangent_and_outer_unit_normal(...)
 void turn_on_warning_for_discontinuous_tangent()
  {
    FaceElement::Ignore_discontinuous_tangent_warning = true;
  }

 /// \short Turn off warning for when there may be discontinuous tangent 
 /// vectors from continuous_tangent_and_outer_unit_normal(...)
 void turn_off_warning_for_discontinuous_tangent()
  {
    FaceElement::Ignore_discontinuous_tangent_warning = false;
  }

 /// \short Compute the tangent vector(s) and the outer unit normal
 /// vector at the specified local coordinate.
 /// In two spatial dimensions, a "tangent direction" is not required.
 /// In three spatial dimensions, a tangent direction is required
 /// (set via set_tangent_direction(...)), and we project the tanent direction
 /// on to the surface. The second tangent vector is taken to be the cross
 /// product of the projection and the unit normal.
 void continuous_tangent_and_outer_unit_normal
  (const Vector<double> &s, Vector<Vector<double> > &tang_vec, 
   Vector<double> &unit_normal) const;

 /// \short Compute the tangent vector(s) and the outer unit normal
 /// vector at the ipt-th integration point. This is a wrapper around
 /// continuous_tangent_and_outer_unit_normal(...) with the integration points
 /// converted into local coordinates.
 void continuous_tangent_and_outer_unit_normal
  (const unsigned &ipt, Vector<Vector<double> > &tang_vec, 
   Vector<double> &unit_normal) const;

 /// \short Compute outer unit normal at the specified local coordinate
 void outer_unit_normal(const Vector<double> &s, 
                        Vector<double> &unit_normal) const;
 
 /// \short Compute outer unit normal at ipt-th integration point
 void outer_unit_normal(const unsigned &ipt, 
                        Vector<double> &unit_normal) const;

 /// Pointer to higher-dimensional "bulk" element
 FiniteElement*& bulk_element_pt() {return Bulk_element_pt;}


 /// Pointer to higher-dimensional "bulk" element (const version)
 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
 
 /// \short Return the pointer to the function that maps the face
 /// coordinate to the bulk coordinate
 CoordinateMappingFctPt &face_to_bulk_coordinate_fct_pt()
  {return Face_to_bulk_coordinate_fct_pt;}

 /// \short Return the pointer to the function that maps the face
 /// coordinate to the bulk coordinate (const version)
 CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt() const
  {return Face_to_bulk_coordinate_fct_pt;}


 /// \short Return the pointer to the function that returns the derivatives
 /// of the bulk coordinates wrt the face coordinates
 BulkCoordinateDerivativesFctPt &bulk_coordinate_derivatives_fct_pt()
  {return Bulk_coordinate_derivatives_fct_pt;}

 /// \short Return the pointer to the function that returns the derivatives
 /// of the bulk coordinates wrt the face coordinates (const version)
 BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt() const
  {return Bulk_coordinate_derivatives_fct_pt;}

#ifdef __clang__
#pragma clang diagnostic pop
#endif
 
 /// \short Return vector of local coordinates in bulk element, 
 /// given the local coordinates in this FaceElement
 Vector<double> local_coordinate_in_bulk(const Vector<double>& s) const;

 /// \short Calculate the vector of local coordinate in the bulk element
 /// given the local coordinates in this FaceElement
 void get_local_coordinate_in_bulk(const Vector<double> &s,
                                   Vector<double> &s_bulk) const;
 
 /// \short Calculate the derivatives of the local coordinates in the
 /// bulk element with respect to the local coordinates in this FaceElement.
 /// In addition return the index of a bulk local coordinate that varies away
 /// from the face.
 void get_ds_bulk_ds_face(const Vector<double> &s,
                          DenseMatrix<double> &dsbulk_dsface,
                          unsigned &interior_direction) const;

 /// \short Return the position type in the "bulk" element that corresponds
 /// to position type i on the FaceElement.
 unsigned& bulk_position_type(const unsigned &i) 
  {return Bulk_position_type[i];}

 /// \short Return the position type in the "bulk" element that corresponds
 /// to the position type i on the FaceElement. Const version
 const unsigned &bulk_position_type(const unsigned &i) const 
  {return Bulk_position_type[i];}

 /// \short Resize the storage for the bulk node numbers
 void bulk_node_number_resize(const unsigned &i) {Bulk_node_number.resize(i);}
 
 /// \short Return the bulk node number that corresponds to the n-th
 /// local node number
 unsigned &bulk_node_number(const unsigned &n) {return Bulk_node_number[n];}
 
 /// \short Return the bulk node number that corresponds to the n-th
 /// local node number (const version)
 const unsigned &bulk_node_number(const unsigned &n) const
  {return Bulk_node_number[n];}

 /// \short Resize the storage for bulk_position_type to i entries.
 void bulk_position_type_resize(const unsigned &i) 
  {Bulk_position_type.resize(i);}

 /// \short Return the number of values originally stored at local node n
 /// (before the FaceElement added additional values to it (if it did))
 unsigned& nbulk_value(const unsigned &n) {return Nbulk_value[n];}

 /// \short Return the number of values originally stored at local node n
 /// (before the FaceElement added additional values to it (if it did))
 /// (const version)
 unsigned nbulk_value(const unsigned &n) const 
  {return Nbulk_value[n];}

 /// \short  Resize the storage for the number of values originally 
 /// stored at the local nodes to i entries.
 void nbulk_value_resize(const unsigned &i) {Nbulk_value.resize(i);}

 /// \short Provide additional storage for a specified number of 
 /// values at the nodes of the FaceElement. (This is needed, for
 /// instance, in free-surface elements, if the non-penetration
 /// condition is imposed by Lagrange multipliers whose values
 /// are only stored at the surface nodes but not in the interior of
 /// the bulk element). \c nadditional_data_values[n] specifies the
 /// number of additional values required at node \c n of the 
 /// FaceElement. \b Note: Since this function is executed separately
 /// for each FaceElement, nodes that are common to multiple elements
 /// might be resized repeatedly. To avoid this, we only allow
 /// a single resize operation by comparing the number of values stored
 /// at each node to the number of values the node had when it 
 /// was simply a member of the associated "bulk"
 /// element. <b>There are cases where this will break!</b> -- e.g. if
 /// a node is common to two FaceElements which require 
 /// additional storage for distinct quantities. Such cases
 /// need to be handled by "hand-crafted" face elements.
 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
  }

 /// Output boundary coordinate zeta
 void output_zeta(std::ostream &outfile, const unsigned& nplot);
 
};



//========================================================================
/// SolidFaceElements combine FaceElements and SolidFiniteElements
/// and overload various functions so they work properly in the
/// FaceElement context
//========================================================================
class SolidFaceElement: public virtual FaceElement,
                        public virtual SolidFiniteElement
{

  public:

 /// \short The "global" intrinsic coordinate of the element when
 /// viewed as part of a geometric object should be given by
 /// the FaceElement representation, by default
 /// This final over-ride is required because both SolidFiniteElements 
 /// and FaceElements overload zeta_nodal
 double zeta_nodal(const unsigned &n, const unsigned &k,           
                          const unsigned &i) const 
  {return FaceElement::zeta_nodal(n,k,i);}     



 /// \short Return i-th FE-interpolated Lagrangian coordinate xi[i] at
 /// local coordinate s. Overloaded from SolidFiniteElement. Note that
 /// the Lagrangian coordinates are those defined in the bulk!
 /// For instance, in a 1D FaceElement that is aligned with
 /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary 
 /// in the FaceElement. This may confuse you if you (wrongly!) believe that 
 /// in a 1D SolidElement there should only a single Lagrangian 
 /// coordinate, namely xi_0!
 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);
  }


 /// \short Compute FE interpolated Lagrangian coordinate vector xi[] at 
 /// local coordinate s as Vector. Overloaded from SolidFiniteElement. Note 
 /// that the Lagrangian coordinates are those defined in the bulk!
 /// For instance, in a 1D FaceElement that is aligned with
 /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary 
 /// in the FaceElement. This may confuse you if you (wrongly!) believe that 
 /// in a 1D SolidElement there should only a single Lagrangian 
 /// coordinate, namely xi_0!
 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);
  }


};



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



//=======================================================================
/// Solid point element
//=======================================================================
 class SolidPointElement : public virtual SolidFiniteElement,
                           public virtual PointElement
{
};


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





//=======================================================================
/// FaceGeometry class definition: This policy class is used to allow 
/// construction of face elements that solve arbitrary equations without 
/// having to tamper with the corresponding "bulk" elements. 
/// The geometrical information for the face element must be specified
/// by each "bulk" element using an explicit specialisation of this class.
//=======================================================================
template<class ELEMENT> 
class FaceGeometry 
{

};


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



//======================================================================
/// Dummy FaceElement for use with purely geometric operations
/// such as mesh generation.
//======================================================================
template <class ELEMENT>
class DummyFaceElement : public virtual FaceGeometry<ELEMENT>, 
public virtual FaceElement
{

public:
 
 /// \short Constructor, which takes a "bulk" element and the 
 /// face index
 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);
  }


 /// \short Constructor 
 DummyFaceElement() : 
  FaceGeometry<ELEMENT>(), FaceElement()
  {  
  }


  /// \short The "global" intrinsic coordinate of the element when
  /// viewed as part of a geometric object should be given by
  /// the FaceElement representation, by default
  double zeta_nodal(const unsigned &n, const unsigned &k,           
                          const unsigned &i) const 
  {return FaceElement::zeta_nodal(n,k,i);}     

 /// Output nodal coordinates
 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;
    }     
  }

 /// Output at n_plot points
 void output(std::ostream &outfile, const unsigned &n_plot)
  {FiniteElement::output(outfile,n_plot);}

 /// C-style output 
 void output(FILE* file_pt)
  {FiniteElement::output(file_pt);}

 /// C_style output at n_plot points
 void output(FILE* file_pt, const unsigned &n_plot)
  {FiniteElement::output(file_pt,n_plot);}

}; 

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


//=========================================================================
/// Base class for elements that can specify a drag and torque
/// (about the origin) -- typically used for immersed particle 
/// computations.
//=========================================================================
class ElementWithDragFunction
{

public:
 
 /// Empty constructor
 ElementWithDragFunction(){}

 /// Empty virtual destructor
 virtual ~ElementWithDragFunction(){}

 /// \short Function that specifies the drag force and the torque about
 /// the origin
 virtual void get_drag_and_torque(Vector<double>& drag_force,
                                  Vector<double>& drag_torque)=0;

};


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


//======================================================================
/// Basic-ified FaceElement, without any of the functionality of
/// of actual FaceElements -- it's just a  surface element of the
/// same geometric type as the FaceGeometry associated with
/// bulk element specified by the template parameter. The element
/// can be used to represent boundaries without actually being 
/// attached to a bulk element. Used mainly during unstructured
/// mesh generation. 
//======================================================================
template <class ELEMENT>
class FreeStandingFaceElement : public virtual FaceGeometry<ELEMENT>
{

public:
 
 /// \short Constructor
 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
   
  }

 /// Access function for the boundary number in bulk mesh
 inline const unsigned& boundary_number_in_bulk_mesh() const
  {return Boundary_number_in_bulk_mesh;}
 
 
 /// Set function for the 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
  }

 /// \short In a FaceElement, the "global" intrinsic coordinate
 /// of the element along the boundary, when viewed as part of
 /// a compound geometric object is specified using the
 /// boundary coordinate defined by the mesh. 
 /// Note: Boundary coordinates will have been set up when
 /// creating the underlying mesh, and their values will have 
 /// been stored at the nodes.
 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:
 
 /// The boundary number in the bulk mesh to which this element is attached
 unsigned Boundary_number_in_bulk_mesh;
 
#ifdef PARANOID
 
 /// \short Has the Boundary_number_in_bulk_mesh been set? Only included if
 /// compiled with PARANOID switched on.
 bool Boundary_number_in_bulk_mesh_has_been_set;

#endif

}; 



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


//======================================================================
/// Pure virtual base class for elements that can be used with
/// PressureBasedSolidLSCPreconditioner.
//======================================================================
class SolidElementWithDiagonalMassMatrix 
{
 
  public:
 
 /// Empty constructor
 SolidElementWithDiagonalMassMatrix(){}
 
 /// Virtual destructor 
 virtual ~SolidElementWithDiagonalMassMatrix(){}
 
 /// Broken copy constructor
 SolidElementWithDiagonalMassMatrix(
  const SolidElementWithDiagonalMassMatrix&) 
  { 
   BrokenCopy::broken_copy(
    "SolidElementWithDiagonalMassMatrix");
  } 
 
 /// Broken assignment operator
 void operator=(const 
                SolidElementWithDiagonalMassMatrix&) 
  {
   BrokenCopy::broken_assign(
    "SolidElementWithDiagonalMassMatrix");
  }
 
 /// \short Get the diagonal of whatever represents the mass matrix
 /// in the specific preconditionable element. For Navier-Stokes
 /// elements this is the velocity mass matrix; for incompressible solids it's 
 /// the mass matrix; ... 
 virtual void get_mass_matrix_diagonal(Vector<double> &mass_diag)=0;
 
};




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



//======================================================================
/// Pure virtual base class for elements that can be used with
/// Navier-Stokes Schur complement preconditioner and provide the diagonal 
/// of their velocity and pressure mass matrices -- needs to be defined here 
/// (in generic)  because this applies to a variety of Navier-Stokes 
/// elements (cartesian, cylindrical polar, ...)
/// that can be preconditioned effectively by the Navier Stokes (!)
/// preconditioners in the (cartesian) Navier-Stokes directory.
//======================================================================
class NavierStokesElementWithDiagonalMassMatrices
{
 
  public:
 
 /// Empty constructor
 NavierStokesElementWithDiagonalMassMatrices(){}
 
 /// Virtual destructor 
 virtual ~NavierStokesElementWithDiagonalMassMatrices(){}
 
 /// Broken copy constructor
 NavierStokesElementWithDiagonalMassMatrices(
  const NavierStokesElementWithDiagonalMassMatrices&) 
  { 
   BrokenCopy::broken_copy(
    "NavierStokesElementWithDiagonalMassMatrices");
  } 
 
 /// Broken assignment operator
 void operator=(const 
                NavierStokesElementWithDiagonalMassMatrices&) 
  {
   BrokenCopy::broken_assign(
    "NavierStokesElementWithDiagonalMassMatrices");
  }
 
 /// \short Compute the diagonal of the velocity/pressure mass matrices.
 /// If which one=0, both are computed, otherwise only the pressure 
 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the 
 /// LSC version of the preconditioner only needs that one)
 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;
 
};




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

} // end of oomph-lib namespace

#endif