[beaminteraction]Dual Hermite shapefunctions for Lagrange Multipliers

[dharinib98]

Description and Context

This is an initial working version of the implementation of dual Hermite shape functions in beaminteraction. The current implementation still requires substantial refactoring to improve code efficiency and maintainability.

Interested Parties

@mayrmt @isteinbrecher @maxfirmbach

Related Issues and Pull Requests

[dharinib98]

@isteinbrecher I moved the dual Hermite implementation to a separate header and addressed most of the other comments. Two points are still open: I still instantiate the rotational variants as well, since the same templated factory path is used for both positional and rotational mortar pairs. Otherwise, configurations with t_hermite_dual can still be generated and lead to linker errors.
The shape-function equations are also still defined through the FEM utility function. Moving them fully to the header complained about some array-bound issues, I am not sure although exactly why. I think both these can be cleaned up once the rotational dual Hermite path is handled properly. Please let me know what you think about the current state of this PR.

[isteinbrecher]

Thanks fort the changes @dharinib98 this is already looking a lot better.

Please have a look at my comments and letzt discuss if there are any questions.

Also, can you add an input file test with the dual shape functions?

[dharinib98]

@isteinbrecher @maxfirmbach @mayrmt This is ready for review again! The core unit test fails in the pipeline here, because of the unit test I added for the dual hermite shape functions. However locally these tests pass for me and I not quite sure why they are problematic here.

[maxfirmbach]

@dharinib98 Looking forward to use this feature too 😉. Looks alright to me ... I just have some minor comments and questions.

[maxfirmbach]

@isteinbrecher Is this what we want here? Of course using a struct forces us to a specific naming scheme, which looks a bit off in the template initialization later on. I don't have a strong opinion on this, just want to know if this would make sense for the other template parameters too in the future.

[isteinbrecher]

I agree @maxfirmbach. I am not sure what the best way is, we could use an alias to define a name like t_hermite_dual after creating the struct.

[maxfirmbach]

Failing tests produce this output:

The difference between Phi(i) and Phi_ref[i] is 0.043972369034521064, which exceeds 1e-10, where
Phi(i) evaluates to -0.052766842841425188,
Phi_ref[i] evaluates to -0.0087944738069041239, and
1e-10 evaluates to 1e-10.

/__w/4C/4C/src/core/fem/tests/utils/4C_fem_general_utils_fem_shapefunctions_test.cpp:56: Failure
The difference between Phi(i) and Phi_ref[i] is 16.761378432263328, which exceeds 1e-10, where
Phi(i) evaluates to -14.366895799082851,
Phi_ref[i] evaluates to 2.3944826331804752, and
1e-10 evaluates to 1e-10.

[  FAILED  ] ElementShapeFunctionsTest.TestDualHermiteLine2 (0 ms)
[ RUN      ] ElementShapeFunctionsTest.TestDualHermiteLine2Length03
/__w/4C/4C/src/core/fem/tests/utils/4C_fem_general_utils_fem_shapefunctions_test.cpp:75: Failure
The difference between Phi(i) and Phi_ref[i] is 0.14657456344840356, which exceeds 1e-10, where
Phi(i) evaluates to -0.17588947613808398,
Phi_ref[i] evaluates to -0.029314912689680415, and
1e-10 evaluates to 1e-10.

/__w/4C/4C/src/core/fem/tests/utils/4C_fem_general_utils_fem_shapefunctions_test.cpp:75: Failure
The difference between Phi(i) and Phi_ref[i] is 55.871261440877753, which exceeds 1e-10, where
Phi(i) evaluates to -47.889652663609503,
Phi_ref[i] evaluates to 7.9816087772682511, and
1e-10 evaluates to 1e-10.

[isteinbrecher]

I agree @maxfirmbach. I am not sure what the best way is, we could use an alias to define a name like t_hermite_dual after creating the struct.

[isteinbrecher]

Thanks for the changes @dharinib98!

@isteinbrecher Will the unit tests actually be enough to check if the ref_length_ is set correctly for now since we have two lengths being checked for in there. Also do we want to add an input file test and a vtu check already in this PR or make it a follow up PR?

We still need the vtu test, because we previously saw that the unit test alone is not enough. Did you verify that the vtu tests produce the correct result?

One more point that we can maybe discuss in person, if we should add a runtime test that the D matrix is actually diagonal when using the dual shape functions.

[mayrmt]

@isteinbrecher I have discussed the idea of testing for a diagonal D with @dharinib98 some time ago, though it is not clear to me, yet, how to achieve this in 4C. It probably must be through solve form of unit test...

Most importantly, I agree that this property should be tested.

[dharinib98]

@isteinbrecher I have discussed the idea of testing for a diagonal D with @dharinib98 some time ago, though it is not clear to me, yet, how to achieve this in 4C. It probably must be through solve form of unit test...

Most importantly, I agree that this property should be tested.

yeah! I am not 100 percent sure either, how this can be done or rather how this can be done correctly. We can discuss about this during the week.

[mayrmt]

@isteinbrecher @dharinib98 A basic question: if we use dual shape functions and we have ensured through unit testing, that the implementation actually satisfies the biorthogonality conditions, why can't we exploit that when assembling D? We know, that there is only one value, so it should be sufficient to compute just this value and put it into the correct location in the matrix. Then, it will be "diagonal" by construction and we do not have to test, if it is diagonal. Furthermore, this avoids to compute many matrix entries, that we expect to evaluate to zero anyways.

What are your thoughts on this?

[dharinib98]

@isteinbrecher @dharinib98 A basic question: if we use dual shape functions and we have ensured through unit testing, that the implementation actually satisfies the biorthogonality conditions, why can't we exploit that when assembling D? We know, that there is only one value, so it should be sufficient to compute just this value and put it into the correct location in the matrix. Then, it will be "diagonal" by construction and we do not have to test, if it is diagonal. Furthermore, this avoids to compute many matrix entries, that we expect to evaluate to zero anyways.

What are your thoughts on this?

@mayrmt The main issue with this would be (since we only deal with positional coupling and not full coupling as yet) to put the computed values into the correct location, since accessing the positional dofs and the lambda translation dofs are not easy with the current tools available. Another personal concern, might be that we would be needing a lot of extra lines of code just for this special case and might make the the current assembly structure that is in place in beaminteraction look a bit messy.

[isteinbrecher]

Thanks for addressing all the points @dharinib98 . From my side only two things are missing, the move of the diagonal check to the beaminteratcion utils and the unittest.

[isteinbrecher]

@isteinbrecher @dharinib98 A basic question: if we use dual shape functions and we have ensured through unit testing, that the implementation actually satisfies the biorthogonality conditions, why can't we exploit that when assembling D? We know, that there is only one value, so it should be sufficient to compute just this value and put it into the correct location in the matrix. Then, it will be "diagonal" by construction and we do not have to test, if it is diagonal. Furthermore, this avoids to compute many matrix entries, that we expect to evaluate to zero anyways.

In theory we could do that. However, with the current code structure, we are able to combine all the different cases into one general implementation. I suspect the overhead of computing the zero entries negligible. Also, if at some point we will want to consider beams that are not fully embedded and then we will need integration again, also for dual shape functions.