diff --git a/applications/GeoMechanicsApplication/custom_elements/compressibility_calculator.cpp b/applications/GeoMechanicsApplication/custom_elements/compressibility_calculator.cpp index 6a865db70e6a..394ef8313975 100644 --- a/applications/GeoMechanicsApplication/custom_elements/compressibility_calculator.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/compressibility_calculator.cpp @@ -38,7 +38,7 @@ Vector CompressibilityCalculator::RHSContribution(const Matrix& rCompressibility Matrix CompressibilityCalculator::LHSContribution(const Matrix& rCompressibilityMatrix) const { - return mInputProvider.GetMatrixScalarFactor()* rCompressibilityMatrix; + return mInputProvider.GetMatrixScalarFactor() * rCompressibilityMatrix; } std::pair CompressibilityCalculator::LocalSystemContribution() diff --git a/applications/GeoMechanicsApplication/custom_elements/permeability_calculator.cpp b/applications/GeoMechanicsApplication/custom_elements/permeability_calculator.cpp index 70c4452b6c34..64c698360706 100644 --- a/applications/GeoMechanicsApplication/custom_elements/permeability_calculator.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/permeability_calculator.cpp @@ -44,11 +44,11 @@ Vector PermeabilityCalculator::RHSContribution(const Matrix& rPermeabilityMatrix Matrix PermeabilityCalculator::CalculatePermeabilityMatrix() const { RetentionLaw::Parameters retention_parameters(mInputProvider.GetElementProperties()); - const auto& r_properties = mInputProvider.GetElementProperties(); - auto integration_coefficients = mInputProvider.GetIntegrationCoefficients(); - const auto shape_function_gradients = mInputProvider.GetShapeFunctionGradients(); - const auto local_dimension = shape_function_gradients[0].size2(); - const Matrix constitutive_matrix = + const auto& r_properties = mInputProvider.GetElementProperties(); + auto integration_coefficients = mInputProvider.GetIntegrationCoefficients(); + const auto shape_function_gradients = mInputProvider.GetShapeFunctionGradients(); + const auto local_dimension = shape_function_gradients[0].size2(); + const Matrix constitutive_matrix = GeoElementUtilities::FillPermeabilityMatrix(r_properties, local_dimension); const auto number_of_nodes = shape_function_gradients[0].size1(); diff --git a/applications/GeoMechanicsApplication/custom_elements/permeability_calculator.h b/applications/GeoMechanicsApplication/custom_elements/permeability_calculator.h index 7f3329ead9db..86ffbee68004 100644 --- a/applications/GeoMechanicsApplication/custom_elements/permeability_calculator.h +++ b/applications/GeoMechanicsApplication/custom_elements/permeability_calculator.h @@ -26,7 +26,7 @@ class PermeabilityCalculator : public ContributionCalculator class InputProvider { public: - InputProvider(std::function GetElementProperties, + InputProvider(std::function GetElementProperties, std::function&()> GetRetentionLaws, std::function GetIntegrationCoefficients, std::function&)> GetNodalValuesOf, @@ -39,9 +39,9 @@ class PermeabilityCalculator : public ContributionCalculator { } - [[nodiscard]] const Properties& GetElementProperties() const; + [[nodiscard]] const Properties& GetElementProperties() const; [[nodiscard]] const std::vector& GetRetentionLaws() const; - [[nodiscard]] Vector GetIntegrationCoefficients() const; + [[nodiscard]] Vector GetIntegrationCoefficients() const; [[nodiscard]] Vector GetNodalValues(const Variable& rVariable) const; [[nodiscard]] Geometry::ShapeFunctionsGradientsType GetShapeFunctionGradients() const; diff --git a/applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.cpp index 69124591620b..2ce29608bcb3 100644 --- a/applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.cpp @@ -15,15 +15,17 @@ // Project includes #include "custom_python/add_custom_utilities_to_python.h" -#include "includes/kratos_parameters.h" -#include "custom_utilities/condition_utilities.hpp" -#include "custom_utilities/element_utilities.hpp" -#include "custom_utilities/interface_element_utilities.h" +#include "custom_utilities/node_utilities.h" namespace Kratos::Python { -void AddCustomUtilitiesToPython(pybind11::module&) { namespace py = pybind11; } +void AddCustomUtilitiesToPython(const pybind11::module& rModule) +{ + pybind11::class_(rModule, "NodeUtilities") + .def("AssignUpdatedVectorVariableToNonFixedComponentsOfNodes", + &NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes); +} } // Namespace Kratos::Python. diff --git a/applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.h b/applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.h index 3a3b6f8df4e9..dd9edb7f77c6 100644 --- a/applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.h +++ b/applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.h @@ -25,7 +25,7 @@ namespace Kratos::Python { -void AddCustomUtilitiesToPython(pybind11::module&); +void AddCustomUtilitiesToPython(const pybind11::module& rModule); } // namespace Kratos::Python. diff --git a/applications/GeoMechanicsApplication/custom_utilities/node_utilities.cpp b/applications/GeoMechanicsApplication/custom_utilities/node_utilities.cpp index e98b14d2a230..7540210fc910 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/node_utilities.cpp +++ b/applications/GeoMechanicsApplication/custom_utilities/node_utilities.cpp @@ -22,7 +22,8 @@ namespace Kratos void NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode, const Variable>& rDestinationVariable, - const array_1d& rNewValues) + const array_1d& rNewValues, + IndexType SolutionStepIndex) { const std::vector components = {"X", "Y", "Z"}; for (const auto& zipped : boost::combine(rNewValues, components)) { @@ -33,9 +34,20 @@ void NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode, if (const auto& component_variable = VariablesUtilities::GetComponentFromVectorVariable(rDestinationVariable.Name(), component); !rNode.IsFixed(component_variable)) { - rNode.FastGetSolutionStepValue(component_variable, 0) = new_value; + rNode.FastGetSolutionStepValue(component_variable, SolutionStepIndex) = new_value; } } } +void NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes( + const ModelPart::NodesContainerType& rNodes, + const Variable>& rDestinationVariable, + const array_1d& rNewValues, + IndexType SolutionStepIndex) +{ + block_for_each(rNodes, [&rDestinationVariable, &rNewValues, SolutionStepIndex](Node& rNode) { + AssignUpdatedVectorVariableToNonFixedComponents(rNode, rDestinationVariable, rNewValues, SolutionStepIndex); + }); +} + } // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_utilities/node_utilities.h b/applications/GeoMechanicsApplication/custom_utilities/node_utilities.h index e0ddc41e6c4f..bcd77544f50e 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/node_utilities.h +++ b/applications/GeoMechanicsApplication/custom_utilities/node_utilities.h @@ -14,6 +14,7 @@ #include "containers/array_1d.h" #include "containers/variable.h" +#include "includes/model_part.h" namespace Kratos { @@ -23,9 +24,17 @@ class Node; class KRATOS_API(GEO_MECHANICS_APPLICATION) NodeUtilities { public: + KRATOS_CLASS_POINTER_DEFINITION(NodeUtilities); + static void AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode, const Variable>& rDestinationVariable, - const array_1d& rNewValues); + const array_1d& rNewValues, + IndexType SolutionStepIndex = 0); + + static void AssignUpdatedVectorVariableToNonFixedComponentsOfNodes(const ModelPart::NodesContainerType& rNodes, + const Variable>& rDestinationVariable, + const array_1d& rNewValues, + IndexType SolutionStepIndex = 0); }; } // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp index cae242401994..4447f27177fc 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp @@ -191,14 +191,6 @@ int KratosGeoSettlement::RunStage(const std::filesystem::path& rWorki std::vector> processes = GetProcesses(project_parameters); std::vector> process_observables(processes.begin(), processes.end()); - if (mpTimeLoopExecutor) { - mpTimeLoopExecutor->SetCancelDelegate(rShouldCancel); - mpTimeLoopExecutor->SetProgressDelegate(rProgressDelegate); - mpTimeLoopExecutor->SetProcessObservables(process_observables); - mpTimeLoopExecutor->SetTimeIncrementor(MakeTimeIncrementor(project_parameters)); - mpTimeLoopExecutor->SetSolverStrategyWrapper(MakeStrategyWrapper(project_parameters, rWorkingDirectory)); - } - for (const auto& process : processes) { process->ExecuteInitialize(); } @@ -208,6 +200,12 @@ int KratosGeoSettlement::RunStage(const std::filesystem::path& rWorki } if (mpTimeLoopExecutor) { + mpTimeLoopExecutor->SetCancelDelegate(rShouldCancel); + mpTimeLoopExecutor->SetProgressDelegate(rProgressDelegate); + mpTimeLoopExecutor->SetProcessObservables(process_observables); + mpTimeLoopExecutor->SetTimeIncrementor(MakeTimeIncrementor(project_parameters)); + mpTimeLoopExecutor->SetSolverStrategyWrapper(MakeStrategyWrapper(project_parameters, rWorkingDirectory)); + // For now, pass a dummy state. THIS PROBABLY NEEDS TO BE REFINED AT SOME POINT! TimeStepEndState start_of_loop_state; start_of_loop_state.convergence_state = TimeStepEndState::ConvergenceState::converged; @@ -341,10 +339,8 @@ std::shared_ptr KratosGeoSettlement::MakeStrategyWrapper(const GetMainModelPart().CloneTimeStep(); if (rProjectParameters["solver_settings"]["reset_displacements"].GetBool()) { - constexpr auto source_index = std::size_t{0}; - constexpr auto destination_index = std::size_t{1}; - RestoreValuesOfNodalVariable(DISPLACEMENT, source_index, destination_index); - RestoreValuesOfNodalVariable(ROTATION, source_index, destination_index); + ResetValuesOfNodalVariable(DISPLACEMENT); + ResetValuesOfNodalVariable(ROTATION); VariableUtils{}.UpdateCurrentToInitialConfiguration(GetComputationalModelPart().Nodes()); } diff --git a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.h b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.h index b3fdf6e24c44..42c84579fc4c 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.h +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.h @@ -19,6 +19,7 @@ #include "includes/kernel.h" #include "includes/kratos_export_api.h" +#include "custom_utilities/node_utilities.h" #include "custom_utilities/process_factory.hpp" #include "geo_mechanics_application.h" #include "linear_solvers_application.h" @@ -73,17 +74,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoSettlement const std::stringstream& rKratosLogBuffer) const; template - void RestoreValuesOfNodalVariable(const TVariableType& rVariable, Node::IndexType SourceIndex, Node::IndexType DestinationIndex) + void ResetValuesOfNodalVariable(const TVariableType& rVariable) { if (!GetComputationalModelPart().HasNodalSolutionStepVariable(rVariable)) return; - VariableUtils{}.SetHistoricalVariableToZero(rVariable, GetComputationalModelPart().Nodes()); - - block_for_each(GetComputationalModelPart().Nodes(), - [&rVariable, SourceIndex, DestinationIndex](auto& node) { - node.GetSolutionStepValue(rVariable, DestinationIndex) = - node.GetSolutionStepValue(rVariable, SourceIndex); - }); + NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes( + GetComputationalModelPart().Nodes(), rVariable, rVariable.Zero(), 0); + NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes( + GetComputationalModelPart().Nodes(), rVariable, rVariable.Zero(), 1); } template diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py index d187b65fa541..8e872345545a 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py @@ -57,14 +57,15 @@ def _CalculateTotalDisplacement(self,node, old_total_displacement): def ResetIfHasNodalSolutionStepVariable(self, variable): if self._GetSolver().main_model_part.HasNodalSolutionStepVariable(variable): - KratosMultiphysics.VariableUtils().SetHistoricalVariableToZero(variable, self._GetSolver().GetComputingModelPart().Nodes) - for node in self._GetSolver().GetComputingModelPart().Nodes: - new_value = node.GetSolutionStepValue(variable, 0) - node.SetSolutionStepValue(variable, 1, new_value) - - def ModifyInitialGeometry(self): - # Overrides the base class. Necessary to let reset_displacements function correctly i.c.w. prescribed displacements/rotations. - # The reset needs to take place befor the Initialize of the processes, as these will set the Dirichlet condition. + zero_vector = Kratos.Array3([0.0, 0.0, 0.0]) + KratosGeo.NodeUtilities.AssignUpdatedVectorVariableToNonFixedComponentsOfNodes( + self._GetSolver().GetComputingModelPart().Nodes, variable, zero_vector, 0) + KratosGeo.NodeUtilities.AssignUpdatedVectorVariableToNonFixedComponentsOfNodes( + self._GetSolver().GetComputingModelPart().Nodes, variable, zero_vector, 1) + + def Initialize(self): + super().Initialize() + self._GetSolver().main_model_part.ProcessInfo[KratosGeo.RESET_DISPLACEMENTS] = self.reset_displacements if self.reset_displacements: self.ResetIfHasNodalSolutionStepVariable(KratosMultiphysics.DISPLACEMENT) diff --git a/applications/GeoMechanicsApplication/tests/truss_with_reset_displacement/ProjectParameters_stage2.json b/applications/GeoMechanicsApplication/tests/truss_with_reset_displacement/ProjectParameters_stage2.json index 42c9b17a1037..f370677acea1 100644 --- a/applications/GeoMechanicsApplication/tests/truss_with_reset_displacement/ProjectParameters_stage2.json +++ b/applications/GeoMechanicsApplication/tests/truss_with_reset_displacement/ProjectParameters_stage2.json @@ -30,7 +30,7 @@ "block_builder": true, "solution_type": "Quasi-Static", "scheme_type": "Newmark", - "reset_displacements": false, + "reset_displacements": true, "newmark_beta": 0.25, "newmark_gamma": 0.5, "newmark_theta": 0.5,