Implementing new finite element in FreeFEM

General Discussion
1

Hello all,

I would like to have some feedback about implementing new finite element in FreeFEM.
How difficult and how long it is to implement new finite element (in particular non conforming) in FreeFEM ?

For anyone who has already done (or tried it), feel free to share your feedback (the type of finite element, the difficulties you have faced with, if it works well or not, …)

Thank you in advance,

Best regards,

Loïc,

2

This depend strongly of DoF of the finite element,
but generally it is not so difficult.

I can help you if you need.

3

Hello @frederichecht

Thank you for you reply,

I would like to implement a non conforming Finite Element which is very popular (in the fluid dynamics domain).

I have explained this finite element in the attached document. I would like to implement the case n=1 and the case n=2. (The case n=2 in priority).

What do you think of the difficulty of this ?

Best regards,

Loïc,
Description_FE_space_FreeFEM.pdf (114.9 KB)

4

it is simple do but a little technical, I will ty to make a canvas.

5

Thank you for your help @frederichecht .

I will try to understand and follow this tutorial: Developers

If you need more details about this finite element, don’t hesitate to let me know.

Best regards,

Loïc,

6

Remark the version P2 + bulle exist

P2pnc piecewise quadratic plus a bubble P3 element with the continuity of the 2 moments on each edge (version 3.59) (need load “Element_P2pnc” is exactly you element in case 1 .

So you can start form this and call it P3pnc for exemple.

7

Thank you,

In the file Element_P2pnc.cpp , I can see the 7 basis functions.
But I cannot understand where the DoF are defined.

Could you help me to understand how the DoF are defined ?

Best regards,

Loïc,

8

you can see in the test testFE-P2pnc.edp
the 2 dof on edge E = (i,j) are
\int_E \lambda_i v and \int_E \lambda_j v where \lambda_i and \lambda_j are barycentric
coordinate

9

I think in the FE definition, these doF are defined there:

for (int i = 0; i < 3; i++) {
  for (int p = 0; p < QFE.n; ++p) {
    R l1 = QFE[p].x, l0 = 1 - QFE[p].x;
    if (oe[i] < 0) {
      swap(l0, l1);
    }

    if (ddd < 3) {
      cout << p << " " << oe[i] << " " << l0 << " " << l1 << endl;
    }

    R p0 = l0, p1 = l1;
    R cc1 = p0 * QFE[p].a;
    R cc0 = p1 * QFE[p].a;
    v[k++] = cc0;
    v[k++] = cc1;
  }
}

for (int p = 0; p < QFK.n; ++p) {
  double w = QFK[p].a;
  R l1 = QFK[p].x, l2 = QFK[p].y, l0 = 1 - l1 - l2;
  R b = 1;
  v[k++] = b * w;
}

ffassert(k == this->pij_alpha.N( ));

}

But in my case I want that my DoF are \int_E v L_k where L_k are the k-th Legendre polynomial. Thus in case n=1, the 2 polynomial should be 1 and \lambda.

I think this is not what P2Nc does, am I wrong ?

Best regards,

Loïc,

10

Yes the degree of freedom are not exactly the same but a the end the space are exactly the same because we have same constraint;
But it is trivial to change

11

Yes sure,

And for example in the case n=2,
we can define the three weight just as P0=1, P1=l_0 and P2= \frac{1}{2} (3 l_0^2-1) ?

And I think I have to change the quadrature formula for the case n=2?

I have also few questions about the code in Element_P2pnc.cpp,

  • What do QFE[p].x, QFE[p].y and QFE[p].a mean ?
  • How the double C1[] is defined and what it does ?
  • What do K.H(i) and K[i] mean ?
  • And what is the function IPJ() ?

Thank you in advance,

Best regards,

Loïc,

12

Bonjour,

le plus simple est de faire une visio,
via zoom , skype, or google meat !

envoyer moi votre mail pour vous envoyer une invitation

mon mail est : Hecht Frédéric <Frederic.Hecht@upmc.fr>

1 Like