rectangular_quadmesh.template.h

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-2022 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//====================================================================
// Header file for a relatively simple Quad Meshe
#ifndef OOMPH_RECTANGULAR_QUADMESH_HEADER
#define OOMPH_RECTANGULAR_QUADMESH_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB headers
#include "../generic/mesh.h"
#include "../generic/quad_mesh.h"
#include "../generic/refineable_quad_mesh.h"
 
namespace oomph
{
  //==========================================================================
  //===========================================================================
  template<class ELEMENT>
  class RectangularQuadMesh : public virtual QuadMeshBase
  {
  protected:
    // Mesh variables
    unsigned Nx;
    unsigned Ny;
    unsigned Np;
 
    double Xmin;
    double Xmax;
 
    double Ymin;
    double Ymax;
 
    bool Xperiodic;
 
    void build_mesh(TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper);
 
    RectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& xmin,
      const double& xmax,
      const double& ymin,
      const double& ymax,
      const bool& periodic_in_x,
      const bool& build,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : Nx(nx),
        Ny(ny),
        Xmin(xmin),
        Xmax(xmax),
        Ymin(ymin),
        Ymax(ymax),
        Xperiodic(periodic_in_x)
    {
      if (build)
      {
        // Call the generic mesh constructor
        build_mesh(time_stepper_pt);
      }
    }
 
  public:
    RectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : Nx(nx),
        Ny(ny),
        Xmin(0.0),
        Xmax(lx),
        Ymin(0.0),
        Ymax(ly),
        Xperiodic(false)
    {
      // Mesh can only be built with 2D Qelements.
      MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
 
      // Call the generic mesh constructor
      build_mesh(time_stepper_pt);
    }
 
    RectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& xmin,
      const double& xmax,
      const double& ymin,
      const double& ymax,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : Nx(nx),
        Ny(ny),
        Xmin(xmin),
        Xmax(xmax),
        Ymin(ymin),
        Ymax(ymax),
        Xperiodic(false)
    {
      // Mesh can only be built with 2D Qelements.
      MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
 
      // Call the generic mesh constructor
      build_mesh(time_stepper_pt);
    }
 
    RectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      const bool& periodic_in_x,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : Nx(nx),
        Ny(ny),
        Xmin(0.0),
        Xmax(lx),
        Ymin(0.0),
        Ymax(ly),
        Xperiodic(periodic_in_x)
    {
      // Mesh can only be built with 2D Qelements.
      MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
 
      // Call the generic mesh constructor
      build_mesh(time_stepper_pt);
    }
 
    RectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& xmin,
      const double& xmax,
      const double& ymin,
      const double& ymax,
      const bool& periodic_in_x,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : Nx(nx),
        Ny(ny),
        Xmin(xmin),
        Xmax(xmax),
        Ymin(ymin),
        Ymax(ymax),
        Xperiodic(periodic_in_x)
    {
      // Mesh can only be built with 2D Qelements.
      MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
 
      // Call the generic mesh constructor
      build_mesh(time_stepper_pt);
    }
 
    const unsigned& nx() const
    {
      // Return the value of Nx
      return Nx;
    }
 
    const unsigned& ny() const
    {
      // Return the value of Ny
      return Ny;
    }
 
    const double x_min() const
    {
      // Return the value of Xmin
      return Xmin;
    }
 
    const double x_max() const
    {
      // Return the value of Xmax
      return Xmax;
    }
 
    const double y_min() const
    {
      // Return the value of Ymin
      return Ymin;
    }
 
    const double y_max() const
    {
      // Return the value of Ymax
      return Ymax;
    }
 
    virtual void element_reorder();
 
    virtual double x_spacing_function(unsigned xelement,
                                      unsigned xnode,
                                      unsigned yelement,
                                      unsigned ynode)
    {
      // Calculate the values of equal increments in nodal values
      double xstep = (Xmax - Xmin) / ((Np - 1) * Nx);
      // Return the appropriate value
      return (Xmin + xstep * ((Np - 1) * xelement + xnode));
    }
 
    virtual double y_spacing_function(unsigned xelement,
                                      unsigned xnode,
                                      unsigned yelement,
                                      unsigned ynode)
    {
      double ystep = (Ymax - Ymin) / ((Np - 1) * Ny);
      // Return the appropriate value
      return (Ymin + ystep * ((Np - 1) * yelement + ynode));
    }
  };
 
 
  //==========================================================================
  //===========================================================================
  template<class ELEMENT>
  class RefineableRectangularQuadMesh
    : public virtual RectangularQuadMesh<ELEMENT>,
      public RefineableQuadMesh<ELEMENT>
  {
  public:
    RefineableRectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(nx, ny, lx, ly, time_stepper_pt)
    {
      // Nodal positions etc. were created in constructor for
      // RectangularMesh<...>. Only need to setup quadtree forest
      this->setup_quadtree_forest();
    }
 
    RefineableRectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      const bool& periodic_in_x,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(
          nx, ny, lx, ly, periodic_in_x, time_stepper_pt)
    {
      // Nodal positions etc. were created in constructor for
      // RectangularMesh<...>. Only need to setup quadtree forest
      this->setup_quadtree_forest();
    }
 
    RefineableRectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& xmin,
      const double& xmax,
      const double& ymin,
      const double& ymax,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(
          nx, ny, xmin, xmax, ymin, ymax, time_stepper_pt)
    {
      // Nodal positions etc. were created in constructor for
      // RectangularMesh<...>. Only need to setup quadtree forest
      this->setup_quadtree_forest();
    }
 
    RefineableRectangularQuadMesh(
      const unsigned& nx,
      const unsigned& ny,
      const double& xmin,
      const double& xmax,
      const double& ymin,
      const double& ymax,
      const bool& periodic_in_x,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(
          nx, ny, xmin, xmax, ymin, ymax, periodic_in_x, time_stepper_pt)
    {
      // Nodal positions etc. were created in constructor for
      // RectangularMesh<...>. Only need to setup quadtree forest
      this->setup_quadtree_forest();
    }
  };
 
 
 
 
  //================================================================
  //================================================================
  template<class ELEMENT>
  class ElasticRectangularQuadMesh
    : public virtual RectangularQuadMesh<ELEMENT>,
      public virtual SolidMesh
  {
  public:
    ElasticRectangularQuadMesh<ELEMENT>(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      const Vector<double>& origin,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(nx, ny, lx, ly, time_stepper_pt)
    {
      // Translate the nodes
      unsigned nnod = nnode();
      for (unsigned j = 0; j < nnod; j++)
      {
        node_pt(j)->x(0) += origin[0];
        node_pt(j)->x(1) += origin[1];
      }
 
      set_lagrangian_nodal_coordinates();
 
      // Setup boundary coordinates
      set_boundary_coordinates(origin);
    }
 
 
    ElasticRectangularQuadMesh<ELEMENT>(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(nx, ny, lx, ly, time_stepper_pt)
    {
      // No shift
      Vector<double> origin(2, 0.0);
 
      set_lagrangian_nodal_coordinates();
 
      // Setup boundary coordinates
      set_boundary_coordinates(origin);
    }
 
 
    ElasticRectangularQuadMesh<ELEMENT>(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      const bool& periodic_in_x,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(
          nx, ny, lx, ly, periodic_in_x, time_stepper_pt)
    {
      // No shift
      Vector<double> origin(2, 0.0);
 
      set_lagrangian_nodal_coordinates();
 
      // Setup boundary coordinates
      set_boundary_coordinates(origin);
    }
 
  private:
    void set_boundary_coordinates(const Vector<double>& origin)
    {
      // 1D vector fo boundary coordinate
      Vector<double> zeta(1);
 
      // Loop over boundaries 0 and 2 where xi_0 is the varying
      // boundary coordinate
      for (unsigned b = 0; b < 3; b += 2)
      {
        // Number of nodes on those boundaries
        unsigned n_nod = nboundary_node(b);
 
        // Loop over the nodes
        for (unsigned i = 0; i < n_nod; i++)
        {
          // Boundary coordinate varies between 0 and L
          zeta[0] = boundary_node_pt(b, i)->xi(0) - origin[0];
          boundary_node_pt(b, i)->set_coordinates_on_boundary(b, zeta);
        }
        Boundary_coordinate_exists[b] = true;
      }
 
 
      // Loop over boundaries 1 and 3 where xi_1 is the varying
      // boundary coordinate
      for (unsigned b = 1; b < 4; b += 2)
      {
        // Number of nodes on those boundaries
        unsigned n_nod = nboundary_node(b);
 
        // Loop over the nodes
        for (unsigned i = 0; i < n_nod; i++)
        {
          // Boundary coordinate varies between +/- H/2
          zeta[0] = boundary_node_pt(b, i)->xi(1) - origin[1] -
                    0.5 * (this->Ymax - this->Ymin);
          boundary_node_pt(b, i)->set_coordinates_on_boundary(b, zeta);
        }
        Boundary_coordinate_exists[b] = true;
      }
    }
  };
 
 
 
 
  //================================================================
  //================================================================
  template<class ELEMENT>
  class ElasticRefineableRectangularQuadMesh
    : public virtual ElasticRectangularQuadMesh<ELEMENT>,
      public RefineableQuadMesh<ELEMENT>
  {
  public:
    ElasticRefineableRectangularQuadMesh<ELEMENT>(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(nx, ny, lx, ly, time_stepper_pt),
        ElasticRectangularQuadMesh<ELEMENT>(nx, ny, lx, ly, time_stepper_pt)
 
    {
      // Nodal positions etc. were created in base class.
      // Only need to setup quadtree forest
      this->setup_quadtree_forest();
    }
 
 
    ElasticRefineableRectangularQuadMesh<ELEMENT>(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      const bool& periodic_in_x,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(
          nx, ny, lx, ly, periodic_in_x, time_stepper_pt),
        ElasticRectangularQuadMesh<ELEMENT>(
          nx, ny, lx, ly, periodic_in_x, time_stepper_pt)
    {
      // Nodal positions etc. were created in base class.
      // Only need to setup quadtree forest
      this->setup_quadtree_forest();
    }
 
 
    ElasticRefineableRectangularQuadMesh<ELEMENT>(
      const unsigned& nx,
      const unsigned& ny,
      const double& lx,
      const double& ly,
      const Vector<double>& origin,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : RectangularQuadMesh<ELEMENT>(nx, ny, lx, ly, time_stepper_pt),
        ElasticRectangularQuadMesh<ELEMENT>(
          nx, ny, lx, ly, origin, time_stepper_pt)
 
    {
      // Nodal positions etc. were created in base class.
      // Only need to setup quadtree forest
      this->setup_quadtree_forest();
    }
  };
 
 
} // namespace oomph
 
#endif