Implement Hessians and Identity Mapping for embedded elements

[henrij22]

This does not correctly work at the moment

This implements two things:

1. Compute Hessians for embedded elements
2. Apply geometry mapping to second derivatives.

Linking also for reference : Ferrite-FEM/Tensors.jl#188

[codecov]

Codecov Report

❌ Patch coverage is 0% with 70 lines in your changes missing coverage. Please review.
✅ Project coverage is 0.00%. Comparing base (f99e206) to head (d3214ec).

Files with missing lines	Patch %	Lines
src/FEValues/FunctionValues.jl	0.00%	59 Missing ⚠️
src/FEValues/GeometryMapping.jl	0.00%	10 Missing ⚠️
src/PointEvalHandler.jl	0.00%	1 Missing ⚠️

❗ There is a different number of reports uploaded between BASE (f99e206) and HEAD (d3214ec). Click for more details.

HEAD has 4 uploads less than BASE

Flag	BASE (f99e206)	HEAD (d3214ec)
	5	1

Additional details and impacted files

@@            Coverage Diff             @@
##           master   #1257       +/-   ##
==========================================
- Coverage   94.19%   0.00%   -94.20%     
==========================================
  Files          40      40               
  Lines        6662    6628       -34     
==========================================
- Hits         6275       0     -6275     
- Misses        387    6628     +6241     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[lijas]

Nice!

Can you also add the hessian equation (for embeded elements) to the docs, similar to how we do it here: https://ferrite-fem.github.io/Ferrite.jl/stable/topics/FEValues/#Identity-mapping

And add a reference to where it is taken from (or perhaps show the derivation :) )

[lijas]

where did the computation of the hessian go?

[henrij22]

Aha, I only added it in above. This code then might not be tested atm. Also, what is the difference between calculate_mapping(gip::ScalarInterpolation) and calculate_mapping(geo_mapping::GeometryMapping). They seem quite similar.

[lijas]

I think we can add tests for the function_hessian also, as we do here:

                coords_nl = [x + rand(x) * 0.01 for x in coords] # add some displacement to nodes
                reinit!(cv, coords_nl)

                _ue_nl = [u_funk(coords_nl[i], V, G, H) for i in 1:n_basefunc_base]
                ue_nl = reinterpret(Float64, _ue_nl)

                for i in 1:getnquadpoints(cv)
                    xqp = spatial_coordinate(cv, i, coords_nl)
                    Hqp, Gqp, Vqp = Tensors.hessian(x -> function_value_from_physical_coord(func_interpol, coords_nl, x, ue_nl), xqp, :all)
                    @test function_value(cv, i, ue_nl) ≈ Vqp
                    @test function_gradient(cv, i, ue_nl) ≈ Gqp
                    if update_hessians
                        @test Ferrite.function_hessian(cv, i, ue_nl) ≈ Hqp
                    end
                end

Can you add a similar test in this @testset for hessians? (I am not sure if function_value_from_physical_coord works with embedded elements)

[lijas]

This might be more difficult than i first thought. For the embedded case, the function_value_from_physical_coord would need some signed distance search to find the parametric coord on the surface

[lijas]

I think this type of construction of SArrays is fine, but I dont know how efficiently the compiler can optimize it, maybe someone else can review this part :)

[henrij22]

You were right, there are a lot of calls in there.
This implementation takes 24 ns on my system, I found a better version which only takes about 2 ns

[henrij22]

Thanks @lijas for the review. I will adress the rest of your points when I have more time :)

[lijas]

I compared the results from function_hessian from the hessian computed with AD:

ip = Lagrange{RefQuadrilateral, 2}()
qr = QuadratureRule{RefQuadrilateral}(2)
cv = CellValues(qr, ip, ip^3; update_hessians = true)

coords = [Vec{3}((2x[1], 0.5x[2], rand())) for x in Ferrite.reference_coordinates(ip)]
coords_nl = [x + rand(x) * 0.01 for x in coords] # Create non-linear geometry
reinit!(cv, coords_nl)
ue = [rand()*0.01 for i in 1:getnbasefunctions(cv)]

for i in 1:getnquadpoints(cv)
    ξqp = cv.qr.points[i]
    xqp = spatial_coordinate(cv, i, coords_nl)
    Hqp, Gqp, Vqp = Tensors.hessian(x -> function_value_from_physical_coord(ip, coords_nl, x, ue_nl), xqp, :all)
    @test function_value(cv, i, ue_nl) ≈ Vqp
    @test function_gradient(cv, i, ue_nl) ≈ Gqp
    @test Ferrite.function_hessian(cv, i, ue_nl) ≈ Hqp
end

The function value and gradient test passes, but the hessian does not. So I suspect that there is some error in you hessian computation.

For this test to work, I needed to change this function to

_solve_helper(A::SMatrix{idim, odim}, b::Vec{idim, T}) where {odim, idim, T} = Vec{odim, T}( inv(A'*A)*(A'*b)) ## pinv(A) * b does not work with ForwardDiff

[henrij22]

@lijas Thanks, okay then I have to take a second look. It's a bit tricky to see where the error could be, its either in the computation of the Hessian (of the geometry) or in the apply_mapping() function.
I think we should have the AD test also somewhere in the test suite.

Edit: I think the error should be in apply_mapping()

[henrij22]

@lijas So my idea to implement this was to use the existing derivation from https://ferrite-fem.github.io/Ferrite.jl/stable/topics/FEValues/#mapping_theory
and just calculate it with SMatrices. However there seems to be an issue that I can't quite solve.
I tried it with einstein summation

@tullio d2Ndx2[i, j] := Jinv[s, j] * d2Ndξ2[s, r] * Jinv[r, i]
@tullio d2Ndx2[i, j] += -dNdξ[r] * (Jinv[r, k] * H[k, p, s] * Jinv[p, i]) * Jinv[s, j]

However it does not pass your test ...
Edit: in here

[lijas]

Hmm, i dont think the equations for the hessian for solids is directly transfereble to embedded elements (the first part in your code snippet maybe is, but not the the second part).

Can I ask, what do you need the embedded hessian for? :)

[henrij22]

Yeah I figured ... unfortuntly :(
Sure, I was implementing shell elements where I needed the second derivative. I can do with only using the Hessian itself, as I do the geometry mapping myself on element level. However I thought if I add the Hessian to Ferrite I can also tackle the geometry mapping in a general way.
If I have some time, I will take a look at the literature if I can find anything on how to do this correctly.

[lijas]

Cool, I guess you are using Kirchhoff-Love theory then if you need the Hessian? How are you getting C^1 continuity?

For fun, I was "vibe-mathing" with ChatGPT, and I think it was able to derive the hessian. Perhaps an approach you can take if you dont find any good reference in the litterature haha.

[henrij22]

@lijas Haha good idea.
Yes, its for KL-shells and shell theories building on top of that. I was using isogeometric shape functions.
Maybe if you are interested more, you can also drop me an email (my GH profile should have it) or DM, then we don't spam this thread haha

[KnutAM]

Maybe if you are interested more, you can also drop me an email (my GH profile should have it) or DM, then we don't spam this thread haha

IMO please keep spamming :D (Or make a separate discussion/issue for shell implementations as this is something the current support is not very good for). These discussions are often super-useful for anyone wanting to implement these things in the future!

[henrij22]

Perfect :D
Thats true, I think for conitnuum elements, Ferrite is excellent, but for embedded elements, it needs some work-arounds.

[KnutAM]

Linking also for reference : Ferrite-FEM/Tensors.jl#188

Just noting that I'm currently trying to push this one over the finish line - so if you want to try stuff with that functionality do feel free to and happy to get feedback on this one!

[henrij22]

@KnutAM thanks for the heads up. The mixed tensors look great!
However I don't really know how to procede here. I've now implemented a transformation I think might be correct. At least for non-curved quadrilaterals I can verify it with the test lijas proposed (not yet pushed), but for curved ones, the test fails, although I am not sure how to test this.

[KnutAM]

I haven't worked much with shells, but I think some key issues (that you might already have hit/solved). If you already have derived things in this direction that you could share, that would be nice and perhaps we can push it in the right direction. I think it could be useful to see the complete formulation that you are aiming at (not necessarily in research but say the KL-formulation for curved shells):

• I think the theory assumes a curvelinear local coordinate system so standard index notation might not work (was googling a bit and found this master thesis https://publications.lib.chalmers.se/records/fulltext/250355/250355.pdf)
• The assumption in the derivations in the docs of diffeomorphic doesn't hold, but I believe this is covered by considering the proper covariant/contravariant bases.

[lijas]

I think the master thesis you linked to uses the gradient/hessian w.r.t surface coordiantes, e.g dN/ds where s is a coordinates along the surface (dim(s) = rdim). This PR implements dN/dX where dim(X) = sdim.

In my experience, shell formulation are often expressed and implemented in terms of dN/ds because the constituve law often requires it, so it would be nice to have this aswell. But I dont know if it is easy to put this in to CellValues directly, or if a new EmbeddedCellValues would be easier.

[KnutAM]

I think the master thesis you linked to uses the gradient/hessian w.r.t surface coordiantes, e.g dN/ds where s is a coordinates along the surface (dim(s) = rdim). This PR implements dN/dX where dim(X) = sdim.

Good point! I'm unsure though how we add in that $\frac{\partial N}{\partial \boldsymbol{x}} \cdot \boldsymbol{n} := 0$ (where $\boldsymbol{n}$ is the surface normal) in the derivations without going through the surface coordinates? AFAIU we only have that $\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{x}} \cdot \boldsymbol{v} = \boldsymbol{v}$ if $\boldsymbol{v} \perp \boldsymbol{n}$? Since the function $\boldsymbol{x}(\boldsymbol{\xi})$ is not unique if we allow variations in $\boldsymbol{x}$ in the normal direction...?

[henrij22]

Hi,
sorry for the late reply. Shells are not my main focus of research, I am mostly using solids. However I wanted to see how well one could implement shell elements in Ferrite. I have implemented the KL shell from [1] and derived formulations for example in [2]. For this formulation we transform the derivatives with the directors and its derivatives in the element formulation itself. So it was enough to not transform the derivatives and just dNd\xi etc from the function values.
I am currently working on implementing a Ferrite Interpolation and Grid based on https://github.com/henrij22/TinyGismo.jl

[1] J. Kiendl, K.-U. Bletzinger, J. Linhard, und R. Wüchner, „Isogeometric shell analysis with Kirchhoff–Love elements“, Computer Methods in Applied Mechanics and Engineering, Bd. 198, Nr. 49–52, S. 3902–3914, Nov. 2009, doi: 10.1016/j.cma.2009.08.013.
[2] B. Oesterle, E. Ramm, und M. Bischoff, „A shear deformable, rotation-free isogeometric shell formulation“, Computer Methods in Applied Mechanics and Engineering, Bd. 307, S. 235–255, Aug. 2016, doi: 10.1016/j.cma.2016.04.015.