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// Copyright (C) 2006-2023 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
{
  //=============================================================================
  /// 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));
    }
  }
 
 
  //=============================================================================
  /// 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);
      }
    }
  }
 
 
  //=============================================================================
  /// 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;
      }
    }
  }
 
 
  //=============================================================================
  /// 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();
  }
 
 
  //=============================================================================
  /// 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;
  }
} // namespace oomph
 
#endif