refineable_young_laplace_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$
//LIC//
//LIC// $LastChangedDate$
//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//====================================================================
#include "refineable_young_laplace_elements.h"


namespace oomph
{


//======================================================================
/// Compute element residual vector taking hanging nodes into account
//======================================================================
void RefineableYoungLaplaceEquations::
fill_in_contribution_to_residuals(Vector<double> &residuals)
{
 
 //Find out how many nodes there are
 unsigned n_node = nnode();
 
 //Set up memory for the shape (test) functions
 Shape psi(n_node);
 DShape dpsidzeta(n_node,2);
 
 //Set the value of n_intpt
 unsigned n_intpt = integral_pt()->nweight();
 
 //Integers to store the local equation numbers
 int local_eqn=0;

 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   //Call the derivatives of the shape (test) functions
   double J = dshape_eulerian_at_knot(ipt,psi,dpsidzeta);
       
   //Premultiply the weights and the Jacobian
   double W = w*J;

   //Calculate local values of the pressure and velocity components
   //Allocate and initialise to zero
   double interpolated_u=0.0;
   Vector<double> interpolated_zeta(2,0.0);
   Vector<double> interpolated_du_dzeta(2,0.0);
   
   //Calculate function value and derivatives:
   //-----------------------------------------
   // Loop over nodes
   for(unsigned l=0;l<n_node;l++) 
    {
     interpolated_u += u(l)*psi(l);
     // Loop over directions
     for(unsigned j=0;j<2;j++)
      {
       interpolated_zeta[j] += nodal_position(l,j)*psi(l);
       interpolated_du_dzeta[j] += u(l)*dpsidzeta(l,j);
      }
    }

   // Allocation and definition of variables necessary for
   // further calculations
   
   /// "Simple" case
   ///--------------
   double nonlinearterm=1.0;
   double sqnorm=0.0;
  
   /// Spine case 
   ///-----------

  // Derivs of position vector w.r.t. global intrinsic coords 
   Vector<Vector<double> > dRdzeta; 
   allocate_vector_of_vectors(2,3,dRdzeta);
   
   // Unnormalised normal
   Vector<double> N_unnormalised(3,0.0);
   
   // Spine and spine basis vectors, entries initialised to zero
   Vector<double> spine_base(3,0.0), spine(3,0.0);

   // Derivative of spine basis vector w.r.t to the intrinsic
   // coordinates: dspine_base[i,j] = j-th component of the deriv.
   // of the spine basis vector w.r.t. to the i-th global intrinsic
   // coordinate
   Vector< Vector<double> > dspine_base; 
   allocate_vector_of_vectors(2,3,dspine_base);

   // Derivative of spine vector w.r.t to the intrinsic
   // coordinates: dspine[i,j] = j-th component of the deriv.
   // of the spine vector w.r.t. to the i-th global intrinsic
   // coordinate
   Vector< Vector<double> > dspine; 
   allocate_vector_of_vectors(2,3,dspine);
   
   // Vector v_\alpha contains the numerator of the variations of the 
   // area element {\cal A}^{1/2} w.r.t. the components of dR/d\zeta_\alpha 
   Vector<double> area_variation_numerator_0(3,0.0);
   Vector<double> area_variation_numerator_1(3,0.0);

   // Vector position
   Vector<double> r(3,0.0);

   //No spines
   //---------
   if (!use_spines())
    {
    
     for (unsigned j=0; j<2; j++) 
      sqnorm += interpolated_du_dzeta[j]*interpolated_du_dzeta[j];
     
     nonlinearterm=1.0/sqrt(1.0+sqnorm);
    }

   //Spines
   //------
   else
    {
     // Get the spines
     get_spine_base(interpolated_zeta, spine_base, dspine_base);
     get_spine(interpolated_zeta, spine, dspine);

     // calculation of dR/d\zeta_\alpha
     for (unsigned alpha=0;alpha<2;alpha++)
      {
       // Product rule for d(u {\bf S} ) / d \zeta_\alpha
       Vector<double> dudzeta_times_spine(3,0.0);
       scalar_times_vector(interpolated_du_dzeta[alpha],
                           spine,dudzeta_times_spine);

       Vector<double> u_times_dspinedzeta(3,0.0);
       scalar_times_vector(interpolated_u,dspine[alpha],u_times_dspinedzeta);

       Vector<double> d_u_times_spine_dzeta(3,0.0);
       vector_sum(dudzeta_times_spine,
                  u_times_dspinedzeta,
                  d_u_times_spine_dzeta);

       // Add derivative of spine base
       vector_sum(d_u_times_spine_dzeta,dspine_base[alpha],dRdzeta[alpha]); 
      }

     /// Get the unnormalized normal
     cross_product(dRdzeta[0],dRdzeta[1],N_unnormalised);
          
     /// Tmp storage 
     Vector<double> v_tmp_1(3,0.0);
     Vector<double> v_tmp_2(3,0.0);

     // Calculation of 
     // |dR/d\zeta_1|^2 dR/d\zeta_0 - <dR/d\zeta_0,dR/d\zeta_1>dR/d\zeta_1
     scalar_times_vector(pow(two_norm(dRdzeta[1]),2), dRdzeta[0], v_tmp_1);
     scalar_times_vector(-1*scalar_product(dRdzeta[0],dRdzeta[1]), 
                         dRdzeta[1], v_tmp_2);
     vector_sum(v_tmp_1,v_tmp_2,area_variation_numerator_0);

     // Calculation of 
     // |dR/d\zeta_0|^2 dR/d\zeta_1 - <dR/d\zeta_0,dR/d\zeta_1>dR/d\zeta_0
     scalar_times_vector(pow(two_norm(dRdzeta[0]),2), dRdzeta[1], v_tmp_1);
     scalar_times_vector(-1*scalar_product(dRdzeta[0],dRdzeta[1]), 
                         dRdzeta[0], v_tmp_2);
     vector_sum(v_tmp_1,v_tmp_2,area_variation_numerator_1);  

     // Global Eulerian cooordinates
     for (unsigned j=0; j<3; j++)
       {
         r[j]=spine_base[j]+interpolated_u*spine[j];
       }
     
    }
   
   
   
   // Assemble residuals
   //-------------------

   // Loop over the nodes for the test functions 
   for(unsigned l=0;l<n_node;l++)
    {
     //Local variables used to store the number of master nodes and the
     //weight associated with the shape function if the node is hanging
     unsigned n_master=1; double hang_weight=1.0;
     
     //If the node is hanging, get the number of master nodes
     if(node_pt(l)->is_hanging())
      {
       n_master = node_pt(l)->hanging_pt()->nmaster();
      }
     //Otherwise there is just one master node, the node itself
     else
      {
       n_master = 1;
      }
     
     //Loop over the master nodes
     for(unsigned m=0;m<n_master;m++)
      {
       //Get the local equation number and hang_weight
       //If the node is hanging
       if(node_pt(l)->is_hanging())
        {
         //Read out the local equation from the master node
         local_eqn = 
          local_hang_eqn(node_pt(l)->hanging_pt()->master_node_pt(m),0);
         //Read out the weight from the master node
         hang_weight = node_pt(l)->hanging_pt()->master_weight(m);
        }
       //If the node is not hanging
       else
        {
         //The local equation number comes from the node itself
         local_eqn = this->u_local_eqn(l);
         //The hang weight is one
         hang_weight = 1.0;
        }
       
       /*IF it's not a boundary condition*/
       if(local_eqn >= 0)
        {
         
         // "simple" calculation case
         if (!use_spines())
          { 
           // Add source term: The curvature
           residuals[local_eqn] += get_kappa()*psi(l)*W*hang_weight; 
           
           // The YoungLaplace bit itself
           for(unsigned k=0;k<2;k++)
            {
             residuals[local_eqn] += nonlinearterm*
              interpolated_du_dzeta[k]*dpsidzeta(l,k)*W*hang_weight;
            }
          }
         
         // Spine calculation case
         else 
          {

            // Calculation of d(u S)/d\zeta_0
            //-------------------------------
            Vector<double> v_tmp_1(3,0.0);
            scalar_times_vector(dpsidzeta(l,0), spine, v_tmp_1);
            
            Vector<double> v_tmp_2(3,0.0);
            scalar_times_vector(psi(l),dspine[0], v_tmp_2);

            Vector<double> d_uS_dzeta0(3,0.0);
            vector_sum(v_tmp_1,v_tmp_2,d_uS_dzeta0);

            // Add contribution to residual 
            residuals[local_eqn] += 
              W*hang_weight*
              scalar_product(area_variation_numerator_0,d_uS_dzeta0)
              /two_norm(N_unnormalised);

            // Calculation of d(u S)/d\zeta_1
            scalar_times_vector(dpsidzeta(l,1), spine, v_tmp_1);
            scalar_times_vector(psi(l),dspine[1], v_tmp_2);
            Vector<double> d_uS_dzeta1(3,0.0);
            vector_sum(v_tmp_1,v_tmp_2,d_uS_dzeta1);

            // Add contribution to residual
            residuals[local_eqn] += 
              W*hang_weight*
              scalar_product(area_variation_numerator_1,d_uS_dzeta1)
              /two_norm(N_unnormalised);    

           // Curvature contribution to the residual  : kappa N S test
           residuals[local_eqn] += W*hang_weight*(get_kappa())*
             scalar_product(N_unnormalised,spine)*psi(l);
          }
        }
       
      } //End of loop over master nodes for residual
     
    } //End of loop over nodes
   
  } // End of loop over integration points

}
 
//====================================================================
// Force build of templates
//====================================================================
template class RefineableQYoungLaplaceElement<2>;
template class RefineableQYoungLaplaceElement<3>;
template class RefineableQYoungLaplaceElement<4>;

}