supermesh projection for KMV elements

[santos-td]

Dear all,

I am trying to project (supermesh projection) a solution obtained with a mesh (m1) onto another mesh (m2). I use KMV elements (Kong-Mulder-Veldhuizen, which are continuous with additional internal bubble functions). The projection works well for KMV P1, P3, P4, and P5 (max order of KMV elements), which means that f1 (a function defined over m1) is numerically equal to f2 (the function projected on m2). However, for KMV P2, f2 seems odd, with error (f1-f2) as high as 30%, at least for the code below.

Would you have any idea about the source of this error so I could try to fix it? Or maybe I am missing something and the projection will only work for Lagrange elements? Any idea would be helpful.

Thank you in advance
Thiago

from firedrake import *

m1 = RectangleMesh(15, 15, 1, 1)
m2 = RectangleMesh(10, 10, 1, 1)

P = 2

V1 = FunctionSpace(m1, 'KMV', P)
V2 = FunctionSpace(m2, 'KMV', P)

f1 = Function(V1).interpolate(Constant(1))
f2 = Function(V2)

#f2 = Projector(f1, V2).project()
f2.project(f1)

File("f1.pvd").write(f1)
File("f2.pvd").write(f2)

[santos-td]

I found out the source of this issue. It is necessary to force the integration rule to be KMV too in both assembly calls, rhs and mixed mass matrix, e.g. (supermeshing.py:194):

qr_x = finat.quadrature.make_quadrature(
V_S_A.finat_element.cell,
V_S_A.ufl_element().degree(),
"KMV"
)
M_SS = assemble(inner(TrialFunction(V_S_A), TestFunction(V_S_B)) * dx(rule=qr_x))

[ksagiyam]

Relevant Slack conversation: https://firedrakeproject.slack.com/archives/C1Q0Y6H8A/p1652626816762209

So it seems that UFL polynomial degree estimate somehow needs to be amended for KMV P2 as it contains a P3 basis.
I think the fix you made is correct and good if we are projecting from KMV to KMV, but I think the degree estimate needs to be fixed for general purposes.

[santos-td]

Hi, thanks for replying. I agree, the degree estimate should consider the degree of the bubble, with is > than the "labeled" degree of the element.

[ksagiyam]

Putting aside degree estimate, I think it is still good to use the designated quadrature rules at

firedrake/firedrake/projection.py

Line 138 in 5c10b51

a = firedrake.inner(u, v)*firedrake.dx

if the target space was KMV as you have done, which indirectly fixes half of the issue, so if you could make a pull request, I'm happy to look at it. I'm not sure if it is ok to use the quad rule for the rhs.

[santos-td]

Thanks for looking at it. I believe there is two alternatives to fix that aside the degree estimate: changing the "degree" label of the KMV elements (please, see my reply to Patrick's comment below); or modifying "project" such that a "rule" or "degree" could be set by the user, e.g., f2.project(f1, degree=5).
Not sure which one is the best, so if you have any thoughts, that would be great

[pefarrell]

My understanding of the "degree" of a finite element is the degree of the minimal complete polynomial space containing the function space of the element. So if a KMV function is (some special kind of) a cubic polynomial, its degree should be 3 (in my opinion). The quadrature estimator would then do the right thing.

[santos-td]

Yea, that is a good point, and I suspect it could be a solution (i.e., shifting the "degree label" of the KMV elements by 1, at least in the degree() method that is called by dx).

Originally the KMV element (or mass lumped element) had a "degree" and an "order". Please, see the figure below. So the implementation of such element in FIAT seems "consistent" with such paper (not sure if consistent from the mathematical point of view): in its constructor , the order is computed considering the degree + bubble. For instance, KMV P2 has degree=2 and order=3.

But then, fiat.CiarletElement.degree() returns self.poly_set.get_embedded_degree() (which is 3 for the case of KMV P2). So, it seems "embedded degree" is the "right one" for KMV.

However, ufl.FiniteElementBase.degree() returns self._degree. And actually there is a "FIXME" there saying that "embedded degree" could be used instead.

So I suspect that the quadrature estimator uses "degree()" from the ufl and not from the fiat (I am still learning how firedrake works in detail actually).

So, a possible fix would be inserting an alias in ufl.FiniteElement such that "degree()" will return the correct degree for KMV, but not sure if this is safe, actually...
Another possibility is modifying "project" such that a "degree" or "rule" could be defined externally by the user, but not sure if it is the best way...

(figure from Journal of Computational Acoustics, Vol. 9, No. 2 (2001) 671-680)

[dham]

I think this is in effect a bug in UFL. I believe that the usual convention is to number finite elements by the space that they span, rather than the space that is needed to span them. However, because UFL uses the degree to estimate quadrature, the latter is needed. This already occurs in UFL for RT elements. Relabelling the element is probably the best workaround for the moment. When we redo UFL under SysGenX we should give elements a degree and a separate required quadrature degree.

[santos-td]

Hi, thanks for looking into it and clarifying those points. I will perform some tests after relabeling the KMV elements, and then make a pull request.

[rckirby]

Relabeling the degree is problematic since it will throw off the degree estimation algorithm in UFL. It’s best to let things ride and just use a custom quadrature when you need the diagonal mass matrix. You can replace project with a call to solve on a very simple L^2 inner product or, even better, interpolation since these coincide for KMV elements with KMV quadrature.

…

On May 20, 2022, at 7:16 AM, santos-td ***@***.***> wrote:  Hi, thanks for looking into it and clarifying those points. I will perform some tests after relabeling the KMV elements, and then make a pull request. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.

[santos-td]

Oh I see, many thanks for this. So, I will stick to an overloaded projection such that I can change the quadrature. Thanks!

[pefarrell]

Err, isn't relabelling the degree exactly what's necessary to get the degree estimation algorithm in UFL to work? (Maybe I'm missing something here.)

The basis functions are e.g. cubic polynomials; they need a quadrature rule to match. In order for UFL to pick the right quadrature rule, we have to tell it that they're cubic.

[dham]

I think @pefarrell is right. This is a wart that we should really fix in a future major version of UFL, but for now it would be right to renumber the element.