dofs in subdomain

[gatling-nrl]

In my application (solving Laplace's equation for electric potential) I have some embedded electrodes. Some of them act as boundary conditions and some act as probes. Each electrode occupies its own subdomain in the mesh. I want to get the dofs in these subdomains to set the boundary conditions, or in the case of the probes, compute the average potential.

I've got a working solution, I just think it's ugly and maybe not that fast. Ideas for improvement? Even though facets use midpoints and so I used midpoints here too, I don't trust == with floats. Hence the np.linalg.norm is acting like a broadcastable isclose.

def get_dofs_in_subdomain(basis, subdomain):
    eps = 1e-8
    fem_mesh = basis.mesh
    subdomain_midpts = list()
    for cell_idx in fem_mesh.subdomains[subdomain]:
        topo = fem_mesh.t.T[cell_idx]
        for edge in [[0,1],[1,2],[2,0]]: # connectivity of one topo element
            subdomain_midpts.append(fem_mesh.p.T[topo[edge]].sum(axis=0)/2.0)
    subdomain_midpts = np.array(subdomain_midpts)            
    facets = fem_mesh.facets_satisfying(
        lambda x: [(np.linalg.norm(subdomain_midpts-x0, axis=1) < eps).any()
            for x0 in x.T]
    )
    return basis.get_dofs(facets)

[gdmcbain]

How about ex26?

scikit-fem/docs/examples/ex26.py

Line 34 in 7500d1c

annulus = np.unique(basis.element_dofs[:, mesh.subdomains['annulus']])

[gatling-nrl]

I think that's a lot prettier! Thanks! Is that common enough to belong in e.g. basis.get_dofs(subdomain='annulus')?

[gdmcbain]

That's an idea, yes.

Or even basis['annulus'], if we were to define the __getitem__ method. (I haven't thought that through at all.)

I'm still not sure about find_dofs versus get_dofs, whether one will eventually be deprecated. There are many more uses of the former in docs/examples but the later additions seem to prefer the latter.

[gatling-nrl]

@gdmcbain So I gave that a try, but they're not quite equivalent. basis.get_dofs returns a DofsView but annulus contains a ndarray after line 34. Later in my code DofsView + ndarray fails where DofsView + DofsView previously succeeded.

Digging around it seems like DofsView is (almost?) always treated like an array and that it can flattened as soon as its gotten and concatenated with other ndarrays in the usual ways. Am I missing something with that?

[gatling-nrl]

Personally, I would find the basis['annulus'] notation confusing because I don't think of the basis as a collection of dofs. I'd intuitively expect basis['annulus'] to return an set of elements or something like that. basis['annulus'].get_dofs() ?

[gatling-nrl]

Since DofsView behaves a lot like ndarray the following seems to work.

bottom_dofs = basis.get_dofs(fem_mesh.facets_satisfying(lambda x: np.isclose(x[1], 0))
sub0_dofs = np.unique(basis.element_dofs[:, fem_mesh.subdomains['subdomain0_name']])
sub1_dofs = np.unique(basis.element_dofs[:, fem_mesh.subdomains['subdomain1_name']])
boundary_dofs = np.hstack([bottom_dofs, sub0_dofs, sub1_dofs])
boundary_values = basis.zeros()
boundary_values[[sub1_dofs]] = 1.0 # e.g. volts
A,b = skfem.enforce(A, D=boundary_dofs, x=boundary_values)
u = skfem.solve(A, b)

[gatling-nrl]

This little bit of syntactic sugar made the app code much more readable

class MyBasis(skfem.Basis):

    def get_dofs(self, *args, **kwargs):
        if 'subdomain' in kwargs:
            mesh = self.mesh
            subdomain = mesh.subdomains[kwargs['subdomain']]
            return np.unique(self.element_dofs[:, subdomain])
        elif 'facets_satisfying' in kwargs:
            facets = self.mesh.facets_satisfying(kwargs['facets_satisfying'])
            return super().get_dofs(facets).flatten()

        else:
            return super().get_dofs(*args, **kwargs).flatten()

as in

basis = MyBasis(skfem.io.from_meshio(mesh), skfem.ElementTriP1())

b_dofs = np.hstack([
    basis.get_dofs(facets_satisfying=lambda x: np.isclose(x[0], 0)), # left
    basis.get_dofs(facets_satisfying=lambda x: np.isclose(x[0], 1)), # right
    basis.get_dofs(facets_satisfying=lambda x: np.isclose(x[1], 0)), # bottom
    basis.get_dofs(subdomain='grounded_electrode'),
    basis.get_dofs(subdomain='source_electrode'),
])

[gatling-nrl]

I'd like to suggest a couple of small changes....

...
      if facets is not None and elements is not None:
          raise Exception('only one of `elements` or `facets` should'
                          'be specified')
    if elements is not None:
        ...
        elif isinstance(elements, (list, tuple)):
            elements = np.hstack([self.get_dofs(elements=e) for e in elements])
        ...

[kinnala]

I don't understand the latter. Why would you want to put DOF indices to elements? Do you mean something like:

elif isinstance(elements, (list, tuple)):
    elements = np.unique(np.array([self.mesh.subdomains[e] for e in elements]).flatten())

in order to support

basis.get_dofs(elements=["left", "right"])

[kinnala]

Support for multiple names is in #778.

[gatling-nrl]

in order to support
basis.get_dofs(elements=["left", "right"])

yes, that's exactly what i was thinking, except I was trying to be more general and also support

basis.get_dofs(elements=["left", lambda x: x[0]>.5])

[kinnala]

That should work after #778, didn't test it though.