export trace

[gdmcbain]

This is a vague feature request; vague because I'm still thinking about how best to formulate the problem, but let's see if I can already make sense.

Say in ex02, the square plate subject to normal pressure, clamped on the left and pinned right and top, that I wanted to export the deflexion x along the {'bottom': m.facets_satisfying(lambda x: x[1] == 0)}. Say I wanted to couple the scikit-fem finite element model to something external. Say, for a definite example, I had a soap-film (Δ w = 0) pinned to the bottom of the plate x[1] == 0 and pinned (w = 0) around the adjacent square 0 < x[0] < 1, 1 < x[1] < 0 with boundary condition on x[1] == 0 prescribed by the deflexion of the plate. That's one-way coupling, but also in other cases two-way coupling might be of interest, and we might use a Lagrange multiplier as suggested for coupling subdomains #163.

I can extract the dofs with

dofs['bottom'] = ib.find_dofs({'bottom': m.facets_satisfying(lambda x: x[1] == 0)})['bottom']
np.vstack([ib.doflocs[0, dofs['bottom'].all()], x[dofs['bottom'].all()]])

or the ordinates

ib.interpolator(x)(ib.doflocs[:, dofs['bottom'].all()])

(Here the dofs aren't the same as the ordinates: the last eight of dofs['bottom'].all() are the dofs['bottom'].facet['u_n'] since x is built with ElementTriMorley.)

But how is the connectivity of dofs['bottom'] to be extracted and exported?

Does it make sense to export some kind of facet mesh corresponding to the FacetBasis?

bb = FacetBasis(ib.mesh, ib.elem, facets=m.facets_satisfying(lambda x: x[1] == 0))

on which the dofs from the ib are also dofs? (I'm not sure what this means when isinstance(ib.elem, ElementTriMorley).)

Or if the external model doesn't have Morley elements, does it make sense to project the trace of x along 'bottom' onto some other facet finite element space, e.g. facets with domain element ElementTriP2 or ElementTriP1 or ElementTriP0? Or onto a lower (dim = ib.elem.dim - 1) dimensional finite element space defined on the part of the boundary; e.g. ElementLineP0 #508, ElementLineP1, or ElementLineP2? If the deflexion of the plate were projected onto ElementTriP2 is there a sense in which the FacetBasis is like ElementLineP2?

So obviously it can all be done by interpolation (ib.interpolator) and reinterpolation on the other side, but are there more direct ways, e.g. involving projection?

Some related issues: #150, #158, #163, #261, #271.

[gdmcbain]

The following extracts something like a MeshLine with ElementLineP1 (actually the corresponding meshio.Mesh).

from skfem import *
import numpy as np
import meshio
from ex02 import m, ib, dofs, x

bottom_facets = m.facets_satisfying(lambda x: x[1] == 0)
dofs['bottom'] = ib.find_dofs({'bottom': bottom_facets})['bottom']
bb = FacetBasis(ib.mesh, ib.elem, facets=bottom_facets)

points = bb.doflocs[0, dofs['bottom'].nodal['u']]
lines_global = bb.element_dofs[
    np.isin(bb.element_dofs,
    dofs['bottom'].nodal['u']).reshape((-1, bb.nelems))
]
lines_local = np.unique(lines_global, return_inverse=True)[1].reshape((-1, bb.nelems))

mesh = meshio.Mesh(np.outer(points, [1, 0, 0]), 
                   [('line', lines_local.T)],
                   point_data={'deflexion': x[dofs['bottom'].nodal['u']]})
mesh.write('bottom.xdmf')

It works, but I don't like it so much as it just plucks the .nodal['u'] dofs from the triangular Morley space rather than projecting onto the lineal P1.

[gdmcbain]

Projecting from ElementLineP1 to ElementLineP0 is easy.

x1 = x[dofs['bottom'].nodal['u']]
mesh = MeshLine(points, lines_local)
l1b = InteriorBasis(mesh, ElementLineP1())
l0b = InteriorBasis(mesh, ElementLineP0())
x0 = project(x1, l1b, l0b)

But I don't see how to project from n dimensions to n − 1.

[kinnala]

I started looking into it. I first extended MeshTri with the method

    def trace(self):
        from ..mesh_line import MeshLine
        ix = self.boundary_facets()
        ix = self.facets[:, ix]
        ixuniq, a = np.unique(ix, return_index=True)
        b = np.zeros(np.max(ix) + 1, dtype=np.int64)
        b[ixuniq] = np.arange(len(ixuniq), dtype=np.int64)
        return MeshLine(self.p[:, ixuniq], b[ix], validate=False)

Then I was able to do

In [5]: me = m.trace()                                                          

In [6]: me                                                                      
Out[6]: One-dimensional mesh with 32 vertices and 32 elements.

In [7]: from skfem.visuals.matplotlib import *                                  

In [8]: plot(me, me.p[0])                                                       
Out[8]: <matplotlib.axes._subplots.AxesSubplot at 0x7f6b3135b970>

In [9]: show() 

This is as far as I got. Now have to do something else for a while.

[kinnala]

Probably better to have this code in Basis as you say.

[gdmcbain]

Whoa. Can a MeshLine represent a planar curve?? Now that's interesting !

I suspect a Mesh.trace would also be of.interest; e.g. for quicker visualization of a Mesh3D. (If a Mesh2D can be a surface in space rather than just in the plane.)

[kinnala]

It kind of can but the resulting mesh does not (yet?) work with many of the other parts of the library. I'm thinking about how to best solve it.

[gdmcbain]

Right, yes, of course. Defining forms on curved manifolds is a big leap. Very interesting.

[kinnala]

I think right now I'd do something like (building on your earlier example)

from skfem import *
import numpy as np
import meshio
from ex02 import m, ib, dofs, x

bottom_facets = m.facets_satisfying(lambda x: x[1] == 0)
dofs['bottom'] = ib.find_dofs({'bottom': bottom_facets})['bottom']
bb = FacetBasis(ib.mesh, ib.elem, facets=bottom_facets, intorder=3)
bb2 = FacetBasis(ib.mesh, ElementTriP1(), facets=bottom_facets, intorder=3)
project(x, bb, bb2, I=bb2.get_dofs(lambda x: x[1] == 0.0).all())

I'll need to check how to build a mesh from the result.

[kinnala]

I see what you mean; matching between the dofs of bb and l1b is kinda nonobvious.

[kinnala]

In case of linears its simple:

from skfem import *
import numpy as np
import meshio
from docs.examples.ex02 import m, ib, dofs, x

bottom_facets = m.facets_satisfying(lambda x: x[1] == 0)
dofs['bottom'] = ib.find_dofs({'bottom': bottom_facets})['bottom']
bb = FacetBasis(ib.mesh, ib.elem, facets=bottom_facets, intorder=3)

# project to linears
bb2 = FacetBasis(ib.mesh, ElementTriP1(), facets=bottom_facets, intorder=3)
I = bb2.get_dofs(lambda x: x[1] == 0.0).all()
y = project(x, bb, bb2, I=I)

# build 1d mesh
ix = np.argsort(bb2.doflocs[0, I])
m = MeshLine(bb2.doflocs[0, I][ix])
basis = InteriorBasis(m, ElementLineP1())
from skfem.visuals.matplotlib import plot
plot(basis, y[ix])

The resulting figure is

However, generalizing this will require some work.

[kinnala]

In general the 1D mesh should be built using bottom_facets, their indexing and coordinates.

[kinnala]

Ugh, I already had it done above but just forgot about it.

[kinnala]

Here is a version utilizing bottom_facets:

from skfem import *
import numpy as np
import meshio
from docs.examples.ex02 import m, ib, dofs, x

bottom_facets = m.facets_satisfying(lambda x: x[1] == 0)
dofs['bottom'] = ib.find_dofs({'bottom': bottom_facets})['bottom']
bb = FacetBasis(ib.mesh, ib.elem, facets=bottom_facets, intorder=3)

# project to linears
bb2 = FacetBasis(ib.mesh, ElementTriP1(), facets=bottom_facets, intorder=3)
I = bb2.get_dofs(lambda x: x[1] == 0.0).all()
y = project(x, bb, bb2, I=I)

# build 1D mesh
ix = ib.mesh.facets[:, bottom_facets]
ixuniq, a = np.unique(ix, return_index=True)
b = np.zeros(np.max(ix) + 1, dtype=np.int64)
b[ixuniq] = np.arange(len(ixuniq), dtype=np.int64)
m = MeshLine(ib.mesh.p[0, ixuniq], b[ix])

# plot
basis = InteriorBasis(m, ElementLineP1())
from skfem.visuals.matplotlib import plot
plot(basis, y)

[kinnala]

Making it more generic:

from skfem import *
import numpy as np
import meshio
from docs.examples.ex02 import m, ib, dofs, x

bottom_facets = m.facets_satisfying(lambda x: x[1] == 0)
dofs['bottom'] = ib.find_dofs({'bottom': bottom_facets})['bottom']
bb = FacetBasis(ib.mesh, ib.elem, facets=bottom_facets, intorder=3)

def trace(basis, facets):
    # project to linears
    fbasis_from = FacetBasis(basis.mesh, basis.elem, facets=facets, intorder=3)
    fbasis_to = FacetBasis(basis.mesh, ElementTriP1(), facets=facets, intorder=3)
    I = fbasis_to.get_dofs(facets).all()
    y = project(x, fbasis_from, fbasis_to, I=I)

    # build 1D mesh
    ix = basis.mesh.facets[:, facets]
    ixuniq, a = np.unique(ix, return_index=True)
    b = np.zeros(np.max(ix) + 1, dtype=np.int64)
    b[ixuniq] = np.arange(len(ixuniq), dtype=np.int64)
    return MeshLine(basis.mesh.p[0, ixuniq], b[ix]), y

# plot
m, y = trace(ib, bottom_facets)
basis = InteriorBasis(m, ElementLineP1())
from skfem.visuals.matplotlib import plot
plot(basis, y)

[kinnala]

Work on this in this PR: #530