Skip to content

[GeoMechanicsApplication] Extend class LineInterfaceElement such that it also includes planar interface elements#13634

Merged
markelov208 merged 28 commits intomasterfrom
geo/13334-extend-LineInterfaceElement
Jul 17, 2025
Merged

[GeoMechanicsApplication] Extend class LineInterfaceElement such that it also includes planar interface elements#13634
markelov208 merged 28 commits intomasterfrom
geo/13334-extend-LineInterfaceElement

Conversation

@markelov208
Copy link
Contributor

@markelov208 markelov208 commented Jul 15, 2025

📝 Description
A brief description of the PR.

Please mark the PR with appropriate tags:

  • renamed LineInterfaceElement as InterfaceElement
  • added pStressStatePolicy in constructors
  • added LumpedIntegrationScheme
  • added LocalSpaceDimension check
  • added Calculate3DRotationMatrixForPlaneGeometry in geometry_utilities
  • added unit tests
  • updated element_hierarchy.svg

Copy link
Contributor

@rfaasse rfaasse left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for adding this well-tested extension of the interface element. I have some suggestions, but nothing major

Copy link
Contributor

@avdg81 avdg81 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Gennady,

Thanks a lot for generalizing the line interface element to an interface element that supports both line and surface geometry. I'm also happy to see that you've added lots of unit tests. Well done! I have several suggestions, but they're all minor. We may want to revisit the test code at some point, but for now this looks good enough to me. Let's try to get this in as soon as possible.

: Element(NewId, rGeometry, rProperties),
mIntegrationScheme(std::make_unique<LobattoIntegrationScheme>(GetGeometry().PointsNumber() / 2)),
mStressStatePolicy(std::make_unique<Line2DInterfaceStressState>())
InterfaceElement::InterfaceElement(IndexType NewId,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm just wondering whether we should name this class GeoInterfaceElement to avoid future name clashes with other applications? Alternatively, we could consider to wrap this class in a Geo namespace. I'm curious how you and @rfaasse feel about this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed. Currently, it is not changed.

Comment on lines +73 to +76
InterfaceElement::InterfaceElement(IndexType NewId,
const GeometryType::Pointer& rGeometry,
std::unique_ptr<StressStatePolicy> pStressStatePolicy)
: Element(NewId, rGeometry), mpStressStatePolicy(std::move(pStressStatePolicy))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just something to consider: perhaps we can use forwarding constructors to avoid some code duplication?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry I've heard about the forwarding constructors but AI does not help. Perhaps, we can discuss it tomorrow.

Comment on lines +216 to +219
template <typename TElementFactory>
InterfaceElement CreateAndInitializeElement(TElementFactory&& rFactory,
const Properties::Pointer& rProperties,
const std::vector<std::pair<std::size_t, array_1d<double, 3>>>& rDisplacements = {})
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I find the first function parameter hard to read (rvalue reference). When I look at its usages, my understanding is that you'd like to pass a function object that builds an interface element given the model and the desired properties. Perhaps removing the && already helps?

I'm also a bit worried about the model instance here, since model will go out of scope as soon as this function returns. Any other objects that keep a pointer or reference to that model, will then have dangling pointers and references, which is very dangerous and hard to debug.
I've discussed this with @rfaasse and we agreed that we can look into that later, in order to not delay merging this PR.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now, removed &&. It still works ;)

Comment on lines +648 to +650
const auto p_geometry = std::make_shared<TriangleInterfaceGeometry3D3Plus3Noded>(nodes);
const InterfaceElement element(0, p_geometry, p_properties,
std::make_unique<Line2DInterfaceStressState>());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also here I'm confused. We're creating a surface interface element (that's what the geometry suggests), but then we pass it a stress state instance that is intended for 2D line interfaces. Perhaps we should not allow that? At least, we shouldn't showcase that in a unit test, in my opinion.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, another mistake. I think we should not allow. Perhaps, a space dimension shall be added to the StressStatePolicy then we can check a correspondence of the geometry and the stress state.

@markelov208 markelov208 requested a review from avdg81 July 16, 2025 19:33
@markelov208 markelov208 requested a review from rfaasse July 16, 2025 19:33
Copy link
Contributor Author

@markelov208 markelov208 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Richard and Anne, thank you very much for the reviews. Not everything is changed, hope we can discuss it.

Copy link
Contributor

@avdg81 avdg81 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Gennady, thanks for processing the review suggestions so swiftly. I feel we're almost there. I have one more question about the function that calculates the expected LHS matrices.

Comment on lines +232 to +261
Matrix CreateExpectedLeftHandSideForTriangleElement(std::size_t dimension,
double NormalStiffness,
double ShearStiffness,
double divider1 = 3.0,
double divider2 = 3.0)
{
Matrix expected_left_hand_side = ZeroMatrix(dimension, dimension);
double value;
std::size_t counter = 0;
std::size_t half_dimension = dimension / 2;
for (std::size_t i = 0; i < half_dimension; i++) {
if (++counter == 3) {
value = NormalStiffness;
counter = 0;
} else {
value = ShearStiffness;
}
if (i < 9) {
value /= divider1;
} else {
value /= divider2;
}
expected_left_hand_side(i, i) = value;
expected_left_hand_side(i, i + half_dimension) = -value;
expected_left_hand_side(i + half_dimension, i) = -value;
expected_left_hand_side(i + half_dimension, i + half_dimension) = value;
}

return expected_left_hand_side;
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I understand that you'd like to reuse one or a few of the expected left hand side matrices, but I'm not sure whether this is the way to go. I don't understand the implementation of this function. How can I check that this makes sense from a physics/FEM point of view? My understanding of @rfaasse's suggestion was to extract a function that returns a matrix with hard-coded numbers, since some of these matrices were repeated.
If we'd like to keep the current implementation, I feel that we need a comment that explains the calculation, perhaps with some references to finite element books/papers?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for the clarification. Now there is a function that provide these values. I put a comment that the values are taken from the element. The theory, it looks 3D case with matrix of 18x18 is consistent with 2D case 8x8 but 3D case of 36x36 is something special.

Comment on lines +430 to +433
const auto prescribed_displacements = PrescribedDisplacements{
{2, array_1d<double, 3>{0.2, 0.5, 0.0}}, {3, array_1d<double, 3>{0.2, 0.5, 0.0}}};
auto element = CreateAndInitializeElement(CreateHorizontalUnitLength2Plus2NodedLineInterfaceElementWithUDofs,
p_properties, prescribed_displacements);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Much clearer. Thank you.

Copy link
Contributor

@rfaasse rfaasse left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for processing the review comments, to me this looks ready to go!

@markelov208 markelov208 merged commit ec79f44 into master Jul 17, 2025
10 checks passed
@markelov208 markelov208 deleted the geo/13334-extend-LineInterfaceElement branch July 17, 2025 12:45
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

[GeoMechanicsApplication] Extend class LineInterfaceElement such that it also includes planar interface elements

3 participants