line_search_contact_strategy.h

Go to the documentation of this file.

// KRATOS    ______            __             __  _____ __                  __                   __
//          / ____/___  ____  / /_____ ______/ /_/ ___// /________  _______/ /___  ___________ _/ /
//         / /   / __ \/ __ \/ __/ __ `/ ___/ __/\__ \/ __/ ___/ / / / ___/ __/ / / / ___/ __ `/ /
//        / /___/ /_/ / / / / /_/ /_/ / /__/ /_ ___/ / /_/ /  / /_/ / /__/ /_/ /_/ / /  / /_/ / /
//        \____/\____/_/ /_/\__/\__,_/\___/\__//____/\__/_/   \__,_/\___/\__/\__,_/_/   \__,_/_/  MECHANICS
//
//  License:         BSD License
//                   license: ContactStructuralMechanicsApplication/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
 
#pragma once
 
// System includes
 
// External includes
 
// Project includes
#include "includes/kratos_parameters.h"
#include "includes/define.h"
#include "includes/model_part.h"
#include "includes/variables.h"
#include "solving_strategies/strategies/implicit_solving_strategy.h"
#include "solving_strategies/strategies/line_search_strategy.h"
#include "utilities/parallel_utilities.h"
#include "utilities/variable_utils.h"
#include "utilities/atomic_utilities.h"
 
/* Convergence criterias */
#include "solving_strategies/convergencecriterias/convergence_criteria.h"
 
/* Default builder and solver */
#include "solving_strategies/builder_and_solvers/residualbased_block_builder_and_solver.h"
 
namespace Kratos
{
 
 
 
 
 
 
template<class TSparseSpace,
         class TDenseSpace, // = DenseSpace<double>,
         class TLinearSolver //= LinearSolver<TSparseSpace,TDenseSpace>
         >
class LineSearchContactStrategy :
    public LineSearchStrategy< TSparseSpace, TDenseSpace, TLinearSolver >
{
public:
    KRATOS_CLASS_POINTER_DEFINITION( LineSearchContactStrategy );
 
    using TConvergenceCriteriaType = ConvergenceCriteria<TSparseSpace, TDenseSpace>;
 
    using SolvingStrategyType = SolvingStrategy<TSparseSpace, TDenseSpace>;
 
    using StrategyBaseType = ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver>;
 
    using NRBaseType = ResidualBasedNewtonRaphsonStrategy<TSparseSpace, TDenseSpace, TLinearSolver>;
 
    using BaseType = LineSearchStrategy<TSparseSpace, TDenseSpace, TLinearSolver>;
 
    using ClassType = LineSearchContactStrategy<TSparseSpace, TDenseSpace, TLinearSolver>;
 
    using TBuilderAndSolverType = typename BaseType::TBuilderAndSolverType;
 
    using TDataType = typename BaseType::TDataType;
 
    using SparseSpaceType = TSparseSpace;
 
    using TSchemeType = typename BaseType::TSchemeType;
 
    using DofsArrayType = typename BaseType::DofsArrayType;
 
    using TSystemMatrixType = typename BaseType::TSystemMatrixType;
 
    using TSystemVectorType = typename BaseType::TSystemVectorType;
 
    using LocalSystemVectorType = typename BaseType::LocalSystemVectorType;
 
    using LocalSystemMatrixType = typename BaseType::LocalSystemMatrixType;
 
    using TSystemMatrixPointerType = typename BaseType::TSystemMatrixPointerType;
 
    using TSystemVectorPointerType = typename BaseType::TSystemVectorPointerType;
 
    using NodesArrayType = typename ModelPart::NodesContainerType;
 
    using ConditionsArrayType = typename ModelPart::ConditionsContainerType;
 
    using IndexType = std::size_t;
 
    explicit LineSearchContactStrategy()
    {
    }
 
    explicit LineSearchContactStrategy(ModelPart& rModelPart, Parameters ThisParameters)
        : BaseType(rModelPart, BaseType::GetDefaultParameters())
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
    }
 
    LineSearchContactStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TLinearSolver::Pointer pNewLinearSolver,
        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
        IndexType MaxIterations = 30,
        bool CalculateReactions = false,
        bool ReformDofSetAtEachStep = false,
        bool MoveMeshFlag = false,
        Parameters ThisParameters =  Parameters(R"({})")
    )
        : BaseType(rModelPart, pScheme, pNewLinearSolver, pNewConvergenceCriteria, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag)
    {
        KRATOS_TRY;
 
        Parameters default_parameters = this->GetDefaultParameters();
 
        ThisParameters.ValidateAndAssignDefaults(default_parameters);
 
        KRATOS_CATCH("");
    }
 
    LineSearchContactStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TLinearSolver::Pointer pNewLinearSolver,
        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
        typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver,
        IndexType MaxIterations = 30,
        bool CalculateReactions = false,
        bool ReformDofSetAtEachStep = false,
        bool MoveMeshFlag = false,
        Parameters ThisParameters =  Parameters(R"({})")
    )
        : BaseType(rModelPart, pScheme, pNewLinearSolver, pNewConvergenceCriteria, pNewBuilderAndSolver, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag )
    {
        KRATOS_TRY;
 
        Parameters default_parameters = this->GetDefaultParameters();
 
        ThisParameters.ValidateAndAssignDefaults(default_parameters);
 
        KRATOS_CATCH("");
    }
 
    ~LineSearchContactStrategy() override
    = default;
 
 
 
    typename SolvingStrategyType::Pointer Create(
        ModelPart& rModelPart,
        Parameters ThisParameters
        ) const override
    {
        return Kratos::make_shared<ClassType>(rModelPart, ThisParameters);
    }
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name" : "line_search_contact_strategy"
        })" );
 
        // Getting base class default parameters
        const Parameters base_default_parameters = BaseType::GetDefaultParameters();
        default_parameters.RecursivelyAddMissingParameters(base_default_parameters);
        return default_parameters;
    }
 
    static std::string Name()
    {
        return "line_search_contact_strategy";
    }
 
 
 
 
    std::string Info() const override
    {
        return "LineSearchContactStrategy";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
protected:
 
 
 
    bool mRecalculateFactor;         // To check if we recalculate or not the scale factor
 
 
 
    void InitializeSolutionStep() override
    {
        BaseType::InitializeSolutionStep();
 
        // TODO: Add something if necessary
    }
 
    void UpdateDatabase(
        TSystemMatrixType& A,
        TSystemVectorType& Dx,
        TSystemVectorType& b,
        const bool MoveMesh
    ) override
    {
        typename TSchemeType::Pointer pScheme = this->GetScheme();
        typename TBuilderAndSolverType::Pointer pBuilderAndSolver = this->GetBuilderAndSolver(); // FIXME: Separate in the parts of LM and displacement
 
        TSystemVectorType aux(b.size()); //TODO: do it by using the space
        TSparseSpace::Assign(aux, 0.5, Dx);
 
        TSystemVectorType DxDisp(b.size());
        TSystemVectorType DxLM(b.size());
        ComputeSplitDx(Dx, DxDisp, DxLM);
 
        // Compute residual without update
        TSparseSpace::SetToZero(b);
        pBuilderAndSolver->BuildRHS(pScheme, BaseType::GetModelPart(), b );
        double roDisp;
        double roLM;
        ComputeMixedResidual(b, roDisp, roLM);
 
        // Compute half step residual
        NRBaseType::UpdateDatabase(A,aux,b,MoveMesh);
        TSparseSpace::SetToZero(b);
        pBuilderAndSolver->BuildRHS(pScheme, BaseType::GetModelPart(), b );
        double rhDisp;
        double rhLM;
        ComputeMixedResidual(b, rhDisp, rhLM);
 
        // Compute full step residual (add another half Dx to the previous half)
        NRBaseType::UpdateDatabase(A,aux,b,MoveMesh);
        TSparseSpace::SetToZero(b);
        pBuilderAndSolver->BuildRHS(pScheme, BaseType::GetModelPart(), b );
        double rfDisp;
        double rfLM;
        ComputeMixedResidual(b, rfDisp, rfLM);
 
        // We compute the parabola
        double XminDisp = 1e-3;
        double XmaxDisp = 1.0;
        double XminLM = 1e-3;
        double XmaxLM = 1.0;
 
        ComputeParabola(XminDisp, XmaxDisp, rfDisp, roDisp, rhDisp);
        ComputeParabola(XminLM, XmaxLM, rfLM, roLM, rhLM);
 
        // Perform final update
        TSparseSpace::Assign(aux,-(1.0 - XmaxDisp), DxDisp);
        TSparseSpace::UnaliasedAdd(aux,-(1.0 - XmaxLM), DxLM);
        NRBaseType::UpdateDatabase(A,aux,b,MoveMesh);
    }
 
    void ComputeSplitDx(
        TSystemVectorType& Dx,
        TSystemVectorType& DxDisp,
        TSystemVectorType& DxLM
        )
    {
        // Now we iterate over all the nodes
        NodesArrayType& r_nodes_array = StrategyBaseType::GetModelPart().Nodes();
        block_for_each(r_nodes_array, [&](Node& rNode) {
            for(auto itDoF = rNode.GetDofs().begin() ; itDoF != rNode.GetDofs().end() ; itDoF++) {
                const int j = (**itDoF).EquationId();
                const std::size_t CurrVar = (**itDoF).GetVariable().Key();
 
                if ((CurrVar == DISPLACEMENT_X) || (CurrVar == DISPLACEMENT_Y) || (CurrVar == DISPLACEMENT_Z)) {
                    DxDisp[j] = Dx[j];
                    DxLM[j] = 0.0;
                } else { // Corresponding with contact
                    DxDisp[j] = 0.0;
                    DxLM[j] = Dx[j];
                }
            }
        });
    }
 
    void ComputeMixedResidual(
        TSystemVectorType& b,
        double& normDisp,
        double& normLM
        )
    {
        // Now we iterate over all the nodes
        NodesArrayType& r_nodes_array = StrategyBaseType::GetModelPart().Nodes();
        block_for_each(r_nodes_array, [&](Node& rNode) {
            for(auto itDoF = rNode.GetDofs().begin() ; itDoF != rNode.GetDofs().end() ; itDoF++) {
                const int j = (**itDoF).EquationId();
                const std::size_t CurrVar = (**itDoF).GetVariable().Key();
 
                if ((CurrVar == DISPLACEMENT_X) || (CurrVar == DISPLACEMENT_Y) || (CurrVar == DISPLACEMENT_Z)) {
                    AtomicAdd(normDisp, b[j] * b[j]);
                } else { // Corresponding with contact
                    AtomicAdd(normLM, b[j] * b[j]);
                }
            }
        });
 
        normDisp = std::sqrt(normDisp);
        normLM = std::sqrt(normLM);
    }
 
    void ComputeParabola(
        double& Xmax,
        double& Xmin,
        const double rf,
        const double ro,
        const double rh
        )
    {
        // Compute optimal (limited to the range 0-1)
        // Parabola is y = a*x^2 + b*x + c -> min/max for
        // x=0   --> r=ro
        // x=1/2 --> r=rh
        // x=1   --> r =
        // c= ro,     b= 4*rh -rf -3*ro,  a= 2*rf - 4*rh + 2*ro
        // max found if a>0 at the position  Xmax = (rf/4 - rh)/(rf - 2*rh);
 
        const double parabole_a = 2 * rf + 2 * ro - 4 * rh;
        const double parabole_b = 4 * rh - rf - 3 * ro;
 
        if( parabole_a > 0.0) //  If parabola has a local minima
        {
            Xmax = -0.5 * parabole_b/parabole_a; // -b / 2a
            if( Xmax > 1.0)
                Xmax = 1.0;
            else if(Xmax < -1.0)
                Xmax = -1.0;
        }
        else // Parabola degenerates to either a line or to have a local max. best solution on either extreme
        {
            if(rf < ro)
                Xmax = 1.0;
            else
                Xmax = Xmin; // Should be zero, but otherwise it will stagnate
        }
    }
 
    void AssignSettings(const Parameters ThisParameters) override
    {
        BaseType::AssignSettings(ThisParameters);
    }
 
 
 
 
    LineSearchContactStrategy(const LineSearchContactStrategy& Other)
    {
    };
 
private:
 
 
 
 
 
 
 
 
}; /* Class LineSearchContactStrategy */
}  // namespace Kratos