axisym_solid_traction_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//    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//====================================================================
//Header file for elements that are used to apply surface loads to 
//the equations of elasticity

#ifndef OOMPH_AXISYMM_SOLID_TRACTION_ELEMENTS_HEADER
#define OOMPH_AXISYMM_SOLID_TRACTION_ELEMENTS_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

// OOMPH-LIB headers
#include "../generic/Qelements.h"
#include "../generic/hermite_elements.h"

namespace oomph
{

//======================================================================
/// A class for elements that allow the imposition of an applied traction
/// in the principle of virtual displacements. 
/// The geometrical information can be read from the FaceGeometry<ELEMENT> 
/// class and and thus, we can be generic enough without the need to have
/// a separate equations class.
//======================================================================
template <class ELEMENT>
class AxisymmetricSolidTractionElement : 
public virtual FaceGeometry<ELEMENT>, public virtual FaceElement
{
private:

 /// Pointer to an imposed traction function
 void (*Traction_fct_pt)(const double& time, const Vector<double> &x, 
                         const Vector<double> &n,
                         Vector<double> &result);

protected:

 /// Return the surface traction force
 void get_traction(const double& time, const Vector<double> &x,
                   const Vector<double> &n,Vector<double> &result) const
  {
   //If the function pointer is zero return zero
   if(Traction_fct_pt == 0)
    {
     //Loop over dimensions and set body forces to zero
     //It's axisymmetric, so we only have "two" dimensions
     for(unsigned i=0;i<2;i++) {result[i] = 0.0;}
    }
   //Otherwise call the function
   else
    {
     (*Traction_fct_pt)(time,x,n,result);
    }
  }

public:

 /// \short Constructor, which takes a "bulk" element and 
 /// the value of the index and its limit
 AxisymmetricSolidTractionElement(FiniteElement* const &element_pt, 
                                  const int &face_index) : 
  FaceGeometry<ELEMENT>(), FaceElement()
  { 
   //Attach the geometrical information to the element. N.B. This function
   //also assigns nbulk_value from the required_nvalue of the bulk element
   element_pt->build_face_element(face_index,this);
 
   //Set the body force function pointer to zero
   Traction_fct_pt = 0;
  }
 
 /// Return the imposed traction pointer
 void (* &traction_fct_pt())(const double&, const Vector<double> &, 
                             const Vector<double> &, Vector<double> &) 
  {return Traction_fct_pt;}
 
 /// Return the residuals
 void fill_in_contribution_to_residuals(Vector<double> &residuals);
 
/// Return the jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals, 
                                   DenseMatrix<double> &jacobian)
  {
   fill_in_contribution_to_residuals(residuals);
   //Call the generic FD jacobian calculation
   FaceGeometry<ELEMENT>::fill_in_jacobian_from_solid_position_by_fd(jacobian);
  }

 /// Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}

 /// Output function: x,y,[z],u,v,[w],p in tecplot format
 void output(std::ostream &outfile, const unsigned &n_plot)
  {FiniteElement::output(outfile,n_plot);}

 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}

 /// Output function: x,y,[z],u,v,[w],p in tecplot format
 void output(FILE* file_pt, const unsigned &n_plot)
  {FiniteElement::output(file_pt,n_plot);}

}; 



/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////


//=======================================================================
/// Return the residuals for the AxisymmetricSolidTractionElements
//=======================================================================
template<class ELEMENT>
void AxisymmetricSolidTractionElement<ELEMENT>::
fill_in_contribution_to_residuals(Vector<double> &residuals)
{
 //Find out how many nodes there are
 unsigned n_node = nnode();
 //Find out how many positional dofs there are
 unsigned n_position_type = nnodal_position_type();

 //Integer to hold the local equation number
 int local_eqn=0;

 //Set up memory for the shape functions
 //The surface is 1D, so we only have one local derivative
 Shape psi(n_node,n_position_type);
 DShape dpsids(n_node,n_position_type,1); 

 //Set the value of n_intpt
 unsigned n_intpt = integral_pt()->nweight();

 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   //Only need to call the local derivatives
   dshape_local_at_knot(ipt,psi,dpsids);

   //Calculate the global position and lagrangian coordinate
   Vector<double> interpolated_x(2,0.0), interpolated_xi(2,0.0);
   //Calculate the global and lagrangian derivtives wrt the local coordinates
   Vector<double> interpolated_dxds(2,0.0), interpolated_dxids(2,0.0);

   //Calculate displacements and derivatives
   for(unsigned l=0;l<n_node;l++) 
    {
     //Loop over positional dofs
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Loop over the number of lagrangian coordinates (2)
       for(unsigned i=0;i<2;i++)
        {
         //Calculate the global position
         interpolated_x[i] += 
          nodal_position_gen(l,bulk_position_type(k),i)*psi(l,k);
         interpolated_xi[i] += 
          this->lagrangian_position_gen(l,bulk_position_type(k),i)*psi(l,k);
         //Calculate the derivatives of the global and lagrangian coordinates
         interpolated_dxds[i] += 
          nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,0);
         interpolated_dxids[i] += 
          this->lagrangian_position_gen(l,bulk_position_type(k),i)
          *dpsids(l,k,0);
        }
      }
    }

   //Now calculate the entries of the deformed surface metric tensor
   //Now find the local deformed metric tensor from the tangent Vectors
   DenseMatrix<double> A(2);
   //The off-diagonal terms are Zero 
   A(0,1) = A(1,0) = 0.0;
   //The diagonal terms are a little complicated
   A(0,0) =  
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1])*
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1]) +
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1])*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);


   A(1,1) =  (interpolated_x[0]*sin(interpolated_xi[1]) +
               interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]));

   //Premultiply the weights and the square-root of the determinant of 
   //the metric tensor
   double W = w*sqrt(A(0,0)*A(1,1));
   
   //Also find the normal -- just the cross product of the metric tensors
   //but I want to express it in terms of e_r and e_theta components
   //N.B. There is an issue at theta = 0,pi, where the normal is e_{r},
   //but given that I never assemble it, should be OK!
   //The minus sign is chosen to ensure that the normal is really outward 
   Vector<double> interpolated_normal(2);
   //Component in the e_{r} direction
   interpolated_normal[0] = -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);
   //Component in the e_{theta} direction
   interpolated_normal[1] =  -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[1]*interpolated_dxids[1] - interpolated_dxds[0]);

   //If we're on the north or south face need to flip normal
   if(s_fixed_value()==-1)
    {
     interpolated_normal[0] *= -1.0;
     interpolated_normal[1] *= -1.0;
    }

   //Now adjust and scale the normal
   double length = 0.0;
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] *= normal_sign();
     length += interpolated_normal[i]*interpolated_normal[i];
    }
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] /= sqrt(length);
    }

   //Now calculate the load
   Vector<double> traction(2);

   //Normal is outwards
   get_traction(time(),interpolated_x,interpolated_normal,traction);

   //=====LOAD TERMS  FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
       
   //Loop over the test functions, nodes of the element
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop of types of dofs
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Loop over the displacement components
       for(unsigned i=0;i<2;i++)
        {
         local_eqn = 
          this->position_local_eqn(l,bulk_position_type(k),i);
         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
           //Add the loading terms to the residuals
           residuals[local_eqn] -= traction[i]*psi(l,k)*W;
          }
        }
      } //End of if not boundary condition
    } //End of loop over shape functions
  } //End of loop over integration points

}


}

#endif