linearised_axisymmetric_fluid_interface_elements.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 (one-dimensional) free surface elements
//Include guards, to prevent multiple includes
#ifndef OOMPH_PATRICKLINEARISED_INTERFACE_ELEMENTS_HEADER
#define OOMPH_PATRICKLINEARISED_INTERFACE_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif
 
// oomph-lib includes
#include "generic.h"
#include "perturbed_spines.h"
 
namespace oomph
{
 
 //========================================================================
 //========================================================================
 class LinearisedAxisymmetricFluidInterfaceElement
  : public virtual FaceElement
 {
   private:
  
  double *Ca_pt;
  
  double *St_pt;
  
  int *Azimuthal_Mode_Number_pt;
  
  static double Default_Physical_Constant_Value;
  
  static int Default_Azimuthal_Mode_Number_Value;
  
  
   protected:
  
  Vector<unsigned> U_index_interface;
  
  Vector<unsigned> Xhat_index_interface;
  
  virtual int kinematic_local_eqn(const unsigned &n,const unsigned &i)=0;
  
   public:
 
  virtual void hijack_kinematic_conditions(const Vector<unsigned> 
                                           &bulk_node_number)=0;
 
 
   protected:
  
  virtual void fill_in_generic_residual_contribution_interface(
   Vector<double> &residuals, 
   DenseMatrix<double> &jacobian, 
   unsigned flag);
  
  virtual void add_additional_residual_contributions(
   Vector<double> &residuals, DenseMatrix<double> &jacobian,
   const unsigned &flag,
   const Shape &psif, const DShape &dpsifds,
   const Vector<double> &interpolated_n, 
   const double &r, const double &W,
   const double &J) {}
  
  double dXhat_dt(const unsigned &n, const unsigned &i) const
  {
   // Get the data's positional timestepper
   TimeStepper* position_time_stepper_pt
    = this->node_pt(n)->position_time_stepper_pt();
   
   // Get the value
   const double xhat_value = this->node_pt(n)->value(Xhat_index_interface[i]);
   
   // Initialise dXhat/dt
   double dXhatdt = 0.0;
 
   // Loop over the timesteps, if there is a non Steady timestepper
   if(!position_time_stepper_pt->is_steady())
    {
     // Number of timsteps (past & present)
     const unsigned n_time = position_time_stepper_pt->ntstorage();
     
     // Add the contributions to the time derivative
     for(unsigned t=0;t<n_time;t++)
      {
       dXhatdt += position_time_stepper_pt->weight(1,t)*xhat_value;
      }
    }
   
   return dXhatdt;
  }
  
 
   public:
  
  LinearisedAxisymmetricFluidInterfaceElement()//: Pext_data_pt(0) 
   {
    // Set all the physical parameter pointers to the default value of zero
    Ca_pt = &Default_Physical_Constant_Value;
    St_pt = &Default_Physical_Constant_Value;
    
    // Set the azimuthal mode number to default (zero)
    Azimuthal_Mode_Number_pt = &Default_Azimuthal_Mode_Number_Value;
   }
  
  virtual double sigma(const Vector<double> &s_local) { return 1.0; }
  
  void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   // Add the residual contributions
   fill_in_generic_residual_contribution_interface(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }
  
  const double &ca() const { return *Ca_pt; }
  
  double* &ca_pt() { return Ca_pt; }
  
  const double &st() const { return *St_pt; }
  
  double* &st_pt() { return St_pt; }
  
  const int &azimuthal_mode_number() const
   { return *Azimuthal_Mode_Number_pt; }
  
  int* &azimuthal_mode_number_pt() { return Azimuthal_Mode_Number_pt; }
  
  double u(const unsigned &n, const unsigned &i)
  {
   return node_pt(n)->value(U_index_interface[i]);
  }
 
  double interpolated_u(const Vector<double> &s, const unsigned &i)
  {
   // Find number of nodes
   const unsigned n_node = nnode();
   
   // Storage for the local shape function
   Shape psi(n_node);
   
   // Get values of shape function at local coordinate s
   this->shape(s,psi);
   
   // Initialise value of u
   double interpolated_u = 0.0;
   
   // Loop over the local nodes and sum
   for(unsigned l=0;l<n_node;l++) { interpolated_u += u(l,i)*psi(l); }
   
   return(interpolated_u);
  }
  
  void output(std::ostream &outfile) { FiniteElement::output(outfile); }
  
  void output(std::ostream &outfile, const unsigned &n_plot);
  
  void output(FILE* file_pt) { FiniteElement::output(file_pt); }
  
  void output(FILE* file_pt, const unsigned &n_plot);
 
};
 
 
 //========================================================================
 //========================================================================
 template<class ELEMENT>
  class PerturbedSpineLinearisedAxisymmetricFluidInterfaceElement : 
 public virtual Hijacked<PerturbedSpineElement<FaceGeometry<ELEMENT> > >,
  public virtual LinearisedAxisymmetricFluidInterfaceElement 
  {
    private:
   
   int kinematic_local_eqn(const unsigned &n, const unsigned &i)
   {
    return this->spine_local_eqn(n,i);
   }
   
   
   void hijack_kinematic_conditions(const Vector<unsigned> &bulk_node_number)
   {
    // Loop over all the passed nodes
    for(Vector<unsigned>::const_iterator it=bulk_node_number.begin();
        it!=bulk_node_number.end();++it)
     {
      // Hijack the spine heights and delete the returned data objects
      delete this->hijack_nodal_spine_value(*it,0);
     }
   }
   
   double dH_dt(const unsigned &n, const unsigned &i) const
   {
    // Get the data's timestepper
    TimeStepper* time_stepper_pt = this->node_pt(n)->time_stepper_pt();
    
    // Upcast from general node to PerturbedSpineNode
    PerturbedSpineNode* perturbed_spine_node_pt =
     dynamic_cast<PerturbedSpineNode*>(this->node_pt(n));
    
    // Initialise dH/dt
    double dHdt = 0.0;
   
    // Loop over the timesteps, if there is a non Steady timestepper
    if(!time_stepper_pt->is_steady())
     {
      // Number of timsteps (past & present)
      const unsigned n_time = time_stepper_pt->ntstorage();
      
      // Add the contributions to the time derivative
      for(unsigned t=0;t<n_time;t++)
       {
        dHdt += time_stepper_pt->weight(1,t)*
         perturbed_spine_node_pt->perturbed_spine_pt()->height(t,i);
       }
     }
    
    return dHdt;
   }
   
    public:
 
   PerturbedSpineLinearisedAxisymmetricFluidInterfaceElement(
    FiniteElement* const &bulk_element_pt, const int &face_index) : 
    Hijacked<PerturbedSpineElement<FaceGeometry<ELEMENT> > >(), 
    LinearisedAxisymmetricFluidInterfaceElement()
     {
      // Attach the geometrical information to the element
      bulk_element_pt->build_face_element(face_index,this);
      
      // Find the index at which the velocity unknowns are stored 
      // from the bulk element
      ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(bulk_element_pt);
      
      // We must have six velocity components
      this->U_index_interface.resize(6);
      for(unsigned i=0;i<6;i++)
       {
        this->U_index_interface[i]=cast_element_pt->u_index_lin_axi_nst(i);
       }
      
      // We must have four components of the perturbations to the
      // nodal positions
      this->Xhat_index_interface.resize(4);
      for(unsigned i=0;i<4;i++)
       {
        this->Xhat_index_interface[i] = 
         cast_element_pt->xhat_index_lin_axi_nst(i);
       }
     }
    
    void fill_in_contribution_to_jacobian(Vector<double> &residuals, 
                                          DenseMatrix<double> &jacobian)
    {
     // Call the generic residuals routine with the flag set to 1
     fill_in_generic_residual_contribution_interface(residuals,jacobian,1);
     // Call the generic routine to handle the spine variables
     // PerturbedSpineElement<FaceGeometry<ELEMENT> >::
     this->fill_in_jacobian_from_geometric_data(jacobian);
    }
    
    void fill_in_generic_residual_contribution_interface(
     Vector<double> &residuals, 
     DenseMatrix<double> &jacobian, 
     unsigned flag);
    
    void output(std::ostream &outfile) {FiniteElement::output(outfile);}
    
    void output(std::ostream &outfile, const unsigned &n_plot)
    {LinearisedAxisymmetricFluidInterfaceElement::output(outfile,n_plot);}
    
    void output(FILE* file_pt) {FiniteElement::output(file_pt);}
    
    void output(FILE* file_pt, const unsigned &n_plot)
    {LinearisedAxisymmetricFluidInterfaceElement::output(file_pt,n_plot);}
 
    void output_interface_position(std::ostream &outfile,
                                   const unsigned &nplot);
    
    void output_perturbation_to_interface(std::ostream &outfile,
                                          const unsigned &nplot);
    
    double interpolated_H(const Vector<double> &s,
                          const unsigned &i) const
    {
     // Determine number of nodes in the element
     const unsigned n_node = nnode();
     
     // Provide storage for local shape functions
     Shape psi(n_node);
     
     // Find values of shape functions
     shape(s,psi);
     
     // Initialise value of H
     double interpolated_H = 0.0;
     
     // Loop over the shape functions and sum
     for(unsigned l=0;l<n_node;l++) 
      {
       // Upcast from general node to PerturbedSpineNode
       PerturbedSpineNode* perturbed_spine_node_pt =
        dynamic_cast<PerturbedSpineNode*>(this->node_pt(l));
       
       // Calculate interpolated i-th component of perturbed spine "heights"
       interpolated_H +=
        perturbed_spine_node_pt->perturbed_spine_pt()->height(i)*psi[l];
      }
     
     return(interpolated_H);
    }
    
  };
 
 
 
 
}
 
#endif