elements.cc

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//====================================================================
//Non-inline member functions for generic elements

#include<float.h>

//oomph-lib includes
#include "elements.h"
#include "timesteppers.h"
#include "integral.h"
#include "shape.h"
#include "oomph_definitions.h"
#include "element_with_external_element.h"

namespace oomph
{

/// Static boolean to suppress warnings about repeated internal
/// data. Defaults to false
bool GeneralisedElement::Suppress_warning_about_repeated_internal_data=false;


/// Static boolean to suppress warnings about repeated external
/// data. Defaults to true
bool GeneralisedElement::Suppress_warning_about_repeated_external_data=true;


///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
//  Functions for generalised elements
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////

//=======================================================================
/// 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()
//=======================================================================
 unsigned GeneralisedElement::add_internal_data(Data* const &data_pt,
                                                const bool &fd)
 {
  //Local cache of numbers of internal and external data
  const unsigned n_internal_data = Ninternal_data;
  const unsigned n_external_data = Nexternal_data;

  //Find out whether the data is already stored in the array
  
  //Loop over the number of internals
  //The internal data are stored at the beginning of the array
  for(unsigned i=0;i<n_internal_data;i++)
   {
    //If the passed pointer is stored in the array
    if(internal_data_pt(i) == data_pt)
     {
#ifdef PARANOID
      if (!Suppress_warning_about_repeated_internal_data)
       {
        oomph_info << std::endl << std::endl;
        oomph_info 
         << "-----------------------------------------------------------------"
         << std::endl
         << "Info: Data is already included in this element's internal storage." 
         << std::endl
         << "It's stored as entry " << i <<  " and I'm not adding it again." 
         << std::endl<< std::endl
         << "Note: You can suppress this message by recompiling without"
         << "\n      PARANOID or setting the boolean \n"
         << "\n      GeneralisedElement::Suppress_warning_about_repeated_internal_data"
         << "\n\n      to true."
         << std::endl
         << "-----------------------------------------------------------------"
         << std::endl << std::endl;
       }
#endif
      //Return the index to the data object
      return i;
     }
   }
  
  //Allocate new storage for the pointers to data
  Data** new_data_pt = new Data*[n_internal_data + n_external_data + 1];

  //Copy the old internal values across to the beginning of the array
  for(unsigned i=0;i<n_internal_data;i++) {new_data_pt[i] = Data_pt[i];}
  
  //Now add the new value to the end of the internal data
  new_data_pt[n_internal_data] = data_pt;
  
  //Copy the external values across
  for(unsigned i=0;i<n_external_data;i++)
   {new_data_pt[n_internal_data + 1 + i] = Data_pt[n_internal_data + i];}
  
  //Delete the storage associated with the previous values
  delete[] Data_pt;
  
  //Set the pointer to the new storage
  Data_pt = new_data_pt;
  
  //Resize the array of boolean flags
  Data_fd.resize(n_internal_data + n_external_data + 1);
  //Shuffle the flags for the external data to the end of the array
  for(unsigned i=n_external_data;i>0;i--)
   {
    Data_fd[n_internal_data + i] = Data_fd[n_internal_data + i-1];
   }
  //Now add the new flag to the end of the internal data
  Data_fd[n_internal_data] = fd;
  
  //Increase the number of internals
  ++Ninternal_data;
  
  //Return the final index to the new internal data
  return n_internal_data;
 }
 
//=======================================================================
/// Add the contents of the queue global_eqn_numbers to
/// the local storage for the equation numbers, which represents the 
/// local-to-global translation scheme. It is essential that the entries
/// are added in order, i.e. from the front.
//=======================================================================
 void GeneralisedElement::add_global_eqn_numbers(
  std::deque<unsigned long> const &global_eqn_numbers,
  std::deque<double*> const &global_dof_pt)
 {
  //Find the number of dofs
  const unsigned n_dof = Ndof;
  //Find the number of additional dofs
  const unsigned n_additional_dof = global_eqn_numbers.size();
  //If there are none, return immediately
  if(n_additional_dof==0) {return;}
  
  //Find the new total number of equation numbers
  const unsigned new_n_dof = n_dof + n_additional_dof;
  //Create storage for all equations, initialised to NULL
  unsigned long *new_eqn_number = new unsigned long[new_n_dof];

  //Copy over the existing values to the start new storage
  for(unsigned i=0;i<n_dof;i++) {new_eqn_number[i] = Eqn_number[i];}
  
  //Set an index to the next position in the new storage
  unsigned index = n_dof;
  //Loop over the queue and add it's entries to our new storage
  for(std::deque<unsigned long>::const_iterator it=global_eqn_numbers.begin();
      it!=global_eqn_numbers.end();++it)
   {
    //Add the value to the storage
    new_eqn_number[index] = *it;
    //Increase the array index
    ++index;
   }
  

  //If a non-empty dof deque has been passed then do stuff
  const unsigned n_additional_dof_pt = global_dof_pt.size();
  if(n_additional_dof_pt > 0)
   {
//If it's size is not the same as the equation numbers complain
#ifdef PARANOID 
    if(n_additional_dof_pt != n_additional_dof)
     {
      std::ostringstream error_stream;
      error_stream 
       << 
       "global_dof_pt is non-empty, yet it does not have the same size\n"
       <<
       "as global_eqn_numbers.\n"
       << 
       "There are " << n_additional_dof << " equation numbers,\n"
       << 
       "but " << n_additional_dof_pt << std::endl;
      
      throw OomphLibError(error_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
#endif

    //Create storge for all dofs initialised to NULL
    double **new_dof_pt = new double*[new_n_dof];
    //Copy over the exisiting values to the start of new storage
    for(unsigned i=0;i<n_dof;i++) {new_dof_pt[i] = Dof_pt[i];}
  
    //Set an index to the next position in the new storage
    unsigned index = n_dof;
    //Loop over the queue and add it's entries to our new storage
    for(std::deque<double*>::const_iterator it=global_dof_pt.begin();
        it!=global_dof_pt.end();++it)
     {
      //Add the value to the storage
      new_dof_pt[index] = *it;
      //Increase the array index
      ++index;
     }

    //Now delete the old storage
    delete[] Dof_pt;
    //Set the pointer to address the new storage
    Dof_pt = new_dof_pt;
   }

  //Now delete the old for the equation numbers storage
  delete[] Eqn_number;
  //Set the pointer to address the new storage
  Eqn_number = new_eqn_number; 
  //Finally update the number of degrees of freedom
  Ndof = new_n_dof;
 }

//========================================================================
/// Empty dense matrix used as a dummy argument to combined
/// residual and jacobian functions in the case when only the residuals
/// are being assembled
//========================================================================
DenseMatrix<double> GeneralisedElement::Dummy_matrix;

//========================================================================
/// Static storage used when pointers to the dofs are being assembled by
/// add_global_eqn_numbers()
//========================================================================
std::deque<double*> GeneralisedElement::Dof_pt_deque;

 
//=========================================================================
/// Default value used as the increment for finite difference calculations
/// of the jacobian matrices
//=========================================================================
 double GeneralisedElement::Default_fd_jacobian_step=1.0e-8;

//==========================================================================
/// Destructor for generalised elements: Wipe internal data. Pointers 
/// to external data get NULLed  but are not deleted because they 
/// are (generally) shared by lots of elements.
//==========================================================================
 GeneralisedElement::~GeneralisedElement()
 {
  //Delete each of the objects stored as internal data
  for(unsigned i=Ninternal_data;i>0;i--)
   {
    //The objects are stored at the beginning of the Data_pt array
    delete Data_pt[i-1];
    Data_pt[i-1] = 0;
   }

  //Now delete the storage for internal and external data
  delete[] Data_pt;
  
  //Now if we have allocated storage for the local equation for
  //the internal and external data, delete it.
  if(Data_local_eqn)
   {
    delete[] Data_local_eqn[0];
    delete[] Data_local_eqn;
   }

  //Delete the storage for the global equation numbers
  delete[] Eqn_number;
 }


//=======================================================================
/// Add a (pointer to an) external data object to the element and
/// return the index required to obtain it from the access
/// function \c external_data_pt()
//=======================================================================
 unsigned GeneralisedElement::add_external_data(Data* const &data_pt,
                                                const bool &fd)
 {
  //Find the numbers of internal and external data
  const unsigned n_internal_data = Ninternal_data;
  const unsigned n_external_data = Nexternal_data;
  //Find out whether the data is already stored in the array
  
  //Loop over the number of externals
  //The external data are stored at the end of the array Data_pt
  for(unsigned i=0;i<n_external_data;i++)
   {
    //If the passed pointer is stored in the array
    if(external_data_pt(i) == data_pt)
     {
#ifdef PARANOID
      if (!Suppress_warning_about_repeated_external_data)
       {
        oomph_info << std::endl << std::endl;
        oomph_info 
         << "-----------------------------------------------------------------"
         << std::endl
         << "Info: Data is already included in this element's external storage." 
         << std::endl
         << "It's stored as entry " << i <<  " and I'm not adding it again" 
         << std::endl << std::endl
         << "Note: You can suppress this message by recompiling without"
         << "\n      PARANOID or setting the boolean \n"
         << "\n      GeneralisedElement::Suppress_warning_about_repeated_external_data"
         << "\n\n      to true."
         << std::endl
         << "-----------------------------------------------------------------"
         << std::endl << std::endl;
       }
#endif
      //Return the index to the data object
      return i;
     }
   }
  
  //Allocate new storage for the pointers to data
  Data** new_data_pt = new Data*[n_internal_data + n_external_data + 1];
  
  //Copy the old internal and external values across to the new array
  for(unsigned i=0;i<(n_internal_data + n_external_data);i++) 
   {new_data_pt[i] = Data_pt[i];}
  
  //Add the new data pointer to the end of the array
  new_data_pt[n_internal_data + n_external_data] = data_pt;
  
  //Delete the storage associated with the previous values
  delete[] Data_pt;
  
  //Set the pointer to the new storage
  Data_pt = new_data_pt;
  
  //Resize the array of boolean flags
  Data_fd.resize(n_internal_data + n_external_data + 1);
  //Now add the new flag to the end of the external data
  Data_fd[n_internal_data + n_external_data] = fd;

  //Increase the number of externals
  ++Nexternal_data;
 
  //Return the final index to the new external data
  return n_external_data;
 }

//========================================================================
/// Flush all the external data, i.e. clear the pointers to external
/// data from the internal storage.
//========================================================================
 void GeneralisedElement::flush_external_data()
 {
  //Get the numbers of internal and external data
  const unsigned n_external_data = Nexternal_data;
  //If there is external data
  if(n_external_data>0)
   {
    //Storage for new data, initialised to NULL
    Data** new_data_pt=0;

    //Find the number of internal data
    const unsigned n_internal_data = Ninternal_data;
    //If there is internal data resize Data_pt and copy over
    //the pointers
    if(n_internal_data  > 0)
     {
      //The new data pointer should only be the size of the internal data
      new_data_pt = new Data*[n_internal_data];
      //Copy over the internal data only
      for(unsigned i=0;i<n_internal_data;i++) {new_data_pt[i] = Data_pt[i];}
     }
    
    //Delete the old storage
    delete[] Data_pt;
    //Set the new storage, this will be NULL if there is no internal data
    Data_pt = new_data_pt;
    //Set the number of externals to zero
    Nexternal_data = 0;

    //Resize the array of boolean flags to the number of internals
    Data_fd.resize(n_internal_data);
   }
 }



//=========================================================================
/// Remove the object addressed by data_pt from the external data array
/// Note that this could mess up the numbering of other external data
//========================================================================
 void GeneralisedElement::flush_external_data(Data* const &data_pt)
 {
  //Get the current numbers of external data
  const unsigned n_external_data = Nexternal_data;
  //Index of the data to be removed
  //We initialise this to be n_external_data, and it will be smaller than
  //n_external_data if the data pointer is found in the array
  unsigned index = n_external_data;
  //Loop over the external data and find the argument
  for(unsigned i=0;i<n_external_data;i++)
   {
    //If we have found the data pointer, set the index and break
    if(external_data_pt(i) == data_pt)
     {
      index=i;
      break;
     }
   }

  //If we have found an index less than Nexternal_data, then we have found
  //the data in the array
  if(index < n_external_data)
   {
    //Initialise a new array to NULL
    Data** new_data_pt = 0;
    
    //Find the number of internals
    const unsigned n_internal_data = Ninternal_data;

    //Find the new number of total data items (one fewer)
    const unsigned n_total_data = n_internal_data + n_external_data - 1;
 
   //Create a new array containing the data items, if we have any
    if(n_total_data > 0) {new_data_pt = new Data*[n_total_data];}

    //Copy over all the internal data
    for(unsigned i=0;i<n_internal_data;i++) {new_data_pt[i] = Data_pt[i];}

    //Now copy over the un-flushed data
    unsigned counter=0;
    for(unsigned i=0;i<n_external_data;i++)
     {
      //If we are not at the deleted index 
      if(i!=index)
       {
        //Copy the undeleted entry into the next entry in our new 
        //array of pointers to Data
        new_data_pt[n_internal_data + counter] = Data_pt[i];
        //Increase the counter
        ++counter;
       }
     }

    //Delete the storage associated with the previous values
    delete[] Data_pt;

    //Set pointers to the new storage, will be NULL if no data left
    Data_pt = new_data_pt;

    //Remove the entry from the array of boolean flags
    Data_fd.erase(Data_fd.begin()+n_internal_data+index);

    //Decrease the number of externals
    --Nexternal_data;
    
    //Issue a warning if there will be external data remaining
    if(Nexternal_data > 1)
     {
      std::ostringstream warning_stream;
      warning_stream << "Data removed from element's external data   "
                     << std::endl
                     << "You may have to update the indices for remaining data "
                     << std::endl
                     << "This can be achieved by using add_external_data()    "
                     << std::endl;
      OomphLibWarning(warning_stream.str(),
                      "GeneralisedElement::flush_external_data()",
                      OOMPH_EXCEPTION_LOCATION);
     }
   }
 }



//==========================================================================
/// This function loops over the internal data of the element and assigns
/// GLOBAL equation numbers to the data objects.
///
/// Pass: 
/// - the current total number of dofs, global_number, which gets
///   incremented. 
/// - Dof_pt, the Vector of pointers to the global dofs 
///   (to which any internal dofs are added).
/// .
//==========================================================================
 void GeneralisedElement::
 assign_internal_eqn_numbers(unsigned long &global_number,
                             Vector<double *> &Dof_pt)
 {
  //Loop over the internal data and assign the equation numbers
  //The internal data are stored at the beginning of the Data_pt array
  for(unsigned i=0;i<Ninternal_data;i++)
   {internal_data_pt(i)->assign_eqn_numbers(global_number,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 GeneralisedElement::describe_dofs(std::ostream& out,
                                        const std::string& current_string) const
 {
  for(unsigned i=0;i<Ninternal_data;i++)
   {
    std::stringstream conversion;
    conversion <<" of Internal Data "<<i<<current_string;
    std::string in(conversion.str());
    internal_data_pt(i)->describe_dofs(out,in);
   }
 }

//==========================================================================
/// \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 GeneralisedElement::
describe_local_dofs(std::ostream& out,const std::string& current_string) const
 {
  // Find the number of internal and external data
  const unsigned n_internal_data = Ninternal_data;
  const unsigned n_external_data = Nexternal_data;

  // Now loop over the internal data and describe local equation numbers
  for(unsigned i=0;i<n_internal_data;i++)
   {
    // Pointer to the internal data
    Data* const data_pt = internal_data_pt(i);
      
    std::stringstream conversion;
    conversion <<" of Internal Data "<<i<<current_string;
    std::string in(conversion.str());
    data_pt->describe_dofs(out,in);
   }   

    
  // Now loop over the external data and assign local equation numbers
  for(unsigned i=0;i<n_external_data;i++)
   {
    // Pointer to the external data
    Data* const data_pt = external_data_pt(i);
          
    std::stringstream conversion;
    conversion<<" of External Data "<<i<<current_string;
    std::string in(conversion.str());
    data_pt->describe_dofs(out,in);  
   }
   
 }

 //==========================================================================
 /// This function loops over the internal data of the element and add 
 /// pointers to their unconstrained values to a map indexed by the global
 /// equation number.
 //==========================================================================
 void GeneralisedElement:: add_internal_value_pt_to_map(
  std::map<unsigned,double*> &map_of_value_pt)
 {
  //Loop over the internal data and add their data to the map
  //The internal data are stored at the beginning of the Data_pt array
  for(unsigned i=0;i<Ninternal_data;i++)
   {internal_data_pt(i)->add_value_pt_to_map(map_of_value_pt);}
 }




#ifdef OOMPH_HAS_MPI
 //=========================================================================
 /// Add all internal data and time history values to the vector in 
 /// the internal storage order
 //=========================================================================
 void GeneralisedElement::add_internal_data_values_to_vector(
  Vector<double> &vector_of_values)
 {
  for(unsigned i=0;i<Ninternal_data;i++)
   {internal_data_pt(i)->add_values_to_vector(vector_of_values);}
 }

 //========================================================================
 /// 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 GeneralisedElement::read_internal_data_values_from_vector(
  const Vector<double> &vector_of_values, 
  unsigned &index)
 {
  for(unsigned i=0;i<Ninternal_data;i++)
   {internal_data_pt(i)->read_values_from_vector(vector_of_values,index);}
 }

 //=========================================================================
 /// Add all equation numbers associated with internal data 
 /// to the vector in the internal storage order
 //=========================================================================
 void GeneralisedElement::add_internal_eqn_numbers_to_vector(
  Vector<long> &vector_of_eqn_numbers)
 {
  for(unsigned i=0;i<Ninternal_data;i++)
   {internal_data_pt(i)->add_eqn_numbers_to_vector(vector_of_eqn_numbers);}
 }

 //=========================================================================
 ///  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 GeneralisedElement::read_internal_eqn_numbers_from_vector(
  const Vector<long> & vector_of_eqn_numbers, unsigned &index)
 {
  for(unsigned i=0;i<Ninternal_data;i++)
   {internal_data_pt(i)->read_eqn_numbers_from_vector(vector_of_eqn_numbers,
                                                      index);}
 }
 
#endif


 //====================================================================
 /// Setup the arrays of local equation numbers for the element. 
 /// In general, this function should not need to be overloaded. Instead
 /// the function assign_all_generic_local_eqn_numbers() will be overloaded
 /// by different types of Element.
 /// That said, the function is virtual so that it 
 /// may be overloaded by the user if they *really* know what they are doing.
 //==========================================================================
 void GeneralisedElement::assign_local_eqn_numbers(
  const bool &store_local_dof_pt) 
 {
  clear_global_eqn_numbers();
  assign_all_generic_local_eqn_numbers(store_local_dof_pt);
  assign_additional_local_eqn_numbers();

  //Check that no global equation numbers are repeated
#ifdef PARANOID

  std::ostringstream error_stream;

  //Loop over the array of equation numbers and add to set to assess
  //uniqueness
  std::map<unsigned,bool> is_repeated;
  std::set<unsigned long> set_of_global_eqn_numbers;
  unsigned old_ndof=0;
  for(unsigned n=0;n<Ndof;++n) 
   {
    set_of_global_eqn_numbers.insert(Eqn_number[n]);
    if (set_of_global_eqn_numbers.size()==old_ndof)
     {
      error_stream << "Repeated global eqn: " << Eqn_number[n] << std::endl;
      is_repeated[Eqn_number[n]]=true;
     }
    old_ndof=set_of_global_eqn_numbers.size();
   }
  
  //If the sizes do not match we have a repeat, throw an error
  if(set_of_global_eqn_numbers.size() != Ndof)
   {
#ifdef OOMPH_HAS_MPI
    error_stream << "Element is ";
    if (!is_halo()) error_stream << "not ";
    error_stream << "a halo element\n\n";
#endif
    error_stream << "\nLocal/lobal equation numbers: " << std::endl;
    for(unsigned n=0;n<Ndof;++n)
     {
      error_stream << "  " << n << "        " << Eqn_number[n] << std::endl;
     }
    
    // It's helpful for debugging purposes to output more about this element
    error_stream << std::endl << " element address is " << this << std::endl;
    
    // Check if the repeated dofs are among the internal Data values
    unsigned nint=this->ninternal_data();
    error_stream << "Number of internal data " << nint << std::endl;
    for (unsigned i=0;i<nint;i++)
     {
      Data* data_pt=this->internal_data_pt(i);
      unsigned nval=data_pt->nvalue();
      for (unsigned j=0;j<nval;j++)
       {
        int eqn_no=data_pt->eqn_number(j);
        error_stream << "Internal dof: " << eqn_no << std::endl;
        if (is_repeated[unsigned(eqn_no)])
         {
          error_stream << "Repeated internal dof: " << eqn_no << std::endl;
         }
       }
     }
    
    // Check if the repeated dofs are among the external Data values
    unsigned next=this->nexternal_data();
    error_stream << "Number of external data " << next << std::endl;
    for (unsigned i=0;i<next;i++)
     {
      Data* data_pt=this->external_data_pt(i);
      unsigned nval=data_pt->nvalue();
      for (unsigned j=0;j<nval;j++)
       {
        int eqn_no=data_pt->eqn_number(j);
        error_stream << "External dof: " << eqn_no << std::endl;
        if (is_repeated[unsigned(eqn_no)])
         {
          error_stream << "Repeated external dof: " << eqn_no;            
          Node* nod_pt=dynamic_cast<Node*>(data_pt);
          if (nod_pt!=0)
           {
            error_stream << " (is a node at: ";
            unsigned ndim=nod_pt->ndim();
            for (unsigned ii=0;ii<ndim;ii++)
             {
              error_stream << nod_pt->x(ii) << " ";
             }
           }
          error_stream << ")\n";
         }
       }
     }
    
    
    // If it's an element with external element check the associated
    // Data
    ElementWithExternalElement* e_el_pt=
     dynamic_cast<ElementWithExternalElement*>(this);
    if (e_el_pt!=0)
     {
      // Check if the repeated dofs are among the external Data values
      {
       unsigned next=e_el_pt->nexternal_interaction_field_data();
       error_stream << "Number of external element data " << next << std::endl;
       Vector<Data*> data_pt(e_el_pt->external_interaction_field_data_pt());
       for (unsigned i=0;i<next;i++)
        {
         unsigned nval=data_pt[i]->nvalue();
         for (unsigned j=0;j<nval;j++)
          {
           int eqn_no=data_pt[i]->eqn_number(j);
           error_stream << "External element dof: " << eqn_no << std::endl;
           if (is_repeated[unsigned(eqn_no)])
            {
             error_stream << "Repeated external element dof: " 
                          << eqn_no;            
             Node* nod_pt=dynamic_cast<Node*>(data_pt[i]);
             if (nod_pt!=0)
              {
               error_stream << " (is a node at: ";
               unsigned ndim=nod_pt->ndim();
               for (unsigned ii=0;ii<ndim;ii++)
                {
                 error_stream << nod_pt->x(ii) << " ";
                }
              }
             error_stream << ")\n";
            }
          }
        }
      }
       
      
      // Check if the repeated dofs are among the external geom Data values
      {
       unsigned next=e_el_pt->nexternal_interaction_geometric_data();
       error_stream << "Number of external element geom data " 
                    << next << std::endl;
       Vector<Data*> data_pt(e_el_pt->
                             external_interaction_geometric_data_pt());
       for (unsigned i=0;i<next;i++)
        {
         unsigned nval=data_pt[i]->nvalue();
         for (unsigned j=0;j<nval;j++)
          {
           int eqn_no=data_pt[i]->eqn_number(j);
           error_stream << "External element geometric dof: " 
                        << eqn_no << std::endl;
           if (is_repeated[unsigned(eqn_no)])
            {
             error_stream << "Repeated external element geometric dof: " 
                          << eqn_no << " " 
                          << data_pt[i]->value(j);            
             Node* nod_pt=dynamic_cast<Node*>(data_pt[i]);
             if (nod_pt!=0)
              {
               error_stream << " (is a node at: ";
               unsigned ndim=nod_pt->ndim();
               for (unsigned ii=0;ii<ndim;ii++)
                {
                 error_stream << nod_pt->x(i) << " ";
                }
               error_stream << ")";
              }
             error_stream << "\n";
            }
          }
        }
      }

     }

    // If it's a FiniteElement then output its nodes
    FiniteElement* f_el_pt=dynamic_cast<FiniteElement*>(this);
    if (f_el_pt!=0)
     {
      unsigned n_node=f_el_pt->nnode();
      for (unsigned n=0;n<n_node;n++)
       {
        Node* nod_pt=f_el_pt->node_pt(n);
        unsigned nval=nod_pt->nvalue();
        for (unsigned j=0;j<nval;j++)
         {
          int eqn_no=nod_pt->eqn_number(j);
          error_stream << "Node " << n 
                       << ": Nodal dof: " 
                       << eqn_no;
          if (eqn_no>=0)
           {
            if (is_repeated[unsigned(eqn_no)])
             {
              error_stream << "Node " << n 
                           << ": Repeated nodal dof: " 
                           << eqn_no;
              if (j>=f_el_pt->required_nvalue(n))
               {
                error_stream << " (resized)";
               }
              error_stream << std::endl;
             }
           }
         }       
        SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
        if (solid_nod_pt!=0)
         {
          Data* data_pt=solid_nod_pt->variable_position_pt();
          unsigned nval=data_pt->nvalue();
          for (unsigned j=0;j<nval;j++)
           {
            int eqn_no=data_pt->eqn_number(j);
            error_stream <<  "Node " << n << ": Positional dof: " 
                         << eqn_no << std::endl;
            if (is_repeated[unsigned(eqn_no)])
             {
              error_stream << "Repeated positional dof: " 
                           << eqn_no << " " << data_pt->value(j) << std::endl;
             }
           }
         }        
       }
      
      // Output nodal coordinates of element
      n_node=f_el_pt->nnode();
      for (unsigned n=0;n<n_node;n++)
       {
        Node* nod_pt=f_el_pt->node_pt(n);
        unsigned n_dim=nod_pt->ndim();
        error_stream << "Node " << n << " at ( ";
        for (unsigned i=0;i<n_dim;i++)
         {
          error_stream << nod_pt->x(i) << " ";
         }
        error_stream << ")" << std::endl;
       }

     }


    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
 }
 

//==========================================================================
/// This function loops over the internal and external data of the element, 
/// adds the GLOBAL equation numbers to the local-to-global look-up scheme and 
/// fills in the look-up schemes for the local equation 
/// numbers.
/// If the boolean argument is true then pointers to the dofs will be 
/// stored in Dof_pt
//==========================================================================
 void GeneralisedElement::assign_internal_and_external_local_eqn_numbers(
  const bool &store_local_dof_pt)
 {
  //Find the number of internal and external data
  const unsigned n_internal_data = Ninternal_data;
  const unsigned n_external_data = Nexternal_data;
  //Find the total number of data
  const unsigned n_total_data = n_internal_data + n_external_data;

  //If there is data
  if(n_total_data  > 0)
   {
    //Find the number of local equations assigned so far
    unsigned local_eqn_number = ndof();

    //We need to find the total number of values stored in all the
    //internal and external data
    //Initialise to the number stored in the first data object
    unsigned n_total_values = Data_pt[0]->nvalue();
    //Loop over the other data and add the number of values stored
    for(unsigned i=1;i<n_total_data;++i)
     {n_total_values += Data_pt[i]->nvalue();}

    //If allocated delete the old storage
    if(Data_local_eqn) 
     {
      delete[] Data_local_eqn[0];
      delete[] Data_local_eqn;
     }
        
    //If there are no values then we are done, null out the storage and
    //return
    if(n_total_values==0) {Data_local_eqn=0; return;}

    //Allocate the storage for the local equation numbers
    //The idea is that we store internal equation numbers followed by
    //external equation numbers

    //Firstly allocate pointers to the rows for each data item
    Data_local_eqn = new int*[n_total_data];
    //Now allocate storage for all the equation numbers
    Data_local_eqn[0] = new int[n_total_values];
    //Set all values to be unclassified
    for(unsigned i=0;i<n_total_values;++i)
     {Data_local_eqn[0][i] = Data::Is_unclassified;}

    //Loop over the remaining rows and set their pointers
    for(unsigned i=1;i<n_total_data;++i)
     {
      //Initially set the pointer to the i-th row to the pointer
      //to the i-1th row
      Data_local_eqn[i] = Data_local_eqn[i-1];
      //Now increase the row pointer by the number of values 
      //stored at the i-1th data object
      Data_local_eqn[i] += Data_pt[i-1]->nvalue();
     }
    
    //A local queue to store the global equation numbers
    std::deque<unsigned long> global_eqn_number_queue;

    //Now loop over the internal data and assign local equation numbers
    for(unsigned i=0;i<n_internal_data;i++)
     {
      //Pointer to the internal data
      Data* const data_pt = internal_data_pt(i);
      //Find the number of values stored at the internal data
      unsigned n_value = data_pt->nvalue();
     
      //Loop over the number of values
      for(unsigned j=0;j<n_value;j++)
       {
        //Get the GLOBAL equation number
        long eqn_number = data_pt->eqn_number(j);
        //If the GLOBAL equation number is positive (a free variable)
        if(eqn_number >= 0)
         {
          //Add the GLOBAL equation number to the queue
          global_eqn_number_queue.push_back(eqn_number);
          //Add pointer to the dof to the queue if required
          if(store_local_dof_pt)
           {
            GeneralisedElement::Dof_pt_deque.push_back(data_pt->value_pt(j));
           }
          //Add the local equation number to the storage scheme
          Data_local_eqn[i][j] = local_eqn_number;
          //Increase the local number
          local_eqn_number++;
         }
        else
         {
          //Set the local scheme to be pinned
          Data_local_eqn[i][j] = Data::Is_pinned;
         } 
       }
     } //End of loop over internal data

    
    //Now loop over the external data and assign local equation numbers
    for(unsigned i=0;i<n_external_data;i++)
     {
      //Pointer to the external data
      Data* const data_pt = external_data_pt(i);
      //Find the number of values stored at the external data
      unsigned n_value = data_pt->nvalue();
     
      //Loop over the number of values
      for(unsigned j=0;j<n_value;j++)
       {
        //Get the GLOBAL equation number
        long eqn_number = data_pt->eqn_number(j);
        //If the GLOBAL equation number is positive (a free variable)
        if(eqn_number >= 0)
         {
          //Add the GLOBAL equation number to the queue
          global_eqn_number_queue.push_back(eqn_number);
          //Add pointer to the dof to the queue if required
          if(store_local_dof_pt)
           {
            GeneralisedElement::Dof_pt_deque.push_back(data_pt->value_pt(j));
           }
          //Add the local equation number to the local scheme
          Data_local_eqn[n_internal_data + i][j] = local_eqn_number;
          //Increase the local number
          local_eqn_number++;
         }
        else
         {
          //Set the local scheme to be pinned
          Data_local_eqn[n_internal_data + i][j] = Data::Is_pinned;
         }
       }
     }

    //Now add our global equations numbers to the internal element storage
    add_global_eqn_numbers(global_eqn_number_queue,
                           GeneralisedElement::Dof_pt_deque);
    //Clear the memory used in the deque
    if(store_local_dof_pt)
     {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}
   }
 }


//============================================================================
/// This function calculates the entries of Jacobian matrix, used in 
/// the Newton method, associated with the internal degrees of freedom.
/// It does this using finite differences, 
/// rather than an analytical formulation, so can be done in total generality.
/// 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 GeneralisedElement::
 fill_in_jacobian_from_internal_by_fd(Vector<double> &residuals,
                                      DenseMatrix<double> &jacobian,
                                      const bool &fd_all_data)
 {
  //Locally cache the number of internal data
  const unsigned n_internal_data = Ninternal_data;

  //If there aren't any internal data, then return straight away
  if(n_internal_data == 0) {return;}

  //Call the update function to ensure that the element is in
  //a consistent state before finite differencing starts
  update_before_internal_fd();

  //Find the number of dofs in the element
  const unsigned n_dof = ndof();

  //Create newres vector
  Vector<double> newres(n_dof);

  //Integer storage for local unknown
  int local_unknown=0;
  
  //Use the default finite difference step
  const double fd_step = Default_fd_jacobian_step;

  //Loop over the internal data
  for(unsigned i=0;i<n_internal_data;i++)
   {
    //If we are doing all finite differences or 
    //the variable is included in the finite difference loop, do it
    if(fd_all_data || internal_data_fd(i))
     {
      //Get the number of value at the internal data
      const unsigned n_value = internal_data_pt(i)->nvalue();
      
      //Loop over the number of values
      for(unsigned j=0;j<n_value;j++)
       {
        //Get the local equation number
        local_unknown = internal_local_eqn(i,j);
        //If it's not pinned
        if(local_unknown >= 0)
         {
          //Get a pointer to the value of the internal data
          double* const value_pt = internal_data_pt(i)->value_pt(j);
          
          //Save the old value of the Internal data
          const double old_var = *value_pt;
          
          //Increment the value of the Internal data
          *value_pt += fd_step;
          
          //Now update any slaved variables
          update_in_internal_fd(i);
          
          //Calculate the new residuals
          get_residuals(newres);
          
          //Do finite differences
          for(unsigned m=0;m<n_dof;m++)
           {
            double sum = (newres[m] - residuals[m])/fd_step;
            //Stick the entry into the Jacobian matrix
            jacobian(m,local_unknown) = sum;
           }
          
          //Reset the Internal data
          *value_pt = old_var;
          
          //Reset any slaved variables
          reset_in_internal_fd(i);
         }
       }
     } //End of finite-differencing for datum (if block)
   }

  //End of finite difference loop
  //Final reset of any slaved data
  reset_after_internal_fd();
 }

//============================================================================
/// This function calculates the entries of Jacobian matrix, used in 
/// the Newton method, associated with the external degrees of freedom.
/// It does this using finite differences, 
/// rather than an analytical formulation, so can be done in total generality.
/// 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 GeneralisedElement::
 fill_in_jacobian_from_external_by_fd(Vector<double> &residuals,
                                      DenseMatrix<double> &jacobian,
                                      const bool &fd_all_data)
 {
  //Locally cache the number of external data
  const unsigned n_external_data = Nexternal_data;
  //If there aren't any external data, then return straight away
  if(n_external_data == 0) {return;}

  //Call the update function to ensure that the element is in
  //a consistent state before finite differencing starts
  update_before_external_fd();

  //Find the number of dofs in the element
  const unsigned n_dof = ndof();

  //Create newres vector
  Vector<double> newres(n_dof);

  //Integer storage for local unknown
  int local_unknown=0;
  
  //Use the default finite difference step
  const double fd_step = Default_fd_jacobian_step;

  //Loop over the external data
  for(unsigned i=0;i<n_external_data;i++)
   {
    //If we are doing all finite differences or 
    //the variable is included in the finite difference loop, do it
    if(fd_all_data || external_data_fd(i))
     {
      //Get the number of value at the external data
      const unsigned n_value = external_data_pt(i)->nvalue();
   
      //Loop over the number of values
      for(unsigned j=0;j<n_value;j++)
       {
        //Get the local equation number
        local_unknown = external_local_eqn(i,j);
        //If it's not pinned
        if(local_unknown >= 0)
         {
          //Get a pointer to the External data value
          double* const value_pt = external_data_pt(i)->value_pt(j);
          
          //Save the old value of the External data
          const double old_var = *value_pt;
          
          //Increment the value of the External data
          *value_pt += fd_step;
          
          //Now update any slaved variables
          update_in_external_fd(i);
          
          //Calculate the new residuals
          get_residuals(newres);
          
          //Do finite differences
          for(unsigned m=0;m<n_dof;m++)
           {
            double sum = (newres[m] - residuals[m])/fd_step;
            //Stick the entry into the Jacobian matrix
            jacobian(m,local_unknown) = sum;
           }
          
          //Reset the External data
          *value_pt = old_var;
          
          //Reset any slaved variables
          reset_in_external_fd(i);
         }
       }
     } //End of finite differencing for datum (if block)
   }

  //End of finite difference loop
  //Final reset of any slaved data
  reset_after_external_fd();
 }

 //=====================================================================
 /// \short Add the elemental contribution to the mass matrix
 /// and the residuals vector. Note that
 /// this function will 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.
 //=====================================================================
 void GeneralisedElement:: 
 fill_in_contribution_to_mass_matrix(Vector<double> &residuals,
                                     DenseMatrix<double> &mass_matrix) 
 {
  std::string error_message =
   "Empty fill_in_contribution_to_mass_matrix() has been called.\n";
  error_message +=
   "This function is called from the default implementation of\n";
  error_message += "get_mass_matrix();\n";
  error_message += 
    "and must calculate the residuals vector and mass matrix ";
  error_message += "without initialising any of their entries.\n\n";
  
  error_message += 
   "If you do not wish to use these defaults, you must overload\n";
  error_message += 
   "get_mass_matrix(), which must initialise the entries\n";
  error_message +=
   "of the residuals vector and mass matrix to zero.\n";
  
  throw 
   OomphLibError(error_message,
                 "GeneralisedElement::fill_in_contribution_to_mass_matrix()",
                 OOMPH_EXCEPTION_LOCATION);
 }

 //=====================================================================
 /// \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. The default is deliberately
 /// broken.
 //=====================================================================
 void GeneralisedElement::fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals,
  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix)
 {
  std::string error_message =
   "Empty fill_in_contribution_to_jacobian_and_mass_matrix() has been ";
  error_message += "called.\n";
  error_message +=
   "This function is called from the default implementation of\n";
  error_message += "get_jacobian_and_mass_matrix();\n";
  error_message += 
    "and must calculate the residuals vector and mass and jacobian matrices ";
  error_message += "without initialising any of their entries.\n\n";
  
  error_message += 
   "If you do not wish to use these defaults, you must overload\n";
  error_message += 
   "get_jacobian_and_mass_matrix(), which must initialise the entries\n";
  error_message +=
   "of the residuals vector, jacobian and mass matrix to zero.\n";
  
  throw 
   OomphLibError(
    error_message,
    "GeneralisedElement::fill_in_contribution_to_jacobian_and_mass_matrix()",
    OOMPH_EXCEPTION_LOCATION);
 }


 //=====================================================================
 /// 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 deliberately broken
 //========================================================================
 void GeneralisedElement::fill_in_contribution_to_dresiduals_dparameter(
  double* const &parameter_pt, Vector<double> &dres_dparam)
 {
  std::string error_message =
   "Empty fill_in_contribution_to_dresiduals_dparameter() has been ";
  error_message += "called.\n";
  error_message +=
   "This function is called from the default implementation of\n";
  error_message += "get_dresiduals_dparameter();\n";
  error_message += 
   "and must calculate the derivatives of the residuals vector with respect\n";
  error_message += "to the passed parameter ";
  error_message += "without initialising any values.\n\n";

  error_message += 
   "If you do not wish to use these defaults, you must overload\n";
  error_message += 
   "get_dresiduals_dparameter(), which must initialise the entries\n";
  error_message +=
   "of the returned vector to zero.\n";
  
  error_message += 
   "This function is intended for use instead of the default (global) \n";
  error_message += 
   "finite-difference routine when analytic expressions are to be used\n";
  error_message += "in continuation and bifurcation tracking problems.\n\n";
  error_message += "This function is only called when the function\n";
  error_message += 
   "Problem::set_analytic_dparameter() has been called in the driver code\n";
  
  throw 
   OomphLibError(
    error_message,
    "GeneralisedElement::fill_in_contribution_to_dresiduals_dparameter()",
    OOMPH_EXCEPTION_LOCATION);
 }

 //======================================================================
 /// 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.
 //========================================================================
 void GeneralisedElement::fill_in_contribution_to_djacobian_dparameter(
  double* const &parameter_pt,
  Vector<double> &dres_dparam, 
  DenseMatrix<double> &djac_dparam)
 {
  std::string error_message =
   "Empty fill_in_contribution_to_djacobian_dparameter() has been ";
  error_message += "called.\n";
  error_message +=
   "This function is called from the default implementation of\n";
  error_message += "get_djacobian_dparameter();\n";
  error_message += 
   "and must calculate the derivatives of residuals vector and jacobian ";
   error_message += "matrix\n";
  error_message += "with respect to the passed parameter ";
  error_message += "without initialising any values.\n\n";

  error_message += 
   "If you do not wish to use these defaults, you must overload\n";
  error_message += 
   "get_djacobian_dparameter(), which must initialise the entries\n";
  error_message +=
   "of the returned vector and matrix to zero.\n\n";
  
  error_message += 
   "This function is intended for use instead of the default (global) \n";
  error_message += 
   "finite-difference routine when analytic expressions are to be used\n";
  error_message += "in continuation and bifurcation tracking problems.\n\n";
  error_message += "This function is only called when the function\n";
  error_message += 
   "Problem::set_analytic_dparameter() has been called in the driver code\n";

  
  throw 
   OomphLibError(
    error_message,
    "GeneralisedElement::fill_in_contribution_to_djacobian_dparameter()",
    OOMPH_EXCEPTION_LOCATION);
 }


 //=====================================================================
 /// \short Add the elemental contribution to the derivative of the 
 /// jacobian matrix,
 /// mass matrix and the residuals vector with respect to a 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. The default is deliberately
 /// broken.
 //=====================================================================
 void GeneralisedElement::
 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)
 {
  std::string error_message =
   "Empty fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter() has";
  error_message += " been called.\n";
  error_message +=
   "This function is called from the default implementation of\n";
  error_message += "get_djacobian_and_dmass_matrix_dparameter();\n";
  error_message += 
    "and must calculate the residuals vector and mass and jacobian matrices ";
  error_message += "without initialising any of their entries.\n\n";
  
  error_message += 
   "If you do not wish to use these defaults, you must overload\n";
  error_message += 
   "get_djacobian_and_dmass_matrix_dparameter(), which must initialise the\n";
    error_message += 
   "entries of the returned vector and  matrices to zero.\n";


  error_message += 
   "This function is intended for use instead of the default (global) \n";
  error_message += 
   "finite-difference routine when analytic expressions are to be used\n";
  error_message += "in continuation and bifurcation tracking problems.\n\n";
  error_message += "This function is only called when the function\n";
  error_message += 
   "Problem::set_analytic_dparameter() has been called in the driver code\n";

  

  throw 
   OomphLibError(
    error_message,
    "GeneralisedElement::fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter()",
    OOMPH_EXCEPTION_LOCATION);
 }



 //========================================================================
 /// 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}
 /// At the moment the dof pointer is passed in to enable
 /// easy calculation of finite difference default
//==========================================================================
 void GeneralisedElement::fill_in_contribution_to_hessian_vector_products(
  Vector<double> const &Y,
  DenseMatrix<double> const &C,
  DenseMatrix<double> &product)
 {
  std::string error_message =
   "Empty fill_in_contribution_to_hessian_vector_products() has been ";
  error_message += "called.\n";
  error_message +=
   "This function is called from the default implementation of\n";
  error_message += "get_hessian_vector_products(); ";
  error_message += " and must calculate the products \n";
  error_message += "of the hessian matrix with the (global) ";
  error_message += "vectors Y and C\n";
  error_message += "without initialising any values.\n\n";

  error_message += 
   "If you do not wish to use these defaults, you must overload\n";
  error_message += 
   "get_hessian_vector_products(), which must initialise the entries\n";
  error_message +=
   "of the returned vector to zero.\n\n";
  
  error_message += 
   "This function is intended for use instead of the default (global) \n";
  error_message += 
   "finite-difference routine when analytic expressions are to be used.\n";
  error_message += "This function is only called when the function\n";
  error_message += 
   "Problem::set_analytic_hessian_products() has been called in the driver code\n";
  
  throw 
   OomphLibError(
    error_message,
    "GeneralisedElement::fill_in_contribution_to_hessian_vector_product()",
    OOMPH_EXCEPTION_LOCATION);
 }

 //========================================================================
 /// Fill in the contribution to the inner products between given
 /// pairs of history values
 //==========================================================================
 void GeneralisedElement::fill_in_contribution_to_inner_products(
  Vector<std::pair<unsigned,unsigned> > const &history_index,
  Vector<double> &inner_product)
 {
  std::string error_message =
   "Empty fill_in_contribution_to_inner_products() has been called.\n";
  error_message +=
   "This function is called from the default implementation of\n";
  error_message += "get_inner_products();\n ";
  error_message += "It must calculate the inner products over the element\n";
  error_message += "of given pairs of history values\n";
  error_message += "without initialision of the output vector.\n\n";

  error_message += 
   "If you do not wish to use these defaults, you must overload\n";
  error_message += 
   "get_inner_products(), which must initialise the entries\n";
  error_message +=
   "of the returned vector to zero.\n\n";
  
  throw 
   OomphLibError(
    error_message,
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
 }

 //========================================================================
 /// \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
 //==========================================================================
 void GeneralisedElement::fill_in_contribution_to_inner_product_vectors(
  Vector<unsigned>  const &history_index,
  Vector<Vector<double> > &inner_product_vector)
 {
  std::string error_message =
   "Empty fill_in_contribution_to_inner_product_vectors() has been called.\n";
  error_message +=
   "This function is called from the default implementation of\n";
  error_message += "get_inner_product_vectors(); ";
  error_message += " and must calculate the vectors \n";
  error_message += "corresponding to the input history values\n ";
  error_message += "that when multiplied by other vectors of history values\n";
  error_message += "return the appropriate dot products.\n\n";
  error_message += "The output vectors must not be initialised.\n\n";

  error_message += 
   "If you do not wish to use these defaults, you must overload\n";
  error_message += 
   "get_inner_products(), which must initialise the entries\n";
  error_message +=
   "of the returned vector to zero.\n\n";
    
  throw 
   OomphLibError(
    error_message,
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
 }





//==========================================================================
/// Self-test: Have all internal values been classified as 
/// pinned/unpinned? Return 0 if OK. 
//==========================================================================
 unsigned GeneralisedElement::self_test()
 {
  // Initialise
  bool passed=true;
 
  //Loop over internal Data
  for(unsigned i=0;i<Ninternal_data;i++)
   {
    if(internal_data_pt(i)->self_test()!=0)
     {
      passed=false;
      oomph_info 
       << "\n ERROR: Failed GeneralisedElement::self_test()!" << std::endl;
      oomph_info 
       << "for internal data object number: " << i << std::endl;
     }
   }

  //Loop over external Data
  for(unsigned i=0;i<Nexternal_data;i++)
   {
    if(external_data_pt(i)->self_test()!=0)
     {
      passed=false;
      oomph_info 
       << "\n ERROR: Failed GeneralisedElement::self_test()!" << std::endl;
      oomph_info 
       << "for external data object number: " << i << std::endl;
     }
   }

  // Return verdict
  if (passed) {return 0;}
  else {return 1;}
 }


//======================================================================
/// Helper namespace for tolerances, number of iterations, etc
/// used in the locate_zeta function in FiniteElement
//======================================================================
namespace Locate_zeta_helpers
{
 /// Convergence tolerance for the newton solver
 double Newton_tolerance = 1.0e-7;      
 
 /// Maximum number of newton iterations
 unsigned Max_newton_iterations = 10;
 
 /// \short Multiplier for (zeta-based) outer radius of element used for
 /// deciding that point is outside element. Set to negative value
 /// to suppress test.
 double Radius_multiplier_for_fast_exit_from_locate_zeta=1.5;

 /// 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 when locating coordinate)
 unsigned N_local_points = 5;
}


///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
//  Functions for finite elements
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////

//======================================================================
 /// 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 FiniteElement::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 node_pt(n)->x(i);
  }



//======================================================================
/// \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 FiniteElement::describe_local_dofs(std::ostream& out,
                                        const std::string& current_string) const
{
 // Call the standard finite element classification.
 GeneralisedElement::describe_local_dofs(out,current_string);
 describe_nodal_local_dofs(out,current_string);
}

//======================================================================
// \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 FiniteElement::
describe_nodal_local_dofs(std::ostream& out,
                          const std::string& current_string) const
{
 //Find the number of nodes
 const unsigned n_node = nnode();
 for(unsigned n=0;n<n_node;n++)
  {
   //Pointer to node
   Node* const nod_pt = node_pt(n);
      
   std::stringstream conversion;
   conversion<<" of Node "<<n<<current_string;
   std::string in(conversion.str());
   nod_pt->describe_dofs(out,in);
  }// End if for n_node
}// End describe_nodal_local_dofs

//========================================================================
/// \short Internal function used to check for singular or negative values
/// of the determinant of the Jacobian of the mapping between local and 
/// global or lagrangian coordinates. Negative jacobians are allowed if the 
/// Accept_negative_jacobian flag is set to true. 
//========================================================================
 void FiniteElement::check_jacobian(const double &jacobian) const
 {
  //First check for a zero jacobian
  if(std::fabs(jacobian) < Tolerance_for_singular_jacobian)
   {
    if (FiniteElement::Suppress_output_while_checking_for_inverted_elements)
     {
      throw OomphLibQuietException();
     }
    else
     {
      std::ostringstream error_stream;
      error_stream 
       << "Determinant of Jacobian matrix is zero --- "
       << "singular mapping!\nThe determinant of the "
       << "jacobian is " << std::fabs(jacobian)
       << " which is smaller than the treshold "
       << Tolerance_for_singular_jacobian <<"\n"
       << "You can change this treshold, by specifying "
       << "FiniteElement::Tolerance_for_singular_jacobian \n"
       << "Here are the nodal coordinates of the inverted element\n" 
       << "in the form \n\n       x,y[,z], hang_status\n\n" 
       << "where hang_status = 1 or 2 for non-hanging or hanging\n" 
       << "nodes respectively (useful for automatic sizing of\n" 
       << "tecplot markers to identify the hanging nodes). \n\n" ;
      unsigned n_dim_nod=node_pt(0)->ndim();
      unsigned n_nod=nnode();
      unsigned hang_count=0;
      for (unsigned j=0;j<n_nod;j++)
       {
        for (unsigned i=0;i<n_dim_nod;i++)
         {
          error_stream << node_pt(j)->x(i) << " ";
         }
        if (node_pt(j)->is_hanging())
         {
          error_stream << " 2";
          hang_count++;
         }
        else
         {
          error_stream << " 1";
         }
        error_stream << std::endl;
       }
      error_stream << std::endl << std::endl;
      if ((Macro_elem_pt!=0)&&(0!=hang_count))
       {
        error_stream 
         << "NOTE:  Offending element is associated with a MacroElement\n"
         << "       AND the element has hanging nodes! \n"
         << "       If an element is thin and highly curved, the \n"
         << "       constraints imposed by\n \n"
         << "       (1) inter-element continuity (imposed by the hanging\n"
         << "           node constraints) and \n\n"
         << "       (2) the need to respect curvilinear domain boundaries\n"
         << "           during mesh refinement (imposed by the element's\n"
         << "           macro element mapping)\n\n"
         << "       may be irreconcilable!  \n \n"
         << "       You may have to re-design your base mesh to avoid \n"
         << "       the creation of thin, highly curved elements during\n"
         << "       the refinement process.\n"
         << std::endl;
       }
      throw OomphLibError(error_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }


  //Now check for negative jacobians, if we're not allowing them (default)
  if((Accept_negative_jacobian==false) && (jacobian < 0.0))
   {

    if (FiniteElement::Suppress_output_while_checking_for_inverted_elements)
     {
      throw OomphLibQuietException();
     }
    else
     {
      std::ostringstream error_stream;
      error_stream 
       << "Negative Jacobian in transform from "
       << "local to global coordinates"
       << std::endl 
       << "       You have an inverted coordinate system!" 
       << std::endl << std::endl;
      error_stream 
       << "Here are the nodal coordinates of the inverted element\n" 
       << "in the form \n\n       x,y[,z], hang_status\n\n" 
       << "where hang_status = 1 or 2 for non-hanging or hanging\n" 
       << "nodes respectively (useful for automatic sizing of\n" 
       << "tecplot markers to identify the hanging nodes). \n\n" ;
      unsigned n_dim_nod=node_pt(0)->ndim();
      unsigned n_nod=nnode();
      unsigned hang_count=0;
      for (unsigned j=0;j<n_nod;j++)
       {
        for (unsigned i=0;i<n_dim_nod;i++)
         {
          error_stream << node_pt(j)->x(i) << " ";
         }
        if (node_pt(j)->is_hanging())
         {
          error_stream << " 2";
          hang_count++;
         }
        else
         {
          error_stream << " 1";
         }
        error_stream << std::endl;
       }
      error_stream << std::endl << std::endl;
      if ((Macro_elem_pt!=0)&&(0!=hang_count))
       {
        error_stream 
         << "NOTE: The inverted element is associated with a MacroElement\n"
         << "       AND the element has hanging nodes! \n"
         << "       If an element is thin and highly curved, the \n"
         << "       constraints imposed by\n \n"
         << "       (1) inter-element continuity (imposed by the hanging\n"
         << "           node constraints) and \n\n"
         << "       (2) the need to respect curvilinear domain boundaries\n"
         << "           during mesh refinement (imposed by the element's\n"
         << "           macro element mapping)\n\n"
         << "       may be irreconcilable!  \n \n"
         << "       You may have to re-design your base mesh to avoid \n"
         << "       the creation of thin, highly curved elements during\n"
         << "       the refinement process.\n"
         << std::endl;
       }
      
      error_stream 
       << std::endl << std::endl
       << "If you believe that inverted elements do not cause any\n" 
       << "problems in your specific application you can \n " 
       << "suppress this test by: " << std::endl
       << "  i) setting the (static) flag " 
       << "FiniteElement::Accept_negative_jacobian to be true" << std::endl;
      error_stream << " ii) switching OFF the PARANOID flag" 
                   << std::endl << std::endl;
      
      throw OomphLibError(error_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
 }

//=========================================================================
/// Internal function that is used to assemble the jacobian of the mapping
/// from local coordinates (s) to the eulerian coordinates (x), given the
/// derivatives of the shape functions. The entire jacobian matrix is 
/// constructed and this function will only work if there are the same number
/// of local coordinates as global coordinates (i.e. for "bulk" elements).
//=========================================================================
void FiniteElement::
assemble_local_to_eulerian_jacobian(const DShape &dpsids,
                                    DenseMatrix<double> &jacobian) const
{
 //Locally cache the elemental dimension
 const unsigned el_dim = dim();
 //The number of shape functions must be equal to the number
 //of nodes (by definition)
 const unsigned n_shape = nnode();
 //The number of shape function types must be equal to the number
 //of nodal position types (by definition)
 const unsigned n_shape_type = nnodal_position_type();

#ifdef PARANOID
 //Check for dimensional compatibility
 if(Elemental_dimension != Nodal_dimension)
  {
   std::ostringstream error_message;
   error_message << "Dimension mismatch" << std::endl;
   error_message << "The elemental dimension: " << Elemental_dimension
                 << " must equal the nodal dimension: " 
                 << Nodal_dimension 
                 << " for the jacobian of the mapping to be well-defined"
                 << std::endl;
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
   

 //Loop over the rows of the jacobian
 for(unsigned i=0;i<el_dim;i++)
  {
   //Loop over the columns of the jacobian
   for(unsigned j=0;j<el_dim;j++)
    {
     //Zero the entry
     jacobian(i,j) = 0.0;
     //Loop over the shape functions
     for(unsigned l=0;l<n_shape;l++)
      {
       for(unsigned k=0;k<n_shape_type;k++)
        {
         //Jacobian is dx_j/ds_i, which is represented by the sum
         //over the dpsi/ds_i of the nodal points X j
         //Call the Non-hanging version of positions
         //This is overloaded in refineable elements
         jacobian(i,j) += raw_nodal_position_gen(l,k,j)*dpsids(l,k,i);
        }
      }
    }
  }
}

//=========================================================================
/// Internal function that is used to assemble the jacobian of second
/// derivatives of the mapping from local coordinates (s) to the 
/// eulerian coordinates (x), given the second derivatives of the 
/// shape functions. 
//=========================================================================
 void FiniteElement::
 assemble_local_to_eulerian_jacobian2(const DShape &d2psids,
                                      DenseMatrix<double> &jacobian2) const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  //Must be equal to the number of nodes and their position types
  //by the definition of the shape function.
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_position_type();
  //Find the number of second derivatives
  const unsigned n_row = N2deriv[el_dim];
 
  //Assemble the "jacobian" (d^2 x_j/ds_i^2) for second derivatives of 
  //shape functions
  //Loop over the rows (number of second derivatives)
  for(unsigned i=0;i<n_row;i++)
   {
    //Loop over the columns (element dimension
    for(unsigned j=0;j<el_dim;j++)
     {
      //Zero the entry
      jacobian2(i,j) = 0.0;
      //Loop over the shape functions
      for(unsigned l=0;l<n_shape;l++)
       {
        //Loop over the shape function types
        for(unsigned k=0;k<n_shape_type;k++)
         {
          //Add the terms to the jacobian entry
          //Call the Non-hanging version of positions
          //This is overloaded in refineable elements
          jacobian2(i,j) += raw_nodal_position_gen(l,k,j)*d2psids(l,k,i);
         }
       }
     }
   }
 }

 //=====================================================================
 /// Assemble the covariant Eulerian base vectors and return them in the
 /// matrix interpolated_G. The derivatives of the shape functions with
 /// respect to the local coordinate should already have been calculated
 /// before calling this function
 //=====================================================================
 void FiniteElement::assemble_eulerian_base_vectors(
  const DShape &dpsids, DenseMatrix<double> &interpolated_G) const
 {
  //Find the number of nodes and position types
  const unsigned n_node = nnode();
  const unsigned n_position_type = nnodal_position_type();
  //Find the dimension of the node and element
  const unsigned n_dim_node = nodal_dimension();
  const unsigned n_dim_element = dim();

  //Loop over the dimensions and compute the entries of the 
  //base vector matrix
  for(unsigned i=0;i<n_dim_element;i++)
   {
    for(unsigned j=0;j<n_dim_node;j++)
     {
      //Initialise the j-th component of the i-th base vector to zero
      interpolated_G(i,j) = 0.0;
      for(unsigned l=0;l<n_node;l++)
       {
        for(unsigned k=0;k<n_position_type;k++)
         {
          interpolated_G(i,j) += raw_nodal_position_gen(l,k,j)*dpsids(l,k,i);
         }
       }
     }
   }
 }


//============================================================================
/// Zero-d specialisation of function to calculate inverse of jacobian mapping
//============================================================================
 template<>
 double FiniteElement::
 invert_jacobian<0>(const DenseMatrix<double> &jacobian,
                    DenseMatrix<double> &inverse_jacobian) const
 {
  //Issue a warning
  oomph_info << "\nWarning: You are trying to invert the jacobian for "
             << "a 'point' element" << std::endl
             << "This makes no sense and is almost certainly an error"
             << std::endl << std::endl;

  //Dummy return
  return(1.0);
 }


//===========================================================================
/// One-d specialisation of function to calculate inverse of jacobian mapping
//===========================================================================
 template<>
 double FiniteElement::
 invert_jacobian<1>(const DenseMatrix<double> &jacobian,
                    DenseMatrix<double> &inverse_jacobian) const
 {
  //Calculate the determinant of the matrix
  const double det = jacobian(0,0);

//Report if Matrix is singular or negative
#ifdef PARANOID
  check_jacobian(det);
#endif 

  //Calculate the inverse --- trivial in 1D
  inverse_jacobian(0,0) = 1.0/jacobian(0,0);
 
  //Return the determinant
  return(det);
 }

//===========================================================================
/// Two-d specialisation of function to calculate inverse of jacobian mapping
//===========================================================================
 template<>
 double FiniteElement::
 invert_jacobian<2>(const DenseMatrix<double> &jacobian,
                    DenseMatrix<double> &inverse_jacobian) const
 {
  //Calculate the determinant of the matrix
  const double det = jacobian(0,0)*jacobian(1,1) - jacobian(0,1)*jacobian(1,0);

//Report if Matrix is singular or negative
#ifdef PARANOID
  check_jacobian(det);
#endif 

  //Calculate the inverset of the 2x2 matrix
  inverse_jacobian(0,0) =  jacobian(1,1)/det;
  inverse_jacobian(0,1) = -jacobian(0,1)/det;
  inverse_jacobian(1,0) = -jacobian(1,0)/det;
  inverse_jacobian(1,1) =  jacobian(0,0)/det;

  //Return the jacobian
  return(det);
 }

//=============================================================================
/// Three-d specialisation of function to calculate inverse of jacobian mapping
//=============================================================================
 template<>
 double FiniteElement::
 invert_jacobian<3>(const DenseMatrix<double> &jacobian,
                    DenseMatrix<double> &inverse_jacobian) const
 {
  //Calculate the determinant of the matrix
  const double det = jacobian(0,0)*jacobian(1,1)*jacobian(2,2) 
   + jacobian(0,1)*jacobian(1,2)*jacobian(2,0) 
   + jacobian(0,2)*jacobian(1,0)*jacobian(2,1) 
   - jacobian(0,0)*jacobian(1,2)*jacobian(2,1) 
   - jacobian(0,1)*jacobian(1,0)*jacobian(2,2) 
   - jacobian(0,2)*jacobian(1,1)*jacobian(2,0); 
 
  //Report if Matrix is singular or negative
#ifdef PARANOID
  check_jacobian(det);
#endif

  //Calculate the inverse of the 3x3 matrix
  inverse_jacobian(0,0) =  (jacobian(1,1)*jacobian(2,2) 
                            - jacobian(1,2)*jacobian(2,1))/det; 
  inverse_jacobian(0,1) = -(jacobian(0,1)*jacobian(2,2) 
                            - jacobian(0,2)*jacobian(2,1))/det; 
  inverse_jacobian(0,2) =  (jacobian(0,1)*jacobian(1,2) 
                            - jacobian(0,2)*jacobian(1,1))/det; 
  inverse_jacobian(1,0) = -(jacobian(1,0)*jacobian(2,2) 
                            - jacobian(1,2)*jacobian(2,0))/det; 
  inverse_jacobian(1,1) =  (jacobian(0,0)*jacobian(2,2) 
                            - jacobian(0,2)*jacobian(2,0))/det; 
  inverse_jacobian(1,2) = -(jacobian(0,0)*jacobian(1,2) 
                            - jacobian(0,2)*jacobian(1,0))/det; 
  inverse_jacobian(2,0) =  (jacobian(1,0)*jacobian(2,1) 
                            - jacobian(1,1)*jacobian(2,0))/det; 
  inverse_jacobian(2,1) = -(jacobian(0,0)*jacobian(2,1) 
                            - jacobian(0,1)*jacobian(2,0))/det; 
  inverse_jacobian(2,2) =  (jacobian(0,0)*jacobian(1,1) 
                            - jacobian(0,1)*jacobian(1,0))/det; 

  //Return the determinant
  return(det);
 }

//========================================================================
/// Template-free interface for inversion of the jacobian of a mapping.
/// This is slightly inefficient, given that it uses a switch statement.
/// It can always be overloaded in specific geometric elements, 
/// for efficiency reasons.
//========================================================================
 double FiniteElement::
 invert_jacobian_mapping(const DenseMatrix<double> &jacobian,
                         DenseMatrix<double> &inverse_jacobian) const
 {
  //Find the spatial dimension of the element
  const unsigned el_dim = dim();
  //Call the appropriate templated function, depending on the 
  //element dimension
  switch(el_dim)
   {
   case 0:
    return invert_jacobian<0>(jacobian,inverse_jacobian);
    break;
   case 1:
    return invert_jacobian<1>(jacobian,inverse_jacobian);
    break;
   case 2:
    return invert_jacobian<2>(jacobian,inverse_jacobian);
    break;
   case 3:
    return invert_jacobian<3>(jacobian,inverse_jacobian);
    break;
    //Catch-all default case: issue warning and die
   default:
    std::ostringstream error_stream;
    error_stream 
     << "Dimension of the element must be 0,1,2 or 3, not " << el_dim 
     << std::endl;
  
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  //Dummy return for Intel compiler
  return 1.0;
 }

//============================================================================
/// Zero-d specialisation of function to calculate the derivative of the
/// jacobian of a mapping with respect to the nodal coordinates X_ij.
//============================================================================
 template<>
 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<0>(
  const DenseMatrix<double> &jacobian,const DShape &dpsids,
  DenseMatrix<double> &djacobian_dX) const
 {
  // Issue a warning
  oomph_info << "\nWarning: You are trying to calculate derivatives of "
             << "a jacobian w.r.t. nodal coordinates for a 'point' "
             << "element." << std::endl
             << "This makes no sense and is almost certainly an error."
             << std::endl << std::endl;
 }

//===========================================================================
/// One-d specialisation of function to calculate the derivative of the
/// jacobian of a mapping with respect to the nodal coordinates X_ij.
//===========================================================================
 template<>
 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<1>(
  const DenseMatrix<double> &jacobian,const DShape &dpsids,
  DenseMatrix<double> &djacobian_dX) const
 {
  // Determine the number of nodes in the element
  const unsigned n_node = nnode();
  
  // Loop over nodes
  for(unsigned j=0;j<n_node;j++)
   {
    djacobian_dX(0,j) = dpsids(j,0);
   }
 }

//===========================================================================
/// Two-d specialisation of function to calculate the derivative of the
/// jacobian of a mapping with respect to the nodal coordinates X_ij.
//===========================================================================
 template<>
 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<2>(
  const DenseMatrix<double> &jacobian,const DShape &dpsids,
  DenseMatrix<double> &djacobian_dX) const
 {
  // Determine the number of nodes in the element
  const unsigned n_node = nnode();

  // Loop over nodes
  for(unsigned j=0;j<n_node;j++)
   {
    // i=0
    djacobian_dX(0,j) = dpsids(j,0)*jacobian(1,1) - dpsids(j,1)*jacobian(0,1);

    // i=1
    djacobian_dX(1,j) = dpsids(j,1)*jacobian(0,0) - dpsids(j,0)*jacobian(1,0);
   }
 }

//=============================================================================
/// Three-d specialisation of function to calculate the derivative of the
/// jacobian of a mapping with respect to the nodal coordinates X_ij.
//=============================================================================
 template<>
 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<3>(
  const DenseMatrix<double> &jacobian,const DShape &dpsids,
  DenseMatrix<double> &djacobian_dX) const
 {
  // Determine the number of nodes in the element
  const unsigned n_node = nnode();

  // Loop over nodes
  for(unsigned j=0;j<n_node;j++)
   {
    // i=0
    djacobian_dX(0,j)
     = dpsids(j,0)*(jacobian(1,1)*jacobian(2,2)
                    - jacobian(1,2)*jacobian(2,1))
     + dpsids(j,1)*(jacobian(0,2)*jacobian(2,1)
                    - jacobian(0,1)*jacobian(2,2))
     + dpsids(j,2)*(jacobian(0,1)*jacobian(1,2)
                    - jacobian(0,2)*jacobian(1,1));
    
    // i=1
    djacobian_dX(1,j)
     = dpsids(j,0)*(jacobian(1,2)*jacobian(2,0)
                    - jacobian(1,0)*jacobian(2,2))
     + dpsids(j,1)*(jacobian(0,0)*jacobian(2,2)
                    - jacobian(0,2)*jacobian(2,0))
     + dpsids(j,2)*(jacobian(0,2)*jacobian(1,0)
                    - jacobian(0,0)*jacobian(1,2));

    // i=2
    djacobian_dX(2,j)
     = dpsids(j,0)*(jacobian(1,0)*jacobian(2,1)
                    - jacobian(1,1)*jacobian(2,0))
     + dpsids(j,1)*(jacobian(0,1)*jacobian(2,0)
                    - jacobian(0,0)*jacobian(2,1))
     + dpsids(j,2)*(jacobian(0,0)*jacobian(1,1)
                    - jacobian(0,1)*jacobian(1,0));
   }
 }

//============================================================================
/// Zero-d specialisation of function to 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$.
//============================================================================
 template<>
 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<0>(
  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
 {
  // Issue a warning
  oomph_info << "\nWarning: You are trying to calculate derivatives of "
             << "eulerian derivatives of shape functions w.r.t. nodal "
             << "coordinates for a 'point' element." << std::endl
             << "This makes no sense and is almost certainly an error."
             << std::endl << std::endl;
 }

//===========================================================================
/// One-d specialisation of function to 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$.
//===========================================================================
 template<>
 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<1>(
  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
 {
  // Find inverse of determinant of jacobian of mapping
  const double inv_det_jac = 1.0/det_jacobian;

  // Determine the number of nodes in the element
  const unsigned n_node = nnode();

  // Loop over the shape functions
  for(unsigned q=0;q<n_node;q++)
   {
    // Loop over the shape functions
    for(unsigned j=0;j<n_node;j++)
     {
      d_dpsidx_dX(0,q,j,0)
       = - djacobian_dX(0,q)*dpsids(j,0)*inv_det_jac*inv_det_jac;
     }
   }
 }

//===========================================================================
/// Two-d specialisation of function to 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$.
//===========================================================================
 template<>
 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<2>(
  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
 {
  // Find inverse of determinant of jacobian of mapping
  const double inv_det_jac = 1.0/det_jacobian;

  // Determine the number of nodes in the element
  const unsigned n_node = nnode();

  // Loop over the spatial dimension (this must be 2)
  for(unsigned p=0;p<2;p++)
   {
    // Loop over the shape functions
    for(unsigned q=0;q<n_node;q++)
     {
      // Loop over the shape functions
      for(unsigned j=0;j<n_node;j++)
       {
        // i=0
        d_dpsidx_dX(p,q,j,0) = - djacobian_dX(p,q)*
         (inverse_jacobian(0,0)*dpsids(j,0)
          + inverse_jacobian(0,1)*dpsids(j,1));
        
        if(p==1)
         {
          d_dpsidx_dX(p,q,j,0)
           += dpsids(j,0)*dpsids(q,1) - dpsids(j,1)*dpsids(q,0);
         }
        d_dpsidx_dX(p,q,j,0) *= inv_det_jac;

        // i=1
        d_dpsidx_dX(p,q,j,1) = - djacobian_dX(p,q)*
         (inverse_jacobian(1,1)*dpsids(j,1)
          + inverse_jacobian(1,0)*dpsids(j,0));
        
        if(p==0)
         {
          d_dpsidx_dX(p,q,j,1)
           += dpsids(j,1)*dpsids(q,0) - dpsids(j,0)*dpsids(q,1);
         }
        d_dpsidx_dX(p,q,j,1) *= inv_det_jac;
       }
     }
   }
 }

//=============================================================================
/// Three-d specialisation of function to 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$.
//=============================================================================
 template<>
 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<3>(
  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
 {
  // Find inverse of determinant of jacobian of mapping
  const double inv_det_jac = 1.0/det_jacobian;
  
  // Determine the number of nodes in the element
  const unsigned n_node = nnode();
  
  // Loop over the spatial dimension (this must be 3)
  for(unsigned p=0;p<3;p++)
   {
    // Loop over the shape functions
    for(unsigned q=0;q<n_node;q++)
     {
      // Loop over the shape functions
      for(unsigned j=0;j<n_node;j++)
       {
        // Terms not multiplied by delta function
        for(unsigned i=0;i<3;i++)
         {
          d_dpsidx_dX(p,q,j,i)
           = - djacobian_dX(p,q)*(inverse_jacobian(i,0)*dpsids(j,0)
                                  + inverse_jacobian(i,1)*dpsids(j,1)
                                  + inverse_jacobian(i,2)*dpsids(j,2));
         }

        // Delta function terms
        switch(p)
         {
         case 0:
          d_dpsidx_dX(p,q,j,1)+=((dpsids(q,2)*jacobian(1,2)
                                  - dpsids(q,1)*jacobian(2,2))*dpsids(j,0)
                                 + (dpsids(q,0)*jacobian(2,2)
                                    - dpsids(q,2)*jacobian(0,2))*dpsids(j,1)
                                 + (dpsids(q,1)*jacobian(0,2)
                                    - dpsids(q,0)*jacobian(1,2))*dpsids(j,2));

          d_dpsidx_dX(p,q,j,2)+=((dpsids(q,1)*jacobian(2,1)
                                  - dpsids(q,2)*jacobian(1,1))*dpsids(j,0)
                                 + (dpsids(q,2)*jacobian(0,1)
                                    - dpsids(q,0)*jacobian(2,1))*dpsids(j,1)
                                 + (dpsids(q,0)*jacobian(1,1)
                                    - dpsids(q,1)*jacobian(0,1))*dpsids(j,2));
          break;

         case 1:
        
          d_dpsidx_dX(p,q,j,0)+=((dpsids(q,1)*jacobian(2,2)
                                  - dpsids(q,2)*jacobian(1,2))*dpsids(j,0)
                                 + (dpsids(q,2)*jacobian(0,2)
                                    - dpsids(q,0)*jacobian(2,2))*dpsids(j,1)
                                 + (dpsids(q,0)*jacobian(1,2)
                                    - dpsids(q,1)*jacobian(0,2))*dpsids(j,2));

          d_dpsidx_dX(p,q,j,2)+=((dpsids(q,2)*jacobian(1,0)
                                  - dpsids(q,1)*jacobian(2,0))*dpsids(j,0)
                                 + (dpsids(q,0)*jacobian(2,0)
                                    - dpsids(q,2)*jacobian(0,0))*dpsids(j,1)
                                 + (dpsids(q,1)*jacobian(0,0)
                                    - dpsids(q,0)*jacobian(1,0))*dpsids(j,2));
          break;
         
         case 2:

          d_dpsidx_dX(p,q,j,0)+=((dpsids(q,2)*jacobian(1,1)
                                  - dpsids(q,1)*jacobian(2,1))*dpsids(j,0)
                                 + (dpsids(q,0)*jacobian(2,1)
                                    - dpsids(q,2)*jacobian(0,1))*dpsids(j,1)
                                 + (dpsids(q,1)*jacobian(0,1)
                                    - dpsids(q,0)*jacobian(1,1))*dpsids(j,2));

          d_dpsidx_dX(p,q,j,1)+=((dpsids(q,1)*jacobian(2,0)
                                  - dpsids(q,2)*jacobian(1,0))*dpsids(j,0)
                                 + (dpsids(q,2)*jacobian(0,0)
                                    - dpsids(q,0)*jacobian(2,0))*dpsids(j,1)
                                 + (dpsids(q,0)*jacobian(1,0)
                                    - dpsids(q,1)*jacobian(0,0))*dpsids(j,2));
          break;
         }
        
        // Divide through by the determinant of the Jacobian mapping
        for(unsigned i=0;i<3;i++)
         {
          d_dpsidx_dX(p,q,j,i) *= inv_det_jac;
         }
       }
     }
   }
 }

//=======================================================================
/// Default value for the number of values at a node
//=======================================================================
 const unsigned FiniteElement::Default_Initial_Nvalue = 0;

//======================================================================
/// \short Default value that is used for the tolerance required when 
/// locating nodes via local coordinates
 const double FiniteElement::Node_location_tolerance = 1.0e-14;

//=======================================================================
/// Set the default tolerance for a singular jacobian
//=======================================================================
 double FiniteElement::Tolerance_for_singular_jacobian = 1.0e-16;

//======================================================================
/// Set the default value of the Accept_negative_jacobian flag to be 
/// false
//======================================================================
 bool FiniteElement::Accept_negative_jacobian = false;


//======================================================================
///  Set default for static boolean to suppress output while checking 
/// for inverted elements
//======================================================================
bool FiniteElement::Suppress_output_while_checking_for_inverted_elements=false;

//========================================================================
/// Static array that holds the number of rows in the second derivative
/// matrix as a function of spatial dimension. In one-dimension, there is
/// only one possible second derivative. In two-dimensions, there are three,
/// the two second derivatives and the mixed derivatives. In three 
/// dimensions there are six.
//=========================================================================
 const unsigned FiniteElement::N2deriv[4]={0,1,3,6};

//==========================================================================
/// Calculate the mapping from local to eulerian 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.
/// The local derivatives are passed as dpsids and the jacobian and 
/// inverse jacobian are returned.
//==========================================================================
 double FiniteElement::
 local_to_eulerian_mapping_diagonal(const DShape &dpsids,
                                    DenseMatrix<double> &jacobian,
                                    DenseMatrix<double> &inverse_jacobian)
  const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  //Equal to the number of nodes and their position types by definition
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_position_type();

#ifdef PARANOID
  //Check for dimension compatibility
 if(Elemental_dimension != Nodal_dimension)
  {
   std::ostringstream error_message;
   error_message << "Dimension mismatch" << std::endl;
   error_message << "The elemental dimension: " << Elemental_dimension
                 << " must equal the nodal dimension: " 
                 << Nodal_dimension 
                 << " for the jacobian of the mapping to be well-defined"
                 << std::endl;
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
  //In this case we assume that there are no cross terms, that is
  //global coordinate 0 is always in the direction of local coordinate 0
 
  //Loop over the coordinates
  for(unsigned i=0;i<el_dim;i++)
   {
    //Zero the jacobian and inverse jacobian entries
    for(unsigned j=0;j<el_dim;j++) 
     {jacobian(i,j) = 0.0; inverse_jacobian(i,j) = 0.0;}
   
    //Loop over the shape functions
    for(unsigned l=0;l<n_shape;l++)
     {
      //Loop over the types of dof
      for(unsigned k=0;k<n_shape_type;k++)
       {
        //Derivatives are always dx_{i}/ds_{i}
        jacobian(i,i) += raw_nodal_position_gen(l,k,i)*dpsids(l,k,i);
       }
     }
   }
 
  //Now calculate the determinant of the matrix
  double det = 1.0;
  for(unsigned i=0;i<el_dim;i++) {det *= jacobian(i,i);}
 
//Report if Matrix is singular, or negative
#ifdef PARANOID
  check_jacobian(det);
#endif
 
  //Calculate the inverse mapping (trivial in this case)
  for(unsigned i=0;i<el_dim;i++) 
   {inverse_jacobian(i,i) = 1.0/jacobian(i,i);}
 
  //Return the value of the Jacobian
  return(det);
 }

//========================================================================
/// Template-free interface calculating the derivative of the jacobian
/// of a mapping with respect to the nodal coordinates X_ij. This is
/// slightly inefficient, given that it uses a switch statement. It can
/// always be overloaded in specific geometric elements, for efficiency
/// reasons.
//========================================================================
void FiniteElement::dJ_eulerian_dnodal_coordinates(
 const DenseMatrix<double> &jacobian,const DShape &dpsids,
 DenseMatrix<double> &djacobian_dX) const
{
 // Determine the spatial dimension of the element
 const unsigned el_dim = dim();
 
#ifdef PARANOID
 // Determine the number of nodes in the element
 const unsigned n_node = nnode();
 
 // Check that djacobian_dX has the correct number of rows (= el_dim)
 if(djacobian_dX.nrow()!=el_dim)
  {
   std::ostringstream error_message;
   error_message << "djacobian_dX must have the same number of rows as the"
                 << "\nspatial dimension of the element.";
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
 // Check that djacobian_dX has the correct number of columns (= n_node)
 if(djacobian_dX.ncol()!=n_node)
  {
   std::ostringstream error_message;
   error_message << "djacobian_dX must have the same number of columns as the"
                 << "\nnumber of nodes in the element.";
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Call the appropriate templated function, depending on the 
 // element dimension
 switch(el_dim)
  {
  case 0:
   dJ_eulerian_dnodal_coordinates_templated_helper<0>(jacobian,dpsids,
                                                      djacobian_dX);
   break;
  case 1:
   dJ_eulerian_dnodal_coordinates_templated_helper<1>(jacobian,dpsids,
                                                      djacobian_dX);
   break;
  case 2:
   dJ_eulerian_dnodal_coordinates_templated_helper<2>(jacobian,dpsids,
                                                      djacobian_dX);
   break;
  case 3:
   dJ_eulerian_dnodal_coordinates_templated_helper<3>(jacobian,dpsids,
                                                      djacobian_dX);
   break;
   // Catch-all default case: issue warning and die
  default:
   std::ostringstream error_stream;
   error_stream 
    << "Dimension of the element must be 0,1,2 or 3, not " << el_dim 
    << std::endl;
   
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
}

//========================================================================
/// Template-free interface calculating 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$ 
/// This function is slightly inefficient, given that it uses a switch
/// statement. It can always be overloaded in specific geometric elements,
/// for efficiency reasons.
//========================================================================
void FiniteElement::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
{
 // Determine the spatial dimension of the element
 const unsigned el_dim = dim();

#ifdef PARANOID
 // Determine the number of nodes in the element
 const unsigned n_node = nnode();
 
 // Check that d_dpsidx_dX is of the correct size
 if(d_dpsidx_dX.nindex1()!=el_dim || d_dpsidx_dX.nindex2()!=n_node
    || d_dpsidx_dX.nindex3()!=n_node || d_dpsidx_dX.nindex4()!=el_dim)
  {
   std::ostringstream error_message;
   error_message << "d_dpsidx_dX must be of the following dimensions:"
                 << "\nd_dpsidx_dX(el_dim,n_node,n_node,el_dim)";
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Call the appropriate templated function, depending on the 
 // element dimension
 switch(el_dim)
  {
  case 0:
   d_dshape_eulerian_dnodal_coordinates_templated_helper<0>(
    det_jacobian,jacobian,djacobian_dX,inverse_jacobian,dpsids,d_dpsidx_dX);
   break;
  case 1:
   d_dshape_eulerian_dnodal_coordinates_templated_helper<1>(
    det_jacobian,jacobian,djacobian_dX,inverse_jacobian,dpsids,d_dpsidx_dX);
   break;
  case 2:
   d_dshape_eulerian_dnodal_coordinates_templated_helper<2>(
    det_jacobian,jacobian,djacobian_dX,inverse_jacobian,dpsids,d_dpsidx_dX);
   break;
  case 3:
   d_dshape_eulerian_dnodal_coordinates_templated_helper<3>(
    det_jacobian,jacobian,djacobian_dX,inverse_jacobian,dpsids,d_dpsidx_dX);
   break;
   // Catch-all default case: issue warning and die
  default:
   std::ostringstream error_stream;
   error_stream 
    << "Dimension of the element must be 0,1,2 or 3, not " << el_dim 
    << std::endl;
   
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
}

//=====================================================================
/// Convert derivatives w.r.t local coordinates to derivatives w.r.t the
/// coordinates used to assemble the inverse jacobian mapping passed as
/// inverse_jacobian. The derivatives passed in dbasis will be
/// modified in this function from dbasisds to dbasisdX.
//======================================================================
 void FiniteElement::transform_derivatives(const DenseMatrix<double> 
                                           &inverse_jacobian,
                                           DShape &dbasis) const
 {
  //Find the number of basis functions and basis function types
  const unsigned n_basis = dbasis.nindex1();
  const unsigned n_basis_type = dbasis.nindex2();
  //Find the dimension of the array (Must be the elemental dimension)
  const unsigned n_dim = dim();
 
  //Allocate temporary (stack) storage of the dimension of the element
  double new_derivatives[n_dim];

  //Loop over the number of basis functions
  for(unsigned l=0;l<n_basis;l++)
   {
    //Loop over type of basis functions
    for(unsigned k=0;k<n_basis_type;k++)
     {
      //Loop over the coordinates
      for(unsigned j=0;j<n_dim;j++)
       {
        //Zero the new transformed derivatives
        new_derivatives[j] = 0.0;
        //Do premultiplication by inverse jacobian
        for(unsigned i=0;i<n_dim;i++)
         {
          new_derivatives[j] += inverse_jacobian(j,i)*dbasis(l,k,i);
         }
       }
      //We now copy the new derivatives into the shape functions
      for(unsigned j=0;j<n_dim;j++) {dbasis(l,k,j) = new_derivatives[j];}
     }
   }
 }

//=====================================================================
/// Convert derivatives w.r.t local coordinates to derivatives w.r.t the
/// coordinates used to assemble the inverse jacobian mapping passed as
/// inverse_jacobian, assuming that the mapping is diagonal. This merely
/// saves a few loops, but is probably worth it.
//======================================================================
 void FiniteElement::transform_derivatives_diagonal(const DenseMatrix<double> 
                                                    &inverse_jacobian,
                                                    DShape &dbasis) const
 {
  //Find the number of basis functions and basis function types
  const unsigned n_basis = dbasis.nindex1();
  const unsigned n_basis_type = dbasis.nindex2();
  //Find the dimension of the array (must be the elemental dimension)
  const unsigned n_dim = dim();
 
  //Loop over the number of basis functions
  for(unsigned l=0;l<n_basis;l++)
   {
    //Loop over type of basis functions
    for(unsigned k=0;k<n_basis_type;k++)
     {
      //Loop over the coordinates
      for(unsigned j=0;j<n_dim;j++)
       {
        dbasis(l,k,j) *= inverse_jacobian(j,j);
       }
     }
   }
 }

//=======================================================================
/// Convert derivatives and second derivatives w.r.t local coordinates to 
/// derivatives w.r.t. the coordinates used to assemble the jacobian,
/// inverse_jacobian and jacobian 2 passed. This must be specialised
/// for each dimension, otherwise it gets very ugly
/// Specialisation to one dimension.
//=======================================================================
 template<>
 void FiniteElement::
 transform_second_derivatives_template<1>(const DenseMatrix<double>
                                          &jacobian,
                                          const DenseMatrix<double> 
                                          &inverse_jacobian,
                                          const DenseMatrix<double> 
                                          &jacobian2,
                                          DShape &dbasis,
                                          DShape &d2basis) const
 {
  //Find the number of basis functions and basis function types
  const unsigned n_basis = dbasis.nindex1();
  const unsigned n_basis_type = dbasis.nindex2();
  
  //The second derivatives are easy, because there can be no mixed terms
  //Loop over number of basis functions
  for(unsigned l=0;l<n_basis;l++) 
   {
    //Loop over number of basis function types
    for(unsigned k=0;k<n_basis_type;k++)
     {
      d2basis(l,k,0) = d2basis(l,k,0)/(jacobian(0,0)*jacobian(0,0))
       //Second term comes from taking d/dx of (dpsi/ds / dx/ds)
       - dbasis(l,k,0)*jacobian2(0,0)/
       (jacobian(0,0)*jacobian(0,0)*jacobian(0,0));
     }
   }

  //Assemble the first derivatives (do this last so that we don't
  //overwrite the dphids before we use it in the above)
  transform_derivatives(inverse_jacobian,dbasis);

 }

//=======================================================================
/// Convert derivatives and second derivatives w.r.t local coordinates to 
/// derivatives w.r.t. the coordinates used to assemble the jacobian,
/// inverse_jacobian and jacobian 2 passed. This must be specialised
/// for each dimension, otherwise it gets very ugly.
/// Specialisation to two spatial dimensions
//=======================================================================
 template<>
 void FiniteElement::
 transform_second_derivatives_template<2>(const DenseMatrix<double>
                                          &jacobian,
                                          const DenseMatrix<double> 
                                          &inverse_jacobian,
                                          const DenseMatrix<double> 
                                          &jacobian2,
                                          DShape &dbasis,
                                          DShape &d2basis) const
 {
  //Find the number of basis functions and basis function types
  const unsigned n_basis = dbasis.nindex1();
  const unsigned n_basis_type = dbasis.nindex2();

  //Calculate the determinant
  const double det = jacobian(0,0)*jacobian(1,1) - jacobian(0,1)*jacobian(1,0);

  //Second derivatives ... the approach taken here is to construct
  //dphidX/ds which can then be used to calculate the second derivatives
  //using the relationship d/dX = inverse_jacobian*d/ds

  double ddetds[2];

  ddetds[0] = jacobian2(0,0)*jacobian(1,1) + jacobian(0,0)*jacobian2(2,1) 
   - jacobian2(0,1)*jacobian(1,0) - jacobian(0,1)*jacobian2(2,0);
  ddetds[1] = jacobian2(2,0)*jacobian(1,1) + jacobian(0,0)*jacobian2(1,1)
   - jacobian2(2,1)*jacobian(1,0) - jacobian(0,1)*jacobian2(1,0);

  //Calculate the derivatives of the inverse jacobian wrt the local coordinates
  double dinverse_jacobiands[2][2][2];

  dinverse_jacobiands[0][0][0] = jacobian2(2,1)/det - 
   jacobian(1,1)*ddetds[0]/(det*det);
  dinverse_jacobiands[0][1][0] = -jacobian2(0,1)/det +
   jacobian(0,1)*ddetds[0]/(det*det);
  dinverse_jacobiands[1][0][0] = -jacobian2(2,0)/det +
   jacobian(1,0)*ddetds[0]/(det*det);
  dinverse_jacobiands[1][1][0] = jacobian2(0,0)/det -
   jacobian(0,0)*ddetds[0]/(det*det);

  dinverse_jacobiands[0][0][1] = jacobian2(1,1)/det - 
   jacobian(1,1)*ddetds[1]/(det*det);
  dinverse_jacobiands[0][1][1] = -jacobian2(2,1)/det +
   jacobian(0,1)*ddetds[1]/(det*det);
  dinverse_jacobiands[1][0][1] = -jacobian2(1,0)/det +
   jacobian(1,0)*ddetds[1]/(det*det);
  dinverse_jacobiands[1][1][1] = jacobian2(2,0)/det -
   jacobian(0,0)*ddetds[1]/(det*det);

  //Set up derivatives of dpsidx wrt local coordinates
  DShape dphidXds0(n_basis,n_basis_type,2), dphidXds1(n_basis,n_basis_type,2);
  
  for(unsigned l=0;l<n_basis;l++)
   {
    for(unsigned k=0;k<n_basis_type;k++)
     {
      for(unsigned j=0;j<2;j++)
       {
        //Note that we can't have an inner loop because of the
        //convention I've chosen for the mixed derivatives!
        dphidXds0(l,k,j) = dinverse_jacobiands[j][0][0]*dbasis(l,k,0)
         + dinverse_jacobiands[j][1][0]*dbasis(l,k,1)
         + inverse_jacobian(j,0)*d2basis(l,k,0)
         + inverse_jacobian(j,1)*d2basis(l,k,2);
       
        dphidXds1(l,k,j) = dinverse_jacobiands[j][0][1]*dbasis(l,k,0)
         + dinverse_jacobiands[j][1][1]*dbasis(l,k,1)
         + inverse_jacobian(j,0)*d2basis(l,k,2)
         + inverse_jacobian(j,1)*d2basis(l,k,1);
       }
     }
   }
  
  //Now calculate the DShape d2phidX
  for(unsigned l=0;l<n_basis;l++)
   {
    for(unsigned k=0;k<n_basis_type;k++)
     {
      //Zero dpsidx 
      for(unsigned j=0;j<3;j++) {d2basis(l,k,j) = 0.0;}
      
      //Do premultiplication by inverse jacobian
      for(unsigned i=0;i<2;i++)
       {
        d2basis(l,k,i) = inverse_jacobian(i,0)*dphidXds0(l,k,i)+
         inverse_jacobian(i,1)*dphidXds1(l,k,i);
       }
      //Calculate mixed derivative term
      d2basis(l,k,2) += inverse_jacobian(0,0)*dphidXds0(l,k,1)
       + inverse_jacobian(0,1)*dphidXds1(l,k,1);
     }
   }
 
  //Assemble the first derivatives second, so that we don't
  //overwrite dphids
  transform_derivatives(inverse_jacobian,dbasis);
 }


//=======================================================================
/// Convert derivatives and second derivatives w.r.t local coordinates to 
/// derivatives w.r.t. the coordinates used to assemble the jacobian,
/// inverse_jacobian and jacobian 2 passed. This must be specialised
/// for each dimension, otherwise it gets very ugly
/// Specialisation to one dimension.
//=======================================================================
 template<>
 void FiniteElement::
 transform_second_derivatives_diagonal<1>(const DenseMatrix<double>
                                          &jacobian,
                                          const DenseMatrix<double> 
                                          &inverse_jacobian,
                                          const DenseMatrix<double> 
                                          &jacobian2,
                                          DShape &dbasis,
                                          DShape &d2basis) const
 {
  FiniteElement::
   transform_second_derivatives_template<1>(jacobian,
                                            inverse_jacobian,
                                            jacobian2,dbasis,d2basis);
 }


//=========================================================================
/// Convert second derivatives w.r.t. local coordinates to 
/// second derivatives w.r.t. the coordinates passed in the tensor
/// coordinate. Specialised to two spatial dimension
//=========================================================================
 template<>
 void FiniteElement::
 transform_second_derivatives_diagonal<2>(const DenseMatrix<double>
                                          &jacobian,
                                          const DenseMatrix<double> 
                                          &inverse_jacobian,
                                          const DenseMatrix<double> 
                                          &jacobian2,
                                          DShape &dbasis,
                                          DShape &d2basis) const
 {
  //Find the number of basis functions and basis function types
  const unsigned n_basis = dbasis.nindex1();
  const unsigned n_basis_type = dbasis.nindex2();
 
  //Again we assume that there are no cross terms and that  coordinate
  //i depends only upon local coordinate i

  //Now calculate the DShape d2phidx
  for(unsigned l=0;l<n_basis;l++)
   {
    for(unsigned k=0;k<n_basis_type;k++)
     {
      //Second derivatives
      d2basis(l,k,0) = d2basis(l,k,0)/(jacobian(0,0)*jacobian(0,0))
       - dbasis(l,k,0)*jacobian2(0,0)/
       (jacobian(0,0)*jacobian(0,0)*jacobian(0,0));

      d2basis(l,k,1) = d2basis(l,k,1)/(jacobian(1,1)*jacobian(1,1))
       - dbasis(l,k,1)*jacobian2(1,1)/
       (jacobian(1,1)*jacobian(1,1)*jacobian(1,1));
     
      d2basis(l,k,2) = d2basis(l,k,2)/(jacobian(0,0)*jacobian(1,1));
     }
   }


  //Assemble the first derivatives
  transform_derivatives_diagonal(inverse_jacobian,dbasis);
 }


//=============================================================================
/// \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. This is a template-free general interface, that should be
/// overloaded for efficiency
//============================================================================
 void FiniteElement::transform_second_derivatives(const DenseMatrix<double>
                                                  &jacobian,
                                                  const DenseMatrix<double> 
                                                  &inverse_jacobian,
                                                  const DenseMatrix<double> 
                                                  &jacobian2,
                                                  DShape &dbasis,
                                                  DShape &d2basis) const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Choose the appropriate function based on the dimension of the element
  switch(el_dim)
   {
   case 1:
    transform_second_derivatives_template<1>(jacobian,
                                             inverse_jacobian,jacobian2,
                                             dbasis,d2basis);
    break;
   case 2:
    transform_second_derivatives_template<2>(jacobian,
                                             inverse_jacobian,jacobian2,
                                             dbasis,d2basis);
    break;
         
   case 3:
    throw OomphLibError("Not implemented yet ... maybe one day",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);

    //transform_second_derivatives_template<3>(dphids,d2phids,jacobian,
    //                                         inverse_jacobian,jacobian2,
    //                                         dphidX,d2phidX);
    break;
    //Catch-all default case: issue warning and die
   default:
    std::ostringstream error_stream;
    error_stream
     << "Dimension of the element must be 1,2 or 3, not " << el_dim
     << std::endl;

    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }


//======================================================================
/// \short The destructor cleans up the memory allocated
/// for storage of pointers to nodes. Internal and external data get
/// wiped by the GeneralisedElement destructor; nodes get
/// killed in mesh destructor. 
//=======================================================================
 FiniteElement::~FiniteElement() 
 {
  //Delete the storage allocated for the pointers to the loca nodes
  delete[] Node_pt;
  
  //Delete the storage allocated for the nodal numbering schemes
  if(Nodal_local_eqn)
   {
    delete[] Nodal_local_eqn[0];
    delete[] Nodal_local_eqn;
   }
 }


 //==============================================================
 /// Get the local fraction of the node j in the element; 
 /// vector sets its own size
 //==============================================================
 void FiniteElement::local_fraction_of_node(const unsigned &j, 
                                            Vector<double> &s_fraction)
 {
  //Default implementation is rather dumb
  //Get the local coordinate and scale by local coordinate range
  local_coordinate_of_node(j,s_fraction);
  unsigned n_coordinates = s_fraction.size();
  for(unsigned i=0;i<n_coordinates;i++)
   {
    s_fraction[i] = (s_fraction[i] - s_min())/(s_max() - s_min());
   }
 }

//=======================================================================
/// Set the spatial integration scheme and also calculate the values of the
/// shape functions and their derivatives w.r.t. the local coordinates,
/// placing the values into storage so that they may be re-used, 
/// without recalculation
//=======================================================================
 void FiniteElement::set_integration_scheme(Integral* const &integral_pt) 
 {
  //Assign the integration scheme
  Integral_pt = integral_pt;
 }

//=========================================================================
/// \short Return the shape function stored at the ipt-th integration
/// point.
//=========================================================================
 void FiniteElement::shape_at_knot(const unsigned &ipt, Shape &psi) const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Storage for the local coordinates of the integration point
  Vector<double> s(el_dim); 
  //Set the local coordinate
  for(unsigned i=0;i<el_dim;i++) {s[i] = integral_pt()->knot(ipt,i);}
  //Get the shape function
  shape(s,psi);
 }

//=========================================================================
/// \short Return the shape function and its derivatives w.r.t. the local
/// coordinates at the ipt-th integration point.
//=========================================================================
 void FiniteElement::dshape_local_at_knot(const unsigned &ipt, Shape &psi,
                                          DShape &dpsids) const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Storage for the local coordinates of the integration point
  Vector<double> s(el_dim); 
  //Set the local coordinate
  for(unsigned i=0;i<el_dim;i++) {s[i] = integral_pt()->knot(ipt,i);}
  //Get the shape function and derivatives
  dshape_local(s,psi,dpsids);
 }

//=========================================================================
/// Calculate the shape function and its first and second derivatives 
/// w.r.t. local coordinates at the ipt-th integration point.
/// Numbering:
/// \b 1D: 
/// d2psids(i,0) = \f$ d^2 \psi_j / d s^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$ 
//=========================================================================
 void FiniteElement::d2shape_local_at_knot(const unsigned &ipt, Shape &psi, 
                                           DShape &dpsids, DShape &d2psids) 
  const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Storage for the local coordinates of the integration point
  Vector<double> s(el_dim); 
  //Set the local coordinate
  for(unsigned i=0;i<el_dim;i++) {s[i] = integral_pt()->knot(ipt,i);}
  //Get the shape function and first and second derivatives
  d2shape_local(s,psi,dpsids,d2psids);
 }

//=========================================================================
/// \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.
/// Most general form of the function, but may be over-loaded, if desired
//=========================================================================
 double FiniteElement::dshape_eulerian(const Vector<double> &s, 
                                       Shape &psi,
                                       DShape &dpsi) const
 {
  //Find the element dimension
  const unsigned el_dim = dim();
 
  //Get the values of the shape functions and their local derivatives
  //Temporarily stored in dpsi
  dshape_local(s,psi,dpsi);
  
  //Allocate memory for the inverse jacobian
  DenseMatrix<double> inverse_jacobian(el_dim);
  //Now calculate the inverse jacobian
  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
  
  //Now set the values of the derivatives to be dpsidx
  transform_derivatives(inverse_jacobian,dpsi);
  //Return the determinant of the jacobian
  return det;
 }
 
//========================================================================
/// \short Compute the geometric shape functions and also first
/// derivatives w.r.t. global coordinates at integration point ipt.
/// Most general form of function, but may be over-loaded if desired
//========================================================================
 double FiniteElement::dshape_eulerian_at_knot(const unsigned &ipt,
                                               Shape &psi, 
                                               DShape &dpsi) const
 {
  //Find the element dimension
  const unsigned el_dim = dim();
 
  //Get the values of the shape function and local derivatives
  //Temporarily store it in dpsi
  dshape_local_at_knot(ipt,psi,dpsi);
 
  //Allocate memory for the inverse jacobian
  DenseMatrix<double> inverse_jacobian(el_dim);
  //Now calculate the inverse jacobian
  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);

  //Now set the values of the derivatives to dpsidx
  transform_derivatives(inverse_jacobian,dpsi);
  //Return the determinant of the jacobian
  return det;
 }


//========================================================================
/// \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 FiniteElement::dJ_eulerian_at_knot(
  const unsigned &ipt,
  Shape &psi, 
  DenseMatrix<double> &djacobian_dX) const
 {
  // Find the element dimension
  const unsigned el_dim = dim();
 
  // Get the values of the shape function and local derivatives
  unsigned nnod=nnode();
  DShape dpsi(nnod,el_dim);
  dshape_local_at_knot(ipt,psi,dpsi);
 
  // Allocate memory for the jacobian and the inverse of the jacobian
  DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);

  // Now calculate the inverse jacobian
  const double det = local_to_eulerian_mapping(dpsi,jacobian,inverse_jacobian);

  // Calculate the derivative of the jacobian w.r.t. nodal coordinates
  // Note: must call this before "transform_derivatives(...)" since this
  // function requires dpsids rather than dpsidx
  dJ_eulerian_dnodal_coordinates(jacobian,dpsi,djacobian_dX);

  // Return the determinant of the jacobian
  return det;
 }


//========================================================================
/// \short Compute the geometric shape functions (psi) and first
/// derivatives w.r.t. global coordinates (dpsidx) at integration point
/// ipt. 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.
/// Most general form of function, but may be over-loaded if desired.
//========================================================================
 double FiniteElement::dshape_eulerian_at_knot(
  const unsigned &ipt,
  Shape &psi, 
  DShape &dpsi,
  DenseMatrix<double> &djacobian_dX,
  RankFourTensor<double> &d_dpsidx_dX) const
 {
  // Find the element dimension
  const unsigned el_dim = dim();
 
  // Get the values of the shape function and local derivatives
  // Temporarily store in dpsi
  dshape_local_at_knot(ipt,psi,dpsi);
 
  // Allocate memory for the jacobian and the inverse of the jacobian
  DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);

  // Now calculate the inverse jacobian
  const double det = local_to_eulerian_mapping(dpsi,jacobian,inverse_jacobian);

  // Calculate the derivative of the jacobian w.r.t. nodal coordinates
  // Note: must call this before "transform_derivatives(...)" since this
  // function requires dpsids rather than dpsidx
  dJ_eulerian_dnodal_coordinates(jacobian,dpsi,djacobian_dX);

  // Calculate the derivative of dpsidx w.r.t. nodal coordinates
  // Note: this function also requires dpsids rather than dpsidx
  d_dshape_eulerian_dnodal_coordinates(det,jacobian,djacobian_dX,
                                       inverse_jacobian,dpsi,d_dpsidx_dX);
                                                            
  // Now set the values of the derivatives to dpsidx
  transform_derivatives(inverse_jacobian,dpsi);

  // Return the determinant of the jacobian
  return det;
 }



//===========================================================================
/// \short Compute the geometric shape functions and also first
/// and second derivatives w.r.t. global coordinates at local coordinate s;
/// Also 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 FiniteElement::d2shape_eulerian(const Vector<double> &s, 
                                        Shape &psi, 
                                        DShape &dpsi, 
                                        DShape &d2psi) const
 {
  //Find the values of the indices of the shape functions
  //Locally cached.
  //Find the element dimension
  const unsigned el_dim = dim();
  //Find the number of second derivatives required
  const unsigned n_deriv = N2deriv[el_dim];
 
  //Get the values of the shape function and local derivatives
  d2shape_local(s,psi,dpsi,d2psi);
 
  //Allocate memory for the jacobian and inverse jacobian 
  DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);
  //Calculate the jacobian and inverse jacobian
  const double det = 
   local_to_eulerian_mapping(dpsi,jacobian,inverse_jacobian);

  //Allocate memory for the jacobian of second derivatives
  DenseMatrix<double> jacobian2(n_deriv,el_dim);
  //Assemble the jacobian of second derivatives
  assemble_local_to_eulerian_jacobian2(d2psi,jacobian2);

  //Now set the value of the derivatives
  transform_second_derivatives(jacobian,inverse_jacobian,
                               jacobian2,dpsi,d2psi);
  //Return the determinant of the mapping
  return det;
 }

//===========================================================================
/// \short Compute the geometric shape functions and also first
/// and second derivatives w.r.t. global coordinates at ipt-th integration
/// point
/// Returns Jacobian of mapping from global to local coordinates.
/// This is the most general version, may be overloaded, if desired.
/// 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 FiniteElement::d2shape_eulerian_at_knot(const unsigned &ipt, 
                                                Shape &psi, 
                                                DShape &dpsi, 
                                                DShape &d2psi) const
 {
  //Find the values of the indices of the shape functions
  //Locally cached
  //Find the element dimension
  const unsigned el_dim = dim();
  //Find the number of second derivatives required
  const unsigned n_deriv = N2deriv[el_dim];

  //Get the values of the shape function and local derivatives
  d2shape_local_at_knot(ipt,psi,dpsi,d2psi);
 
  //Allocate memory for the jacobian and inverse jacobian 
  DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);
  //Calculate the jacobian and inverse jacobian
  const double det = 
   local_to_eulerian_mapping(dpsi,jacobian,inverse_jacobian);

  //Allocate memory for the jacobian of second derivatives
  DenseMatrix<double> jacobian2(n_deriv,el_dim);
  //Assemble the jacobian of second derivatives
  assemble_local_to_eulerian_jacobian2(d2psi,jacobian2);

  //Now set the value of the derivatives
  transform_second_derivatives(jacobian,inverse_jacobian,
                               jacobian2,dpsi,d2psi);
  //Return the determinant of the mapping
  return det;
 }



//==========================================================================
/// This function loops over the nodal data of the element, adds the
/// GLOBAL equation numbers to the local-to-global look-up scheme and 
/// fills in the Nodal_local_eqn look-up scheme for the local equation 
/// numbers
/// If the boolean argument is true then pointers to the dofs will be 
/// stored in Dof_pt
//==========================================================================
 void FiniteElement::assign_nodal_local_eqn_numbers(
  const bool &store_local_dof_pt)
 {
  //Find the number of nodes
  const unsigned n_node = nnode();
  //If there are nodes
  if(n_node > 0)
   {
    //Find the number of local equations assigned so far
    unsigned local_eqn_number = ndof();
   
    //We need to find the total number of values stored at the node
    //Initialise to the number of values stored at the first node
    unsigned n_total_values = node_pt(0)->nvalue();
    //Loop over the other nodes and add the values stored
    for(unsigned n=1;n<n_node;n++)
     {n_total_values += node_pt(n)->nvalue();}

    //If allocated delete the old storage
    if(Nodal_local_eqn)
     {
      delete[] Nodal_local_eqn[0];
      delete[] Nodal_local_eqn;
     }
    
    //If there are no values, we are done, null out the storage and
    //return
    if(n_total_values==0) {Nodal_local_eqn=0; return;}

    //Resize the storage for the nodal local equation numbers
    //Firstly allocate pointers to rows for each node
    Nodal_local_eqn = new int*[n_node];
    //Now allocate storage for the equation numbers
    Nodal_local_eqn[0] = new int[n_total_values];
    //initially all local equations are unclassified
    for(unsigned i=0;i<n_total_values;i++)
     {Nodal_local_eqn[0][i] = Data::Is_unclassified;}
    
    //Loop over the remaining rows and set their pointers
    for(unsigned n=1;n<n_node;++n)
     {
      //Initially set the pointer to the i-th row to the pointer
      //to the i-1th row
      Nodal_local_eqn[n] = Nodal_local_eqn[n-1];
      //Now increase the row pointer by the number of values 
      //stored at the i-1th node
      Nodal_local_eqn[n] += Node_pt[n-1]->nvalue();
     }
    

    //A local queue to store the global equation numbers
    std::deque<unsigned long> global_eqn_number_queue;

    //Now loop over the nodes again and assign local equation numbers
    for(unsigned n=0;n<n_node;n++)
     {
      //Pointer to node
      Node* const nod_pt = node_pt(n);

      //Find the number of values stored at the node
      unsigned n_value = nod_pt->nvalue();
     
      //Loop over the number of values
      for(unsigned j=0;j<n_value;j++)
       {
        //Get the GLOBAL equation number
        long eqn_number = nod_pt->eqn_number(j);
        //If the GLOBAL equation number is positive (a free variable)
        if(eqn_number >= 0)
         {
          //Add the GLOBAL equation number to the queue
          global_eqn_number_queue.push_back(eqn_number);
          //Add pointer to the dof to the queue if required
          if(store_local_dof_pt)
           {
            GeneralisedElement::Dof_pt_deque.push_back(
             nod_pt->value_pt(j));
           }
          //Add the local equation number to the local scheme
          Nodal_local_eqn[n][j] = local_eqn_number;
          //Increase the local number
          local_eqn_number++;
         }
        else
         {
          //Set the local scheme to be pinned
          Nodal_local_eqn[n][j] = Data::Is_pinned;
         }
       }
     }
    
    //Now add our global equations numbers to the internal element storage
    add_global_eqn_numbers(global_eqn_number_queue,
                           GeneralisedElement::Dof_pt_deque);
    //Clear the memory used in the deque
    if(store_local_dof_pt)
     {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}
   }
 }


//============================================================================
/// This function calculates the entries of Jacobian matrix, used in 
/// the Newton method, associated with the nodal degrees of freedom.
/// It does this using finite differences, 
/// rather than an analytical formulation, so can be done in total generality.
//==========================================================================
 void FiniteElement::
 fill_in_jacobian_from_nodal_by_fd(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
 {
  //Find the number of nodes
  const unsigned n_node = nnode();
  //If there aren't any nodes, then return straight awayy
  if(n_node == 0) {return;}

  //Call the update function to ensure that the element is in
  //a consistent state before finite differencing starts
  update_before_nodal_fd();

  //Find the number of dofs in the element
  const unsigned n_dof = ndof();
  //Create newres vector
  Vector<double> newres(n_dof);

  //Integer storage for local unknown
  int local_unknown=0;
  
  //Use the default finite difference step
  const double fd_step = Default_fd_jacobian_step;

  //Loop over the nodes
  for(unsigned n=0;n<n_node;n++)
   {
    //Get the number of values stored at the node
    const unsigned n_value = node_pt(n)->nvalue();
   
    //Loop over the number of values
    for(unsigned i=0;i<n_value;i++)
     {
      //Get the local equation number
      local_unknown = nodal_local_eqn(n,i);
      //If it's not pinned
      if(local_unknown >= 0)
       {
        //Store a pointer to the nodal data value
        double* const value_pt = node_pt(n)->value_pt(i);
        
        //Save the old value of the Nodal data
        const double old_var = *value_pt;
       
        //Increment the value of the Nodal data
        *value_pt += fd_step;

        //Now update any slaved variables
        update_in_nodal_fd(i);
        
        //Calculate the new residuals
        get_residuals(newres);
       
        //Do finite differences
        for(unsigned m=0;m<n_dof;m++)
         {
          double sum = (newres[m] - residuals[m])/fd_step;
          //Stick the entry into the Jacobian matrix
          jacobian(m,local_unknown) = sum;
         }
       
        //Reset the Nodal data
        *value_pt = old_var;

        //Reset any slaved variables
        reset_in_nodal_fd(i);
       }
     }
   }
  
  //End of finite difference loop
  //Final reset of any slaved data
  reset_after_nodal_fd();
 }




//=======================================================================
/// 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}
////=======================================================================
void FiniteElement::get_dresidual_dnodal_coordinates(
 RankThreeTensor<double>& dresidual_dnodal_coordinates)
{
 // Number of nodes
 unsigned n_nod=nnode();
 
 // If the element has no nodes (why??!!) return straightaway
 if (n_nod==0) return;
 
 // Get dimension from first node
 unsigned dim_nod=node_pt(0)->ndim();
 
 // Number of dofs
 unsigned n_dof=ndof();
 
 // Get reference residual
 Vector<double> res(n_dof);
 Vector<double> res_pls(n_dof);
 get_residuals(res);
  
  // FD step 
  double eps_fd=GeneralisedElement::Default_fd_jacobian_step;

  // Do FD loop
  for (unsigned j=0;j<n_nod;j++)
   {
    // Get node
    Node* nod_pt=node_pt(j);

    // Loop over coordinate directions
    for (unsigned i=0;i<dim_nod;i++)
     {
      // Make backup
      double backup=nod_pt->x(i);

      // Do FD step. No node update required as we're
      // attacking the coordinate directly...
      nod_pt->x(i)+=eps_fd;

      // Perform auxiliary node update function
      nod_pt->perform_auxiliary_node_update_fct();
  
      // Get advanced residual
      get_residuals(res_pls);
      
      // Fill in FD entries [Loop order is "wrong" here as l is the
      // slow index but this is in a function that's costly anyway
      // and gives us the fastest loop outside where these tensor
      // is actually used.]
      for (unsigned l=0;l<n_dof;l++)
       {
        dresidual_dnodal_coordinates(l,i,j)=(res_pls[l]-res[l])/eps_fd;
       }

      // Reset coordinate. No node update required as we're
      // attacking the coordinate directly...
      nod_pt->x(i)=backup;

      // Perform auxiliary node update function
      nod_pt->perform_auxiliary_node_update_fct();
      
     }
   }
  
 }


//===============================================================
/// Return the number of the node located at *node_pt 
/// if this node is in the element, else return \f$ -1 \f$ 
//===============================================================
 int FiniteElement::get_node_number(Node* const &global_node_pt) const
 {
   //Initialise the number to -1
   int number = -1;
   //Find the number of nodes
   unsigned n_node = nnode();
#ifdef PARANOID
   {
    //Error check that node does not appear in element more than once
    unsigned count=0;
    //Storage for the local node numbers of the element
    std::vector<int> local_node_number;
    //Loop over the nodes
    for(unsigned i=0;i<n_node;i++)
     {
      //If the node is present increase the counter
      //and store the local node number
      if(node_pt(i)==global_node_pt) 
       {
        ++count;
        local_node_number.push_back(i);
       }
     }
    
    //If the node appears more than once, complain
    if(count > 1)
     {
      std::ostringstream error_stream;
      error_stream << "Node " << global_node_pt << " appears " 
                   << count << " times in an element." << std::endl
                   << "In positions: ";
      for(std::vector<int>::iterator it=local_node_number.begin();
          it!=local_node_number.end();++it)
       {
        error_stream << *it << " ";
       }
      error_stream << std::endl
                   << "That seems very odd." << std::endl;
      
      throw OomphLibError(error_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif
   
   //Loop over the nodes
   for(unsigned i=0;i<n_node;i++)
    {
     //If the passed node pointer is present in the element
     //set number to be its local node number
     if(node_pt(i)==global_node_pt)
      {      
       number=i;
       break;
      }
    }
   
   //Return the node number
   return number;
 }
 

 //==========================================================================
 /// \short If there is a node at the local coordinate, s, return the pointer 
 /// to  the node. If not return 0. Note that this is a default, brute
 /// force implementation, can almost certainly be made more efficient for
 /// specific elements.
 //==========================================================================
 Node* FiniteElement::get_node_at_local_coordinate(const Vector<double> &s) const
  {
   //Locally cache the tolerance
   const double tol = Node_location_tolerance;
   Vector<double> s_node;
   //Locally cache the member data
   const unsigned el_dim = Elemental_dimension;
   const unsigned n_node = Nnode;
   //Loop over the nodes
   for(unsigned n=0;n<n_node;n++)
    {
     bool Match = true;
     //Find the local coordinate of the node
     local_coordinate_of_node(n,s_node);
     for(unsigned i=0;i<el_dim;i++)
      {
       //Calculate the difference between coordinates
       //and if it's bigger than our tolerance 
       //break out of the (inner)loop
       if(std::fabs(s[i] - s_node[i]) > tol)
        {
         Match = false;
         break;
        }
      }
     //If we haven't complained then we have a match
     if(Match) {return node_pt(n);}
    }
   //If we get here, we have no match
   return 0;
  }
   
//======================================================================
/// 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. NOTE: All computed in terms of zeta!
//======================================================================
 void FiniteElement::get_centre_of_gravity_and_max_radius_in_terms_of_zeta(
  Vector<double>& cog, double& max_radius) const
 {
  // Initialise
  cog.resize(Elemental_dimension);
  max_radius=0.0;

  // Get cog
  unsigned nnod=nnode();
  for (unsigned j=0;j<nnod;j++)
   {
    for (unsigned i=0;i<Elemental_dimension;i++)
     {
      cog[i]+=zeta_nodal(j,0,i);
     }
   }
  for (unsigned i=0;i<Elemental_dimension;i++)
   {
    cog[i]/=double(nnod);
   }

  // Get max distance
  for (unsigned j=0;j<nnod;j++)
   {
    double dist_squared=0.0;
    for (unsigned i=0;i<Elemental_dimension;i++)
     {
      dist_squared+=(cog[i]-zeta_nodal(j,0,i))*(cog[i]-zeta_nodal(j,0,i));
     }
    if (dist_squared>max_radius) max_radius=dist_squared;
   }
  max_radius=sqrt(max_radius);
 }

//======================================================================
/// Return FE interpolated coordinate x[i] at local coordinate s
//======================================================================
 double FiniteElement::interpolated_x(const Vector<double> &s, 
                                      const unsigned &i) const
 {
  //Find the number of nodes
  const unsigned n_node = nnode();
  //Find the number of positional types
  const unsigned n_position_type = nnodal_position_type();
  //Assign storage for the local shape function
  Shape psi(n_node,n_position_type);
  //Find the values of shape function
  shape(s,psi);

  //Initialise value of x
  double interpolated_x = 0.0;
  //Loop over the local nodes
  for(unsigned l=0;l<n_node;l++) 
   {
    //Loop over the number of dofs
    for(unsigned k=0;k<n_position_type;k++)
     {
      interpolated_x += nodal_position_gen(l,k,i)*psi(l,k);
     }
   }
   
  return(interpolated_x);
 }

//=========================================================================
/// Return FE interpolated coordinate x[i] at local coordinate s
/// at previous timestep t (t=0: present; t>0: previous timestep)
//========================================================================
 double FiniteElement::interpolated_x(const unsigned &t,
                                      const Vector<double> &s, 
                                      const unsigned &i) const
 {
  //Find the number of nodes
  const unsigned n_node = nnode();
  //Find the number of positional types
  const unsigned n_position_type = nnodal_position_type();

  //Assign storage for the local shape function
  Shape psi(n_node,n_position_type);
  //Find the values of shape function
  shape(s,psi);

  //Initialise value of x
  double interpolated_x = 0.0;
  //Loop over the local nodes
  for(unsigned l=0;l<n_node;l++) 
   {
    //Loop over the number of dofs
    for(unsigned k=0;k<n_position_type;k++)
     {
      interpolated_x += nodal_position_gen(t,l,k,i)*psi(l,k);
     }
   }
   
  return(interpolated_x);
 }

//=======================================================================
/// Return FE interpolated position x[] at local coordinate s as Vector
//=======================================================================
 void FiniteElement::interpolated_x(const Vector<double> &s, Vector<double> &x)
  const
 {
   //Find the number of nodes
  const unsigned n_node = nnode();
  //Find the number of positional types
  const unsigned n_position_type = nnodal_position_type();
  //Find the dimension stored in the node
  const unsigned nodal_dim = nodal_dimension();

  //Assign storage for the local shape function
  Shape psi(n_node,n_position_type);
  //Find the values of shape function
  shape(s,psi);

  //Loop over the dimensions
  for(unsigned i=0;i<nodal_dim;i++)
   {
    //Initilialise value of x[i] to zero
    x[i] = 0.0;
    //Loop over the local nodes
    for(unsigned l=0;l<n_node;l++) 
     {
      //Loop over the number of dofs
      for(unsigned k=0;k<n_position_type;k++)
       {
        x[i] += nodal_position_gen(l,k,i)*psi(l,k);
       }
     }
   }
 }

//==========================================================================
/// Return FE interpolated position x[] at local coordinate s
/// at previous timestep t as Vector (t=0: present; t>0: previous timestep)
//==========================================================================
 void FiniteElement::interpolated_x(const unsigned &t, 
                                    const Vector<double> &s, 
                                    Vector<double>& x) const
 {
  //Find the number of nodes
  const unsigned n_node = nnode();
  //Find the number of positional types
  const unsigned n_position_type = nnodal_position_type();
  //Find the dimensions of the nodes
  const unsigned nodal_dim = nodal_dimension();

  //Assign storage for the local shape function
  Shape psi(n_node,n_position_type);
  //Find the values of shape function
  shape(s,psi);

  //Loop over the dimensions
  for(unsigned i=0;i<nodal_dim;i++)
   {
    //Initilialise value of x[i] to zero
    x[i] = 0.0;
    //Loop over the local nodes
    for(unsigned l=0;l<n_node;l++) 
     {
      //Loop over the number of dofs
      for(unsigned k=0;k<n_position_type;k++)
       {
        x[i] += nodal_position_gen(t,l,k,i)*psi(l,k);
       }
     }
   }
 }

//========================================================================
/// \short Calculate the determinant of the 
/// Jacobian of the mapping between local and global
/// coordinates at the position. Works directly from the base vectors
/// without assuming that coordinates match spatial dimension. Will
/// be overloaded in FaceElements, in which the elemental dimension does
/// not match the spatial dimension. WARNING: this is always positive
/// and cannot be used to check if the element is inverted, say!
//========================================================================
 double FiniteElement::J_eulerian(const Vector<double> &s) const
 {
  //Find the number of nodes and position types
  const unsigned n_node = nnode();
  const unsigned n_position_type = nnodal_position_type();
  //Find the dimension of the node and element
  const unsigned n_dim_node = nodal_dimension();
  const unsigned n_dim_element = dim();

  //Set up dummy memory for the shape functions
  Shape psi(n_node,n_position_type);
  DShape dpsids(n_node,n_position_type,n_dim_element);
  //Get the shape functions and local derivatives
  dshape_local(s,psi,dpsids);
  
  //Right calculate the base vectors
  DenseMatrix<double> interpolated_G(n_dim_element,n_dim_node);
  assemble_eulerian_base_vectors(dpsids,interpolated_G);
 
  //Calculate the metric tensor of the element
  DenseMatrix<double> G(n_dim_element,n_dim_element,0.0);
  for(unsigned i=0;i<n_dim_element;i++)
   {
    for(unsigned j=0;j<n_dim_element;j++)
     {
      for(unsigned k=0;k<n_dim_node;k++) 
       {
        G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
       }
     }
   }

  //Calculate the determinant of the metric tensor
  double det = 0.0;
  switch(n_dim_element)
   {
   case 0:
    throw OomphLibError("Cannot calculate J_eulerian() for point element\n",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
    break;
   case 1:
    det = G(0,0);
    break;
   case 2:
    det = G(0,0)*G(1,1) - G(0,1)*G(1,0);
    break;
   case 3:
    det = G(0,0)*G(1,1)*G(2,2) + G(0,1)*G(1,2)*G(2,0) + G(0,2)*G(1,0)*G(2,1)
     - G(0,0)*G(1,2)*G(2,1) - G(0,1)*G(1,0)*G(2,2) - G(0,2)*G(1,1)*G(2,0);
    break;
   default:
    oomph_info << "More than 3 dimensions in J_eulerian()" << std::endl;
    break;
   }

  //Return the Jacobian (square-root of the determinant of the metric tensor)
  return sqrt(det);
 }

//========================================================================
/// \short Compute the Jacobian of the mapping between the local and global
/// coordinates at the ipt-th integration point
//========================================================================
 double FiniteElement::J_eulerian_at_knot(const unsigned &ipt)
  const
 {
  //Find the number of nodes
  const unsigned n_node = nnode();
  //Find the number of position types
  const unsigned n_position_type = nnodal_position_type();
  //Find the dimension of the node and element
  const unsigned n_dim_node = nodal_dimension();
  const unsigned n_dim_element = dim();

  //Set up dummy memory for the shape functions
  Shape psi(n_node,n_position_type);
  DShape dpsids(n_node,n_position_type,n_dim_element);
  //Get the shape functions and local derivatives at the knot
  //This call may use the stored versions, which is why this entire
  //function doesn't just call J_eulerian(s), after reading out s from
  //the knots.
  dshape_local_at_knot(ipt,psi,dpsids);

  //Right calculate the base vectors
  DenseMatrix<double> interpolated_G(n_dim_element,n_dim_node);
  assemble_eulerian_base_vectors(dpsids,interpolated_G);
 
  //Calculate the metric tensor of the element
  DenseMatrix<double> G(n_dim_element,n_dim_element,0.0);
  for(unsigned i=0;i<n_dim_element;i++)
   {
    for(unsigned j=0;j<n_dim_element;j++)
     {
      for(unsigned k=0;k<n_dim_node;k++) 
       {
        G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
       }
     }
   }

  //Calculate the determinant of the metric tensor
  double det = 0.0;
  switch(n_dim_element)
   {
   case 0:
    throw OomphLibError("Cannot calculate J_eulerian() for point element\n",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
    break;
   case 1:
    det = G(0,0);
    break;
   case 2:
    det = G(0,0)*G(1,1) - G(0,1)*G(1,0);
    break;
   case 3:
    det = G(0,0)*G(1,1)*G(2,2) + G(0,1)*G(1,2)*G(2,0) + G(0,2)*G(1,0)*G(2,1)
     - G(0,0)*G(1,2)*G(2,1) - G(0,1)*G(1,0)*G(2,2) - G(0,2)*G(1,1)*G(2,0);
    break;
   default:
    oomph_info << "More than 3 dimensions in J_eulerian()" << std::endl;
    break;
   }
 
  //Return the Jacobian (square-root of the determinant of the metric tensor)
  return sqrt(det);
 }

//========================================================================
/// \short Check that Jacobian of mapping between local and Eulerian
/// coordinates at all integration points is positive.
//========================================================================
void FiniteElement::check_J_eulerian_at_knots(bool& passed) const
 {

  // Bypass check deep down in the guts...
  bool backup=FiniteElement::Accept_negative_jacobian;
  FiniteElement::Accept_negative_jacobian=true;

  // Let's be optimistic...
  passed=true;

  //Find the number of nodes
  const unsigned n_node = nnode();

  //Find the number of position types
  const unsigned n_position_type = nnodal_position_type();

  //Find the dimension of the element
  const unsigned n_dim_element = dim();

  //Set up dummy memory for the shape functions
  Shape psi(n_node,n_position_type);
  DShape dpsi(n_node,n_dim_element);

  unsigned nintpt=integral_pt()->nweight();
  for (unsigned ipt=0;ipt<nintpt;ipt++)
   {
    double jac=dshape_eulerian_at_knot(ipt,psi,dpsi);

    // Are we dead yet?
    if (jac<=0.0)
     {
      passed=false;
      
      // Reset
      FiniteElement::Accept_negative_jacobian=backup;

      return;
     }
   }

  // Reset
  FiniteElement::Accept_negative_jacobian=backup;

 }

//====================================================================
/// Calculate the size of the element (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 FiniteElement::size() const
 {
  //Initialise the sum to zero
  double sum = 0.0;
 
  //Loop over the integration points
  const unsigned n_intpt = integral_pt()->nweight();
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {
    //Get the integral weight
    double w = integral_pt()->weight(ipt);
    //Get the value of the Jacobian of the mapping to global coordinates
    double J = J_eulerian_at_knot(ipt);

    //Add the product to the sum
    sum += w*J;
   }
 
  //Return the answer
  return(sum);
 }

//==================================================================
/// Integrate Vector-valued time-dep function over element
//==================================================================
 void FiniteElement::
 integrate_fct(FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt,
               const double& time,
               Vector<double>& integral)
 {
  //Initialise all components of integral Vector and setup integrand vector
  const unsigned ncomponents=integral.size();
  Vector<double> integrand(ncomponents);
  for (unsigned i=0;i<ncomponents;i++) {integral[i]=0.0;}

  // Figure out the global (Eulerian) spatial dimension of the
  // element 
  const unsigned n_dim_eulerian = nodal_dimension();
  
  // Allocate Vector of global Eulerian coordinates
  Vector<double> x(n_dim_eulerian);

  // Set the value of n_intpt
  const unsigned n_intpt = integral_pt()->nweight();

  // Vector of local coordinates
  const unsigned n_dim = dim();
  Vector<double> s(n_dim);

  //Loop over the integration points
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {
    //Assign the values of s 
    for(unsigned i=0;i<n_dim;i++) {s[i] = integral_pt()->knot(ipt,i);}

    //Assign the values of the global Eulerian coordinates
    for(unsigned i=0;i<n_dim_eulerian;i++) {x[i] = interpolated_x(s,i);}

    //Get the integral weight
    double w = integral_pt()->weight(ipt);
     
    //Get Jacobian of mapping
    double J = J_eulerian(s);
     
    // Evaluate the integrand at the knot points
    integrand_fct_pt(time,x,integrand);

    //Add to components of integral Vector
    for (unsigned i=0;i<ncomponents;i++) {integral[i]+=integrand[i]*w*J;}
   }
 }

//==================================================================
/// Integrate Vector-valued function over element
//==================================================================
 void FiniteElement::
 integrate_fct(FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt,
               Vector<double>& integral)
 {
  //Initialise all components of integral Vector
  const unsigned ncomponents=integral.size();
  Vector<double> integrand(ncomponents);
  for (unsigned i=0;i<ncomponents;i++) {integral[i]=0.0;}

  // Figure out the global (Eulerian) spatial dimension of the
  // element by checking the Eulerian dimension of the nodes
  const unsigned n_dim_eulerian = nodal_dimension();
  
  // Allocate Vector of global Eulerian coordinates
  Vector<double> x(n_dim_eulerian);

  // Set the value of n_intpt
  const unsigned n_intpt = integral_pt()->nweight();

  // Vector of local coordinates
  const unsigned n_dim = dim();
  Vector<double> s(n_dim);

  //Loop over the integration points
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {
    //Assign the values of s 
    for(unsigned i=0;i<n_dim;i++) {s[i] = integral_pt()->knot(ipt,i);}

    //Assign the values of the global Eulerian coordinates
    for(unsigned i=0;i<n_dim_eulerian;i++) {x[i] = interpolated_x(s,i);}

    //Get the integral weight
    double w = integral_pt()->weight(ipt);
     
    //Get Jacobian of mapping
    double J = J_eulerian(s);
     
    // Evaluate the integrand at the knot points
    integrand_fct_pt(x,integrand);

    //Add to components of integral Vector
    for (unsigned i=0;i<ncomponents;i++) {integral[i]+=integrand[i]*w*J;}
   }
 }

//==========================================================================
/// Self-test: Have all internal values been classified as 
/// pinned/unpinned? Has pointer to spatial integration scheme
/// been set? Return 0 if OK. 
//==========================================================================
 unsigned FiniteElement::self_test()
 {
  // Initialise
  bool passed=true;

  if(GeneralisedElement::self_test()!=0) {passed=false;}

  // Check that pointer to spatial integration scheme has been set
  if(integral_pt()==0)
   {
    passed=false;

    OomphLibWarning(
     "Pointer to spatial integration scheme has not been set.",
     "FiniteElement::self_test()",
     OOMPH_EXCEPTION_LOCATION);
   }

  //If the dimension of the element is zero (point element), there
  //is not jacobian
  const unsigned dim_el = dim();

  if(dim_el> 0)
   {
    // Loop over integration points to check sign of Jacobian
    //-------------------------------------------------------
   
    //Set the value of n_intpt
    const unsigned n_intpt = integral_pt()->nweight();
 
    //Set the Vector to hold local coordinates
    Vector<double> s(dim_el);

    //Find the number of local nodes
    const unsigned n_node = nnode();
    const unsigned n_position_type = nnodal_position_type();
    //Set up memory for the shape and test functions
    Shape psi(n_node,n_position_type);
    DShape dpsidx(n_node,dim_el);
            
    // Jacobian
    double jacobian;
   

    // Two ways of testing for negative Jacobian for non-FaceElements
    unsigned ntest=1; 
   
    // For FaceElements checking the Jacobian via dpsidx doesn't
    // make sense
    FiniteElement* tmp_pt=const_cast<FiniteElement*>(this);
    FaceElement* face_el_pt=dynamic_cast<FaceElement*>(tmp_pt);
    if (face_el_pt==0)
     {
      ntest=2;
     }
   
    // For now overwrite -- the stuff above fails for Bretherton.
    // Not sure why.
    ntest=1;

    //Loop over the integration points
    for(unsigned ipt=0;ipt<n_intpt;ipt++)
     {
      //Assign values of s
      for(unsigned i=0;i<dim_el;i++)
       {
        s[i] = integral_pt()->knot(ipt,i);
       }
     

      // Do tests
      for (unsigned test=0;test<ntest;test++)
       {
   
        switch (test)
         {
         
         case 0:
         
          // Get the jacobian from the mapping between local and Eulerian
          // coordinates
          jacobian = J_eulerian(s);
         
          break;
         
         case 1:

          //Call the geometrical shape functions and derivatives  
          //This also computes the Jacobian by a slightly different 
          //method
          jacobian = dshape_eulerian_at_knot(ipt,psi,dpsidx);

          break;

         default:

          throw OomphLibError("Never get here",
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
         }

       
        //Check for a singular jacobian
        if(std::fabs(jacobian) < 1.0e-16)
         {
          std::ostringstream warning_stream;
          warning_stream << "Determinant of Jacobian matrix is zero at ipt "
                         << ipt << std::endl;
          OomphLibWarning(warning_stream.str(),
                          "FiniteElement::self_test()",
                          OOMPH_EXCEPTION_LOCATION);
          passed = false;
          //Skip the next test
          continue;
         }
       
        //Check sign of Jacobian
        if ((Accept_negative_jacobian==false) && (jacobian < 0.0))
         {
          std::ostringstream warning_stream;
          warning_stream << "Jacobian negative at integration point ipt=" 
                         << ipt << std::endl;
          warning_stream 
           << "If you think that this is what you want you may: " << std::endl
           << "set the (static) flag " 
           << "FiniteElement::Accept_negative_jacobian to be true" << std::endl;
         
          OomphLibWarning(warning_stream.str(),
                          "FiniteElement::self_test()",
                          OOMPH_EXCEPTION_LOCATION);
          passed=false;
         }
       
       } // end of loop over two tests

     }
   } //End of non-zero dimension check
   
  // Return verdict
  if (passed) {return 0;}
  else {return 1;}
 }




//=======================================================================
/// Return the t-th time-derivative of the 
/// i-th FE-interpolated Eulerian coordinate at 
/// local coordinate s.
//=======================================================================
 double FiniteElement::interpolated_dxdt(const Vector<double> &s, 
                                         const unsigned &i,
                                         const unsigned &t_deriv)
 {
  //Find the number of nodes and positions (locally cached)
  const unsigned n_node = nnode();
  const unsigned n_position_type = nnodal_position_type();
  // Get shape functions: specify # of nodes, # of positional dofs
  Shape psi(n_node,n_position_type);
  shape(s,psi);
 
  // Initialise 
  double drdt=0.0;
 
  // Assemble time derivative
  //Loop over nodes
  for(unsigned l=0;l<n_node;l++) 
   {
    //Loop over types of dof
    for(unsigned k=0;k<n_position_type;k++)
     {
      drdt+=dnodal_position_gen_dt(t_deriv,l,k,i)*psi(l,k);  
     }
   }
  return drdt;
 }



//=======================================================================
/// Compute t-th time-derivative of the
/// FE-interpolated Eulerian coordinate vector at
/// local coordinate s. 
//=======================================================================
 void FiniteElement::interpolated_dxdt(const Vector<double> &s, 
                                       const unsigned &t_deriv,
                                       Vector<double>& dxdt)
 {
  //Find the number of nodes and positions (locally cached)
  const unsigned n_node = nnode();
  const unsigned n_position_type = nnodal_position_type();
  const unsigned nodal_dim = nodal_dimension();

  // Get shape functions: specify # of nodes, # of positional dofs
  Shape psi(n_node,n_position_type);
  shape(s,psi);
 
  // Loop over directions
  for (unsigned i=0;i<nodal_dim;i++)
   {
    // Initialise 
    dxdt[i]=0.0;
    
    // Assemble time derivative
    //Loop over nodes
    for(unsigned l=0;l<n_node;l++) 
     {
      //Loop over types of dof
      for(unsigned k=0;k<n_position_type;k++)
       {
        dxdt[i]+=dnodal_position_gen_dt(t_deriv,l,k,i)*psi(l,k);  
       }
     }
   }
 }

//============================================================================
/// 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 FiniteElement::interpolated_zeta(const Vector<double> &s,  
                                       Vector<double> &zeta) const       
{
 //If there is a macro element use it
 if(Macro_elem_pt!=0) {this->get_x_from_macro_element(s,zeta);}
 //Otherwise interpolate zeta_nodal using the shape functions
 else
  {
   //Find the number of nodes
   const unsigned n_node = this->nnode();
   //Find the number of positional types
   const unsigned n_position_type = this->nnodal_position_type();
   //Storage for the shape functions
   Shape psi(n_node,n_position_type);
   //Get the values of the shape functions at the local coordinate s
   this->shape(s,psi);
   
   //Find the number of coordinates
   const unsigned ncoord = this->dim();
   //Initialise the value of zeta to zero
   for(unsigned i=0;i<ncoord;i++) {zeta[i] = 0.0;}

   //Add the contributions from each nodal dof to the interpolated value
   //of zeta.
   for(unsigned l=0;l<n_node;l++)
    {
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Locally cache the value of the shape function
       const double psi_ = psi(l,k);
       for(unsigned i=0;i<ncoord;i++)
        {
         zeta[i] += this->zeta_nodal(l,k,i)*psi_;
        }
      }
    }
  }
}

//==========================================================================
/// 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. 
/// This is achieved in generality by using Newton's method to find the value 
/// of the local coordinate, s, such that
/// interpolated_zeta(s) is equal 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.
/// Setting the optional bool argument to true means that the coordinate
/// argument "s" is used as the initial guess. (Default is false).
//=========================================================================
void FiniteElement::locate_zeta(const Vector<double> &zeta,
                                GeomObject*& geom_object_pt, Vector<double> &s,
                                const bool& use_coordinate_as_initial_guess)
{
 
 //Find the number of coordinates, the dimension of the element
 //This must be the same for the local and intrinsic global coordinate
 unsigned ncoord = this->dim();
 
 // Fast return based on centre of gravity and max. radius of any nodal
 // point 
 if (Locate_zeta_helpers::Radius_multiplier_for_fast_exit_from_locate_zeta>0.0)
  {
   Vector<double> cog(ncoord);
   double max_radius=0.0;
   get_centre_of_gravity_and_max_radius_in_terms_of_zeta(cog,max_radius);
   
   // Get radius
   double radius=0.0;
   for (unsigned i=0;i<ncoord;i++)
    {
     radius+=(cog[i]-zeta[i])*(cog[i]-zeta[i]);
    }
   radius=sqrt(radius);
   if (radius>
       Locate_zeta_helpers::Radius_multiplier_for_fast_exit_from_locate_zeta*
       max_radius)
    {   
     geom_object_pt=0;
     return;    
    }
  }
  
 //Assign storage for the vector and matrix used in Newton's method
 Vector<double> dx(ncoord,0.0);
 DenseDoubleMatrix jacobian(ncoord,ncoord,0.0);
 
 // // Make a list of (equally-spaced) local coordinates inside the element
 // unsigned n_list=Locate_zeta_helpers::N_local_points; 

 // double list_space=(1.0/(double(n_list)-1.0))*(s_max()-s_min());
 
 // Possible initial guesses for Newton's method
 Vector<Vector<double> > s_list;
 
 // If the boolean argument use_coordinate_as_initial_guess was set
 // to true then we don't need to initialise s
 if (!use_coordinate_as_initial_guess)
  {
   //Vector of local coordinates
   Vector<double> s_c(ncoord);
   
   // Loop over plot points
   unsigned num_plot_points=nplot_points(Locate_zeta_helpers::N_local_points);
   for (unsigned iplot=0;iplot<num_plot_points;iplot++)
    {
     // Get local coordinates of plot point
     get_s_plot(iplot,Locate_zeta_helpers::N_local_points,s_c);
     s_list.push_back(s_c);
    }
  }
 
 //Counter for the number of Newton steps
 unsigned count=0;
 
 //Control flag for the Newton loop
 bool keep_going=true;
 
 //Storage for the interpolated value of x
 Vector<double> inter_x(ncoord);
 
 // Ifwe have specified the coordinate already
 if (use_coordinate_as_initial_guess)
  {
   //Get the value of x at the initial guess
   this->interpolated_zeta(s,inter_x);
   
   //Set up the residuals
   for(unsigned i=0;i<ncoord;i++) {dx[i] = zeta[i] - inter_x[i];}
  }
 else
  {
   // Find the smallest residual from the list of coordinates made earlier
   double my_min_resid=DBL_MAX;
   Vector<double> s_local(ncoord);
   Vector<double> work_x(ncoord);
   Vector<double> work_dx(ncoord);
   
   unsigned n_list_coord=s_list.size();
   
   for (unsigned i_coord=0; i_coord<n_list_coord; i_coord++)
    {
     for (unsigned i=0;i<ncoord;i++)
      {
       s_local[i]=s_list[i_coord][i];
      }
     // get_x for this coordinate
     this->interpolated_zeta(s_local,work_x);
     
     // calculate residuals
     for(unsigned i=0;i<ncoord;i++)
      {
       work_dx[i] = zeta[i] - work_x[i];
      }
     
     double maxres =
      std::fabs(*std::max_element(work_dx.begin(),work_dx.end(),
                                  AbsCmp<double>()));
     
     // test against previous residuals
     if (maxres<my_min_resid)
      {
       my_min_resid=maxres;
       dx=work_dx;
       inter_x=work_x;
       s=s_local;
      }
     
    }
  }
 
 //Main Newton Loop
 do    // start of do while loop
  {
   //Increase loop counter
   count++;
   
   //Bail out if necessary (without an error for now...)
   if(count > Locate_zeta_helpers::Max_newton_iterations)
    {
     keep_going=false;
     continue;
    }
   
   //If it's the first time round the loop, check the initial residuals
   if(count==1)
    {
     double maxres =
      std::fabs(*std::max_element(dx.begin(),dx.end(),AbsCmp<double>()));
     
     //If it's small enough exit
     if(maxres < Locate_zeta_helpers::Newton_tolerance)
      {
       keep_going=false;
       continue;
      }
    }
   
   //Is there a macro element? If so, assemble the Jacobian by FD-ing
   if (macro_elem_pt()!=0)
    {
     // Assemble jacobian on the fly by finite differencing
     Vector<double> work_s=s;
     Vector<double> r=inter_x; // i.e. the result of previous call to get_x
     
     // Finite difference step
     double fd_step=GeneralisedElement::Default_fd_jacobian_step;
     
     // Storage for calculated r from incremented s
     Vector<double> work_r(ncoord,0.0);
     
     // Loop over s coordinates
     for (unsigned i=0; i<ncoord; i++)
      {
       // Increment work_s by a small amount
       work_s[i]+=fd_step;
       
       // Calculate work_r from macro element
       this->interpolated_zeta(work_s,work_r);
       
       // Loop over r to fill Jacobian
       for (unsigned j=0; j<ncoord; j++)
        {
         jacobian(j,i)=-(work_r[j]-r[j])/fd_step;
        }
       
       // Reset work_s
       work_s[i]=s[i];
       
      }
     
    }
   else // no macro element, so compute Jacobian with shape functions etc.
    {
     //Compute the entries of the Jacobian matrix
     unsigned n_node = this->nnode();
     unsigned n_position_type = this->nnodal_position_type();
     Shape psi(n_node,n_position_type);
     DShape dpsids(n_node,n_position_type,ncoord);
     
     //Get the local shape functions and their derivatives
     dshape_local(s,psi,dpsids);
     
     //Calculate the values of dxds
     DenseMatrix<double> interpolated_dxds(ncoord,ncoord,0.0);
     
// MH: No longer needed
//          //This implementation will only work for n_position_type=1
//          //since the function nodal_position_gen does not yet exist
// #ifdef PARANOID
//          if (n_position_type!=1)
//           {
//            std::ostringstream error_stream;
//            error_stream << "This implementation does not exist yet;\n"
//                         << "it currently uses raw_nodal_position_gen\n"
//                         << "which does not take hangingness into account\n"
//                         << "It will work if n_position_type=1\n";
//            throw OomphLibError(error_stream.str(),
//         OOMPH_CURRENT_FUNCTION,
//                                OOMPH_EXCEPTION_LOCATION);
//           }
// #endif
     
     // Loop over the nodes
     for(unsigned l=0;l<n_node;l++)
      {
       // Loop over position type even though it should be 1; the
       // functionality for n_position_type>1 will exist in the future
       for(unsigned k=0;k<n_position_type;k++)
        {
         // Add the contribution from the nodal coordinates to the matrix
         for(unsigned i=0;i<ncoord;i++)
          {
           for(unsigned j=0;j<ncoord;j++)
            {
             interpolated_dxds(i,j) +=
              this->zeta_nodal(l,k,i)*dpsids(l,k,j);
            }
          }
        }
      }
     
     //The entries of the Jacobian matrix are merely dresiduals/ds
     //i.e. \f$ -dx/ds \f$
     for(unsigned i=0;i<ncoord;i++)
      {
       for(unsigned j=0;j<ncoord;j++)
        {
         jacobian(i,j) = - interpolated_dxds(i,j);
        }
      }
     
    }
   
   //Now solve the damn thing
   try
    {
     jacobian.solve(dx);
    }
   catch(OomphLibError &error)
    {
     // I've caught the error so shut up!
     error.disable_error_message();
#ifdef PARANOID
     oomph_info << "Error in linear solve for "
                << "FiniteElement::locate_zeta()"
                << std::endl;
     oomph_info << "Should not affect the result!" << std::endl;
#endif
    }
   
   //Add the correction to the local coordinates
   for(unsigned i=0;i<ncoord;i++) {s[i] -= dx[i];}
   
   //Get the new residuals
   this->interpolated_zeta(s,inter_x);
   for(unsigned i=0;i<ncoord;i++) {dx[i] = zeta[i] - inter_x[i];}
   
   //Get the maximum residuals
   double maxres =
    std::fabs(*std::max_element(dx.begin(),dx.end(),AbsCmp<double>()));
   
   //If we have converged jump straight to the test at the end of the loop
   if(maxres < Locate_zeta_helpers::Newton_tolerance)
    {
     keep_going=false;
     continue;
    }
  }
 while(keep_going);
 
 //Test whether the local coordinates are valid or not
 bool valid=local_coord_is_valid(s);

 // If not valid, experimentally push back into element
 // and see if the result is still valid (within the Newton tolerance)
 if (!valid)
  {
   move_local_coord_back_into_element(s);

   // Check residuals again
   this->interpolated_zeta(s,inter_x);
   for(unsigned i=0;i<ncoord;i++) {dx[i] = zeta[i] - inter_x[i];}
   
   //Get the maximum residuals
   double maxres =
    std::fabs(*std::max_element(dx.begin(),dx.end(),AbsCmp<double>()));
   
   // Are we still OK?
   if(maxres > Locate_zeta_helpers::Newton_tolerance)
    {
     // oomph_info 
     // << "Pushing back inside has violated the Newton tolerance: max_res = "
     // << maxres << std::endl;
     geom_object_pt=0;
     return;
    }
  }
 
 // It is also possible now that it may not have converged "correctly", 
 // i.e. count is greater than Max_newton_iterations
 if (count > Locate_zeta_helpers::Max_newton_iterations)
  {
   // Don't trust the current answer, return null
   geom_object_pt=0;
   return;
  }
  
 //Otherwise the required point is located in "this" element:
 geom_object_pt = this;
}



//=======================================================================
/// Loop over all nodes in the element and update their positions
/// using each node's (algebraic) update function
//=======================================================================
 void FiniteElement::node_update()
 {
  const unsigned n_node = nnode();
  for(unsigned n=0;n<n_node;n++) {node_pt(n)->node_update();}
 }

//======================================================================
/// 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 (a) the pointer to a Data object
/// and (b) 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. Note that the geometric data, 
/// which includes the positions
/// of SolidNodes, is treated separately by the function 
/// \c identify_geometric_data()
//======================================================================
void FiniteElement::identify_field_data_for_interactions(
 std::set<std::pair<Data*,unsigned> > &paired_field_data)
{
 //Loop over all internal data
 const unsigned n_internal = this->ninternal_data();
 for(unsigned n=0;n<n_internal;n++)
  {
   //Cache the data pointer
   Data* const dat_pt = this->internal_data_pt(n);
   //Find the number of data values stored in the data object
   const unsigned n_value = dat_pt->nvalue();
   //Add the index of each data value and the pointer to the set
   //of pairs
   for(unsigned i=0;i<n_value;i++)
    {
     paired_field_data.insert(std::make_pair(dat_pt,i));
    }
  }

 //Loop over all the nodes
 const unsigned n_node = this->nnode();
 for(unsigned n=0;n<n_node;n++)
  {
   //Cache the node pointer
   Node* const nod_pt = this->node_pt(n);
   //Find the number of values stored at the node
   const unsigned n_value = nod_pt->nvalue();
   //Add the index of each data value and the pointer to the set
   //of pairs
   for(unsigned i=0;i<n_value;i++) 
    {
     paired_field_data.insert(std::make_pair(nod_pt,i));
    }
  }
}

void FiniteElement::build_face_element(const int &face_index, 
                                       FaceElement* face_element_pt)
{
 //Set the nodal dimension
 face_element_pt->set_nodal_dimension(nodal_dimension());

 //Set the pointer to the orginal "bulk" element
 face_element_pt->bulk_element_pt()=this;

#ifdef OOMPH_HAS_MPI
 // Pass on non-halo proc ID
 face_element_pt->set_halo(Non_halo_proc_ID);
#endif

 // Set the face index
 face_element_pt->face_index() = face_index;

 // Get number of bulk nodes on a face of this element
 const unsigned nnode_face = nnode_on_face();

 // Set the function pointer for face coordinate to bulk coordinate
 // mapping
 face_element_pt->face_to_bulk_coordinate_fct_pt() =
  face_to_bulk_coordinate_fct_pt(face_index);

 // Set the function pointer for the derivative of the face coordinate to
 // bulk coordinate mapping
 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
  bulk_coordinate_derivatives_fct_pt(face_index);

 // Resize storage for the number of values originally stored each of the
 // face element's nodes.
 face_element_pt->nbulk_value_resize(nnode_face);

 // Resize storage for the bulk node numbers corresponding to the face
 // element's nodes.
 face_element_pt->bulk_node_number_resize(nnode_face);

 // Copy bulk_node_numbers and nbulk_values
 for(unsigned i=0; i<nnode_face; i++)
  {
   // Find the corresponding bulk node's number
   unsigned bulk_number = get_bulk_node_number(face_index, i);

   // Assign the pointer and number into the face element
   face_element_pt->node_pt(i) = node_pt(bulk_number);
   face_element_pt->bulk_node_number(i) = bulk_number;

   // Set the number of values originally stored at this node
   face_element_pt->nbulk_value(i) = required_nvalue(bulk_number); 
  }

 // Set the outer unit normal sign
 face_element_pt->normal_sign() = face_outer_unit_normal_sign(face_index);

}

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

// Initialise the static variable. 
// Do not ignore warning for discontinuous tangent vectors.
bool FaceElement::Ignore_discontinuous_tangent_warning = false;


//========================================================================
/// Output boundary coordinate zeta
//========================================================================
void FaceElement::output_zeta(std::ostream &outfile, const unsigned& nplot)
{
 //Vector of local coordinates
 unsigned n_dim=dim();
 Vector<double> s(n_dim);
 
 // Tecplot header info
 outfile << tecplot_zone_string(nplot);
 
 // Loop over plot points
 unsigned num_plot_points=nplot_points(nplot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   // Get local coordinates of plot point
   get_s_plot(iplot,nplot,s);
   
   // Spatial coordinates are one higher
   for(unsigned i=0;i<n_dim+1;i++) 
    {
     outfile << interpolated_x(s,i) << " ";
    }
   
   // Boundary coordinate
   Vector<double> zeta(n_dim);
   interpolated_zeta(s,zeta);
   for(unsigned i=0;i<n_dim;i++) 
    {
     outfile << zeta[i] << " ";
    }
   outfile << std::endl;   
  }
 
 // Write tecplot footer (e.g. FE connectivity lists)
 write_tecplot_zone_footer(outfile,nplot);
}


//========================================================================
/// \short Calculate the determinant of the 
/// Jacobian of the mapping between local and global
/// coordinates at the position s. Overloaded from FiniteElement.
//========================================================================
 double FaceElement::J_eulerian(const Vector<double> &s) const
 {

  //Find out the sptial dimension of the element
  unsigned n_dim_el = this->dim();

  // Bail out if we're in a point element -- not sure what
  // J_eulerian actually is, but this is harmless 
  if (n_dim_el==0) return 1.0;
   
  //Find out how many nodes there are
  unsigned n_node = nnode();
  
  //Find out how many positional dofs there are
  unsigned n_position_type = this->nnodal_position_type();
  
  //Find out the dimension of the node
  unsigned n_dim = this->nodal_dimension();
    
  //Set up memory for the shape functions
  Shape psi(n_node,n_position_type);
  DShape dpsids(n_node,n_position_type,n_dim_el); 

  //Only need to call the local derivatives
  dshape_local(s,psi,dpsids);
  
  //Also calculate the surface Vectors (derivatives wrt local coordinates)
  DenseMatrix<double> interpolated_A(n_dim_el,n_dim,0.0);   
  
  //Calculate positions and derivatives
  for(unsigned l=0;l<n_node;l++) 
   {
    //Loop over positional dofs
    for(unsigned k=0;k<n_position_type;k++)
     {
      //Loop over coordinates
      for(unsigned i=0;i<n_dim;i++)
       {
        //Loop over LOCAL derivative directions, to calculate the tangent(s)
        for(unsigned j=0;j<n_dim_el;j++)
         {
          interpolated_A(j,i) += 
           nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,j);
         }
       }
     }
   }
  //Now find the local deformed metric tensor from the tangent Vectors
  DenseMatrix<double> A(n_dim_el,n_dim_el,0.0);
  for(unsigned i=0;i<n_dim_el;i++)
   {
    for(unsigned j=0;j<n_dim_el;j++)
     {
      //Take the dot product
      for(unsigned k=0;k<n_dim;k++)
       { 
        A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
       }
     }
   }
  
  //Find the determinant of the metric tensor
  double Adet =0.0;
  switch(n_dim_el) 
   {
   case 1:
    Adet = A(0,0);
    break;
   case 2:
    Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
    break;
   default:
    throw 
     OomphLibError("Wrong dimension in FaceElement",
                   "FaceElement::J_eulerian()",
                   OOMPH_EXCEPTION_LOCATION);
   }
  
  // Return 
  return sqrt(Adet);
 }


//========================================================================
/// \short Compute the Jacobian of the mapping between the local and global
/// coordinates at the ipt-th integration point. Overloaded from 
/// FiniteElement.
//========================================================================
 double FaceElement::J_eulerian_at_knot(const unsigned &ipt)
  const
 {
  //Find the number of nodes
  const unsigned n_node = nnode();

  //Find the number of position types
  const unsigned n_position_type = nnodal_position_type();

  //Find the dimension of the node and element
  const unsigned n_dim = nodal_dimension();
  const unsigned n_dim_el = dim();

  //Set up dummy memory for the shape functions
  Shape psi(n_node,n_position_type);
  DShape dpsids(n_node,n_position_type,n_dim_el);

  //Get the shape functions and local derivatives at the knot
  //This call may use the stored versions, which is why this entire
  //function doesn't just call J_eulerian(s), after reading out s from
  //the knots.
  dshape_local_at_knot(ipt,psi,dpsids);

  //Also calculate the surface Vectors (derivatives wrt local coordinates)
  DenseMatrix<double> interpolated_A(n_dim_el,n_dim,0.0);   
  
  //Calculate positions and derivatives
  for(unsigned l=0;l<n_node;l++) 
   {
    //Loop over positional dofs
    for(unsigned k=0;k<n_position_type;k++)
     {
      //Loop over coordinates
      for(unsigned i=0;i<n_dim;i++)
       {
        //Loop over LOCAL derivative directions, to calculate the tangent(s)
        for(unsigned j=0;j<n_dim_el;j++)
         {
          interpolated_A(j,i) += 
           nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,j);
         }
       }
     }
   }

  //Now find the local deformed metric tensor from the tangent Vectors
  DenseMatrix<double> A(n_dim_el,n_dim_el,0.0);
  for(unsigned i=0;i<n_dim_el;i++)
   {
    for(unsigned j=0;j<n_dim_el;j++)
     {
      //Take the dot product
      for(unsigned k=0;k<n_dim;k++)
       { 
        A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
       }
     }
   }
  
  //Find the determinant of the metric tensor
  double Adet =0.0;
  switch(n_dim_el) 
   {
   case 1:
    Adet = A(0,0);
    break;
   case 2:
    Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
    break;
   default:
    throw 
     OomphLibError("Wrong dimension in FaceElement",
                   "FaceElement::J_eulerian_at_knot()",
                   OOMPH_EXCEPTION_LOCATION);
   }
  
  // Return 
  return sqrt(Adet);
 }

//========================================================================
/// \short Check that Jacobian of mapping between local and Eulerian
/// coordinates at all integration points is positive.
//========================================================================
 void FaceElement::check_J_eulerian_at_knots(bool& passed) const
 {
  // Let's be optimistic...
  passed=true;

  //Find the number of nodes
  const unsigned n_node = nnode();

  //Find the number of position types
  const unsigned n_position_type = nnodal_position_type();

  //Find the dimension of the node and element
  const unsigned n_dim = nodal_dimension();
  const unsigned n_dim_el = dim();

  //Set up dummy memory for the shape functions
  Shape psi(n_node,n_position_type);
  DShape dpsids(n_node,n_position_type,n_dim_el);


  unsigned nintpt=integral_pt()->nweight();
  for (unsigned ipt=0;ipt<nintpt;ipt++)
   {
    
    //Get the shape functions and local derivatives at the knot
    //This call may use the stored versions, which is why this entire
    //function doesn't just call J_eulerian(s), after reading out s from
    //the knots.
    dshape_local_at_knot(ipt,psi,dpsids);
    
    //Also calculate the surface Vectors (derivatives wrt local coordinates)
    DenseMatrix<double> interpolated_A(n_dim_el,n_dim,0.0);   
    
    //Calculate positions and derivatives
    for(unsigned l=0;l<n_node;l++) 
     {
      //Loop over positional dofs
      for(unsigned k=0;k<n_position_type;k++)
       {
        //Loop over coordinates
        for(unsigned i=0;i<n_dim;i++)
         {
          //Loop over LOCAL derivative directions, to calculate the tangent(s)
          for(unsigned j=0;j<n_dim_el;j++)
           {
            interpolated_A(j,i) += 
             nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,j);
           }
         }
       }
     }
    
    //Now find the local deformed metric tensor from the tangent Vectors
    DenseMatrix<double> A(n_dim_el,n_dim_el,0.0);
    for(unsigned i=0;i<n_dim_el;i++)
     {
      for(unsigned j=0;j<n_dim_el;j++)
       {
        //Take the dot product
        for(unsigned k=0;k<n_dim;k++)
         { 
          A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
         }
       }
     }
    
    //Find the determinant of the metric tensor
    double Adet =0.0;
    switch(n_dim_el) 
     {
     case 1:
      Adet = A(0,0);
      break;
     case 2:
      Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
      break;
     default:
      throw 
       OomphLibError("Wrong dimension in FaceElement",
                     "FaceElement::J_eulerian_at_knot()",
                     OOMPH_EXCEPTION_LOCATION);
     }
    

    // Are we dead yet?
    if (Adet<=0.0)
     {
      passed=false;
      return;
     }
   }

 }

//=======================================================================
/// \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 FaceElement::continuous_tangent_and_outer_unit_normal
(const Vector<double> &s, Vector<Vector<double> > &tang_vec, 
 Vector<double> &unit_normal) const
{
 //Find the spatial dimension of the FaceElement
 const unsigned element_dim = dim();

 //Find the overall dimension of the problem 
 //(assume that it's the same for all nodes)
 const unsigned spatial_dim = nodal_dimension();
 
#ifdef PARANOID
 //Check the number of local coordinates passed
 if(s.size()!=element_dim)
  {
   std::ostringstream error_stream;
   error_stream
    << "Local coordinate s passed to outer_unit_normal() has dimension " 
    << s.size() << std::endl
    << "but element dimension is " << element_dim << std::endl;
   
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

 // Check that if the Tangent_direction_pt is set, then
 // it is the same length as the spatial dimension.
 if(Tangent_direction_pt != 0 && Tangent_direction_pt->size()!=spatial_dim)
  {
   std::ostringstream error_stream;
   error_stream
    << "Tangent direction vector has dimension " 
    << Tangent_direction_pt->size() << std::endl
    << "but spatial dimension is " << spatial_dim << std::endl;

   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
 
 //Check the dimension of the normal vector
 if(unit_normal.size()!=spatial_dim)
  {
   std::ostringstream error_stream;
   error_stream
    << "Unit normal passed to outer_unit_normal() has dimension " 
    << unit_normal.size() << std::endl
    << "but spatial dimension is " << spatial_dim << std::endl;
   
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }


 // The number of tangent vectors.
 unsigned ntangent_vec = tang_vec.size();

 // For the tangent vector, 
 // if element_dim = 2, tang_vec is a 2D Vector of length 3.
 // if element_dim = 1, tang_vec is a 1D Vector of length 2.
 // if element_dim = 0, tang_vec is a 1D Vector of length 1.
 switch(element_dim)
  {
   //Point element, derived from a 1D element.
  case 0:
  {
   // Check that tang_vec is a 1D vector.
   if(ntangent_vec != 1)
    {
     std::ostringstream error_stream;
     error_stream
      << "This is a 0D FaceElement, we need one tangent vector.\n" 
      << "You have given me " << tang_vec.size() << " vector(s).\n";
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
  }
  break;

  //Line element, derived from a 2D element.
  case 1:
  {
   // Check that tang_vec is a 1D vector.
   if(ntangent_vec != 1)
    {
     std::ostringstream error_stream;
     error_stream
      << "This is a 1D FaceElement, we need one tangent vector.\n" 
      << "You have given me " << tang_vec.size() << " vector(s).\n";
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
  }
  break;

  //Plane element, derived from 3D element.
  case 2:
  {
   // Check that tang_vec is a 2D vector.
   if(ntangent_vec != 2)
    {
     std::ostringstream error_stream;
     error_stream
      << "This is a 2D FaceElement, we need two tangent vectors.\n" 
      << "You have given me " << tang_vec.size() << " vector(s).\n";
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
  }
  break;
  // There are no other elements!
  default:

   throw OomphLibError(
    "Cannot have a FaceElement with dimension higher than 2",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
   break;
  }

 // Check the lengths of each sub vectors.
 for(unsigned vec_i = 0; vec_i < ntangent_vec; vec_i++)
  {
   if(tang_vec[vec_i].size() != spatial_dim)
    {
     std::ostringstream error_stream;
     error_stream
      << "This problem has " << spatial_dim << " spatial dimensions.\n" 
      << "But the " << vec_i << " tangent vector has size " 
      << tang_vec[vec_i].size() << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
  }

#endif   


 //Now let's consider the different element dimensions
 switch(element_dim)
  {

   //Point element, derived from a 1D element.
  case 0:
  {
   std::ostringstream error_stream;
   error_stream
    << "It is unclear how to calculate a tangent and normal vectors for "
    << "a point element.\n"; 
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
  break;

  //Line element, derived from a 2D element, in this case
  //the normal is a mess of cross products
  //We need an interior direction, so we must find the local
  //derivatives in the BULK element
  case 1:
  {
   //Find the number of nodes in the bulk element
   const unsigned n_node_bulk = Bulk_element_pt->nnode();
   //Find the number of position types in the bulk element
   const unsigned n_position_type_bulk = 
    Bulk_element_pt->nnodal_position_type();
    
   //Construct the local coordinate in the bulk element
   Vector<double> s_bulk(2);
   //Get the local coordinates in the bulk element
   get_local_coordinate_in_bulk(s,s_bulk);
    
   //Allocate storage for the shape functions and their derivatives wrt
   //local coordinates
   Shape psi(n_node_bulk,n_position_type_bulk);
   DShape dpsids(n_node_bulk,n_position_type_bulk,2);
   //Get the value of the shape functions at the given local coordinate
   Bulk_element_pt->dshape_local(s_bulk,psi,dpsids);
 
   //Calculate all derivatives of the spatial coordinates 
   //wrt local coordinates
   DenseMatrix<double> interpolated_dxds(2,spatial_dim,0.0);
    
   //Loop over all parent nodes
   for(unsigned l=0;l<n_node_bulk;l++)
    {
     //Loop over all position types in the bulk
     for(unsigned k=0;k<n_position_type_bulk;k++)
      {
       //Loop over derivative direction
       for(unsigned j=0;j<2;j++)
        {
         //Loop over coordinate directions
         for(unsigned i=0;i<spatial_dim;i++)
          {
           //Compute the spatial derivative
           interpolated_dxds(j,i) += 
            Bulk_element_pt->nodal_position_gen(l,k,i)*dpsids(l,k,j);
          }
        }
      }
    }
    
   //Initialise the tangent, interior tangent and normal vectors to zero
   //The idea is that even if the element is in a two-dimensional space,
   //the normal cannot be calculated without embedding the element in three
   //dimensions, in which case, the tangent and interior tangent will have
   //zero z-components.
   Vector<double> tangent(3,0.0), interior_tangent(3,0.0), normal(3,0.0);
    
   //We must get the relationship between the coordinate along the face
   //and the local coordinates in the bulk element
   //We must also find an interior direction
   DenseMatrix<double> dsbulk_dsface(2,1,0.0);
   unsigned interior_direction=0;
   get_ds_bulk_ds_face(s,dsbulk_dsface,interior_direction);
   //Load in the values for the tangents
   for(unsigned i=0;i<spatial_dim;i++)
    {
     //Tangent to the face is the derivative wrt to the face coordinate
     //which is calculated using dsbulk_dsface and the chain rule
     tangent[i] = interpolated_dxds(0,i)*dsbulk_dsface(0,0)
      + interpolated_dxds(1,i)*dsbulk_dsface(1,0);
     //Interior tangent to the face is the derivative in the interior 
     //direction
     interior_tangent[i] = interpolated_dxds(interior_direction,i);
    }

   //Now the (3D) normal to the element is the interior tangent 
   //crossed with the tangent
   VectorHelpers::cross(interior_tangent,tangent,normal);

   //We find the line normal by crossing the element normal with the tangent
   Vector<double> full_normal(3);
   VectorHelpers::cross(normal,tangent,full_normal);
    
   //Copy the appropriate entries into the unit normal
   //Two or Three depending upon the spatial dimension of the system
   for(unsigned i=0;i<spatial_dim;i++) {unit_normal[i] = full_normal[i];}

   //Finally normalise unit normal and multiply by the Normal_sign
   double length = VectorHelpers::magnitude(unit_normal);
   for(unsigned i=0;i<spatial_dim;i++) 
    {unit_normal[i] *= Normal_sign/length;}

   // Create the tangent vector
   tang_vec[0][0]=-unit_normal[1];
   tang_vec[0][1]=unit_normal[0]; 
  }
  break;

  //Plane element, derived from 3D element, in this case the normal
  //is just the cross product of the two surface tangents
  //We assume, therefore, that we have three spatial coordinates
  //and two surface coordinates
  //Then we need only to get the derivatives wrt the local coordinates
  //in this face element
  case 2:
  {
#ifdef PARANOID
   //Check that we actually have three spatial dimensions
   if(spatial_dim != 3)
    {
     std::ostringstream error_stream;
     error_stream << "There are only " << spatial_dim
                  << "coordinates at the nodes of this 2D FaceElement,\n"
                  << "which must have come from a 3D Bulk element\n";
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   //Find the number of nodes in the element
   const unsigned n_node  = this->nnode();
   //Find the number of position types
   const unsigned n_position_type = this->nnodal_position_type();
    
   //Allocate storage for the shape functions and their derivatives wrt
   //local coordinates
   Shape psi(n_node,n_position_type);
   DShape dpsids(n_node,n_position_type,2);
   //Get the value of the shape functions at the given local coordinate
   this->dshape_local(s,psi,dpsids);
 
   //Calculate all derivatives of the spatial coordinates 
   //wrt local coordinates
   Vector<Vector<double> > interpolated_dxds(2,Vector<double>(3,0));
    
   //Loop over all nodes
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop over all position types
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Loop over derivative directions
       for(unsigned j=0;j<2;j++)
        {
         //Loop over coordinate directions
         for(unsigned i=0;i<3;i++)
          {
           //Compute the spatial derivative
           //Remember that we need to translate the position type
           //to its location in the bulk node
           interpolated_dxds[j][i] += 
            this->nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,j);
          }
        }
      }
    }

   // We now take the cross product of the two normal vectors
   VectorHelpers::cross(interpolated_dxds[0],interpolated_dxds[1],
                        unit_normal);

   // Finally normalise unit normal
   double normal_length = VectorHelpers::magnitude(unit_normal);

   for(unsigned i=0;i<spatial_dim;i++) 
    {unit_normal[i] *= Normal_sign/normal_length;}

   // Next we create the continuous tangent vectors!
   if(Tangent_direction_pt != 0)
    // There is a general direction that the first tangent vector should follow.
    {

     // Project the Tangent direction vector onto the surface.
     // Project Tangent_direction_pt D onto the plane P defined by
     // T1 and T2:
     // proj_P(D) = proj_T1(D) + proj_T2(D), where D is Tangent_direction_pt,
     // recall that proj_u(v) = (u.v)/(u.u) * u

     // Get the direction vector. The vector is NOT copied! :)
     Vector<double> &direction_vector = *Tangent_direction_pt;

#ifdef PARANOID
     // Check that the angle between the direction vector and the normal
     // is not less than 10 degrees
     if(VectorHelpers::angle(direction_vector,unit_normal) < 
        (10.0 * MathematicalConstants::Pi/180.0))
      {
       std::ostringstream err_stream;
       err_stream << "The angle between the unit normal and the \n"
                  << "direction vector is less than 10 degrees.";
       throw OomphLibError(err_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
#endif

     // Calculate the two scalings, (u.v) / (u.u)
     double t1_scaling 
      = VectorHelpers::dot(interpolated_dxds[0],direction_vector) / 
      VectorHelpers::dot(interpolated_dxds[0],interpolated_dxds[0]);

     double t2_scaling 
      = VectorHelpers::dot(interpolated_dxds[1],direction_vector) / 
      VectorHelpers::dot(interpolated_dxds[1],interpolated_dxds[1]);

     // Finish off the projection.
     tang_vec[0][0] = t1_scaling*interpolated_dxds[0][0] 
      + t2_scaling*interpolated_dxds[1][0];
     tang_vec[0][1] = t1_scaling*interpolated_dxds[0][1] 
      + t2_scaling*interpolated_dxds[1][1];
     tang_vec[0][2] = t1_scaling*interpolated_dxds[0][2] 
      + t2_scaling*interpolated_dxds[1][2];

     // The second tangent vector is the cross product of
     // tang_vec[0] and the normal vector N.
     VectorHelpers::cross(tang_vec[0],unit_normal,tang_vec[1]);

     // Normalise the tangent vectors 
     for(unsigned vec_i=0; vec_i<2; vec_i++)
      {
       // Get the length...
       double tang_length = VectorHelpers::magnitude(tang_vec[vec_i]);

       for(unsigned dim_i=0;dim_i<spatial_dim;dim_i++) 
        {tang_vec[vec_i][dim_i] /= tang_length;}
      }
    }
   else
    // A direction for the first tangent vector is not supplied, we do our best
    // to keep it continuous. But if the normal vector is aligned with the 
    // z axis, then the tangent vector may "flip", even within elements.
    // This is because we check this by comparing doubles.
    // The code is copied from impose_parallel_outflow_element.h,
    // NO modification is done except for a few variable name changes.
    {
     double a,b,c;
     a=unit_normal[0];
     b=unit_normal[1];
     c=unit_normal[2]; 

     if(a!=0.0 || b!=0.0)
      {
       double a_sq=a*a;
       double b_sq=b*b;
       double c_sq=c*c;
       
       tang_vec[0][0]=-b/sqrt(a_sq+b_sq);
       tang_vec[0][1]= a/sqrt(a_sq+b_sq);
       tang_vec[0][2]=0;
          
       double z=(a_sq +b_sq)
        /sqrt(a_sq*c_sq +b_sq*c_sq +(a_sq +b_sq)*(a_sq +b_sq)) ;
        
       tang_vec[1][0]=-(a*c*z)/(a_sq + b_sq) ;
       tang_vec[1][1]= -(b*c*z)/(a_sq + b_sq);
       tang_vec[1][2]= z;
       // NB : we didn't use the fact that N is normalized,
       // that's why we have these unsimplified formulas
      }
     else if (c!=0.0)
      {
       if(!FaceElement::Ignore_discontinuous_tangent_warning)
        {
         std::ostringstream warn_stream;
         warn_stream << "I have detected that your normal vector is [0,0,1].\n"
                     << "Since this function compares against 0.0, you may\n"
                     << "get discontinuous tangent vectors.";
         OomphLibWarning(warn_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
        }

       tang_vec[0][0]= 1.0;
       tang_vec[0][1]= 0.0;       
       tang_vec[0][2]= 0.0;
          
       tang_vec[1][0]=0.0;
       tang_vec[1][1]= 1.0;       
       tang_vec[1][2]= 0.0;
      }
     else
      {
       throw 
        OomphLibError("You have a zero normal vector!! ",
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
      }
    }

  }
  break;

  default:

   throw OomphLibError(
    "Cannot have a FaceElement with dimension higher than 2",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
   break;
  }
   
}

//=======================================================================
/// \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 FaceElement::continuous_tangent_and_outer_unit_normal
(const unsigned &ipt, Vector<Vector<double> > &tang_vec,
 Vector<double> &unit_normal) const
{
 //Find the dimension of the element
 const unsigned element_dim = dim();
 //Find the local coordinates of the ipt-th integration point
 Vector<double> s(element_dim);
 for(unsigned i=0;i<element_dim;i++) {s[i] = integral_pt()->knot(ipt,i);}
 //Call the outer unit normal function
 continuous_tangent_and_outer_unit_normal(s, tang_vec,unit_normal);
}

//=======================================================================
/// Compute the outer unit normal at the specified local coordinate
//=======================================================================
 void FaceElement::outer_unit_normal(const Vector<double> &s,
                                     Vector<double> &unit_normal) const
 {
  //Find the spatial dimension of the FaceElement
  const unsigned element_dim = dim();

  //Find the overall dimension of the problem 
  //(assume that it's the same for all nodes)
  const unsigned spatial_dim = nodal_dimension();
 
#ifdef PARANOID
  //Check the number of local coordinates passed
  if(s.size()!=element_dim)
   {
    std::ostringstream error_stream;
    error_stream
     << "Local coordinate s passed to outer_unit_normal() has dimension " 
     << s.size() << std::endl
     << "but element dimension is " << element_dim << std::endl;
   
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 
  //Check the dimension of the normal vector
  if(unit_normal.size()!=spatial_dim)
   {
    std::ostringstream error_stream;
    error_stream
     << "Unit normal passed to outer_unit_normal() has dimension " 
     << unit_normal.size() << std::endl
     << "but spatial dimension is " << spatial_dim << std::endl;
   
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif   

/*  //The spatial dimension of the bulk element will be element_dim+1
  const unsigned bulk_dim = element_dim + 1;

  //Find the number of nodes in the bulk element
  const unsigned n_node_bulk = Bulk_element_pt->nnode();
  //Find the number of position types in the bulk element
  const unsigned n_position_type_bulk = 
   Bulk_element_pt->nnodal_position_type();

  //Construct the local coordinate in the bulk element
  Vector<double> s_bulk(bulk_dim);
  //Set the value of the bulk coordinate that is fixed on the face
  //s_bulk[s_fixed_index()] = s_fixed_value();

  //Set the other bulk coordinates
  //for(unsigned i=0;i<element_dim;i++) {s_bulk[bulk_s_index(i)] = s[i];} 

  get_local_coordinate_in_bulk(s,s_bulk);

  //Allocate storage for the shape functions and their derivatives wrt
  //local coordinates
  Shape psi(n_node_bulk,n_position_type_bulk);
  DShape dpsids(n_node_bulk,n_position_type_bulk,bulk_dim);
  //Get the value of the shape functions at the given local coordinate
  Bulk_element_pt->dshape_local(s_bulk,psi,dpsids);
 
  //Calculate all derivatives of the spatial coordinates wrt local coordinates
  DenseMatrix<double> interpolated_dxds(bulk_dim,spatial_dim);
  //Initialise to zero
  for(unsigned j=0;j<bulk_dim;j++)
   {for(unsigned i=0;i<spatial_dim;i++) {interpolated_dxds(j,i) = 0.0;}}

  //Loop over all parent nodes
  for(unsigned l=0;l<n_node_bulk;l++)
   {
    //Loop over all position types in the bulk
    for(unsigned k=0;k<n_position_type_bulk;k++)
     {
      //Loop over derivative direction
      for(unsigned j=0;j<bulk_dim;j++)
       {
        //Loop over coordinate directions
        for(unsigned i=0;i<spatial_dim;i++)
         {
          //Compute the spatial derivative
          interpolated_dxds(j,i) += 
           Bulk_element_pt->nodal_position_gen(l,k,i)*dpsids(l,k,j);
         }
       }
     }
     }*/

  //Now let's consider the different element dimensions
  switch(element_dim)
   {
    //Point element, derived from a 1D element, in this case
    //the normal is merely the tangent to the bulk element
    //and there is only one free coordinate in the bulk element
    //Hence we will need to calculate the derivatives wrt the
    //local coordinates in the BULK element.
   case 0:
   {
    //Find the number of nodes in the bulk element
    const unsigned n_node_bulk = Bulk_element_pt->nnode();
    //Find the number of position types in the bulk element
    const unsigned n_position_type_bulk = 
     Bulk_element_pt->nnodal_position_type();

    //Construct the local coordinate in the bulk element
    Vector<double> s_bulk(1);

    //Get the local coordinates in the bulk element
    get_local_coordinate_in_bulk(s,s_bulk);
    
    //Allocate storage for the shape functions and their derivatives wrt
    //local coordinates
    Shape psi(n_node_bulk,n_position_type_bulk);
    DShape dpsids(n_node_bulk,n_position_type_bulk,1);
    //Get the value of the shape functions at the given local coordinate
    Bulk_element_pt->dshape_local(s_bulk,psi,dpsids);
 
    //Calculate all derivatives of the spatial coordinates wrt 
    //local coordinates
    DenseMatrix<double> interpolated_dxds(1,spatial_dim,0.0);
    
    //Loop over all parent nodes
    for(unsigned l=0;l<n_node_bulk;l++)
     {
      //Loop over all position types in the bulk
      for(unsigned k=0;k<n_position_type_bulk;k++)
       {
        //Loop over coordinate directions
        for(unsigned i=0;i<spatial_dim;i++)
         {
          //Compute the spatial derivative
          interpolated_dxds(0,i) += 
             Bulk_element_pt->nodal_position_gen(l,k,i)*dpsids(l,k,0);
         }
       }
     }

    //Now the unit normal is just the derivative of the position vector
    //with respect to the single coordinate
    for(unsigned i=0;i<spatial_dim;i++)
     {unit_normal[i] = interpolated_dxds(0,i);}
   }
    break;

    //Line element, derived from a 2D element, in this case
    //the normal is a mess of cross products
    //We need an interior direction, so we must find the local
    //derivatives in the BULK element
   case 1:
   {
    //Find the number of nodes in the bulk element
    const unsigned n_node_bulk = Bulk_element_pt->nnode();
    //Find the number of position types in the bulk element
    const unsigned n_position_type_bulk = 
     Bulk_element_pt->nnodal_position_type();
    
    //Construct the local coordinate in the bulk element
    Vector<double> s_bulk(2);
    //Get the local coordinates in the bulk element
    get_local_coordinate_in_bulk(s,s_bulk);
    
    //Allocate storage for the shape functions and their derivatives wrt
    //local coordinates
    Shape psi(n_node_bulk,n_position_type_bulk);
    DShape dpsids(n_node_bulk,n_position_type_bulk,2);
    //Get the value of the shape functions at the given local coordinate
    Bulk_element_pt->dshape_local(s_bulk,psi,dpsids);
 
    //Calculate all derivatives of the spatial coordinates 
    //wrt local coordinates
    DenseMatrix<double> interpolated_dxds(2,spatial_dim,0.0);
    
    //Loop over all parent nodes
    for(unsigned l=0;l<n_node_bulk;l++)
     {
      //Loop over all position types in the bulk
      for(unsigned k=0;k<n_position_type_bulk;k++)
       {
        //Loop over derivative direction
        for(unsigned j=0;j<2;j++)
         {
          //Loop over coordinate directions
          for(unsigned i=0;i<spatial_dim;i++)
           {
            //Compute the spatial derivative
            interpolated_dxds(j,i) += 
             Bulk_element_pt->nodal_position_gen(l,k,i)*dpsids(l,k,j);
           }
         }
       }
     }
    
    //Initialise the tangent, interior tangent and normal vectors to zero
    //The idea is that even if the element is in a two-dimensional space,
    //the normal cannot be calculated without embedding the element in three
    //dimensions, in which case, the tangent and interior tangent will have
    //zero z-components.
    Vector<double> tangent(3,0.0), interior_tangent(3,0.0), normal(3,0.0);
    
    //We must get the relationship between the coordinate along the face
    //and the local coordinates in the bulk element
    //We must also find an interior direction
    DenseMatrix<double> dsbulk_dsface(2,1,0.0);
    unsigned interior_direction=0;
    get_ds_bulk_ds_face(s,dsbulk_dsface,interior_direction);
    //Load in the values for the tangents
    for(unsigned i=0;i<spatial_dim;i++)
     {
      //Tangent to the face is the derivative wrt to the face coordinate
      //which is calculated using dsbulk_dsface and the chain rule
      tangent[i] = interpolated_dxds(0,i)*dsbulk_dsface(0,0)
       + interpolated_dxds(1,i)*dsbulk_dsface(1,0);
      //Interior tangent to the face is the derivative in the interior 
      //direction
      interior_tangent[i] = interpolated_dxds(interior_direction,i);
     }

    //Now the (3D) normal to the element is the interior tangent 
    //crossed with the tangent
    VectorHelpers::cross(interior_tangent,tangent,normal);
   
    //We find the line normal by crossing the element normal with the tangent
    Vector<double> full_normal(3);
    VectorHelpers::cross(normal,tangent,full_normal);

    //Copy the appropriate entries into the unit normal
    //Two or Three depending upon the spatial dimension of the system
    for(unsigned i=0;i<spatial_dim;i++) {unit_normal[i] = full_normal[i];}
   }
   break;

   //Plane element, derived from 3D element, in this case the normal
   //is just the cross product of the two surface tangents
   //We assume, therefore, that we have three spatial coordinates
   //and two surface coordinates
   //Then we need only to get the derivatives wrt the local coordinates
   //in this face element
   case 2:
   {
#ifdef PARANOID
    //Check that we actually have three spatial dimensions
    if(spatial_dim != 3)
     {
      std::ostringstream error_stream;
      error_stream << "There are only " << spatial_dim
                   << "coordinates at the nodes of this 2D FaceElement,\n"
                   << "which must have come from a 3D Bulk element\n";
      throw OomphLibError(error_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
#endif

    //Find the number of nodes in the element
    const unsigned n_node  = this->nnode();
    //Find the number of position types
    const unsigned n_position_type = this->nnodal_position_type();
    
    //Allocate storage for the shape functions and their derivatives wrt
    //local coordinates
    Shape psi(n_node,n_position_type);
    DShape dpsids(n_node,n_position_type,2);
    //Get the value of the shape functions at the given local coordinate
    this->dshape_local(s,psi,dpsids);
 
    //Calculate all derivatives of the spatial coordinates 
    //wrt local coordinates
    Vector<Vector<double> > interpolated_dxds(2,Vector<double>(3,0));
    
    //Loop over all nodes
    for(unsigned l=0;l<n_node;l++)
     {
      //Loop over all position types
      for(unsigned k=0;k<n_position_type;k++)
       {
        //Loop over derivative directions
        for(unsigned j=0;j<2;j++)
         {
          //Loop over coordinate directions
          for(unsigned i=0;i<3;i++)
           {
            //Compute the spatial derivative
            //Remember that we need to translate the position type
            //to its location in the bulk node
            interpolated_dxds[j][i] += 
             this->nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,j);
           }
         }
       }
     }

    //We now take the cross product of the two normal vectors
    VectorHelpers::cross(interpolated_dxds[0],interpolated_dxds[1],
                         unit_normal);
   }
    break;

   default:

    throw OomphLibError(
     "Cannot have a FaceElement with dimension higher than 2",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
    break;
   }
   
  //Finally normalise unit normal
  double length = VectorHelpers::magnitude(unit_normal);
  for(unsigned i=0;i<spatial_dim;i++) 
   {unit_normal[i] *= Normal_sign/length;}
 }

//=======================================================================
/// Compute the outer unit normal at the ipt-th integration point
//=======================================================================
 void FaceElement::outer_unit_normal(const unsigned &ipt,
                                     Vector<double> &unit_normal) const
 {
  //Find the dimension of the element
  const unsigned element_dim = dim();
  //Find the local coordiantes of the ipt-th integration point
  Vector<double> s(element_dim);
  for(unsigned i=0;i<element_dim;i++) {s[i] = integral_pt()->knot(ipt,i);}
  //Call the outer unit normal function
  outer_unit_normal(s,unit_normal);
 }



//=======================================================================
/// Return vector of local coordinates in bulk element, 
/// given the local coordinates in this FaceElement
//=======================================================================
Vector<double> FaceElement::local_coordinate_in_bulk(
  const Vector<double>& s) const
 {
  //Find the dimension of the bulk element
  unsigned dim_bulk = Bulk_element_pt->dim();

  // Vector of local coordinates in bulk element
  Vector<double> s_bulk(dim_bulk);

  //Use the function pointer if it is set
  if(Face_to_bulk_coordinate_fct_pt)
   {
    //Call the translation function
    (*Face_to_bulk_coordinate_fct_pt)(s,s_bulk);
   }
  else
   {
    throw OomphLibError("Face_to_bulk_coordinate mapping not set",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  
  // Return it
  return s_bulk;
 }



//=======================================================================
/// Calculate the  vector of local coordinates in bulk element, 
/// given the local coordinates in this FaceElement
//=======================================================================
void FaceElement::get_local_coordinate_in_bulk(
  const Vector<double>& s, Vector<double> &s_bulk) const
 {
  //Use the function pointer if it is set
  if(Face_to_bulk_coordinate_fct_pt)
   {
    //Call the translation function
    (*Face_to_bulk_coordinate_fct_pt)(s,s_bulk);
   }
  //Otherwise use the existing (not general) interface
  else
   {
        throw OomphLibError("Face_to_bulk_coordinate mapping not set",
                            OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }


//=======================================================================
///  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 FaceElement::get_ds_bulk_ds_face(const Vector<double> &s,
                                       DenseMatrix<double> &dsbulk_dsface,
                                       unsigned &interior_direction) const
 {
  //Use the function pointer if it is set
  if(Bulk_coordinate_derivatives_fct_pt)
   {
    //Call the translation function
    (*Bulk_coordinate_derivatives_fct_pt)(s,dsbulk_dsface,interior_direction);
   }
  //Otherwise throw an error
  else
   {
    throw OomphLibError(
     "No function for derivatives of bulk coords wrt face coords set",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
 }



///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
//  Functions for elastic general elements
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////



//==================================================================
/// Calculate the L2 norm of the displacement u=R-r to overload the
/// compute_norm function in the GeneralisedElement base class
//==================================================================
 void SolidFiniteElement::compute_norm(double& el_norm)
 {
  // Initialise el_norm to 0.0
  el_norm=0.0;
  
  unsigned n_dim=dim();

  // Vector of local coordinates
  Vector<double> s(n_dim);
  
  // Displacement vector, Lagrangian coordinate vector and Eulerian
  // coordinate vector respectively
  Vector<double> disp(n_dim,0.0);
  Vector<double> xi(n_dim,0.0);
  Vector<double> x(n_dim,0.0);
  
  // Find out how many nodes there are in the element
  unsigned n_node=this->nnode();
  
  // Construct a shape function
  Shape psi(n_node);
  
  // Get the number of integration points
  unsigned n_intpt=this->integral_pt()->nweight();
  
  // Loop over the integration points
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {     
    // Assign values of s
    for(unsigned i=0;i<n_dim;i++)
     {
      s[i]=this->integral_pt()->knot(ipt,i);
     }
    
    // Get the integral weight
    double w=this->integral_pt()->weight(ipt);
    
    // Get jacobian of mapping
    double J=this->J_eulerian(s);
    
    // Premultiply the weights and the Jacobian
    double W=w*J;
    
    // Get the Lagrangian and Eulerian coordinate at this point, respectively
    this->interpolated_xi(s,xi);
    this->interpolated_x(s,x);
    
    // Calculate the displacement vector, u=R-r=x-xi
    for (unsigned idim=0;idim<n_dim;idim++)
     {
      disp[idim]=x[idim]-xi[idim];
     }
    
    // Add to norm
    for (unsigned ii=0;ii<n_dim;ii++)
     {
      el_norm+=(disp[ii]*disp[ii])*W;
     }
   }    
 } // End of compute_norm(...)



//=========================================================================
/// \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 SolidFiniteElement::
describe_local_dofs(std::ostream& out,const std::string& current_string) const
{
 // Call the standard finite element description.
 FiniteElement::describe_local_dofs(out,current_string);
 describe_solid_local_dofs(out,current_string);
}


//=========================================================================
/// Internal function that is used to assemble the jacobian of the mapping
/// from local coordinates (s) to the lagrangian coordinates (xi), given the
/// derivatives of the shape functions. 
//=========================================================================
 void SolidFiniteElement::
 assemble_local_to_lagrangian_jacobian(const DShape &dpsids,
                                       DenseMatrix<double> &jacobian) const
 {
  //Find the the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  //We shall ASSUME (ENFORCE) that Lagrangian coordinates must 
  //be interpolated through the nodes
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_lagrangian_type();

#ifdef PARANOID
 //Check for dimensional compatibility
 if(el_dim != Lagrangian_dimension)
  {
   std::ostringstream error_message;
   error_message << "Dimension mismatch" << std::endl;
   error_message << "The elemental dimension: " << el_dim
                 << " must equal the nodal Lagrangian dimension: " 
                 << Lagrangian_dimension 
                 << " for the jacobian of the mapping to be well-defined"
                 << std::endl;
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

  //Loop over the rows of the jacobian
  for(unsigned i=0;i<el_dim;i++)
   {
    //Loop over the columns of the jacobian
    for(unsigned j=0;j<el_dim;j++)
     {
      //Zero the entry
      jacobian(i,j) = 0.0;
      //Loop over the shape functions
      for(unsigned l=0;l<n_shape;l++)
       {
        for(unsigned k=0;k<n_shape_type;k++)
         {
          //Jacobian is dx_j/ds_i, which is represented by the sum
          //over the dpsi/ds_i of the nodal points X j
          //Call the Non-hanging version of positions
          //This is overloaded in refineable elements
          jacobian(i,j) += raw_lagrangian_position_gen(l,k,j)*dpsids(l,k,i);
         }
       }
     }
   }
 }

//=========================================================================
/// Internal function that is used to assemble the jacobian of second
/// derivatives of the the mapping from local coordinates (s) to the 
/// lagrangian coordinates (xi), given the second derivatives of the 
/// shape functions. 
//=========================================================================
 void SolidFiniteElement::
 assemble_local_to_lagrangian_jacobian2(const DShape &d2psids,
                                        DenseMatrix<double> &jacobian2) const
 {
  //Find the the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  //We ENFORCE that Lagrangian coordinates must be interpolated through
  //the nodes
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_lagrangian_type();
  //Find the number of second derivatives
  const unsigned n_row = N2deriv[el_dim];
 
  //Assemble the "jacobian" (d^2 x_j/ds_i^2) for second derivatives of 
  //shape functions
  //Loop over the rows (number of second derivatives)
  for(unsigned i=0;i<n_row;i++)
   {
    //Loop over the columns (element dimension
    for(unsigned j=0;j<el_dim;j++)
     {
      //Zero the entry
      jacobian2(i,j) = 0.0;
      //Loop over the shape functions
      for(unsigned l=0;l<n_shape;l++)
       {
        //Loop over the shape function types
        for(unsigned k=0;k<n_shape_type;k++)
         {
          //Add the terms to the jacobian entry
          //Call the Non-hanging version of positions
          //This is overloaded in refineable elements
          jacobian2(i,j) += raw_lagrangian_position_gen(l,k,j)*d2psids(l,k,i);
         }
       }
     }
   }
 }

//============================================================================
/// \short Destructor for SolidFiniteElement:
//============================================================================
 SolidFiniteElement:: ~SolidFiniteElement() 
 {
  //Delete the storage allocated for the positional local equations
  delete[] Position_local_eqn;
 }


//==========================================================================
/// Calculate the mapping from local to lagrangian 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.
/// The local derivatives are passed as dpsids and the jacobian and 
/// inverse jacobian are returned.
//==========================================================================
 double SolidFiniteElement::
 local_to_lagrangian_mapping_diagonal(const DShape &dpsids,
                                      DenseMatrix<double> &jacobian,
                                      DenseMatrix<double> &inverse_jacobian)
  const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
 //We shall ASSUME (ENFORCE) that Lagrangian coordinates must 
  //be interpolated through the nodes
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_lagrangian_type();
 
  //In this case we assume that there are no cross terms, that is
  //global coordinate 0 is always in the direction of local coordinate 0

  //Loop over the coordinates
  for(unsigned i=0;i<el_dim;i++)
   {
    //Zero the jacobian and inverse jacobian entries
    for(unsigned j=0;j<el_dim;j++) 
     {jacobian(i,j) = 0.0; inverse_jacobian(i,j) = 0.0;}
   
    //Loop over the shape functions
    for(unsigned l=0;l<n_shape;l++)
     {
      //Loop over the types of dof
      for(unsigned k=0;k<n_shape_type;k++)
       {
        //Derivatives are always dx_{i}/ds_{i}
        jacobian(i,i) += raw_lagrangian_position_gen(l,k,i)*dpsids(l,k,i);
       }
     }
   }
 
  //Now calculate the determinant of the matrix
  double det = 1.0;
  for(unsigned i=0;i<el_dim;i++) {det *= jacobian(i,i);}
 
//Report if Matrix is singular, or negative
#ifdef PARANOID
  check_jacobian(det);
#endif
 
  //Calculate the inverse mapping (trivial in this case)
  for(unsigned i=0;i<el_dim;i++) 
   {inverse_jacobian(i,i) = 1.0/jacobian(i,i);}
 
  //Return the value of the Jacobian
  return(det);
 }

//========================================================================
/// 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.
/// General case, may be overloaded
//========================================================================
 double SolidFiniteElement::
 dshape_lagrangian(const Vector<double> &s, Shape &psi, DShape &dpsi) const
 {
  //Find the element dimension
  const unsigned el_dim = dim();
 
  //Get the values of the shape function and local derivatives
  //Temporarily stored in dpsi
  dshape_local(s,psi,dpsi);
 
  //Allocate memory for the inverse jacobian
  DenseMatrix<double> inverse_jacobian(el_dim);
  //Now calculate the inverse jacobian
  const double det = local_to_lagrangian_mapping(dpsi,inverse_jacobian);
 
  //Now set the values of the derivatives to be dpsidxi
  transform_derivatives(inverse_jacobian,dpsi);
  //Return the determinant of the jacobian
  return det;
 }

//========================================================================= 
/// Compute the geometric shape functions and also first
/// derivatives w.r.t. Lagrangian coordinates at integration point ipt.
/// Most general form of function, but may be over-loaded if desired
//========================================================================
 double SolidFiniteElement::dshape_lagrangian_at_knot(const unsigned &ipt,
                                                      Shape &psi, 
                                                      DShape &dpsi) const
 {
  //Find the element dimension
  const unsigned el_dim = dim();
 
  //Shape function for the local derivatives
  //Again we ASSUME (insist) that the lagrangian coordinates
  //are interpolated through the nodes
  //Get the values of the shape function and local derivatives
  dshape_local_at_knot(ipt,psi,dpsi);
 
  //Allocate memory for the inverse jacobian
  DenseMatrix<double> inverse_jacobian(el_dim);
  //Now calculate the inverse jacobian
  const double det = local_to_lagrangian_mapping(dpsi,inverse_jacobian);

  //Now set the values of the derivatives
  transform_derivatives(inverse_jacobian,dpsi);
  //Return the determinant of the jacobian
  return det;
 }

//========================================================================
/// \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 SolidFiniteElement::
 d2shape_lagrangian(const Vector<double> &s, Shape &psi, 
                    DShape &dpsi, DShape &d2psi) const
 {
  //Find the element dimension
  const unsigned el_dim = dim();
  //Find the number of second derivatives required
  const unsigned n_deriv = N2deriv[el_dim];

  //Get the values of the shape function and local derivatives
  d2shape_local(s,psi,dpsi,d2psi);
 
  //Allocate memory for the jacobian and inverse jacobian 
  DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);
  //Calculate the jacobian and inverse jacobian
  const double det = 
   local_to_lagrangian_mapping(dpsi,jacobian,inverse_jacobian);

  //Allocate memory for the jacobian of second derivatives
  DenseMatrix<double> jacobian2(n_deriv,el_dim);
  //Assemble the jacobian of second derivatives
  assemble_local_to_lagrangian_jacobian2(d2psi,jacobian2);

  //Now set the value of the derivatives
  transform_second_derivatives(jacobian,inverse_jacobian,
                               jacobian2,dpsi,d2psi);
  //Return the determinant of the mapping
  return det;
 }

//========================================================================== 
/// \short Compute 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 \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 SolidFiniteElement::d2shape_lagrangian_at_knot(const unsigned &ipt, 
                                                       Shape &psi, 
                                                       DShape &dpsi, 
                                                       DShape &d2psi) const
 {
  //Find the values of the indices of the shape functions
  //Find the element dimension
  const unsigned el_dim = dim();
  //Find the number of second derivatives required
  const unsigned n_deriv = N2deriv[el_dim];
 
  //Get the values of the shape function and local derivatives
  d2shape_local_at_knot(ipt,psi,dpsi,d2psi);
 
  //Allocate memory for the jacobian and inverse jacobian 
  DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);
  //Calculate the jacobian and inverse jacobian
  const double det = 
   local_to_lagrangian_mapping(dpsi,jacobian,inverse_jacobian);

  //Allocate memory for the jacobian of second derivatives
  DenseMatrix<double> jacobian2(n_deriv,el_dim);
  //Assemble the jacobian of second derivatives
  assemble_local_to_lagrangian_jacobian2(d2psi,jacobian2);

  //Now set the value of the derivatives
  transform_second_derivatives(jacobian,inverse_jacobian,
                               jacobian2,dpsi,d2psi);
  //Return the determinant of the mapping
  return det;
 }

//============================================================================
/// \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 SolidFiniteElement::
describe_solid_local_dofs(std::ostream& out,
                          const std::string& current_string) const
{
 //Find the number of nodes
 const unsigned n_node = nnode();
 //Loop over the nodes
 for(unsigned n=0;n<n_node;n++)
  {
   //Cast to a solid node
   SolidNode* cast_node_pt = static_cast<SolidNode*>(node_pt(n));
   std::stringstream conversion;
   conversion<<" of Solid Node "<<n<<current_string;
   std::string in(conversion.str());
   cast_node_pt->describe_dofs(out,in);
  }
}

//============================================================================
/// Assign local equation numbers for the solid equations in the element.
//  This can be done at a high level assuming, as I am, that the equations will
//  always be formulated in terms of nodal positions. 
/// If the boolean argument is true then pointers to the dofs will be 
/// stored in Dof_pt.
//============================================================================
 void SolidFiniteElement::assign_solid_local_eqn_numbers(
  const bool &store_local_dof_pt)
 {
  //Find the number of nodes
  const unsigned n_node = nnode();
 
  //Check there are nodes!
  if(n_node > 0)
   {
    //Find the number of position types and dimensions of the nodes
    //Local caching
    const unsigned n_position_type = nnodal_position_type();
    const unsigned nodal_dim = nodal_dimension();
     
    //Delete the existing storage
    delete[] Position_local_eqn;
    //Resize the storage for the positional equation numbers
    Position_local_eqn = new int[n_node*n_position_type*nodal_dim];
    
    //A local queue to store the global equation numbers
    std::deque<unsigned long> global_eqn_number_queue;

    //Get the number of dofs so far, this must be outside both loops
    //so that both can use it
    unsigned local_eqn_number = ndof();
    
    //Loop over the nodes
    for(unsigned n=0;n<n_node;n++)
     {
      //Cast to a solid node
      SolidNode* cast_node_pt = static_cast<SolidNode*>(node_pt(n));

      //Loop over the number of position dofs
      for(unsigned j=0;j<n_position_type;j++)
       {
        //Loop over the dimension of each node
        for(unsigned k=0;k<nodal_dim;k++)
         {
          //Get equation number
          //Note eqn_number is long !
          long eqn_number = cast_node_pt->position_eqn_number(j,k);
          //If equation_number positive add to array
          if(eqn_number >= 0)
           {
            //Add to global array
            global_eqn_number_queue.push_back(eqn_number);
            //Add pointer to the dof to the queue if required
            if(store_local_dof_pt)
             {
              GeneralisedElement::Dof_pt_deque.push_back(
               &(cast_node_pt->x_gen(j,k)));
             }

            //Add to look-up scheme
            Position_local_eqn[(n*n_position_type + j)*nodal_dim + k] = 
             local_eqn_number;
            //Increment the local equation number
            local_eqn_number++;
           }
          else
           {
            Position_local_eqn[(n*n_position_type + j)*nodal_dim + k] = 
             Data::Is_pinned;
           }
         }
       }
     } //End of loop over nodes

    //Now add our global equations numbers to the internal element storage
    add_global_eqn_numbers(global_eqn_number_queue,
                           GeneralisedElement::Dof_pt_deque);
    //Clear the memory used in the deque
    if(store_local_dof_pt)
     {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}

   } //End of the case when there are nodes
 }


//============================================================================
/// This function calculates the entries of Jacobian matrix, used in 
/// the Newton method, associated with the elastic problem in which the
/// nodal position is a variable. It does this using finite differences, 
/// rather than an analytical formulation, so can be done in total generality.
//==========================================================================
 void SolidFiniteElement::
 fill_in_jacobian_from_solid_position_by_fd(Vector<double> &residuals,
                                            DenseMatrix<double> &jacobian)
 {

  // Flag to indicate if we use first or second order FD
  // bool use_first_order_fd=false;

  //Find the number of nodes
  const unsigned n_node = nnode();

  //If there aren't any nodes, then return straight away
  if(n_node == 0) {return;}

  //Call the update function to ensure that the element is in
  //a consistent state before finite differencing starts
  update_before_solid_position_fd();

  //Get the number of position dofs and dimensions at the node
  const unsigned n_position_type = nnodal_position_type();
  const unsigned nodal_dim = nodal_dimension();

  //Find the number of dofs in the element
  const unsigned n_dof = this->ndof();

  //Create newres vectors
  Vector<double> newres(n_dof);
  // Vector<double> newres_minus(n_dof);
  
  //Integer storage for local unknown
  int local_unknown=0;
 
  //Should probably give this a more global scope
  const double fd_step = Default_fd_jacobian_step;
 
  //Loop over the nodes
  for(unsigned n=0;n<n_node;n++)
   {
    //Loop over position dofs
    for(unsigned k=0;k<n_position_type;k++)
     {
      //Loop over dimension
      for(unsigned i=0;i<nodal_dim;i++)
       {
        //If the variable is free
        local_unknown = position_local_eqn(n,k,i);
        if(local_unknown >= 0)
         {
          //Store a pointer to the (generalised) Eulerian nodal position
          double* const value_pt = &(node_pt(n)->x_gen(k,i));

          //Save the old value of the (generalised) Eulerian nodal position
          const double old_var = *value_pt;
           
          //Increment the (generalised) Eulerian nodal position
          *value_pt += fd_step;
          
          // Perform any auxialiary node updates
          node_pt(n)->perform_auxiliary_node_update_fct();

          // Update any other slaved variables
          update_in_solid_position_fd(n);

          //Calculate the new residuals
          get_residuals(newres);
         
//          if (use_first_order_fd)
           {
            //Do forward finite differences
            for(unsigned m=0;m<n_dof;m++)
             {
              //Stick the entry into the Jacobian matrix
              jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
             }
           }
//           else
//            {
//             //Take backwards step for the  (generalised) Eulerian nodal 
//             // position
//             node_pt(n)->x_gen(k,i) = old_var-fd_step;
           
//             //Calculate the new residuals at backward position
//             get_residuals(newres_minus);

//             //Do central finite differences
//             for(unsigned m=0;m<n_dof;m++)
//              {
//               //Stick the entry into the Jacobian matrix
//               jacobian(m,local_unknown) = 
//                (newres[m] - newres_minus[m])/(2.0*fd_step);
//              }
//            }

          //Reset the (generalised) Eulerian nodal position
          *value_pt = old_var;


          // Perform any auxialiary node updates
          node_pt(n)->perform_auxiliary_node_update_fct();
          
          // Reset any other slaved variables
          reset_in_solid_position_fd(n);
         }
       }
     }
   }

  //End of finite difference loop
  //Final reset of any slaved data
  reset_after_solid_position_fd();
 }

//============================================================================
/// Return i-th FE-interpolated Lagrangian coordinate at
/// local coordinate s.
//============================================================================
 double SolidFiniteElement::interpolated_xi(const Vector<double> &s, 
                                            const unsigned &i) const
 { 
  //Find the number of nodes
  const unsigned n_node = nnode();

  //Find the number of lagrangian types from the first node
  const unsigned n_lagrangian_type = nnodal_lagrangian_type();

  //Assign storage for the local shape function
  Shape psi(n_node,n_lagrangian_type);

  //Find the values of shape function
  shape(s,psi);

  //Initialise value of xi
  double interpolated_xi = 0.0;

  //Loop over the local nodes
  for(unsigned l=0;l<n_node;l++) 
   {
    //Loop over the number of dofs
    for(unsigned k=0;k<n_lagrangian_type;k++)
     {
      interpolated_xi += lagrangian_position_gen(l,k,i)*psi(l,k);
     }
   }
   
  return(interpolated_xi);
 }

//============================================================================
/// Compute FE-interpolated Lagrangian coordinate vector xi[] at
/// local coordinate s.
//============================================================================
 void SolidFiniteElement::interpolated_xi(const Vector<double> &s, 
                                          Vector<double> &xi) const
 { 
  //Find the number of nodes
  const unsigned n_node = nnode();

  //Find the number of lagrangian types from the first node
  const unsigned n_lagrangian_type = nnodal_lagrangian_type();

  //Assign storage for the local shape function
  Shape psi(n_node,n_lagrangian_type);

  //Find the values of shape function
  shape(s,psi);

  //Read out the number of lagrangian coordinates from the node
  const unsigned n_lagrangian = lagrangian_dimension();

  //Loop over the number of lagrangian coordinates
  for(unsigned i=0;i<n_lagrangian;i++)
   {
    //Initialise component to zero
    xi[i] = 0.0;

    //Loop over the local nodes
    for(unsigned l=0;l<n_node;l++) 
     {
      //Loop over the number of dofs
      for(unsigned k=0;k<n_lagrangian_type;k++)
       {
        xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
       }
     }
   }
 }





//============================================================================
/// Compute derivatives of FE-interpolated Lagrangian coordinates xi
/// with respect to local coordinates: dxids[i][j]=dxi_i/ds_j
//============================================================================
void SolidFiniteElement::interpolated_dxids(const Vector<double> &s, 
                                            DenseMatrix<double> &dxids) const
{ 
 //Find the number of nodes
 const unsigned n_node = nnode();
 
 //Find the number of lagrangian types from the first node
 const unsigned n_lagrangian_type = nnodal_lagrangian_type();
 
 // Dimension of the element =  number of local coordinates
 unsigned el_dim=dim();
 
 //Assign storage for the local shape function
 Shape psi(n_node,n_lagrangian_type);
 DShape dpsi(n_node,n_lagrangian_type,el_dim);
 
 //Find the values of shape function and its derivs w.r.t. to local coords
 dshape_local(s,psi,dpsi);
 
 //Read out the number of lagrangian coordinates from the node
 const unsigned n_lagrangian = lagrangian_dimension();
 
 //Loop over the number of lagrangian and local coordinates 
 for(unsigned i=0;i<n_lagrangian;i++)
  {
   for(unsigned j=0;j<el_dim;j++)
    {
     //Initialise component to zero
     dxids(i,j) = 0.0;
     
     //Loop over the local nodes
     for(unsigned l=0;l<n_node;l++) 
      {
       //Loop over the number of dofs
       for(unsigned k=0;k<n_lagrangian_type;k++)
        {
         dxids(i,j) += lagrangian_position_gen(l,k,i)*dpsi(l,k,j);
        }
      }
    }
  }
}

//=======================================================================
/// Add jacobian and residuals for consistent assignment of 
/// initial "accelerations" in Newmark scheme. Jacobian is the mass matrix.
//=======================================================================
void SolidFiniteElement::
fill_in_jacobian_for_newmark_accel(DenseMatrix<double> &jacobian)
{
 
#ifdef PARANOID
 // Check that we're computing the real residuals, not the
 // ones corresponding to the assignement of a prescribed displacement
 if ((Solid_ic_pt!=0)||(!Solve_for_consistent_newmark_accel_flag))
  {
   std::ostringstream error_stream;
   error_stream << "Can only call fill_in_jacobian_for_newmark_accel()\n" 
                << "With Solve_for_consistent_newmark_accel_flag:" 
                << Solve_for_consistent_newmark_accel_flag << std::endl;
   error_stream << "Solid_ic_pt: " << Solid_ic_pt << std::endl;

   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 //Find the number of nodes
 const unsigned n_node = nnode();
 const unsigned n_position_type = nnodal_position_type();
 const unsigned nodal_dim = nodal_dimension();

 //Set the number of Lagrangian coordinates
 const unsigned n_lagrangian = dim();
  
 //Integer storage for local equation and unknown
 int local_eqn=0, local_unknown=0;

 // Set up memory for the shape functions:
 // # of nodes, # of positional dofs
 Shape psi(n_node,n_position_type);

 // # of nodes, # of positional dofs, # of lagrangian coords (for deriv)
 // Not needed but they come for free when compute the Jacobian.
 DShape dpsidxi(n_node,n_position_type,n_lagrangian);

 //Set # of integration points
 const unsigned n_intpt = integral_pt()->nweight();

 //Set the Vector to hold local coordinates
 Vector<double> s(n_lagrangian);

 // Get positional timestepper from first nodal point
 TimeStepper* timestepper_pt= node_pt(0)->position_time_stepper_pt();

#ifdef PARANOID
 // Of course all this only works if we're actually using a
 // Newmark timestepper!
 if (timestepper_pt->type()!="Newmark")
  {
   std::ostringstream error_stream;
   error_stream
    << "Assignment of Newmark accelerations obviously only works\n"
    << "for Newmark timesteppers. You've called me with " 
    << timestepper_pt->type() << std::endl;

   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // "Acceleration" is the last history value:
 const unsigned ntstorage=timestepper_pt->ntstorage();
 const double accel_weight=timestepper_pt->weight(2,ntstorage-1);

 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
     
   //Assign values of s
   for(unsigned i=0;i<n_lagrangian;i++) {s[i] = integral_pt()->knot(ipt,i);}
   
   //Get the integral weight
   double w = integral_pt()->weight(ipt);

   // Jacobian of mapping between local and Lagrangian coords and shape
   // functions
   double J = dshape_lagrangian(s,psi,dpsidxi);

   // Get Lagrangian coordinate
   Vector<double> xi(n_lagrangian);
   interpolated_xi(s,xi);

   // Get multiplier for inertia terms
   double factor=multiplier(xi);

   //Premultiply the weights and the Jacobian
   double W = w*J;

   //Loop over the nodes
   for(unsigned l0=0;l0<n_node;l0++)
    {
     //Loop over the positional dofs: 'Type': 0: displacement; 1: deriv
     for(unsigned k0=0;k0<n_position_type;k0++)
      {
       // Loop over Eulerian directions
       for(unsigned i0=0;i0<nodal_dim;i0++)
        {
         local_eqn = position_local_eqn(l0,k0,i0);
         //If it's a degree of freedom, add contribution
         if(local_eqn >= 0)
          {
           //Loop over the nodes
           for(unsigned l1=0;l1<n_node;l1++)
            {
             //Loop over the positional dofs: 'Type': 0: displacement; 
             //1: deriv
             for(unsigned k1=0;k1<n_position_type;k1++)
              {
               // Loop over Eulerian directions
               unsigned i1=i0;
               {
                local_unknown = position_local_eqn(l1,k1,i1);
                //If it's a degree of freedom, add contribution
                if(local_unknown >= 0)
                 {
                  // Add contribution: Mass matrix, multiplied by
                  // weight for "acceleration" in Newmark scheme
                  // and general multiplier 
                  jacobian(local_eqn,local_unknown) +=
                   factor*accel_weight*psi(l0,k0)*psi(l1,k1)*W;
                 }
               }
              }
            }
          }
        }
      }
    }
   
  } //End of loop over the integration points

}





//=======================================================================
/// 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 SolidFiniteElement::fill_in_generic_jacobian_for_solid_ic(
 Vector<double> &residuals,
 DenseMatrix<double> &jacobian,
 const unsigned& flag)
{
 //Find the number of nodes and position types
 const unsigned n_node = nnode();
 const unsigned n_position_type = nnodal_position_type();

 //Set the dimension of the global coordinates
 const unsigned nodal_dim = nodal_dimension();
 
 //Find the number of lagragian types from the first node
 const unsigned n_lagrangian_type = nnodal_lagrangian_type();

 //Set the number of lagrangian coordinates
 const unsigned n_lagrangian = dim();
  
 //Integer storage for local equation number
 int local_eqn=0;
 int local_unknown=0;

 // # of nodes, # of positional dofs
 Shape psi(n_node,n_position_type);
 
 // # of nodes, # of positional dofs, # of lagrangian coords (for deriv)
 // not needed but they come for free when we compute the Jacobian
 //DShape dpsidxi(n_node,n_position_type,n_lagrangian);
 
 //Set # of integration points
 const unsigned n_intpt = integral_pt()->nweight();
 
 //Set the Vector to hold local coordinates
 Vector<double> s(n_lagrangian);
 
 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {

   //Assign values of s
   for(unsigned i=0;i<n_lagrangian;i++) 
    {s[i] = integral_pt()->knot(ipt,i);}

   //Get the integral weight
   double w = integral_pt()->weight(ipt);

   // Shape fcts
   shape(s,psi);

   // Get Lagrangian coordinate
   Vector<double> xi(n_lagrangian,0.0);

   //Loop over the number of lagrangian coordinates
   for(unsigned i=0;i<n_lagrangian;i++)
    {
     //Loop over the local nodes
     for(unsigned l=0;l<n_node;l++) 
      {
       //Loop over the number of dofs
       for(unsigned k=0;k<n_lagrangian_type;k++)
        {
         xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
        }
      }
    }

   // Get initial condition
   Vector<double> drdt_ic(nodal_dim);
   Solid_ic_pt->geom_object_pt()->
    dposition_dt(xi,Solid_ic_pt->ic_time_deriv(),drdt_ic);
   
   // Weak form of assignment of initial guess
   
   //Loop over the number of node
   for (unsigned l=0;l<n_node;l++)
    {
     //Loop over the type of degree of freedom
     for (unsigned k=0;k<n_position_type;k++)
      {
       //Loop over the coordinate directions
       for (unsigned i=0;i<nodal_dim;i++)
        {
         local_eqn = position_local_eqn(l,k,i);

         //If it's not a boundary condition
         if(local_eqn >= 0)
          {
           // Note we're ignoring the mapping between local and
           // global Lagrangian coords -- doesn't matter here;
           // we're just enforcing a slightly different
           // weak form.
           residuals[local_eqn] += (interpolated_x(s,i)-drdt_ic[i])*
            psi(l,k)*w;


           // Do Jacobian too?
           if (flag==1)
            {     

             //Loop over the number of node
             for (unsigned ll=0;ll<n_node;ll++)
              {                
               //Loop over the type of degree of freedom
               for (unsigned kk=0;kk<n_position_type;kk++)
                {    

                 // Only diagonal term
                 unsigned ii=i;
                                 
                 local_unknown = position_local_eqn(ll,kk,ii);
                 
                 //If it's not a boundary condition
                 if(local_unknown >= 0)
                  {
                   // Note we're ignoring the mapping between local and
                   // global Lagrangian coords -- doesn't matter here;
                   // we're just enforcing a slightly different
                   // weak form.
                   jacobian(local_eqn,local_unknown) +=
                    psi(ll,kk)*psi(l,k)*w;
                  }
                 else
                  {
                   oomph_info 
                    << "WARNING: You should really free all Data" 
                    << std::endl
                    << "         before setup of initial guess" 
                    << std::endl
                    << "ll, kk, ii " << ll << " " << kk << " " << ii 
                    << std::endl;
                  }
                }
              }
             
            }
          }
         else
          {
           oomph_info 
            << "WARNING: You should really free all Data" << std::endl
            << "         before setup of initial guess" << std::endl
            << "l, k, i " << l << " " << k << " " << i << std::endl;
          }
        }
      }
    }
  
  } //End of loop over the integration points
}


//===============================================================
/// Return the geometric shape function at the local coordinate s
//===============================================================
 void PointElement::shape(const Vector<double> &s, Shape &psi) const
  {
   // Single shape function always has value 1
   psi[0]=1.0;
  }
 
//=======================================================================
/// Assign the static Default_integration_scheme
//=======================================================================
PointIntegral PointElement::Default_integration_scheme;


}