FacetBasis constructed around a closed loop doesn't have consistent normal

[gatling-nrl]

mesh = skfem.MeshTri().refined(1)
def is_facet_on_loop(x):
    return (abs(x[0]-.3)+abs(x[1]-.3)) < .3
    
loop_facets = mesh.facets_satisfying(is_facet_on_loop)
fbasis_p1 = skfem.FacetBasis(mesh, skfem.ElementTriP1(), facets=loop_facets)

@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]  # arbitrarily, to suppress exception


fig, ax = plt.subplots(figsize=(5,5))
skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
skfem.asm(test1, fbasis_p1)  # to add quadrature points and arrows to plot
ax.set_xlabel('x[0]'); ax.set_ylabel('x[1]')

I'm going to guess this is a very non-trivial problem to solve, but it would be awesome if FacetBasis automatically made the normal go consistently in or out when it had a closed loop of facets... perhaps FacetBasis(..., outward=True) and outward could default to None preserving current behavior. But if outward was True or False, FacetBasis would try to comply and raise if the supplied facets didn't actually form a loop.

[gatling-nrl]

So one brute force solution, after you know you have a closed loop, would be to https://en.wikipedia.org/wiki/Point_in_polygon test the arrowheads and flip the failing ones. I assume this extends to 3D...

Ugh... there are a lot of edge cases and degenerate loops. Does a figure 8 have an outward normal? "It depends." seems to be the answer :( and figuring out whether a set of facets close a loop doesn't look trivial to me at all.

So I guess there's a lot of pathological geometries and that's annoying, but I think it would be really helpful even if it started out only working for convex shapes.

[gatling-nrl]

Ok, for convex loops, it is possible to write

mesh = skfem.MeshTri().refined(1)
def is_facet_on_loop(x):
    return (abs(x[0]-.3)+abs(x[1]-.3)) < .3
    
loop_facets = mesh.facets_satisfying(is_facet_on_loop)
fbasis_p1 = skfem.FacetBasis(mesh, skfem.ElementTriP1(), facets=loop_facets)

@skfem.Functional
def test1(w):
    interior_pt = (.3, .3)
    pts = w.x.reshape(2,-1)
    for p, n in zip(pts.T, 0.2*w.n.reshape(2,-1).T):
        # test whether w.n was pointing in or out for this point
        # ONLY WORKS FOR CONVEX LOOPS
        pos = np.linalg.norm(interior_pt-(p+n))
        neg = np.linalg.norm(interior_pt-(p-n))
        if pos < neg:
            # w.n was pointing in for this p, so flip it
            n *= -1
        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]  # arbitrarily, to suppress exception


fig, ax = plt.subplots(figsize=(5,5))
skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
skfem.asm(test1, fbasis_p1)  # to add quadrature points and arrows to plot
ax.set_xlabel('x[0]'); ax.set_ylabel('x[1]')

[gdmcbain]

I did have a vague idea on something like this the other day when I was thinking about electrodes and subdomains. I didn't get to try it out yet, but what if we had a characteristic function for the interior, i.e. one inside, zero outside. Then the sign of the component of the gradient in the direction of the normal should tell which way the normal points?

[gatling-nrl]

Hmmm... so give FacetBasis a P0 type projection and insist that it has normal vectors in the direction of positive gradient. Dude, that's an amazing idea. Because it could raise if the gradient was zero. In that event, one of the facets wasn't on the boundary between "in" and "out". That's incredibly flexible too, handling lots of the pathological cases. You could have disjoint regions and your idea will still produce "outward" normals... at least outwards in the sense that I think of the words. Math papers on geometry are hard to read.

I love your idea.

[gatling-nrl]

If I can dream up some syntax.... this seems neat:

interior = CellBasis(mesh, ElementTriP0()).project(is_interior)
arbitrary_boundary = FacetBasis(mesh, ElementTriP1(), interior_cells=interior)

Let FacetBasis do the work of finding the facets that bound interior_cells.

But why have two lines when you can have one?

arbitrary_boundary = FacetBasis(mesh, ElementTriP1(), interior_cells=is_interior)

I don't like the name "interior_cells" very much. But this feels really extremely natrual to write.

[gatling-nrl]

I guess it's not going to be quite as simple as I thought.... the gradient is all zeros.

mesh = skfem.MeshTri().refined(1)
basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
interior = basis_p0.zeros()
interior[basis_p0.get_dofs(elements=lambda x: abs(x[0]-.3)+abs(x[1]-.3)<.1)] = 1
fig, ax = plt.subplots(figsize=(6,6))
skfem.visuals.matplotlib.plot(basis_p0, interior, ax=ax, shading='gouraud')
skfem.visuals.matplotlib.draw(mesh, ax=ax)
print(basis_p0.interpolate(interior).grad)

[[[0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]]

 [[0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]]]

[gdmcbain]

Ah, the characteristic function might need to be continuous, like ElementTriP1.

[kinnala]

I wouldn't advice using BoundaryFacetBasis to do integration on internal facets because you can access only one of the traces. There is another class InteriorFacetBasis for that purpose. And I think the suggestion of @gdmcbain is a good one.

[gatling-nrl]

I didn't even realize till just now FacetBasis==BoundaryFacetBasis. And the distinction between InteriorFacetBasis and BoundaryFacetBasis is quite lost on me... This must be an implementation based distinction? I don't understand the significance of making those two classes, although I would bet it makes implementation of solutions to certain specific problems easier?

Are my labels right in the picture below?

[gatling-nrl]

Mathematically we have the concept of a "trace". The value of the trace depends on the side which is used to approach the discontinuity. On external boundary the trace is unique because nothing happens outside of the domain. On internal facets you choose the direction of the trace as follows: InteriorFacetBasis(..., side=0) or side=1.

[gatling-nrl]

So the use case I am focused on is this: I want to draw the blue line, get an outward facing normal in w.n and integrate dot(interpolation(thing).grad, w.n) over the blue line.

Frustrations I am encountering so far are that gmsh makes labeling the blue line really extremely aggravating and skfem doesn't want to give me an outward normal for it even if I got the labeling to work or resort to np.logical_whatever().

[gatling-nrl]

From an outside user perspective, and especially not knowing the inner workings of implementing FEM all that well, the distinction between BoundaryFacetBasis and InteriorFacetBasis is very, very confusing. My three drawings above all look fundamentally the same to me.

It is even more so something I don't understand because they apparently are all the same mathematical "basis" (and the same "basis" as something that seems higher dimensional, CellBasis) and apparently only differ in capabilities of interrogating the space. Again, looking from the outside in, and judging from the names, I thought they were different mathematical basis sets. Now I suspect the evolution of their existence has more to do with implementation and solving specific domain problems... which is why understanding the implementation is a prerequisite to understanding how to use the various basis objects.

This is rather frustrating for an outside user (me). But I think also understandable, because I also suspect there are significant implications to performance when making these kinds of design choices. Too many times I find good computer science abstractions are just in direct opposition of computations efficiency, and I'm not at all surprised if that is the case here. Am I wrong about this?

Sorry... its late for me and I just want to understand how to identify and integrate around the blue line, dotted with the outward normal. (I intend to make everything inside the blue line a dirichlet boundary condition too, if that matters.)

[gatling-nrl]

how to identify ... the blue line

This also has been a huge source of aggravation. Neither gmsh nor skfem make this a simple task, in the general case. For my present efforts, the blue line will be the closest approximation to some circle. In gmsh I can easily define the subdomain that is the disk enclosed by the circle. But getting the perimeter of the circle, at least with occ has been impossible so far. And skfem doesn't have a facets=the_ones_closest_to_this_path way of interrogation the mesh.

[gdmcbain]

Ah, the characteristic function might need to be continuous, like ElementTriP1.

But seeing as we might not want that for the purposes of defining a discontinous conductivity as in ex17 (and while we're at it just use the conductivity instead of a zero–one characteristic function), what about a Raviart–Thomas gradient of the P0 conductivity/characteristic function?

[kinnala]

From an outside user perspective, and especially not knowing the inner workings of implementing FEM all that well, the distinction between BoundaryFacetBasis and InteriorFacetBasis is very, very confusing. My three drawings above all look fundamentally the same to me.

InteriorFacetBasis was introduced for the specific purpose of evaluating a posteriori error estimators in adaptive mesh refinement. Later it was extended to support interior penalty/discontinuous Galerkin methods. The whole point is that in error estimators the value you are integrating is different from one element to another and you want to minimize that by adaptive refinement.

[kinnala]

Sorry... its late for me and I just want to understand how to identify and integrate around the blue line, dotted with the outward normal. (I intend to make everything inside the blue line a dirichlet boundary condition too, if that matters.)

I can give you an example soon which uses the idea of @gdmcbain using the indicator function.

[kinnala]

Are my labels right in the picture below?

Yes, but the main difference is missing: for interior facets (not on the boundary) you have in general two values, one from each of the elements. It is just coincidence that what you are integrating right now has unique value.

[kinnala]

So this is the indicator function that @gdmcbain is proposing:

The nice thing is that -grad(ind) will give the direction for outward normal.

Edit: Updated the figure to match the blue line above.

[kinnala]

Are my labels right in the picture below?

... It is just coincidence that what you are integrating right now has unique value.

I was incorrect. What you are integrating does not have a unique value because gradient has different value inside and outside of the blue line. So which of these you want to use? How about on the boundary? You are making a mistake by thinking that this is a simple matter. :-)

[kinnala]

Took me a while to resolve because there was a bug in COOData.todefault() (apparently accidentally pushed the bug fix directly to master - I hope it is correct - I need to sleep a bit more).

However, after fixing the bug, the idea of @gdmcbain can be implemented like this:

import numpy as np
from skfem import *
from skfem.helpers import dot, grad
from skfem.visuals.matplotlib import *

m = MeshTri().refined(1)
basis_p1 = Basis(m, ElementTriP1())

# create indicator function for the domain enclosed by the blue line
ind = basis_p1.zeros()
ind[basis_p1.get_dofs(elements=lambda x: abs(x[0]) + abs(x[1]) < 1.)] = 1
ind[m.nodes_satisfying(lambda x: (x[0] == 0) * (x[1] == 0))] = 0

# integrating over a combination of boundary and interior facets,
# you need all facets and all traces you might ever imagine:
bases = [
    InteriorFacetBasis(m, ElementTriP1(), side=0),
    InteriorFacetBasis(m, ElementTriP1(), side=1),
    FacetBasis(m, ElementTriP1()),
]

@Functional
def trace(w):
    # outward normal, normalize
    n = -grad(w['ind'])
    length = np.sqrt(dot(n, n))
    # divide by length
    # special case: length can be zero because ind is sometimes constant
    n = np.divide(n, length, out=np.zeros_like(n), where=~np.isclose(length, 0.))
    # use ind to check that we are in the correct part of the domain
    return np.isclose(w['ind'], 1.) * dot(grad(w['thing']), n)

# use ind as "thing" for testing, integral should be like this:
# - 3 * 1 / (sqrt(0.5) / 2) * sqrt(0.5)
print(asm(trace, bases, ind=ind, thing=ind))

fig, ax = plt.subplots(figsize=(6,6))
plot(basis_p1, ind, ax=ax, shading='gouraud')
draw(m, ax=ax)

show()

This now uses the gradient just outside of the boundary on interior facets.

[kinnala]

arbitrary_boundary = FacetBasis(mesh, ElementTriP1(), interior_cells=is_interior)

So you would like to initialize in a single call a basis object which allows integrating both, over boundary facets and over interior facets, as defined by is_interior. Moreover, on interior facets it would use (by default) the values of the function just outside of the boundary.

The call would be something like:

sfbasis = SubdomainFacetBasis(mesh, ElementTriP1(), elements=is_interior)

Definitely possible to do such a thing.

(There is a caveat: This might give surprising results if the finite element is discontinuous because on interior facets it uses values "outside" of the subdomain and on boundary facets it uses values "inside" of the subdomain.)

This functionality would however need to be inside a new class because it does not fit into the BoundaryFacetBasis/InteriorFacetBasis distinction too well. It doesn't matter though, we already have MortarFacetBasis which is neither of these, we can have more *FacetBasis classes if there is a proper need for them.

[gdmcbain]

Thanks very much @kinnala for implementing my sketchy idea of a ℙ₁ approximation of the indicator function. Meanwhile I've followed up on my revision of that to try harder to get a meaningful weak gradient of the true discontinuous indicator by using ElementTriRT0 but got stuck trying to pass a scalar coefficient to a vector form #823.

[gdmcbain]

I sorted that out easily enough and calculated what I think should be the RT0 gradient of the P0 indicator but am now pondering how to render it #824 to see whether it makes any sense.

[gatling-nrl]

which allows integrating both, over boundary facets and over interior facets

The more I think about this, the more I have a hard time justifying anything but interior facets. I have not thought of any use cases for the blue boundary above including part of the domain covered by the boundary conditions.

Edit:

(There is a caveat: This might give surprising results if the finite element is discontinuous because on interior facets it uses values "outside" of the subdomain and on boundary facets it uses values "inside" of the subdomain.)

With this in mind, I especially wonder if only interior facets should ever be in the blue boundary.

[kinnala]

I included some internal changes in #822 which allows more easily implementing SubdomainFacetBasis (or similar). There is now an undocumented kwarg _tind in BoundaryFacetBasis which would eventually be used by SubdomainFacetBasis to pass a list of elements corresponding to the list of facets on the blue boundary. Then SubdomainFacetBasis constructor would need to find the facets on the blue boundary (passed to super() as facets=...) and next the elements just outside/inside of the boundary (passed to super() as _tind=...).

[gatling-nrl]

SubdomainFacetBasis should inherit BoundaryFacetBasis or AbstractBasis?

[kinnala]

It could be pretty similar to the implementation of InteriorFacetBasis (inherit from BoundaryFacetBasis):

class SubdomainFacetBasis(BoundaryFacetBasis):

    def __init__(self,
                 mesh: Mesh,
                 elem: Element,
                 mapping: Optional[Mapping] = None,
                 intorder: Optional[int] = None,
                 quadrature: Optional[Tuple[ndarray, ndarray]] = None,
                 elements = None,
                 side: int = 0,  # 0 outside, 1 inside?
                 dofs: Optional[Dofs] = None):

        if side not in (0, 1):
            raise ValueError("'side' must be 0 or 1.")

        if elements is None:
            raise ValueError("Subdomain must be specified.")

        elements = self._normalize_elements(elements)
        facets = ...  # find facets on the outer boundary of 'elements'
        tind = ...  # find element indices just outside/inside of 'facets' depending on side

        super(InteriorFacetBasis, self).__init__(mesh,
                                                 elem,
                                                 mapping=mapping,
                                                 intorder=intorder,
                                                 quadrature=quadrature,
                                                 facets=facets,
                                                 _tind=tind,
                                                 dofs=dofs)

[gatling-nrl]

elements just outside/inside of the boundary

Do you think it is a good idea to do just outside/inside depending on boundary or not, or should it refuse to include boundary facets?

[kinnala]

_tind can be any combination of outside/inside elements, but is there any use case? Inside elements would work for any facets but outside elements obviously do not exist for boundary facets. If you don't need a combination of inside/outside elements, then we can simply raise error if an element outside of the boundary facet is requested. If you feel that it would be useful in your use case, then we can fall back to using inside elements on boundaries.

[gatling-nrl]

Outward normals using only an indicator function projected into ElementTriP0!

import numpy as np
import matplotlib.pyplot as plt
import skfem
import skfem.visuals.matplotlib

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

basis_p0 = skfem.CellBasis(mesh, skfem.ElementTriP0())
facet_basis_p0_s0 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP0(), side=0)
facet_basis_p0_s1 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP0(), side=1)

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


def plot_normals(ax):
    arrow_style = dict(width=.0, length_includes_head=True, head_width=.025, head_length=.025, facecolor='k', fill=True)
    arrow_length = .05
    def _form(w):
        # arrays are shaped to group by facet, but we don't care for that here so reshape(2,-1)
        ii0_s0 = w.ind_p0_s0.reshape(-1)
        ii0_s1 = w.ind_p0_s1.reshape(-1)        
        pp = w.x.reshape(2,-1).T  
        nn = w.n.reshape(2,-1).T
        for p, n, i0_s0, i0_s1  in zip(pp, nn, ii0_s0, ii0_s1):
            if i0_s0 > 0:
                ax.arrow(*p, *(arrow_length*n), **arrow_style)
            elif  i0_s1 > 0:
                ax.arrow(*p, *(-arrow_length*n), **arrow_style)  # NOTE THE MINUS SIGN!!!
        return w.x[0] # must return *something* of the right shape
    return skfem.Functional(_form)


fig, ax = plt.subplots(1, 1, figsize=(6,6))
ax.set_title('ind_p0 with outward normals')
skfem.visuals.matplotlib.plot(basis_p0, ind_p0, ax=ax, cmap='Set3', shading='gouraud')
skfem.visuals.matplotlib.draw(mesh, ax=ax)
skfem.asm(
    plot_normals(ax),
    [facet_basis_p0_s0, facet_basis_p0_s1],
    ind_p0_s0=facet_basis_p0_s0.interpolate(ind_p0),
    ind_p0_s1=facet_basis_p0_s1.interpolate(ind_p0),
)

Yeah, yeah, I need to get rid of that for loop and wrap this up in an AbstractBasis class. Oh, and this is using #849, so that's good.

By the way, you said

I don't think this is really an issue. I don't see when you would want to iterate like that. It just never makes any sense.

and I understand where you were coming from... but doing stuff like this is very instructive and I'm glad that I can.

[gatling-nrl]

Improved logic, no more "if" statement.

import numpy as np
import matplotlib.pyplot as plt
import skfem
import skfem.visuals.matplotlib

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

basis_p0 = skfem.CellBasis(mesh, skfem.ElementTriP0())
facet_basis_p0_s0 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP0(), side=0)
facet_basis_p0_s1 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP0(), side=1)

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


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


fig, ax = plt.subplots(1, 1, figsize=(6,6))
ax.set_title('ind_p0 with outward normals')
skfem.visuals.matplotlib.plot(basis_p0, ind_p0, ax=ax, cmap='Set3', shading='gouraud')
skfem.visuals.matplotlib.draw(mesh, ax=ax)
skfem.asm(
    plot_normals(ax),
    [facet_basis_p0_s0, facet_basis_p0_s1],
    ind_p0_s0=facet_basis_p0_s0.interpolate(ind_p0),
    ind_p0_s1=facet_basis_p0_s1.interpolate(ind_p0),
)

[gatling-nrl]

Using #851

import numpy as np
import matplotlib.pyplot as plt
import skfem
import skfem.visuals.matplotlib
from skfem.assembly.basis.subdomain_facet_basis import SubdomainFacetBasis

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

def elem_query(x):
    return np.logical_or(
        np.linalg.norm(x.T-[0.35,0.35], axis=1) < .05,
        np.linalg.norm(x.T-[0.65,0.65], axis=1) < .1,
    )

### Just for visualization, not required by SubdomainFacetBasis
cell_basis_p0 = skfem.CellBasis(mesh, skfem.ElementTriP0())
ind_p0 = cell_basis_p0.zeros()
ind_p0[cell_basis_p0.get_dofs(elements=elem_query)] = 1
###

basis = SubdomainFacetBasis(mesh, skfem.ElementTriP1(), elements=elem_query, side=0)

def plot_normals(ax):
    arrow_style = dict(width=.0, length_includes_head=True, head_width=.025, head_length=.025, facecolor='k', fill=True)
    arrow_length = .05
    def _form(w):
        # arrays are shaped to group by facet, but we don't care for that here so reshape(2,-1)
        pp = w.x.reshape(2,-1).T  
        nn = w.n.reshape(2,-1).T
        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(1, 1, figsize=(6,6))
ax.set_title('ind_p0 with outward normals')
skfem.visuals.matplotlib.plot(cell_basis_p0, ind_p0, ax=ax, cmap='Set3', shading='gouraud')
skfem.visuals.matplotlib.draw(mesh, ax=ax)

skfem.asm(
    plot_normals(ax),
    basis,
)

[gdmcbain]

I see, yes, just using a difference of the indicator functions defined on the two sides is much simpler than my suggestion of the Raviart–Thomas gradient! Nice. Maybe one would have been enough, seeing whether the indicator is one or zero? But I see that this is subsequently subsumed by #851, which I'll look into next.