simpler_cubic_mesh.h

#ifndef OOMPH_SIMPLER_CUBIC_MESH_HEADER
#define OOMPH_SIMPLER_CUBIC_MESH_HEADER

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

//Include the OOMPH-LIB header files
#include "../../src/generic/mesh.h"
#include "../../src/generic/matrices.h"
#include "../../src/generic/brick_mesh.h"
#include "../../src/generic/refineable_brick_mesh.h"

namespace oomph
{

  template <class ELEMENT> class SimplerCubicMesh;

  namespace simpler_mesh_helpers
  {

    /// Helper function to compare equality for two floating point
    /// positions.
    inline bool fp_equal(const double& a, const double& b)
    {
      return (std::abs(a - b) < 1e-13);
    }

    /// Helper function to compare equality for two vectors (or vectors of
    /// vectors or ...) of floating point positions.
    template<typename T>
    inline bool fp_equal(const Vector<T>& a,
                         const Vector<T>& b)
    {
      if(a.size() != b.size()) return false;

      const unsigned ni = a.size();
      for(unsigned i=0; i<ni; i++)
        {
          if(! fp_equal(a[i], b[i]))
            {
              return false;
            }
        }

      return true;
    }

    /// Check if neighbour_node is in the correct position to be a node of
    /// element ele_pt. If so insert it, otherwise throw an error.
    inline void insert_as_corresponding_node(FiniteElement* ele_pt,
                                             const Vector<double>& ele_length,
                                             const Vector<double>& ele_x,
                                             Node* neighbour_node)
      {
        // Get position of the neighbouring node
        Vector<double> neighbour_node_x(3);
        neighbour_node->position(neighbour_node_x);

        // Check that the node is inside the element to within fp tol
#ifdef PARANOID
        for(unsigned i=0; i<3; i++)
          {
            if((neighbour_node_x[i] < (ele_x[i] - 1e-13))
               || (neighbour_node_x[i] > (ele_x[i] + ele_length[i] + 1e-13)))
              {
                std::string err = "Node is outside the element:";
                err += "\nnode at " + to_string(neighbour_node_x);
                err += "\nelement at " + to_string(ele_x);
                err += "\nwith size " + to_string(ele_length);
                throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                                    OOMPH_EXCEPTION_LOCATION);
              }
          }
#endif


        for(unsigned nd=0; nd<ele_pt->nnode(); nd++)
          {

            // Calculate position of this node
            Vector<double> s_fraction, node_x(3, 0.0);;
            ele_pt->local_fraction_of_node(nd, s_fraction);
            for(unsigned i=0; i<3; i++)
              {
                node_x[i] = ele_x[i] + ele_length[i]*(s_fraction[i]);
              }

            // If positions are the same then insert it.
            if(fp_equal(node_x, neighbour_node_x))
              {
                // Unless a node is already in the slot (probably from
                // already handling a previous element), in which case
                // check that they are the same node.
                if(ele_pt->node_pt(nd) == 0)
                  {
                    ele_pt->node_pt(nd) = neighbour_node;
                  }
                else
                  {
#ifdef PARANOID
                    if(ele_pt->node_pt(nd) != neighbour_node)
                      {
                        std::string err = "Node mismatch!";
                        err += "\nnode at " + to_string(neighbour_node_x);
                        err += "\nelement at " + to_string(ele_x);
                        err += "\nwith size " + to_string(ele_length);
                        err += "\nold node pt " + to_string(ele_pt->node_pt(nd));
                        err += "\nnew node pt " + to_string(neighbour_node);
                        throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                                            OOMPH_EXCEPTION_LOCATION);
                      }
#endif
                  }

                // succeeded
                return;
              }
           }

        std::string err = "Failed to insert node into element.";
        throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }

    /// Get the neighbouring element to element ijk in coordinate direction
    /// "direction" (in increasing coordinate direction if up is true,
    /// decreasing otherwires). Lookup needs a cubic array of elements.
    inline FiniteElement* neighbouring_element
    (const Vector<unsigned>& ijk,
     bool up, const unsigned& direction,
     const Vector<Vector<Vector <FiniteElement*> > > elements)
    {
#ifdef PARANOID
      if(direction > 2)
        {
          std::string err = "Direction must be 0, 1 or 2.";
          throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif

      // Get number of elements in each direction. Assumed that e.g. all
      // entries of elements[i] have same number of entries, which they
      // should.
      Vector<unsigned> nijk(3, 0);
      nijk[0] = elements.size();
      nijk[1] = elements[0].size();
      nijk[2] = elements[0][0].size();

#ifdef PARANOID
      for(unsigned j=0; j<3; j++)
        {
          if(ijk[j] >= nijk[j])
            {
              std::string err = "ijk out of range in direction " + to_string(j);
              throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
        }
#endif

      // Get next element indices, use int in case it goes negative.
      Vector<int> neighbour_ijk;
      neighbour_ijk.insert(neighbour_ijk.begin(), ijk.begin(), ijk.end());

      if(up) {neighbour_ijk[direction]++;}
      else  {neighbour_ijk[direction]--;}

      // No next element in this direction
      if((neighbour_ijk[direction] >= int(nijk[direction])) ||
         (neighbour_ijk[direction] < 0))
        {
          return 0;
        }
      // All ok
      else
        {
          return elements[neighbour_ijk[0]][neighbour_ijk[1]][neighbour_ijk[2]];
        }
    }

    /// Get a list of nodes which are on the neighbouring face of the
    /// neighbouring element to element ijk in coordinate direction
    /// "direction" (in increasing coordinate direction if up is true,
    /// decreasing otherwires). Lookup needs a cubic array of elements.
    inline Vector<Node*> neighbouring_nodes
    (const Vector<unsigned>& ijk, bool up,
     const unsigned& direction,
     const Vector<Vector<Vector <FiniteElement*> > > elements)
    {
      // Get the element
      FiniteElement* ele_pt = neighbouring_element(ijk, up, direction, elements);

      Vector<Node*> nodes;

      // If element exists
      if(ele_pt != 0)
        {
          // The number of the face index: "face with x constant" = 1, y =
          // 2, z = 3.
          int num = direction + 1;

          // The sign of the face index depends on which direction we are
          // looking from the initial element. If we are looking "upward"
          // (in whichever coordinate) then we want the bottom face of the
          // neighbouring element, which has negative sign.
          int sign = up ? -1 : 1;

          int face_index = sign * num;

          // Get the nodes on face
          for(unsigned i=0; i<ele_pt->nnode_on_face(); i++)
            {
              nodes.push_back
                (ele_pt->node_pt(ele_pt->get_bulk_node_number(face_index, i)));
            }
        }

      return nodes;
    }


    /// In element ele_pt create any nodes whose pointers are currently
    /// null.
    template<class ELEMENT>
    inline Vector<Node*> simple_fill_in_nodes(FiniteElement* ele_pt,
                                              const Vector<double>& ele_x,
                                              const Vector<double>& ele_length,
                                              TimeStepper* time_stepper_pt,
                                              const SimplerCubicMesh<ELEMENT>* mesh_pt)
    {
      Vector<Node*> new_nodes;

      const unsigned n_node = ele_pt->nnode();
      for(unsigned nd=0; nd<n_node; nd++)
        {
          // If node exists then do nothing, otherwise create it.
          if(ele_pt->node_pt(nd) == 0)
            {

              // Calculate node position
              Vector<double> s_fraction, node_x(3, 0.0);;
              ele_pt->local_fraction_of_node(nd, s_fraction);
              for(unsigned i=0; i<3; i++)
                {
                  node_x[i] = ele_x[i] + ele_length[i]*(s_fraction[i]);
                }

              // If position is on boundary make a boundary node, otherwise
              // make a normal node.
              if(mesh_pt->is_on_boundary(node_x))
                {
                  ele_pt->construct_boundary_node(nd, time_stepper_pt);
                }
              else
                {
                  ele_pt->construct_node(nd, time_stepper_pt);
                }

              // Fill in position of node
              for(unsigned i=0; i<3; i++)
                {
                  ele_pt->node_pt(nd)->x(i) = node_x[i];
                }

              // Add to list of new nodes
              new_nodes.push_back(ele_pt->node_pt(nd));
            }
        }

      return new_nodes;
    }

  }

  using namespace simpler_mesh_helpers;


  //=======================================================================
  /// Simpler cubic 3D Brick mesh class.
  //=======================================================================
  template <class ELEMENT>
  class SimplerCubicMesh : public virtual RefineableBrickMesh<ELEMENT>
  {

  public:

    /// \short Constructor: Pass number of elements in the x, y, and z directions,
    /// and the corresponding dimensions. Assume that the back lower left corner
    /// is located at (0,0,0)
    /// Timestepper defaults to Steady.
    SimplerCubicMesh(const unsigned &nx, const unsigned &ny, const unsigned &nz,
                     const double &lx, const double &ly, const double &lz,
                     TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
      Nx(nx), Ny(ny), Nz(nz), Xmin(0.0), Xmax(lx), Ymin(0.0), Ymax(ly),
      Zmin(0.0), Zmax(lz)
    {
      // Mesh can only be built with 3D Qelements.
      MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3);

      //Call the generic build function
      build_mesh(time_stepper_pt);

      // Setup octree forest (for refinement)
      this->setup_octree_forest();
    }

    /// \short Constructor: Pass the number of elements in the x,y and z directions
    /// and the correspoding minimum and maximum values of the coordinates in
    /// each direction
    SimplerCubicMesh(const unsigned &nx, const unsigned &ny, const unsigned &nz,
                     const double &xmin, const double &xmax, const double &ymin,
                     const double &ymax, const double &zmin,const double &zmax,
                     TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
      Nx(nx), Ny(ny), Nz(nz), Xmin(xmin), Xmax(xmax), Ymin(ymin), Ymax(ymax),
      Zmin(zmin), Zmax(zmax)
    {
      // Mesh can only be built with 3D Qelements.
      MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3);

      //Call the generic mesh constructor
      build_mesh(time_stepper_pt);

      // Setup octree forest (for refinement)
      this->setup_octree_forest();
    }


    /// Access function for number of elements in x directions
    const unsigned &nx() const {return Nx;}

    /// Access function for number of elements in y directions
    const unsigned &ny() const {return Ny;}

    /// Access function for number of elements in y directions
    const unsigned &nz() const {return Nx;}

    /// Check if position is on one of the boundaries
    bool is_on_boundary(const Vector<double>& node_x) const
    {
      // Check vs boundaries
      return ((fp_equal(node_x[0], Xmin))
              || (fp_equal(node_x[0], Xmax))
              || (fp_equal(node_x[1], Ymin))
              || (fp_equal(node_x[1], Ymax))
              || (fp_equal(node_x[2], Zmin))
              || (fp_equal(node_x[2], Zmax))
              );
    }

  protected:

    /// Number of elements in x direction
    unsigned Nx;

    /// Number of elements in y direction
    unsigned Ny;

    /// Number of elements in y direction
    unsigned Nz;

    /// Minimum value of x coordinate
    double Xmin;

    /// Maximum value of x coordinate
    double Xmax;

    /// Minimum value of y coordinate
    double Ymin;

    /// Minimum value of y coordinate
    double Ymax;

    /// Minimum value of z coordinate
    double Zmin;

    /// Maximum value of z coordinate
    double Zmax;

    /// Generic mesh construction function: contains all the hard work
    void build_mesh(TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper);

    void build_mesh2(TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper);

    void add_to_boundaries(Node* node_pt)
    {
      // Get location
      Vector<double> x(3, 0.0);
      node_pt->position(x);

      // Boundary numbering the same as in the image for simple cubic mesh
      // on the oomph-lib website.

      // Check vs boundaries
      if(fp_equal(x[0], Xmin)) {this->add_boundary_node(4, node_pt);}
      else if(fp_equal(x[0], Xmax)) {this->add_boundary_node(2, node_pt);}

      if(fp_equal(x[1], Ymin)) {this->add_boundary_node(1, node_pt);}
      else if(fp_equal(x[1], Ymax)) {this->add_boundary_node(3, node_pt);}

      if(fp_equal(x[2], Zmin)) {this->add_boundary_node(0, node_pt);}
      else if(fp_equal(x[2], Zmax)) {this->add_boundary_node(5, node_pt);}
    }

  };


  template <class ELEMENT>
  void SimplerCubicMesh<ELEMENT>::build_mesh(TimeStepper* time_stepper_pt)
    {
      // Mesh can only be built with 3D Qelements.
      MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3);

      //Set the number of boundaries
      this->set_nboundary(6);

      // Get nnode1d using a dummy element
      unsigned nn1d;
      {
        FiniteElement* dummy_pt = new ELEMENT;
        nn1d = dummy_pt->nnode_1d();
        delete dummy_pt; dummy_pt = 0;
      }

      // Calculate the size of each element
      Vector<double> ele_length(3, 0.0);
      ele_length[0] = (Xmax - Xmin)/double(Nx);
      ele_length[1] = (Ymax - Ymin)/double(Ny);
      ele_length[2] = (Zmax - Zmin)/double(Nz);

      // Create storage for elements and element locations. The first is a
      // 3d array Nx x Ny x Nz of null pointers.
      Vector<Vector<Vector <FiniteElement*> > > eles_pt
        (Nx, Vector<Vector<FiniteElement*> >
         (Ny, Vector<FiniteElement*>
          (Nz, 0)));


      // Initialise sizes of class variable vectors so that we aren't
      // reallocating inside the loop.
      this->Element_pt.reserve(Nx*Ny*Nz);
      this->Node_pt.reserve((1 + (nn1d-1)*Nx)*(1 + (nn1d-1)*Ny)*(1 + (nn1d-1)*Nz));


      // Loop over x, y, z directions creating elements and nodes
      for(unsigned ix=0; ix < Nx; ix++)
        {
          for(unsigned iy=0; iy < Ny; iy++)
            {
              for(unsigned iz=0; iz < Nz; iz++)
                {
                  // Create the element itself
                  FiniteElement* ele_pt = new ELEMENT;
                  eles_pt[ix][iy][iz] = ele_pt;

                  // Copy ijk information to vector (c++11 would make this
                  // cleaner...)
                  Vector<unsigned> ijk(3);
                  ijk[0] = ix; ijk[1] = iy; ijk[2] = iz;

                  // Calculate position of the elements "lowest" (in all
                  // directions) corner.
                  Vector<double> ele_corner(3, 0.0);
                  ele_corner[0] = Xmin + (Xmax - Xmin)*(double(ix)/double(Nx));
                  ele_corner[1] = Ymin + (Ymax - Ymin)*(double(iy)/double(Ny));
                  ele_corner[2] = Zmin + (Zmax - Zmin)*(double(iz)/double(Nz));

                  // Lookup any existing nodes from neighbouring elements
                  for(unsigned dir=0; dir<3; dir++)
                    {
                      // Get neighbouring nodes from element "below"
                      Vector<Node*> neighbour_nodes =
                        neighbouring_nodes(ijk, false, dir, eles_pt);

                      // Add nodes from element above
                      {
                        Vector<Node*> top_neighbour_nodes =
                          neighbouring_nodes(ijk, true, dir, eles_pt);

                        neighbour_nodes.insert(neighbour_nodes.end(),
                                               top_neighbour_nodes.begin(),
                                               top_neighbour_nodes.end());
                      }


                      // Brute force search for node slot which has the
                      // correct position and insert node pt. Assuming that
                      // nnode1d is sensibly small this is cheap. May be
                      // slow for spectral elements.
                      const unsigned ni = neighbour_nodes.size();
                      for(unsigned i=0; i<ni; i++)
                        {
                          insert_as_corresponding_node(ele_pt,
                                                       ele_length,
                                                       ele_corner,
                                                       neighbour_nodes[i]);
                        }
                    }

                  // Fill in any remaining nodes, and get pointers to them.
                  Vector<Node*> new_nodes =
                    simple_fill_in_nodes(ele_pt, ele_corner, ele_length,
                                         time_stepper_pt, this);

                  // Store the new element and node pointers in the class
                  // variables
                  this->Element_pt.push_back(ele_pt);
                  this->Node_pt.insert(this->Node_pt.end(), new_nodes.begin(),
                                 new_nodes.end());
                }
            }
        }

      // Now loop over the nodes and set up boundary info based on location
      for(unsigned i=0; i<this->Node_pt.size(); i++)
        {
          this->add_to_boundaries(this->Node_pt[i]);
        }

      // Setup lookup scheme that establishes which elements are located
      // on the mesh boundaries
      this->setup_boundary_element_info();


#ifdef PARANOID
      if(this->Element_pt.size() != Nx*Ny*Nz)
        {
          std::string err = "Wrong Element_pt size.";
          throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }

      // //O(N^2) time debugging code
      // const unsigned ni = Node_pt.size();
      // for(unsigned i=0; i<ni; i++)
      //   {
      //     const unsigned nj = Node_pt.size();
      //     for(unsigned j=0; j<nj; j++)
      //       {
      //         // Don't compare same two nodes twice
      //         if(i < j)
      //           {
      //             // Check no duplicate pointers
      //             if(Node_pt[j] == Node_pt[i])
      //             {
      //               std::string err = "Multiple copies of node pt ";
      //               err += to_string(Node_pt[i]);
      //               err += " " + to_string(Node_pt[j]);
      //               throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
      //                                   OOMPH_EXCEPTION_LOCATION);
      //             }

      //             // Check no duplicate positions
      //             Vector<double> xi(3), xj(3);
      //             Node_pt[j]->position(xj);
      //             Node_pt[i]->position(xi);
      //             if(fp_equal(xi, xj))
      //               {
      //                 std::string err = "Two nodes in same place: ";
      //                 err += to_string(Node_pt[i]);
      //                 err += " " + to_string(Node_pt[j]);
      //                 err += "\nat position " + to_string(xi);
      //                 err += "\nand " + to_string(xj);
      //                 throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
      //                                     OOMPH_EXCEPTION_LOCATION);
      //               }
      //           }
      //       }
      //   }

      if(this->Node_pt.size() != (1 + (nn1d-1)*Nx)*(1 + (nn1d-1)*Ny)*(1 + (nn1d-1)*Nz))
        {
          std::string err = "Wrong Node_pt size";
          throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif

    }

}

#endif