young_laplace_contact_angle_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 "young_laplace_contact_angle_elements.h"
#include "young_laplace_elements.h"
#include "refineable_young_laplace_elements.h"


namespace oomph
{

//===========================================================================
/// Constructor, takes the pointer to the "bulk" element  and the 
/// index of the face to which the element is attached.
//===========================================================================
 template<class ELEMENT>
 YoungLaplaceContactAngleElement<ELEMENT>::
 YoungLaplaceContactAngleElement(FiniteElement* const &bulk_el_pt, 
                                 const int &face_index) : 
  FaceGeometry<ELEMENT>(), FaceElement()
 { 
  
  // Let the bulk element build the FaceElement, i.e. setup the pointers 
  // to its nodes (by referring to the appropriate nodes in the bulk
  // element), etc.
  bulk_el_pt->build_face_element(face_index,this);
  
  // Initialise the prescribed contact angle pointer to zero
  Prescribed_cos_gamma_pt = 0;
  
#ifdef PARANOID
  // Extract the dimension of the problem from the dimension of 
  // the first node
  unsigned dim = node_pt(0)->ndim();
  if (dim!=2) 
   {
    throw OomphLibError(
     "YoungLaplaceContactAngleElement only work with 2D nodes",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
 }





//================================================================
/// Compute the element's contribution to the residual vector 
//================================================================
 template<class ELEMENT>
 void YoungLaplaceContactAngleElement<ELEMENT>::
 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 functions
  Shape psi(n_node);
  
  //Number of integration points
  unsigned n_intpt = integral_pt()->nweight();
  
  // Set the Vector to hold local coordinates
  // Note: We need the coordinate itself below for the evaluation
  // of the contact line vectors even though we use the *at_knot
  // version for the various shape-function-related functions
  Vector<double> s(1); 
  
  //Integers to hold the local equation and unknown numbers
  int local_eqn=0;
  
  //Loop over the integration points
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {
    //Assign value of s
    s[0] = integral_pt()->knot(ipt,0);
   
   //Get the integral weight
   double w = integral_pt()->weight(ipt); 

   //Find the shape functions
   shape_at_knot(ipt,psi); 

   //Get the prescribed cos_gamma
   double cos_gamma=prescribed_cos_gamma();

   // Get the various contact line vectors
   Vector<double> tangent(3);
   Vector<double> normal(3);
   Vector<double> spine(3);
   double norm_of_drds;
   contact_line_vectors(s,tangent,normal,spine,norm_of_drds);

   // Get beta factor:

   // Cross product of spine and tangent to contact line is
   // the wall normal
   Vector<double> wall_normal(3);
   ELEMENT::cross_product(spine,tangent,wall_normal);
 
   // Normalise
   double norm=ELEMENT::two_norm(wall_normal);
   for (unsigned i=0;i<3;i++) wall_normal[i]/=norm;
   
   // Take cross product with tangent to get the normal to the
   // contact line parallel to wall
   Vector<double> normal_to_contact_line_parallel_to_wall(3);
   ELEMENT::cross_product(tangent,wall_normal,
                          normal_to_contact_line_parallel_to_wall);

   double beta=ELEMENT::scalar_product(spine,
                                     normal_to_contact_line_parallel_to_wall);

                                    
   //Now add to the appropriate equations
   
   //Loop over the test functions
   for(unsigned l=0;l<n_node;l++)
    {
     local_eqn = u_local_eqn(l);
     /*IF it's not a boundary condition*/
     if(local_eqn >= 0)
      {
       //Add to residual:
       residuals[local_eqn] -= beta*cos_gamma*psi[l]*norm_of_drds*w;
      }
    }
  }

}




//========================================================================
/// Get the actual contact angle
//========================================================================
 template<class ELEMENT>
 double YoungLaplaceContactAngleElement<ELEMENT>::actual_cos_contact_angle(
  const Vector<double>& s)
 {  
  // Get pointer to bulk element
  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
  
  double cos_gamma=0.0;
  
  // Spine case
  if ( bulk_elem_pt->use_spines() )
   {
    // Get various contact line vectors
    Vector<double> tangent(3);
    Vector<double> normal(3);
    Vector<double> spine(3);
    double norm_of_drds;
    contact_line_vectors(s,tangent,normal,spine,norm_of_drds);
    
    // Get the wall normal: Both the tangent to the contact line
    // and the spine vector are tangential to the wall:
    Vector<double> wall_normal(3); 
    Vector<double> axe_ez(3,0.0);
    axe_ez[2]=1.0;
    ELEMENT::cross_product(axe_ez,tangent,wall_normal);
    
    // Normalise
    double norm=0.0;
    for (unsigned i=0;i<3;i++) norm+=wall_normal[i]*wall_normal[i];
    for (unsigned i=0;i<3;i++)
     { 
      wall_normal[i]/=sqrt(norm);
     }
    
    // Get the auxiliary unit vector that's normal to 
    // the contact line and tangent to the wall
    Vector<double> aux(3);
    ELEMENT::cross_product(tangent,wall_normal,aux);
    
    // Cos of contact angle is dot product with wall normal
    cos_gamma=ELEMENT::scalar_product(aux,normal);
    
   }
  
  // Cartesian case
  else
   {
    
    // Get local coordinates in bulk element by copy construction
    Vector<double> s_bulk(local_coordinate_in_bulk(s));
    
    // Number of nodes in bulk element
    unsigned nnode_bulk=bulk_elem_pt->nnode();
    
    
#ifdef PARANOID
    // Dimension of (= number of local coordinates in) bulk element 
    unsigned dim_bulk=bulk_elem_pt->dim();

    if (dim_bulk!=2) 
     {
      throw OomphLibError(
       "YoungLaplaceContactAngleElements only work with 2D bulk elements",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
#endif
    
    //Set up memory for the shape functions
    Shape psi(nnode_bulk);
    DShape dpsidzeta(nnode_bulk,2);
    
    //Call the derivatives of the shape and test functions
    //in the bulk -- must do this via s because this point
    //is not an integration point the bulk element!
    (void)bulk_elem_pt->dshape_eulerian(s_bulk,psi,dpsidzeta);
    
    // Get the gradient at s
    Vector<double> gradient_u(2,0.0);
    
    //Calculate function value and derivatives:
    //-----------------------------------------
    // Loop over nodes
    for(unsigned l=0;l<nnode_bulk;l++) 
     {
      // Loop over directions
      for(unsigned j=0;j<2;j++)
       { gradient_u[j]+= bulk_elem_pt->u(l)*dpsidzeta(l,j); }
     }
    
    // Get the outer unit normal to boundary
    Vector<double> outer_normal(2,0.0);
    outer_unit_normal(s,outer_normal);
    
    // Compute the cosinus of the angle
    double gradient_norm_2=
     ELEMENT::two_norm(gradient_u)*
     ELEMENT::two_norm(gradient_u);
    cos_gamma=ELEMENT::scalar_product(gradient_u,outer_normal)/
     sqrt(1+gradient_norm_2);
    
   }                                                     
    
 return cos_gamma;

}



//========================================================================
/// Get unit tangent and normal to contact line and the spine itself (this 
/// allows the wall normal to be worked out by a cross product). Final
/// argument gives the norm of dR/ds where R is the vector to the
/// contact line and s the local coordinate in the element.
//========================================================================
template<class ELEMENT>
void YoungLaplaceContactAngleElement<ELEMENT>::contact_line_vectors(
 const Vector<double>& s, 
 Vector<double>& tangent, 
 Vector<double>& normal,
 Vector<double>& spine,
 double& norm_of_drds)
{
 
 // Get pointer to bulk element
 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
 
 // Dimension of (= number of local coordinates in) bulk element 
 unsigned dim_bulk=bulk_elem_pt->dim();
 
#ifdef PARANOID
 if (dim_bulk!=2) 
  {
   throw OomphLibError(
    "YoungLaplaceContactAngleElements only work with 2D bulk elements",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Which local coordinate is const in bulk element?
 unsigned s_fixed_index_in_bulk;
 
 DenseMatrix<double> ds_bulk_ds_face(2,1);
 this->get_ds_bulk_ds_face(s,ds_bulk_ds_face,
                           s_fixed_index_in_bulk);
 
 // Get local coordinates in bulk element by copy construction
 Vector<double> s_bulk(local_coordinate_in_bulk(s));
 
 // Number of nodes in bulk element
 unsigned nnode_bulk=bulk_elem_pt->nnode();

 // Get the shape functions and their derivatives w.r.t. 
 // to the local coordinates in the bulk element
 Shape psi_bulk(nnode_bulk);
 DShape dpsi_bulk(nnode_bulk,dim_bulk);
 bulk_elem_pt->dshape_local(s_bulk,psi_bulk,dpsi_bulk);
 
 // Displacement along spine
 double interpolated_u=0.0;

 // Derivative of u w.r.t. tangential and pseudo-normal bulk coordinate
 double interpolated_du_ds_tangent=0.0;
 double interpolated_du_ds_pseudo_normal=0.0;
 
 // Global intrinsic coordinates of current point for evaluation of
 // spine and spine base
 Vector<double> interpolated_zeta(dim_bulk,0.0);
 

 // Derivatives of zeta (in bulk) w.r.t. to tangent and pseudo-normal
 // local coordinate
 Vector<double> interpolated_dzeta_ds_tangent(dim_bulk);
 Vector<double> interpolated_dzeta_ds_pseudo_normal(dim_bulk);
 
 // Loop over nodes                                                           
 for(unsigned l=0;l<nnode_bulk;l++)
  {                    
   interpolated_u+=bulk_elem_pt->u(l)*psi_bulk(l);
   
   // Chain rule for tangent derivative
   double aux=0.0;
   for (unsigned i=0;i<dim_bulk;i++)
    {
     aux+=dpsi_bulk(l,i)*ds_bulk_ds_face(i,0);
    }
   interpolated_du_ds_tangent+=bulk_elem_pt->u(l)*aux;
   
   // Straight derivative w.r.t. one of the bulk coordinates
   // that's fixed along the face
   interpolated_du_ds_pseudo_normal += 
    bulk_elem_pt->u(l)*dpsi_bulk(l,s_fixed_index_in_bulk);
   
   // Loop over directions                                                    
   for(unsigned j=0;j<dim_bulk;j++)
    {                    
     interpolated_zeta[j]+=bulk_elem_pt->nodal_position(l,j)*psi_bulk(l);
     
     // Chain rule for tangent derivative
     double aux=0.0;
     for (unsigned i=0;i<dim_bulk;i++)
      {
       aux+=dpsi_bulk(l,i)*ds_bulk_ds_face(i,0);
      }
     interpolated_dzeta_ds_tangent[j] +=
      bulk_elem_pt->nodal_position(l,j)*aux;
     
     // Straight derivative w.r.t. one of the bulk coordinates
     // that's fixed along the face
     interpolated_dzeta_ds_pseudo_normal[j] +=
      bulk_elem_pt->nodal_position(l,j)*dpsi_bulk(l,s_fixed_index_in_bulk);
     
    }
  }                                                                         
 
 // Auxiliary vector (tangent to non-fixed bulk coordinate but
 // not necessarily normal to contact line)
 Vector<double> aux_vector(3);
 double tang_norm=0.0;
 double aux_norm=0.0;
 
 if ( bulk_elem_pt->use_spines() )
  {
   // Get the spine and spine base vector at this point from the bulk element
   Vector<double> spine_base(3,0.0);
   Vector< Vector<double> > dspine;
   ELEMENT::allocate_vector_of_vectors(2,3,dspine);
   Vector< Vector<double> > dspine_base;
   ELEMENT::allocate_vector_of_vectors(2,3,dspine_base);
   bulk_elem_pt->get_spine(interpolated_zeta,spine,dspine);
   bulk_elem_pt->get_spine_base(interpolated_zeta,spine_base,dspine_base);
   
   // Derivative of spine and spine base w.r.t. local coordinate in 
   // FaceElement:
   Vector<double> dspine_ds_tangent(3,0.0);
   Vector<double> dspine_base_ds_tangent(3,0.0);
   Vector<double> dspine_ds_pseudo_normal(3,0.0);
   Vector<double> dspine_base_ds_pseudo_normal(3,0.0);
   for (unsigned i=0;i<3;i++)
    {
     
     dspine_ds_tangent[i]+=
      dspine[0][i]*interpolated_dzeta_ds_tangent[0]+
      dspine[1][i]*interpolated_dzeta_ds_tangent[1];
     
     dspine_base_ds_tangent[i]+=
      dspine_base[0][i]*interpolated_dzeta_ds_tangent[0]+
      dspine_base[1][i]*interpolated_dzeta_ds_tangent[1];
          
     dspine_ds_pseudo_normal[i]+=
      dspine[0][i]*interpolated_dzeta_ds_pseudo_normal[0]+
      dspine[1][i]*interpolated_dzeta_ds_pseudo_normal[1];
     
     dspine_base_ds_pseudo_normal[i]+=
      dspine_base[0][i]*interpolated_dzeta_ds_pseudo_normal[0]+
      dspine_base[1][i]*interpolated_dzeta_ds_pseudo_normal[1];
    }
   
   // Auxiliary vector (tangent to non-fixed bulk coordinate but
   // not necessarily normal to contact line)
   for (unsigned i=0;i<3;i++)
    {
     tangent[i]=
      dspine_base_ds_tangent[i]+interpolated_du_ds_tangent*spine[i]+
      interpolated_u*dspine_ds_tangent[i];
     tang_norm+=tangent[i]*tangent[i];
     
     
     aux_vector[i]=
      dspine_base_ds_pseudo_normal[i]+
      interpolated_du_ds_pseudo_normal*spine[i]+
      interpolated_u*dspine_ds_pseudo_normal[i];
     aux_norm+=aux_vector[i]*aux_vector[i];
    }
   
  }
 
 // Cartesian case
 else
  {
   for(unsigned i=0;i<2;i++)
    {
     tangent[i]= interpolated_dzeta_ds_tangent[i];
     tang_norm+=tangent[i]*tangent[i];
     aux_vector[i]= interpolated_dzeta_ds_pseudo_normal[i];
     aux_norm+=aux_vector[i]*aux_vector[i];
     
    }
   
   tangent[2]=interpolated_du_ds_tangent; 
   tang_norm+=tangent[2]*tangent[2];
   aux_vector[2]=interpolated_du_ds_pseudo_normal; 
   aux_norm+=aux_vector[2]*aux_vector[2];
   
  }
    
 // Norm of the rate of change
 norm_of_drds=sqrt(tang_norm);
 
 // Normalise
 double tang_norm_fact=1.0/sqrt(tang_norm);
 double aux_norm_fact=1.0/sqrt(aux_norm);
 for (unsigned i=0;i<3;i++)
   {
     tangent[i]*=tang_norm_fact;
     aux_vector[i]*=aux_norm_fact;
   }
 
 
 // Normal to meniscus is the cross product between the
 // two contact line vectors:
 Vector<double> meniscus_normal(3);
 ELEMENT::cross_product(tangent,aux_vector,meniscus_normal);
 
 // Calculate the norm
 double men_norm_fact=0.0;
 for (unsigned i=0;i<3;i++)
   {
     men_norm_fact+=meniscus_normal[i]*meniscus_normal[i];
   }
 
 // Normalise and adjust direction
 double sign=-double(normal_sign());
 for (unsigned i=0;i<3;i++)
  {
   meniscus_normal[i]*=sign/sqrt(men_norm_fact);
  }
 
 // The actual (bi) normal to the contact line is given
 // by another cross product
 ELEMENT::cross_product(meniscus_normal,tangent,normal);
 
}





//============================================================
// Build the required elements
//============================================================
template class YoungLaplaceContactAngleElement<QYoungLaplaceElement<2> >;
template class YoungLaplaceContactAngleElement<QYoungLaplaceElement<3> >;
template class YoungLaplaceContactAngleElement<QYoungLaplaceElement<4> >;

template class YoungLaplaceContactAngleElement<RefineableQYoungLaplaceElement<2> >;
template class YoungLaplaceContactAngleElement<RefineableQYoungLaplaceElement<3> >;
template class YoungLaplaceContactAngleElement<RefineableQYoungLaplaceElement<4> >;

}