residualbased_newton_raphson_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 "contact_structural_mechanics_application_variables.h"
#include "includes/kratos_parameters.h"
#include "includes/define.h"
#include "includes/model_part.h"
#include "includes/variables.h"
 
/* Strategies */
#include "solving_strategies/strategies/residualbased_newton_raphson_strategy.h"
 
/* Utilities */
#include "utilities/variable_utils.h"
#include "utilities/color_utilities.h"
#include "utilities/math_utils.h"
#include "custom_python/process_factory_utility.h"
#include "custom_utilities/contact_utilities.h"
 
namespace Kratos {
 
 
 
 
 
 
template<class TSparseSpace,
         class TDenseSpace, // = DenseSpace<double>,
         class TLinearSolver //= LinearSolver<TSparseSpace,TDenseSpace>
         >
class ResidualBasedNewtonRaphsonContactStrategy :
    public ResidualBasedNewtonRaphsonStrategy< TSparseSpace, TDenseSpace, TLinearSolver >
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( ResidualBasedNewtonRaphsonContactStrategy );
 
    using SolvingStrategyType = SolvingStrategy<TSparseSpace, TDenseSpace>;
 
    using StrategyBaseType = ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver>;
 
    using BaseType = ResidualBasedNewtonRaphsonStrategy<TSparseSpace, TDenseSpace, TLinearSolver>;
 
    using ClassType = ResidualBasedNewtonRaphsonContactStrategy<TSparseSpace, TDenseSpace, TLinearSolver>;
 
    using TConvergenceCriteriaType = ConvergenceCriteria<TSparseSpace, TDenseSpace>;
 
    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 ElementsArrayType = typename ModelPart::ElementsContainerType;
 
    using ConditionsArrayType = typename ModelPart::ConditionsContainerType;
 
    using ProcessesListType = typename ProcessFactoryUtility::Pointer;
 
    using IndexType = std::size_t;
 
    explicit ResidualBasedNewtonRaphsonContactStrategy()
    {
    }
 
    explicit ResidualBasedNewtonRaphsonContactStrategy(ModelPart& rModelPart, Parameters ThisParameters)
        : BaseType(rModelPart),
          mpMyProcesses(nullptr),
          mpPostProcesses(nullptr)
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
 
        // Auxiliary assign
        mConvergenceCriteriaEchoLevel = BaseType::mpConvergenceCriteria->GetEchoLevel();
    }
 
    ResidualBasedNewtonRaphsonContactStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
        typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver,
        IndexType MaxIterations = 30,
        bool CalculateReactions = false,
        bool ReformDofSetAtEachStep = false,
        bool MoveMeshFlag = false,
        Parameters ThisParameters =  Parameters(R"({})"),
        ProcessesListType pMyProcesses = nullptr,
        ProcessesListType pPostProcesses = nullptr
        )
        : BaseType(rModelPart, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag ),
        mThisParameters(ThisParameters),
        mpMyProcesses(pMyProcesses),
        mpPostProcesses(pPostProcesses)
    {
        KRATOS_TRY;
 
        mConvergenceCriteriaEchoLevel = pNewConvergenceCriteria->GetEchoLevel();
 
        Parameters default_parameters = GetDefaultParameters();
        mThisParameters.ValidateAndAssignDefaults(default_parameters);
 
        KRATOS_CATCH("");
    }
 
    ResidualBasedNewtonRaphsonContactStrategy(
        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"({})"),
        ProcessesListType pMyProcesses = nullptr,
        ProcessesListType pPostProcesses = nullptr
    )
        : BaseType(rModelPart, pScheme, pNewLinearSolver, pNewConvergenceCriteria, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag),
        mThisParameters(ThisParameters),
        mpMyProcesses(pMyProcesses),
        mpPostProcesses(pPostProcesses)
    {
        KRATOS_TRY;
 
        mConvergenceCriteriaEchoLevel = pNewConvergenceCriteria->GetEchoLevel();
 
        Parameters default_parameters = GetDefaultParameters();
        mThisParameters.ValidateAndAssignDefaults(default_parameters);
 
        KRATOS_CATCH("");
    }
 
    ResidualBasedNewtonRaphsonContactStrategy(
        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"({})"),
        ProcessesListType pMyProcesses = nullptr,
        ProcessesListType pPostProcesses = nullptr
        )
        : ResidualBasedNewtonRaphsonStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(rModelPart, pScheme, pNewLinearSolver, pNewConvergenceCriteria, pNewBuilderAndSolver, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag ),
        mThisParameters(ThisParameters),
        mpMyProcesses(pMyProcesses),
        mpPostProcesses(pPostProcesses)
    {
        KRATOS_TRY;
 
        mConvergenceCriteriaEchoLevel = pNewConvergenceCriteria->GetEchoLevel();
 
        Parameters default_parameters = GetDefaultParameters();
        mThisParameters.ValidateAndAssignDefaults(default_parameters);
 
        KRATOS_CATCH("");
    }
 
    ~ResidualBasedNewtonRaphsonContactStrategy() override
    = default;
 
 
 
    typename SolvingStrategyType::Pointer Create(
        ModelPart& rModelPart,
        Parameters ThisParameters
        ) const override
    {
        return Kratos::make_shared<ClassType>(rModelPart, ThisParameters);
    }
 
    void Predict() override
    {
        KRATOS_TRY
 
        // Auxiliary zero array
        const array_1d<double, 3> zero_array = ZeroVector(3);
 
        // Set to zero the weighted gap
        ModelPart& r_model_part = StrategyBaseType::GetModelPart();
        NodesArrayType& r_nodes_array = r_model_part.GetSubModelPart("Contact").Nodes();
        const bool frictional = r_model_part.Is(SLIP);
 
        // We predict contact pressure in case of contact problem
        if (r_nodes_array.begin()->SolutionStepsDataHas(WEIGHTED_GAP)) {
            VariableUtils().SetVariable(WEIGHTED_GAP, 0.0, r_nodes_array);
            if (frictional) {
                VariableUtils().SetVariable(WEIGHTED_SLIP, zero_array, r_nodes_array);
            }
 
            // Compute the current gap
            ContactUtilities::ComputeExplicitContributionConditions(r_model_part.GetSubModelPart("ComputingContact"));
 
            // We predict a contact pressure
            ProcessInfo& r_process_info = r_model_part.GetProcessInfo();
            const std::size_t step = r_process_info[STEP];
 
            if (step == 1) {
                block_for_each(r_nodes_array, [&](Node& rNode) {
                    noalias(rNode.Coordinates()) += rNode.FastGetSolutionStepValue(DISPLACEMENT);
                });
            } else {
                block_for_each(r_nodes_array, [&](Node& rNode) {
                    noalias(rNode.Coordinates()) += (rNode.FastGetSolutionStepValue(DISPLACEMENT) - rNode.FastGetSolutionStepValue(DISPLACEMENT, 1));
                });
            }
        }
 
//         BaseType::Predict();  // NOTE: May cause problems in dynamics!!!
//
//         // Set to zero the weighted gap // NOTE: This can be done during the search if the predict is deactivated
//         ModelPart& r_model_part = StrategyBaseType::GetModelPart();
//         NodesArrayType& nodes_array = r_model_part.GetSubModelPart("Contact").Nodes();
//
//         // We predict contact pressure in case of contact problem
//         if (nodes_array.begin()->SolutionStepsDataHas(WEIGHTED_GAP)) {
//             VariableUtils().SetVariable(WEIGHTED_GAP, 0.0, nodes_array);
//
//             // Compute the current gap
//             ContactUtilities::ComputeExplicitContributionConditions(r_model_part.GetSubModelPart("ComputingContact"));
//
//             // We predict a contact pressure
//             ProcessInfo& r_process_info = r_model_part.GetProcessInfo();
//             const double initial_penalty_parameter = r_process_info[INITIAL_PENALTY];
//
//             // We iterate over the nodes
//             const bool is_components = nodes_array.begin()->SolutionStepsDataHas(LAGRANGE_MULTIPLIER_CONTACT_PRESSURE) ? false : true;
//
//             block_for_each(r_nodes_array, [&initial_penalty_parameter, &is_components](Node& rNode) {
//                 const double current_gap = rNode.FastGetSolutionStepValue(WEIGHTED_GAP);
//                 const double penalty = rNode.Has(INITIAL_PENALTY) ? rNode.GetValue(INITIAL_PENALTY) : initial_penalty_parameter;
//                 if (current_gap < 0.0) {
//                     rNode.Set(ACTIVE, true);
//                     if (is_components) {
//                         rNode.FastGetSolutionStepValue(LAGRANGE_MULTIPLIER_CONTACT_PRESSURE) = penalty * current_gap;
//                     } else {
//                         const array_1d<double, 3>& normal = rNode.FastGetSolutionStepValue(NORMAL);
//                         rNode.FastGetSolutionStepValue(VECTOR_LAGRANGE_MULTIPLIER) = penalty * current_gap * normal;
//                     }
//                 }
//             });
//         }
 
        KRATOS_CATCH("")
    }
 
    void Initialize() override
    {
        KRATOS_TRY;
 
        BaseType::Initialize();
        mFinalizeWasPerformed = false;
 
        // Initializing NL_ITERATION_NUMBER
        ModelPart& r_model_part = StrategyBaseType::GetModelPart();
        ProcessInfo& r_process_info = r_model_part.GetProcessInfo();
        r_process_info[NL_ITERATION_NUMBER] = 1;
 
        KRATOS_CATCH("");
    }
 
    double Solve() override
    {
        this->Initialize();
        this->InitializeSolutionStep();
        this->Predict();
        this->SolveSolutionStep();
        this->FinalizeSolutionStep();
 
        // TODO: Add something if necessary
 
        return 0.0;
    }
 
    void InitializeSolutionStep() override
    {
        BaseType::mpConvergenceCriteria->SetEchoLevel(0);
        BaseType::InitializeSolutionStep();
        BaseType::mpConvergenceCriteria->SetEchoLevel(mConvergenceCriteriaEchoLevel);
 
        mFinalizeWasPerformed = false;
    }
 
    void FinalizeSolutionStep() override
    {
        KRATOS_TRY;
 
        if (mFinalizeWasPerformed == false) {
            BaseType::FinalizeSolutionStep();
 
            // To avoid compute twice the FinalizeSolutionStep
            mFinalizeWasPerformed = true;
        }
 
        KRATOS_CATCH("");
    }
 
    bool SolveSolutionStep() override
    {
        KRATOS_TRY;
 
//         bool is_converged = BaseType::SolveSolutionStep(); // FIXME: Requires to separate the non linear iterations
 
//         bool is_converged = BaseSolveSolutionStep(); // Direct solution
        bool is_converged = false;
 
        // Getting model part
        ModelPart& r_model_part = StrategyBaseType::GetModelPart();
 
        if (r_model_part.IsNot(INTERACTION)) {
            // We get the system
            TSystemMatrixType& A = *BaseType::mpA;
            TSystemVectorType& Dx = *BaseType::mpDx;
            TSystemVectorType& b = *BaseType::mpb;
 
            // We get the process info
            ProcessInfo& r_process_info = r_model_part.GetProcessInfo();
 
            int inner_iteration = 0;
            while (!is_converged && inner_iteration < mThisParameters["inner_loop_iterations"].GetInt()) {
                ++inner_iteration;
 
                if (mConvergenceCriteriaEchoLevel > 0 && StrategyBaseType::GetModelPart().GetCommunicator().MyPID() == 0 ) {
                    std::cout << std::endl << BOLDFONT("Simplified semi-smooth strategy. INNER ITERATION: ") << inner_iteration;;
                }
 
                // We solve one loop
                r_process_info[NL_ITERATION_NUMBER] = 1;
                r_process_info[INNER_LOOP_ITERATION] = inner_iteration;
                is_converged = BaseSolveSolutionStep();
 
                // We check the convergence
                BaseType::mpConvergenceCriteria->SetEchoLevel(0);
                is_converged = BaseType::mpConvergenceCriteria->PostCriteria(r_model_part, BaseType::GetBuilderAndSolver()->GetDofSet(), A, Dx, b);
                BaseType::mpConvergenceCriteria->SetEchoLevel(mConvergenceCriteriaEchoLevel);
 
                if (mConvergenceCriteriaEchoLevel > 0 && StrategyBaseType::GetModelPart().GetCommunicator().MyPID() == 0 ) {
                    if (is_converged) std::cout << BOLDFONT("Simplified semi-smooth strategy. INNER ITERATION: ") << BOLDFONT(FGRN("CONVERGED")) << std::endl;
                    else std::cout << BOLDFONT("Simplified semi-smooth strategy. INNER ITERATION: ") << BOLDFONT(FRED("NOT CONVERGED")) << std::endl;
                }
            }
        } else {
            // We compute the base loop
            r_model_part.GetProcessInfo()[INNER_LOOP_ITERATION] = 1;
            is_converged = BaseSolveSolutionStep();
        }
 
        if (mThisParameters["adaptative_strategy"].GetBool()) {
            if (!is_converged) {
                is_converged = AdaptativeStep();
            }
        }
 
        return is_converged;
 
        KRATOS_CATCH("");
    }
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name"                                : "newton_raphson_contact_strategy",
            "adaptative_strategy"                 : false,
            "split_factor"                        : 10.0,
            "max_number_splits"                   : 3,
            "inner_loop_iterations"               : 5
        })" );
 
        // 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 "newton_raphson_contact_strategy";
    }
 
 
 
 
 
protected:
 
 
 
    Parameters mThisParameters;        
 
    // ADAPTATIVE STRATEGY PARAMETERS
    bool mFinalizeWasPerformed;        
    ProcessesListType mpMyProcesses;   
    ProcessesListType mpPostProcesses; 
 
    // OTHER PARAMETERS
    int mConvergenceCriteriaEchoLevel; 
 
 
 
    void AssignSettings(const Parameters ThisParameters) override
    {
        BaseType::AssignSettings(ThisParameters);
 
        // Copy the parameters
        mThisParameters = ThisParameters;
    }
 
    bool BaseSolveSolutionStep()
    {
        KRATOS_TRY;
 
        // Pointers needed in the solution
        ModelPart& r_model_part = StrategyBaseType::GetModelPart();
        ProcessInfo& r_process_info = r_model_part.GetProcessInfo();
        typename TSchemeType::Pointer p_scheme = BaseType::GetScheme();
        typename TBuilderAndSolverType::Pointer p_builder_and_solver = BaseType::GetBuilderAndSolver();
        auto& r_dof_set = p_builder_and_solver->GetDofSet();
 
        TSystemMatrixType& rA = *BaseType::mpA;
        TSystemVectorType& rDx = *BaseType::mpDx;
        TSystemVectorType& rb = *BaseType::mpb;
 
        // Initializing the parameters of the Newton-Raphson cicle
        IndexType iteration_number = 1;
        r_process_info[NL_ITERATION_NUMBER] = iteration_number;
 
        bool is_converged = false;
        bool residual_is_updated = false;
        p_scheme->InitializeNonLinIteration(r_model_part, rA, rDx, rb);
        BaseType::mpConvergenceCriteria->InitializeNonLinearIteration(r_model_part, r_dof_set, rA, rDx, rb);
        is_converged = BaseType::mpConvergenceCriteria->PreCriteria(r_model_part, r_dof_set, rA, rDx, rb);
 
        // We do a geometry check before solve the system for first time
        if (mThisParameters["adaptative_strategy"].GetBool()) {
            if (CheckGeometryInverted()) {
                KRATOS_WARNING("Element inverted") << "INVERTED ELEMENT BEFORE FIRST SOLVE"  << std::endl;
                r_process_info[STEP] -= 1; // We revert one step in the case that the geometry is already broken before start the computing
                return false;
            }
        }
 
        // Function to perform the building and the solving phase.
        if (StrategyBaseType::mRebuildLevel > 1 || StrategyBaseType::mStiffnessMatrixIsBuilt == false) {
            TSparseSpace::SetToZero(rA);
            TSparseSpace::SetToZero(rDx);
            TSparseSpace::SetToZero(rb);
 
            p_builder_and_solver->BuildAndSolve(p_scheme, r_model_part, rA, rDx, rb);
        } else {
            TSparseSpace::SetToZero(rDx); //Dx=0.00;
            TSparseSpace::SetToZero(rb);
 
            p_builder_and_solver->BuildRHSAndSolve(p_scheme, r_model_part, rA, rDx, rb);
        }
 
        // Debugging info
        BaseType::EchoInfo(iteration_number);
 
        // Updating the results stored in the database
        UpdateDatabase(rA, rDx, rb, StrategyBaseType::MoveMeshFlag());
 
        // We now check the geometry
        if (mThisParameters["adaptative_strategy"].GetBool()) {
            if (CheckGeometryInverted()) {
                KRATOS_WARNING("Element inverted") << "INVERTED ELEMENT DURING DATABASE UPDATE" << std::endl;
                r_process_info[STEP] -= 1; // We revert one step in the case that the geometry is already broken before start the computing
                return false;
            }
        }
 
        p_scheme->FinalizeNonLinIteration(r_model_part, rA, rDx, rb);
        BaseType::mpConvergenceCriteria->FinalizeNonLinearIteration(r_model_part, r_dof_set, rA, rDx, rb);
 
        if (is_converged) {
            // Initialisation of the convergence criteria
            BaseType::mpConvergenceCriteria->InitializeSolutionStep(r_model_part, r_dof_set, rA, rDx, rb);
 
            if (BaseType::mpConvergenceCriteria->GetActualizeRHSflag()) {
                TSparseSpace::SetToZero(rb);
 
                p_builder_and_solver->BuildRHS(p_scheme, r_model_part, rb);
            }
 
            is_converged = BaseType::mpConvergenceCriteria->PostCriteria(r_model_part, r_dof_set, rA, rDx, rb);
        }
 
        // Iteration Cicle... performed only for NonLinearProblems
        while (is_converged == false && iteration_number++<BaseType::mMaxIterationNumber) {
            //setting the number of iteration
            r_process_info[NL_ITERATION_NUMBER] = iteration_number;
 
            p_scheme->InitializeNonLinIteration(r_model_part, rA, rDx, rb);
            BaseType::mpConvergenceCriteria->InitializeNonLinearIteration(r_model_part, r_dof_set, rA, rDx, rb);
 
            is_converged = BaseType::mpConvergenceCriteria->PreCriteria(r_model_part, r_dof_set, rA, rDx, rb);
 
            //call the linear system solver to find the correction mDx for the
            //it is not called if there is no system to solve
            if (SparseSpaceType::Size(rDx) != 0) {
                if (StrategyBaseType::mRebuildLevel > 1 || StrategyBaseType::mStiffnessMatrixIsBuilt == false ) {
                    if( BaseType::GetKeepSystemConstantDuringIterations() == false) {
                        //A = 0.00;
                        TSparseSpace::SetToZero(rA);
                        TSparseSpace::SetToZero(rDx);
                        TSparseSpace::SetToZero(rb);
 
                        p_builder_and_solver->BuildAndSolve(p_scheme, r_model_part, rA, rDx, rb);
                    }
                    else {
                        TSparseSpace::SetToZero(rDx);
                        TSparseSpace::SetToZero(rb);
 
                        p_builder_and_solver->BuildRHSAndSolve(p_scheme, r_model_part, rA, rDx, rb);
                    }
                }
                else {
                    TSparseSpace::SetToZero(rDx);
                    TSparseSpace::SetToZero(rb);
 
                    p_builder_and_solver->BuildRHSAndSolve(p_scheme, r_model_part, rA, rDx, rb);
                }
            } else {
                KRATOS_WARNING("No DoFs") << "ATTENTION: no free DOFs!! " << std::endl;
            }
 
            // Debugging info
            BaseType::EchoInfo(iteration_number);
 
            // Updating the results stored in the database
            UpdateDatabase(rA, rDx, rb, StrategyBaseType::MoveMeshFlag());
 
            // We now check the geometry
            if (mThisParameters["adaptative_strategy"].GetBool()) {
                if (CheckGeometryInverted()) {
                    KRATOS_WARNING("Element inverted") << "INVERTED ELEMENT DURING DATABASE UPDATE" << std::endl;
                    r_process_info[STEP] -= 1; // We revert one step in the case that the geometry is already broken before start the computing
                    return false;
                }
            }
 
            p_scheme->FinalizeNonLinIteration(r_model_part, rA, rDx, rb);
            BaseType::mpConvergenceCriteria->FinalizeNonLinearIteration(r_model_part, r_dof_set, rA, rDx, rb);
 
            residual_is_updated = false;
 
            if (is_converged) {
 
                if (BaseType::mpConvergenceCriteria->GetActualizeRHSflag()) {
                    TSparseSpace::SetToZero(rb);
 
                    p_builder_and_solver->BuildRHS(p_scheme, r_model_part, rb);
                    residual_is_updated = true;
                    //std::cout << "mb is calculated" << std::endl;
                }
 
                is_converged = BaseType::mpConvergenceCriteria->PostCriteria(r_model_part, r_dof_set, rA, rDx, rb);
            }
        }
 
        // Plots a warning if the maximum number of iterations is exceeded
        if (iteration_number >= BaseType::mMaxIterationNumber && r_model_part.GetCommunicator().MyPID() == 0)
            MaxIterationsExceeded();
 
        // Recalculate residual if needed
        // (note that some convergence criteria need it to be recalculated)
        if (residual_is_updated == false) {
            // NOTE:
            // The following part will be commented because it is time consuming
            // and there is no obvious reason to be here. If someone need this
            // part please notify the community via mailing list before uncommenting it.
            // Pooyan.
 
            //    TSparseSpace::SetToZero(mb);
            //    p_builder_and_solver->BuildRHS(p_scheme, r_model_part, mb);
        }
 
        // Calculate reactions if required
        if (BaseType::mCalculateReactionsFlag)
            p_builder_and_solver->CalculateReactions(p_scheme, r_model_part, rA, rDx, rb);
 
        return is_converged;
 
        KRATOS_CATCH("");
    }
 
    bool AdaptativeStep()
    {
        KRATOS_TRY;
 
        bool is_converged = false;
        // Plots a warning if the maximum number of iterations is exceeded
        if (mpMyProcesses == nullptr && StrategyBaseType::mEchoLevel > 0)
            KRATOS_WARNING("No python processes") << "If you have not implemented any method to recalculate BC or loads in function of time, this strategy will be USELESS" << std::endl;
 
        if (mpPostProcesses == nullptr && StrategyBaseType::mEchoLevel > 0)
            KRATOS_WARNING("No python post processes") << "If you don't add the postprocesses and the time step if splitted you won't postprocess that steps" << std::endl;
 
        ModelPart& r_model_part = StrategyBaseType::GetModelPart();
        ProcessInfo& r_process_info = r_model_part.GetProcessInfo();
 
        const double original_delta_time = r_process_info[DELTA_TIME]; // We save the delta time to restore later
 
        int split_number = 0;
 
        // We iterate until we reach the convergence or we split more than desired
        while (is_converged == false && split_number <= mThisParameters["max_number_splits"].GetInt()) {
            // Expliting time step as a way to try improve the convergence
            split_number += 1;
            double aux_delta_time, current_time;
            const double aux_time = SplitTimeStep(aux_delta_time, current_time);
            current_time += aux_delta_time;
 
            bool inside_the_split_is_converged = false;
            IndexType inner_iteration = 0;
            while (current_time <= aux_time) {
                inner_iteration += 1;
                r_process_info[STEP] += 1;
 
                if (inner_iteration == 1) {
                    if (StrategyBaseType::MoveMeshFlag())
                        UnMoveMesh();
 
                    NodesArrayType& r_nodes_array = r_model_part.Nodes();
                    block_for_each(r_nodes_array, [&](Node& rNode) {
                        rNode.OverwriteSolutionStepData(1, 0);
//                         rNode.OverwriteSolutionStepData(2, 1);
                    });
 
                    r_process_info.SetCurrentTime(current_time); // Reduces the time step
 
                    FinalizeSolutionStep();
                } else {
                    NodesArrayType& r_nodes_array = r_model_part.Nodes();
                    block_for_each(r_nodes_array, [&](Node& rNode) {
                        rNode.CloneSolutionStepData();
                    });
 
                    r_process_info.CloneSolutionStepInfo();
                    r_process_info.ClearHistory(r_model_part.GetBufferSize());
                    r_process_info.SetAsTimeStepInfo(current_time); // Sets the new time step
                }
 
                // We execute the processes before the non-linear iteration
                if (mpMyProcesses != nullptr)
                    mpMyProcesses->ExecuteInitializeSolutionStep();
 
                if (mpPostProcesses != nullptr)
                    mpPostProcesses->ExecuteInitializeSolutionStep();
 
                // In order to initialize again everything
                BaseType::mInitializeWasPerformed = false;
                mFinalizeWasPerformed = false;
 
                // We repeat the solve with the new DELTA_TIME
                this->Initialize();
                this->InitializeSolutionStep();
                this->Predict();
                inside_the_split_is_converged = BaseType::SolveSolutionStep();
                this->FinalizeSolutionStep();
 
                // We execute the processes after the non-linear iteration
                if (mpMyProcesses != nullptr)
                    mpMyProcesses->ExecuteFinalizeSolutionStep();
 
                if (mpPostProcesses != nullptr)
                    mpPostProcesses->ExecuteFinalizeSolutionStep();
 
                if (mpMyProcesses != nullptr)
                    mpMyProcesses->ExecuteBeforeOutputStep();
 
                if (mpPostProcesses != nullptr)
                    mpPostProcesses->PrintOutput();
 
                if (mpMyProcesses != nullptr)
                    mpMyProcesses->ExecuteAfterOutputStep();
 
                current_time += aux_delta_time;
            }
 
            if (inside_the_split_is_converged)
                is_converged = true;
        }
 
        // Plots a warning if the maximum number of iterations and splits are exceeded
        if (is_converged == false)
            MaxIterationsAndSplitsExceeded();
 
        // Restoring original DELTA_TIME
        r_process_info[DELTA_TIME] = original_delta_time;
 
        return is_converged;
 
        KRATOS_CATCH("");
    }
 
    void UpdateDatabase(
        TSystemMatrixType& A,
        TSystemVectorType& Dx,
        TSystemVectorType& b,
        const bool MoveMesh
        ) override
    {
        BaseType::UpdateDatabase(A,Dx,b,MoveMesh);
 
        // TODO: Add something if necessary
    }
 
    bool CheckGeometryInverted()
    {
        ModelPart& r_model_part = StrategyBaseType::GetModelPart();
        ProcessInfo& r_process_info = r_model_part.GetProcessInfo();
        bool inverted_element = false;
 
        ElementsArrayType& elements_array = r_model_part.Elements();
 
        // NOT OMP
        for(int i = 0; i < static_cast<int>(elements_array.size()); ++i) {
            auto it_elem = elements_array.begin() + i;
            auto& geom = it_elem->GetGeometry();
            if (geom.DeterminantOfJacobian(0) < 0.0) {
                if (mConvergenceCriteriaEchoLevel > 0) {
                    KRATOS_WATCH(it_elem->Id())
                    KRATOS_WATCH(geom.DeterminantOfJacobian(0))
                }
                return true;
            }
 
            // We check now the deformation gradient
            std::vector<Matrix> deformation_gradient_matrices;
            it_elem->CalculateOnIntegrationPoints( DEFORMATION_GRADIENT, deformation_gradient_matrices, r_process_info);
 
            for (IndexType i_gp = 0; i_gp  < deformation_gradient_matrices.size(); ++i_gp) {
                const double det_f = MathUtils<double>::Det(deformation_gradient_matrices[i_gp]);
                if (det_f < 0.0) {
                    if (mConvergenceCriteriaEchoLevel > 0) {
                        KRATOS_WATCH(it_elem->Id())
                        KRATOS_WATCH(det_f)
                    }
                    return true;
                }
            }
        }
 
        return inverted_element;
    }
 
    double SplitTimeStep(
        double& AuxDeltaTime,
        double& CurrentTime
        )
    {
        KRATOS_TRY;
 
        const double aux_time = StrategyBaseType::GetModelPart().GetProcessInfo()[TIME];
        AuxDeltaTime = StrategyBaseType::GetModelPart().GetProcessInfo()[DELTA_TIME];
        CurrentTime = aux_time - AuxDeltaTime;
 
        StrategyBaseType::GetModelPart().GetProcessInfo()[TIME] = CurrentTime; // Restore time to the previous one
        AuxDeltaTime /= mThisParameters["split_factor"].GetDouble();
        StrategyBaseType::GetModelPart().GetProcessInfo()[DELTA_TIME] = AuxDeltaTime; // Change delta time
 
        CoutSplittingTime(AuxDeltaTime, aux_time);
 
        return aux_time;
 
        KRATOS_CATCH("");
    }
 
    void UnMoveMesh()
    {
        KRATOS_TRY;
 
        if (StrategyBaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(DISPLACEMENT_X) == false)
            KRATOS_ERROR << "It is impossible to move the mesh since the DISPLACEMENT var is not in the model_part. Either use SetMoveMeshFlag(False) or add DISPLACEMENT to the list of variables" << std::endl;
 
        NodesArrayType& r_nodes_array = StrategyBaseType::GetModelPart().Nodes();
        block_for_each(r_nodes_array, [&](Node& rNode) {
            noalias(rNode.Coordinates()) = rNode.GetInitialPosition().Coordinates();
            noalias(rNode.Coordinates()) += rNode.FastGetSolutionStepValue(DISPLACEMENT, 1);
        });
 
        KRATOS_CATCH("");
    }
 
    void CoutSolvingProblem()
    {
        if (mConvergenceCriteriaEchoLevel != 0) {
            std::cout << "STEP: " << StrategyBaseType::GetModelPart().GetProcessInfo()[STEP] << "\t NON LINEAR ITERATION: " << StrategyBaseType::GetModelPart().GetProcessInfo()[NL_ITERATION_NUMBER] << "\t TIME: " << StrategyBaseType::GetModelPart().GetProcessInfo()[TIME] << "\t DELTA TIME: " << StrategyBaseType::GetModelPart().GetProcessInfo()[DELTA_TIME]  << std::endl;
        }
    }
 
    void CoutSplittingTime(
        const double AuxDeltaTime,
        const double AuxTime
        )
    {
        if (mConvergenceCriteriaEchoLevel > 0 && StrategyBaseType::GetModelPart().GetCommunicator().MyPID() == 0 ) {
            const double Time = StrategyBaseType::GetModelPart().GetProcessInfo()[TIME];
            std::cout.precision(4);
            std::cout << "|----------------------------------------------------|" << std::endl;
            std::cout << "|     " << BOLDFONT("SPLITTING TIME STEP") << "                            |" << std::endl;
            std::cout << "| " << BOLDFONT("COMING BACK TO TIME: ") << std::scientific << Time << "                    |" << std::endl;
            std::cout << "| " << BOLDFONT("      NEW TIME STEP: ") << std::scientific << AuxDeltaTime << "                    |" << std::endl;
            std::cout << "| " << BOLDFONT("         UNTIL TIME: ") << std::scientific << AuxTime << "                    |" << std::endl;
            std::cout << "|----------------------------------------------------|" << std::endl;
        }
    }
 
    void MaxIterationsExceeded() override
    {
        if (mConvergenceCriteriaEchoLevel > 0 && StrategyBaseType::GetModelPart().GetCommunicator().MyPID() == 0 ) {
            std::cout << "|----------------------------------------------------|" << std::endl;
            std::cout << "|        " << BOLDFONT(FRED("ATTENTION: Max iterations exceeded")) << "          |" << std::endl;
            std::cout << "|----------------------------------------------------|" << std::endl;
        }
    }
 
    void MaxIterationsAndSplitsExceeded()
    {
        if (mConvergenceCriteriaEchoLevel > 0 && StrategyBaseType::GetModelPart().GetCommunicator().MyPID() == 0 ) {
            std::cout << "|----------------------------------------------------|" << std::endl;
            std::cout << "|        " << BOLDFONT(FRED("ATTENTION: Max iterations exceeded")) << "          |" << std::endl;
            std::cout << "|        " << BOLDFONT(FRED("   Max number of splits exceeded  ")) << "          |" << std::endl;
            std::cout << "|----------------------------------------------------|" << std::endl;
        }
    }
 
 
 
 
    ResidualBasedNewtonRaphsonContactStrategy(const ResidualBasedNewtonRaphsonContactStrategy& Other)
    {
    };
 
 
}; /* Class ResidualBasedNewtonRaphsonContactStrategy */
}  // namespace Kratos