From 9c0e16cdddc87346e1575d4a53784072ca68df47 Mon Sep 17 00:00:00 2001 From: JonasHill Date: Fri, 22 Nov 2024 14:49:02 +0100 Subject: [PATCH 1/8] added Update method for rotational dofs --- .../structural_mechanics_bossak_scheme.h | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h index 5e0196bf9585..5cf32f582a15 100644 --- a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h +++ b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h @@ -121,6 +121,38 @@ class StructuralMechanicsBossakScheme }); } } + void Update( + ModelPart& rModelPart, + DofsArrayType& rDofSet, + TSystemMatrixType& rA, + TSystemVectorType& rDx, + TSystemVectorType& rb + ) override + { + KRATOS_TRY + + // Call the base Update method + BossakBaseType::Update(rModelPart, rDofSet, rA, rDx, rb); + + // Updating angular time derivatives if they are available (nodally for efficiency) + if(rModelPart.HasNodalSolutionStepVariable(ROTATION)) + { + block_for_each(rModelPart.Nodes(), array_1d(), [&](Node& rNode, array_1d& rDeltaRotationTLS){ + noalias(rDeltaRotationTLS) = rNode.FastGetSolutionStepValue(ROTATION) - rNode.FastGetSolutionStepValue(ROTATION, 1); + + array_1d& r_current_angular_velocity = rNode.FastGetSolutionStepValue(ANGULAR_VELOCITY); + const array_1d& r_previous_angular_velocity = rNode.FastGetSolutionStepValue(ANGULAR_VELOCITY, 1); + + array_1d& r_current_angular_acceleration = rNode.FastGetSolutionStepValue(ANGULAR_ACCELERATION); + const array_1d& r_previous_angular_acceleration = rNode.FastGetSolutionStepValue(ANGULAR_ACCELERATION, 1); + + UpdateAngularVelocity(r_current_angular_velocity, rDeltaRotationTLS, r_previous_angular_velocity, r_previous_angular_acceleration); + UpdateAngularAcceleration(r_current_angular_acceleration, rDeltaRotationTLS, r_previous_angular_velocity, r_previous_angular_acceleration); + }); + } + + KRATOS_CATCH("") + } void InitializeSolutionStep( ModelPart &rModelPart, @@ -261,6 +293,26 @@ class StructuralMechanicsBossakScheme } } + inline void UpdateAngularVelocity( + array_1d& rCurrentAngularVelocity, + const array_1d& rDeltaRotation, + const array_1d& rPreviousAngularVelocity, + const array_1d& rPreviousAngularAcceleration + ) + { + noalias(rCurrentAngularVelocity) = (this->mBossak.c1 * rDeltaRotation - this->mBossak.c4 * rPreviousAngularVelocity - this->mBossak.c5 * rPreviousAngularAcceleration); + } + + inline void UpdateAngularAcceleration( + array_1d& rCurrentAngularAcceleration, + const array_1d& rDeltaRotation, + const array_1d& rPreviousAngularVelocity, + const array_1d& rPreviousAngularAcceleration + ) + { + noalias(rCurrentAngularAcceleration) = (this->mBossak.c0 * rDeltaRotation - this->mBossak.c2 * rPreviousAngularVelocity - this->mBossak.c3 * rPreviousAngularAcceleration); + } + ///@} private: ///@name Private Member Variables From 261303c40ad9c18490371024f70b7f28b8976450 Mon Sep 17 00:00:00 2001 From: JonasHill Date: Fri, 22 Nov 2024 16:01:23 +0100 Subject: [PATCH 2/8] include and typename fixings --- .../structural_mechanics_bossak_scheme.h | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h index 5cf32f582a15..a2fb4ac7ea30 100644 --- a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h +++ b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h @@ -17,6 +17,8 @@ // External includes // Project includes +#include "includes/variables.h" +#include "structural_mechanics_application_variables.h" #include "includes/cfd_variables.h" //TODO: For OSS_SWITCH (think about a better place for this variable) #include "processes/calculate_nodal_area_process.h" @@ -54,6 +56,9 @@ class StructuralMechanicsBossakScheme /// Class type for the scheme using ClassType = StructuralMechanicsBossakScheme; + /// Type for the dofs array within BossakBaseType + using DofsArrayType = typename BossakBaseType::DofsArrayType; + /// Type for the system matrix within BossakBaseType using TSystemMatrixType = typename BossakBaseType::TSystemMatrixType; @@ -121,13 +126,14 @@ class StructuralMechanicsBossakScheme }); } } + void Update( - ModelPart& rModelPart, - DofsArrayType& rDofSet, - TSystemMatrixType& rA, - TSystemVectorType& rDx, - TSystemVectorType& rb - ) override + ModelPart& rModelPart, + DofsArrayType& rDofSet, + TSystemMatrixType& rA, + TSystemVectorType& rDx, + TSystemVectorType& rb + ) override { KRATOS_TRY From 2999b9d92e07d405e44846cad807361981d1bfe9 Mon Sep 17 00:00:00 2001 From: JonasHill Date: Mon, 25 Nov 2024 20:18:15 +0100 Subject: [PATCH 3/8] added Predict method for rotational dofs --- .../structural_mechanics_bossak_scheme.h | 60 +++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h index a2fb4ac7ea30..dc4ac961963f 100644 --- a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h +++ b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h @@ -68,6 +68,9 @@ class StructuralMechanicsBossakScheme /// Pointer type for the BaseType using BaseTypePointer = typename BaseType::Pointer; + /// Component type as 'double' + using ComponentType = double; + ///@} ///@name Life Cycle ///@{ @@ -127,6 +130,63 @@ class StructuralMechanicsBossakScheme } } + void Predict( + ModelPart& rModelPart, + DofsArrayType& rDofSet, + TSystemMatrixType& rA, + TSystemVectorType& rDx, + TSystemVectorType& rb + ) override + { + KRATOS_TRY + + // Call the base Predict method + BossakBaseType::Predict(rModelPart, rDofSet, rA, rDx, rb); + + const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo(); + const double delta_time = r_current_process_info[DELTA_TIME]; + + // Predicting angular time derivatives if they are available (nodally for efficiency) + if(rModelPart.HasNodalSolutionStepVariable(ROTATION) && rModelPart.Nodes().size() > 0) + { + const auto it_node_begin = rModelPart.Nodes().begin(); + + // Getting position + KRATOS_ERROR_IF_NOT(it_node_begin->HasDofFor(ROTATION_X)) << "StructuralMechanicsBossakScheme:: ROTATION is not added" << std::endl; + const int rot_pos = it_node_begin->GetDofPosition(ROTATION_X); + + // Getting dimension + const std::size_t dimension = r_current_process_info.Has(DOMAIN_SIZE) ? r_current_process_info.GetValue(DOMAIN_SIZE) : 3; + + //Auxiliar variable + const std::array*, 3> rot_components = {&ROTATION_X, &ROTATION_Y, &ROTATION_Z}; + + block_for_each(rModelPart.Nodes(), array_1d(), [&](Node& rNode, array_1d& rDeltaRotationTLS){ + //Predicting + const array_1d& r_previous_rotation = rNode.FastGetSolutionStepValue(ROTATION, 1); + const array_1d& r_previous_angular_velocity = rNode.FastGetSolutionStepValue(ANGULAR_VELOCITY, 1); + const array_1d& r_previous_angular_acceleration = rNode.FastGetSolutionStepValue(ANGULAR_ACCELERATION, 1); + array_1d& r_current_rotation = rNode.FastGetSolutionStepValue(ROTATION); + array_1d& r_current_angular_velocity = rNode.FastGetSolutionStepValue(ANGULAR_VELOCITY); + array_1d& r_current_angular_acceleration = rNode.FastGetSolutionStepValue(ANGULAR_ACCELERATION); + + for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) { + if (!rNode.GetDof(*rot_components[i_dim], rot_pos + i_dim).IsFixed()) { + r_current_rotation[i_dim] = r_previous_rotation[i_dim] + delta_time * r_previous_angular_velocity[i_dim] + 0.5 * std::pow(delta_time, 2) * r_previous_angular_acceleration[i_dim]; + } + } + + // Updating + noalias(rDeltaRotationTLS) = r_current_rotation - r_previous_rotation; + UpdateAngularVelocity(r_current_angular_velocity, rDeltaRotationTLS, r_previous_angular_velocity, r_previous_angular_acceleration); + UpdateAngularAcceleration(r_current_angular_acceleration, rDeltaRotationTLS, r_previous_angular_velocity, r_previous_angular_acceleration); + + }); + } + + KRATOS_CATCH("") + } + void Update( ModelPart& rModelPart, DofsArrayType& rDofSet, From cf081cd80f900a341673fb7230a720175d0a867b Mon Sep 17 00:00:00 2001 From: JonasHill Date: Mon, 25 Nov 2024 20:53:58 +0100 Subject: [PATCH 4/8] adopted InitializeSolutionStep for rotational dofs --- .../structural_mechanics_bossak_scheme.h | 50 ++++++++++++++++++- 1 file changed, 48 insertions(+), 2 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h index dc4ac961963f..600b5f737b2a 100644 --- a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h +++ b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h @@ -146,7 +146,7 @@ class StructuralMechanicsBossakScheme const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo(); const double delta_time = r_current_process_info[DELTA_TIME]; - // Predicting angular time derivatives if they are available (nodally for efficiency) + // Updating angular time derivatives if they are available (nodally for efficiency) if(rModelPart.HasNodalSolutionStepVariable(ROTATION) && rModelPart.Nodes().size() > 0) { const auto it_node_begin = rModelPart.Nodes().begin(); @@ -226,15 +226,61 @@ class StructuralMechanicsBossakScheme TSystemVectorType &Dx, TSystemVectorType &b) override { + KRATOS_TRY + // Call the base InitializeSolutionStep method BossakBaseType::InitializeSolutionStep(rModelPart, A, Dx, b); - // Update the NODAL_AREA (note that this strictly required only in the updated Lagrangian case) + // The current process info const auto& r_current_process_info = rModelPart.GetProcessInfo(); + + // Updating angular time derivatives if they are available (nodally for efficiency) + if(rModelPart.HasNodalSolutionStepVariable(ROTATION)) + { + const auto it_node_begin = rModelPart.Nodes().begin(); + + // Getting dimension + const std::size_t dimension = r_current_process_info.Has(DOMAIN_SIZE) ? r_current_process_info.GetValue(DOMAIN_SIZE) : 3; + + // Getting position + const int angular_vel_pos = it_node_begin->HasDofFor(ANGULAR_VELOCITY_X) ? static_cast(it_node_begin->GetDofPosition(ANGULAR_VELOCITY_X)) : -1; + const int angular_acc_pos = it_node_begin->HasDofFor(ANGULAR_ACCELERATION_X) ? static_cast(it_node_begin->GetDofPosition(ANGULAR_ACCELERATION_X)) : -1; + + std::array fixed = {false, false, false}; + const std::array*, 3> rot_components = {&ROTATION_X, &ROTATION_Y, &ROTATION_Z}; + const std::array*, 3> ang_vel_components = {&ANGULAR_VELOCITY_X, &ANGULAR_VELOCITY_Y, &ANGULAR_VELOCITY_Z}; + const std::array*, 3> ang_acc_components = {&ANGULAR_ACCELERATION_X, &ANGULAR_ACCELERATION_Y, &ANGULAR_ACCELERATION_Z}; + + block_for_each(rModelPart.Nodes(), fixed, [&](Node& rNode, std::array& rFixedTLS){ + for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) { + rFixedTLS[i_dim] = false; + } + + if (angular_acc_pos > -1) { + for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim){ + if(rNode.GetDof(*ang_acc_components[i_dim], angular_acc_pos + i_dim).IsFixed()){ + rNode.Fix(*rot_components[i_dim]); + rFixedTLS[i_dim] = true; + } + } + } + if (angular_vel_pos > -1) { + for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim){ + if(rNode.GetDof(*ang_vel_components[i_dim], angular_vel_pos + i_dim).IsFixed() && !rFixedTLS[i_dim]){ + rNode.Fix(*rot_components[i_dim]); + } + } + } + }); + } + + // Update the NODAL_AREA (note that this strictly required only in the updated Lagrangian case) const bool oss_switch = r_current_process_info.Has(OSS_SWITCH) ? r_current_process_info[OSS_SWITCH] : false; if (oss_switch && mUpdateNodalArea) { CalculateNodalAreaProcess(rModelPart).Execute(); } + + KRATOS_CATCH("") } void FinalizeNonLinIteration( From 7396068cc5283334b3fe31470b7234f0f8085ef4 Mon Sep 17 00:00:00 2001 From: JonasHill Date: Tue, 26 Nov 2024 09:48:34 +0100 Subject: [PATCH 5/8] adapted Update method for angular time derivative fixing --- .../structural_mechanics_bossak_scheme.h | 46 ++++++++++++++++--- 1 file changed, 40 insertions(+), 6 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h index 600b5f737b2a..12b53335f25c 100644 --- a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h +++ b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h @@ -135,7 +135,7 @@ class StructuralMechanicsBossakScheme DofsArrayType& rDofSet, TSystemMatrixType& rA, TSystemVectorType& rDx, - TSystemVectorType& rb + TSystemVectorType& rb ) override { KRATOS_TRY @@ -147,21 +147,37 @@ class StructuralMechanicsBossakScheme const double delta_time = r_current_process_info[DELTA_TIME]; // Updating angular time derivatives if they are available (nodally for efficiency) - if(rModelPart.HasNodalSolutionStepVariable(ROTATION) && rModelPart.Nodes().size() > 0) + if(rModelPart.HasNodalSolutionStepVariable(ROTATION)) { const auto it_node_begin = rModelPart.Nodes().begin(); // Getting position KRATOS_ERROR_IF_NOT(it_node_begin->HasDofFor(ROTATION_X)) << "StructuralMechanicsBossakScheme:: ROTATION is not added" << std::endl; const int rot_pos = it_node_begin->GetDofPosition(ROTATION_X); + const int ang_vel_pos = it_node_begin->HasDofFor(ANGULAR_VELOCITY_X) ? static_cast(it_node_begin->GetDofPosition(ANGULAR_VELOCITY_X)) : -1; + const int ang_acc_pos = it_node_begin->HasDofFor(ANGULAR_ACCELERATION_X) ? static_cast(it_node_begin->GetDofPosition(ANGULAR_ACCELERATION_X)) : -1; // Getting dimension const std::size_t dimension = r_current_process_info.Has(DOMAIN_SIZE) ? r_current_process_info.GetValue(DOMAIN_SIZE) : 3; //Auxiliar variable + array_1d delta_rotation; + std::array predicted = {false, false, false}; const std::array*, 3> rot_components = {&ROTATION_X, &ROTATION_Y, &ROTATION_Z}; + const std::array*, 3> ang_vel_components = {&ANGULAR_VELOCITY_X, &ANGULAR_VELOCITY_Y, &ANGULAR_VELOCITY_Z}; + const std::array*, 3> ang_acc_components = {&ANGULAR_ACCELERATION_X, &ANGULAR_ACCELERATION_Y, &ANGULAR_ACCELERATION_Z}; + + typedef std::tuple, std::array> TLSContainerType; + TLSContainerType aux_TLS = std::make_tuple(delta_rotation, predicted); + + block_for_each(rModelPart.Nodes(), aux_TLS, [&](Node& rNode, TLSContainerType& rAuxTLS){ + auto& r_delta_rotation = std::get<0>(rAuxTLS); + auto& r_predicted = std::get<1>(rAuxTLS); + + for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) { + r_predicted[i_dim] = false; + } - block_for_each(rModelPart.Nodes(), array_1d(), [&](Node& rNode, array_1d& rDeltaRotationTLS){ //Predicting const array_1d& r_previous_rotation = rNode.FastGetSolutionStepValue(ROTATION, 1); const array_1d& r_previous_angular_velocity = rNode.FastGetSolutionStepValue(ANGULAR_VELOCITY, 1); @@ -170,6 +186,24 @@ class StructuralMechanicsBossakScheme array_1d& r_current_angular_velocity = rNode.FastGetSolutionStepValue(ANGULAR_VELOCITY); array_1d& r_current_angular_acceleration = rNode.FastGetSolutionStepValue(ANGULAR_ACCELERATION); + if(ang_acc_pos > -1) { + for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) { + if (!rNode.GetDof(*ang_acc_components[i_dim], ang_acc_pos + i_dim).IsFixed()) { + r_delta_rotation[i_dim] = (r_current_angular_acceleration[i_dim] + this->mBossak.c3 * r_previous_angular_acceleration[i_dim] + this->mBossak.c2 * r_previous_angular_velocity[i_dim])/(this->mBossak.c0); + r_current_rotation[i_dim] = r_previous_rotation[i_dim] + r_delta_rotation[i_dim]; + r_predicted[i_dim] = true; + } + } + } + if(ang_vel_pos > -1) { + for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) { + if (!rNode.GetDof(*ang_vel_components[i_dim], ang_vel_pos + i_dim).IsFixed() && !r_predicted[i_dim]) { + r_delta_rotation[i_dim] = (r_current_angular_velocity[i_dim] + this->mBossak.c4 * r_previous_angular_velocity[i_dim] + this->mBossak.c5 * r_previous_angular_acceleration[i_dim])/(this->mBossak.c1); + r_current_rotation[i_dim] = r_previous_rotation[i_dim] + r_delta_rotation[i_dim]; + r_predicted[i_dim] = true; + } + } + } for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) { if (!rNode.GetDof(*rot_components[i_dim], rot_pos + i_dim).IsFixed()) { r_current_rotation[i_dim] = r_previous_rotation[i_dim] + delta_time * r_previous_angular_velocity[i_dim] + 0.5 * std::pow(delta_time, 2) * r_previous_angular_acceleration[i_dim]; @@ -177,9 +211,9 @@ class StructuralMechanicsBossakScheme } // Updating - noalias(rDeltaRotationTLS) = r_current_rotation - r_previous_rotation; - UpdateAngularVelocity(r_current_angular_velocity, rDeltaRotationTLS, r_previous_angular_velocity, r_previous_angular_acceleration); - UpdateAngularAcceleration(r_current_angular_acceleration, rDeltaRotationTLS, r_previous_angular_velocity, r_previous_angular_acceleration); + noalias(r_delta_rotation) = r_current_rotation - r_previous_rotation; + UpdateAngularVelocity(r_current_angular_velocity, r_delta_rotation, r_previous_angular_velocity, r_previous_angular_acceleration); + UpdateAngularAcceleration(r_current_angular_acceleration, r_delta_rotation, r_previous_angular_velocity, r_previous_angular_acceleration); }); } From ac24be84845dcdf07f9216c20b1ba94da64747d2 Mon Sep 17 00:00:00 2001 From: JonasHill Date: Tue, 26 Nov 2024 09:57:16 +0100 Subject: [PATCH 6/8] bugfix --- .../custom_schemes/structural_mechanics_bossak_scheme.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h index 12b53335f25c..5b1da57d0fff 100644 --- a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h +++ b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h @@ -188,18 +188,20 @@ class StructuralMechanicsBossakScheme if(ang_acc_pos > -1) { for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) { - if (!rNode.GetDof(*ang_acc_components[i_dim], ang_acc_pos + i_dim).IsFixed()) { + if (rNode.GetDof(*ang_acc_components[i_dim], ang_acc_pos + i_dim).IsFixed()) { r_delta_rotation[i_dim] = (r_current_angular_acceleration[i_dim] + this->mBossak.c3 * r_previous_angular_acceleration[i_dim] + this->mBossak.c2 * r_previous_angular_velocity[i_dim])/(this->mBossak.c0); r_current_rotation[i_dim] = r_previous_rotation[i_dim] + r_delta_rotation[i_dim]; + std::cout << "Updating rotation by angular acceleration" << std::endl; r_predicted[i_dim] = true; } } } if(ang_vel_pos > -1) { for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) { - if (!rNode.GetDof(*ang_vel_components[i_dim], ang_vel_pos + i_dim).IsFixed() && !r_predicted[i_dim]) { + if (rNode.GetDof(*ang_vel_components[i_dim], ang_vel_pos + i_dim).IsFixed() && !r_predicted[i_dim]) { r_delta_rotation[i_dim] = (r_current_angular_velocity[i_dim] + this->mBossak.c4 * r_previous_angular_velocity[i_dim] + this->mBossak.c5 * r_previous_angular_acceleration[i_dim])/(this->mBossak.c1); r_current_rotation[i_dim] = r_previous_rotation[i_dim] + r_delta_rotation[i_dim]; + std::cout << "Updating rotation by angular velocity" << std::endl; r_predicted[i_dim] = true; } } From dcb7c31de2e7d3fdbc3a234da56291a24e0e2bb7 Mon Sep 17 00:00:00 2001 From: JonasHill Date: Tue, 26 Nov 2024 12:37:02 +0100 Subject: [PATCH 7/8] deleted std::cout message --- .../custom_schemes/structural_mechanics_bossak_scheme.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h index 5b1da57d0fff..b49cd9fe8088 100644 --- a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h +++ b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h @@ -191,7 +191,6 @@ class StructuralMechanicsBossakScheme if (rNode.GetDof(*ang_acc_components[i_dim], ang_acc_pos + i_dim).IsFixed()) { r_delta_rotation[i_dim] = (r_current_angular_acceleration[i_dim] + this->mBossak.c3 * r_previous_angular_acceleration[i_dim] + this->mBossak.c2 * r_previous_angular_velocity[i_dim])/(this->mBossak.c0); r_current_rotation[i_dim] = r_previous_rotation[i_dim] + r_delta_rotation[i_dim]; - std::cout << "Updating rotation by angular acceleration" << std::endl; r_predicted[i_dim] = true; } } @@ -201,7 +200,6 @@ class StructuralMechanicsBossakScheme if (rNode.GetDof(*ang_vel_components[i_dim], ang_vel_pos + i_dim).IsFixed() && !r_predicted[i_dim]) { r_delta_rotation[i_dim] = (r_current_angular_velocity[i_dim] + this->mBossak.c4 * r_previous_angular_velocity[i_dim] + this->mBossak.c5 * r_previous_angular_acceleration[i_dim])/(this->mBossak.c1); r_current_rotation[i_dim] = r_previous_rotation[i_dim] + r_delta_rotation[i_dim]; - std::cout << "Updating rotation by angular velocity" << std::endl; r_predicted[i_dim] = true; } } From 54b7fa8176f870075a71eb7a6edf0ca64232cadf Mon Sep 17 00:00:00 2001 From: JonasHill Date: Tue, 26 Nov 2024 13:21:22 +0100 Subject: [PATCH 8/8] added doxygen for Predict and Update methods --- .../structural_mechanics_bossak_scheme.h | 26 ++++++++++++++----- 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h index b49cd9fe8088..595b2791a7f1 100644 --- a/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h +++ b/applications/StructuralMechanicsApplication/custom_strategies/custom_schemes/structural_mechanics_bossak_scheme.h @@ -130,6 +130,15 @@ class StructuralMechanicsBossakScheme } } + /** + * @brief Apply the predictor. + * @details @f[ x_{k+1} = x_{k} + \dot x_{k} \Delta t + \frac{1}{2} \ddot x_{k} \Delta t^2 @f] + * @param rModelPart @ref ModelPart to update. + * @param rDofSet Set of all primary variables. + * @param rA Left hand side matrix. + * @param rDx Primary variable updates. + * @param rb Right hand side vector. + */ void Predict( ModelPart& rModelPart, DofsArrayType& rDofSet, @@ -147,8 +156,7 @@ class StructuralMechanicsBossakScheme const double delta_time = r_current_process_info[DELTA_TIME]; // Updating angular time derivatives if they are available (nodally for efficiency) - if(rModelPart.HasNodalSolutionStepVariable(ROTATION)) - { + if(rModelPart.HasNodalSolutionStepVariable(ROTATION)){ const auto it_node_begin = rModelPart.Nodes().begin(); // Getting position @@ -221,6 +229,14 @@ class StructuralMechanicsBossakScheme KRATOS_CATCH("") } + /** @brief Update state variables and its time derivatives. + * @details @f[ x_{n+1}^{k+1} = x_{n+1}^k+ \Delta x @f] + * @param rModelPart @ref ModelPart to update. + * @param rDofSet Set of all primary variables. + * @param rA Left hand side matrix. + * @param rDx Primary variable updates. + * @param rb Right hand side vector. + */ void Update( ModelPart& rModelPart, DofsArrayType& rDofSet, @@ -235,8 +251,7 @@ class StructuralMechanicsBossakScheme BossakBaseType::Update(rModelPart, rDofSet, rA, rDx, rb); // Updating angular time derivatives if they are available (nodally for efficiency) - if(rModelPart.HasNodalSolutionStepVariable(ROTATION)) - { + if(rModelPart.HasNodalSolutionStepVariable(ROTATION)){ block_for_each(rModelPart.Nodes(), array_1d(), [&](Node& rNode, array_1d& rDeltaRotationTLS){ noalias(rDeltaRotationTLS) = rNode.FastGetSolutionStepValue(ROTATION) - rNode.FastGetSolutionStepValue(ROTATION, 1); @@ -269,8 +284,7 @@ class StructuralMechanicsBossakScheme const auto& r_current_process_info = rModelPart.GetProcessInfo(); // Updating angular time derivatives if they are available (nodally for efficiency) - if(rModelPart.HasNodalSolutionStepVariable(ROTATION)) - { + if(rModelPart.HasNodalSolutionStepVariable(ROTATION)){ const auto it_node_begin = rModelPart.Nodes().begin(); // Getting dimension