contact_domain_condition.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosContactMechanicsApplication $
//   Created by:          $Author:              JMCarbonell $
//   Last modified by:    $Co-Author:                       $
//   Date:                $Date:                  July 2016 $
//   Revision:            $Revision:                    0.0 $
//
//
 
#if !defined(KRATOS_CONTACT_DOMAIN_CONDITION_H_INCLUDED )
#define  KRATOS_CONTACT_DOMAIN_CONDITION_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "includes/constitutive_law.h"
#include "includes/checks.h"
#include "includes/element.h"
#include "includes/condition.h"
#include "includes/kratos_flags.h"
 
#include "custom_utilities/contact_domain_utilities.hpp"
#include "custom_utilities/mesh_data_transfer_utilities.hpp"
 
namespace Kratos
{
 
 
 
 
class KRATOS_API(CONTACT_MECHANICS_APPLICATION) ContactDomainCondition
    : public Condition
{
public:
 
 
  typedef ConstitutiveLaw ConstitutiveLawType;
  typedef ConstitutiveLawType::Pointer ConstitutiveLawPointerType;
  typedef GeometryData::IntegrationMethod IntegrationMethod;
 
  typedef Node NodeType;
  typedef Geometry<NodeType> GeometryType;
  typedef Element::ElementType ElementType;
 
 
  typedef ContactDomainUtilities::PointType             PointType;
  typedef ContactDomainUtilities::SurfaceVector     SurfaceVector;
  typedef ContactDomainUtilities::SurfaceScalar     SurfaceScalar;
  typedef ContactDomainUtilities::BaseLengths         BaseLengths;
 
  typedef GlobalPointersVector<Node >       NodeWeakPtrVectorType;
  typedef GlobalPointersVector<Element>     ElementWeakPtrVectorType;
  typedef GlobalPointersVector<Condition> ConditionWeakPtrVectorType;
 
  typedef ContactDomainUtilities::ScalarBaseType   ScalarBaseType;
  typedef ContactDomainUtilities::SurfaceBase         SurfaceBase;
 
protected:
 
  typedef struct
  {
    //Geometrical surface tangent gaps:
    ScalarBaseType   CurrentGap;     //tangential gap
 
    //Contact constraint parameters
    double   Multiplier;            //Lagrange Multipliyer tangent
    double   Penalty;               //Penalty Parameter tangent
 
    //Variables of the contact domain elements
    Vector          dN_dt;      //Discrete variacion of the shape function  in the current tangent direction
    std::vector<Vector >       Tsigma;
 
    ScalarBaseType   CurrentTensil;
    double           GapSign;
 
  } ContactSurfaceParameters;
 
 
 
  typedef struct
  {
    ContactSurfaceParameters A;
    ContactSurfaceParameters B;
 
    SurfaceBase    CovariantBase;
    SurfaceBase    ContravariantBase;
 
    //geometrical variables
    double ReferenceArea;
    double CurrentArea;
 
    double FactorArea;
 
    double EquivalentHeigh;
 
    double ElementSize;
 
  } ContactTangentParameters;
 
 
  typedef struct
  {
    Flags           Options;               //calculation options
 
    //Geometrical gaps:
    SurfaceScalar   CurrentGap;             //normal and tangential gap
 
    //The stabilization factor or penalty factor
    SurfaceScalar   ContactFactor;
 
    //Friction:
    double          FrictionCoefficient;   //total friction coeffitient mu
    double          TangentialGapSign;     //sign or direction of the tangential gap
 
    SurfaceScalar   Multiplier;            //Lagrange Multipliyer normal and tangent
    SurfaceScalar   Penalty;               //Penalty Parameter normal and tangent
 
 
    //variables of the contact element 2D
    Vector          dN_dn;      //Discrete variacion of the shape function  in the current normal direction
    Vector          dN_dt;      //Discrete variacion of the shape function  in the current tangent direction
    Vector          dN_drn;     //Discrete variacion of the shape function  in the reference normal direction
 
    std::vector<Vector >       Nsigma;
    std::vector<Vector >       Tsigma;
 
    //Geometric variables
    SurfaceVector        CurrentSurface;
 
    std::vector<BaseLengths>   CurrentBase;    //Current Base Lengths variables
    std::vector<BaseLengths>   ReferenceBase;  //Reference Base Lengths variables
 
    //Resultant mechanical tractions
    SurfaceScalar              CurrentTensil;  //Tangential and Normal modulus of the traction vector components
 
 
    //Tangent parameters used in 3D
    ContactTangentParameters Tangent;
 
 
  } ContactParameters;
 
 
 
  typedef struct
  {
    ScalarBaseType   PreviousGapA;
    ScalarBaseType   PreviousGapB;
 
    SurfaceBase    CovariantBase;
    SurfaceBase    ContravariantBase;
 
  } ContactTangentVariables;
 
 
  struct ConditionVariables
  {
    double  detF;
    double  detJ;
    Vector  StrainVector;
    Vector  StressVector;
    Vector  N;
    Matrix  F;
    Matrix  DN_DX;
    Matrix  ConstitutiveMatrix;
 
    ContactParameters Contact;
 
    //Axisymmetric
    double  CurrentRadius;
    double  ReferenceRadius;
 
    //variables including all integration points
    const GeometryType::ShapeFunctionsGradientsType* pDN_De;
    const Matrix* pNcontainer;
 
    GeometryType::JacobiansType J;
    GeometryType::JacobiansType j;
    Matrix  DeltaPosition;
 
    void SetShapeFunctionsGradients(const GeometryType::ShapeFunctionsGradientsType &rDN_De)
    {
      pDN_De=&rDN_De;
    };
 
    void SetShapeFunctions(const Matrix& rNcontainer)
    {
      pNcontainer=&rNcontainer;
    };
 
 
    const GeometryType::ShapeFunctionsGradientsType& GetShapeFunctionsGradients()
    {
      return *pDN_De;
    };
 
    const Matrix& GetShapeFunctions()
    {
      return *pNcontainer;
    };
 
  };
 
 
  struct ContactVariables
  {
 
    //Iteration counter:
    int             IterationCounter;     //the number of the step iteration
 
    //The stabilization parameter and penalty parameter
    double          StabilizationFactor;
    double          PenaltyFactor;
 
    //Geometrical gaps:
    SurfaceScalar        PreviousGap;     //effective normal and tangential gap in previous time step configuration
 
    //Geometric variables
    SurfaceVector        PreStepSurface;
    SurfaceVector        ReferenceSurface;
 
    //Tangent variables in 3D
    ContactTangentVariables  Tangent;
 
    //Contact condition conectivities
    std::vector<unsigned int> nodes;
    std::vector<unsigned int> order;
    std::vector<unsigned int> slaves;
 
    //Resultant mechanical tractions
    PointType       TractionVector;       //Traction Vector in the reference configuration
 
    //Pointer Variables
    GeometryType*         mpMasterGeometry;
    ElementType*          mpMasterElement;
    Condition*            mpMasterCondition;
    NodeType*             mpMasterNode;
 
    void SetMasterGeometry  (GeometryType& rGeometry){ mpMasterGeometry = &rGeometry; }
    void SetMasterElement   (ElementType& rElement){ mpMasterElement = &rElement; }
    void SetMasterCondition (ConditionType& rCondition){ mpMasterCondition = &rCondition; }
    void SetMasterNode      (NodeType& rNode){ mpMasterNode = &rNode; }
 
    GeometryType& GetMasterGeometry()   { return (*mpMasterGeometry); }
    ElementType& GetMasterElement()     { return (*mpMasterElement); }
    ConditionType& GetMasterCondition() { return (*mpMasterCondition); }
    NodeType& GetMasterNode()           { return (*mpMasterNode); }
 
 
  };
 
 
 
 
  struct LocalSystemComponents
  {
   private:
 
    //for calculation local system with compacted LHS and RHS
    MatrixType *mpLeftHandSideMatrix;
    VectorType *mpRightHandSideVector;
 
    //for calculation local system with LHS and RHS components
    std::vector<MatrixType> *mpLeftHandSideMatrices;
    std::vector<VectorType> *mpRightHandSideVectors;
 
    //LHS variable components
    const std::vector< Variable< MatrixType > > *mpLeftHandSideVariables;
 
    //RHS variable components
    const std::vector< Variable< VectorType > > *mpRightHandSideVariables;
 
 
   public:
 
    //calculation flags
    Flags  CalculationFlags;
 
    void SetLeftHandSideMatrix( MatrixType& rLeftHandSideMatrix ) { mpLeftHandSideMatrix = &rLeftHandSideMatrix; };
    void SetLeftHandSideMatrices( std::vector<MatrixType>& rLeftHandSideMatrices ) { mpLeftHandSideMatrices = &rLeftHandSideMatrices; };
    void SetLeftHandSideVariables(const std::vector< Variable< MatrixType > >& rLeftHandSideVariables ) { mpLeftHandSideVariables = &rLeftHandSideVariables; };
 
    void SetRightHandSideVector( VectorType& rRightHandSideVector ) { mpRightHandSideVector = &rRightHandSideVector; };
    void SetRightHandSideVectors( std::vector<VectorType>& rRightHandSideVectors ) { mpRightHandSideVectors = &rRightHandSideVectors; };
    void SetRightHandSideVariables(const std::vector< Variable< VectorType > >& rRightHandSideVariables ) { mpRightHandSideVariables = &rRightHandSideVariables; };
 
 
    MatrixType& GetLeftHandSideMatrix() { return *mpLeftHandSideMatrix; };
    std::vector<MatrixType>& GetLeftHandSideMatrices() { return *mpLeftHandSideMatrices; };
    const std::vector< Variable< MatrixType > >& GetLeftHandSideVariables() { return *mpLeftHandSideVariables; };
 
    VectorType& GetRightHandSideVector() { return *mpRightHandSideVector; };
    std::vector<VectorType>& GetRightHandSideVectors() { return *mpRightHandSideVectors; };
    const std::vector< Variable< VectorType > >& GetRightHandSideVariables() { return *mpRightHandSideVariables; };
 
  };
 
 
public:
 
  KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION( ContactDomainCondition );
 
 
 
 
  ContactDomainCondition() : Condition() {};
 
  ContactDomainCondition(IndexType NewId, GeometryType::Pointer pGeometry);
 
  ContactDomainCondition(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties);
 
  ContactDomainCondition(ContactDomainCondition const& rOther);
 
 
  virtual ~ContactDomainCondition();
 
 
  ContactDomainCondition& operator=(ContactDomainCondition const& rOther);
 
 
 
  Condition::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes, PropertiesType::Pointer pProperties) const override;
 
  Condition::Pointer Clone(IndexType NewId,
                           NodesArrayType const& ThisNodes) const override;
 
  //************* GETTING METHODS
 
  IntegrationMethod GetIntegrationMethod() const override;
 
  void GetDofList(DofsVectorType& rConditionalDofList, const ProcessInfo& rCurrentProcessInfo) const override;
 
  void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo& rCurrentProcessInfo) const override;
 
  void GetValuesVector(Vector& rValues, int Step = 0) const override;
 
  void GetFirstDerivativesVector(Vector& rValues, int Step = 0) const override;
 
  void GetSecondDerivativesVector(Vector& rValues, int Step = 0) const override;
 
 
 
  //on integration points:
  //SET
  void SetValuesOnIntegrationPoints(const Variable<double>& rVariable, const std::vector<double>& rValues, const ProcessInfo& rCurrentProcessInfo) override;
  void SetValuesOnIntegrationPoints(const Variable<Vector>& rVariable, const std::vector<Vector>& rValues, const ProcessInfo& rCurrentProcessInfo) override;
 
  void SetValuesOnIntegrationPoints(const Variable<Matrix>& rVariable, const std::vector<Matrix>& rValues, const ProcessInfo& rCurrentProcessInfo) override;
 
 
 
  //************* STARTING - ENDING  METHODS
 
  void Initialize(const ProcessInfo& CurrentProcessInfo) override;
 
 
  void InitializeSolutionStep(const ProcessInfo& CurrentProcessInfo) override;
 
  void InitializeNonLinearIteration(const ProcessInfo& CurrentProcessInfo) override;
 
 
  void FinalizeNonLinearIteration(const ProcessInfo& CurrentProcessInfo) override;
 
 
  void FinalizeSolutionStep(const ProcessInfo& CurrentProcessInfo) override;
 
 
 
  //************* COMPUTING  METHODS
 
  void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector,const ProcessInfo& rCurrentProcessInfo) override;
 
  void CalculateLocalSystem(std::vector< MatrixType >& rLeftHandSideMatrices,
                            const std::vector< Variable< MatrixType > >& rLHSVariables,
                            std::vector< VectorType >& rRightHandSideVectors,
                            const std::vector< Variable< VectorType > >& rRHSVariables,
                            const ProcessInfo& rCurrentProcessInfo);
 
  void CalculateRightHandSide(VectorType& rRightHandSideVector,const ProcessInfo& rCurrentProcessInfo) override;
 
 
  void CalculateRightHandSide(std::vector< VectorType >& rRightHandSideVectors,
                              const std::vector< Variable< VectorType > >& rRHSVariables,
                              const ProcessInfo& rCurrentProcessInfo);
 
  void CalculateLeftHandSide (MatrixType& rLeftHandSideMatrix, const ProcessInfo& rCurrentProcessInfo) override;
 
 
  virtual void AddExplicitContribution(const VectorType& rRHSVector,
                                       const Variable<VectorType>& rRHSVariable,
                                       const Variable<array_1d<double,3> >& rDestinationVariable,
                                       const ProcessInfo& rCurrentProcessInfo) override;
 
 
  //on integration points:
  void CalculateOnIntegrationPoints(const Variable<double>& rVariable, std::vector<double>& rOutput, const ProcessInfo& rCurrentProcessInfo) override;
 
  void CalculateOnIntegrationPoints(const Variable<Vector>& rVariable, std::vector<Vector>& rOutput, const ProcessInfo& rCurrentProcessInfo) override;
 
  void CalculateOnIntegrationPoints(const Variable<Matrix >& rVariable, std::vector< Matrix >& rOutput, const ProcessInfo& rCurrentProcessInfo) override;
 
 
  //************************************************************************************
  //************************************************************************************
  int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
  //std::string Info() const;
 
 
 
  //      virtual String Info() const;
 
  //      virtual void PrintInfo(std::ostream& rOStream) const;
 
  //      virtual void PrintData(std::ostream& rOStream) const;
 
protected:
 
  IntegrationMethod mThisIntegrationMethod;
 
  ContactVariables      mContactVariables;
 
  ContactDomainUtilities  mContactUtilities;
 
 
  virtual void SetMasterGeometry()
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
  };
 
 
  virtual void CalculateContactFactor(const ProcessInfo& rCurrentProcessInfo)
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
  };
 
  virtual void CalculatePreviousGap()
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
 
    };
 
 
  virtual void InitializeConditionVariables (ConditionVariables& rVariables,
                         const ProcessInfo& rCurrentProcessInfo);
 
  virtual void CalculateConditionSystem(LocalSystemComponents& rLocalSystem,
                                        const ProcessInfo& rCurrentProcessInfo);
  void ClearNodalForces ();
 
 
  void CalculateNodalForces (const ProcessInfo& CurrentProcessInfo);
 
 
  void ClearMasterElementNodalForces(ElementType& rMasterElement);
 
 
  void SetContactIntegrationVariable (Vector & rContactVariable,
                                      std::vector<Vector> & rMasterVariables,
                                      const unsigned int& rPointNumber);
 
  virtual void CalculateKinematics(ConditionVariables& rVariables,
                                   const ProcessInfo& rCurrentProcessInfo,
                                   const unsigned int& rPointNumber);
 
  Matrix& CalculateDeltaPosition(Matrix & rDeltaPosition);
 
  virtual void CalculateExplicitFactors(ConditionVariables& rVariables,
                                        const ProcessInfo& rCurrentProcessInfo)
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
 
    };
  virtual void CalculateDomainShapeN(ConditionVariables& rVariables)
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
 
    };
 
 
 
  virtual void CalculateRelativeVelocity(ConditionVariables& rVariables,
                                         PointType & TangentVelocity,
                                         const ProcessInfo& rCurrentProcessInfo);
 
  virtual void CalculateRelativeDisplacement(ConditionVariables& rVariables,
                                             PointType & TangentDisplacement,
                                             const ProcessInfo& rCurrentProcessInfo);
 
  virtual PointType & CalculateCurrentTangent(PointType &rTangent)
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
 
    };
 
 
  virtual void CalculateFrictionCoefficient(ConditionVariables& rVariables,
                                            const PointType & TangentVelocity);
 
 
  virtual double& CalculateIntegrationWeight(double& rIntegrationWeight);
 
 
  virtual void InitializeSystemMatrices(MatrixType& rLeftHandSideMatrix,
                                        VectorType& rRightHandSideVector,
                                        Flags& rCalculationFlags);
 
 
  void CalculatePerturbedLeftHandSide (MatrixType& rLeftHandSideMatrix,
                                       const ProcessInfo& rCurrentProcessInfo);
 
  virtual void CalculateAndAddLHS(LocalSystemComponents& rLocalSystem,
                                  ConditionVariables& rVariables,
                                  double& rIntegrationWeight);
 
  virtual void CalculateAndAddRHS(LocalSystemComponents& rLocalSystem,
                                  ConditionVariables& rVariables,
                                  double& rIntegrationWeight);
 
 
  virtual void CalculateAndAddKuug(MatrixType& rLeftHandSideMatrix,
                                   ConditionVariables& rVariables,
                                   double& rIntegrationWeight);
 
  virtual void CalculateContactStiffness (double &Kcont,ConditionVariables& rVariables,
                                          unsigned int& ndi,unsigned int& ndj,
                                          unsigned int& idir,unsigned int& jdir)
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
 
    };
 
  virtual void CalculateAndAddContactForces(VectorType& rRightHandSideVector,
                                            ConditionVariables& rVariables,
                                            double& rIntegrationWeight);
 
 
 
  virtual void CalculateNormalForce       (double &F,ConditionVariables& rVariables,
                                           unsigned int& ndi,unsigned int& idir)
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
  };
  virtual void CalculateTangentStickForce (double &F,ConditionVariables& rVariables,
                                           unsigned int& ndi,unsigned int& idir)
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
  };
 
  virtual void CalculateTangentSlipForce  (double &F,ConditionVariables& rVariables,
                                           unsigned int& ndi,unsigned int& idir)
  {
    KRATOS_THROW_ERROR( std::invalid_argument, "Calling base class in contact domain", "" )
  };
 
 
  friend class Serializer;
 
  void save(Serializer& rSerializer) const override;
 
  void load(Serializer& rSerializer) override;
 
 
private:
 
 
}; // Class ContactDomainCondition
 
 
/*  inline std::istream& operator >> (std::istream& rIStream,
    ContactDomainCondition& rThis);
*/
/*  inline std::ostream& operator << (std::ostream& rOStream,
    const ContactDomainCondition& rThis)
    {
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
    }*/
 
} // namespace Kratos.
#endif // KRATOS_CONTACT_DOMAIN_CONDITION_H_INCLUDED  defined