SubdomainFacetBasis implementation, how do I...

[gatling-nrl]

SubdomainFacetBasis is supposed to give you at least a P1 (but ideally any) Basis on a set of facets comprising an internal boundary of the domain. (Facets on the exterior boundary aren't allowed.) The normals are supposed to be consistently outward. To do this I'm subclassing BoundaryFacetBasis and there are three important steps I'm needing help figuring out. (1) I have working, but its hacky. (2) is working and probably okay. (3) is just broken and I don't know how to fix it. Here they are:

1. How to get the facets on the boundary.

If the boundary is defined by a zero/one indicator function (1 in a cell , then facets with discontinuous traces are on the boundary and the rest aren't. I used asm to make this happen. As a by-product, the sign of the difference of these traces indicates the whether the skfem normal is "inward" or "outward", so you can just multiply this difference with the self.normals and they'll all get pointed the right way (outward). This method is robust, does not require contiguous subdomains and automatically excludes actual boundary facets.

However, it uses asm and Functional. This creates an import loop, so that's bad. It is also not really the purpose of asm. On the other hand, this works. So I feel like getting rid of asm and doing it right is going to involve rewriting mesh traversal code that already exists and that I don't know how to write. This is the main reason I wondered in a different discussion if a formal iterator is something we should have.

Nevertheless, this implementation works, in that it does found the subdomain boundary and which facets need to be "flipped".

        def get_boundary(facet_list, flip_list):
            def _form(w):
                for i, indicator in zip(
                    np.nonzero(mesh.f2t[1] != -1)[0], w.ind_p0_s0 - w.ind_p0_s1
                ):
                    if indicator[0] != 0 and i not in facet_list:
                        facet_list.append(i)
                        flip_list.append(indicator)
                return w.x[0]  # must return *something* of the right shape

            return Functional(_form)

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

        ind_p0 = basis_p0.zeros()
        ind_p0[basis_p0.get_dofs(elements=elements)] = 1

        facet_list = list()  # type: List[int]
        flip_list = list()  # type: List[ndarray]
        asm(
            get_boundary(facet_list, flip_list),
            [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),
        )
 
        flips = np.array([a[0] for a in flip_list])
        facets = np.array(facet_list)

2. Flipping the normals.

Since in the previous step I ended up with an array of integers that are 1 for already correct normals and -1 for normals pointing the wrong way, "fixing" the normals was as easy as

if side==1:
    flips*= -1

self.normals *= flip_list

Except this broke after #849. A work around here is to insist

self.normals = DiscreteField(flip_list * self.normals)

But I want to make sure this is in the spirit of the changes in #849. It seems like all the things in DiscreteField are the result of linear operations and at least scalar multiplication and addition should be supported because they can just be applied to everything.

Also I wonder if ufunc(DiscreteField, array) and ufunc(DiscreteField, DiscreteField) should return DiscreteField with the extra attribute stripped?

3. How to I get the right trace values?

Here is where the implementation is just broken. I know which trace needs to be used on which facet. I thought I could write

tind = np.hstack([
    mesh.f2t[0, facets[flips>0]],
    mesh.f2t[1, facets[flips<0]],
])

But it didn't seem to work, and I'm not too sure what is tind, f2t, or AffineMapping.... these things are getting way down in the guts of implementation. I think this is another reason for the formal iterator I mention earlier, but the biggest reason why I don't write that myself: I'm not even sure where to start. It might be that this idea of iteration is just me not understanding what some of the lower level functions do. There names are pretty terse. I'll try to add some docstrings in #851 too.

[kinnala]

Here are some ideas:

In [1]: from skfem import *

In [2]: m = MeshTri().refined(3)

In [3]: tind = m.elements_satisfying(lambda x: (x[0] - 0.5) ** 2 + (x[1] - 0.5)
   ...: **2 < 0.25**2)

In [4]: tind
Out[4]: 
array([  7,   9,  11,  18,  20,  22,  31,  39,  54,  62,  63,  71,  73,
        75,  82,  84,  86,  94, 103, 105, 107, 114, 116, 118, 126, 127])

In [6]: import numpy as np

In [7]: all_facets, counts = np.unique(m.t2f[:, tind], return_counts=True)

In [9]: facets = all_facets[counts == 1]  # facets on subdomain boundary?

In [10]: m.f2t[:, facets]  # triangles on both sides of subdomain boundary?
Out[10]: 
array([[ 82,  62,  84,  94,  89,  31,  91,  55,  50,  41,  52,  43, 126,
         95],
       [ 90,  70,  92,  78,  73,  47,  75,  63, 114, 105, 116, 107,  30,
        127]])

Finally you check if first or second row is in the given elements set.

[kinnala]

No need to modify self.normals, it should come automatically if you set _tind correctly.

[kinnala]

Here is something I wrote but did not try if it works:

from skfem import *
import numpy as np


class SubdomainFacetBasis(BoundaryFacetBasis):

    def __init__(self,
                 mesh,
                 elem,
                 elements,
                 mapping=None,
                 intorder=None,
                 quadrature=None,
                 side: int = 0,
                 dofs=None):
        self.mesh = mesh  # required by _normalize_elements
        elements = self._normalize_elements(elements)

        all_facets, counts = np.unique(m.t2f[:, elements], return_counts=True)
        facets = all_facets[counts == 1]
        tind = m.f2t[:, facets].flatten()
        ix = np.in1d(tind, elements)
        super().__init__(
            mesh,
            elem,
            mapping=mapping,
            intorder=intorder,
            quadrature=quadrature,
            facets=facets,
            _tind=tind[~ix] if side == 1 else tind[ix],
            dofs=dofs,
        )


m = MeshTri().refined(3)
sfbasis = SubdomainFacetBasis(m, ElementTriP1(), m.elements_satisfying(lambda x: (x[0] - 0.5) ** 2 + (x[1] - 0.5)**2 < 0.25**2))

[kinnala]

Alright, seems there was a "bug" remaining in BoundaryFacetBasis and it did not construct the normals based on _tind even though that was how I originally meant it to work. I'll need to look into it.

[kinnala]

My take on this is here: #861

[kinnala]

• Mesh depends on (Element, Dofs, Mapping) triplet.
• Basis depends on Mesh, and possibly different (Element, Dofs, Mapping) triplet.

The first triplet (Element, Dofs, Mapping) inside Mesh defines the geometrical shape of the Mesh (e.g., boundary could be quadratic and not piecewise linear if you use ElementTriP2, or the mesh could be periodic if you use ElementTriP0).

The second quadruplet (Mesh, Element, Dofs, Mapping) defines the finite element basis that you use to discretize PDE's.

Originally we supported piecewise linear meshes, essentially (ElementTriP1, Dofs, MappingAffine), and the support for arbitrary Element classes was added only in 3.0.0.

[kinnala]

Examples that use these features in a nontrivial manner:

• https://scikit-fem.readthedocs.io/en/latest/listofexamples.html#example-42-periodic-meshes
• https://scikit-fem.readthedocs.io/en/latest/listofexamples.html#example-31-curved-elements

[kinnala]

Let me add that 3.0.0 was a huge refactoring and it took me really several tries to make it backwards-compatible to even some extent. Moreover, the features are still somewhat experimental, e.g., I'd like to be able to solve FEM problems on a surface but it is not currently possible because of some shortcomings inside Dofs and missing stuff in Mapping.

[gatling-nrl]

Now we don't need self.mesh = mesh before a call to normalize*.

Ok, I see it now. You need to call self._normalize_elements() before super().__init__() but self._normalize_elements() was depending on running after super().__init__()

[gatling-nrl]

Digging deeper I realize part of the problem. The way you have them written now, _normalize_facets and _normalize_elmements are defined in AbstractBasis but they actually belong to Mesh. You can tell because you got rid of all the references to self in both functions.

[gatling-nrl]

Looking now at your code. What is tind/_tind and f2t/t2f?

[gatling-nrl]

I'm working on adding some stuff to docstrings. I thought the t had something to do with topology/connectivity. Triangles makes a lot more sense.

[kinnala]

I started with only triangles and tetrahedra, however, now it's pretty much everything and I think "topology" sort of generalizes the idea as described here:

• https://mfem.org/mesh-format-v1.x/
• http://www.nutils.org/en/stable/tutorial/#topology-vs-geometry

[gatling-nrl]

While we're on this subject, do you have any links for reading on AffineMappings and IsoparametricMappings, or is this something I can reasonably ignore right now? It isn't something I've ever studied before, but it is where the normals seem to be getting calculated so I wasn't sure how important understanding it was to SubdomainFacetBasis. I know nothing about these mappings.

[kinnala]

Those are pretty crucial to how finite element method is programmed using a reference element. I think pretty much any FEM book/course material will cover affine mappings. Isoparametric mapping uses the finite element basis itself to construct the mapping and, while common, it might not be covered in detail. I'm not sure what I would recommend though. I read tens of books so I'm not a good judge anymore.

[gatling-nrl]

Ah! find is just a variable for "facet indicies"!! .. I hope...

[gatling-nrl]

Is

mesh.f2t[:, facets].flatten('F')

equivalent to

mesh.f2t[:, facets].T.reshape(-1)

just more efficient?

Edit:
They are very nearly identical. However, per the numpy docs, .reshape tries to avoid making a copy, preferring instead to make a view. .flatten always makes a copy.

A = np.empty([2000,400,200])
%time a = A.T.reshape(-1)

Wall time: 3.51 s

%time b = A.flatten('F')

Wall time: 3.45 s

(a==b).all()

True

[gatling-nrl]

I see, so the way you intend this to work is that

        super().__init__(
            mesh,
            elem,
            mapping=mapping,
            intorder=intorder,
            quadrature=quadrature,
            # the next three arguments are 1d array of the same length.
            facets=facets,  # this first one specifies which facets are in the basis
            _tind=_tind,  # this one specifies which triangle to use for the trace values
            _tind_normals=_tind,  # this one specifies which triangle to point the normal out of
            dofs=dofs,
        )

Edit: the comments were wrong, but they're correct now.