Skip to content

[GeoMechanicsApplication] Improvements for the nodal extrapolation process #14068

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We鈥檒l occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
avdg81 merged 3 commits into from
Dec 18, 2025
Merged

[GeoMechanicsApplication] Improvements for the nodal extrapolation process #14068

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b578a17
Removed override of member `Execute`
avdg81 Dec 16, 2025
8b5f254
Removed an erroneously added inner loop over model parts
avdg81 Dec 16, 2025
2d489ac
Removed calls to `ExecuteFinalize` from the unit tests
avdg81 Dec 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 11 additions & 16 deletions ...pplication/custom_processes/geo_extrapolate_integration_point_values_to_nodes_process.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,12 +60,6 @@ void GeoExtrapolateIntegrationPointValuesToNodesProcess::FillVariableLists(const
}
}

void GeoExtrapolateIntegrationPointValuesToNodesProcess::Execute()
{
ExecuteBeforeSolutionLoop();
ExecuteFinalizeSolutionStep();
}

void GeoExtrapolateIntegrationPointValuesToNodesProcess::ExecuteBeforeSolutionLoop()
{
InitializeAverageVariablesForElements();
Expand Down Expand Up @@ -132,10 +126,11 @@ void GeoExtrapolateIntegrationPointValuesToNodesProcess::ExecuteFinalizeSolution
CacheExtrapolationMatricesForElements();

for (const auto& r_model_part : mrModelParts) {
block_for_each(r_model_part.get().Elements(), [this](Element& rElement) {
block_for_each(r_model_part.get().Elements(),
[this, &r_process_info = r_model_part.get().GetProcessInfo()](Element& rElement) {
if (rElement.IsActive()) {
AddIntegrationPointContributionsForAllVariables(
rElement, GetCachedExtrapolationMatrixFor(rElement));
rElement, GetCachedExtrapolationMatrixFor(rElement), r_process_info);
}
});
}
Expand Down Expand Up @@ -173,25 +168,25 @@ const Matrix& GeoExtrapolateIntegrationPointValuesToNodesProcess::GetCachedExtra
}

void GeoExtrapolateIntegrationPointValuesToNodesProcess::AddIntegrationPointContributionsForAllVariables(
Element& rElement, const Matrix& rExtrapolationMatrix) const
Element& rElement, const Matrix& rExtrapolationMatrix, const ProcessInfo& rProcessInfo) const
{
const auto integration_points_number = GeoElementUtilities::GetNumberOfIntegrationPointsOf(rElement);

for (const auto p_var : mDoubleVariables) {
AddIntegrationContributionsToNodes(rElement, *p_var, rExtrapolationMatrix,
integration_points_number, AtomicAdd<double>);
integration_points_number, rProcessInfo, AtomicAdd<double>);
}
for (const auto p_var : mArrayVariables) {
AddIntegrationContributionsToNodes(rElement, *p_var, rExtrapolationMatrix,
integration_points_number, AtomicAdd<double, 3>);
AddIntegrationContributionsToNodes(rElement, *p_var, rExtrapolationMatrix, integration_points_number,
rProcessInfo, AtomicAdd<double, 3>);
}
for (const auto p_var : mVectorVariables) {
AddIntegrationContributionsToNodes(rElement, *p_var, rExtrapolationMatrix,
integration_points_number, AtomicAddVector<Vector, Vector>);
AddIntegrationContributionsToNodes(rElement, *p_var, rExtrapolationMatrix, integration_points_number,
rProcessInfo, AtomicAddVector<Vector, Vector>);
}
for (const auto p_var : mMatrixVariables) {
AddIntegrationContributionsToNodes(rElement, *p_var, rExtrapolationMatrix,
integration_points_number, AtomicAddMatrix<Matrix, Matrix>);
AddIntegrationContributionsToNodes(rElement, *p_var, rExtrapolationMatrix, integration_points_number,
rProcessInfo, AtomicAddMatrix<Matrix, Matrix>);
}
}

Expand Down
33 changes: 15 additions & 18 deletions ...sApplication/custom_processes/geo_extrapolate_integration_point_values_to_nodes_process.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoExtrapolateIntegrationPointValues

~GeoExtrapolateIntegrationPointValuesToNodesProcess() override;

void Execute() override;
void ExecuteBeforeSolutionLoop() override;
void ExecuteFinalizeSolutionStep() override;
void ExecuteFinalize() override;
Expand Down Expand Up @@ -96,27 +95,23 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoExtrapolateIntegrationPointValues
const Variable<T>& rVariable,
const Matrix& rExtrapolationMatrix,
SizeType NumberOfIntegrationPoints,
const ProcessInfo& rProcessInfo,
const U& rAtomicAddOperation) const
{
auto& r_geometry = rElement.GetGeometry();
std::vector<T> values_on_integration_points(NumberOfIntegrationPoints);
for (const auto& r_model_part : mrModelParts) {
rElement.CalculateOnIntegrationPoints(rVariable, values_on_integration_points,
r_model_part.get().GetProcessInfo());
Comment thread
markelov208 marked this conversation as resolved.

for (IndexType iNode = 0; iNode < r_geometry.PointsNumber(); ++iNode) {
// We first initialize the source, which we need to do by getting the first value,
// because we don't know the size of the dynamically allocated Vector/Matrix
T source = rExtrapolationMatrix(iNode, 0) * values_on_integration_points[0];
for (IndexType i_gauss_point = 1;
i_gauss_point < values_on_integration_points.size(); ++i_gauss_point) {
source += rExtrapolationMatrix(iNode, i_gauss_point) *
values_on_integration_points[i_gauss_point];
}
source /= r_geometry[iNode].GetValue(mrAverageVariable);

rAtomicAddOperation(r_geometry[iNode].FastGetSolutionStepValue(rVariable), source);
rElement.CalculateOnIntegrationPoints(rVariable, values_on_integration_points, rProcessInfo);

for (IndexType iNode = 0; iNode < r_geometry.PointsNumber(); ++iNode) {
// We first initialize the source, which we need to do by getting the first value,
// because we don't know the size of the dynamically allocated Vector/Matrix
T source = rExtrapolationMatrix(iNode, 0) * values_on_integration_points[0];
for (IndexType i_gauss_point = 1; i_gauss_point < values_on_integration_points.size(); ++i_gauss_point) {
source += rExtrapolationMatrix(iNode, i_gauss_point) * values_on_integration_points[i_gauss_point];
}
source /= r_geometry[iNode].GetValue(mrAverageVariable);

rAtomicAddOperation(r_geometry[iNode].FastGetSolutionStepValue(rVariable), source);
}
}

Expand All @@ -125,7 +120,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoExtrapolateIntegrationPointValues
[[nodiscard]] bool ExtrapolationMatrixIsCachedFor(const Element& rElement) const;
[[nodiscard]] const Matrix& GetCachedExtrapolationMatrixFor(const Element& rElement) const;

void AddIntegrationPointContributionsForAllVariables(Element& rElement, const Matrix& rExtrapolationMatrix) const;
void AddIntegrationPointContributionsForAllVariables(Element& rElement,
const Matrix& rExtrapolationMatrix,
const ProcessInfo& rProcessInfo) const;
};

} // namespace Kratos.
60 changes: 54 additions & 6 deletions ...tests/custom_processes/test_geo_extrapolate_integration_point_values_to_nodes_process.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,10 @@
#include "tests/cpp_tests/geo_mechanics_fast_suite.h"

#include <numbers>
#include <string>

using namespace std::numbers;
using namespace std::string_literals;

namespace Kratos::Testing
{
Expand Down Expand Up @@ -155,7 +157,8 @@ KRATOS_TEST_CASE_IN_SUITE(TestExtrapolationProcess_ExtrapolatesCorrectlyForConst
})");

GeoExtrapolateIntegrationPointValuesToNodesProcess process(model, parameters);
process.Execute();
process.ExecuteBeforeSolutionLoop();
process.ExecuteFinalizeSolutionStep();

for (auto& node : model_part.Nodes()) {
KRATOS_EXPECT_NEAR(node.FastGetSolutionStepValue(HYDRAULIC_HEAD), 1.0, 1e-6);
Expand Down Expand Up @@ -186,7 +189,8 @@ KRATOS_TEST_CASE_IN_SUITE(TestExtrapolationProcess_ExtrapolatesCorrectlyForTwoCo
})");

GeoExtrapolateIntegrationPointValuesToNodesProcess process(model, parameters);
process.Execute();
process.ExecuteBeforeSolutionLoop();
process.ExecuteFinalizeSolutionStep();

std::vector<double> expected_values = {1.0, 1.5, 1.5, 1.0, 2.0, 2.0};
std::vector<double> actual_values;
Expand Down Expand Up @@ -223,7 +227,48 @@ KRATOS_TEST_CASE_IN_SUITE(TestExtrapolationProcess_ExtrapolatesCorrectlyForLinea
"list_of_variables" : ["HYDRAULIC_HEAD"]
})");
GeoExtrapolateIntegrationPointValuesToNodesProcess process(model, parameters);
process.Execute();
process.ExecuteBeforeSolutionLoop();
process.ExecuteFinalizeSolutionStep();

std::vector<double> expected_values = {-1, 0, 1, -1, -1, 1};
std::vector<double> actual_values;
actual_values.reserve(model_part.Nodes().size());
std::transform(model_part.Nodes().begin(), model_part.Nodes().end(), std::back_inserter(actual_values),
[](const auto& node) { return node.FastGetSolutionStepValue(HYDRAULIC_HEAD); });

KRATOS_EXPECT_VECTOR_NEAR(actual_values, expected_values, 1e-6)
}

KRATOS_TEST_CASE_IN_SUITE(TestExtrapolationProcess_ExtrapolatesCorrectlyForLinearFields_EvenIfUnrelatedEmptyModelPartsAreSupplied,
KratosGeoMechanicsFastSuiteWithoutKernel)
{
// This test uses the following two-element system.
// 4------3------6
// | El1 | El2 |
// 1------2------5
//
Model model;
auto& model_part = CreateModelPartWithTwoStubElements(model);
model.CreateModelPart("foo"s);
model.CreateModelPart("bar"s);

// Linear field in x between -1 and 1
dynamic_cast<StubElementForNodalExtrapolationTest&>(model_part.Elements()[1]).mIntegrationDoubleValues = {
-inv_sqrt3, inv_sqrt3, inv_sqrt3, -inv_sqrt3};

// Linear field in y between -1 and 1
dynamic_cast<StubElementForNodalExtrapolationTest&>(model_part.Elements()[2]).mIntegrationDoubleValues = {
-inv_sqrt3, -inv_sqrt3, inv_sqrt3, inv_sqrt3};

const auto parameters = Parameters(R"(
{
"model_part_name_list" : ["MainModelPart", "foo", "bar"],
"echo_level" : 0,
"list_of_variables" : ["HYDRAULIC_HEAD"]
})");
GeoExtrapolateIntegrationPointValuesToNodesProcess process(model, parameters);
process.ExecuteBeforeSolutionLoop();
process.ExecuteFinalizeSolutionStep();

std::vector<double> expected_values = {-1, 0, 1, -1, -1, 1};
std::vector<double> actual_values;
Expand Down Expand Up @@ -263,7 +308,8 @@ KRATOS_TEST_CASE_IN_SUITE(TestExtrapolationProcess_ExtrapolatesMatrixCorrectlyFo
})");

GeoExtrapolateIntegrationPointValuesToNodesProcess process(model, parameters);
process.Execute();
process.ExecuteBeforeSolutionLoop();
process.ExecuteFinalizeSolutionStep();

std::vector<Matrix> expected_values = {ScalarMatrix(3, 3, -1), ScalarMatrix(3, 3, 0),
ScalarMatrix(3, 3, 1), ScalarMatrix(3, 3, -1),
Expand Down Expand Up @@ -308,7 +354,8 @@ KRATOS_TEST_CASE_IN_SUITE(TestExtrapolationProcess_ExtrapolatesVectorCorrectlyFo
})");

GeoExtrapolateIntegrationPointValuesToNodesProcess process(model, parameters);
process.Execute();
process.ExecuteBeforeSolutionLoop();
process.ExecuteFinalizeSolutionStep();

std::vector<Vector> expected_values = {ScalarVector(6, -1), ScalarVector(6, 0),
ScalarVector(6, 1), ScalarVector(6, -1),
Expand Down Expand Up @@ -353,7 +400,8 @@ KRATOS_TEST_CASE_IN_SUITE(TestExtrapolationProcess_ExtrapolatesArrayCorrectlyFor
})");

GeoExtrapolateIntegrationPointValuesToNodesProcess process(model, parameters);
process.Execute();
process.ExecuteBeforeSolutionLoop();
process.ExecuteFinalizeSolutionStep();

std::vector<Vector> expected_values = {ScalarVector(3, -1), ScalarVector(3, 0),
ScalarVector(3, 1), ScalarVector(3, -1),
Expand Down
Loading