text: 
  
  // System includes
  
  // External includes
  
  // Project includes
  #include "custom_elements/solid_elements/axisymmetric_updated_lagrangian_element.hpp"
  #include "solid_mechanics_application_variables.h"
  
  
  namespace Kratos
  {
  
  //******************************CONSTRUCTOR*******************************************
  //************************************************************************************
  
  AxisymmetricUpdatedLagrangianElement::AxisymmetricUpdatedLagrangianElement( IndexType NewId, GeometryType::Pointer pGeometry )
      : LargeDisplacementElement( NewId, pGeometry )
  {
      //DO NOT ADD DOFS HERE!!!
  }
  
  
  //******************************CONSTRUCTOR*******************************************
  //************************************************************************************
  
  AxisymmetricUpdatedLagrangianElement::AxisymmetricUpdatedLagrangianElement( IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties )
      : LargeDisplacementElement( NewId, pGeometry, pProperties )
  {
      mThisIntegrationMethod = GetGeometry().GetDefaultIntegrationMethod();
      //mThisIntegrationMethod = GeometryData::GI_GAUSS_1;
      //mThisIntegrationMethod = GeometryData::GI_GAUSS_2;
  }
  
  
  //******************************COPY CONSTRUCTOR**************************************
  //************************************************************************************
  
  AxisymmetricUpdatedLagrangianElement::AxisymmetricUpdatedLagrangianElement( AxisymmetricUpdatedLagrangianElement const& rOther)
      :LargeDisplacementElement(rOther)
      ,mDeformationGradientF0(rOther.mDeformationGradientF0)
      ,mDeterminantF0(rOther.mDeterminantF0)
  {
  }
  
  
  //*******************************ASSIGMENT OPERATOR***********************************
  //************************************************************************************
  
  AxisymmetricUpdatedLagrangianElement&  AxisymmetricUpdatedLagrangianElement::operator=(AxisymmetricUpdatedLagrangianElement const& rOther)
  {
      LargeDisplacementElement::operator=(rOther);
  
      mDeformationGradientF0.clear();
      mDeformationGradientF0.resize(rOther.mDeformationGradientF0.size());
  
      for(unsigned int i=0; i<mConstitutiveLawVector.size(); i++)
      {
          mDeformationGradientF0[i] = rOther.mDeformationGradientF0[i];
      }
  
      mDeterminantF0 = rOther.mDeterminantF0;
  
      return *this;
  }
  
  
  //*********************************OPERATIONS*****************************************
  //************************************************************************************
  
  Element::Pointer AxisymmetricUpdatedLagrangianElement::Create( IndexType NewId, NodesArrayType const& rThisNodes, PropertiesType::Pointer pProperties ) const
  {
      return Kratos::make_intrusive< AxisymmetricUpdatedLagrangianElement >(NewId, GetGeometry().Create(rThisNodes), pProperties);
  }
  
  
  //************************************CLONE*******************************************
  //************************************************************************************
  
  Element::Pointer AxisymmetricUpdatedLagrangianElement::Clone( IndexType NewId, NodesArrayType const& rThisNodes ) const
  {
  
      AxisymmetricUpdatedLagrangianElement NewElement(NewId, GetGeometry().Create( rThisNodes ), pGetProperties() );
  
      //-----------//
  
      NewElement.mThisIntegrationMethod = mThisIntegrationMethod;
  
      if ( NewElement.mConstitutiveLawVector.size() != mConstitutiveLawVector.size() )
        {
  	NewElement.mConstitutiveLawVector.resize(mConstitutiveLawVector.size());
  
  	if( NewElement.mConstitutiveLawVector.size() != NewElement.GetGeometry().IntegrationPointsNumber() )
  	  KRATOS_THROW_ERROR( std::logic_error, "constitutive law not has the correct size ", NewElement.mConstitutiveLawVector.size() );
        }
  
  
      for(unsigned int i=0; i<mConstitutiveLawVector.size(); i++)
        {
  	NewElement.mConstitutiveLawVector[i] = mConstitutiveLawVector[i]->Clone();
        }
  
      //-----------//
  
      if ( NewElement.mDeformationGradientF0.size() != mDeformationGradientF0.size() )
        NewElement.mDeformationGradientF0.resize(mDeformationGradientF0.size());
  
      for(unsigned int i=0; i<mDeformationGradientF0.size(); i++)
      {
          NewElement.mDeformationGradientF0[i] = mDeformationGradientF0[i];
      }
  
      NewElement.mDeterminantF0 = mDeterminantF0;
  
      NewElement.SetData(this->GetData());
      NewElement.SetFlags(this->GetFlags());
  
      return Kratos::make_intrusive< AxisymmetricUpdatedLagrangianElement >(NewElement);
  }
  
  
  //*******************************DESTRUCTOR*******************************************
  //************************************************************************************
  
  AxisymmetricUpdatedLagrangianElement::~AxisymmetricUpdatedLagrangianElement()
  {
  }
  
  //************************************************************************************
  //************************************************************************************
  
  
  //*********************************SET DOUBLE VALUE***********************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::SetValuesOnIntegrationPoints( const Variable<double>& rVariable,
          const std::vector<double>& rValues,
          const ProcessInfo& rCurrentProcessInfo )
  {
  
    if (rVariable == DETERMINANT_F){
  
      const unsigned int& integration_points_number = mConstitutiveLawVector.size();
  
  
      for ( unsigned int PointNumber = 0;  PointNumber < integration_points_number; PointNumber++ )
        {
  	mDeterminantF0[PointNumber] = rValues[PointNumber];
  	mConstitutiveLawVector[PointNumber]->SetValue(rVariable, rValues[PointNumber], rCurrentProcessInfo);
        }
  
    }
    else{
  
      LargeDisplacementElement::SetValuesOnIntegrationPoints( rVariable, rValues, rCurrentProcessInfo );
  
    }
  
  
  }
  
  
  //************************************************************************************
  //************************************************************************************
  
  //**********************************GET DOUBLE VALUE**********************************
  //************************************************************************************
  
  
  void AxisymmetricUpdatedLagrangianElement::GetValueOnIntegrationPoints( const Variable<double>& rVariable,
          std::vector<double>& rValues,
          const ProcessInfo& rCurrentProcessInfo )
  {
  
    if (rVariable == DETERMINANT_F){
  
      const unsigned int& integration_points_number = mConstitutiveLawVector.size();
  
      if ( rValues.size() != integration_points_number )
        rValues.resize( integration_points_number );
  
      for ( unsigned int PointNumber = 0;  PointNumber < integration_points_number; PointNumber++ )
        {
  	rValues[PointNumber] = mDeterminantF0[PointNumber];
        }
  
    }
    else{
  
      LargeDisplacementElement::GetValueOnIntegrationPoints( rVariable, rValues, rCurrentProcessInfo );
  
    }
  
  }
  
  
  //************* STARTING - ENDING  METHODS
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::Initialize()
  {
      KRATOS_TRY
  
      LargeDisplacementElement::Initialize();
  
      SizeType integration_points_number = GetGeometry().IntegrationPointsNumber( mThisIntegrationMethod );
  
      //Resize historic deformation gradient
      if ( mDeformationGradientF0.size() != integration_points_number )
          mDeformationGradientF0.resize( integration_points_number );
  
      if ( mDeterminantF0.size() != integration_points_number )
          mDeterminantF0.resize( integration_points_number, false );
  
      for ( unsigned int PointNumber = 0; PointNumber < integration_points_number; PointNumber++ )
      {
          mDeterminantF0[PointNumber] = 1;
          mDeformationGradientF0[PointNumber] = identity_matrix<double> (3);
      }
  
  
      KRATOS_CATCH( "" )
  }
  
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::InitializeElementData (ElementDataType & rVariables, const ProcessInfo& rCurrentProcessInfo)
  {
      KRATOS_TRY
  
      const SizeType number_of_nodes  = GetGeometry().size();
      const SizeType dimension        = GetGeometry().WorkingSpaceDimension();
      const unsigned int voigt_size      = 4;
  
      rVariables.Initialize(voigt_size, dimension, number_of_nodes);
  
      rVariables.F.resize(3,3,false);
      rVariables.F = IdentityMatrix(3);
      rVariables.F0.resize(3,3,false);
      rVariables.F0 = IdentityMatrix(3);
  
      rVariables.H.resize(3,3,false);
  
      //set variables including all integration points values
  
      //reading shape functions
      rVariables.SetShapeFunctions(GetGeometry().ShapeFunctionsValues( mThisIntegrationMethod ));
  
      //reading shape functions local gradients
      rVariables.SetShapeFunctionsGradients(GetGeometry().ShapeFunctionsLocalGradients( mThisIntegrationMethod ));
  
      //set process info
      rVariables.SetProcessInfo(rCurrentProcessInfo);
  
      //calculating the current jacobian from cartesian coordinates to parent coordinates for all integration points [dx_n+1/d£]
      rVariables.j = GetGeometry().Jacobian( rVariables.j, mThisIntegrationMethod );
  
      //Calculate Delta Position
      ElementUtilities::CalculateDeltaPosition(rVariables.DeltaPosition,this->GetGeometry());
  
      //calculating the reference jacobian from cartesian coordinates to parent coordinates for all integration points [dx_n/d£]
      rVariables.J = GetGeometry().Jacobian( rVariables.J, mThisIntegrationMethod, rVariables.DeltaPosition );
  
      KRATOS_CATCH( "" )
  }
  
  ////************************************************************************************
  ////************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::FinalizeStepVariables( ElementDataType & rVariables, const double& rPointNumber )
  {
      //update internal (historical) variables
      mDeterminantF0[rPointNumber]         = rVariables.detF * rVariables.detF0;
      noalias(mDeformationGradientF0[rPointNumber]) = prod(rVariables.F, rVariables.F0);
  }
  
  
  //************* COMPUTING  METHODS
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::CalculateAndAddLHS(LocalSystemComponents& rLocalSystem, ElementDataType& rVariables, double& rIntegrationWeight)
  {
  
      double IntegrationWeight = rIntegrationWeight * 2.0 * Globals::Pi * rVariables.CurrentRadius;
      if ( this->GetProperties().Has( THICKNESS ) )
        IntegrationWeight /= GetProperties()[THICKNESS];
  
      //contributions to stiffness matrix calculated on the reference config
  
      LargeDisplacementElement::CalculateAndAddLHS( rLocalSystem, rVariables, IntegrationWeight );
  
      //KRATOS_WATCH( rLeftHandSideMatrix )
  }
  
  
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::CalculateAndAddRHS(LocalSystemComponents& rLocalSystem, ElementDataType& rVariables, Vector& rVolumeForce, double& rIntegrationWeight)
  {
      double IntegrationWeight = rIntegrationWeight * 2.0 * Globals::Pi * rVariables.CurrentRadius;
      if ( this->GetProperties().Has( THICKNESS ) )
        IntegrationWeight /= GetProperties()[THICKNESS];
  
      //contribution to external forces
  
      LargeDisplacementElement::CalculateAndAddRHS( rLocalSystem, rVariables, rVolumeForce, IntegrationWeight );
  
      //KRATOS_WATCH( rRightHandSideVector )
  }
  
  
  
  //************************************CALCULATE TOTAL MASS****************************
  //************************************************************************************
  
  double& AxisymmetricUpdatedLagrangianElement::CalculateTotalMass( double& rTotalMass, const ProcessInfo& rCurrentProcessInfo )
  {
      KRATOS_TRY
  
      //Compute the Volume Change acumulated:
      ElementDataType Variables;
      this->InitializeElementData(Variables,rCurrentProcessInfo);
  
      const GeometryType::IntegrationPointsArrayType& integration_points = GetGeometry().IntegrationPoints( mThisIntegrationMethod );
  
      rTotalMass = 0;
      //reading integration points
      for ( unsigned int PointNumber = 0; PointNumber < integration_points.size(); PointNumber++ )
        {
  	//compute element kinematics
  	this->CalculateKinematics(Variables,PointNumber);
  
  	//getting informations for integration
          Variables.IntegrationWeight = Variables.detJ * integration_points[PointNumber].Weight();
  
  	//compute point volume change
  	double PointVolumeChange = 0;
  	PointVolumeChange = this->CalculateVolumeChange( PointVolumeChange, Variables );
  
  	rTotalMass += PointVolumeChange * GetProperties()[DENSITY] * 2.0 * Globals::Pi * Variables.CurrentRadius * Variables.IntegrationWeight;
  
        }
  
      return rTotalMass;
  
      KRATOS_CATCH( "" )
  }
  
  
  
  //*********************************COMPUTE KINEMATICS*********************************
  //************************************************************************************
  
  
  void AxisymmetricUpdatedLagrangianElement::CalculateKinematics(ElementDataType& rVariables,
          const double& rPointNumber)
  
  {
      KRATOS_TRY
  
      //Get the parent coodinates derivative [dN/d£]
      const GeometryType::ShapeFunctionsGradientsType& DN_De = rVariables.GetShapeFunctionsGradients();
  
      //Get the shape functions for the order of the integration method [N]
      const Matrix& Ncontainer = rVariables.GetShapeFunctions();
  
      //Parent to reference configuration
      rVariables.StressMeasure = ConstitutiveLaw::StressMeasure_Cauchy;
  
      //Calculating the inverse of the jacobian and the parameters needed [d£/dx_n]
      Matrix InvJ;
      MathUtils<double>::InvertMatrix( rVariables.J[rPointNumber], InvJ, rVariables.detJ);
  
      //Compute cartesian derivatives [dN/dx_n]
      noalias( rVariables.DN_DX ) = prod( DN_De[rPointNumber], InvJ );
  
      //Set Shape Functions Values for this integration point
      noalias(rVariables.N) = matrix_row<const Matrix>( Ncontainer, rPointNumber);
  
      //Calculate IntegrationPoint radius
      CalculateRadius (rVariables.CurrentRadius, rVariables.ReferenceRadius, rVariables.N);
  
      //Current Deformation Gradient [dx_n+1/dx_n]
      CalculateDeformationGradient (rVariables.F, rVariables.DN_DX, rVariables.DeltaPosition, rVariables.CurrentRadius, rVariables.ReferenceRadius);
  
      //Determinant of the deformation gradient F
      rVariables.detF  = MathUtils<double>::Det(rVariables.F);
  
      //Calculating the inverse of the jacobian and the parameters needed [d£/dx_n+1]
      Matrix Invj;
      MathUtils<double>::InvertMatrix( rVariables.j[rPointNumber], Invj, rVariables.detJ); //overwrites detJ
  
      //Compute cartesian derivatives [dN/dx_n+1]
      noalias(rVariables.DN_DX) = prod( DN_De[rPointNumber], Invj ); //overwrites DX now is the current position dx
  
      //Determinant of the Deformation Gradient F0
      rVariables.detF0 = mDeterminantF0[rPointNumber];
      noalias(rVariables.F0)  = mDeformationGradientF0[rPointNumber];
  
      //Compute the deformation matrix B
      CalculateDeformationMatrix(rVariables.B, rVariables.DN_DX, rVariables.N, rVariables.CurrentRadius);
  
  
      KRATOS_CATCH( "" )
  }
  
  
  
  //*************************COMPUTE AXYSIMMETRIC RADIUS********************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::CalculateRadius(double & rCurrentRadius,
  							   double & rReferenceRadius,
  							   const Vector& rN)
  
  
  {
  
      KRATOS_TRY
  
      const SizeType number_of_nodes  = GetGeometry().PointsNumber();
  
      const SizeType dimension  = GetGeometry().WorkingSpaceDimension();
  
      rCurrentRadius=0;
      rReferenceRadius=0;
  
      if ( dimension == 2 )
      {
          for ( SizeType i = 0; i < number_of_nodes; i++ )
          {
              //Displacement from the reference to the current configuration
              // array_1d<double, 3 > & CurrentDisplacement  = GetGeometry()[i].FastGetSolutionStepValue(DISPLACEMENT);
              // array_1d<double, 3 > & PreviousDisplacement = GetGeometry()[i].FastGetSolutionStepValue(DISPLACEMENT,1);
              // array_1d<double, 3 > DeltaDisplacement      = CurrentDisplacement-PreviousDisplacement;
  	    // array_1d<double, 3 > & CurrentPosition      = GetGeometry()[i].Coordinates();
  	    // array_1d<double, 3 > ReferencePosition      = CurrentPosition - DeltaDisplacement;
  
              // rCurrentRadius   += CurrentPosition[0]*rN[i];
              // rReferenceRadius += ReferencePosition[0]*rN[i];
              //std::cout<<" node "<<i<<" -> DeltaDisplacement : "<<DeltaDisplacement<<std::endl;
  
  	    rCurrentRadius   += rN[i] * GetGeometry()[i].X();
  	    rReferenceRadius += rN[i] * (GetGeometry()[i].X() + GetGeometry()[i].FastGetSolutionStepValue(DISPLACEMENT_X,1) - GetGeometry()[i].FastGetSolutionStepValue(DISPLACEMENT_X));
  
          }
      }
  
  
      if ( dimension == 3 )
      {
          std::cout<<" AXISYMMETRIC case and 3D is not possible "<<std::endl;
      }
  
      KRATOS_CATCH( "" )
  }
  
  //*************************COMPUTE DEFORMATION GRADIENT*******************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::CalculateDeformationGradient(Matrix&  rF,
                                                                          const Matrix& rDN_DX,
                                                                          const Matrix& rDeltaPosition,
                                                                          const double & rCurrentRadius,
                                                                          const double & rReferenceRadius)
  {
      KRATOS_TRY
  
      const SizeType number_of_nodes  = GetGeometry().PointsNumber();
      const SizeType dimension        = GetGeometry().WorkingSpaceDimension();
  
      rF = identity_matrix<double> ( 3 );
  
      if( dimension == 2 )
      {
  
          for ( SizeType i = 0; i < number_of_nodes; i++ )
          {
              rF ( 0 , 0 ) += rDeltaPosition(i,0)*rDN_DX ( i , 0 );
              rF ( 0 , 1 ) += rDeltaPosition(i,0)*rDN_DX ( i , 1 );
              rF ( 1 , 0 ) += rDeltaPosition(i,1)*rDN_DX ( i , 0 );
              rF ( 1 , 1 ) += rDeltaPosition(i,1)*rDN_DX ( i , 1 );
          }
  
          rF ( 2 , 2 ) = rCurrentRadius/rReferenceRadius;
      }
      else if( dimension == 3)
      {
  
          std::cout<<" AXISYMMETRIC case and 3D is not possible "<<std::endl;
      }
      else
      {
  
          KRATOS_THROW_ERROR( std::invalid_argument, "something is wrong with the dimension", "" );
  
      }
  
      KRATOS_CATCH( "" )
  }
  
  
  
  
  //************************************************************************************
  //************************************************************************************
  
  
  void AxisymmetricUpdatedLagrangianElement::CalculateDeformationMatrix(Matrix& rB,
                                                                        const Matrix& rDN_DX,
                                                                        const Vector& rN,
                                                                        const double & rCurrentRadius)
  {
      KRATOS_TRY
  
      const SizeType number_of_nodes  = GetGeometry().PointsNumber();
      const SizeType dimension        = GetGeometry().WorkingSpaceDimension();
  
      rB.clear(); //set all components to zero
  
      if( dimension == 2 )
      {
  
          for ( SizeType i = 0; i < number_of_nodes; i++ )
          {
              unsigned int index = 2 * i;
  
              rB( 0, index + 0 ) = rDN_DX( i, 0 );
              rB( 1, index + 1 ) = rDN_DX( i, 1 );
              rB( 2, index + 0 ) = rN[i]/rCurrentRadius;
              rB( 3, index + 0 ) = rDN_DX( i, 1 );
              rB( 3, index + 1 ) = rDN_DX( i, 0 );
  
          }
  
      }
      else if( dimension == 3 )
      {
  
          std::cout<<" AXISYMMETRIC case and 3D is not possible "<<std::endl;
  
      }
      else
      {
  
          KRATOS_THROW_ERROR( std::invalid_argument, "something is wrong with the dimension", "" );
  
      }
  
      KRATOS_CATCH( "" )
  }
  
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::CalculateGreenLagrangeStrain(const Matrix& rF,
          Vector& rStrainVector )
  {
      KRATOS_TRY
  
      const SizeType dimension   = GetGeometry().WorkingSpaceDimension();
  
      //Right Cauchy-Green Calculation
      Matrix C ( 3, 3 );
      noalias ( C )  = prod( trans( rF ), rF );
  
      if( dimension == 2 )
      {
  
          //Green Lagrange Strain Calculation
          if ( rStrainVector.size() != 4 ) rStrainVector.resize( 4, false );
  
          rStrainVector[0] = 0.5 * ( C( 0, 0 ) - 1.00 );
  
          rStrainVector[1] = 0.5 * ( C( 1, 1 ) - 1.00 );
  
          rStrainVector[2] = 0.5 * ( C( 2, 2 ) - 1.00 );
  
          rStrainVector[3] = C( 0, 1 ); // xy
  
      }
      else if( dimension == 3 )
      {
  
          std::cout<<" AXISYMMETRIC case and 3D is not possible "<<std::endl;
  
      }
      else
      {
  
          KRATOS_THROW_ERROR( std::invalid_argument, "something is wrong with the dimension", "" );
  
      }
  
      KRATOS_CATCH( "" )
  }
  
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::CalculateAlmansiStrain(const Matrix& rF,
          Vector& rStrainVector )
  {
      KRATOS_TRY
  
      const SizeType dimension  = GetGeometry().WorkingSpaceDimension();
  
      //Left Cauchy-Green Calculation
      Matrix LeftCauchyGreen(rF.size1(),rF.size1());
      noalias(LeftCauchyGreen) = prod( rF, trans( rF ) );
  
      //Calculating the inverse of the jacobian
      Matrix InverseLeftCauchyGreen(rF.size1(),rF.size1());
      double det_b=0;
      MathUtils<double>::InvertMatrix( LeftCauchyGreen, InverseLeftCauchyGreen, det_b);
  
      if( dimension == 2 )
      {
  
          //Almansi Strain Calculation
          if ( rStrainVector.size() != 4 ) rStrainVector.resize( 4, false );
  
          rStrainVector[0] = 0.5 * (  1.00 - InverseLeftCauchyGreen( 0, 0 ) );
  
          rStrainVector[1] = 0.5 * (  1.00 - InverseLeftCauchyGreen( 1, 1 ) );
  
          rStrainVector[2] = 0.5 * ( 1.00 - InverseLeftCauchyGreen( 2, 2 ) );
  
          rStrainVector[3] = - InverseLeftCauchyGreen( 0, 1 ); // xy
  
      }
      else if( dimension == 3 )
      {
  
          std::cout<<" AXISYMMETRIC case and 3D is not possible "<<std::endl;
  
      }
      else
      {
  
          KRATOS_THROW_ERROR( std::invalid_argument, "something is wrong with the dimension", "" )
  
      }
  
  
      KRATOS_CATCH( "" )
  }
  
  
  
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::CalculateAndAddKuug(MatrixType& rK,
          ElementDataType& rVariables,
          double& rIntegrationWeight)
  
  {
      KRATOS_TRY
  
      const SizeType number_of_nodes  = GetGeometry().size();
  
      //Matrix Kh = rK;
  
      // axisymmetric geometric matrix
  
      double alpha1 = 0;
      double alpha2 = 0;
      double alpha3 = 0;
  
      unsigned int indexi = 0;
      unsigned int indexj = 0;
  
      for ( SizeType i = 0; i < number_of_nodes; i++ )
      {
          indexj =0;
          for ( SizeType j = 0; j < number_of_nodes; j++ )
          {
              alpha1 = rVariables.DN_DX(j,0) * ( rVariables.DN_DX(i,0) * rVariables.StressVector[0] + rVariables.DN_DX(i,1) * rVariables.StressVector[3] );
              alpha2 = rVariables.DN_DX(j,1) * ( rVariables.DN_DX(i,0) * rVariables.StressVector[3] + rVariables.DN_DX(i,1) * rVariables.StressVector[1] );
              alpha3 = rVariables.N[i] * rVariables.N[j] * rVariables.StressVector[2] * (1.0/rVariables.CurrentRadius*rVariables.CurrentRadius);
  
              rK(indexi,indexj)     += (alpha1 + alpha2 + alpha3) * rIntegrationWeight ;
              rK(indexi+1,indexj+1) += (alpha1 + alpha2) * rIntegrationWeight ;
  
              indexj+=2;
          }
  
          indexi+=2;
  
      }
  
      //std::cout<<std::endl;
      //std::cout<<" Kgeo "<<rK-Kh<<std::endl;
  
      KRATOS_CATCH( "" )
  }
  
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::GetHistoricalVariables( ElementDataType& rVariables, const double& rPointNumber )
  {
      LargeDisplacementElement::GetHistoricalVariables(rVariables,rPointNumber);
  
      //Deformation Gradient F0
      rVariables.detF0 = mDeterminantF0[rPointNumber];
      rVariables.F0    = mDeformationGradientF0[rPointNumber];
  
      rVariables.CurrentRadius = rVariables.ReferenceRadius;
  }
  
  
  //************************************CALCULATE VOLUME CHANGE*************************
  //************************************************************************************
  
  double& AxisymmetricUpdatedLagrangianElement::CalculateVolumeChange( double& rVolumeChange, ElementDataType& rVariables )
  {
      KRATOS_TRY
  
      rVolumeChange = 1.0 / (rVariables.detF * rVariables.detF0);
  
      return rVolumeChange;
  
      KRATOS_CATCH( "" )
  }
  
  
  //************************************************************************************
  //************************************************************************************
  
  int AxisymmetricUpdatedLagrangianElement::Check( const ProcessInfo& rCurrentProcessInfo )
  {
      KRATOS_TRY
  
      // Perform base element checks
      int ErrorCode = 0;
      ErrorCode = LargeDisplacementElement::Check(rCurrentProcessInfo);
  
      // Check that the element nodes contain all required SolutionStepData and Degrees of freedom
      for(SizeType i=0; i<this->GetGeometry().size(); ++i)
        {
  	// Nodal data
  	Node<3> &rNode = this->GetGeometry()[i];
  	KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(DISPLACEMENT,rNode);
  	//KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(VOLUME_ACCELERATION,rNode);
  
  	// Nodal dofs
  	KRATOS_CHECK_DOF_IN_NODE(DISPLACEMENT_X,rNode);
  	KRATOS_CHECK_DOF_IN_NODE(DISPLACEMENT_Y,rNode);
  	if( rCurrentProcessInfo[SPACE_DIMENSION] == 3)
  	  KRATOS_CHECK_DOF_IN_NODE(DISPLACEMENT_Z,rNode);
        }
  
      return ErrorCode;
  
      KRATOS_CATCH( "" );
  }
  
  //************************************************************************************
  //************************************************************************************
  
  void AxisymmetricUpdatedLagrangianElement::save( Serializer& rSerializer ) const
  {
      KRATOS_SERIALIZE_SAVE_BASE_CLASS( rSerializer, LargeDisplacementElement )
      rSerializer.save("DeformationGradientF0",mDeformationGradientF0);
      rSerializer.save("DeterminantF0",mDeterminantF0);
  }
  
  void AxisymmetricUpdatedLagrangianElement::load( Serializer& rSerializer )
  {
      KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, LargeDisplacementElement )
      rSerializer.load("DeformationGradientF0",mDeformationGradientF0);
      rSerializer.load("DeterminantF0",mDeterminantF0);
  }
  
  
  } // Namespace Kratos
  