hermite_element_quad_mesh.template.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//====================================================================
// the mesh assembly functions for the hermite element quad mesh
#ifndef OOMPH_HERMITE_QUAD_MESH_TEMPLATE_CC
#define OOMPH_HERMITE_QUAD_MESH_TEMPLATE_CC


// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

//oomph-lib includes
#include "hermite_element_quad_mesh.template.h"


namespace oomph
{


//=============================================================================
/// \short sets the generalised position of the node (i.e. - x_i, dx_i/ds_0, 
/// dx_i/ds_1 & d2x_i/ds_0ds_1 for i = 1,2). Takes the x,y coordinates of the 
/// node from which its position can be determined.
//=============================================================================
template<class ELEMENT>
void HermiteQuadMesh<ELEMENT>::set_position_of_node(const unsigned& node_num_x,
                                                    const unsigned& node_num_y,
                                                    Node* node_pt)
{
 
 // get the generalised macro element position of the node
 DenseMatrix<double> m_gen(4,2);
 generalised_macro_element_position_of_node(node_num_x,node_num_y,m_gen);
 
 // copy the macro element coordinates
 Vector<double> node_macro_position(2);
 node_macro_position[0] = m_gen(0,0);
 node_macro_position[1] = m_gen(0,1);
 
 // get the global position of the nodes
 Vector<double> x_node(2);
 Domain_pt->macro_element_pt(0)->macro_map(node_macro_position,x_node);
 
 // get the jacobian of the mapping from macro to eulerian coordinates
 DenseMatrix<double> jacobian_MX(2,2);
 Domain_pt->macro_element_pt(0)->
  assemble_macro_to_eulerian_jacobian(node_macro_position,jacobian_MX);
 
 
 // get the jacobian2 of the mapping from macro to eulerian coordinates
 DenseMatrix<double> jacobian2_MX(3,2);
 Domain_pt->macro_element_pt(0)->
  assemble_macro_to_eulerian_jacobian2(node_macro_position,jacobian2_MX);   
 
 // set x_0
 node_pt->x_gen(0,0) = x_node[0];
 // set x_1
 node_pt->x_gen(0,1) = x_node[1];
 
 // set dxi/dsj for i = 0,1 and j = 0,1
 // computation : dxi/dsj = dm0/dsj*dxi/dm0 + dm1/dsj*dxi/dm1
 for (unsigned i = 0; i < 2; i++)
  {
   for (unsigned j = 0; j < 2; j++)
    {
     node_pt->x_gen(j+1,i) = m_gen(j+1,0)*jacobian_MX(0,i) 
      + m_gen(j+1,1)*jacobian_MX(1,i);
    }
  }
 
 // set d2xi/ds0ds1 for i = 0,1
 // d2xi/ds0ds1 = d2m0/ds0ds1*dxi/dm0 
 //               + d2m1/ds0ds1*dxi/dm1
 //               + (dm0/ds1*d2xi/dm0^2 + dm1/ds1*d2xi/dm0dm1)*dm0/ds0
 //               + (dm0/ds1*d2xi/dm0dm1 + dm1/ds1*d2xi/dm1^2)*dm1/ds0   
 for (unsigned i = 0; i < 2; i++)
  {
   node_pt->x_gen(3,i) = m_gen(3,0)*jacobian_MX(0,i) 
    + m_gen(3,1)*jacobian_MX(1,i) 
    + m_gen(1,0)*(m_gen(2,0)*jacobian2_MX(0,i) + 
                  m_gen(2,1)*jacobian2_MX(2,i))
    + m_gen(1,1)*(m_gen(2,0)*jacobian2_MX(2,i) + 
                  m_gen(2,1)*jacobian2_MX(1,i));
  }
}


 //=============================================================================
/// \short sets the generalised position of the node (i.e. - x_i, dx_i/ds_0, 
/// dx_i/ds_1 & d2x_i/ds_0ds_1 for i = 1,2). Takes the x,y coordinates of the 
/// node from which its position can be determined. Also sets coordinates
/// on boundary vector for the node to be the generalised position of the node
/// in macro element coordinates
//=============================================================================
template<class ELEMENT>
void HermiteQuadMesh<ELEMENT>::
set_position_of_boundary_node(const unsigned& node_num_x,
                              const unsigned& node_num_y, 
                              BoundaryNode<Node>* node_pt)
{
 
 // get the generalised macro element position of the node
 DenseMatrix<double> m_gen(4,2,0.0);
 generalised_macro_element_position_of_node(node_num_x,node_num_y,m_gen);
 
 // copy the macro element coordinates
 Vector<double> node_macro_position(2);
 node_macro_position[0] = m_gen(0,0);
 node_macro_position[1] = m_gen(0,1);
 
 // get the global position of the nodes
 Vector<double> x_node(2);
 Domain_pt->macro_element_pt(0)->macro_map(node_macro_position,x_node);
 
 // get the jacobian of the mapping from macro to eulerian coordinates
 DenseMatrix<double> jacobian_MX(2,2);
 Domain_pt->macro_element_pt(0)->
  assemble_macro_to_eulerian_jacobian(node_macro_position,jacobian_MX);
 
 // get the jacobian2 of the mapping from macro to eulerian coordinates
 DenseMatrix<double> jacobian2_MX(3,2);
 Domain_pt->macro_element_pt(0)->
  assemble_macro_to_eulerian_jacobian2(node_macro_position,jacobian2_MX); 
 
 // set x_0
 node_pt->x_gen(0,0) = x_node[0];
 // set x_1
 node_pt->x_gen(0,1) = x_node[1];
 
 // set dxi/dsj for i = 0,1 and j = 0,1
 // computation : dxi/dsj = dm0/dsj*dxi/dm0 + dm1/dsj*dxi/dm1
 for (unsigned i = 0; i < 2; i++)
  {
   for (unsigned j = 0; j < 2; j++)
    {
     node_pt->x_gen(j+1,i) = m_gen(j+1,0)*jacobian_MX(0,i) 
      + m_gen(j+1,1)*jacobian_MX(1,i);
    }
  }
 
 // set d2xi/ds0ds1 for i = 0,1
 // d2xi/ds0ds1 = d2m0/ds0ds1*dxi/dm0 
 //               + d2m1/ds0ds1*dxi/dm1
 //               + (dm0/ds1*d2xi/dm0^2 + dm1/ds1*d2xi/dm0dm1)*dm0/ds0
 //               + (dm0/ds1*d2xi/dm0dm1 + dm1/ds1*d2xi/dm1^2)*dm1/ds0   
 for (unsigned i = 0; i < 2; i++)
  {
   node_pt->x_gen(3,i) = m_gen(3,0)*jacobian_MX(0,i) 
    + m_gen(3,1)*jacobian_MX(1,i) 
    + m_gen(1,0)*(m_gen(2,0)*jacobian2_MX(0,i) + 
                  m_gen(2,1)*jacobian2_MX(2,i))
    + m_gen(1,1)*(m_gen(2,0)*jacobian2_MX(2,i) + 
                    m_gen(2,1)*jacobian2_MX(1,i));
  }
 
 
 // pass generalised macro element position to vector to pass to boundary
 // node
 for (unsigned b = 0; b < 4; b++)
  {
   if (node_pt->is_on_boundary(b))
    { 
     Vector<double> boundary_macro_position(2,0);
     boundary_macro_position[0] = m_gen(0,b%2);
     boundary_macro_position[1] = m_gen(1+b%2,b%2);
     node_pt->set_coordinates_on_boundary(b,boundary_macro_position);
    }
  }
}


//=============================================================================
/// \short computes the generalised position of the node at position 
/// (node_num_x, node_num_y) in the macro element coordinate scheme.
///     index 0 of m_gen : 0 - m_i
///                        1 - dm_i/ds_0
///                        2 - dm_i/ds_1
///                        3 - d2m_i/ds_0ds_1    (where i is index 1 of m_gen)
//=============================================================================
template<class ELEMENT>
void HermiteQuadMesh<ELEMENT>::generalised_macro_element_position_of_node
(const unsigned& node_num_x,const unsigned& node_num_y, 
 DenseMatrix<double>& m_gen)
{  
   // obtain position of node
   Vector<double> s_macro_N(2);
   macro_coordinate_position(node_num_x,node_num_y,s_macro_N);

   // set position of node in macro element coordinates
   m_gen(0,0) = s_macro_N[0];
   m_gen(0,1) = s_macro_N[1];

   // if on left hand edge
   if (node_num_x == 0)
    {

     // first dm0ds0 and dm1ds0

     // get position of right node
     Vector<double> s_macro_R(2);
     macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);

     // compute dm0ds0
     m_gen(1,0) = 0.5*(s_macro_R[0]-s_macro_N[0]);

     // compute dm1ds0
     m_gen(1,1) = 0.5*(s_macro_R[1]-s_macro_N[1]);

     // now d2m0/ds0ds1 and d2m1/ds0ds1

     // if lower left hand corner
     if (node_num_y == 0)
      {

       // get position of upper right node
       Vector<double> s_macro_UR(2);
       macro_coordinate_position(node_num_x+1,node_num_y+1,s_macro_UR);
       
       // get position of upper node
       Vector<double> s_macro_U(2);
       macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);       

       // compute d2m0/ds0ds1
       m_gen(3,0) = 0.5 * (0.5*(s_macro_UR[0] - s_macro_U[0]) - m_gen(1,0));

        // compute d2m1/ds0ds1
       m_gen(3,1) = 0.5 * (0.5*(s_macro_UR[1] - s_macro_U[1]) - m_gen(1,1));
      }

     // else if upper left hand corner
     else if (node_num_y == Nelement[1])
      {

       // get position of lower right node
       Vector<double> s_macro_DR(2);
       macro_coordinate_position(node_num_x+1,node_num_y-1,s_macro_DR);
       
       // get position of lower node
       Vector<double> s_macro_D(2);
       macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);    

       // compute d2m0/ds0ds1
       m_gen(3,0) = 0.5 * (m_gen(1,0) - 0.5*(s_macro_DR[0]-s_macro_D[0]));

       // compute d2m0/ds0ds1
       m_gen(3,1) = 0.5 * (m_gen(1,1) - 0.5*(s_macro_DR[1]-s_macro_D[1]));
      }

     // else left hand edge (not corner)
     else
      {
       
       // get position of lower right node
       Vector<double> s_macro_DR(2);
       macro_coordinate_position(node_num_x+1,node_num_y-1,s_macro_DR);
       
       // get position of lower node
       Vector<double> s_macro_D(2);
       macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);    

       // get position of upper right node
       Vector<double> s_macro_UR(2);
       macro_coordinate_position(node_num_x+1,node_num_y+1,s_macro_UR);
       
       // get position of upper node
       Vector<double> s_macro_U(2);
       macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);  

       // compute d2m0/ds0ds1
       m_gen(3,0) = 0.5*std::min(m_gen(1,0)-0.5*(s_macro_DR[0]-s_macro_D[0]),
                                 0.5*(s_macro_UR[0]-s_macro_U[0])-m_gen(1,0));

       // compute d2m1/ds0ds1
       m_gen(3,1) = 0.5*std::min(m_gen(1,1)-0.5*(s_macro_DR[1]-s_macro_D[1]),
                                 0.5*(s_macro_UR[1]-s_macro_U[1])-m_gen(1,1));
      }
    }

   // if on right hand edge
   else if (node_num_x == Nelement[0])
    {

     // first dm0ds0 and dm1ds0

     // get position of left node
     Vector<double> s_macro_L(2);
     macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);

     // compute dm0ds0
     m_gen(1,0) = 0.5*(s_macro_N[0]-s_macro_L[0]);

     // compute dm1ds0
     m_gen(1,1) = 0.5*(s_macro_N[1]-s_macro_L[1]);

     // now d2m0/ds0ds1 and d2m1/ds0ds1

     // if lower right hand corner
     if (node_num_y == 0)
      {

       // get position of upper left node
       Vector<double> s_macro_UL(2);
       macro_coordinate_position(node_num_x-1,node_num_y+1,s_macro_UL);   
       
       // get position of upper node
       Vector<double> s_macro_U(2);
       macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);       

       // compute d2m0/ds0ds1
       m_gen(3,0) = 0.5 * (0.5*(s_macro_U[0] - s_macro_UL[0]) - m_gen(1,0));

        // compute d2m1/ds0ds1
       m_gen(3,1) = 0.5 * (0.5*(s_macro_U[1] - s_macro_UL[1]) - m_gen(1,1));
      }

     // else if upper right hand corner
     else if (node_num_y == Nelement[1])
      {

       // get position of lower left node
       Vector<double> s_macro_DL(2);
       macro_coordinate_position(node_num_x-1,node_num_y-1,s_macro_DL);
       
       // get position of lower node
       Vector<double> s_macro_D(2);
       macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);    

       // compute d2m0/ds0ds1
       m_gen(3,0) = 0.5 * (m_gen(1,0) - 0.5*(s_macro_D[0]-s_macro_DL[0]));

       // compute d2m0/ds0ds1
       m_gen(3,1) = 0.5 * (m_gen(1,1) - 0.5*(s_macro_D[1]-s_macro_DL[1]));
      }

     // else left hand edge (not corner)
     else
      {
       

       // get position of lower left node
       Vector<double> s_macro_DL(2);
       macro_coordinate_position(node_num_x-1,node_num_y-1,s_macro_DL);
       
       // get position of lower node
       Vector<double> s_macro_D(2);
       macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);    

       // get position of upper left node
       Vector<double> s_macro_UL(2);
       macro_coordinate_position(node_num_x-1,node_num_y+1,s_macro_UL);   
       
       // get position of upper node
       Vector<double> s_macro_U(2);
       macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);  

       // compute d2m0/ds0ds1
       m_gen(3,0) = 0.5*std::min(m_gen(1,0)-0.5*(s_macro_D[0]-s_macro_DL[0]),
                                 0.5*(s_macro_U[0]-s_macro_UL[0])-m_gen(1,0));

       // compute d2m0/ds0ds1
       m_gen(3,1) = 0.5*std::min(m_gen(1,1)-0.5*(s_macro_D[1]-s_macro_DL[1]),
                                 0.5*(s_macro_U[1]-s_macro_UL[1])-m_gen(1,1));
      }
    }
   // 
   else
    {
     // get position of left node
     Vector<double> s_macro_L(2);
     macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);

     // get position of right node
     Vector<double> s_macro_R(2);    
     macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);

     // compute dm0/ds0
     m_gen(1,0) = 0.5*std::min(s_macro_N[0]-s_macro_L[0],
                              s_macro_R[0]-s_macro_N[0]);

     //compute dm1/ds0
     Vector<double> node_space(2);
     node_space[0] = s_macro_N[1]-s_macro_L[1];
     node_space[1] = s_macro_R[1]-s_macro_N[1];
     // set nodal dof
     if (node_space[0] > 0 && node_space[1] > 0)
      m_gen(1,1)  = 0.5 * std::min(node_space[0],node_space[1]);
     else if (node_space[0] < 0 && node_space[1] < 0)
      m_gen(1,1)  = 0.5 * std::max(node_space[0],node_space[1]);
     else
      m_gen(1,1) = 0;
    }

   // if node in lower row
   if (node_num_y == 0)
    {

     // now dm0ds1 and dm1ds1

     // get position of upper node
     Vector<double> s_macro_U(2);
     macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);

     // compute dm0ds1
     m_gen(2,0) = 0.5*(s_macro_U[0]-s_macro_N[0]);

     // compute dm1ds1
     m_gen(2,1) = 0.5*(s_macro_U[1]-s_macro_N[1]);
    
     // and if not corner node
     if (node_num_x > 0 && node_num_x < Nelement[0])
      {

       // now dm0/ds0ds1 and dm1/ds0ds1

       // get position of upper left node
       Vector<double> s_macro_UL(2);
       macro_coordinate_position(node_num_x-1,node_num_y+1,s_macro_UL);       
       
       // get position of left node
       Vector<double> s_macro_L(2);
       macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);
       
       // get position of right node
       Vector<double> s_macro_R(2);
       macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);
       
       // get position of upper right node
       Vector<double> s_macro_UR(2);
       macro_coordinate_position(node_num_x+1,node_num_y+1,s_macro_UR);
       
       // set dm0/ds0ds1
       m_gen(3,0) = 0.5*std::min(m_gen(2,0)-0.5*(s_macro_UL[0]-s_macro_L[0]),
                                 0.5*(s_macro_UR[0]-s_macro_R[0])-m_gen(2,0));

       // set dm1/ds0ds1
       m_gen(3,1) = 0.5*std::min(m_gen(2,1)-0.5*(s_macro_UL[1]-s_macro_L[1]),
                                 0.5*(s_macro_UR[1]-s_macro_R[1])-m_gen(2,1));

      }
    }

   // else if node in upper row
   else if (node_num_y == Nelement[1])
    {

     // now dm0ds1 and dm1ds1

     // get position of lower node
     Vector<double> s_macro_D(2);
     macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);

     // compute dm0ds1
     m_gen(2,0) = 0.5*(s_macro_N[0]-s_macro_D[0]);

     // compute dm1ds1
     m_gen(2,1) = 0.5*(s_macro_N[1]-s_macro_D[1]);
    
     // and if not corner node
     if (node_num_x > 0 && node_num_x < Nelement[0])
      {

       // now dm0/ds0ds1 and dm1/ds0ds1

       // get position of lower left node
       Vector<double> s_macro_DL(2);
       macro_coordinate_position(node_num_x-1,node_num_y-1,s_macro_DL);
       
       // get position of left node
       Vector<double> s_macro_L(2);
       macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);
       
       // get position of right node
       Vector<double> s_macro_R(2);
       macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);
       
       // get position of lower right node
       Vector<double> s_macro_DR(2);
       macro_coordinate_position(node_num_x+1,node_num_y-1,s_macro_DR);
       
       // set dm0/ds0ds1
       m_gen(3,0) = 0.5*std::min(m_gen(2,0)-0.5*(s_macro_L[0]-s_macro_DL[0]),
                                 0.5*(s_macro_R[0]-s_macro_DR[0])-m_gen(2,0));

       // set dm1/ds0ds1
       m_gen(3,1) = 0.5*std::min(m_gen(2,1)-0.5*(s_macro_L[1]-s_macro_DL[1]),
                                 0.5*(s_macro_R[1]-s_macro_DR[1])-m_gen(2,1));
      }
    }
   else
    {
    
     // get position of upper node
     Vector<double> s_macro_U(2);
     macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);

     // get position of lower node
     Vector<double> s_macro_D(2);
     macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);

     //compute dm0/ds1
     Vector<double> node_space(2);
     node_space[0] = s_macro_N[0]-s_macro_D[0];
     node_space[1] = s_macro_U[0]-s_macro_N[0];        
     // set nodal coordinate
     if (node_space[0] > 0 && node_space[1] > 0)
      m_gen(2,0)  = 0.5 * std::min(node_space[0],node_space[1]);
     else if (node_space[0] < 0 && node_space[1] < 0)
      m_gen(2,0)  = 0.5 * std::max(node_space[0],node_space[1]);
     else
      m_gen(2,0) = 0;

     // compute dm1/ds1
     m_gen(2,1) = 0.5*std::min(s_macro_N[1]-s_macro_D[1],
                               s_macro_U[1]-s_macro_N[1]);

     // for interior nodes
     if (node_num_x > 0 && node_num_x < Nelement[0])
      {

       // get position of left node
       Vector<double> s_macro_L(2);
       macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);
       
       // get position of upper left node
       Vector<double> s_macro_UL(2);
       macro_coordinate_position(node_num_x-1,node_num_y+1,s_macro_UL);   
       
       // get position of upper right node
       Vector<double> s_macro_UR(2);
       macro_coordinate_position(node_num_x+1,node_num_y+1,s_macro_UR);
       
       // get position of right node
       Vector<double> s_macro_R(2);
       macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);
       
       // get position of lower right node
       Vector<double> s_macro_DR(2);
       macro_coordinate_position(node_num_x+1,node_num_y-1,s_macro_DR);
       
       // get position of lower left node
       Vector<double> s_macro_DL(2);
       macro_coordinate_position(node_num_x-1,node_num_y-1,s_macro_DL);

       Vector<double> node_space(2);
       // comute dm0/ds0ds1 wrt to node above and below
       node_space[0] = m_gen(1,0) - 0.5*std::min(s_macro_D[0]-s_macro_DL[0],
                                          s_macro_DR[0]-s_macro_D[0]);
       node_space[1] = 0.5*std::min(s_macro_U[0]-s_macro_UL[0],
                              s_macro_UR[0]-s_macro_U[0]) - m_gen(1,0);
       // set nodal dof
       if (node_space[0] > 0 && node_space[1] > 0)
        m_gen(3,0) = 0.5 * std::min(node_space[0],node_space[1]);
       else if (node_space[0] < 0 && node_space[1] < 0)
        m_gen(3,0) = 0.5 * std::max(node_space[0],node_space[1]);
       else
        m_gen(3,0) = 0;  
       
       // comute dm1/ds0ds1 wrt node left and right
       node_space[0] = m_gen(2,1) - 0.5*std::min(s_macro_L[0]-s_macro_DL[0],
                                          s_macro_UL[0]-s_macro_L[0]);
       node_space[1] = 0.5*std::min(s_macro_R[0]-s_macro_DR[0],
                              s_macro_UR[0]-s_macro_R[0]) - m_gen(2,1);
       // set nodal dof
       if (node_space[0] > 0 && node_space[1] > 0)
        m_gen(3,1)  = 0.5 * std::min(node_space[0],node_space[1]);
       else if (node_space[0] < 0 && node_space[1] < 0)
        m_gen(3,1) = 0.5 * std::max(node_space[0],node_space[1]);
       else
        m_gen(3,1) = 0;  
      }
    }
  }


//=============================================================================
/// \short Generic mesh construction function to build the mesh
//=============================================================================
template<class ELEMENT>
void HermiteQuadMesh<ELEMENT>::build_mesh(TimeStepper* time_stepper_pt)
{
 // Mesh can only be built with 2D QHermiteElements.
 MeshChecker::assert_geometric_element<QHermiteElementBase,ELEMENT>(2);

 //Set the number of boundaries
 set_nboundary(4);
 
 //Allocate the store for the elements
 Element_pt.resize(Nelement[0]*Nelement[1]); 
 
 //Allocate the memory for the first element
 Element_pt[0] = new ELEMENT;
 finite_element_pt(0)->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
 
 //Can now allocate the store for the nodes 
 Node_pt.resize((1 + Nelement[0])*(1 + Nelement[1]));
 
 //Set up geometrical data
 unsigned long node_count=0;
 
 //Now assign the topology
 //Boundaries are numbered 0 1 2 3 from the bottom proceeding anticlockwise
 //Pinned value are denoted by an integer value 1
 //Thus if a node is on two boundaries, ORing the values of the 
 //boundary conditions will give the most restrictive case (pinning)
 
 
 // we first create the lowest row of elements
 // ##########################################
 // ##########################################
 
 
 // FIRST ELEMENT - lower left hand corner
 // **************************************
 
 
 // LOWER LEFT HAND NODE
 
 //Set the corner node
 //Allocate memory for the node
 Node_pt[node_count] = 
  finite_element_pt(0)->construct_boundary_node(0,time_stepper_pt);
 
 //Push the node back onto boundaries
 add_boundary_node(0,Node_pt[node_count]);
 add_boundary_node(3,Node_pt[node_count]);
 
 // set the position of the boundary node
 set_position_of_boundary_node
  (0,0,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
 
 //Increment the node number
 node_count++;
 
 // LOWER RIGHT HAND NODE
 
 //Allocate memory for the node
 Node_pt[node_count] = 
  finite_element_pt(0)->construct_boundary_node(1,time_stepper_pt);
 
 //Push the node back onto boundaries
 add_boundary_node(0,Node_pt[node_count]);
 
 //If we only have one column then the RHS node is on the right boundary
 if(Nelement[0]==1)
  {add_boundary_node(1,Node_pt[node_count]);}
 
 // set the position of the node
 set_position_of_boundary_node
  (1,0,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
 
 //Increment the node number
 node_count++;
 
 // UPPER LEFT HAND NODE
 
 //Allocate memory for the node
 Node_pt[node_count] = 
  finite_element_pt(0)->construct_boundary_node(2,time_stepper_pt);
 
 //Push the node back onto boundaries
 add_boundary_node(3,Node_pt[node_count]);
 
 //If we only have one row, then the top node is on the top boundary
 if(Nelement[1]==1)
  {add_boundary_node(2,Node_pt[node_count]);}
 
 // set the position of the boundary node
 set_position_of_boundary_node
  (0,1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
 
 //Increment the node number 
 node_count++;
 
 // UPPER RIGHT NODE
 
 //Allocate the memory for the node
 Node_pt[node_count] = 
  finite_element_pt(0)->construct_node(3,time_stepper_pt);
 
 //If we only have one column then the RHS node is on the right boundary
 if(Nelement[0]==1)
  {add_boundary_node(1,Node_pt[node_count]);}
 //If we only have one row, then the top node is on the top boundary
 if(Nelement[1]==1)
  {add_boundary_node(2,Node_pt[node_count]);}
 
 // if the node is a boundary node
 if (Nelement[0]==1 || Nelement[1]==1)
  {
   // set the position of the boundary node
   set_position_of_boundary_node
    (1,1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
  }
 else
  {
   // set the position of the node
   set_position_of_node(1,1,Node_pt[node_count]);
  }

 //Increment the node number
 node_count++;
 
 // END OF FIRST ELEMENT

 
 // CENTRE OF FIRST ROW OF ELEMENTS
 // *******************************
 
 
 //Now loop over the first row of elements, apart from final element
 for(unsigned j=1;j<(Nelement[0]-1);j++)
  {
   //Allocate memory for new element
   Element_pt[j] = new ELEMENT;
   finite_element_pt(j)->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
   
   // LOWER LEFT NODE
   
   //lower left hand node column of nodes is same as lower right hand node
   //in neighbouring element
   finite_element_pt(j)->node_pt(0) = 
    finite_element_pt(j-1)->node_pt(1);
   
   // LOWER RIGHT NODE
   
   //Allocate memory for the nodes
   Node_pt[node_count] = 
    finite_element_pt(j)->construct_boundary_node(1,time_stepper_pt);
   
   //Push the node back onto boundaries
   add_boundary_node(0,Node_pt[node_count]);
   
   // set the position of the boundary node
   set_position_of_boundary_node
    (j+1,0,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
   
   //Increment the node number
   node_count++;
   
   // UPPER LEFT NODE
   
   //First column of nodes is same as neighbouring element
   finite_element_pt(j)->node_pt(2) = 
    finite_element_pt(j-1)->node_pt(3);
   
   // UPPER RIGHT NODE
   
   //Allocate memory for the nodes
   Node_pt[node_count] = 
    finite_element_pt(j)->construct_node(3,time_stepper_pt);
   
   //If we only have one row, then the top node is on the top boundary
   if(Nelement[0]==1)
    {add_boundary_node(2,Node_pt[node_count]);}
   
   // set the position of the (boundary) node
   if (Nelement[0]==1)
    set_position_of_boundary_node
     (j+1,1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
   else
    set_position_of_node(j+1,1,Node_pt[node_count]);
   
   //Increment the node number
   node_count++;
  } 
 
 // FINAL ELEMENT IN FIRST ROW (lower right corner element)
 // **************************
 
 
 //Only allocate if there is more than one element in the row
 if(Nelement[0] > 1)
  {
   //Allocate memory for element
   Element_pt[Nelement[0]-1] = new ELEMENT;
   finite_element_pt(Nelement[0]-1)->
    set_macro_elem_pt(Domain_pt->macro_element_pt(0));
   
   // LOWER LEFT NODE
   
   //First column of nodes is same as neighbouring element
   finite_element_pt(Nelement[0]-1)->node_pt(0) = 
    finite_element_pt(Nelement[0]-2)->node_pt(1);
   
   // LOWER RIGHT NODE
   
   //If we have an Xperiodic mesh then the final node is the first node in
   //the row
   if(Xperiodic==true)
    {
     //Note that this is periodic in the x-direction   
     finite_element_pt(Nelement[0]-1)->node_pt(1) = 
      finite_element_pt(0)->node_pt(0);
    }
   //Otherwise it's a new final node
   else
    {
     //Allocate memory for the node
     Node_pt[node_count] = 
      finite_element_pt(Nelement[0]-1)->
      construct_boundary_node(1,time_stepper_pt);
    }
   
   //Push the node back onto boundaries
   add_boundary_node(0,Node_pt[node_count]);
   add_boundary_node(1,Node_pt[node_count]);

   // set the position of the boundary node
   set_position_of_boundary_node
    (Nelement[0],0,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));

   // if not periodic mesh
   if(Xperiodic==false)
    {
     node_count++;
    }
   
   // UPPER LEFT NODE
   
   // same as upper right node in element to left
   finite_element_pt(Nelement[0]-1)->node_pt(2) = 
    finite_element_pt(Nelement[0]-2)->node_pt(3);
   
   //If we only have one row, then the top node is on the top boundary
   if(Nelement[1]==1)
    {add_boundary_node(2,Node_pt[node_count]);}   

   // set the position of the (boundary) node
   if (Nelement[1]==1)
    set_position_of_boundary_node
     (Nelement[0]-1,1,
      dynamic_cast<BoundaryNode<Node> *>
      (finite_element_pt(Nelement[0]-2)->node_pt(3)));

   // UPPER RIGHT NODE   
   
   //If we have an Xperiodic mesh then the nodes in the final column are 
   //those in the first column
   if(Xperiodic==true)
    { 
     //Note that these are periodic in the x-direction
     finite_element_pt(Nelement[0]-1)->node_pt(3) = 
      finite_element_pt(0)->node_pt(2);
    }
   //Otherwise we need new nodes for the final column
   else
    {
     //Allocate memory for node
     Node_pt[node_count] = finite_element_pt(Nelement[0]-1)->
      construct_boundary_node(3,time_stepper_pt);
    }   
   
   //If we only have one row, then the top node is on the top boundary
   if(Nelement[1]==1)
    {add_boundary_node(2,Node_pt[node_count]);}
   
   // boundary 1 node
   add_boundary_node(1,Node_pt[node_count]);
   
   // set the position of the boundary node
   set_position_of_boundary_node
    (Nelement[0],1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));

   //Increment the node number
   if(Xperiodic==false)
    {   
     //Increment the node number
     node_count++;
    }
  }
 
 // now create all remaining central rows
 // #####################################
 // #####################################
 
 
 //Loop over remaining element rows in the mesh
 for(unsigned i=1;i<(Nelement[1]-1);i++)
  {
   
   
   // set the first element in the row
   // ********************************
   
   
   //Allocate memory for element
   Element_pt[Nelement[0]*i] = new ELEMENT;
   finite_element_pt(Nelement[0]*i)
    ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));

   //The first row of nodes is copied from the element below
   for(unsigned l=0;l<2;l++)
    {  
     finite_element_pt(Nelement[0]*i)->node_pt(l) = 
      finite_element_pt(Nelement[0]*(i-1))->node_pt(2 + l);
    }
   
   // UPPER LEFT HAND NODE
   
   //Allocate memory for node
   Node_pt[node_count] = finite_element_pt(Nelement[0]*i)->
    construct_boundary_node(2,time_stepper_pt);
   
   //Push the node back onto boundaries
   add_boundary_node(3,Node_pt[node_count]);
   
   // set the position of the boundary node
   set_position_of_boundary_node
    (0,i+1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));

   //Increment the node number
   node_count++;
   
   // UPPER RIGHT HAND NODE
   
   //If we only have one column, the node could be on the boundary
   if(Nelement[0]==1)
    {
     Node_pt[node_count] = finite_element_pt(Nelement[0]*i)->
      construct_boundary_node(3,time_stepper_pt);
    }
   else
    {
     Node_pt[node_count] = 
      finite_element_pt(Nelement[0]*i)->
      construct_node(3,time_stepper_pt);
    }
   
   //If we only have one column then the RHS node is on the 
   //right boundary
   if(Nelement[0]==1)
    {add_boundary_node(1,Node_pt[node_count]);}
   
   // set the position of the (boundary) node
   if (Nelement[0]==1)
    set_position_of_boundary_node
     (1,i+1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
   else
    set_position_of_node(1,i+1,Node_pt[node_count]);

   //Increment the node number
   node_count++;

     
   // loop over central elements in row
   // *********************************
   
   //Now loop over the rest of the elements in the row, apart from the last
   for(unsigned j=1;j<(Nelement[0]-1);j++)
    {
     //Allocate memory for new element
     Element_pt[Nelement[0]*i+j] = new ELEMENT;
     finite_element_pt(Nelement[0]*i+j)->
      set_macro_elem_pt(Domain_pt->macro_element_pt(0));

     // LOWER LEFT AND RIGHT NODES
     
     //The first row is copied from the lower element
     for(unsigned l=0;l<2;l++)
      {  
       finite_element_pt(Nelement[0]*i+j)->node_pt(l) = 
        finite_element_pt(Nelement[0]*(i-1)+j)->node_pt(2 + l);
      }
     
     // UPPER LEFT NODE
     
     //First column is same as neighbouring element
     finite_element_pt(Nelement[0]*i+j)->node_pt(2) = 
      finite_element_pt(Nelement[0]*i+(j-1))->node_pt(3);
     
     // UPPER RIGHT NODE
     
     //Allocate memory for the nodes
     Node_pt[node_count] = 
      finite_element_pt(Nelement[0]*i+j)->
      construct_node(3,time_stepper_pt);

     // set position of node
     set_position_of_node(j+1,i+1,Node_pt[node_count]);
     
     //Increment the node number
     node_count++;
    } 
   
   
   // final element in row
   // ********************
   
   
   //Only if there is more than one column
   if(Nelement[0] > 1)
    {
     //Allocate memory for element
     Element_pt[Nelement[0]*i+Nelement[0]-1] = new ELEMENT;
     finite_element_pt(Nelement[0]*i+Nelement[0]-1)->
      set_macro_elem_pt(Domain_pt->macro_element_pt(0));

     // LOWER LEFT AND RIGHT NODES
     
     //The first row is copied from the lower element
     for(unsigned l=0;l<2;l++)
      {  
       finite_element_pt(Nelement[0]*i+Nelement[0]-1)->node_pt(l) = 
        finite_element_pt(Nelement[0]*(i-1)
                          +Nelement[0]-1)->node_pt(2 + l);
      }
     
     // UPPER LEFT NODE
     
     // First node is same as neighbouring element
     finite_element_pt(Nelement[0]*i+Nelement[0]-1)->node_pt(2) = 
      finite_element_pt(Nelement[0]*i+Nelement[0]-2)->node_pt(3);
     
     // UPPER RIGHT NODE
     
     //If we have an Xperiodic mesh then this node is the same 
     //as the first node
     if(Xperiodic==true)
      {
       //Allocate memory for node, periodic in x-direction
     finite_element_pt(Nelement[0]*i+Nelement[0]-1)->node_pt(3) = 
      finite_element_pt(Nelement[0]*i)->node_pt(2);
      }
     //Otherwise allocate memory for a new node
     else
      {
       Node_pt[node_count] = 
        finite_element_pt(Nelement[0]*i+Nelement[0]-1)->
        construct_boundary_node(3,time_stepper_pt);
      }
     
     //Push the node back onto boundaries
     add_boundary_node(1,Node_pt[node_count]);
     
     // set position of boundary node
     set_position_of_boundary_node
      (Nelement[0],i+1,
       dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));

     //Increment the node number
     node_count++;
    }
  }
 

 // final element row
 // #################
 // #################
 
 
 // only if there is more than one row
 if(Nelement[1] > 1)
  {
   
   // first element in upper row (upper left corner)
   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   
   //Allocate memory for element
   Element_pt[Nelement[0]*(Nelement[1]-1)] = new ELEMENT;
   finite_element_pt(Nelement[0]*(Nelement[1]-1))->
    set_macro_elem_pt(Domain_pt->macro_element_pt(0));
   
   // LOWER LEFT AND LOWER RIGHT NODES
   
   //The first row of nodes is copied from the element below
   for(unsigned l=0;l<2;l++)
    {  
     finite_element_pt(Nelement[0]*(Nelement[1]-1))->node_pt(l) 
      = finite_element_pt(Nelement[0]*(Nelement[1]-2))->node_pt(2 + l);
    }
   
   // UPPER LEFT NODE
   
   //Allocate memory for node
   Node_pt[node_count] = finite_element_pt(Nelement[0]*(Nelement[1]-1))->
    construct_boundary_node(2,time_stepper_pt);
   
   //Push the node back onto boundaries
   add_boundary_node(2,Node_pt[node_count]);
   add_boundary_node(3,Node_pt[node_count]);
   
   // set position of boundary node
   set_position_of_boundary_node
    (0,Nelement[1],dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
   
   //Increment the node number
   node_count++;
   
   // UPPER RIGHT NODE
   
   //Allocate memory for the node
   Node_pt[node_count] = 
    finite_element_pt(Nelement[0]*(Nelement[1]-1))->
    construct_boundary_node(3,time_stepper_pt);
   
   //Push the node back onto boundaries
   add_boundary_node(2,Node_pt[node_count]);
   
   // set position of boundary node
   set_position_of_boundary_node
    (1,Nelement[1],dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));

   //Increment the node number
   node_count++;
   
   
   // loop over central elements of upper element row
   // ***********************************************
   
   
   //Now loop over the rest of the elements in the row, apart from the last
   for(unsigned j=1;j<(Nelement[0]-1);j++)
    {
     
     //Allocate memory for element
     Element_pt[Nelement[0]*(Nelement[1]-1)+j] = new ELEMENT;
     finite_element_pt(Nelement[0]*(Nelement[1]-1)+j)->
      set_macro_elem_pt(Domain_pt->macro_element_pt(0));
     
     // LOWER LEFT AND LOWER RIGHT NODES
     
     //The first row is copied from the lower element
     for(unsigned l=0;l<2;l++)
      {  
       finite_element_pt(Nelement[0]*(Nelement[1]-1)+j)->node_pt(l) = 
        finite_element_pt(Nelement[0]*(Nelement[1]-2)+j)->node_pt(2 + l);
      }
     
     // UPPER LEFT NODE
     
     //First column is same as neighbouring element
     finite_element_pt(Nelement[0]*(Nelement[1]-1)+j)->node_pt(2) 
      = finite_element_pt(Nelement[0]*(Nelement[1]-1)+(j-1))->node_pt(3);
     
     // UPPER RIGHT NODE
     
     //Allocate memory for node
     Node_pt[node_count] = 
      finite_element_pt(Nelement[0]*(Nelement[1]-1)+j)->
      construct_boundary_node(3,time_stepper_pt);
     
     //Push the node back onto boundaries
     add_boundary_node(2,Node_pt[node_count]);
     
     // set position of boundary node
     set_position_of_boundary_node
      (j+1,Nelement[1],
       dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
     
     //Increment the node number   
     node_count++;
     
    } //End of loop over central elements in row
   
   
   // final element in row (upper right corner)
   // *****************************************
   
   //Only if there is more than one column
   if(Nelement[0] > 1)
    {
     
     //Allocate memory for element
     Element_pt[Nelement[0]*(Nelement[1]-1)+Nelement[0]-1] = new ELEMENT;
     finite_element_pt(Nelement[0]*(Nelement[1]-1)+Nelement[0]-1)->
      set_macro_elem_pt(Domain_pt->macro_element_pt(0));
     
     // LOWER LEFT AND LOWER RIGHT NODES
     
     //The first row is copied from the lower element
     for(unsigned l=0;l<2;l++)
      {  
       finite_element_pt(Nelement[0]*(Nelement[1]-1)+Nelement[0]-1)->
        node_pt(l) = finite_element_pt(Nelement[0]*(Nelement[1]-2)+
                                       Nelement[0]-1)->node_pt(2 + l);
      }
     
     // UPPER LEFT NODE
     
     //First column is same as neighbouring element
     finite_element_pt(Nelement[0]*(Nelement[1]-1)+Nelement[0]-1)
      ->node_pt(2) = finite_element_pt(Nelement[0]*(Nelement[1]-1)+
                                       Nelement[0]-2)->node_pt(3);
     
     // UPPER RIGHT NODE
     
     //If we have an Xperiodic mesh, the node must be copied from 
     //the first column
     if(Xperiodic==true)
      {
       //Allocate memory for node, periodic in x-direction     
       finite_element_pt(Nelement[0]*(Nelement[1]-1)
                         +Nelement[0]-1)->node_pt(3) = 
        finite_element_pt(Nelement[0]*(Nelement[1]-1))->node_pt(2);
      }
     //Otherwise, allocate new memory for node
     else
      {
       Node_pt[node_count] = finite_element_pt(Nelement[0]*(Nelement[1]-1)+
                                               Nelement[0]-1)->
        construct_boundary_node(3,time_stepper_pt);
      }
     
     //Push the node back onto boundaries
     add_boundary_node(1,Node_pt[node_count]);
     add_boundary_node(2,Node_pt[node_count]);
     
     // set position of boundary node
     set_position_of_boundary_node
      (Nelement[0],Nelement[1],
       dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));

     //Increment the node number
     node_count++;         
    }
  }
 
 // Setup boundary element lookup schemes
 setup_boundary_element_info(); 
}



//=============================================================================
/// \short Setup lookup schemes which establish which elements are located next
/// to which boundaries (Doc to outfile if it's open).
/// Specific version for HermiteQuadMesh to ensure that the order of the 
/// elements in Boundary_element_pt matches the actual order along the 
/// boundary. This is required when hijacking the BiharmonicElement to apply
/// the BiharmonicFluidBoundaryElement in 
/// BiharmonicFluidProblem::impose_traction_free_edge(...)
//================================================================
template <class ELEMENT>
void HermiteQuadMesh<ELEMENT>::
setup_boundary_element_info(std::ostream &outfile)
{

 bool doc=false;
 if (outfile) doc=true;
 
 // Number of boundaries
 unsigned nbound=nboundary();
 
 // Wipe/allocate storage for arrays
 Boundary_element_pt.clear();
 Face_index_at_boundary.clear();
 Boundary_element_pt.resize(nbound);
 Face_index_at_boundary.resize(nbound);
 
 // Temporary vector of sets of pointers to elements on the boundaries: 
 Vector<Vector<FiniteElement*> > vector_of_boundary_element_pt;
 vector_of_boundary_element_pt.resize(nbound);
 
 // Matrix map for working out the fixed local coord for elements on boundary
 MapMatrixMixed<unsigned,FiniteElement*,Vector<int>* > 
  boundary_identifier;
 

 // Loop over elements
 //-------------------
 unsigned nel=nelement();
 for (unsigned e=0;e<nel;e++)
  {
   // Get pointer to element
   FiniteElement* fe_pt=finite_element_pt(e);
   
   if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
   
   // Only include 2D elements! Some meshes contain interface elements too.
   if (fe_pt->dim()==2)
    {
     // Loop over the element's nodes and find out which boundaries they're on
     // ----------------------------------------------------------------------
     unsigned nnode_1d=fe_pt->nnode_1d();
     
     // Loop over nodes in order
     for (unsigned i0=0;i0<nnode_1d;i0++)
      {
       for (unsigned i1=0;i1<nnode_1d;i1++)
        {
         // Local node number
         unsigned j=i0+i1*nnode_1d;
         
         // Get pointer to vector of boundaries that this
         // node lives on
         std::set<unsigned>* boundaries_pt=0;
         fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);

         // If the node lives on some boundaries....
         if (boundaries_pt!=0)
          {
           // Loop over boundaries
           //unsigned nbound=(*boundaries_pt).size();
           for (std::set<unsigned>::iterator it=boundaries_pt->begin();
                it != boundaries_pt->end();++it)
            {
             // Add pointer to finite element to set for the appropriate 
             // boundary -- storage in set makes sure we don't count elements
             // multiple times
             unsigned temp_size = vector_of_boundary_element_pt[*it].size();
             bool temp_flag = true;
             for (unsigned temp_i = 0; temp_i < temp_size; temp_i++)
              {
               if (vector_of_boundary_element_pt[*it][temp_i] == fe_pt)
                temp_flag = false;
              }
             if (temp_flag)
              vector_of_boundary_element_pt[*it].push_back(fe_pt);
             
             // For the current element/boundary combination, create
             // a vector that stores an indicator which element boundaries
             // the node is located (boundary_identifier=-/+1 for nodes
             // on the left/right boundary; boundary_identifier=-/+2 for nodes
             // on the lower/upper boundary. We determine these indices
             // for all corner nodes of the element and add them to a vector
             // to a vector. This allows us to decide which face of the element
             // coincides with the boundary since the (quad!) element must 
             // have exactly two corner nodes on the boundary.
             if (boundary_identifier(*it,fe_pt)==0)
              {
               boundary_identifier(*it,fe_pt)=new Vector<int>;
              }
             
             // Are we at a corner node?
             if (((i0==0)||(i0==nnode_1d-1))&&((i1==0)||(i1==nnode_1d-1)))
              {
               // Create index to represent position relative to s_0
               (*boundary_identifier(*it,fe_pt)).
                push_back(1*(2*i0/(nnode_1d-1)-1));               
               
               // Create index to represent position relative to s_1
               (*boundary_identifier(*it,fe_pt)).
                push_back(2*(2*i1/(nnode_1d-1)-1));
              } 

            }
          }
         //else
         // {
         //  oomph_info << "...does not live on any boundaries " << std::endl;
         // }
  
        }
      }
    }
    //else
    //{
     //oomph_info << "Element " << e << " does not qualify" << std::endl;
    //}
  }

 
 // Now copy everything across into permanent arrays
 //-------------------------------------------------

 // Note: vector_of_boundary_element_pt contains all elements
 // that have (at least) one corner node on a boundary -- can't copy 
 // them across into Boundary_element_pt
 // yet because some of them might have only one node on the
 // the boundary in which case they don't qualify as 
 // boundary elements!

 // Loop over boundaries
 //---------------------
 for (unsigned i=0;i<nbound;i++)
  {
   // Number of elements on this boundary (currently stored in a set)
   unsigned nel=vector_of_boundary_element_pt[i].size();
    
   // Allocate storage for the coordinate identifiers
   Face_index_at_boundary[i].resize(nel);
   
   // Loop over elements that have at least one corner node on this boundary
   //-----------------------------------------------------------------------
   unsigned e_count=0;
   typedef Vector<FiniteElement*>::iterator IT;
   for (IT it=vector_of_boundary_element_pt[i].begin();
        it!=vector_of_boundary_element_pt[i].end();
        it++)
    {    
     // Recover pointer to element
     FiniteElement* fe_pt=*it;

     // Initialise count for boundary identiers (-2,-1,1,2)
     std::map<int,unsigned> count;
     
     // Loop over coordinates
     for (unsigned ii=0;ii<2;ii++)
      {
       // Loop over upper/lower end of coordinates
       for (int sign=-1;sign<3;sign+=2)
        {
         count[(ii+1)*sign]=0;
        }
      }
     
     // Loop over boundary indicators for this element/boundary
     unsigned n_indicators=(*boundary_identifier(i,fe_pt)).size();
     for (unsigned j=0;j<n_indicators;j++)
      {
       count[(*boundary_identifier(i,fe_pt))[j] ]++;
      }
     delete boundary_identifier(i,fe_pt);

     // Determine the correct boundary indicator by checking that it 
     // occurs twice (since two corner nodes of the element's boundary
     // need to be located on the domain boundary
     int indicator=-10;

     //Check that we're finding exactly one boundary indicator
     bool found=false;

     // Loop over coordinates
     for (unsigned ii=0;ii<2;ii++)
      {
       // Loop over upper/lower end of coordinates
       for (int sign=-1;sign<3;sign+=2)
        {
         if (count[(ii+1)*sign]==2)
          {
           // Check that we haven't found multiple boundaries
           if (found)
            {
             std::string error_message=
              "Trouble: Multiple boundary identifiers!\n";
             error_message += 
              "Elements should only have at most 2 corner ";
             error_message += 
              "nodes on any one boundary.\n";

             throw OomphLibError(
              error_message,
              OOMPH_CURRENT_FUNCTION,
              OOMPH_EXCEPTION_LOCATION);
            }
           found=true;
           indicator=(ii+1)*sign;
          }
        }
      }
     
     // Element has exactly two corner nodes on boundary
     if (found)
      {
       // Add to permanent storage
       Boundary_element_pt[i].push_back(fe_pt);

       // Now convert boundary indicator into information required
       // for FaceElements
       switch (indicator)
        {
         
        case -2:
         
         // s_1 is fixed at -1.0:
         Face_index_at_boundary[i][e_count]= -2;
         break;
         
        case -1:
         
         // s_0 is fixed at -1.0:
         Face_index_at_boundary[i][e_count] = -1;
         break;
         
         
        case 1:
         
         // s_0 is fixed at 1.0:
         Face_index_at_boundary[i][e_count] = 1;
         break;
         
        case 2:
         
         // s_1 is fixed at 1.0:
         Face_index_at_boundary[i][e_count] = 2;
         break;
         
        default:

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

       // Increment counter
       e_count++;
      }
     
    }
  }
 


 // Doc?
 //-----
 if (doc)
  {
   // Loop over boundaries
   for (unsigned i=0;i<nbound;i++)
    {
     unsigned nel=Boundary_element_pt[i].size();
     outfile << "Boundary: " << i
             << " is adjacent to " << nel
             << " elements" << std::endl;
     
     // Loop over elements on given boundary
     for (unsigned e=0;e<nel;e++)
      {
       FiniteElement* fe_pt=Boundary_element_pt[i][e];
       outfile << "Boundary element:" <<  fe_pt
               << " Face index of element along boundary is " 
               <<  Face_index_at_boundary[i][e] << std::endl;
      }
    }
  }
 

 // Lookup scheme has now been setup yet
 Lookup_for_elements_next_boundary_is_setup=true;

}
}

#endif