Projecting on ElementDG

[HelgeGehring]

Regarding https://github.com/kinnala/scikit-fem/discussions/895#discussioncomment-4044517 : I've made a minimal example showing a problem projecting on ElementDG:

import numpy as np
import matplotlib.pyplot as plt

from skfem import *

fig, axs = plt.subplots(1, 3, subplot_kw={'aspect': 1})

mesh = MeshTri().refined(4)

basis = Basis(mesh, ElementVector(ElementTriP1()))
x = basis.project(lambda x: np.array((x[1], -x[0])))
basis.plot(x, ax=axs[0])

basis2 = basis.with_element(ElementVector(ElementTriP1()))
x2 = basis2.project(basis.interpolate(x))
(x2_x, basis2_x), (x2_y, basis2_y) = basis2.split(x2)
basis2_x.plot(x2_x, ax=axs[1], shading='gouraud')

basis3 = basis.with_element(ElementVector(ElementDG(ElementTriP1())))
x3 = basis3.project(basis.interpolate(x))
(x3_x, basis3_x), (x3_y, basis3_y) = basis3.split(x3)
basis3_x.plot(x3_x, ax=axs[2], shading='gouraud')

plt.show()

[kinnala]

Thanks! Perhaps it's related to #829.

[kinnala]

Yep, I think it's the same issue as #829. Projection is just fine, split method just does the splitting wrong.

[kinnala]

Now I see that the split logic is completely wrong for ElementVector if there are multiple DOF's related to the same topological entity.

[kinnala]

This looks better:

[kinnala]

Can you confirm that it works now?

[HelgeGehring]

Looks way better!

But for my application it doesn't seem to work so good, it seems as if some elements are mirrored

[HelgeGehring]

import numpy as np
import matplotlib.pyplot as plt

from skfem import *

fig, axs = plt.subplots(1, 3, subplot_kw={'aspect': 1})

mesh = MeshTri().refined(4)

basis = Basis(mesh, ElementTriN0())
x = basis.project(lambda x: np.array((x[1], -x[0]*0)))
basis.plot(x, ax=axs[0])

basis2 = basis.with_element(ElementVector(ElementTriP1()))
x2 = basis2.project(basis.interpolate(x))
(x2_x, basis2_x), (x2_y, basis2_y) = basis2.split(x2)
basis2_x.plot(x2_x, ax=axs[1], shading='gouraud')

basis3 = basis.with_element(ElementVector(ElementDG(ElementTriP1())))
x3 = basis3.project(basis.interpolate(x))
(x3_x, basis3_x), (x3_y, basis3_y) = basis3.split(x3)
basis3_x.plot(x3_x, ax=axs[2], shading='gouraud')

plt.show()

[kinnala]

Are you sure this simply isn't how the field looks like? AFAIK ElementTriN0 is not C0 continuous and only tangential components are continuous over the edges.

[kinnala]

I mean, by projecting to ElementTriP1 you are forcing the x-component to be C0 continuous which is not necessarily the case in ElementTriN0?

[HelgeGehring]

From my understanding ElementTriN0 doesn't force the normal components to be continuous, but is also okay with that :D So, for this example I don't see any reason why the ElementTriN0 should do something not continuous as I project something continuous on it.

From the physical point: If there's an interface, the field will have a discontinuity there and I also want to see the discontinuity in the colorpolot. But if there's not discontinuity I don't see any reason why the field should have a discontinuity.

Thus, I have as I expect a problem at the interface (for the mode simulations, not for this simplified example) if I project on ElementTriP1(), therefore the hope is that ElementDG(ElementTriP1()) let's things be continuous within the element and also in the areas where no abrupt discontinuity is.

[HelgeGehring]

Also, from the physical point of view: There's no reason why the field should be aware of the discretization. If there's a change of epsilon or something - nice, I'm happy about discontinuities in the plot. But if nothing changes in the environment, the border of the element should also not be visible in the plot.

[kinnala]

I'm just saying that if you look at the definition of the basis functions it's rather obvious that it won't be continuous. Like piecewise constant basis, no matter what continuous function you start with, the result won't be continuous because of the basis.

[HelgeGehring]

Okay, then It's probably the wrong approach for me, and this issue is solved :) Would there be a (plotable) way to let the Elements be continuous within subdomains, but disconnected at the boundaries?

[HelgeGehring]

Or would it be solved by using a higher order Nedelec Element? https://defelement.com/elements/examples/triangle-Nedelec-2.html

[kinnala]

ElementTriN1 would most likely fix your issue because it will increase the accuracy of the method significantly.

However, implementing such element is not as easy as it sounds. It's easy to copy paste those basis functions given in the reference but since this H(curl) is still slightly unexplored territory in scikit-fem, there are no examples to follow and it requires some thinking to get them written in the correct order while some of them must be multiplied by -1 so that the orientation is correct. As you might imagine, 8 basis functions in correct order leads already to quite many ways it can go wrong if done randomly. :-)

Thus, another option is to just refine more.

[HelgeGehring]

ElementTriiN1 would probably also increase the accuracy of the simulation itself, right? I'm still unsure about the trade-offs between more elements and higher order elements, but I'd guess that ElementTriiN1 would also help a lot for the simulation as the field is not constant within the element?

I can feel the pain of implementing it :D But the question would be: is there a good way to test the implementation? But here I don't have (yet :) ) the understanding how to do that - would it be an option to just project a large number of functions on the basis and then compare the values on the quadrature points? Using for example sympy we could also automatically differentiate those equations and compare also curl/... 🤔 I'd guess if we have a good way of testing them automatically extensively, implementing those elements wouldn't be a pain?

[kinnala]

I just implemented ElementTriN1 and accidentally pushed to master -_-. I was lucky and it seems to work on the first try.

[kinnala]

We have more automatic testing for other element types but ElementHcurl does not have an automatic convergence test and other tests such as nodality of the DOFs.

[kinnala]

BTW I think I will rename these elements as follows: ElementTriN0 -> ElementTriN1 and ElementTriN1 -> ElementTriN2 so that the number matches the polynomial degree, similar to what we have for other elements.

[HelgeGehring]

Uijj, cool, thanks! Looking forward to try it!

The renaming sounds good, I was already wondering why it didn't match the degrees in defelement.com and thought there was some other naming scheme behind it.

[HelgeGehring]

I'm still wondering if it would somehow make sense to have elements which are only within the subdomains continuous for this interpolation. Would that be just a question of implementing it, or is there a fundamental reason why this is a bad idea?

[kinnala]

It's certainly possible. However, this library has been written, due to simpler vectorization, assuming that every finite element is the same. Removing continuity only for some edges is against this assumption and probably more easily done by splitting the mesh to multiple Mesh objects and treating them separately and not as subdomains.

[HelgeGehring]

Okay, thanks!