refineable_line_element.cc

Go to the documentation of this file.

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

#include<algorithm>

// oomph-lib headers
#include "mesh.h"
#include "algebraic_elements.h"
#include "macro_element_node_update_element.h"
#include "refineable_line_element.h"

namespace oomph
{

 //==========================================================================
 /// Setup static matrix for coincidence between son nodal points and father
 /// boundaries:
 ///
 /// Father_bound[nnode_1d](nnode_son,son_type) = {L/R/OMEGA}
 ///
 /// so that node nnode_son in element of type son_type lies on boundary/
 /// vertex Father_bound[nnode_1d](nnode_son,son_type) in its father
 /// element. If the node doesn't lie on a boundary the value is OMEGA.
 //==========================================================================
 void RefineableQElement<1>::setup_father_bounds()
 {
  using namespace BinaryTreeNames;
  
  // Find the number of nodes along a 1D edge (which is the number of nodes
  // in the element for a 1D element!)
  const unsigned n_node = nnode_1d();
  
  // Allocate space for the boundary information
  Father_bound[n_node].resize(n_node,2);
  
  // Initialise: By default points are not on the boundary
  for(unsigned n=0;n<n_node;n++)
   {
    for(unsigned ison=0;ison<2;ison++)
     {
      Father_bound[n_node](n,ison) = Tree::OMEGA;
     }
   }
  
  // Left-hand son:
  // --------------

  // L node (0) is the L node of the parent
  Father_bound[n_node](0,L) = L;

  // Other boundary is in the interior
  
  // Right-hand son:
  // ---------------

  // R node (n_node-1) is the R node of the parent
  Father_bound[n_node](n_node-1,R) = R;

  // Other boundary is in the interior
 }
 
 //==========================================================================
 /// If a neighbouring element has already created a node at a position
 /// corresponding to the local fractional position within the present
 /// element, s_fraction, return a pointer to that node. If not, return
 /// NULL (0). If the node is periodic the flag is_periodic will be true.
 //==========================================================================
 Node* RefineableQElement<1>::
 node_created_by_neighbour(const Vector<double> &s_fraction, 
                           bool &is_periodic) 
 {  
  using namespace BinaryTreeNames;

  // Initialise edge of current element on which node lies
  int edge_in_current = OMEGA;
  
  // Determine the edge of the current element on which the node lies
  if(s_fraction[0]==0.0) { edge_in_current = L; }
  if(s_fraction[0]==1.0) { edge_in_current = R; }
  
  // If the node does not lie on an edge then there is no neighbour:
  // return NULL
  if(edge_in_current==OMEGA) { return 0; }

  // Allocate storage for edge in neighbouring element
  int edge_in_neighbour;

  // Allocate storage for difference in size between current and
  // neighbouring element
  int diff_level;

  // Allocate storage for local coordinate of node in neighbouring tree
  Vector<double> s_in_neighbour(1);
  
  // Allocate storage for flag indicating if the node is not in the same
  // binary tree
  bool in_neighbouring_tree;

  // Allocate storage for the pointer to the neighbouring element
  // (using its binary tree representation)
  BinaryTree* neighbour_pt;

  // Find pointer to neighbouring element along the edge in question and
  // calculate the local coordinate of the node within that element
  // s_in_neighbour
  neighbour_pt = binary_tree_pt()->
   gteq_edge_neighbour(edge_in_current,s_in_neighbour,edge_in_neighbour,
                       diff_level,in_neighbouring_tree);
  
  // If a neighbour exists...
  if(neighbour_pt!=0)
   {
    // ...check whether its nodes have been created yet
    if(neighbour_pt->object_pt()->nodes_built())
     {
      // If they have, find the node in question in the neighbour
      Node* neighbour_node_pt = 
       neighbour_pt->object_pt()->get_node_at_local_coordinate(s_in_neighbour);
      
      // If there is no node at this position, there is a problem, since in
      // a 1D element (whose nodes have already been built) there should
      // ALWAYS be a node at each edge of the element.
      if(neighbour_node_pt==0)
       {
        std::string error_message =
         "Problem: an element claims to have had its nodes built, yet\n";
        error_message +=
         "it is missing (a least) a node at its edge.\n";
        throw OomphLibError(error_message,
                            OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
       }
      // Otherwise, carry on
      else
       {
        // Now work out whether it's a periodic boundary (this is only
        // (possible if we have moved into a neighbouring tree)
        if(in_neighbouring_tree)
         {
          // Return whether the neighbour is periodic 
          is_periodic = 
           binary_tree_pt()->root_pt()->is_neighbour_periodic(edge_in_current);
         }
        // Return the pointer to the neighbouring node
        return neighbour_node_pt;
       }
     }
   }
  // Node not found, return null
  return 0;
 }
 
 //==========================================================================
 /// Build the element by doing the following:
 /// - Give it nodal positions (by establishing the pointers to its nodes).
 ///   In the process create new nodes where required (i.e. if they don't
 ///   already exist in the father element and have not already been created
 ///   whilst building new neighbouring elements). Node building involves
 ///   the following steps:
 ///   - Get the nodal position from the father element.
 ///   - Establish the time-history of the newly created nodal point
 ///     (its coordinates and previous values) consistent with the father's
 ///     history.
 ///   - Add the new node to the mesh itself.
 ///   - Doc newly created nodes in "new_nodes.dat" stored in the directory
 ///     of the DocInfo object (only if it's open!).
 ///   - NOTE: Unlike in higher-dimensional elements, in 1D it is impossible
 ///     for newly-created nodes to be on a mesh boundary, since any boundary
 ///     nodes must exist in the initial (coarse) mesh. Therefore it is not
 ///     necessary to add any nodes to the mesh's boundary node storage
 ///     schemes, and we always create normal "bulk" nodes.
 /// - Once the element has a full complement of nodes, excute the element-
 ///   specific further_build() (empty by default -- must be overloaded for
 ///   specific elements). This deals with any build operations that are not
 ///   included in the generic process outlined above. For instance, in
 ///   Crouzeix Raviart elements we need to initialise the internal pressure
 ///   pressure values in a manner consistent with the pressure distribution
 ///   in the father element.
 //==========================================================================
 void RefineableQElement<1>::build(Mesh*& mesh_pt, 
                                   Vector<Node*>& new_node_pt,
                                   bool& was_already_built,
                                   std::ofstream &new_nodes_file)
 {
  using namespace BinaryTreeNames;
  
  // Find the number of nodes along a 1D edge (which is the number of nodes
  // in the element for a 1D element!)
  const unsigned n_node = nnode_1d();
  
  // Check whether static father_bound needs to be created
  if(Father_bound[n_node].nrow()==0) { setup_father_bounds(); }
  
  // Pointer to the current element's father (in binary tree impersonation)
  BinaryTree* father_pt
   = dynamic_cast<BinaryTree*>(binary_tree_pt()->father_pt());
  
  // What type of son is the current element?
  // Ask its binary tree representation...
  const int son_type = Tree_pt->son_type();
  
  // Has somebody built the current element already?
  // Check this by determining whether or not the first node has been built
  
  // If this element has not already been built...
  if (!nodes_built())
   {
#ifdef PARANOID
    if (father_pt==0)
     {
      std::string error_message =
       "Something fishy here: I have no father and yet \n";
      error_message +=
       "I have no nodes. Who has created me then?!\n";
      
      throw OomphLibError(error_message,
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
#endif
    
    // Indicate status
    was_already_built = false;
    
    // Return pointer to father element
    RefineableQElement<1>* father_el_pt
     = dynamic_cast<RefineableQElement<1>*>(father_pt->object_pt());
    
    // Timestepper should be the same for all nodes in father element.
    // Use it create timesteppers for new nodes.
    TimeStepper* time_stepper_pt=father_el_pt->node_pt(0)->time_stepper_pt();
    
    // Determine number of history values (including present)
    const unsigned ntstorage = time_stepper_pt->ntstorage();
    
    // Currently we can't handle the case of generalised coordinates
    // since we haven't established how they should be interpolated.
    // Buffer this case:
    if (father_el_pt->node_pt(0)->nposition_type()!=1)
     {
      throw OomphLibError(
       "Can't handle generalised nodal positions (yet).",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
    
    // Set up vertex coordinates in the father element:
    // ------------------------------------------------    

    // Allocate storage for the vertex coordinates s_left and s_right of
    // the current element as viewed by this element's father (i.e. s_left[0]
    // stores the local coordinate within the father element at which the
    // node on the current element's left-hand edge is located. Likewise,
    // s_right[0] stores the local coordinate within the father element at
    // which the node on the current element's right-hand edge is located).
    Vector<double> s_left(1);
    Vector<double> s_right(1);

    // In order to set up the vertex coordinates we need to know which
    // type of son the current element is
    switch(son_type)
     {
     case L:
      s_left[0] = -1.0;
      s_right[0] = 0.0;
      break;
      
     case R:
      s_left[0] = 0.0;
      s_right[0] = 1.0;
      break;
     }
    
    // Pass macro element pointer on to sons and set coordinates in
    // macro element
    // hierher why can I see this?
    if(father_el_pt->Macro_elem_pt!=0)
     {
      set_macro_elem_pt(father_el_pt->Macro_elem_pt);
      
      s_macro_ll(0)=      father_el_pt->s_macro_ll(0)+
       0.5*(s_left[0]+1.0)*(father_el_pt->s_macro_ur(0)-
                            father_el_pt->s_macro_ll(0));
      s_macro_ur(0)=      father_el_pt->s_macro_ll(0)+
       0.5*(s_right[0]+1.0)*(father_el_pt->s_macro_ur(0)-
                             father_el_pt->s_macro_ll(0));
     }

    // If the father element hasn't been generated yet, we're stuck...
    if(father_el_pt->node_pt(0)==0)
     {
      throw OomphLibError(
       "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
    // Otherwise, carry on
    else
     {
      // Allocate storage for the location of a node in terms of the local
      // coordinate of the father element
      Vector<double> s_in_father(1);
      
      // Allocate storage for the fractional position (in the current
      // element) of the node in the direction of s[0]
      Vector<double> s_fraction(1);

      // Loop over all nodes in the element
      for(unsigned n=0;n<n_node;n++)
       {
        // Get the fractional position (in the current element) of the node
        // in the direction of s[0]
        s_fraction[0] = local_one_d_fraction_of_node(n,0);
        
        // Evaluate the local coordinate of the node in the father element
        s_in_father[0] = s_left[0] + (s_right[0]-s_left[0])*s_fraction[0];
        
        // Initialise flag: So far, this node hasn't been built or copied yet
        bool node_done = false;
        
        // Get the pointer to the node in the father (returns NULL if there
        // is no node at this position)
        Node* created_node_pt
         = father_el_pt->get_node_at_local_coordinate(s_in_father);
        
        // Does this node already exist in father element?
        // -----------------------------------------------

        // If it does...
        if(created_node_pt!=0)
         {
          // ...copy node across
          node_pt(n) = created_node_pt;
          
          // Make sure that we update the values at the node so that they
          // are consistent with the present representation. This is only
          // needed for mixed interpolation where the value at the father
          // could now become active.

          // Loop over all history values
          for(unsigned t=0;t<ntstorage;t++)
           {
            // Get values from father element
            // Note: get_interpolated_values() sets Vector size itself
            Vector<double> prev_values;
            father_el_pt->get_interpolated_values(t,s_in_father,prev_values);

            // Find the minimum number of values (either those stored at the
            // node, or those returned by the function)
            unsigned n_val_at_node = created_node_pt->nvalue();
            unsigned n_val_from_function = prev_values.size(); 

            // Use the ternary conditional operator here
            unsigned n_var = n_val_at_node < n_val_from_function ?
             n_val_at_node : n_val_from_function;

            // Assign the values that we can
            for(unsigned k=0;k<n_var;k++)
             {
              created_node_pt->set_value(t,k,prev_values[k]);
             }
           }
          
          // Indicate that node has been created by copy
          node_done = true;
         }

        // Node does not exist in father element but might already
        // -------------------------------------------------------
        // have been created by neighbouring elements
        // ------------------------------------------

        else
         {
          // Was the node created by one of its neighbours? Whether or not
          // the node lies on an edge can be calculated from the fractional
          // position.
          bool is_periodic = false;
          created_node_pt = node_created_by_neighbour(s_fraction,is_periodic);
          
          // If the node was so created...
          if(created_node_pt!=0)
           {
            // ...assign the pointer
            node_pt(n) = created_node_pt;
            
            // Indicate that node has been created by copy
            node_done = true; 
            
            // In a 1D mesh there is no way that a periodic node (which must
            // be on a boundary) can exist without being part of the initial
            // (coarse) mesh. Therefore issue an error message if
            // node_created_by_neighbour(...) returns `is_periodic==true'.
#ifdef PARANOID
            if(is_periodic)
             {
              std::string error_message =
               "node_created_by_neighbour returns a node which claims\n";
              error_message += 
               "to be periodic. In a 1D mesh any periodic nodes must exist\n";
              error_message +=
               "in the initial (coarse) mesh.";
              
              throw OomphLibError(error_message,
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
             }
#endif
           }
         }
        
        // Node does not exist in father element or neighbouring element
        // -------------------------------------------------------------
        
        // If the node has not been built anywhere ---> build it here
        if(!node_done)
         {
          // In a 1D mesh any node which lies on the boundary must exist in
          // the initial (coarse) mesh, so any newly-built nodes cannot be
          // boundary nodes. Therefore we always create a normal "bulk" node.
          
          // Create node and set the pointer to it from the element 
          created_node_pt = construct_node(n,time_stepper_pt);
          
          // Add to vector of new nodes
          new_node_pt.push_back(created_node_pt);
          
          // Now we set the position and values at the newly created node.
          // In the first instance use macro element or FE representation
          // to create past and present nodal positions.
          // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC ELEMENTS AS NOT
          // ALL OF THEM NECESSARILY IMPLEMENT NONTRIVIAL NODE UPDATE
          // FUNCTIONS. CALLING THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL
          // LEAVE THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
          // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL NOT ASSIGN SENSIBLE
          // INITIAL POSITIONS!)
          
          // Loop over history values
          for(unsigned t=0;t<ntstorage;t++)
           {
            // Allocate storage for the previous position of the node
            Vector<double> x_prev(1);
            
            // Get position from father element -- this uses the macro
            // element representation if appropriate.
            father_el_pt->get_x(t,s_in_father,x_prev);
            
            // Set the previous position of the new node
            created_node_pt->x(t,0) = x_prev[0];

            // Allocate storage for the previous values at the node
            // NOTE: the size of this vector is equal to the number of values
            // (e.g. 3 velocity components and 1 pressure, say)
            // associated with each node and NOT the number of history values
            // which the node stores!
            Vector<double> prev_values;
            
            // Get values from father element
            // Note: get_interpolated_values() sets Vector size itself.
            father_el_pt->get_interpolated_values(t,s_in_father,prev_values);
            
            // Determine the number of values at the new node
            const unsigned n_value = created_node_pt->nvalue();
            
            // Loop over all values and set the previous values
            for(unsigned v=0;v<n_value;v++)
             {
              created_node_pt->set_value(t,v,prev_values[v]);
             }
           } // End of loop over history values
          
          // Add new node to mesh
          mesh_pt->add_node_pt(created_node_pt);
         
         } // End of case where we build the node ourselves
        
        // Check if the element is an algebraic element
        AlgebraicElementBase* alg_el_pt=
         dynamic_cast<AlgebraicElementBase*>(this);             
        
        // If it is, set up node position (past and present) from algebraic
        // node position (past and present) from algebraic node update
        // node update function. This over-writes previous assingments that
        // were made based on the macro-element/FE representation.
        // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE NODE IS A
        // MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN THE SAME NODAL
        // POSITIONS BUT WE NEED TO ADD THE REMESH INFO FOR *ALL* ROOT
        // ELEMENTS!
        if (alg_el_pt!=0)
         {
          // Build algebraic node update info for new node. This sets up
          // the node update data for all node update functions that are
          // shared by all nodes in the father element.
          alg_el_pt->setup_algebraic_node_update(node_pt(n),s_in_father,
                                                 father_el_pt);
         }
        
        // If we have built the node and we are documenting our progress,
        // write the (hopefully consistent position) to the outputfile
        if((!node_done) && (new_nodes_file.is_open()))
         {
          new_nodes_file << node_pt(n)->x(0) << std::endl;
         } 
       } // End of loop over all nodes in element
      

      // If the element is a MacroElementNodeUpdateElement, set the update
      // parameters for the current element's nodes -- all this needs is
      // the vector of (pointers to the) geometric objects that affect the
      // MacroElement-based node update. This is the same as that in the
      // father element
      MacroElementNodeUpdateElementBase* father_m_el_pt
       = dynamic_cast<MacroElementNodeUpdateElementBase*>(father_el_pt);
      if(father_m_el_pt!=0)
       {
        // Get Vector of geometric objects from father (construct Vector
        // via copy operation)
        Vector<GeomObject*> geom_object_pt(father_m_el_pt->geom_object_pt());
        
        // Cast current element to MacroElementNodeUpdateElement:
        MacroElementNodeUpdateElementBase* m_el_pt
         = dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
        
#ifdef PARANOID
        if (m_el_pt==0)
         {
          std::string error_message =
           "Failed to cast to MacroElementNodeUpdateElementBase*\n";
          error_message += 
           "Strange -- if the father is a MacroElementNodeUpdateElement\n";
          error_message += "the son should be too....\n";
          
          throw OomphLibError(error_message,
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
         }
#endif
        // Build update info by passing Vector of geometric objects:
        // This sets the current element to be the update element for all
        // of the element's nodes. This is reversed if the element is
        // ever un-refined in the father element's rebuild_from_sons()
        // function which overwrites this assignment to avoid nasty
        // segmentation faults that occur when a node tries to update
        // itself via an element that no longer exists...
        m_el_pt->set_node_update_info(geom_object_pt);
       }
      
#ifdef OOMPH_HAS_MPI
      // Pass on non-halo proc ID
      Non_halo_proc_ID=
       tree_pt()->father_pt()->object_pt()->non_halo_proc_ID();
#endif
      
      // Is the new element an ElementWithMovingNodes? 
      ElementWithMovingNodes* aux_el_pt
       = dynamic_cast<ElementWithMovingNodes*>(this);
      
      // Pass down the information re the method for the evaluation
      // of the shape derivatives
      if(aux_el_pt!=0)
       {
        ElementWithMovingNodes* aux_father_el_pt
         = dynamic_cast<ElementWithMovingNodes*>(father_el_pt);
        
#ifdef PARANOID
        if(aux_father_el_pt==0)
         {
          std::string error_message =
           "Failed to cast to ElementWithMovingNodes*\n";
          error_message += 
           "Strange -- if the son is a ElementWithMovingNodes\n";
          error_message += "the father should be too....\n";
          
          throw OomphLibError(error_message,
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
         }
#endif
      
        //If evaluating the residuals by finite differences in the father
        //continue to do so in the child
        if(aux_father_el_pt
           ->are_dresidual_dnodal_coordinates_always_evaluated_by_fd())
         {
          aux_el_pt->
           enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
         }
        
        aux_el_pt->method_for_shape_derivs()=
         aux_father_el_pt->method_for_shape_derivs();
        
        //If bypassing the evaluation of fill_in_jacobian_from_geometric_data
        //continue to do so
        if(aux_father_el_pt
           ->is_fill_in_jacobian_from_geometric_data_bypassed())
         {
          aux_el_pt->enable_bypass_fill_in_jacobian_from_geometric_data();
         }
       }
      
      
      // Now do further build (if any)
      further_build();
      
     } // Sanity check: Father element has been generated
    
   } // End of if element has not already been built
  
  // If element has already been built, say so
  else
   {
    was_already_built = true;
   }
 }
 
 //==========================================================================
 ///  Print corner nodes, use colour (default "BLACK")
 //==========================================================================
 void RefineableQElement<1>::output_corners(
  std::ostream& outfile,const std::string& colour) const
 {
  // Allocate storage for local coordinate s
  Vector<double> s(1);

  // Allocate storage for global coordinate of element vertex
  Vector<double> vertex(1);
  
  outfile << "ZONE I=2,J=2, C=" << colour << std::endl;
  
  // Left-hand vertex
  s[0]=-1.0;
  get_x(s,vertex);
  outfile << vertex[0]  << " " << Number << std::endl;
  
  // Right-hand vertex
  s[0]=1.0;
  get_x(s,vertex);
  outfile << vertex[0]  << " " << Number << std::endl;
  
  outfile << "TEXT  CS = GRID, X = " << vertex[0] <<
   ", HU = GRID, H = 0.01, AN = MIDCENTER, T =\"" 
          << Number << "\"" << std::endl;      
 }
 
 //==========================================================================
 /// Check inter-element continuity of 
 /// - nodal positions
 /// - (nodally) interpolated function values
 //==========================================================================
 void RefineableQElement<1>::check_integrity(double& max_error)
 {
  using namespace BinaryTreeNames;
  
  // Calculate number of nodes
  const unsigned n_node = nnode_1d();
  
  // Number of timesteps (including present) for which continuity is 
  // to be checked
  const unsigned n_time=1;

  // Initialise errors
  max_error=0.0;
  Vector<double> max_error_x(1,0.0);
  double max_error_val = 0.0;
  
  // Initialise vector of element edges in the current element
  Vector<int> edge_in_current(2);
  edge_in_current[0] = L;
  edge_in_current[1] = R;
  
  // Loop over both edges
  for(unsigned edge_counter=0;edge_counter<2;edge_counter++)
   {
    // Allocate storage for the fractional position (in the current
    // element) of the node in the direction of s[0]
    Vector<double> s_fraction(1);

    // Allocate storage for the location of a node in terms of the local
    // coordinate of the current element
    Vector<double> s_in_current(1);

    // Allocate storage for the location of a node in terms of the local
    // coordinate of the neighbouring element
    Vector<double> s_in_neighbour(1);

    // Allocate storage for edge in neighbouring element
    int edge_in_neighbour;

    // Allocate storage for difference in size between current and
    // neighbouring element
    int diff_level;

    // Allocate storage for flag indicating if the node is not in the same
    // binary tree
    bool in_neighbouring_tree;

    // Calculate the local coordinate and the local coordinate as viewed
    // from the neighbour
    
    // Allocate storage for the pointer to the neighbouring element
    // (using its binary tree representation)
    BinaryTree* neighbour_pt;

    // Find pointer to neighbouring element along the edge in question and
    // calculate the local coordinate of the node within that element,
    // s_in_neighbour
    neighbour_pt = binary_tree_pt()->
     gteq_edge_neighbour(edge_in_current[edge_counter],s_in_neighbour,
                         edge_in_neighbour,diff_level,in_neighbouring_tree);

    // If neighbour exists and has existing nodes
    if((neighbour_pt!=0) && (neighbour_pt->object_pt()->nodes_built()))
     {
      // Need to exclude periodic nodes from this check.
      // There are only periodic nodes if we are in a neighbouring tree...
      bool is_periodic = false;
      if(in_neighbouring_tree)
       {
        // Is it periodic?
        is_periodic = tree_pt()->root_pt()
         ->is_neighbour_periodic(edge_in_current[edge_counter]);
       }
      
      // Allocate storage for pointer to the local node
      Node* local_node_pt = 0;
      
      switch(edge_counter)
       {
       case 0:
        // Local fraction of node
        s_fraction[0] = 0.0; 
        // Get pointer to local node
        local_node_pt = node_pt(0);
        break;
        
       case 1:
        // Local fraction of node
        s_fraction[0] = 1.0; 
        // Get pointer to local node
        local_node_pt = node_pt(n_node-1);
        break;
       }
      
      // Evaluate the local coordinate of the node in the current element
      s_in_current[0] = -1.0 + 2.0*s_fraction[0];

      // NOTE: We have already calculated the local coordinate of the node
      // in the neighbouring element above when calling gteq_edge_neighbour()
      
      // Loop over timesteps
      for(unsigned t=0;t<n_time;t++)
       {
        // Allocate storage for the nodal position in the neighbouring element
        Vector<double> x_in_neighbour(1);

        // Get the nodal position from the neighbouring element
        neighbour_pt->object_pt()->interpolated_x(t,s_in_neighbour,
                                                  x_in_neighbour);
        
        // Check error only if the node is NOT periodic
        if(is_periodic==false)
         {
          // Find the spatial error
          double err = std::fabs(local_node_pt->x(t,0) - x_in_neighbour[0]);
          
          // If it's bigger than our tolerance, say so
          if(err>1e-9)
           {
            oomph_info << "errx " << err << " " << t << " " 
                       << local_node_pt->x(t,0) 
                       << " " <<  x_in_neighbour[0]<< std::endl;
            
            oomph_info << "at " <<  local_node_pt->x(0) << std::endl;
           }
          
          // If it's bigger than the previous max error, it is the
          // new max error!
          if(err>max_error_x[0]) { max_error_x[0]=err; }
         }
        // Allocate storage for the nodal values in the neighbouring element
        Vector<double> values_in_neighbour;        

        // Get the values from neighbouring element. NOTE: Number of values
        // gets set by routine (because in general we don't know how many
        // interpolated values a certain element has.
        neighbour_pt->object_pt()
         ->get_interpolated_values(t,s_in_neighbour,values_in_neighbour);

        // Allocate storage for the nodal values in the current element
        Vector<double> values_in_current;

        // Get the values in current element
        get_interpolated_values(t,s_in_current,values_in_current);
        
        // Now figure out how many continuously interpolated values there are
        const unsigned num_val
         = neighbour_pt->object_pt()->ncont_interpolated_values();
        
        // Check error
        for(unsigned ival=0;ival<num_val;ival++)
         {
          // Find the spatial error
          double err = std::fabs(values_in_current[ival]
                                - values_in_neighbour[ival]);

          // If it's bigger than our tolerance, say so           
          if (err>1.0e-10)
           {
            oomph_info <<  local_node_pt->x(0) << "\n# " 
                       << "erru " << err << " " << ival << " " 
                       << get_node_number(local_node_pt) << " "
                       << values_in_current[ival]
                       << " " << values_in_neighbour[ival] << std::endl;
           }

          // If it's bigger than the previous max error, it is the
          // new max error!       
          if (err>max_error_val) { max_error_val=err; }          
         }
       } // End of loop over timesteps
     } // End of if neighbour exists and has existing nodes
   } // End of loop over edges
  
  // Update max_error if necessary
  max_error = max_error_x[0];
  if(max_error_val>max_error) { max_error = max_error_val; }
  
  // Output max_error information
  if (max_error>1e-9)
   {
    oomph_info << "\n#------------------------------------ \n#Max error " ;
    oomph_info << max_error_x[0] 
               << " " << max_error_val << std::endl;
    oomph_info << "#------------------------------------ \n " << std::endl; 
   }
 }
 
 //==========================================================================
 /// Static matrix for coincidence between son nodal points and father
 /// boundaries
 //==========================================================================
 std::map<unsigned,DenseMatrix<int> > RefineableQElement<1>::Father_bound;

} // End of namespace