Using Project, Interpolate, and Functional

[gatling-nrl]

I'm still using Functional in unintended ways to better understand how to fit the data structures in skfem together. In this scenario I have a function that I want to project into TriP0 and TriP1. I also want these projections on CellBasis and on InteriorFacetBasis.

mesh = skfem.MeshTri().refined(1)

basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
facet_basis_p1_s0 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP1(), side=0)
facet_basis_p0_s0 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP0(), side=0)
facet_basis_p1_s1 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP1(), side=1)
facet_basis_p0_s1 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP0(), side=1)

ind_p0 = basis_p0.zeros()  # returns a "Projection" in P0 space
print('ind_p0.shape =', ind_p0.shape)
ind_p1 = basis_p1.zeros() # returns a "Projection" in P1 space
print('ind_p1.shape =', ind_p1.shape)


def elem_query(x):
    return np.linalg.norm(x.T-[0.4,0.4], axis=1) < .2
ind_p0[basis_p0.get_dofs(elements=elem_query)] = 1
ind_p1[basis_p1.get_dofs(elements=elem_query)] = 1


ind_p0_interpolated_by_basis_p0 = basis_p0.interpolate(ind_p0)
print('ind_p0_interpolated_by_basis_p0.value.shape =', ind_p0_interpolated_by_basis_p0.value.shape)

ind_p1_interpolated_by_basis_p1 = basis_p1.interpolate(ind_p1)
print('ind_p1_interpolated_by_basis_p1.value.shape =', ind_p1_interpolated_by_basis_p1.value.shape)

ind_p0_interpolated_by_facet_basis_p0_s0 = facet_basis_p0_s0.interpolate(ind_p0)
print('ind_p0_interpolated_by_facet_basis_p0_s0.value.shape =', ind_p0_interpolated_by_facet_basis_p0_s0.value.shape)

ind_p1_interpolated_by_facet_basis_p1_s0 = facet_basis_p1_s0.interpolate(ind_p1)
print('ind_p1_interpolated_by_facet_basis_p1_s0.value.shape =', ind_p1_interpolated_by_facet_basis_p1_s0.value.shape)

ind_p0.shape = (8,)
ind_p1.shape = (9,)
ind_p0_interpolated_by_basis_p0.value.shape = (8, 3)
ind_p1_interpolated_by_basis_p1.value.shape = (8, 3)
ind_p0_interpolated_by_facet_basis_p0_s0.value.shape = (8, 2)
ind_p1_interpolated_by_facet_basis_p1_s0.value.shape = (8, 2)

def functional_testing(ax):
    arrow_style = dict(width=.0, length_includes_head=True, head_width=.025, head_length=.025, facecolor='k', fill=True)
    arrow_length = .1
    def _form(w):
        print('w.x.shape = ', w.x.shape)
        print('w.n.shape = ', w.n.shape)
        print('w.ind_p0.shape = ', w.ind_p0.shape)
        print('w.ind_p1.shape = ', w.ind_p1.shape)
        print('--')
        # arrays are shaped to group by facet, but we don't care for that here so reshape(2,-1)
        ii0 = w.ind_p0.reshape(2,-1).T
        ii1 = w.ind_p1.reshape(2,-1).T
        gg = grad(w['ind_p1']).reshape(2,-1).T
        pp = w.x.reshape(2,-1).T  
        nn = w.n.reshape(2,-1).T
        print('pp.shape=',pp.shape)
        print('nn.shape=',nn.shape)
        print('gg.shape=',gg.shape)
        print('ii0.shape=',ii0.shape)
        print('ii1.shape=',ii1.shape)
        for p, n, in zip(pp, nn):
            ax.arrow(*p, *(arrow_length*n), **arrow_style)
        return w.x[0] # must return *something* of the right shape
    return skfem.Functional(_form)

fig, ax = plt.subplots(figsize=(8,8))
skfem.visuals.matplotlib.draw(mesh, ax=ax)
skfem.asm(
    functional_testing(ax),
    facet_basis,
    ind_p1=ind_p1,
    ind_p0=facet_basis_p1_s0.project(facet_basis_p0_s0.interpolate(ind_p0)),
)

w.x.shape =  (2, 8, 2)
w.n.shape =  (2, 8, 2)
w.ind_p0.shape =  (8, 2)
w.ind_p1.shape =  (8, 2)
--
pp.shape= (16, 2)
nn.shape= (16, 2)
gg.shape= (16, 2)
ii0.shape= (8, 2)
ii1.shape= (8, 2)

Things are starting to work, but the syntax is really awkward. Particularly

asm(
    ...
    facet_basis_p1_s1,
    ...
    ind_p0=facet_basis_p1_s0.project(facet_basis_p0_s0.interpolate(ind_p0)),

ind_p0 wasn't created with facet_basis_p0_s0, so why am I using it to interpolate? I think the answer is that facet_basis_p0_s0 shares basis functions with basis_p0 and it shares quadrature points with facet_basis_p1_s0 and so it has to be the bridge between them.

This kind of makes sense in retrospect, but its very awkward. Is there is a better way I'm not thinking of?

Also, do confusing things start happening if you mix the side argument of the facet basis in the three places (interpolation, projection, and asm argument)?

[gatling-nrl]

Hmm, things have gone wrong. I should have realized ind_p0=facet_basis_p1_s0.project was going to mess with things.

def plot_annotations(ax):
    def _functional(w):
        pts = w.x.reshape(2,-1).T
        vals = w['ind_p0'].value.reshape(-1)
        for p, v in zip(pts, vals):
            ax.plot(p[0], p[1], 'or')
            ax.annotate(f'{v:.1f}', p, size=15)
        return w['ind_p1']
    return skfem.Functional(_functional)


fig, ax = plt.subplots(figsize=(8,8))
skfem.visuals.matplotlib.draw(mesh, ax=ax)
skfem.visuals.matplotlib.plot(basis_p0, ind_p0, ax=ax, cmap='Blues', vmax=1.5, shading='gouraud')
skfem.asm(
    plot_annotations(ax),
    facet_basis_p1_s1,
    ind_p1=ind_p1,
    ind_p0=facet_basis_p1_s0.project(facet_basis_p0_s0.interpolate(ind_p0)),
)

ind_p0 is supposed to be 0 or 1 everywhere. I understand that ind_p0=facet_basis_p1_s0.project mangled it, but what am I supposed to do instead?

[kinnala]

Isn't

asm(
    ...
    facet_basis_p1_s1,
    ...
    ind_p0=facet_basis_p0_s0.interpolate(ind_p0)),

enough?

[kinnala]

These keyword arguments to asm are in general of type DiscreteField but if they are not then basis.interpolate is called behind the scenes because quite often people work with only one basis.

[gatling-nrl]

Ah, okay... so if I give asm a keyword argument of type ndarray, it assumes that is a projection (GridFunction) and calls interpolate using the second argument to asm. If I give asm a keyword argument of type DiscreteField, it assumes this has values at the same quadrature points as the second argument to asm.

Is there any loss of information (other than floating point rounding errors) moving back and forth between GridFunction and DiscretreField?

basis_p0 = basis_p1.with_element(...)
x = basis_p0.project(some_callable)
x_df = basis_p0.interpolate(x)
y = basis_p1.project(x_df)
y_df = basis_p1.interpolate(y)
x2 = basis_p0.project(y_df)
np.isclose(x, x2) # True?
x == x2 # False?

If I were to go the other way though, from P1 down to P0 an then back up to P1, there would be loss of information, correct?

By the way, your suggestion worked perfectly.

skfem.asm(
    plot_annotations(ax),
    facet_basis_p1_s1,
    ind_p1=ind_p1,
    ind_p0=facet_basis_p0_s1.interpolate(ind_p0),
)

[kinnala]

For all x (with correct shape) it should hold that array_almost_equal(x, basis.project(basis.interpolate(x))). However, as you say, using different spaces there will be loss of information.

[gatling-nrl]

Thanks! #816 is making way more sense now.