ex13: why is ElementTriP2() critical?

[gatling-nrl]

In ex13, if you change elements = ElementTriP2() to elements = ElementTriP1(), then conductance changes drastically:

elements = ElementTriP1() -> {'ground': -0.42608425018596274, 'positive': 0.45793999701984645}
elements = ElementTriP2() -> {'ground': -0.44141575677566447, 'positive': 0.44141575677563377}

I thought that, starting from ElementTriP1, instead of going to P2 I could decrease lcar. So I did, absurdly so:

lcar = .001

But still...

elements = ElementTriP1() -> {'ground': -0.4411121178299844, 'positive': 0.4414304277445507}

So things got better, I guess, but not by that much... and the computation time penalty was quite painful. (Also, I cannot even report P2 results for this extreme lcar because the python kernel couldn't cope when using juypter.)

So my question is... why is P2 critical in this case? Is it critical in all cases where flux measurements are the goal?

[gdmcbain]

How does the domain-integral calculation of conductance

scikit-fem/docs/examples/ex13.py

Line 51 in 71eb03e

conductance = {'skfem': u @ A @ u,

do with P1? Better? (I'm away from the terminal for a few days.) It should be more robust; see McBain et al (2018).

[gdmcbain]

As general advice, I would recommend the highest order elements available for problems with nice smooth solutions like this. I don't think the relative performance of P1 and P2 that you report is unusual for an analytic solution with all singularities outside the domain.

Note that (as later revealed in the study written up in 2018), this example is a bit lucky in its smothness in that the junction of Dirichlet and Neumann boundaries just so happens to have the Neumann boundaries lying along lines of current; were the angles other than 90°, there would be singularities at these corners. P1 might not compare so unfavourably then.