channel_spine_mesh.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//====================================================================
#ifndef OOMPH_CHANNEL_SPINE_MESH_HEADER
#define OOMPH_CHANNEL_SPINE_MESH_HEADER
 
// oomph-lib includes
#include "../generic/spines.h"
#include "rectangular_quadmesh.template.h"
 
namespace oomph
{
  //======================================================================
  //======================================================================
  template<class ELEMENT>
  class ChannelSpineMesh : public RectangularQuadMesh<ELEMENT>, public SpineMesh
  {
  public:
    ChannelSpineMesh(const unsigned& nx0,
                     const unsigned& nx1,
                     const unsigned& nx2,
                     const unsigned& ny,
                     const double& lx0,
                     const double& lx1,
                     const double& lx2,
                     const double& h,
                     GeomObject* wall_pt,
                     TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper);
 
    ChannelSpineMesh(const unsigned& nx0,
                     const unsigned& nx1,
                     const unsigned& nx2,
                     const unsigned& ny,
                     const double& lx0,
                     const double& lx1,
                     const double& lx2,
                     const double& h,
                     GeomObject* wall_pt,
                     const bool& periodic_in_x,
                     TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper);
 
    FiniteElement*& left_element_pt(const unsigned long& i)
    {
      return Left_element_pt[i];
    }
 
    FiniteElement*& centre_element_pt(const unsigned long& i)
    {
      return Centre_element_pt[i];
    }
 
    FiniteElement*& right_element_pt(const unsigned long& i)
    {
      return Right_element_pt[i];
    }
 
    unsigned long nleft() const
    {
      return Left_element_pt.size();
    }
 
    unsigned long ncentre() const
    {
      return Centre_element_pt.size();
    }
 
    unsigned long nright() const
    {
      return Right_element_pt.size();
    }
 
    unsigned long nbulk() const
    {
      unsigned long Nbulk = Left_element_pt.size() + Centre_element_pt.size() +
                            Right_element_pt.size();
      return Nbulk;
    }
 
    void element_reorder();
 
    virtual void spine_node_update(SpineNode* spine_node_pt)
    {
      // Get spine node's fraction along the spine
      double W = spine_node_pt->fraction();
 
      // Get local coordinates
      Vector<double> s_wall(1);
      s_wall[0] = spine_node_pt->spine_pt()->geom_parameter(0);
 
      // Get position vector to wall
      Vector<double> position(2);
      spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_wall, position);
 
      // Set the value of y
      spine_node_pt->x(1) = this->Ymin + W * position[1];
    }
 
    virtual double x_spacing_function(unsigned xelement,
                                      unsigned xnode,
                                      unsigned yelement,
                                      unsigned ynode)
    {
      // Calculate the values of equal increments in nodal values in left region
      double xstep1 = Lx0 / ((this->Np - 1) * Nx0);
      // Calculate the values of equal increments in nodal values in centre
      // region
      double xstep2 = Lx1 / ((this->Np - 1) * Nx1);
      // Calculate the values of equal increments in nodal values in right
      // region
      double xstep3 = Lx2 / ((this->Np - 1) * Nx2);
 
      // left region
      if (xelement < Nx0)
      {
        // Return the appropriate value
        return (this->Xmin + xstep1 * ((this->Np - 1) * xelement + xnode));
      }
      // centre region
      else if (xelement < Nx0 + Nx1)
      {
        // Return the appropriate value
        return (Lx0 + xstep2 * ((this->Np - 1) * (xelement - Nx0) + xnode));
      }
      // right region
      else if (xelement < Nx0 + Nx1 + Nx2)
      {
        // Return the appropriate value
        return (Lx0 + Lx1 +
                xstep3 * ((this->Np - 1) * (xelement - Nx0 - Nx1) + xnode));
      }
      else
      {
        throw OomphLibError("Should not have got here",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      // Dummy return to keep compiler from barking
      return 0.0;
    }
 
    Spine*& left_spine_pt(const unsigned long& i)
    {
#ifdef PARNOID
      if (i > Nleft_spine)
      {
        throw OomphLibError("Arguemnt out of range",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Spine_pt[i];
    }
 
    Spine*& centre_spine_pt(const unsigned long& i)
    {
      if (i > Ncentre_spine)
      {
        throw OomphLibError("Arguemnt out of range",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else
      {
        return Spine_pt[Nleft_spine + i];
      }
    }
 
    Spine*& right_spine_pt(const unsigned long& i)
    {
      if (i > Nright_spine)
      {
        throw OomphLibError("Arguemnt out of range",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else
      {
        return Spine_pt[Nleft_spine + Ncentre_spine - 1 + i];
      }
    }
 
    unsigned nleft_spine()
    {
      return Nleft_spine;
    }
 
    unsigned ncentre_spine()
    {
      return Ncentre_spine;
    }
 
    unsigned nright_spine()
    {
      return Nright_spine;
    }
 
    GeomObject* wall_pt()
    {
      return Wall_pt;
    }
 
    GeomObject* straight_wall_pt()
    {
      return Straight_wall_pt;
    }
 
  protected:
    Vector<FiniteElement*> Left_element_pt;
 
    Vector<FiniteElement*> Centre_element_pt;
 
    Vector<FiniteElement*> Right_element_pt;
 
    virtual void build_channel_spine_mesh(TimeStepper* time_stepper_pt);
 
    unsigned Nx0;
 
    unsigned Nx1;
 
    unsigned Nx2;
 
    double Lx0;
 
    double Lx1;
 
    double Lx2;
 
    unsigned Nleft_spine;
 
    unsigned Ncentre_spine;
 
    unsigned Nright_spine;
 
    GeomObject* Wall_pt;
 
    GeomObject* Straight_wall_pt;
  };
 
} // namespace oomph
 
#endif