refineable_linearised_navier_stokes_elements.cc

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//    Version 1.0; svn revision $LastChangedRevision: 1097 $
//LIC//
//LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
//LIC// 
//LIC// Copyright (C) 2006-2016 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//====================================================================
// Non-inline functions for the refineable linearised axisymmetric
// Navier-Stokes elements

// oomph-lib includes
#include "refineable_linearised_navier_stokes_elements.h"

namespace oomph
{

 //======================================================================= 
 /// Compute the residuals for the refineable linearised axisymmetric
 /// Navier--Stokes equations; flag=1(or 0): do (or don't) compute the
 /// Jacobian as well.
 //======================================================================= 
 void RefineableLinearisedNavierStokesEquations::
 fill_in_generic_residual_contribution_linearised_nst(
  Vector<double> &residuals, 
  DenseMatrix<double> &jacobian, 
  DenseMatrix<double> &mass_matrix,
  unsigned flag)
 {
   // Get the time from the first node in the element
  const double time = this->node_pt(0)->time_stepper_pt()->time();

  // Determine number of nodes in the element
  const unsigned n_node = nnode();
 
  // Determine how many pressure values there are associated with
  // a single pressure component
  const unsigned n_pres = npres_linearised_nst();

  const unsigned n_veloc = 4*DIM;
  
  // Get the nodal indices at which the velocity is stored
  unsigned u_nodal_index[n_veloc];
  for(unsigned i=0;i<n_veloc;++i)
   { 
    u_nodal_index[i] = u_index_linearised_nst(i);
   }
  
  // Which nodal values represent the two pressure components?
  // (Negative if pressure is not based on nodal interpolation).
  Vector<int> p_index(2);
  for(unsigned i=0;i<2;i++)
   {
    p_index[i] = this->p_index_linearised_nst(i);
   }

  // Local array of booleans that are true if the l-th pressure value is
  // hanging (avoid repeated virtual function calls)
  bool pressure_dof_is_hanging[n_pres];

  // If the pressure is stored at a node
  if(p_index[0] >= 0)
   {
    // Read out whether the pressure is hanging
    for(unsigned l=0;l<n_pres;++l)
     {
      pressure_dof_is_hanging[l] = 
       pressure_node_pt(l)->is_hanging(p_index[0]);
     }
   }
  // Otherwise the pressure is not stored at a node and so cannot hang
  else
   {
    for(unsigned l=0;l<n_pres;++l)
     {
      pressure_dof_is_hanging[l] = false;
     }
   }


  // Set up memory for the fluid shape and test functions
  Shape psif(n_node), testf(n_node);
  DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
  
  // Set up memory for the pressure shape and test functions
  Shape psip(n_pres), testp(n_pres);
  
  // Determine number of integration points
  const unsigned n_intpt = integral_pt()->nweight();
  
  // Set up memory for the vector to hold local coordinates (two dimensions)
  Vector<double> s(DIM);

  // Get physical variables from the element
  // (Reynolds number must be multiplied by the density ratio)
  const double scaled_re = re()*density_ratio();
  const double scaled_re_st = re_st()*density_ratio();
  const double visc_ratio = viscosity_ratio();

  const double eval_real = lambda();
  const double eval_imag = omega();

  const std::complex<double> eigenvalue(eval_real,eval_imag);
  
  // Integers used to store the local equation numbers
  int local_eqn = 0;

  // Local storage for pointers to hang info objects
  HangInfo *hang_info_pt=0;
  
  // Loop over the integration points
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {
    // Assign values of the local coordinates s
    for(unsigned i=0;i<DIM;i++) { s[i] = integral_pt()->knot(ipt,i); }

    // Get the integral weight
    const double w = integral_pt()->weight(ipt);
    
    // Calculate the fluid shape and test functions, and their derivatives
    // w.r.t. the global coordinates
    const double J = 
     dshape_and_dtest_eulerian_at_knot_linearised_nst(ipt,psif,dpsifdx,
                                                          testf,dtestfdx);
    
    // Calculate the pressure shape and test functions
    pshape_linearised_nst(s,psip,testp);
    
    // Premultiply the weights and the Jacobian of the mapping between
    // local and global coordinates
    const double W = w*J;
    
    // Allocate storage for the position and the derivative of the 
    // mesh positions w.r.t. time
    Vector<double> interpolated_x(DIM,0.0);
    //Vector<double> mesh_velocity(2,0.0);

    // Allocate storage for the velocity components (six of these)
    // and their derivatives w.r.t. time
    Vector<std::complex< double> > interpolated_u(DIM);
    //Vector<double> dudt(6,0.0);
    //Allocate storage for the eigen function normalisation
    Vector<std::complex<double> >interpolated_u_normalisation(DIM);
    for(unsigned i=0;i<DIM;++i)
     {
      interpolated_u[i].real(0.0); interpolated_u[i].imag(0.0);
      interpolated_u_normalisation[i].real(0.0);
      interpolated_u_normalisation[i].imag(0.0);
     }
    
    // Allocate storage for the pressure components (two of these
    std::complex<double>  interpolated_p(0.0,0.0);
    std::complex<double>  interpolated_p_normalisation(0.0,0.0);
    
    // Allocate storage for the derivatives of the velocity components
    // w.r.t. global coordinates (r and z)    
    Vector<Vector<std::complex<double> > >interpolated_dudx(DIM);
    for(unsigned i=0;i<DIM;++i)
     {
      interpolated_dudx[i].resize(DIM);
      for(unsigned j=0;j<DIM;++j)
       {
        interpolated_dudx[i][j].real(0.0);
        interpolated_dudx[i][j].imag(0.0);
       }
     }
    
    // Calculate pressure at the integration point
    // -------------------------------------------

    // Loop over pressure degrees of freedom (associated with a single
    // pressure component) in the element
    for(unsigned l=0;l<n_pres;l++)
     {
      // Cache the shape function
      const double psip_ = psip(l);

      //Get the complex nodal pressure values
      const std::complex<double>
       p_value(this->p_linearised_nst(l,0),this->p_linearised_nst(l,1));

      // Add contribution
      interpolated_p += p_value*psip_;

      //Repeat for normalisation
      const std::complex<double>
       p_norm_value(this->p_linearised_nst(l,2),this->p_linearised_nst(l,3));
      interpolated_p_normalisation += p_norm_value*psip_;
     }
    // End of loop over the pressure degrees of freedom in the element
   
    // Calculate velocities and their derivatives at the integration point
    // -------------------------------------------------------------------

    // Loop over the element's nodes
    for(unsigned l=0;l<n_node;l++) 
     {
      // Cache the shape function
      const double psif_ = psif(l);

      // Loop over the DIM coordinate directions
      for(unsigned i=0;i<DIM;i++)
       {
        interpolated_x[i] += this->nodal_position(l,i)*psif_;
       }
      
     // Loop over the DIM complex velocity components
     for(unsigned i=0;i<DIM;i++)
      {
       // Get the value
       const std::complex<double>
        u_value(this->nodal_value(l,u_nodal_index[2*i+0]),
                this->nodal_value(l,u_nodal_index[2*i+1]));

       // Add contribution
       interpolated_u[i] += u_value*psif_;

       // Add contribution to dudt
       //dudt[i] += du_dt_linearised_nst(l,i)*psif_;

       // Loop over two coordinate directions (for derivatives)
       for(unsigned j=0;j<DIM;j++)
        {
         interpolated_dudx[i][j] += u_value*dpsifdx(l,j);
        }
       
       //Interpolate the normalisation function
       const std::complex<double>
        normalisation_value(this->nodal_value(l,u_nodal_index[2*(DIM+i)]),
                            this->nodal_value(l,
                                              u_nodal_index[2*(DIM+i)+1]));
       interpolated_u_normalisation[i] += normalisation_value*psif_;
      }
     } // End of loop over the element's nodes

    // Get the mesh velocity if ALE is enabled
    /*if(!ALE_is_disabled)
     {
      // Loop over the element's nodes
      for(unsigned l=0;l<n_node;l++) 
       {
        // Loop over the two coordinate directions
        for(unsigned i=0;i<2;i++)
         {
          mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif(l);
         }
       }
       }*/

    // Get velocities and their derivatives from base flow problem
    // -----------------------------------------------------------

    // Allocate storage for the velocity components of the base state
    // solution (initialise to zero)
    Vector<double> base_flow_u(DIM,0.0);

    // Get the user-defined base state solution velocity components
    get_base_flow_u(time,ipt,interpolated_x,base_flow_u);

    // Allocate storage for the derivatives of the base state solution's
    // velocity components w.r.t. global coordinate (r and z)
    // N.B. the derivatives of the base flow components w.r.t. the
    // azimuthal coordinate direction (theta) are always zero since the
    // base flow is axisymmetric
    DenseMatrix<double> base_flow_dudx(DIM,DIM,0.0);

    // Get the derivatives of the user-defined base state solution
    // velocity components w.r.t. global coordinates
    get_base_flow_dudx(time,ipt,interpolated_x,base_flow_dudx);


   //MOMENTUM EQUATIONS
   //------------------

    // Number of master nodes
    unsigned n_master = 1;
    
    // Storage for the weight of the shape function
    double hang_weight = 1.0;
    
   //Loop over the test functions
    for(unsigned l=0;l<n_node;l++)
     {
      //Local boolean to indicate whether the node is hanging
      bool is_node_hanging = node_pt(l)->is_hanging();

      if(is_node_hanging)
       {
        hang_info_pt = node_pt(l)->hanging_pt();
        
        //Read out number of master nodes from hanging data
        n_master = hang_info_pt->nmaster();
       }
      //Otherwise the node is its own master
      else
       {
        n_master = 1;
       }
      
      //Loop over the master nodes
      for(unsigned m=0;m<n_master;m++)
       {
        //Loop over the velocity components
        for(unsigned i=0;i<DIM;i++)
         {
          //Assemble the residuals
          //Time dependent term
          std::complex<double> residual_contribution
           = -scaled_re_st*eigenvalue*interpolated_u[i]*testf[l]*W;
          //Pressure term
          residual_contribution += interpolated_p*dtestfdx(l,i)*W;
          //Viscous terms
          for(unsigned k=0;k<DIM;++k)
         {
          residual_contribution -= visc_ratio*
           (interpolated_dudx[i][k] + Gamma[i]*interpolated_dudx[k][i])*
           dtestfdx(l,k)*W;
         }
        
        //Advective terms
        for(unsigned k=0;k<DIM;++k)
         {
          residual_contribution -= scaled_re*(
           base_flow_u[k]*interpolated_dudx[i][k]
           + interpolated_u[k]*base_flow_dudx(i,k))*testf[l]*W;
         }
        
        //Now separate real and imaginary parts

        if(is_node_hanging)
         {
          local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
                                           u_nodal_index[2*i]);
          hang_weight = hang_info_pt->master_weight(m);
         }
        //If node is not hanging number or not
        else
         {
          local_eqn = nodal_local_eqn(l,u_nodal_index[2*i]);
          hang_weight = 1.0;
         }
          
        if(local_eqn >= 0)
         {
          residuals[local_eqn] += residual_contribution.real()*hang_weight;
         }


        if(is_node_hanging)
         {
          local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
                                           u_nodal_index[2*i+1]);
          hang_weight = hang_info_pt->master_weight(m);
         }
        //If node is not hanging number or not
        else
         {
          local_eqn = nodal_local_eqn(l,u_nodal_index[2*i+1]);
          hang_weight = 1.0;
         }
        if(local_eqn >= 0)
         {
          residuals[local_eqn] += residual_contribution.imag()*hang_weight;
         }
        
      
        //CALCULATE THE JACOBIAN
        /*if(flag)
          {
          //Loop over the velocity shape functions again
          for(unsigned l2=0;l2<n_node;l2++)
          { 
          //Loop over the velocity components again
          for(unsigned i2=0;i2<DIM;i2++)
          {
          //If at a non-zero degree of freedom add in the entry
          local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
          if(local_unknown >= 0)
          {
          //Add contribution to Elemental Matrix
          jacobian(local_eqn,local_unknown) 
          -= visc_ratio*Gamma[i]*dpsifdx(l2,i)*dtestfdx(l,i2)*W;
          
          //Extra component if i2 = i
          if(i2 == i)
          {      
          //Loop over velocity components
          for(unsigned k=0;k<DIM;k++)
          { 
          jacobian(local_eqn,local_unknown)
          -= visc_ratio*dpsifdx(l2,k)*dtestfdx(l,k)*W;
          }
          }
          
          //Now add in the inertial terms
          jacobian(local_eqn,local_unknown)
          -= scaled_re*psif[l2]*interpolated_dudx(i,i2)*testf[l]*W;
          
          //Extra component if i2=i
          if(i2 == i)
          {
          //Add the mass matrix term (only diagonal entries)
          //Note that this is positive because the mass matrix
          //is taken to the other side of the equation when
          //formulating the generalised eigenproblem.
          if(flag==2)
          {
          mass_matrix(local_eqn,local_unknown) +=
          scaled_re_st*psif[l2]*testf[l]*W;
          }
          
          //du/dt term
          jacobian(local_eqn,local_unknown)
          -= scaled_re_st*
          node_pt(l2)->time_stepper_pt()->weight(1,0)*
          psif[l2]*testf[l]*W;  
          
          //Loop over the velocity components
          for(unsigned k=0;k<DIM;k++)
          {
          double tmp=scaled_re*interpolated_u[k];
          if (!ALE_is_disabled) tmp-=scaled_re_st*mesh_velocity[k];
          jacobian(local_eqn,local_unknown) -=
          tmp*dpsifdx(l2,k)*testf[l]*W;
          }
          }
          
          }
          }
          }
          
          //Now loop over pressure shape functions
          //This is the contribution from pressure gradient
          for(unsigned l2=0;l2<n_pres;l2++)
          {
          //If we are at a non-zero degree of freedom in the entry
          local_unknown = p_local_eqn(l2);
          if(local_unknown >= 0)
          {
          jacobian(local_eqn,local_unknown)
          += psip[l2]*dtestfdx(l,i)*W;
          }
          }
          } //End of Jacobian calculation
          
          }*/ //End of if not boundary condition statement
        
         } //End of loop over dimension
       } //End of loop over master nodes
     } //End of loop over shape functions
   
    
   
    //CONTINUITY EQUATION
    //-------------------
    
    //Loop over the shape functions
    for(unsigned l=0;l<n_pres;l++)
     {
      if(pressure_dof_is_hanging[l])
       {
        hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index[0]);
        n_master = hang_info_pt->nmaster();
       }
      else
       {
        n_master=1;
       }

      //Loop over the master nodes
      for(unsigned m=0;m<n_master;++m)
       {
        //Assemble the residuals
        std::complex<double> residual_contribution
         = interpolated_dudx[0][0];
        for(unsigned k=1;k<DIM;++k)
         {
          residual_contribution += interpolated_dudx[k][k];
         }

        if(pressure_dof_is_hanging[l])
         {
          local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
                                           p_index[0]);
          hang_weight = hang_info_pt->master_weight(m);
         }
        else
         {
          local_eqn = this->p_local_eqn(l,0);
          hang_weight = 1.0;
         }
        
      //If not a boundary conditions
        if(local_eqn >= 0)
         {
          residuals[local_eqn] +=
           residual_contribution.real()*testp[l]*W*hang_weight;
         }

        if(pressure_dof_is_hanging[l])
         {
          local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
                                           p_index[1]);
          hang_weight = hang_info_pt->master_weight(m);
         }
        else
         {
          local_eqn = this->p_local_eqn(l,1);
          hang_weight = 1.0;
         }

      //If not a boundary conditions
      if(local_eqn >= 0)
       {
        residuals[local_eqn] +=
         residual_contribution.imag()*testp[l]*W*hang_weight;
       }

       } //End of loop over master nodes
     } //End of loop over l

    //Normalisation condition. Leave this alone because there is
    //no test function involved.
    std::complex<double> residual_contribution =
     interpolated_p_normalisation*interpolated_p;
    for(unsigned k=0;k<DIM;++k)
     {
      residual_contribution +=
       interpolated_u_normalisation[k]*interpolated_u[k];
     }

    local_eqn = this->eigenvalue_local_eqn(0);
    if(local_eqn >= 0)
     {
      residuals[local_eqn] += residual_contribution.real()*W;
     }

    local_eqn = this->eigenvalue_local_eqn(1);
    if(local_eqn >= 0)
     {
      residuals[local_eqn] += residual_contribution.imag()*W;
     }
    
   } // End of loop over the integration points

 } // End of fill_in_generic_residual_contribution_linearised_nst
 
} // End of namespace oomph