FacetBasis.trace()

[kinnala]

This aims to address #514 by adding a helper method to FacetBasis. FacetBasis.trace will create a trace mesh and restrict a solution vector to that trace mesh. Optionally it performs L2 projection onto a new finite element basis on the trace mesh.

Because we do not have Mesh classes defined for embedded meshes (surfaces etc.), FacetBasis.trace will return a pair (p, t), i.e. an array of nodes and their connectivity. It is up to the user to, e.g., save them to an external format or project them to lower dimensional space so that we can initialize one of the existing Mesh classes.

Here is an usage example:

from skfem import *
from skfem.visuals.matplotlib import *
m = MeshQuad()
m.refine(3)
basis = FacetBasis(m, ElementQuad1(), facets=m.facets_satisfying(lambda x: x[1]==0.0))
p, t, y = basis.trace(m.p[0], ElementQuad0())
plot(InteriorBasis(MeshLine(p[0], t), ElementLineP0()), y)
show()

Fixes #514

[kinnala]

@gdmcbain Do you feel it addresses #514 ?

[gdmcbain]

Oh, no. Sorry, I don't think I'll be able to review this before next year. It looks excellent, reading through the code. I'm very pleased at how concise it is; that's a good testament to the design of what's gone before, I think. I won't have time to test it before the break.

Would it be cleaner to omit the optional preliminary projection? The user can just do that themselves if they need to, can't they? Actually I envisage projecting afterwards, which should be much cheaper since the trace is smaller than the domain. In the motivating application, I solve the Navier–Stokes equation with quadrilateral Taylor–Hood elements and want to (i) extract the velocity along a part of the boundary, (ii) project it to piecewise constant for use as an inlet velocity boundary condition in a finite volume code (OpenFOAM). Projecting the entire field to ElementVectorH1(ElementQuad0) seems wasteful, but this might be premature optimisation.

[kinnala]

Would it be cleaner to omit the optional preliminary projection? The user can just do that themselves if they need to, can't they?

Sure, but note that many times there is no trivially corresponding boundary element, e.g., in the ElementTriMorley example, and then projection may be required.

Actually I envisage projecting afterwards, which should be much cheaper since the trace is smaller than the domain.

If you read carefully, the L2 projection is performed only on the boundary.

In the motivating application, I solve the Navier–Stokes equation with quadrilateral Taylor–Hood elements and want to (i) extract the velocity along a part of the boundary, (ii) project it to piecewise constant for use as an inlet velocity boundary condition in a finite volume code (OpenFOAM). Projecting the entire field to ElementVectorH1(ElementQuad0) seems wasteful, but this might be premature optimisation.

See above.

[kinnala]

In principle, each Element could carry a reference to the "best possible" boundary element, e.g., ElementTriMorley would reference ElementTriDG(ElementTriP2()) which would transform to something like ElementLineDG(ElementLineP2()) (I think?). Thus, one wouldn't be required to supply this optional argument even in complicated cases. What do you think about that?

[gdmcbain]

Thanks for the further comment on Morley. (Actually the use of Morley in the example was accidental; it was just that ex02 was the first example with a nonzero trace and it uses Morley.) I was puzzled at the idea that ‘there is no trivially corresponding boundary element’ to the Morley space but am not really too familiar with the element. My successful use of it in ex18 was just by analogy: it works for the clamped plate therefore it works for the creeping streamfunction since they're governed by the same biharmonic equation and boundary conditions—but those do involve the trace… So I went looking for a reference and found a very good one on scicomp.stackexchange.com. Thanks!

I do like the idea that

each Element could carry a reference to the "best possible" boundary element,

or just a NotImplemented for exotic cases, either to be filled in or avoided by preliminary projection.

And on

If you read carefully, the L2 projection is performed only on the boundary.

Sorry, yes, this is indeed obvious.

I'll get back to this properly in the new year but thanks again meanwhile for the explanations of Morley.

[kinnala]

I concluded that this will not work if the provided elem has more than one interior DOF. In such a case interpolation/linear combination of basis functions is required to restrict to the boundary. I added an exception which tries to detect such cases.

[kinnala]

I made a major change to how trace was implemented in the previous two commits. FacetBasis.trace will return InteriorBasis (instead of p and t) and the projected solution vector.

User provides a function which defines how the points on the boundary of a d-dimensional mesh are projected onto d-1-dimensional space. This is used to initialize the d-1 dimensional mesh.

User can also request a "target element" to be used in the trace mesh. Not all target elements are supported because there may not always be a one-to-one correspondence between facet DOFs of d-dimensional mesh and DOFs of d-1-dimensional mesh.

[gdmcbain]

O. K. The example from 2020-12-23 is modified to:

from skfem import *
from skfem.visuals.matplotlib import *

m = MeshQuad()
m.refine(3)
basis = FacetBasis(m, ElementQuad1(), facets=m.facets_satisfying(lambda x: x[1] == 0.0))
tb, y = basis.trace(m.p[0], lambda x: x[0])
plot(tb, y)
show()

and produces

which looks good. I'll try it out on the motivating OpenFOAM example and report back. I won't worry about ex02 as it's not immediately obvious what to do with Morley elements and they're not of immediate interest.

[kinnala]

What kind of discrepancy? The trace values are different?

[gdmcbain]

Yes, slightly different values, about 10%, which is weird. Sorry for such an incomplete report. That was just on closing yesterday. I'm sure the error is on my side but today I'll craft a minimal working example; yesterday's was part of an extended version of ex27, 660 lines of code with lots of irrelevance.

[gdmcbain]

Sorry, I haven't been able to reproduce the discrepancy in an MWE. I'm even more sure that the error isn't in the new FacetBasis.trace code but haven't quite been able to pin it down.

The MWE is something like

from skfem import *

import numpy as np


def f(y: np.ndarray) -> np.ndarray:
    return y * (1 - y)


mesh = MeshQuad.init_tensor(*np.tile(np.linspace(0, 1, 9), [2, 1]))
mesh.define_boundary("left", lambda x: x[0] == 0)
ib = InteriorBasis(mesh, ElementQuad2())
fb = FacetBasis(ib.mesh, ib.elem, facets=mesh.boundaries["left"])

y = fb.doflocs[1, fb.get_dofs(fb.find).nodal['u']]
left_mesh = MeshLine(y)
left_basis = InteriorBasis(left_mesh, ElementLineP0())

u_old = project(lambda y: f(y[0]), basis_to=left_basis)

_, u_new = fb.trace(lambda x: f(x[1]), lambda x: x[1])
print(np.linalg.norm(u_new - u_old))

except that that works perfectly!

I'm pretty sure the PR is fine. Sorry for this terrible bug report!

[gdmcbain]

O. K., I've finally found the cause of the discrepancy and as expected it was completely in the extraneous code. (I had inferred the limits of the patch of boundary from the coordinates of the points passed to fun in skfem.utils.project; wrong because they're not doflocs but quadrature points and usually don't include the termini!) Easily fixed in the extraneous code using, e.g.,

mesh.p[:, mesh.facets[:, mesh.boundaries["patch"]]]

or basis.doflocs.

I'll now look at the actual PR.

[gdmcbain]

Would it be more Pythonic to have the DEFAULT_TARGET as an attribute or property of self.mesh? And EAFP over LBYL? Define as NotImplemented or similar for Mesh then override in subclasses? The resulting default error message should be almost as obvious.

Could the key use self.mesh.brefdom?

{'line': ElementLineP0, 'tri': ElementTriP0, 'quad': ElementQuad0}

(and worry about wedges #411 &c. later…)