Geo/14021 interface fluid body flow calculator by WPK4FEM · Pull Request #14215 · KratosMultiphysics/Kratos

WPK4FEM · 2026-02-18T08:36:39Z

📝 Description
Use of permeability and fluid body flow calculator on interface element

🆕 Changelog

Added calculators
Added plane strain and three dimensional unit tests at different orientations.

…ationPoint for Pw_element.

# Conflicts: # applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.h # applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp

…permeability matrix

Added use of fluid body flow calculator and the projected gravity vector needed for it.

I'm missing a - sing somewhere, because now PermeabilityFlow and FluidBodyFlow are equal where they should be opposite.

…d body flow vector. Now they are in an unexpected order for me.

Added unit test for interface at 90 deg with gravity.

additional tests

rfaasse

Thank you for connecting the fluid body flow and permeability calculators to the interface! I did pay slightly less attention to the functions I expect will also be merged back from master (due to the overlap with the coupling) since they've been reviewed before. I didn't see any major issues, just have some minor suggestions/questions.

applications/GeoMechanicsApplication/custom_elements/Pw_element.cpp

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.cpp

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.h

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.cpp

rfaasse · 2026-02-18T12:40:57Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+        {12, array_1d<double, 3>{0.0, 0.0, 0.0}}, {13, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {14, array_1d<double, 3>{0.0, 0.0, 0.0}}, {15, array_1d<double, 3>{0.0, 0.0, 0.0}}};
+
+    // verbeteren, dit zijn indexen en geen knoopnummers


I think this comment is no longer accurate

This was corrected and the comment removed.

rfaasse · 2026-02-18T12:41:28Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+        {12, array_1d<double, 3>{0.0, 0.0, 0.0}}, {13, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {14, array_1d<double, 3>{0.0, 0.0, 0.0}}, {15, array_1d<double, 3>{0.0, 0.0, 0.0}}};
+
+    // verbeteren, dit zijn indexen en geen knoopnummers


Same here, I think the comment can be removed

rfaasse · 2026-02-18T12:42:22Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+    const auto prescribed_displacements = PrescribedDisplacements{
+        {0, array_1d<double, 3>{0.0, 0.0, 0.0}},  {1, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {2, array_1d<double, 3>{0.0, 0.0, 0.0}},  {3, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {4, array_1d<double, 3>{0.0, 0.0, 0.0}},  {5, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {6, array_1d<double, 3>{0.0, 0.0, 0.0}},  {7, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {8, array_1d<double, 3>{0.0, 0.0, 0.0}},  {9, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {10, array_1d<double, 3>{0.0, 0.0, 0.0}}, {11, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {12, array_1d<double, 3>{0.0, 0.0, 0.0}}, {13, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {14, array_1d<double, 3>{0.0, 0.0, 0.0}}, {15, array_1d<double, 3>{0.0, 0.0, 0.0}}};


We could probably leave this part out (I believe displacements are initialized as zeros by default)

It fills the 1st of 3 optional arguments, I wanted to fill it clearly when the following optional arguments are needed.

rfaasse · 2026-02-18T12:44:49Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+
+    for (const auto& [idx, disp] : prescribed_displacements) {


We might be able to make this for-loop a simple (probably templated) utility within this test file, such that we don't need to repeat the for-loop this often 👍

Added a utility function SetVariableOnGeometry in the anonymous namespace and used it throughout this cpp file

rfaasse · 2026-02-18T12:48:22Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+    constexpr auto number_of_pw_dofs       = std::size_t{8};
+    const auto     expected_u_block_vector = Vector{number_of_u_dofs, 0.0};
+    AssertUBlockVectorIsNear(actual_right_hand_side, expected_u_block_vector, number_of_u_dofs, number_of_pw_dofs);
+    // The fluid body flow of a vertical interface subjected to vertical gravity should be zero


Based on the title of this test, I would expect this to be described as a horizontal interface

Suggested change

// The fluid body flow of a vertical interface subjected to vertical gravity should be zero

// The fluid body flow of a horizontal interface subjected to vertical gravity should be zero

Changed, thank you for the sharp observation.

…-fluid-body-flow-calculator Also resolved several merge conflicts.

…and out.

avdg81

Hi Wijtze Pieter,
Thanks a lot for all the hard work that you've done to include the fluid body flow and the permeability contribution. I have several suggestions, but all of them are minor. It's perfectly fine if you'd rather pick up some of the suggestions in a separate PR. I'll leave that up to you.

avdg81 · 2026-02-18T16:06:06Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.cpp

+#include "contribution_calculators/fluid_body_flow_calculator.hpp"
+#include "contribution_calculators/permeability_calculator.hpp"


These two includes have also been added to the header file of the U-Pw interface element. Therefore, it's not needed to repeat them here.

Removed all contribution_calculators includes from the cpp file.

avdg81 · 2026-02-18T16:08:58Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.cpp

    const auto number_of_dofs = GetDofs().size();
-    rLeftHandSideMatrix       = ZeroMatrix{number_of_dofs, number_of_dofs};
+    const auto ignore_undrained =
+        this->GetProperties().Has(IGNORE_UNDRAINED) ? this->GetProperties()[IGNORE_UNDRAINED] : false;


Nitpicking: this-> is redundant here.

Suggested change

this->GetProperties().Has(IGNORE_UNDRAINED) ? this->GetProperties()[IGNORE_UNDRAINED] : false;

GetProperties().Has(IGNORE_UNDRAINED) ? GetProperties()[IGNORE_UNDRAINED] : false;

In addition, we may consider to put this logic in a local helper function, e.g. GetIgnoreUndrained. We could start by adding it to the anonymous namespace. If we later find we also need it elsewhere, we can move it to a utility function class.

Tried a function with the properties as argument.

avdg81 · 2026-02-18T16:18:01Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.cpp


+void UPwInterfaceElement::CalculateAndAssignPermeabilityMatrix(Element::MatrixType& rLeftHandSideMatrix)
+{
+    switch (GetGeometry().PointsNumber()) {


I would suggest to be specific about the geometry that we mean: either the one for the displacement field (GetDisplacementGeometry) or the one for the water pressure field (GetWaterPressureGeometry). GetGeometry refers to the displacement geometry. I'm not sure whether that's what you intended...?

Furthermore, I would suggest to rewrite the switch statement using a map, as was done before for other member function, e.g. CalculateAndAssignUPCouplingMatrix.

Now using the water pressure geometry and did it in the form of a map.

avdg81 · 2026-02-18T16:24:52Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.cpp

+    case 6:
+        CalculateAndAssemblePermeabilityFlowVector<6>(rRightHandSideVector);


This suggests that the interface geometry (either a 3+3 line or a 3+3 surface) is irrelevant here; it's just the number of nodes that matters. Am I correct in saying that?

Its just the number of Pw D.o.F.

avdg81 · 2026-02-18T16:38:55Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.cpp

+    r_water_pressure_mid_geometry.ShapeFunctionsLocalGradients(shape_functions_local_gradient, rIntegrationPoint);
+    // local derivative
+    shape_functions_local_gradient /= r_water_pressure_mid_geometry.DeterminantOfJacobian(rIntegrationPoint);


This code suggests that ShapeFunctionsLocalGradients returns gradients that include the determinants of the Jacobian matrices, which apparently is undesirable here. I was just wondering whether class Geometry perhaps offers API that returns those gradients instead, just so we don't need to correct for that here?

Looked in e.g. line_2d_3.h, but did not find such a gradient.

avdg81 · 2026-02-27T10:12:24Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+        {12, array_1d<double, 3>{0.0, 0.0, 0.0}}, {13, array_1d<double, 3>{0.0, 0.0, 0.0}},
+        {14, array_1d<double, 3>{0.0, 0.0, 0.0}}, {15, array_1d<double, 3>{0.0, 0.0, 0.0}}};
+
+    // verbeteren, dit zijn indexen en geen knoopnummers


Does this (Dutch) comment still apply?

No. Removed it.

avdg81 · 2026-02-27T10:19:40Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+    const auto prescribed_volume_accelerations = std::vector<std::pair<std::size_t, array_1d<double, 3>>>{
+        {0, array_1d<double, 3>{0.0, 0.0, -10.0}},  {1, array_1d<double, 3>{0.0, 0.0, -10.0}},
+        {2, array_1d<double, 3>{0.0, 0.0, -10.0}},  {3, array_1d<double, 3>{0.0, 0.0, -10.0}},
+        {4, array_1d<double, 3>{0.0, 0.0, -10.0}},  {5, array_1d<double, 3>{0.0, 0.0, -10.0}},
+        {6, array_1d<double, 3>{0.0, 0.0, -10.0}},  {7, array_1d<double, 3>{0.0, 0.0, -10.0}},
+        {8, array_1d<double, 3>{0.0, 0.0, -10.0}},  {9, array_1d<double, 3>{0.0, 0.0, -10.0}},
+        {10, array_1d<double, 3>{0.0, 0.0, -10.0}}, {11, array_1d<double, 3>{0.0, 0.0, -10.0}},
+        {12, array_1d<double, 3>{0.0, 0.0, -10.0}}, {13, array_1d<double, 3>{0.0, 0.0, -10.0}},
+        {14, array_1d<double, 3>{0.0, 0.0, -10.0}}, {15, array_1d<double, 3>{0.0, 0.0, -10.0}}};


Instead of defining the 3D acceleration vector many times, we could define is just once (make it a static const variable, e.g. gravity_acceleration_3d). We could make it even nicer by having a nonmember factory function that receives a list of node IDs and a vector to be associated with each node, which then produces these prescribed volume accelerations. Perhaps that is a bit too much for now, so feel free to pick it up in a separate PR.

Did so in many places.

avdg81 · 2026-02-27T10:20:34Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+    // The fluid body flow of a vertical interface subjected to vertical gravity should be zero
+    const auto expected_p_block_vector = Vector{number_of_pw_dofs, 0.0};
+
+    KRATOS_INFO("Actual RHS") << actual_right_hand_side << std::endl;


I think this is a left-over from a debugging session?

Yes, removed it.

avdg81 · 2026-02-27T10:21:55Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+    // The fluid body flow of a vertical interface subjected to vertical gravity should be zero
+    const auto expected_p_block_vector = Vector{number_of_pw_dofs, 0.0};
+
+    KRATOS_INFO("Actual RHS") << actual_right_hand_side << std::endl;


I assume this print statement is no longer needed.

Indeed, removed it.

avdg81 · 2026-02-27T10:23:40Z

...ions/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_U_Pw_interface_element.cpp

+    constexpr auto number_of_pw_dofs       = std::size_t{8};
+    const auto     expected_u_block_vector = Vector{number_of_u_dofs, 0.0};
+    AssertUBlockVectorIsNear(actual_right_hand_side, expected_u_block_vector, number_of_u_dofs, number_of_pw_dofs);
+    auto       flow = (1.0 - 0.5 * sqrt(3.0)) * 62.5;


Nitpicking:

Suggested change

auto flow = (1.0 - 0.5 * sqrt(3.0)) * 62.5;

const auto flow = (1.0 - 0.5 * sqrt3) * 62.5;

made const.

…-fluid-body-flow-calculator' into geo/14021-interface-fluid-body-flow-calculator

…gs const )

…ings const )

…gs const )

…const )

…st )

Compute components once and reuse.

…it anyway. ( Review comment Anne )

-most format statements -use of ublas assignment -left behind print statements.

Gave repeated values a name. Used AssertRhsVectorBlocksAreNear ( once, more to come )

Gave more repeated values a name. Used more AssertRhsVectorBlocksAreNear ( not yet possible on limited precision regression values )

Renamed types Corrected quadrilateral typo's

…Near Gave some values machine precision. Name more repeatedly used values.

WPK4FEM added 16 commits January 29, 2026 13:12

Limited the scope of J_container in CalculateProjectedGravityAtIntegr…

159b2ff

…ationPoint for Pw_element.

1st building version with PermeabilityCalculator

7574fb2

# Conflicts: # applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.h # applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp

building version, not tried yet

da4eb10

Correct orientation of PwB matrix. Use interface version of material …

155d84b

…permeability matrix

Formatting.

0edb224

Indicate per test which calculators to use.

980c33d

More steps towards testing permeability flow.

eae21cd

Changed to assemble for RHS vectors as contributions will overlap.

70b8145

Added use of fluid body flow calculator and the projected gravity vector needed for it.

SonarCloud remark redundant return type.

da6a077

Closer, but no correct yet.

b082d5f

I'm missing a - sing somewhere, because now PermeabilityFlow and FluidBodyFlow are equal where they should be opposite.

Interchanging the sings in the PwBMatrix changes the sign of the flui…

4032b09

…d body flow vector. Now they are in an unexpected order for me.

Preparation for ignore undrained. This still needs unit tests.

0da0559

Avoid ublas assignment

56877da

Project gravity to the conductivity directions.

4403a63

Added unit test for interface at 90 deg with gravity.

correction 3D rotation.

a9345a5

additional tests

formatting

76bf9c7

WPK4FEM requested a review from avdg81 February 18, 2026 08:36

WPK4FEM requested a review from a team as a code owner February 18, 2026 08:36

WPK4FEM linked an issue Feb 18, 2026 that may be closed by this pull request

[GeoMechanicsApplication] Calculator for permeability contributions for the U Pw interface #14022

Open

WPK4FEM self-assigned this Feb 18, 2026

rfaasse reviewed Feb 18, 2026

View reviewed changes

WPK4FEM linked an issue Feb 20, 2026 that may be closed by this pull request

[GeoMechanicsApplication] Calculator for fluid body flow in U Pw interface element #14021

Closed

avdg81 and others added 4 commits February 26, 2026 11:51

Merge remote-tracking branch 'origin/master' into geo/14021-interface…

452f9d4

…-fluid-body-flow-calculator Also resolved several merge conflicts.

Removed copy of function.

0e5b605

Moving similar functionality closer together to let the similarity st…

7304b45

…and out.

Simplified the calculation of the projected gravity a bit

f54a696

avdg81 force-pushed the geo/14021-interface-fluid-body-flow-calculator branch from 452f9d4 to f54a696 Compare February 27, 2026 08:39

avdg81 reviewed Feb 27, 2026

View reviewed changes

WPK4FEM added 2 commits March 2, 2026 11:41

Merge remote-tracking branch 'refs/remotes/origin/geo/14021-interface…

c9a1ad7

…-fluid-body-flow-calculator' into geo/14021-interface-fluid-body-flow-calculator

Reuse nodal value getter with variable name argument

645ab69

WPK4FEM added 17 commits March 2, 2026 14:32

Use map for the stiffness force vector part ( made the necessary thin…

a8c0f23

…gs const )

Use map for the permeability flow vector part ( made the necessary th…

35daff1

…ings const )

Use map for the fluid body flow vector part ( made the necessary thin…

973f815

…gs const )

Use map for the permeability matrix part ( made the necessary things …

e5e5c90

…const )

Use map for the stiffness matrix part ( made the necessary things con…

9e27827

…st )

Corrected comments in test U Pw interface tests

5911826

These includes are in the header file, so not necessary in the .cpp

d8c4df5

Compute components once and reuse.

IGNORE_UNDRAINED through a getter.

3c69d41

Removed geometry argument from CalculatePwBMatrix as we can approach …

f682c03

…it anyway. ( Review comment Anne )

const geometry&

d805995

Removed:

6837aa1

-most format statements -use of ublas assignment -left behind print statements.

Replaced square root of 3

d07d5db

Made flow const

c423f8e

Utility function for setting variables on geometry.

7686f14

Gave repeated values a name. Used AssertRhsVectorBlocksAreNear ( once, more to come )

More use of utility function for setting variables on geometry.

fb2e949

Gave more repeated values a name. Used more AssertRhsVectorBlocksAreNear ( not yet possible on limited precision regression values )

More use of SetVariableOnGeometry

08a674e

Renamed types Corrected quadrilateral typo's

More use of AssertRHSVectorBlocksAreNear and AssertLHSMatrixBlocksAre…

fff28fb

…Near Gave some values machine precision. Name more repeatedly used values.

WPK4FEM requested review from avdg81 and rfaasse March 3, 2026 12:29

avdg81 added the GeoMechanics Issues related to the GeoMechanicsApplication label Mar 3, 2026

avdg81 added this to GeoMechanicsApplication (Deltares) Mar 3, 2026

github-project-automation bot moved this to In progress in GeoMechanicsApplication (Deltares) Mar 3, 2026

	// The fluid body flow of a vertical interface subjected to vertical gravity should be zero
	// The fluid body flow of a horizontal interface subjected to vertical gravity should be zero

		#include "contribution_calculators/fluid_body_flow_calculator.hpp"
		#include "contribution_calculators/permeability_calculator.hpp"

	this->GetProperties().Has(IGNORE_UNDRAINED) ? this->GetProperties()[IGNORE_UNDRAINED] : false;
	GetProperties().Has(IGNORE_UNDRAINED) ? GetProperties()[IGNORE_UNDRAINED] : false;

		case 6:
		CalculateAndAssemblePermeabilityFlowVector<6>(rRightHandSideVector);

	auto flow = (1.0 - 0.5 * sqrt(3.0)) * 62.5;
	const auto flow = (1.0 - 0.5 * sqrt3) * 62.5;

Conversation

WPK4FEM commented Feb 18, 2026

Uh oh!

rfaasse left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

avdg81 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development