Implementing a new Finite Element with Rotational degrees of Freedom

[Lukas Johannimloh]

Hello everyone,

As part of my Master's Thesis I need to implement a Triangular Element with 2 translational and 1 rotational degree of freedom at each corner, so 9 DoFs per element

The element is published and explained in:

A family of triangular and tetrahedral elements with rotational degrees of freedom
Christian Bucher, 2017
https://www.researchgate.net/publication/321196381_A_family_of_triangular_and_tetrahedral_elements_with_rotational_degrees_of_freedom

I created a script where I implemented a CST (via FESystem with two SimplexP) to make sure the mesh import, solver, and export of the results to ParaView work. The next plan was to then implement the Bucher triangle and use it in the existing script.

My current problem is that I don't know how to do that, 

The first idea was to implement only the rotational DoF and then combine it with two SimplexP elements, but my current understanding is that this would not be possible with FESystem due to the interaction between the rotation and translation 
I spent the last few days reading a lot of documentation and trying to understand how other elements work. Sadly, I am slowly running out of time and cannot continue like this until it eventually works, so I need some help

When I asked AI, it told me that I only need a header file, but I am quite skeptical since all the other elements also have .cc and .inst,in files, so I don't know why that should be not neccesarry for me?

What files should I take a look at, or how should I proceed to make progress again?

Any other things to consider or look at?

Thanks a lot beforehand, and I hope I'm not wasting your time with stupid questions. In a month from now, I will be far deeper into the topic but I don't have time to wait untill then

[Wolfgang Bangerth]

On 7/22/26 05:58, Lukas Johannimloh wrote:
>
> My current problem is that I don't know how to do that,
>
> The first idea was to implement only the rotational DoF and then combine it
> with two SimplexP elements, but my current understanding is that this would
> not be possible with FESystem due to the interaction between the rotation and
> translation
> I spent the last few days reading a lot of documentation and trying to
> understand how other elements work. Sadly, I am slowly running out of time and
> cannot continue like this until it eventually works, so I need some help
>
> When I asked AI, it told me that I only need a header file, but I am quite
> skeptical since all the other elements also have .cc and .inst,in files, so I
> don't know why that should be not neccesarry for me?
>
> What files should I take a look at, or how should I proceed to make progress
> again?

Whether you put things into a header file only, or split into .h/.cc files is
a secondary consideration. The starting point for implementing finite elements
is always this:
* How do the shape functions look like on each cell? For example, are they
polynomials and if so what are they?
* How are shape functions transformed from the reference cell to the real cell?
* What is the connectivity between the shape functions on neighboring cells.
For example, are they continuous? (That's easy.) Do they need to be continuous
differentiable? (That's hard.)

If you could summarize the answers to these questions, we might be able to
point you to an element that is similar and already implemented. That could
then serve as the basis for your own work.

Best
W.

[Lukas Johannimloh]

First of all, thanks a lot,

* How do the shape functions look like on each cell? For example, are they
polynomials and if so what are they?

The shapefunctions are based on areal coordinates:

λₖ = Aₖ / A

The shape function consists of the linearly intertpolated displacement of the corners

uᵤ^(k) = λₖ,       uᵥ^(k) = 0,

vᵤ^(k) = 0,         vᵥ^(k) = λₖ.


and the displacement due to the rotation of the corners

uᵣ^(k)(x, y) = −(y − y^(k)) λₖ²

vᵣ^(k)(x, y) = (x − x^(k)) λₖ²

therefore the total displacement for a random point within the element is

u(x, y) = Σₖ= λₖ * u^(k) − (y − y^(k)) λₖ² * r^(k),

v(x, y) = Σₖ= λₖ * v^(k) + (x − x^(k)) λₖ² * r^(k).

This results in an cubic interpolation of the shapefunction

* How are shape functions transformed from the reference cell to the real cell?

Not comletely sure if I understand the question correct, but I would say that the shape functions are transformed due to the Areal coordinates since the coordinates of the corners of the real cell are within it

λ₁(x, y) = (x₂y₃ − x₃y₂ + x y₂ − x y₃ + x₃y − x₂y) / (2A),

λ₂(x, y) = (−x₁y₃ + x₃y₁ + x y₃ − x y₁ − x₃y + x₁y) / (2A),

λ₃(x, y) = 1 − λ₁ − λ₂,

where 2A = |x₂(y₃ − y₁) + x₁(y₂ − y₃) + x₃(y₁ − y₂)|.

* What is the connectivity between the shape functions on neighboring cells.
For example, are they continuous? (That's easy.) Do they need to be continuous
differentiable? (That's hard.)

the corners for the neighboring cells are continuous C⁰, but as far as I am aware not continuously differentiable C¹

Thanks a lot again, hope it helps
Lukas