Skip to content

Commit

Permalink
[GeoMechanicsApplication] Fix ordering issue for reset displacement (#…
Browse files Browse the repository at this point in the history 
…12845)

Only reset variables if their dofs are not fixed and do the reset after the ExecuteInitialize of processes
  • Loading branch information
@rfaasse
rfaasse authored Nov 14, 2024
1 parent cab1cb3 commit fb91210
Show file tree
Hide file tree
Showing 11 changed files with 65 additions and 47 deletions.
2 changes: 1 addition & 1 deletion applications/GeoMechanicsApplication/custom_elements/compressibility_calculator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Matrix, Vector> CompressibilityCalculator::LocalSystemContribution()
Expand Down
10 changes: 5 additions & 5 deletions applications/GeoMechanicsApplication/custom_elements/permeability_calculator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down
6 changes: 3 additions & 3 deletions applications/GeoMechanicsApplication/custom_elements/permeability_calculator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ class PermeabilityCalculator : public ContributionCalculator
class InputProvider
{
public:
InputProvider(std::function<const Properties&()> GetElementProperties,
InputProvider(std::function<const Properties&()> GetElementProperties,
std::function<const std::vector<RetentionLaw::Pointer>&()> GetRetentionLaws,
std::function<Vector()> GetIntegrationCoefficients,
std::function<Vector(const Variable<double>&)> GetNodalValuesOf,
Expand All @@ -39,9 +39,9 @@ class PermeabilityCalculator : public ContributionCalculator
{
}

[[nodiscard]] const Properties& GetElementProperties() const;
[[nodiscard]] const Properties& GetElementProperties() const;
[[nodiscard]] const std::vector<RetentionLaw::Pointer>& GetRetentionLaws() const;
[[nodiscard]] Vector GetIntegrationCoefficients() const;
[[nodiscard]] Vector GetIntegrationCoefficients() const;
[[nodiscard]] Vector GetNodalValues(const Variable<double>& rVariable) const;
[[nodiscard]] Geometry<Node>::ShapeFunctionsGradientsType GetShapeFunctionGradients() const;

Expand Down
12 changes: 7 additions & 5 deletions applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_<NodeUtilities>(rModule, "NodeUtilities")
.def("AssignUpdatedVectorVariableToNonFixedComponentsOfNodes",
&NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes);
}

} // Namespace Kratos::Python.
2 changes: 1 addition & 1 deletion applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
namespace Kratos::Python
{

void AddCustomUtilitiesToPython(pybind11::module&);
void AddCustomUtilitiesToPython(const pybind11::module& rModule);

} // namespace Kratos::Python.

Expand Down
16 changes: 14 additions & 2 deletions applications/GeoMechanicsApplication/custom_utilities/node_utilities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@ namespace Kratos

void NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode,
const Variable<array_1d<double, 3>>& rDestinationVariable,
const array_1d<double, 3>& rNewValues)
const array_1d<double, 3>& rNewValues,
IndexType SolutionStepIndex)
{
const std::vector<std::string> components = {"X", "Y", "Z"};
for (const auto& zipped : boost::combine(rNewValues, components)) {
Expand All @@ -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<array_1d<double, 3>>& rDestinationVariable,
const array_1d<double, 3>& rNewValues,
IndexType SolutionStepIndex)
{
block_for_each(rNodes, [&rDestinationVariable, &rNewValues, SolutionStepIndex](Node& rNode) {
AssignUpdatedVectorVariableToNonFixedComponents(rNode, rDestinationVariable, rNewValues, SolutionStepIndex);
});
}

} // namespace Kratos
11 changes: 10 additions & 1 deletion applications/GeoMechanicsApplication/custom_utilities/node_utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

#include "containers/array_1d.h"
#include "containers/variable.h"
#include "includes/model_part.h"

namespace Kratos
{
Expand All @@ -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<array_1d<double, 3>>& rDestinationVariable,
const array_1d<double, 3>& rNewValues);
const array_1d<double, 3>& rNewValues,
IndexType SolutionStepIndex = 0);

static void AssignUpdatedVectorVariableToNonFixedComponentsOfNodes(const ModelPart::NodesContainerType& rNodes,
const Variable<array_1d<double, 3>>& rDestinationVariable,
const array_1d<double, 3>& rNewValues,
IndexType SolutionStepIndex = 0);
};

} // namespace Kratos
20 changes: 8 additions & 12 deletions applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,14 +191,6 @@ int KratosGeoSettlement::RunStage(const std::filesystem::path& rWorki
std::vector<std::shared_ptr<Process>> processes = GetProcesses(project_parameters);
std::vector<std::weak_ptr<Process>> 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();
}
Expand All @@ -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;
Expand Down Expand Up @@ -341,10 +339,8 @@ std::shared_ptr<StrategyWrapper> 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());
}
Expand Down
14 changes: 6 additions & 8 deletions applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -73,17 +74,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoSettlement
const std::stringstream& rKratosLogBuffer) const;

template <typename TVariableType>
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 <typename ProcessType>
Expand Down
17 changes: 9 additions & 8 deletions applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion ...GeoMechanicsApplication/tests/truss_with_reset_displacement/ProjectParameters_stage2.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down

0 comments on commit fb91210

Please sign in to comment.