hermite_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 functions for the QHermiteElement classes

//oomph-lib header
#include "hermite_elements.h"
#include "Qelement_face_coordinate_translation_schemes.h"

namespace oomph
{


////////////////////////////////////////////////////////////////
//       1D Hermite elements
////////////////////////////////////////////////////////////////

//=======================================================================
/// Assign the static Default_integration_scheme
//=======================================================================
template<unsigned DIM>
Gauss<DIM,3> QHermiteElement<DIM>::Default_integration_scheme;

//=======================================================================
/// Shape function for specific QHermiteElement<1>
//=======================================================================
template <>
void QHermiteElement<1>::shape(const Vector<double> &s, Shape &psi) 
 const
{
 //Local storage
 double Psi[2][2];
 //Call the OneDimensional Shape functions
 OneDimHermite::shape(s[0],Psi);

 //Loop over the number of nodes
 for(unsigned l=0;l<2;l++)
  {
   //Loop over the number of dofs
   for(unsigned k=0;k<2;k++)
    {
     psi(l,k) = Psi[l][k];
    }
  }
}

//=======================================================================
/// Derivatives of shape functions for specific  QHermiteElement<1>
//=======================================================================
template <>
void QHermiteElement<1>::dshape_local(const Vector<double> &s, 
                                      Shape &psi, 
                                      DShape &dpsids) const
{
 //Local storage
 double Psi[2][2], DPsi[2][2];
 //Call the OneDimensional Shape functions
 OneDimHermite::shape(s[0],Psi);
 OneDimHermite::dshape(s[0],DPsi);

 //Loop over number of nodes
 for(unsigned l=0;l<2;l++) 
  {
   //Loop over number of dofs
   for(unsigned k=0;k<2;k++)
    {
     psi(l,k) = Psi[l][k];
     dpsids(l,k,0) = DPsi[l][k];
    }
  }
}

//=======================================================================
/// Derivatives and second derivatives of shape functions for specific  
/// QHermiteElement<1> 
/// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
//=======================================================================
template<>
void QHermiteElement<1>::d2shape_local(const Vector<double> &s, 
                                       Shape &psi, 
                                       DShape &dpsids, 
                                       DShape &d2psids) const
{
 //Local storage
 double Psi[2][2], DPsi[2][2], D2Psi[2][2];
 //Call the OneDimensional Shape functions
 OneDimHermite::shape(s[0],Psi);
 OneDimHermite::dshape(s[0],DPsi);
 OneDimHermite::d2shape(s[0],D2Psi);
 
 //Loop over number of nodes
 for(unsigned l=0;l<2;l++) 
  {
   //Loop over number of dofs
   for(unsigned k=0;k<2;k++)
    {
     psi(l,k) = Psi[l][k];
     dpsids(l,k,0) = DPsi[l][k];
     d2psids(l,k,0) = D2Psi[l][k];
    }
  }
}



//=======================================================================
/// The output function for general 1D QHermiteElements
//=======================================================================
template <>
void QHermiteElement<1>::output(std::ostream &outfile)
{
 //Tecplot header info 
 outfile << "ZONE I=" << 2 << std::endl;

 //Find the dimension of the node
 unsigned n_dim = this->nodal_dimension();
   
 //Loop over element nodes
 for(unsigned l=0;l<2;l++)
  {
   //Loop over the dimensions and output the position
   for(unsigned i=0;i<n_dim;i++)
    {
     outfile << node_pt(l)->x(i) << " ";
    }

   //Find the number of types of dof stored at each node
   unsigned n_position_type = node_pt(l)->nposition_type();
   //Loop over the additional positional dofs
   for(unsigned k=1;k<n_position_type;k++)
    {
     for(unsigned i=0;i<n_dim;i++)
      {
       outfile << node_pt(l)->x_gen(k,i) << " ";
      }
    }
   
   //Find out how many data values at the node
   unsigned initial_nvalue = node_pt(l)->nvalue();
   //Lopp over the data and output whether pinned or not
   for(unsigned i=0;i<initial_nvalue;i++)
    {
     outfile << node_pt(l)->is_pinned(i) << " ";
    }
   outfile << std::endl;
  }
  outfile << std::endl;
}

//=======================================================================
/// The output function for n_plot points in each coordinate direction
//=======================================================================
template <>
void QHermiteElement<1>::output(std::ostream &outfile, const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(1);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << std::endl;

 //Find the dimension of the first node
 unsigned n_dim = this->nodal_dimension();

 //Loop over plot points
 for(unsigned l=0;l<n_plot;l++)
  {
   s[0] = -1.0 + l*2.0/(n_plot-1);

   //Output the coordinates
   for (unsigned i=0;i<n_dim;i++)
    {
     outfile << interpolated_x(s,i) << " " ;
    }
   outfile <<  std::endl;
  }
  outfile << std::endl;
}




//=======================================================================
/// The C-style output function for general 1D QHermiteElements
//=======================================================================
template <>
void QHermiteElement<1>::output(FILE* file_pt)
{
 //Tecplot header info 
 fprintf(file_pt,"ZONE I=2\n");

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();
   
 //Loop over element nodes
 for(unsigned l=0;l<2;l++)
  {
   //Loop over the dimensions and output the position
   for(unsigned i=0;i<n_dim;i++)
    {
     fprintf(file_pt,"%g ",node_pt(l)->x(i));
    }

   //Find the number of types of dof stored at each node
   unsigned n_position_type = node_pt(l)->nposition_type();
   //Loop over the additional positional dofs
   for(unsigned k=1;k<n_position_type;k++)
    {
     for(unsigned i=0;i<n_dim;i++)
      {
       fprintf(file_pt,"%g ",node_pt(l)->x_gen(k,i));
      }
    }
   
   //Find out how many data values at the node
   unsigned initial_nvalue = node_pt(l)->nvalue();
   //Lopp over the data and output whether pinned or not
   for(unsigned i=0;i<initial_nvalue;i++)
    {
     fprintf(file_pt,"%i ",node_pt(l)->is_pinned(i));
    }
   fprintf(file_pt,"\n");
  }
 fprintf(file_pt,"\n");
}


//=======================================================================
/// The C-style output function for n_plot points in each coordinate direction
//=======================================================================
template <>
void QHermiteElement<1>::output(FILE* file_pt, const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(1);

 //Tecplot header info 
 fprintf(file_pt,"ZONE I=%i \n",n_plot);

 //Find the dimension of the first node
 unsigned n_dim = this->nodal_dimension();

 //Loop over element nodes
 for(unsigned l=0;l<n_plot;l++)
  {
   s[0] = -1.0 + l*2.0/(n_plot-1);
   //Output the coordinates
   for (unsigned i=0;i<n_dim;i++)
    {
     fprintf(file_pt,"%g ",interpolated_x(s,i));
    }
   fprintf(file_pt,"\n");
  }
}



//===========================================================
/// Function to setup geometrical information for lower-dimensional 
/// FaceElements (single node elements)
//===========================================================
template<>
void QHermiteElement<1>::build_face_element(const int &face_index,
                                            FaceElement *face_element_pt)
{
 // Set the nodal dimension from the "bulk"
 face_element_pt->set_nodal_dimension(node_pt(0)->ndim());

 // Set the pointer to the "bulk" element
 face_element_pt->bulk_element_pt()=this;
 
#ifdef OOMPH_HAS_MPI
 // Pass on non-halo ID
 face_element_pt->set_halo(Non_halo_proc_ID);
#endif

 // Resize the storage for the original number of values at the (one and only)
 // node of the face element.
 face_element_pt->nbulk_value_resize(1);

 // Resize the storage for the bulk node number corresponding to the (one
 // and only) node of the face element
 face_element_pt->bulk_node_number_resize(1);

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

 //Now set up the node pointer
 //The convention is that the "normal", the coordinate, should always point
 //out of the element
 switch(face_index)
  {
   //Bottom, normal sign is negative (coordinate points into element)
  case(-1):
   face_element_pt->node_pt(0) = node_pt(0);
   face_element_pt->bulk_node_number(0) = 0;
   face_element_pt->normal_sign() = -1;
   
   //Set the pointer to the function that determines the bulk coordinates
   //in the face element
   face_element_pt->face_to_bulk_coordinate_fct_pt() = 
    &QElement1FaceToBulkCoordinates::face0;

   //Set the pointer to the function that determines the mapping of
   //derivatives
   face_element_pt->bulk_coordinate_derivatives_fct_pt() =
    &QElement1BulkCoordinateDerivatives::faces0;

   //Set the number of values stored when the node is part of the "bulk"
   //element.
   face_element_pt->nbulk_value(0) = required_nvalue(0);
   break;
   
   //Top, normal sign is positive (coordinate points out of element)
  case(1):
   face_element_pt->node_pt(0) = node_pt(1);
   face_element_pt->bulk_node_number(0) = 1;
   face_element_pt->normal_sign() = +1;

   //Set the pointer to the function that determines the bulk coordinates
   //in the face element
   face_element_pt->face_to_bulk_coordinate_fct_pt() = 
    &QElement1FaceToBulkCoordinates::face1;

   //Set the pointer to the function that determines the mapping of
   //derivatives
   face_element_pt->bulk_coordinate_derivatives_fct_pt() =
    &QElement1BulkCoordinateDerivatives::faces0;

   //Set the number of values stored when the node is part of the "bulk"
   //element.
   face_element_pt->nbulk_value(0) = required_nvalue(1);
   break;
   
   //Other cases
  default:
   std::ostringstream error_message;
   error_message << "Face_index should only take "
                 << "the values +/-1, not " << face_index << std::endl;

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

////////////////////////////////////////////////////////////////
//       2D Hermite elements
////////////////////////////////////////////////////////////////


//=======================================================================
/// Shape function for specific QHermiteElement<2>
//=======================================================================
template<>
void QHermiteElement<2>::shape(const Vector<double> &s, Shape &psi) 
 const
{
 //Local storage
 double Psi[2][2][2];

 //Call the OneDimensional Shape functions
 OneDimHermite::shape(s[0],Psi[0]);
 OneDimHermite::shape(s[1],Psi[1]);

 //Set up the two dimensional shape functions
 //Set up the functions at corner 0 
 //psi_0 = 1 when s1 = 0, s2 = 0
 psi(0,0) = Psi[0][0][0]*Psi[1][0][0];
 //dpsi_0/ds1 = 1 when s1 = 0, s2 = 0
 psi(0,1) = Psi[0][0][1]*Psi[1][0][0];
 //dpsi_0/ds2 = 1 when s1 = 0, s2 = 0
 psi(0,2) = Psi[0][0][0]*Psi[1][0][1];
 //dpsi_0/ds2ds1 = 1 when s1 = 0, s2 = 0
 psi(0,3) = Psi[0][0][1]*Psi[1][0][1];
    
 //Set up the functions at corner 1
 //psi_1 = 1 when s1 = 1, s2 = 0
 psi(1,0) = Psi[0][1][0]*Psi[1][0][0]; 
 //dpsi_1/ds1 = 1 when s1 = 1, s2 = 0
 psi(1,1) = Psi[0][1][1]*Psi[1][0][0];
 //dpsi_1/ds2 = 1 when s1 = 1, s2 = 0
 psi(1,2) = Psi[0][1][0]*Psi[1][0][1];
 //dpsi_1/ds1ds2  = 1 when s1 = 1, s2 = 0
 psi(1,3) = Psi[0][1][1]*Psi[1][0][1];

 //Set up the functions at the corner 2
 //psi_2 = 1 when s1 = 0, s2 = 1
 psi(2,0) = Psi[0][0][0]*Psi[1][1][0];
 //dpsi_2/ds1 = 1 when s1 = 0, s2 = 1
 psi(2,1) = Psi[0][0][1]*Psi[1][1][0];
 //dpsi_2/ds2 = 1 when s1 = 0, s2 = 1
 psi(2,2) = Psi[0][0][0]*Psi[1][1][1];
 //dpsi_2/ds2ds1 = 1 when s1 = 0, s2 = 1
 psi(2,3) = Psi[0][0][1]*Psi[1][1][1];
    
 //Set up the functions at corner 3
 //psi_3 = 1 when s1 = 1, s2 = 1
 psi(3,0) = Psi[0][1][0]*Psi[1][1][0]; 
 //dpsi_3/ds1 = 1 when s1 = 1, s2 = 1
 psi(3,1) = Psi[0][1][1]*Psi[1][1][0];
 //dpsi_3/ds2 = 1 when s1 = 1, s2 = 1
 psi(3,2) = Psi[0][1][0]*Psi[1][1][1];
 //dpsi_3/ds1ds2  = 1 when s1 = 1, s2 = 1
 psi(3,3) = Psi[0][1][1]*Psi[1][1][1];
}

//=======================================================================
/// Derivatives of shape functions for specific QHermiteElement<2>
//=======================================================================
template <>
void QHermiteElement<2>::dshape_local(const Vector<double> &s,
                                      Shape &psi, 
                                      DShape &dpsids) const
{
 //Local storage
 double Psi[2][2][2];
 double DPsi[2][2][2];
 //Call the OneDimensional Shape functions
 OneDimHermite::shape(s[0],Psi[0]);
 OneDimHermite::shape(s[1],Psi[1]);
 OneDimHermite::dshape(s[0],DPsi[0]);
 OneDimHermite::dshape(s[1],DPsi[1]);

 
 //Set up the two dimensional shape functions
 //Set up the functions at corner 0 
 //psi_0 = 1 when s1 = 0, s2 = 0
 psi(0,0) = Psi[0][0][0]*Psi[1][0][0];
 //dpsi_0/ds1 = 1 when s1 = 0, s2 = 0
 psi(0,1) = Psi[0][0][1]*Psi[1][0][0];
 //dpsi_0/ds2 = 1 when s1 = 0, s2 = 0
 psi(0,2) = Psi[0][0][0]*Psi[1][0][1];
 //dpsi_0/ds2ds1 = 1 when s1 = 0, s2 = 0
 psi(0,3) = Psi[0][0][1]*Psi[1][0][1];
    
 //Set up the functions at corner 1
 //psi_1 = 1 when s1 = 1, s2 = 0
 psi(1,0) = Psi[0][1][0]*Psi[1][0][0]; 
 //dpsi_1/ds1 = 1 when s1 = 1, s2 = 0
 psi(1,1) = Psi[0][1][1]*Psi[1][0][0];
 //dpsi_1/ds2 = 1 when s1 = 1, s2 = 0
 psi(1,2) = Psi[0][1][0]*Psi[1][0][1];
 //dpsi_1/ds1ds2  = 1 when s1 = 1, s2 = 0
 psi(1,3) = Psi[0][1][1]*Psi[1][0][1];

 //Set up the functions at the corner 2
 //psi_2 = 1 when s1 = 0, s2 = 1
 psi(2,0) = Psi[0][0][0]*Psi[1][1][0];
 //dpsi_2/ds1 = 1 when s1 = 0, s2 = 1
 psi(2,1) = Psi[0][0][1]*Psi[1][1][0];
 //dpsi_2/ds2 = 1 when s1 = 0, s2 = 1
 psi(2,2) = Psi[0][0][0]*Psi[1][1][1];
 //dpsi_2/ds2ds1 = 1 when s1 = 0, s2 = 1
 psi(2,3) = Psi[0][0][1]*Psi[1][1][1];
    
 //Set up the functions at corner 3
 //psi_3 = 1 when s1 = 1, s2 = 1
 psi(3,0) = Psi[0][1][0]*Psi[1][1][0]; 
 //dpsi_3/ds1 = 1 when s1 = 1, s2 = 1
 psi(3,1) = Psi[0][1][1]*Psi[1][1][0];
 //dpsi_3/ds2 = 1 when s1 = 1, s2 = 1
 psi(3,2) = Psi[0][1][0]*Psi[1][1][1];
 //dpsi_3/ds1ds2  = 1 when s1 = 1, s2 = 1
 psi(3,3) = Psi[0][1][1]*Psi[1][1][1];

 //FIRST DERIVATIVES 

 //D/Ds[0]

 //Set up the functions at corner 0 
 dpsids(0,0,0) = DPsi[0][0][0]*Psi[1][0][0];
 dpsids(0,1,0) = DPsi[0][0][1]*Psi[1][0][0];
 dpsids(0,2,0) = DPsi[0][0][0]*Psi[1][0][1];
 dpsids(0,3,0) = DPsi[0][0][1]*Psi[1][0][1];
    
 //Set up the functions at corner 1
 dpsids(1,0,0) = DPsi[0][1][0]*Psi[1][0][0]; 
 dpsids(1,1,0) = DPsi[0][1][1]*Psi[1][0][0];
 dpsids(1,2,0) = DPsi[0][1][0]*Psi[1][0][1];
 dpsids(1,3,0) = DPsi[0][1][1]*Psi[1][0][1];

 //Set up the functions at the corner 2
 dpsids(2,0,0) = DPsi[0][0][0]*Psi[1][1][0];
 dpsids(2,1,0) = DPsi[0][0][1]*Psi[1][1][0];
 dpsids(2,2,0) = DPsi[0][0][0]*Psi[1][1][1];
 dpsids(2,3,0) = DPsi[0][0][1]*Psi[1][1][1];
    
 //Set up the functions at corner 3
 dpsids(3,0,0) = DPsi[0][1][0]*Psi[1][1][0]; 
 dpsids(3,1,0) = DPsi[0][1][1]*Psi[1][1][0];
 dpsids(3,2,0) = DPsi[0][1][0]*Psi[1][1][1];
 dpsids(3,3,0) = DPsi[0][1][1]*Psi[1][1][1];

 //D/Ds[1]

 //Set up the functions at corner 0 
 dpsids(0,0,1) = Psi[0][0][0]*DPsi[1][0][0];
 dpsids(0,1,1) = Psi[0][0][1]*DPsi[1][0][0];
 dpsids(0,2,1) = Psi[0][0][0]*DPsi[1][0][1];
 dpsids(0,3,1) = Psi[0][0][1]*DPsi[1][0][1];
    
 //Set up the functions at corner 1
 dpsids(1,0,1) = Psi[0][1][0]*DPsi[1][0][0]; 
 dpsids(1,1,1) = Psi[0][1][1]*DPsi[1][0][0];
 dpsids(1,2,1) = Psi[0][1][0]*DPsi[1][0][1];
 dpsids(1,3,1) = Psi[0][1][1]*DPsi[1][0][1];

 //Set up the functions at the corner 2
 dpsids(2,0,1) = Psi[0][0][0]*DPsi[1][1][0];
 dpsids(2,1,1) = Psi[0][0][1]*DPsi[1][1][0];
 dpsids(2,2,1) = Psi[0][0][0]*DPsi[1][1][1];
 dpsids(2,3,1) = Psi[0][0][1]*DPsi[1][1][1];
    
 //Set up the functions at corner 3
 dpsids(3,0,1) = Psi[0][1][0]*DPsi[1][1][0]; 
 dpsids(3,1,1) = Psi[0][1][1]*DPsi[1][1][0];
 dpsids(3,2,1) = Psi[0][1][0]*DPsi[1][1][1];
 dpsids(3,3,1) = Psi[0][1][1]*DPsi[1][1][1];

}

//======================================================================
/// Second derivatives of the shape functions wrt local coordinates. 
/// 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$ 
//======================================================================
template<>
void QHermiteElement<2>::d2shape_local(const Vector<double> &s,
                                       Shape &psi,
                                       DShape &dpsids,
                                       DShape &d2psids) const
{
 //Local storage
 double Psi[2][2][2];
 double DPsi[2][2][2];
 double D2Psi[2][2][2];

 //Call the OneDimensional Shape functions
 OneDimHermite::shape(s[0],Psi[0]);
 OneDimHermite::shape(s[1],Psi[1]);
 OneDimHermite::dshape(s[0],DPsi[0]);
 OneDimHermite::dshape(s[1],DPsi[1]);
 OneDimHermite::d2shape(s[0],D2Psi[0]);
 OneDimHermite::d2shape(s[1],D2Psi[1]);
 
 //Set up the two dimensional shape functions
 //Set up the functions at corner 0 
 //psi_0 = 1 when s1 = 0, s2 = 0
 psi(0,0) = Psi[0][0][0]*Psi[1][0][0];
 //dpsi_0/ds1 = 1 when s1 = 0, s2 = 0
 psi(0,1) = Psi[0][0][1]*Psi[1][0][0];
 //dpsi_0/ds2 = 1 when s1 = 0, s2 = 0
 psi(0,2) = Psi[0][0][0]*Psi[1][0][1];
 //dpsi_0/ds2ds1 = 1 when s1 = 0, s2 = 0
 psi(0,3) = Psi[0][0][1]*Psi[1][0][1];
    
 //Set up the functions at corner 1
 //psi_1 = 1 when s1 = 1, s2 = 0
 psi(1,0) = Psi[0][1][0]*Psi[1][0][0]; 
 //dpsi_1/ds1 = 1 when s1 = 1, s2 = 0
 psi(1,1) = Psi[0][1][1]*Psi[1][0][0];
 //dpsi_1/ds2 = 1 when s1 = 1, s2 = 0
 psi(1,2) = Psi[0][1][0]*Psi[1][0][1];
 //dpsi_1/ds1ds2  = 1 when s1 = 1, s2 = 0
 psi(1,3) = Psi[0][1][1]*Psi[1][0][1];

 //Set up the functions at the corner 2
 //psi_2 = 1 when s1 = 0, s2 = 1
 psi(2,0) = Psi[0][0][0]*Psi[1][1][0];
 //dpsi_2/ds1 = 1 when s1 = 0, s2 = 1
 psi(2,1) = Psi[0][0][1]*Psi[1][1][0];
 //dpsi_2/ds2 = 1 when s1 = 0, s2 = 1
 psi(2,2) = Psi[0][0][0]*Psi[1][1][1];
 //dpsi_2/ds2ds1 = 1 when s1 = 0, s2 = 1
 psi(2,3) = Psi[0][0][1]*Psi[1][1][1];
    
 //Set up the functions at corner 3
 //psi_3 = 1 when s1 = 1, s2 = 1
 psi(3,0) = Psi[0][1][0]*Psi[1][1][0]; 
 //dpsi_3/ds1 = 1 when s1 = 1, s2 = 1
 psi(3,1) = Psi[0][1][1]*Psi[1][1][0];
 //dpsi_3/ds2 = 1 when s1 = 1, s2 = 1
 psi(3,2) = Psi[0][1][0]*Psi[1][1][1];
 //dpsi_3/ds1ds2  = 1 when s1 = 1, s2 = 1
 psi(3,3) = Psi[0][1][1]*Psi[1][1][1];

 //FIRST DERIVATIVES
 
 //D/Ds[0]

 //Set up the functions at corner 0 
 dpsids(0,0,0) = DPsi[0][0][0]*Psi[1][0][0];
 dpsids(0,1,0) = DPsi[0][0][1]*Psi[1][0][0];
 dpsids(0,2,0) = DPsi[0][0][0]*Psi[1][0][1];
 dpsids(0,3,0) = DPsi[0][0][1]*Psi[1][0][1];
    
 //Set up the functions at corner 1
 dpsids(1,0,0) = DPsi[0][1][0]*Psi[1][0][0]; 
 dpsids(1,1,0) = DPsi[0][1][1]*Psi[1][0][0];
 dpsids(1,2,0) = DPsi[0][1][0]*Psi[1][0][1];
 dpsids(1,3,0) = DPsi[0][1][1]*Psi[1][0][1];

 //Set up the functions at the corner 2
 dpsids(2,0,0) = DPsi[0][0][0]*Psi[1][1][0];
 dpsids(2,1,0) = DPsi[0][0][1]*Psi[1][1][0];
 dpsids(2,2,0) = DPsi[0][0][0]*Psi[1][1][1];
 dpsids(2,3,0) = DPsi[0][0][1]*Psi[1][1][1];
    
 //Set up the functions at corner 3
 dpsids(3,0,0) = DPsi[0][1][0]*Psi[1][1][0]; 
 dpsids(3,1,0) = DPsi[0][1][1]*Psi[1][1][0];
 dpsids(3,2,0) = DPsi[0][1][0]*Psi[1][1][1];
 dpsids(3,3,0) = DPsi[0][1][1]*Psi[1][1][1];

 // D/Ds[1]

 //Set up the functions at corner 0 
 dpsids(0,0,1) = Psi[0][0][0]*DPsi[1][0][0];
 dpsids(0,1,1) = Psi[0][0][1]*DPsi[1][0][0];
 dpsids(0,2,1) = Psi[0][0][0]*DPsi[1][0][1];
 dpsids(0,3,1) = Psi[0][0][1]*DPsi[1][0][1];
    
 //Set up the functions at corner 1
 dpsids(1,0,1) = Psi[0][1][0]*DPsi[1][0][0]; 
 dpsids(1,1,1) = Psi[0][1][1]*DPsi[1][0][0];
 dpsids(1,2,1) = Psi[0][1][0]*DPsi[1][0][1];
 dpsids(1,3,1) = Psi[0][1][1]*DPsi[1][0][1];

 //Set up the functions at the corner 2
 dpsids(2,0,1) = Psi[0][0][0]*DPsi[1][1][0];
 dpsids(2,1,1) = Psi[0][0][1]*DPsi[1][1][0];
 dpsids(2,2,1) = Psi[0][0][0]*DPsi[1][1][1];
 dpsids(2,3,1) = Psi[0][0][1]*DPsi[1][1][1];
    
 //Set up the functions at corner 3
 dpsids(3,0,1) = Psi[0][1][0]*DPsi[1][1][0]; 
 dpsids(3,1,1) = Psi[0][1][1]*DPsi[1][1][0];
 dpsids(3,2,1) = Psi[0][1][0]*DPsi[1][1][1];
 dpsids(3,3,1) = Psi[0][1][1]*DPsi[1][1][1];

 //SECOND DERIVATIVES
 //Convention: index 0 is d^2/ds[0]^2, 
 //            index 1 is d^2/ds[1]^2,
 //            index 2 is the mixed derivative

 //D^2/Ds[0]^2

 //Set up the functions at corner 0 
 d2psids(0,0,0) = D2Psi[0][0][0]*Psi[1][0][0];
 d2psids(0,1,0) = D2Psi[0][0][1]*Psi[1][0][0];
 d2psids(0,2,0) = D2Psi[0][0][0]*Psi[1][0][1];
 d2psids(0,3,0) = D2Psi[0][0][1]*Psi[1][0][1];
    
 //Set up the functions at corner 1
 d2psids(1,0,0) = D2Psi[0][1][0]*Psi[1][0][0]; 
 d2psids(1,1,0) = D2Psi[0][1][1]*Psi[1][0][0];
 d2psids(1,2,0) = D2Psi[0][1][0]*Psi[1][0][1];
 d2psids(1,3,0) = D2Psi[0][1][1]*Psi[1][0][1];

 //Set up the functions at the corner 2
 d2psids(2,0,0) = D2Psi[0][0][0]*Psi[1][1][0];
 d2psids(2,1,0) = D2Psi[0][0][1]*Psi[1][1][0];
 d2psids(2,2,0) = D2Psi[0][0][0]*Psi[1][1][1];
 d2psids(2,3,0) = D2Psi[0][0][1]*Psi[1][1][1];
    
 //Set up the functions at corner 3
 d2psids(3,0,0) = D2Psi[0][1][0]*Psi[1][1][0]; 
 d2psids(3,1,0) = D2Psi[0][1][1]*Psi[1][1][0];
 d2psids(3,2,0) = D2Psi[0][1][0]*Psi[1][1][1];
 d2psids(3,3,0) = D2Psi[0][1][1]*Psi[1][1][1];

 // D^2/Ds[1]^2

 //Set up the functions at corner 0 
 d2psids(0,0,1) = Psi[0][0][0]*D2Psi[1][0][0];
 d2psids(0,1,1) = Psi[0][0][1]*D2Psi[1][0][0];
 d2psids(0,2,1) = Psi[0][0][0]*D2Psi[1][0][1];
 d2psids(0,3,1) = Psi[0][0][1]*D2Psi[1][0][1];
    
 //Set up the functions at corner 1
 d2psids(1,0,1) = Psi[0][1][0]*D2Psi[1][0][0]; 
 d2psids(1,1,1) = Psi[0][1][1]*D2Psi[1][0][0];
 d2psids(1,2,1) = Psi[0][1][0]*D2Psi[1][0][1];
 d2psids(1,3,1) = Psi[0][1][1]*D2Psi[1][0][1];

 //Set up the functions at the corner 2
 d2psids(2,0,1) = Psi[0][0][0]*D2Psi[1][1][0];
 d2psids(2,1,1) = Psi[0][0][1]*D2Psi[1][1][0];
 d2psids(2,2,1) = Psi[0][0][0]*D2Psi[1][1][1];
 d2psids(2,3,1) = Psi[0][0][1]*D2Psi[1][1][1];
    
 //Set up the functions at corner 3
 d2psids(3,0,1) = Psi[0][1][0]*D2Psi[1][1][0]; 
 d2psids(3,1,1) = Psi[0][1][1]*D2Psi[1][1][0];
 d2psids(3,2,1) = Psi[0][1][0]*D2Psi[1][1][1];
 d2psids(3,3,1) = Psi[0][1][1]*D2Psi[1][1][1];

 //D^2/Ds[0]Ds[1]

 //Set up the functions at corner 0 
 d2psids(0,0,2) = DPsi[0][0][0]*DPsi[1][0][0];
 d2psids(0,1,2) = DPsi[0][0][1]*DPsi[1][0][0];
 d2psids(0,2,2) = DPsi[0][0][0]*DPsi[1][0][1];
 d2psids(0,3,2) = DPsi[0][0][1]*DPsi[1][0][1];
    
 //Set up the functions at corner 1
 d2psids(1,0,2) = DPsi[0][1][0]*DPsi[1][0][0]; 
 d2psids(1,1,2) = DPsi[0][1][1]*DPsi[1][0][0];
 d2psids(1,2,2) = DPsi[0][1][0]*DPsi[1][0][1];
 d2psids(1,3,2) = DPsi[0][1][1]*DPsi[1][0][1];

 //Set up the functions at the corner 2
 d2psids(2,0,2) = DPsi[0][0][0]*DPsi[1][1][0];
 d2psids(2,1,2) = DPsi[0][0][1]*DPsi[1][1][0];
 d2psids(2,2,2) = DPsi[0][0][0]*DPsi[1][1][1];
 d2psids(2,3,2) = DPsi[0][0][1]*DPsi[1][1][1];
    
 //Set up the functions at corner 3
 d2psids(3,0,2) = DPsi[0][1][0]*DPsi[1][1][0]; 
 d2psids(3,1,2) = DPsi[0][1][1]*DPsi[1][1][0];
 d2psids(3,2,2) = DPsi[0][1][0]*DPsi[1][1][1];
 d2psids(3,3,2) = DPsi[0][1][1]*DPsi[1][1][1];
}



//===========================================================
/// The output function for QHermiteElement<2,ORDER>
//===========================================================
template<>
void QHermiteElement<2>::output(std::ostream &outfile)
{
 //Tecplot header info 
 outfile << "ZONE I=" << 2 << ", J=" << 2 << std::endl;

 //Find the dimension of the node
 unsigned n_dim = this->nodal_dimension();

 //Loop over element nodes
 for(unsigned l2=0;l2<2;l2++)
  {
   for(unsigned l1=0;l1<2;l1++)
    {
     unsigned l = l2*2 + l1;

     //Loop over the dimensions and output the position
     for(unsigned i=0;i<n_dim;i++)
      {
       outfile << node_pt(l)->x(i) << " ";
      }

     //Find out number of types of dof stored at each node
     unsigned n_position_type = node_pt(l)->nposition_type();
     //Loop over the additional positional dofs
     for(unsigned k=1;k<n_position_type;k++)
      {
       for(unsigned i=0;i<n_dim;i++)
        {
         outfile << node_pt(l)->x_gen(k,i) << " ";
        }
      }

     //Find out how many data values at the node
     unsigned initial_nvalue = node_pt(l)->nvalue();
     //Loop over the data and output whether pinned or not
     for(unsigned i=0;i<initial_nvalue;i++)
      {
       outfile << node_pt(l)->is_pinned(i) << " ";
      }
     outfile << std::endl;
    }
  }
  outfile << std::endl;
}

//=======================================================================
///The output function for n_plot points in each coordinate direction
//=======================================================================
template<>
void QHermiteElement<2>::output(std::ostream &outfile, const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(2);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;

 //Find the dimension of the first node
 unsigned n_dim = this->nodal_dimension();

 //Loop over plot points
 for(unsigned l2=0;l2<n_plot;l2++)
  {
   s[1] = -1.0 + l2*2.0/(n_plot-1);
   
   for(unsigned l1=0;l1<n_plot;l1++)
    {
     s[0] = -1.0 + l1*2.0/(n_plot-1);

     //Output the coordinates
     for (unsigned i=0;i<n_dim;i++)
      {
       outfile << interpolated_x(s,i) << " " ;
      }
    }
  }
 outfile << std::endl;
}









//===========================================================
/// The C-style output function for QHermiteElement<2,ORDER>
//===========================================================
template<>
void QHermiteElement<2>::output(FILE* file_pt)
{
 //Tecplot header info 
 fprintf(file_pt,"ZONE I=2, J=2");
 
 //Find the dimension of the node
 unsigned n_dim = this->nodal_dimension();

 //Loop over element nodes
 for(unsigned l2=0;l2<2;l2++)
  {
   for(unsigned l1=0;l1<2;l1++)
    {
     unsigned l = l2*2 + l1;

     //Loop over the dimensions and output the position
     for(unsigned i=0;i<n_dim;i++)
      {
       fprintf(file_pt,"%g ",node_pt(l)->x(i));
      }

     //Find out number of types of dof stored at each node
     unsigned n_position_type = node_pt(l)->nposition_type();
     //Loop over the additional positional dofs
     for(unsigned k=1;k<n_position_type;k++)
      {
       for(unsigned i=0;i<n_dim;i++)
        {
         fprintf(file_pt,"%g ",node_pt(l)->x_gen(k,i));
        }
      }

     //Find out how many data values at the node
     unsigned initial_nvalue = node_pt(l)->nvalue();
     //Loop over the data and output whether pinned or not
     for(unsigned i=0;i<initial_nvalue;i++)
      {
       fprintf(file_pt,"%i ",node_pt(l)->is_pinned(i));
      }
     fprintf(file_pt,"\n");
    }
  }
  fprintf(file_pt,"\n");
}

//=======================================================================
///The C-style output function for n_plot points in each coordinate direction
//=======================================================================
template<>
void QHermiteElement<2>::output(FILE* file_pt, const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(2);

 //Tecplot header info 
 fprintf(file_pt,"ZONE I=%i, J=%i \n",n_plot,n_plot);

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();

 //Loop over plot points
 for(unsigned l2=0;l2<n_plot;l2++)
  {
   s[1] = -1.0 + l2*2.0/(n_plot-1);
   for(unsigned l1=0;l1<n_plot;l1++)
    {
     s[0] = -1.0 + l1*2.0/(n_plot-1);

     //Output the coordinates
     for (unsigned i=0;i<n_dim;i++)
      {
       fprintf(file_pt,"%g ",interpolated_x(s,i));
      }
    }
  }
 fprintf(file_pt,"\n");
}










//=======================================================================
/// Function to setup geometrical information for lower-dimensional 
/// FaceElements (of type QHermiteElement<1>).
//=======================================================================
template<>
void QHermiteElement<2>::build_face_element(const int &face_index,
                                            FaceElement *face_element_pt)
{
 // Set the nodal dimension from the "bulk"
 face_element_pt->set_nodal_dimension(node_pt(0)->ndim());

 // Set the pointer to the "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

 //Resize the bulk_position_type translation scheme to the number of position
 //types in the 1D element: 2 (position and slope)
 face_element_pt->bulk_position_type_resize(2);
 
 // Resize the storage for the original number of values at 
 // the two nodes of the FaceElement
 face_element_pt->nbulk_value_resize(2);
 
 // Resize the storage for the bulk node numbers coressponding to the
 // two nodes of the FaceElement
 face_element_pt->bulk_node_number_resize(2);
 
 // Set the face index in the face element
 // The faces are
 //                       +1    East       
 //                       -1    West
 //                       +2    North
 //                       -2    South
 
 //Set the face index in the face element
 face_element_pt->face_index() = face_index;
 
 //Now set up the node pointers
 //The convention here is that interior_tangent X tangent X tangent
 //is the OUTWARD normal
 //IMPORTANT NOTE: Need to ensure that numbering is consistent here
 //i.e. node numbers increase in positive x and y directions as they should
 //If not, normals will be inward. 
 switch(face_index)
  {
   unsigned bulk_number;
   //West face, normal sign is positive
  case(-1):
   //Set the pointer to the bulk coordinate translation scheme
   face_element_pt->face_to_bulk_coordinate_fct_pt() = 
    &QElement2FaceToBulkCoordinates::face0;
   
   //Set the pointer to the derivative mappings
   face_element_pt->bulk_coordinate_derivatives_fct_pt() = 
    &QElement2BulkCoordinateDerivatives::faces0;
   
   for(unsigned i=0;i<2;i++) 
    {
     bulk_number = i*2;
     face_element_pt->node_pt(i) = node_pt(bulk_number);
     face_element_pt->bulk_node_number(i) = bulk_number;
     face_element_pt->normal_sign() = 1;
     //Set the number of values originally stored at this node
     face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
     //Set the position type for the slope, which is in the s[1] direction, 
     //so is position_type 2 in the "bulk" element
     face_element_pt->bulk_position_type(1) = 2;
    }
   break;
   //South face, normal sign is positive
  case(-2):
   //Set the pointer to the bulk coordinate translation scheme
   face_element_pt->face_to_bulk_coordinate_fct_pt() = 
    &QElement2FaceToBulkCoordinates::face1;
   
   //Set the pointer to the derivative mappings
   face_element_pt->bulk_coordinate_derivatives_fct_pt() = 
    &QElement2BulkCoordinateDerivatives::faces1;
   
   for(unsigned i=0;i<2;i++) 
    {
     bulk_number = i;
     face_element_pt->node_pt(i) = node_pt(bulk_number);
     face_element_pt->bulk_node_number(i) = bulk_number;
     face_element_pt->normal_sign() = 1;
     //Set the number of values originally stored at this node
     face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
     //Set the position type for the slope, which is in the s[0] direction,
      //so is position_type 1 in "bulk" element
     face_element_pt->bulk_position_type(1) = 1;
    }
   break;
   //East face, normal sign is negative
  case(1):
   //Set the pointer to the bulk coordinate translation scheme
   face_element_pt->face_to_bulk_coordinate_fct_pt() = 
    &QElement2FaceToBulkCoordinates::face2;
   
   //Set the pointer to the derivative mappings
   face_element_pt->bulk_coordinate_derivatives_fct_pt() = 
    &QElement2BulkCoordinateDerivatives::faces0;
   
   for(unsigned i=0;i<2;i++) 
     {
      bulk_number = 2*i + 1;
      face_element_pt->node_pt(i) = node_pt(bulk_number);
      face_element_pt->bulk_node_number(i) = bulk_number;
      face_element_pt->normal_sign() = -1;
      //Set the number of values originally stored at this node
      face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
      //Set the position type for the slope, which is in the s[1] direction,
      //so is position_type 2 in the bulk element
      face_element_pt->bulk_position_type(1) = 2;  
     }
   break;
   //North face, normal sign is negative
  case(2):
   //Set the pointer to the bulk coordinate translation scheme
   face_element_pt->face_to_bulk_coordinate_fct_pt() = 
    &QElement2FaceToBulkCoordinates::face3;
   
   //Set the pointer to the derivative mappings
   face_element_pt->bulk_coordinate_derivatives_fct_pt() = 
    &QElement2BulkCoordinateDerivatives::faces1;
   
   for(unsigned i=0;i<2;i++)
    {
     bulk_number = 2 + i;
     face_element_pt->node_pt(i) = node_pt(bulk_number);
     face_element_pt->bulk_node_number(i) = bulk_number;
     face_element_pt->normal_sign() = -1;
     //Set the number of values originally stored at this node
     face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
     //Set the position type for the slope, which is in the s[0] direction,
     //so is position_type 1 in the "bulk" element
     face_element_pt->bulk_position_type(1) = 1;  
    }
   break;
   
   //Now cover the other cases
  default:
   std::ostringstream error_message;
   error_message << "Face index should only take the values +/- 1 or +/- 2,"
                 << " not " << face_index << std::endl;
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);

   }
 }

/////////////////////////////////////////////////////////////////////
//     1D SolidQHermiteElements
//////////////////////////////////////////////////////////////////////


//=====================================================================
/// Overload the output function
//====================================================================
template<unsigned DIM>
void SolidQHermiteElement<DIM>::output(std::ostream &outfile) 
{QHermiteElement<DIM>::output(outfile);}


//=======================================================================
/// The output function for n_plot points in each coordinate direction
/// for the 1D element
//=======================================================================
template <>
void SolidQHermiteElement<1>::output(std::ostream &outfile, 
                                    const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(1);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << std::endl;

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();

 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over plot points 
 for(unsigned l=0;l<n_plot;l++)
  {
   s[0] = -1.0 + l*2.0/(n_plot-1);

   //Output the Eulerian coordinates
   for (unsigned i=0;i<n_dim;i++)
    {
     outfile << interpolated_x(s,i) << " " ;
    }

   //Output the Lagrangian coordinates
   for(unsigned i=0;i<n_lagr;i++)
    {
     outfile << interpolated_xi(s,i) << " " ;
    }
   outfile << std::endl;
  }
}



//=====================================================================
/// Overload the C-style output function
//====================================================================
template<unsigned DIM>
void SolidQHermiteElement<DIM>::output(FILE* file_pt)
{QHermiteElement<DIM>::output(file_pt);}


//=======================================================================
/// The C-style output function for n_plot points in each coordinate direction
/// for the 1D element
//=======================================================================
template <>
void SolidQHermiteElement<1>::output(FILE* file_pt,
                                    const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(1);

 //Tecplot header info 
 fprintf(file_pt,"ZONE I=%i\n",n_plot);

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();
 
 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
 
 //Loop over plot points 
 for(unsigned l=0;l<n_plot;l++)
  {
   s[0] = -1.0 + l*2.0/(n_plot-1);
   
   //Output the Eulerian coordinates
   for (unsigned i=0;i<n_dim;i++)
    {
     fprintf(file_pt,"%g ",interpolated_x(s,i));
    }
   
   //Output the Lagrangian coordinates
   for(unsigned i=0;i<n_lagr;i++)
    {
     fprintf(file_pt,"%g ",interpolated_xi(s,i));
    }
   fprintf(file_pt,"\n");
  }
}


///////////////////////////////////////////////////////////////////////////
//   2D SolidQHermiteElements
//////////////////////////////////////////////////////////////////////////


//=====================================================================
/// The output function for any number of points per element
//=====================================================================
template<>
void SolidQHermiteElement<2>::output(std::ostream &outfile, const unsigned &n_p)
{
 //Local variables
 Vector<double> s(2);

 //Tecplot header info 
 outfile << "ZONE I=" << n_p << ", J=" << n_p << std::endl;

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();

 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over element nodes
 for(unsigned l2=0;l2<n_p;l2++)
  {
   s[1] = -1.0 + l2*2.0/(n_p-1);
   for(unsigned l1=0;l1<n_p;l1++)
    {
     s[0] = -1.0 + l1*2.0/(n_p-1);
     
     //Output the Eulerian coordinates
     for (unsigned i=0;i<n_dim;i++)
      {
       outfile << interpolated_x(s,i) << " " ;
      }
     
     //Output the Lagrangian coordinates
     for(unsigned i=0;i<n_lagr;i++)
      {
       outfile << interpolated_xi(s,i) << " " ;
      }
     outfile << std::endl;
    }
  }
 outfile << std::endl;
}


//=====================================================================
/// The C-style output function for any number of points per element
//=====================================================================
template<>
void SolidQHermiteElement<2>::output(FILE* file_pt, const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(2);

 //Tecplot header info 
 fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();

 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over element nodes
 for(unsigned l2=0;l2<n_plot;l2++)
  {
   s[1] = -1.0 + l2*2.0/(n_plot-1);
   for(unsigned l1=0;l1<n_plot;l1++)
    {
     s[0] = -1.0 + l1*2.0/(n_plot-1);
     
     //Output the Eulerian coordinates
     for (unsigned i=0;i<n_dim;i++)
      {
       fprintf(file_pt,"%g ",interpolated_x(s,i));
      }
     
     //Output the Lagrangian coordinates
     for(unsigned i=0;i<n_lagr;i++)
      {
       fprintf(file_pt,"%g ",interpolated_xi(s,i));
      }
     fprintf(file_pt,"\n");
    }
  }
 fprintf(file_pt,"\n");
}



//=======================================================================
/// Function to setup geometrical information for lower-dimensional 
/// FaceElements for the solid hermite elements. We need to 
/// construct the basic  element and then sort out the Lagrangian 
/// coordinates
//=======================================================================
template<unsigned DIM>
void SolidQHermiteElement<DIM>::
build_face_element(const int &face_index,
                   FaceElement *face_element_pt)
{
 //Build the standard non-solid FaceElement
 QHermiteElement<DIM>::
  build_face_element(face_index,face_element_pt);
 
 //Set the Lagrangian dimension from the first node of the present element
 dynamic_cast<SolidFiniteElement*>(face_element_pt)->
  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
}

//==================================================================
/// Force build of templates
//==================================================================
template class QHermiteElement<1>;
template class QHermiteElement<2>;
template class DiagQHermiteElement<1>;
template class DiagQHermiteElement<2>;


template class SolidQHermiteElement<1>;
template class SolidQHermiteElement<2>;
template class SolidDiagQHermiteElement<1>;
template class SolidDiagQHermiteElement<2>;

}