bretherton_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_BRETHERTON_SPINE_MESH_HEADER
#define OOMPH_BRETHERTON_SPINE_MESH_HEADER
 
 
// The mesh
#include "single_layer_spine_mesh.template.h"
#include "simple_rectangular_quadmesh.template.h"
 
namespace oomph
{
  //============================================================================
  //============================================================================
  template<class ELEMENT, class INTERFACE_ELEMENT>
  class BrethertonSpineMesh : public SingleLayerSpineMesh<ELEMENT>
  {
  public:
    BrethertonSpineMesh(
      const unsigned& nx1,
      const unsigned& nx2,
      const unsigned& nx3,
      const unsigned& nh,
      const unsigned& nhalf,
      const double& h,
      GeomObject* lower_wall_pt,
      GeomObject* upper_wall_pt,
      const double& zeta_start,
      const double& zeta_transition_start,
      const double& zeta_transition_end,
      const double& zeta_end,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper);
 
 
    FiniteElement*& interface_element_pt(const unsigned long& i)
    {
      return Interface_element_pt[i];
    }
 
    unsigned long ninterface_element() const
    {
      return Interface_element_pt.size();
    }
 
    FiniteElement*& bulk_element_pt(const unsigned long& i)
    {
      return Bulk_element_pt[i];
    }
 
    unsigned long nbulk() const
    {
      return Bulk_element_pt.size();
    }
 
 
    unsigned nfree_surface_spines()
    {
      unsigned np = this->finite_element_pt(0)->nnode_1d();
      return 2 * (Nx1 + Nx2 + Nhalf) * (np - 1);
    }
 
 
    double find_distance_to_free_surface(GeomObject* const& fs_geom_object_pt,
                                         Vector<double>& initial_zeta,
                                         const Vector<double>& spine_base,
                                         const Vector<double>& spine_end);
 
    void reposition_spines(const double& zeta_lo_transition_start,
                           const double& zeta_lo_transition_end,
                           const double& zeta_up_transition_start,
                           const double& zeta_up_transition_end);
 
    void pin_all_spines()
    {
      unsigned n_spine = this->nspine();
      for (unsigned i = 0; i < n_spine; i++)
      {
        this->spine_pt(i)->spine_height_pt()->pin(0);
      }
    }
 
    void spine_node_update(SpineNode* spine_node_pt)
    {
      unsigned id = spine_node_pt->node_update_fct_id();
      switch (id)
      {
        case 0:
          spine_node_update_film_lower(spine_node_pt);
          break;
 
        case 1:
          spine_node_update_horizontal_transition_lower(spine_node_pt);
          break;
 
        case 2:
          spine_node_update_vertical_transition_lower(spine_node_pt);
          break;
 
        case 3:
          spine_node_update_vertical_transition_upper(spine_node_pt);
          break;
 
        case 4:
          spine_node_update_horizontal_transition_upper(spine_node_pt);
          break;
 
        case 5:
          spine_node_update_film_upper(spine_node_pt);
          break;
 
        case 6:
          spine_node_update_channel(spine_node_pt);
          break;
 
        default:
          std::ostringstream error_message;
          error_message << "Incorrect spine update id " << id << std::endl;
 
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
    }
 
 
    ELEMENT* control_element_pt()
    {
      return Control_element_pt;
    }
 
 
    double spine_centre_fraction() const
    {
      return *Spine_centre_fraction_pt;
    }
 
    void set_spine_centre_fraction_pt(double* const& fraction_pt)
    {
      Spine_centre_fraction_pt = fraction_pt;
    }
 
  protected:
    void spine_node_update_film_lower(SpineNode* spine_node_pt)
    {
      // Get fraction along the spine
      double w = spine_node_pt->fraction();
      // Get spine height
      double h = spine_node_pt->h();
 
      // Get wall coordinate
      Vector<double> s_lo(1);
      s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
 
      // Get position vector to wall
      Vector<double> r_wall_lo(2);
      spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_wall_lo);
 
      // Update the nodal position
      spine_node_pt->x(0) = r_wall_lo[0];
      spine_node_pt->x(1) = r_wall_lo[1] + w * h;
    }
 
 
    void spine_node_update_horizontal_transition_lower(SpineNode* spine_node_pt)
    {
      // Get fraction along the spine
      double w = spine_node_pt->fraction();
      // Get spine height
      double h = spine_node_pt->h();
 
      // Get wall coordinate
      Vector<double> s_lo(1);
      s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
 
      // Get position vector to wall
      Vector<double> r_wall_lo(2);
      spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_wall_lo);
 
      // Get the spine centre
      Vector<double> s_transition_lo(1), s_transition_up(1);
      s_transition_lo[0] = spine_node_pt->spine_pt()->geom_parameter(1);
      s_transition_up[0] = spine_node_pt->spine_pt()->geom_parameter(2);
      Vector<double> r_transition_lo(2), r_transition_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_transition_lo,
                                                             r_transition_lo);
      spine_node_pt->spine_pt()->geom_object_pt(2)->position(s_transition_up,
                                                             r_transition_up);
 
      Vector<double> spine_centre(2);
      // Horizontal position is always halfway between the walls
      spine_centre[0] = 0.5 * (r_transition_lo[0] + r_transition_up[0]);
      // Vertical spine centre is given by fraction between the walls
      spine_centre[1] =
        r_transition_lo[1] +
        spine_centre_fraction() * (r_transition_up[1] - r_transition_lo[1]);
 
      // Get vector twoards spine origin
      Vector<double> N(2);
      N[0] = spine_centre[0] - r_wall_lo[0];
      N[1] = spine_centre[1] - r_wall_lo[1];
      double inv_length = 1.0 / sqrt(N[0] * N[0] + N[1] * N[1]);
      // Update the nodal position
      spine_node_pt->x(0) = r_wall_lo[0] + w * h * N[0] * inv_length;
      spine_node_pt->x(1) = r_wall_lo[1] + w * h * N[1] * inv_length;
    }
 
 
    void spine_node_update_vertical_transition_lower(SpineNode* spine_node_pt)
    {
      // Get fraction along the spine
      double w = spine_node_pt->fraction();
      // Get spine height
      double h = spine_node_pt->h();
 
      // Get local coordinates
      Vector<double> s_lo(1), s_up(1);
      s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
      s_up[0] = spine_node_pt->spine_pt()->geom_parameter(1);
 
      // Get position vector to wall
      Vector<double> r_lo(2), r_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_lo);
      spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_up, r_up);
 
      // Get fraction along vertical line across
      double vertical_fraction = spine_node_pt->spine_pt()->geom_parameter(2);
 
      // Origin of spine
      Vector<double> S0(2);
      S0[0] = r_lo[0] + vertical_fraction * (r_up[0] - r_lo[0]);
      S0[1] = r_lo[1] + vertical_fraction * (r_up[1] - r_lo[1]);
 
 
      // Get the spine centre
      Vector<double> s_transition_lo(1), s_transition_up(1);
      s_transition_lo[0] = spine_node_pt->spine_pt()->geom_parameter(3);
      s_transition_up[0] = spine_node_pt->spine_pt()->geom_parameter(4);
      Vector<double> r_transition_lo(2), r_transition_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(2)->position(s_transition_lo,
                                                             r_transition_lo);
      spine_node_pt->spine_pt()->geom_object_pt(3)->position(s_transition_up,
                                                             r_transition_up);
 
      Vector<double> spine_centre(2);
      // Horizontal position is always halfway between the walls
      spine_centre[0] = 0.5 * (r_transition_lo[0] + r_transition_up[0]);
      // Vertical spine centre is given by fraction between the walls
      spine_centre[1] =
        r_transition_lo[1] +
        spine_centre_fraction() * (r_transition_up[1] - r_transition_lo[1]);
 
      // Get vector towards origin
      Vector<double> N(2);
      N[0] = spine_centre[0] - S0[0];
      N[1] = spine_centre[1] - S0[1];
 
      double inv_length = 1.0 / sqrt(N[0] * N[0] + N[1] * N[1]);
      // Update the nodal position
      spine_node_pt->x(0) = S0[0] + w * h * N[0] * inv_length;
      spine_node_pt->x(1) = S0[1] + w * h * N[1] * inv_length;
    }
 
 
    void spine_node_update_vertical_transition_upper(SpineNode* spine_node_pt)
    {
      // Get fraction along the spine
      double w = spine_node_pt->fraction();
      // Get spine height
      double h = spine_node_pt->h();
 
      // Get local coordinates
      Vector<double> s_lo(1), s_up(1);
      s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
      s_up[0] = spine_node_pt->spine_pt()->geom_parameter(1);
 
      // Get position vector to wall
      Vector<double> r_lo(2), r_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_lo);
      spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_up, r_up);
 
      // Get fraction along vertical line across
      double vertical_fraction = spine_node_pt->spine_pt()->geom_parameter(2);
 
      // Origin of spine
      Vector<double> S0(2);
      S0[0] = r_lo[0] + vertical_fraction * (r_up[0] - r_lo[0]);
      S0[1] = r_lo[1] + vertical_fraction * (r_up[1] - r_lo[1]);
 
 
      // Get the spine centre
      Vector<double> s_transition_lo(1), s_transition_up(1);
      s_transition_lo[0] = spine_node_pt->spine_pt()->geom_parameter(3);
      s_transition_up[0] = spine_node_pt->spine_pt()->geom_parameter(4);
      Vector<double> r_transition_lo(2), r_transition_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(2)->position(s_transition_lo,
                                                             r_transition_lo);
      spine_node_pt->spine_pt()->geom_object_pt(3)->position(s_transition_up,
                                                             r_transition_up);
 
      Vector<double> spine_centre(2);
      // Horizontal position is always halfway between the walls
      spine_centre[0] = 0.5 * (r_transition_lo[0] + r_transition_up[0]);
      // Vertical spine centre is given by fraction between the walls
      spine_centre[1] =
        r_transition_lo[1] +
        spine_centre_fraction() * (r_transition_up[1] - r_transition_lo[1]);
 
      // Get vector towards origin
      Vector<double> N(2);
      N[0] = spine_centre[0] - S0[0];
      N[1] = spine_centre[1] - S0[1];
 
      double inv_length = 1.0 / sqrt(N[0] * N[0] + N[1] * N[1]);
      // Update the nodal position
      spine_node_pt->x(0) = S0[0] + w * h * N[0] * inv_length;
      spine_node_pt->x(1) = S0[1] + w * h * N[1] * inv_length;
    }
 
 
    void spine_node_update_horizontal_transition_upper(SpineNode* spine_node_pt)
    {
      // Get fraction along the spine
      double w = spine_node_pt->fraction();
      // Get spine height
      double h = spine_node_pt->h();
 
      // Get wall coordinate
      Vector<double> s_up(1);
      s_up[0] = spine_node_pt->spine_pt()->geom_parameter(0);
 
      // Get position vector to wall
      Vector<double> r_wall_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_up, r_wall_up);
 
      // Get the spine centre
      Vector<double> s_transition_lo(1), s_transition_up(1);
      s_transition_lo[0] = spine_node_pt->spine_pt()->geom_parameter(1);
      s_transition_up[0] = spine_node_pt->spine_pt()->geom_parameter(2);
      Vector<double> r_transition_lo(2), r_transition_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_transition_lo,
                                                             r_transition_lo);
      spine_node_pt->spine_pt()->geom_object_pt(2)->position(s_transition_up,
                                                             r_transition_up);
 
      Vector<double> spine_centre(2);
      // Horizontal position is always halfway between the walls
      spine_centre[0] = 0.5 * (r_transition_lo[0] + r_transition_up[0]);
      // Vertical spine centre is given by fraction between the walls
      spine_centre[1] =
        r_transition_lo[1] +
        spine_centre_fraction() * (r_transition_up[1] - r_transition_lo[1]);
 
      // Get vector towards origin
      Vector<double> N(2);
      N[0] = spine_centre[0] - r_wall_up[0];
      N[1] = spine_centre[1] - r_wall_up[1];
      double inv_length = 1.0 / sqrt(N[0] * N[0] + N[1] * N[1]);
 
      // Update the nodal position
      spine_node_pt->x(0) = r_wall_up[0] + w * h * N[0] * inv_length;
      spine_node_pt->x(1) = r_wall_up[1] + w * h * N[1] * inv_length;
    }
 
 
    void spine_node_update_film_upper(SpineNode* spine_node_pt)
    {
      // Get fraction along the spine
      double w = spine_node_pt->fraction();
      // Get spine height
      double h = spine_node_pt->h();
 
      // Get wall coordinate
      Vector<double> s_up(1);
      s_up[0] = spine_node_pt->spine_pt()->geom_parameter(0);
 
      // Get position vector to wall
      Vector<double> r_wall_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_up, r_wall_up);
 
      // Update the nodal position
      spine_node_pt->x(0) = r_wall_up[0];
      spine_node_pt->x(1) = r_wall_up[1] - w * h;
    }
 
 
    void spine_node_update_channel(SpineNode* spine_node_pt)
    {
      // Get fraction along the spine
      double w = spine_node_pt->fraction();
 
      // Get upper and lower local coordinates
      Vector<double> s_lo(1), s_up(1);
      s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
      s_up[0] = spine_node_pt->spine_pt()->geom_parameter(1);
 
      // Get position vector to lower wall
      Vector<double> r_lo(2), r_up(2);
      spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_lo);
      spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_up, r_up);
 
      // Update the nodal position
      spine_node_pt->x(0) = r_lo[0] + w * (r_up[0] - r_lo[0]);
      spine_node_pt->x(1) = r_lo[1] + w * (r_up[1] - r_lo[1]);
    }
 
 
    void initial_element_reorder();
 
    unsigned Nx1;
 
    unsigned Nx2;
 
    unsigned Nx3;
 
    unsigned Nhalf;
 
    unsigned Nh;
 
    double H;
 
    GeomObject* Upper_wall_pt;
 
    GeomObject* Lower_wall_pt;
 
    double Zeta_start;
 
    double Zeta_end;
 
    double Zeta_transition_start;
 
    double Zeta_transition_end;
 
    double* Spine_centre_fraction_pt;
 
    double Default_spine_centre_fraction;
 
    ELEMENT* Control_element_pt;
 
    Vector<FiniteElement*> Bulk_element_pt;
 
    Vector<FiniteElement*> Interface_element_pt;
  };
 
} // namespace oomph
 
#endif