[GeoMechanicsApplication] Add a minimal 3+3 line interface element #12674

rfaasse · 2024-09-16T09:38:55Z

📝 Description
This PR adds a minimal implementation for a line interface element (which could be 2+2 or 3+3 depending on the geometry passed on during construction). The functionality (e.g. for the LHS/RHS calculations and the other functions) is unit tested.

…ero matrix with the correct size

…geometry is the same as the underlying mid geometry and changed the tests accordingly

…-element

… for a line interface element

…x for a line interface element

Moved the application-wide tolerance values for testing to a dedicated namespace.

Don't increment a local variable at function call time. Increment first, then pass the updated value to the function.

To better reflect what the test case actually covers.

…or a line interface element

So far we expect a zero vector (this will soon change). Also renamed a local variable.

…imal-3_3-line-interface-element

So far, the function under test doesn't calculate anything. It simply returns a default-constructed vector.

The calculation follows the formulation found in any textbook about the finite element method

…-element

…es counterclockwise

…AIN` and `CAUCHY_STRESS_VECTOR` might be overwritten by the base class

…-interface-element' into geo/12649-add-a-minimal-3_3-line-interface-element

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

avdg81 · 2024-09-16T11:56:03Z

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

+    KRATOS_EXPECT_EQ(degrees_of_freedom.size(), 8);
+    const std::vector<double> expected_dof_values = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
+    for (std::size_t i = 0; i < degrees_of_freedom.size(); i++) {
+        KRATOS_EXPECT_DOUBLE_EQ(degrees_of_freedom[i]->GetSolutionStepValue(), expected_dof_values[i]);
+    }


As discussed before, perhaps we'd better check the degrees of freedom themselves, rather than their associated values (since that should be covered by the unit tests for the degrees of freedom class itself).

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

avdg81 · 2024-09-16T12:05:35Z

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

+
+    // Assert
+    Vector expected_cauchy_stress{2};
+    expected_cauchy_stress <<= 10.0, 2.0;


This is the traction vector using the local axis system. Is that we intend?

This is intended, but let's discuss with @WPK4FEM before merging

For me the tractions ( I also think its clearer to call the thing traction in the test, even if we are forced to use functions for strain and Cauchy stress ) are clearer in the local axes than in a global system.

For our continuum elements the local and the global axes system coincide, for interface elements that would be very uncomfortable. Please keep the interface element integration point results local.

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

…ts instead of pointers

mnabideltares

@rfaasse and @avdg81
Thank you for adding interface element. Here some of my observations:

mnabideltares · 2024-09-17T08:35:36Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.h

+
+    using Element::GeometryType;
+    using Element::PropertiesType;
+    using BaseType = Element;


Why to not remove "using BaseType = Element" and use "Element" everywhere in the function, to be more clear.

Initially we did it to be able to insert a class into the hierarchy more easily. In that case some of these Element:: should remain as they are, and some should be changed to a new base class. However, we had some discussion about it before already if we really need it and after your comment decided to indeed remove BaseType, good suggestion!

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

mnabideltares · 2024-09-17T08:42:27Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.h

+    using Element::PropertiesType;
+    using BaseType = Element;
+
+    LineInterfaceElement() = default;


As this is a peurly header file, preferly "default" to be in the cpp, to be consistent with the other constructors. Also to be more clear in reading the code.

The argument here is that (since it's default only, this doesn't hold when there is a real implementation), in this way we don't have to go to the source file just to see it's a default constructor, but I think it's a matter of preference.

applications/GeoMechanicsApplication/custom_elements/line_interface_element.h

mnabideltares · 2024-09-17T09:53:49Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+                                           const Geometry<GeometricalObject::NodeType>::Pointer& rGeometry,
+                                           const Properties::Pointer& rProperties)
+    : Element(NewId, rGeometry, rProperties),
+      mIntegrationScheme(std::make_unique<LobattoIntegrationScheme>(GetGeometry().PointsNumber() / 2)),


Here, the integration scheme is defined based on the number of the points. For now, it is enough, but for further development, it is not flexible. In other elements, we define "GetIntegrationMethod" to be able define the scheme based on the desired order. For example, consider different order element. However, this further investigation/definition can be in a different PR.

The problem with GetIntegrationMethod is that the enum it returns doesn't have entries for Lobatto integration, which is why we cannot directly use it.

mnabideltares · 2024-09-17T09:59:53Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+void LineInterfaceElement::CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo&)
+{
+    // Currently, the left-hand side matrix only includes the stiffness matrix. In the future, it
+    // will also include water pressure contributions and coupling terms.


It is for my curiousity. Do you then set the pressure as zero? As far as I know, there is no solver for U available in Geo, and our solvers are Pw and U-Pw (and T)

We'll have to double check if this method works when looking at integration tests (first indication is it does, but we'll have to see if that holds), but we're just not giving Pw as degrees of freedom in these interfaces, meaning the interfaces shouldn't change water pressure at all. It basically works as a 'diff order' element, but the pressure geometry is set to empty.

mnabideltares · 2024-09-17T10:43:54Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+{
+    // Currently, the right-hand side only includes the internal force vector. In the future, it
+    // will also include water pressure contributions and coupling terms.
+    const auto local_b_matrices = CalculateLocalBMatricesAtIntegrationPoints();


Calculation of B-Matrix is done twice, once for left hand, then for right hand, by means of "CalculateLocalBMatricesAtIntegrationPoints()". It is more efficient to calculate it once for the left hand side, then use it for the right hand side.

Since we're using the CalculateRightHandSide function from the Element interface, we cannot change the signature of the function. I agree we do double work now, so we could look into caching certain lists if we see it's a performance bottleneck. For now, I think it would add complexity to the element that might have limited value until we solve more prominent performance problems.

mnabideltares · 2024-09-17T10:46:18Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+    const auto local_b_matrices = CalculateLocalBMatricesAtIntegrationPoints();
+    const auto relative_displacements = CalculateRelativeDisplacementsAtIntegrationPoints(local_b_matrices);
+    const auto tractions = CalculateTractionsAtIntegrationPoints(relative_displacements);
+    const auto integration_coefficients = CalculateIntegrationCoefficients();


"CalculateIntegrationCoefficients();" is done for left as well as, right hand side. Keep the values in an array reuse them

Valid point, but as per my previous comment, I'd rather choose readability for now and optimize performance later.

mnabideltares · 2024-09-17T10:59:14Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+Element::Pointer LineInterfaceElement::Create(IndexType               NewId,
+                                              const NodesArrayType&   rNodes,
+                                              PropertiesType::Pointer pProperties) const
+{
+    return Create(NewId, this->GetGeometry().Create(rNodes), pProperties);
+}
+
+Element::Pointer LineInterfaceElement::Create(IndexType               NewId,
+                                              GeometryType::Pointer   pGeometry,
+                                              PropertiesType::Pointer pProperties) const
+{
+    return make_intrusive<LineInterfaceElement>(NewId, pGeometry, pProperties);
+}


Preferly Creat to be after the constructors to be consistent with the other element classes design

Agreed, moved

WPK4FEM

Dear Anne and Richard,
Thank you for the element setup. Not much to complain about from my side.
Wijtze Pieter

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

WPK4FEM · 2024-09-18T07:54:19Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+        // For interface elements, the shape function gradients are not defined, since these are
+        // non-continuum elements. Therefore, we pass an empty matrix.


Shape function gradients are always defined, they are not used for interface elements.

Adjusted the comment, thanks for the clarification!

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

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

rfaasse

Thank you for the thorough reviews! I think I processed/answered to all of them, please let me know if you think any other changes are necessary

applications/GeoMechanicsApplication/custom_elements/line_interface_element.h

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

rfaasse · 2024-09-18T13:43:29Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+                                           const Geometry<GeometricalObject::NodeType>::Pointer& rGeometry,
+                                           const Properties::Pointer& rProperties)
+    : Element(NewId, rGeometry, rProperties),
+      mIntegrationScheme(std::make_unique<LobattoIntegrationScheme>(GetGeometry().PointsNumber() / 2)),


The problem with GetIntegrationMethod is that the enum it returns doesn't have entries for Lobatto integration, which is why we cannot directly use it.

rfaasse · 2024-09-18T13:45:50Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+void LineInterfaceElement::CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo&)
+{
+    // Currently, the left-hand side matrix only includes the stiffness matrix. In the future, it
+    // will also include water pressure contributions and coupling terms.


We'll have to double check if this method works when looking at integration tests (first indication is it does, but we'll have to see if that holds), but we're just not giving Pw as degrees of freedom in these interfaces, meaning the interfaces shouldn't change water pressure at all. It basically works as a 'diff order' element, but the pressure geometry is set to empty.

rfaasse · 2024-09-18T13:49:50Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

+{
+    // Currently, the right-hand side only includes the internal force vector. In the future, it
+    // will also include water pressure contributions and coupling terms.
+    const auto local_b_matrices = CalculateLocalBMatricesAtIntegrationPoints();


Since we're using the CalculateRightHandSide function from the Element interface, we cannot change the signature of the function. I agree we do double work now, so we could look into caching certain lists if we see it's a performance bottleneck. For now, I think it would add complexity to the element that might have limited value until we solve more prominent performance problems.

rfaasse · 2024-09-19T07:15:01Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.h

+
+    using Element::GeometryType;
+    using Element::PropertiesType;
+    using BaseType = Element;


Initially we did it to be able to insert a class into the hierarchy more easily. In that case some of these Element:: should remain as they are, and some should be changed to a new base class. However, we had some discussion about it before already if we really need it and after your comment decided to indeed remove BaseType, good suggestion!

rfaasse · 2024-09-19T07:16:20Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.h

+    using Element::PropertiesType;
+    using BaseType = Element;
+
+    LineInterfaceElement() = default;


The argument here is that (since it's default only, this doesn't hold when there is a real implementation), in this way we don't have to go to the source file just to see it's a default constructor, but I think it's a matter of preference.

applications/GeoMechanicsApplication/custom_elements/line_interface_element.h

rfaasse · 2024-09-19T07:19:27Z

applications/GeoMechanicsApplication/custom_elements/line_interface_element.h

+    void Initialize(const ProcessInfo& rCurrentProcessInfo) override;
+
+    Element::Pointer Create(IndexType NewId, const NodesArrayType& rNodes, PropertiesType::Pointer pProperties) const override;
+    Element::Pointer Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override;


Good catch! Added a check that calls the base class check, verifies that the constitutive laws are there and calls check on them.

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

WPK4FEM

Dear Anne and Richard,
To me this looks quite good to go, I don't find things that are really worth putting substantial effort in first anymore.
Thank you, Wijtze Pieter

WPK4FEM · 2024-09-19T09:18:01Z

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

+    expected_right_hand_side <<= 0.33333333, 1.66666667, 0.33333333, 1.66666667, 1.33333333,
+        6.66666667, -0.33333333, -1.66666667, -0.33333333, -1.66666667, -1.33333333, -6.66666667;


Should we use this or write 1/3, 5/3, 1/3, 5/3, 4/3, 19/3, -1/3 et.c

Similarly maybe at line 429, 463 and the similar things where square roots of 3 are involved

Guess this is a matter of preference.

For this set of 'n/3' I completely agree, if I remember correctly the 429 and 463 examples, the analytical solution was not a lot more clear than just the numbers generated by Annes python 'check' script.

I'll wait for the sonarcloud results here: if I need to make changes because of that, I'll change this line, otherwise I'd propose to merge this and avoid another pipeline wait.

mnabideltares

It looks OK to me.

rfaasse and others added 30 commits August 27, 2024 14:38

First commit, added line interface element deriving from Element

ccf38a0

Added create function + unit test

f4844ae

Added create function with node input + unit test

6d764c6

Implemented the function to return the dofs

a669049

Added function for the equation ids

bfddbc5

Minor refactoring: added private GetDofs function

1fdca88

Removed some redundant sets

d2a0f2e

Added implementation for CalculateLeftHandSide which just returns a z…

cbcb60b

…ero matrix with the correct size

Fixed signed/unsigned issue breaking the GCC pipeline

5002b52

Made a fix to make sure the working space dimension of the interface …

37c7718

…geometry is the same as the underlying mid geometry and changed the tests accordingly

Removed unused variable in test to fix pipeline build

e9a2113

Merge branch 'master' into geo/12649-add-a-minimal-3_3-line-interface…

8107cd5

…-element

Started to extend unit test for calculating the left hand side matrix…

3dd885d

… for a line interface element

Finished to extend unit test for calculating the left hand side matri…

1586535

…x for a line interface element

Removed a redundant assignment

e8955ca

Replaced a raw for loop by a transform operation

801583c

Fixed a compile-time error when using Unity builds

68cd292

Moved the application-wide tolerance values for testing to a dedicated namespace.

Fixed two code smells found by SonarQube

2767aee

Don't increment a local variable at function call time. Increment first, then pass the updated value to the function.

Renamed a unit test case

ffc9e65

To better reflect what the test case actually covers.

Added a first unit test for checking the right hand side calculated f…

d7b7b6f

…or a line interface element

Test for the expected RHS vector

37ae3ab

So far we expect a zero vector (this will soon change). Also renamed a local variable.

Merge remote-tracking branch 'origin/master' into geo/12649-add-a-min…

f1c8ced

…imal-3_3-line-interface-element

Added a skeleton test case for calculating the internal force vector

0c745dc

So far, the function under test doesn't calculate anything. It simply returns a default-constructed vector.

Properly implemented function CalculateInternalForceVector

f041353

The calculation follows the formulation found in any textbook about the finite element method

Added rotation matrix with a test for the left-hand side

bf36825

Fixed memory issue (don't use auto for matrices)

1e00062

Merge branch 'master' into geo/12649-add-a-minimal-3_3-line-interface…

83490e3

…-element

Modified a test case to use a line interface element rotated 30 degre…

43454ba

…es counterclockwise

The constitutive laws at the integration points can now be retrieved

588f071

Check the cloned constitutive laws that are kept by the element

81981c7

rfaasse and others added 5 commits September 16, 2024 10:54

Added a test for 3 plus 3 noded interface element

298d241

Fixed a potential issue, where the integration point results for `STR…

806666e

…AIN` and `CAUCHY_STRESS_VECTOR` might be overwritten by the base class

Merge remote-tracking branch 'origin/geo/12649-add-a-minimal-3_3-line…

aab721e

…-interface-element' into geo/12649-add-a-minimal-3_3-line-interface-element

Merge remote-tracking branch 'origin/geo/12649-add-a-minimal-3_3-line…

884059b

…-interface-element' into geo/12649-add-a-minimal-3_3-line-interface-element

Removed clang-tidy warnings

a6b5a03

rfaasse linked an issue Sep 16, 2024 that may be closed by this pull request

[GeoMechanicsApplication] Add a minimal 3+3 line interface element #12649

Closed

rfaasse assigned rfaasse and avdg81 Sep 16, 2024

avdg81 reviewed Sep 16, 2024

View reviewed changes

rfaasse added 5 commits September 16, 2024 16:10

Clean up

6723b5f

Changed helper functions to return the concrete line interface elemen…

7ebfc51

…ts instead of pointers

Renamed some variables and fixed const-correctness of tests

b1f9147

Fix 3+3 noded test and changed scope of doflist test slightly

565d397

Reverted unintended change to test

4c705c1

rfaasse added the GeoMechanics Issues related to the GeoMechanicsApplication label Sep 16, 2024

rfaasse requested review from mnabideltares and markelov208 September 16, 2024 15:20

rfaasse marked this pull request as ready for review September 16, 2024 15:20

Fixed two clang-tidy/sonar cloud issues

3e58bb5

mnabideltares reviewed Sep 17, 2024

View reviewed changes

WPK4FEM reviewed Sep 18, 2024

View reviewed changes

Processed review comments

34ff383

rfaasse commented Sep 19, 2024

View reviewed changes

rfaasse requested review from mnabideltares and WPK4FEM September 19, 2024 07:26

Follow rule of 5 to fix noexcept code smell

1575958

WPK4FEM approved these changes Sep 19, 2024

View reviewed changes

mnabideltares approved these changes Sep 19, 2024

View reviewed changes

rfaasse merged commit 1677ea8 into master Sep 19, 2024
9 of 11 checks passed

rfaasse deleted the geo/12649-add-a-minimal-3_3-line-interface-element branch September 19, 2024 13:43

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[GeoMechanicsApplication] Add a minimal 3+3 line interface element #12674

[GeoMechanicsApplication] Add a minimal 3+3 line interface element #12674

rfaasse commented Sep 16, 2024 •

edited

Loading

avdg81 Sep 16, 2024

rfaasse Sep 19, 2024

avdg81 Sep 16, 2024

rfaasse Sep 16, 2024

WPK4FEM Sep 17, 2024

mnabideltares left a comment

mnabideltares Sep 17, 2024

rfaasse Sep 19, 2024

mnabideltares Sep 17, 2024

rfaasse Sep 19, 2024

mnabideltares Sep 17, 2024

rfaasse Sep 18, 2024

mnabideltares Sep 17, 2024

rfaasse Sep 18, 2024

mnabideltares Sep 17, 2024

rfaasse Sep 18, 2024

mnabideltares Sep 17, 2024

rfaasse Sep 18, 2024

mnabideltares Sep 17, 2024

rfaasse Sep 18, 2024

WPK4FEM left a comment

WPK4FEM Sep 18, 2024

rfaasse Sep 18, 2024

rfaasse left a comment

rfaasse Sep 18, 2024

rfaasse Sep 18, 2024

rfaasse Sep 18, 2024

rfaasse Sep 19, 2024

rfaasse Sep 19, 2024

rfaasse Sep 19, 2024

WPK4FEM left a comment

WPK4FEM Sep 19, 2024

rfaasse Sep 19, 2024

mnabideltares left a comment

		// For interface elements, the shape function gradients are not defined, since these are
		// non-continuum elements. Therefore, we pass an empty matrix.

		expected_right_hand_side <<= 0.33333333, 1.66666667, 0.33333333, 1.66666667, 1.33333333,
		6.66666667, -0.33333333, -1.66666667, -0.33333333, -1.66666667, -1.33333333, -6.66666667;

[GeoMechanicsApplication] Add a minimal 3+3 line interface element #12674

[GeoMechanicsApplication] Add a minimal 3+3 line interface element #12674

Conversation

rfaasse commented Sep 16, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mnabideltares left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

WPK4FEM left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rfaasse left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

WPK4FEM left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mnabideltares left a comment

Choose a reason for hiding this comment

rfaasse commented Sep 16, 2024 •

edited

Loading