[Core] Initialize bossak constants in Predict

applications/ContactMechanicsApplication/custom_strategies/schemes/residual_based_bossak_displacement_rotation_scheme.hpp

@@ -542,6 +542,35 @@ class ResidualBasedBossakDisplacementRotationScheme: public Scheme<TSparseSpace,
 
         KRATOS_TRY
 
+        double DeltaTime = CurrentProcessInfo[DELTA_TIME];
+
+        if (DeltaTime == 0)
+            KRATOS_THROW_ERROR(std::logic_error, "detected delta_time = 0 in the Solution Scheme ... check if the time step is created correctly for the current model part", "" )
+
+        //Set Newmark coefficients
+
+        if( mNewmark.static_dynamic != 0 ) {
+          CurrentProcessInfo[NEWMARK_BETA]  = mNewmark.beta;
+          CurrentProcessInfo[NEWMARK_GAMMA] = mNewmark.gamma;
+          CurrentProcessInfo[BOSSAK_ALPHA]  = mAlpha.m;
+          CurrentProcessInfo[COMPUTE_DYNAMIC_TANGENT] = true;
+        }
+
+        //initializing Newmark constants
+        mNewmark.deltatime = DeltaTime;
+
+        mNewmark.c0 = ( mNewmark.gamma / ( DeltaTime * mNewmark.beta ) );
+        mNewmark.c1 = ( 1.0 / ( DeltaTime * DeltaTime * mNewmark.beta ) );
+
+        mNewmark.c2 = ( DeltaTime * ( 1.0 - mNewmark.gamma ) );
+        mNewmark.c3 = ( DeltaTime * mNewmark.gamma );
+        mNewmark.c4 = ( DeltaTime / mNewmark.beta );
+        mNewmark.c5 = ( DeltaTime * DeltaTime * ( 0.5 - mNewmark.beta ) / mNewmark.beta );
+
+        //extra
+        mNewmark.c6 = ( 1.0 / (mNewmark.beta * DeltaTime) );
+        mNewmark.c7 = ( 0.5 / (mNewmark.beta) - 1.0 );
+
         //std::cout << " Prediction " << std::endl;
 
         //double DeltaTime = r_model_part.GetProcessInfo()[DELTA_TIME];

applications/FluidDynamicsApplication/custom_strategies/schemes/residualbased_predictorcorrector_velocity_bossak_scheme_turbulent.h

@@ -357,6 +357,20 @@ namespace Kratos {
             // if (rModelPart.GetCommunicator().MyPID() == 0)
             //     std::cout << "prediction" << std::endl;
 
+            ProcessInfo& CurrentProcessInfo = rModelPart.GetProcessInfo();
+
+            double DeltaTime = CurrentProcessInfo[DELTA_TIME];
+            KRATOS_ERROR_IF(DeltaTime < 1.0e-12) << "Detected delta_time = 0 in the Bossak scheme. Check if the time step is created correctly for the current model part" << std::endl;
+
+            //initializing constants
+            ma0 = 1.0 / (mGammaNewmark * DeltaTime);
+            ma1 = DeltaTime * mBetaNewmark / mGammaNewmark;
+            ma2 = (-1 + mGammaNewmark) / mGammaNewmark;
+            ma3 = DeltaTime;
+            ma4 = pow(DeltaTime, 2)*(-2.0 * mBetaNewmark + 1.0) / 2.0;
+            ma5 = pow(DeltaTime, 2) * mBetaNewmark;
+            mam = (1.0 - mAlphaBossak) / (mGammaNewmark * DeltaTime);
+
             int NumThreads = ParallelUtilities::GetNumThreads();
             OpenMPUtils::PartitionVector NodePartition;
             OpenMPUtils::DivideInPartitions(rModelPart.Nodes().size(), NumThreads, NodePartition);

applications/PfemSolidMechanicsApplication/custom_strategies/schemes/residual_based_bossak_scheme.hpp

@@ -237,6 +237,27 @@ class ResidualBasedBossakScheme: public Scheme<TSparseSpace,TDenseSpace>
         TSystemVectorType& b
     ) override
     {
+        ProcessInfo CurrentProcessInfo= r_model_part.GetProcessInfo();
+
+        double DeltaTime = CurrentProcessInfo[DELTA_TIME];
+
+        if (DeltaTime == 0) {
+            std::cout<<" WARNING: detected delta_time = 0 in the Solution Scheme "<<std::endl;
+            std::cout<<" DELTA_TIME set to 1 considering a Quasistatic step with one step only "<<std::endl;
+            std::cout<<" PLEASE : check if the time step is created correctly for the current model part "<<std::endl;
+
+            CurrentProcessInfo[DELTA_TIME] = 1;
+            DeltaTime = CurrentProcessInfo[DELTA_TIME];
+        }
+
+        //initializing Newmark constants
+        mNewmark.c0 = ( 1.0 / (mNewmark.beta * DeltaTime * DeltaTime) );
+        mNewmark.c1 = ( mNewmark.gamma / (mNewmark.beta * DeltaTime) );
+        mNewmark.c2 = ( 1.0 / (mNewmark.beta * DeltaTime) );
+        mNewmark.c3 = ( 0.5 / (mNewmark.beta) - 1.0 );
+        mNewmark.c4 = ( (mNewmark.gamma / mNewmark.beta) - 1.0  );
+        mNewmark.c5 = ( DeltaTime * 0.5 * ( ( mNewmark.gamma / mNewmark.beta ) - 2 ) );
+
         //std::cout << " Prediction " << std::endl;
         array_1d<double, 3 > DeltaDisplacement;
 
@@ -681,7 +702,7 @@ class ResidualBasedBossakScheme: public Scheme<TSparseSpace,TDenseSpace>
         Condition::Pointer rCurrentCondition,
         LocalSystemVectorType& RHS_Contribution,
         Element::EquationIdVectorType& EquationId,
-        ProcessInfo& CurrentProcessInfo) 
+        ProcessInfo& CurrentProcessInfo)
     {
         KRATOS_TRY
 
@@ -723,7 +744,7 @@ class ResidualBasedBossakScheme: public Scheme<TSparseSpace,TDenseSpace>
     virtual void GetElementalDofList(
         Element::Pointer rCurrentElement,
         Element::DofsVectorType& ElementalDofList,
-        ProcessInfo& CurrentProcessInfo) 
+        ProcessInfo& CurrentProcessInfo)
     {
         rCurrentElement->GetDofList(ElementalDofList, CurrentProcessInfo);
     }

applications/RANSApplication/custom_strategies/bossak_relaxation_scalar_scheme.h

@@ -96,6 +96,27 @@ class BossakRelaxationScalarScheme
     ///@name Operations
     ///@{
 
+    void Predict(
+        ModelPart& rModelPart,
+        SystemMatrixType& rA,
+        SystemVectorType& rDx,
+        SystemVectorType& rb) override
+    {
+        BaseType::Predict(rModelPart, rA, rDx, rb);
+
+        const double delta_time = rModelPart.GetProcessInfo()[DELTA_TIME];
+
+        KRATOS_ERROR_IF(delta_time < std::numeric_limits<double>::epsilon())
+            << "detected delta_time = 0 in the Bossak Scheme ... "
+               "check if the time step is created correctly for "
+               "the current model part.";
+
+        BossakRelaxationScalarScheme::CalculateBossakConstants(
+            mBossak, mAlphaBossak, delta_time);
+
+        rModelPart.GetProcessInfo()[BOSSAK_ALPHA] = mBossak.Alpha;
+    }
+
     void InitializeSolutionStep(
         ModelPart& rModelPart,
         SystemMatrixType& rA,

applications/ULFapplication/custom_strategies/schemes/residualbased_predictorcorrector_bossak_scheme.h

@@ -283,6 +283,18 @@ class ResidualBasedPredictorCorrectorBossakScheme : public Scheme<TSparseSpace,
         array_1d<double, 3 > DeltaDisp;
         double DeltaTime = r_model_part.GetProcessInfo()[DELTA_TIME];
 
+        if (DeltaTime == 0)
+            KRATOS_THROW_ERROR(std::logic_error, "detected delta_time = 0 in the Bossak Scheme ... check if the time step is created correctly for the current model part", "");
+
+        //initializing constants
+        ma0 = 1.0 / (mBetaNewmark * pow(DeltaTime, 2));
+        ma1 = mGammaNewmark / (mBetaNewmark * DeltaTime);
+        ma2 = 1.0 / (mBetaNewmark * DeltaTime);
+        ma3 = 1.0 / (2.0 * mBetaNewmark) - 1.0;
+        ma4 = mGammaNewmark / mBetaNewmark - 1.0;
+        ma5 = DeltaTime * 0.5 * (mGammaNewmark / mBetaNewmark - 2.0);
+        mam = (1.0 - mAlphaBossak) / (mBetaNewmark * pow(DeltaTime, 2));
+
 
         for (ModelPart::NodeIterator i = r_model_part.NodesBegin();
                 i != r_model_part.NodesEnd(); ++i)

kratos/solving_strategies/schemes/residual_based_adjoint_bossak_scheme.h

@@ -174,6 +174,23 @@ class ResidualBasedAdjointBossakScheme : public Scheme<TSparseSpace, TDenseSpace
         KRATOS_CATCH("");
     }
 
+    void Predict(
+        ModelPart& rModelPart,
+        DofsArrayType& rDofSet,
+        SystemMatrixType& A,
+        SystemVectorType& Dx,
+        SystemVectorType& b) override
+    {
+        KRATOS_TRY;
+
+        const auto& r_current_process_info = rModelPart.GetProcessInfo();
+        mBossak = CalculateBossakConstants(mBossak.Alpha, GetTimeStep(r_current_process_info));
+
+        this->CalculateNodeNeighbourCount(rModelPart);
+
+        KRATOS_CATCH("");
+    }
+
     void InitializeSolutionStep(
         ModelPart& rModelPart,
         SystemMatrixType& rA,

kratos/solving_strategies/schemes/residual_based_bossak_displacement_scheme.hpp

@@ -236,6 +236,14 @@ class ResidualBasedBossakDisplacementScheme
         const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo();
         const double delta_time = r_current_process_info[DELTA_TIME];
 
+        // Initializing Bossak constants
+        mBossak.c0 = ( 1.0 / (mBossak.beta * delta_time * delta_time) );
+        mBossak.c1 = ( mBossak.gamma / (mBossak.beta * delta_time) );
+        mBossak.c2 = ( 1.0 / (mBossak.beta * delta_time) );
+        mBossak.c3 = ( 0.5 / (mBossak.beta) - 1.0 );
+        mBossak.c4 = ( (mBossak.gamma / mBossak.beta) - 1.0  );
+        mBossak.c5 = ( delta_time * 0.5 * ( ( mBossak.gamma / mBossak.beta ) - 2.0 ) );
+
         // Updating time derivatives (nodally for efficiency)
         if (rModelPart.Nodes().size() > 0) {
             const auto it_node_begin = rModelPart.Nodes().begin();