render ElementTriRT0 field

[gdmcbain]

I think after #823 I've solved my problem of calculatring a weak ElementTriRT0 gradient of an ElementTriP0 field like the indicator function of a subdomain, but I don't know how to sensibly plot the answer to see whether it makes any sense.

from skfem import *
from skfem.helpers import div, dot

import ex17


@BilinearForm
def vector_mass(sigma, tau, w):
    return dot(sigma, tau)


@LinearForm
def linf(tau, w):
    return - w["phi"] * div(tau)


mesh, conductivity = ex17.basis0.mesh.to_meshtri(ex17.conductivity)
basis = CellBasis(mesh, ElementTriRT0())

indicator = (conductivity == conductivity.max()).astype(float)
M = vector_mass.assemble(basis)
rhs = linf.assemble(basis, phi=basis.with_element(ElementTriP0()).interpolate(indicator))
sigma = solve(M, rhs)

There's an intriguing use of side-effects in a Functional by @gatling-nrl on Friday, which looks like it could be adapted to work (possibly with a special quadrature rule supported on the midpoints of the facets) but I wonder whether there's a less devious way.

[kinnala]

You should be able to do Basis.refinterp and then apply plt.quiver or similar. Should there be such a functionality in skfem.visuals.matplotlib?

[kinnala]

Would this make sense? I cannot now think clearly that what zeroing out these DOFs means.

x0[mref.boundary_nodes()] = 0
x1[mref.boundary_nodes()] = 0

draw(mesh).quiver(*mref.p, x0, x1, minshaft=1, minlength=0)
show()

[kinnala]

Oh yeah, in the current master you can do this to get rid of few lines:

basis1 = basis.with_element(ElementTriDG(ElementTriP1()))
sigmai = basis.interpolate(sigma)
Y = [basis1.project(sigmai.value[d]) for d in [0, 1]]

[gdmcbain]

I don't think ElementTriRT0 has unique value: the normal component is continuous but not the entire field.

Of course this is completely correct, that being the whole point of H(div) and RT0; sorry, my question was wrong or at least needs precising. What I wanted was to render the degrees of freedom in a low-level way so that I could see exactly what was in the discrete algebraic representation of the vector-field, so very like figure 3.2(a) of Brezzi & Fortin (1991), figure 35 of Braess (2001) or figure 1.9 of Ern & Guermond (2004), which show an arrow pointing normally out of each of the three sides of a triangle. This would be the analogue of the markers in the plot of an ElementLineP1 function in ex15.

In concentrating thus on the normal part and looking at these three cited figures, I ignored and then forgot the tangential part. This was in the context of #821, viz. finding a consistently oriented unit normal to an interface between subdomains, in which the tangential part is irrelevant and should tend to zero as the mesh is refined. (I hope, I haven't checked that.) In fact, Ern & Guermond (figure 1.22) show a typical example of the

global shape functions associated with the [lowest order] Raviart–Thomas…finite elements in dimension 2

and point out that

the normal…component of the Raviart–Thomas…global shape function is continuous across the interface, but sinced the triangles are not isosceles, the tangential…component is not antisymmetric.

So of course my suggestion yesterday that the double-valuedness at the facets is an artifact of interpolation is completely wrong. Thanks for pointing that out.

The complete script, including the most recent suggestions is in e1d8de5.

I'll have a think about whether plotting just the degree-of-freedom arrows is advisable or not.

Another idea: the RT0 can be cast as Lagrange finite elements (Ern & Guermond 2004, Remark 1.40(ii)), so plotting the arrows for those degrees of freedom might render nicely simultaneously the discrete representation and the interpolated field.

[gatling-nrl]

I wonder whether it would be better to use something more like @gatling-nrl 's clever Functional, referred to above (except that I still hesitate to abuse asm to perform this side-effect).

I've been away from this for a bit now, and coming back to it, I think the trick with the Functional is useful to know, but probably should be discouraged... it is very much a side effect, and without adequate comments in the code it would be quite confusing, I think. You may have already come to this conclusion... I'm still catching up on where things are.

On the other hand, it does really expose exactly how the data structures in skfem are organized and how to manipulate them. That is a very good thing. Nevertheless, example code really shouldn't be mischievous.

[gdmcbain]

Exactly, on all points.

The good news is that yesterday I worked out how to get at the feet (.global_coordinates()) and directions (.normals) of the normals to facets using a FacetBasis, in stepping through the assembling of your Functional and the construction of the default_parameters(). So the Functional was educational but won't be required.

I plan to tie this all together, plotting side-by-side yesterday's intrinsic DoF-based and the earlier refinterp representations. I think both are useful. I probably won't use the piecewise-constant ex17.conductivity field but perhaps something like a two-dimensional version of ex37 ‘Mixed Poisson equation’.

[gdmcbain]

Tangential idea: in the example, the quadrilateral mesh from ex17 was split into a triangular mesh

mesh, conductivity = ex17.basis0.mesh.to_meshtri(ex17.conductivity)

as we don't have ElementQuadRT0 yet. Would that be a worthwhile element? Reference: Brezzi & Fortin (1991, §III.3.2 ‘Rectangular approximations of H(div; K)’)

[kinnala]

Sure!

[gdmcbain]

O. K., i've taken this up in #831.

[gdmcbain]

That was implemented in #832, merged, and then released in 5.2.0; hence adapting the above without recourse to MeshTri:

[gdmcbain]

That was with

basisRT = basis0.with_element(ElementQuadRT0())

indicator = (conductivity == conductivity.max()).astype(float)
M = vector_mass.assemble(basisRT)
rhs = linf.assemble(basisRT, phi=basis0.interpolate(indicator))
sigma = solve(M, rhs)

# render via refinterp

basis1 = basis0.with_element(ElementQuadDG(ElementQuad1()))
sigmai = basisRT.interpolate(sigma)
Y = [basis1.project(sigmai.value[d]) for d in range(basis1.elem.dim)]
((mref, x0), (_, x1)) = [basis1.refinterp(y, 1) for y in Y]

ax = draw(basis1.mesh)
plot(basis1.mesh, indicator, ax=ax).quiver(*mref.p, x0, x1)

[gdmcbain]

The normals and their feet are attributes of a FacetBasis, so it's not necessary to make recourse to a false assembly:

interfacial_facets = np.intersect1d(
    *(basis0.mesh.t2f[:, s] for s in basis0.mesh.subdomains.values())
)
fbasis = FacetBasis(
    basis0.mesh,
    basis0.elem,
    facets=interfacial_facets,
    quadrature=(np.full((1, 1), 0.5), np.ones(1)),
)

plot(basis0.mesh, indicator, ax=draw(basis0.mesh)).quiver(
    *fbasis.global_coordinates().value.squeeze(),
    *fbasis.normals.value.squeeze() * sigma[interfacial_facets]
).axes.get_figure().savefig(Path(__file__).with_suffix(".png"))

[gatling-nrl]

Am I understanding ElementTriRT0 correctly, that it is a basis such that the divergence of functions projected in this space will always be in ElementTriP0?

And this is useful, for example, in an electromagnetics problem where div(electric_field) = piecewise_constant_charge_distribution?

And in this example electric_field would use ElementTriRT0 and piecewise_constant_charge_distribution would use ElementTriP0?

[gdmcbain]

Am I understanding ElementTriRT0 correctly, that it is a basis such that the divergence of functions projected in this space will always be in ElementTriP0?

Yes, exactly; as Raviart & Thomas (1977, §3) put it more generally: ‘for any integer k ≥ 0 we shall associate a space Q_K of vector-valued functions q_K ∈ H(div; K) such that : (i) div q is a polynomial of degree ≤ k ; (ii) the restriction of q ⋅ ν_K to any side K' of K is a polynomial of degree ≤ k.’ These days Q_k is usually written ℝ𝕋_k, although Arnold's Periodic Table of the Finite Elements appears to adopt a different convention calling our ℝ𝕋_k RT_{k + 1} (probably for a very good reason as his k is shared by other elements in the same column).

And this is useful, for example, in an electromagnetics problem where div(electric_field) = piecewise_constant_charge_distribution?

Yes, Raviart–Thomas elements are popular for the Poisson equation when the flux is of as much interest as the potential.

And in this example electric_field would use ElementTriRT0 and piecewise_constant_charge_distribution would use ElementTriP0?

I think you also need a ElementTriP0 potential the gradient of which is equal to the electric field. (Or is it minus the gradient?) This equation is integrated by parts so that the piecewise constant potential doesn't have to be differentiated; this brings in the divergence of the field.

This is what's done in ex37 ‘Mixed Poisson equation’ in a cube with uniform charge distribution and isopotential boundaries.

An interesting feature of the mixed formulation is that the natural and essential boundary conditions are swapped; whereas for Lagrangian Pk elements the homogeneous natural condition is insulation, for RT it's grounding. This results from the aforementioned integration by parts of the statement that the flux is the (negative?) gradient of the potential.

(Other quibble besides the sign: I'm conflating field and flux. The latter probably involves some coefficient which multiplies the gradient. This is of course very important when that coefficient varies in space, as in ex17 ‘Insulated wire’.)

[gatling-nrl]

Over in #821 you said "calculated what I think should be the RT0 gradient of the P0 indicator"

I hadn't heard of Raviart-Thomas elements at the time and assumed it was going to be some kind of weak derivative of a piecewise constant function. But now I've learned a little more about it and realize I don't understand at all what you're doing. Somehow grad(scalar_function_in_P0)==vector_function_in_RT0? And this vector_function_in_RT0 is "well-behaved"?

[gdmcbain]

Yes, I think your original assumption is correct: in the usual application of Raviart–Thomas elements, the RT0 vector field is the (weak) gradient of the P0 scalar. This is exactly what I did here:

https://github.com/gdmcbain/scikit-fem/blob/e1d8de53183831c7d0dea17ab01aba41ae695fe7/docs/examples/ex43.py#L11-L29

In long-hand, something like:

• σ = grad φ
• (σ, τ) = (grad φ, τ), ∀ τ
• = [k n, τ]_∂Ω − (φ, div τ)

and in the example k = 0 in the boundary integral because we know the indicated subdomain is strictly interior. This leaves (σ, τ) = − (φ, div τ); hence

@BilinearForm
def vector_mass(sigma, tau, w):
    return dot(sigma, tau)

@LinearForm
def linf(tau, w):
    return -w["phi"] * div(tau)
...
M = vector_mass.assemble(basis)
rhs = linf.assemble(
    basis, phi=basis.with_element(ElementTriP0()).interpolate(indicator)
)
sigma = solve(M, rhs)

And this vector_function_in_RT0 is "well-behaved"?

Depends what you mean by ‘well-behaved’! The vector fields of RT0 do have some surprises, like the components tangential to facets not being continuous. This bit me above. But for sure they can be used in mixed discretizations of the Poisson equation, as in ex37. I think maybe in a way they have the very least amount of well-behavedness to make such a problem soluble and this has a certain appeal.

[gatling-nrl]

By the way, you can _^(ab)use get_dofs like this too...

def plot_query_points(x, ax, style, label):
    ax.plot(x[0], x[1], style, label=label)

plt.subplots(figsize=(5,5))
skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
basis_p1.get_dofs(facets=lambda x: plot_query_points(x, plt.gca(), 'or', 'facets'))
basis_p1.get_dofs(elements=lambda x: plot_query_points(x, plt.gca(), 'ob', 'elements'))
plt.legend()

lambda x: plt.plot is the printf of 2021!

[gdmcbain]

Ha! That is very very clever! Thank you.

Yes, on looking in the sources I see that get_dofs calls .mesh.facets_satisyfing or mesh.elements_satisfying which calculate the midpoints!

scikit-fem/skfem/mesh/mesh.py

Lines 299 to 300 in 2ad3763

	midp = [np.sum(self.p[itr, self.facets], axis=0) / self.facets.shape[0]
	for itr in range(self.dim())]

(But why doesn't this use mean? I've proposed #840 for that.)

[gatling-nrl]

I think these kinds of visualizations of what asm and get_dofs are doing is quite valuable and I want to work them into the docs, but I'm still not sure where to put this kind of mid-level under-the-hood stuff so newcomers (like myself) can find it. If we wrote an skfem book, it feels like... chapter 5? 1: Installation, hello fem! 2: Weak form, solving poissons equation. 3: visualizations of solutions, 4: saving, loading solutions, 5: skfem, a closer look at the api. 6: beyond ElementTriP1, all about elements, 7: beyond poissons equation, other diffeqs, mixed boundary conditions, etc, etc, etc.

[gdmcbain]

Very valuable yes.

That list of chapters is pretty good; about right, I think. What would I add? Looking back to when I started with scikit-fem, the things that I immediately had to add in order to do anything useful were:

• load meshes with tagged subdomains and boundaries (ex13)

scikit-fem/examples/ex13.py

Line 40 in 366392b

geom.add_physical_line(lines[-1], 'ground')

• postprocessing (ex12)

• inhomogeneous and particularly coordinate-dependent boundary data (ex14)

• 1-D problems (ex15) (conceptually easy, but a recurring source of corner-condition bugs)

• inhomogeneous and particularly coordinate-dependent coefficients (ex16)

• subdomain-dependent coefficients (ex17)

• mixed formulations, e.g. velocity-pressure for Stokes (ex18)

• unsteadiness (ex19)

I've started using quite a lot of finite element packages and really a lot of the introductions omit almost all of these (not the last two). They're kind of simple things but impossible if you don't see how to do it in the new package. Anyway most of these are solved now in scikti-fem :) though I might add to the list the consistent orientation of an interface #821; hopefully that's getting close.

But yeah maybe these would all go in ch. 7; that'd be reasonable.

[gdmcbain]

So I have written a script adapting ex37 from Tet to Quad (and flipping the sign of the source term for convenience) and plotted side-by-side the RT0 fluxes and refinterp-ed vectors. d866851/docs/examples/ex37quad.py.

I'm thinking of the flux density as minus the gradient but wanted to have the arrows pointing inwards for a tidier figure; hence the flip in sign of the solution. The centre is now a minimum, attracting flux, with the gradient pointing outwards.

It's working O. K. but doesn't look quite right in one respect: the RT0 DoF facet-fluxes look exactly as I would expect them to but some of the interpolated vectors aren't pointing the right way, I mean the ones that are one square in diagonally from the corners.

[kinnala]

It's working O. K. but doesn't look quite right in one respect: the RT0 DoF facet-fluxes look exactly as I would expect them to but some of the interpolated vectors aren't pointing the right way, I mean the ones that are one square in diagonally from the corners.

What do you think is the right way?

Because to me they look just fine.

[gdmcbain]

I mean they should all point inwards, towards the centre, in both plots, with the RT0 fluxes subject to the constraint that they're normal to facets. The refinterp-ed vectors mostly do this except for the four one diagonal in from the corners. It could be a plotting artifact; I haven't checked the numerical values of y; the arrows are small and smudgy but when I used show and zoomed in they didn't look right.

Ah, no, I see what's happened now! At the nodes of the MeshQuad the field is multivalued and what i saw as a smudgy arrow pointing outwards to the corner is actually a superposition of at least three arrows which might well be pointing inwards. Yes, I do believe they are. O. K., sorry about that. Thanks.