FacetBasis questions

[gatling-nrl]

Suppose I want to integrate

along a line from (0.5, 0) to (0.5, 1). So I try

mesh = skfem.MeshTri().refined(1)
basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
vert_line_facets = mesh.facets_satisfying(lambda x: np.isclose(x[0], 0.5))
vert_line_basis_p1 = skfem.FacetBasis(mesh, basis_p1.elem, facets=vert_line_facets)

def f(x):
    return 2*abs(x[0]+x[1]-1)

projection_p1 = skfem.utils.projection(f, basis_p1)

@skfem.Functional
def test1(w):
    return w['u']

skfem.asm(test1, vert_line_basis_p1 , u=basis_p1.interpolate(projection_p1))

ValueError: Quadrature mismatch: 'u' should have same number of integration points as the basis object.

I think this makes sense, because how could basis_p1 and vert_line_basis_p1 have the same quadrature points... they aren't even the same dimensionality. But still there must be some way to get projection_p1 onto the quadrature points of vert_line_basis_p1.

So I just tried switching out the basis used to do the interpolating:

skfem.asm(test1, fbasis_p1, u=vert_line_basis_p1.interpolate(projection_p1))

0.4999999999999999

But how is vert_line_basis_p1 able to understand a projection created with basis_p1? Is it just some lucky coincidence in this toy problem? What is the right way to interpolate projection_p1 for use in the @Functional test1?

[kinnala]

Actually, basis_p1.zeros() and vert_line_basis_p1.zeros() are compatible by design. The number and the ordering of DOFs is matching. So are the matrices if assembled using these bases. This is useful, e.g., when you have boundary forms in the problem. You are able to sum the resulting matrices. Quadrature points are different but it does not matter at the DOFs level.

[kinnala]

In fact, basis_p1 and vert_line_basis_p1 represent the same basis functions. Only integration happens over a different domain because the quadrature points are different.

[gatling-nrl]

Ok... so going from P1 to P0 I'm changing basis functions, and I use .with_element() to keep the quadrature points the same so I can work with projections from both basis sets in the same Forms.

Going from CellBasis_P1 to FacetBasis_P1, I am not changing the basis set at all, but am on purpose changing the quadrature points so I can integrate over a different domain.

Is it wasteful to build FacetBasis from scratch if I already have a CellBasis of the same "kind"?

[gatling-nrl]

is CellBasis(...,quadrature=points_along_the_facets) equivalent to FacetBasis(...)?

[kinnala]

Not really because for both the most expensive part is evaluating the global basis at a large number of points. Now these points are (in general) completely different so there isn't much that can be shared. With some special quadrature rule you might have the boundary quadrature to be a subset of domain quadrature but I didn't consider that at all.

[kinnala]

is CellBasis(...,quadrature=points_along_the_facets) equivalent to FacetBasis(...)?

It is not because FacetBasis uses quadrature points on the boundary reference domain (line) and CellBasis uses quadrature points on the triangle. Moreover, quadrature kwarg is used to specify the quadrature points in the reference triangle/line. You could technically take any quadrature rule from https://github.com/nschloe/quadpy and pass it to quadrature if you make sure that the reference domain is the same. Usually it's some simple scaling to change from one reference domain to another.

[gatling-nrl]

mesh = skfem.MeshTri().refined(1)
def is_my_facet(x):
    return x[0]>-1
#     return np.logical_and(
#         np.isclose(x[0], 0.75),
#         np.isclose(x[1], 0.75),
#     )
    
basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
my_facet = mesh.facets_satisfying(is_my_facet)
fbasis_p1 = skfem.FacetBasis(mesh, basis_p1.elem, facets=my_facet)

@skfem.Functional
def test1(w):
    pts = w.x.reshape(2,-1)
    for p, n in zip(pts.T, 0.2*w.n.reshape(2,-1).T):
        plt.arrow(*p, *n, width=.0, length_includes_head=True, head_width=.025, head_length=.05, facecolor='k', fill=True)
    plt.plot(pts[0], pts[1], 'or')
    return w.x[0]


fig, ax = plt.subplots(figsize=(5,5))
skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
ax.set_xlabel('x[0]'); ax.set_ylabel('x[1]')
ax.set_title('test1, w.x')
skfem.asm(test1, fbasis_p1)

I'm not sure how to make sense of w.n.

[kinnala]

There is no other rule than its outward on the boundary. Inside it has arbitrary direction.

[gatling-nrl]

Why not just

super().__init__(...)

instead of

scikit-fem/skfem/assembly/basis/boundary_facet_basis.py

Line 59 in d5aba7e

super(BoundaryFacetBasis, self).__init__(mesh,

and

scikit-fem/skfem/assembly/basis/interior_facet_basis.py

Line 40 in d5aba7e

super(InteriorFacetBasis, self).__init__(mesh,

Asking because I want to make sure I'm not missing something subtle about SubdomainFacetBasis.

[kinnala]

I think that's some remnant of Python 2 era.