How do I move data from one basis/quadrature to another

[gatling-nrl]

cb_p1 = skfem.CellBasis(mesh, skfem.ElementTriP1())
cb_p0 = cb_p1.with_element(skfem.ElementTriP0())
ifb_p1 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP1())

sigma = cb_p0.zeros()
sigma[cb_p0.get_dofs(elements='a')] = 1
sigma[cb_p0.get_dofs(elements='b')] = 2
sigma[cb_p0.get_dofs(elements='c')] = 3

u = cb_p1.zeros()

@skfem.Functional
def func1(w):
    return w.u

@skfem.Functional
def func2(w):
    return w.sigma


func1.assemble(cb_p1, u=u)  # ok
func2.assemble(cb_p0, sigma=sigma) # ok

func1.assemble(ifb_p1, u=u) # ok, and half of what_i_want (see below)
# func2.assemble(ifb_p1, sigma=sigma)  # ValueError: Input array has wrong size.
# func2.assemble(ifb_p1, sigma=cb_p0.interpolate(sigma))  # ValueError: Quadrature mismatch
func2.assemble(ifb_p1, sigma=cb_p1.project(cb_p0.interpolate(sigma)))  # No error, but inexact sigma (see plot)

# also tried cb_dg = skfem.CellBasis(mesh, skfem.ElementDG(skfem.ElementTriP1()))
# in place of cb_p0, but cb_dg.zeros().shape != cb_p1.zeros().shape
# I guess because cb_dg and cb_p1 and different basises

fig, ax = plt.subplots(1, 2, figsize=(10,5))
skfem.visuals.matplotlib.plot(cb_p0, sigma, ax=ax[0], cmap='viridis', shading='gouraud')
skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
skfem.visuals.matplotlib.plot(cb_p1, cb_p1.project(cb_p0.interpolate(sigma)), ax=ax[1], cmap='viridis', shading='gouraud')
skfem.visuals.matplotlib.draw(mesh, ax=ax[1])

@skfem.Functional
def what_i_want(w):
    return dot(w.n, sigma*grad(u))

what_i_want.assemble(ifb_p1, sigma=sigma, u=u)  # ValueError: Input array has wrong size.

How can I get the piecewise constant sigma values into what_i_want?

[kinnala]

Try always doing an explicit call to Basis.interpolate() when passing kwargs to Form.assemble(). It will never work automagically if you combine two different bases.

[gatling-nrl]

what_i_want.assemble(ifb_p1, sigma=cb_p1.project(cb_p0.interpolate(sigma)))  # No error, but inexact sigma (see plot)

This was the closest that I could get, but it doesn't bring the exact values into the form. I understand why, but not how to fix it.

[kinnala]

Output of interpolate is DiscreteField and assemble kwargs should be of type DiscreteField unless you want to use the implicit interpolate inside assemble.

So it is something like this:

ifb_p1 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP1())
ifb_p0 = ifb_p1.with_element(skfem.ElementTriP1())
what_i_want.assemble(ifb_p1, sigma=ifb_p0.interpolate(sigma), u=ifb_p1.interpolate(u))

[gatling-nrl]

Oh, its obvious now. I didn't think about making an ifb_p0. That solves everything. Thanks!

[gatling-nrl]

@kinnala

Supposing I have some continuous f(x) that I write in python with

def f(x):
    return np.sin(x[0])  # for example

Then when I write

cb1 = CellBasis(mesh, ElementTriP1())
g = cb1.project(f)

g is the closest L2 approximation using the basis cb1.

Then when I write

h = cb1.interpolate(g)

h is samples of g (not f!), sampled at the quadrature points defined by cb1.

When I'm working with CellBasis and FacetBasis, they share the same mathematical basis, so I can treat the projection objects as belonging to either.

So I can write

cb1 = CellBasis(mesh, ElementTriP1())
fb1 = InteriorFacetBasis(mesh, ElementTriP1())  # same basis, different quadrature

g = cb1.project(f)
hc = cb1.interpolate(g)
hf = fb1.interpolate(g)  # also works, because g is equally valid in cb1 and fb1.

When working with lots of basis, and moving things from P1 and P0 into the same Functional for integration, I'm going to need to use cell + facet basises just like this, one of each, to get everything into a DiscreteField. Care has to be taken to make sure the DiscreteFields are using the same quadrature, even if they come from different basis classes.

Forgive me asking questions like this over and over, I keep getting index errors, size mismatch errors, etc, and they all keep coming back to messing up/ mixing up the basis and/or quadrature points.

But is my code/ideas in this comment correct?

[kinnala]

Yes everything's correct.