optional subset of elements in InteriorBasis

[gdmcbain]

It is already possible to solve a PDE in a subdomain #159 but:

I'm surprised that InteriorBasis.init does not have an optional argument for specifying the set of elements.

This would be consistent with the optional facets argument of FacetBasis.

[kinnala]

A preliminary version is now available in https://github.com/kinnala/scikit-fem/tree/subdomain-assembly.

But after starting to modify ex26 to conform to this, I became unsure whether this is something that is actually wanted there.

If you set now the argument 'elements' to something other than None, you still get a matrix with same size as the original assembly but just DOF's related to other elements ignored in the assembly phase (rows and columns are zero).

[gdmcbain]

Right. I think that this is correct behaviour for an elements keyword-argument to InteriorBasis that is consistent with the existing facets keyword-argument of FacetBasis. If the interior or facet basis passed to asm is restricted to a subdomain of elements or subboundary of facets, it should still return a sparse matrix of the same (square) shape, just with different .nnz. This is just the same as, say in ex26,

asm(laplace, FacetBasis(mesh, basis.elem, facets=mesh.boundaries['convection']))

is (683, 683) with nnz==624 while

asm(laplace, FacetBasis(mesh, basis.elem, facets=mesh.boundaries['convection'][:9]))

is (683, 683) with nnz==75. (Not that the latter matrix has much physical significance.)

I think that #166 does correctly fix the concern from #159:

Later it did occur to me that there was a bit of wasted assembly which could be avoided with such an optional argument.

I find that on branch subdomain-assembly, I can get the same correct solution for the temperature field in ex26 with

wire_basis = InteriorBasis(mesh, basis.elem, elements=mesh.subdomains['wire'])
L = asm(laplace, wire_basis)
f = asm(unit_load, wire_basis)

[gdmcbain]

But after starting to modify ex26 to conform to this, I became unsure whether this is something that is actually wanted there.

If you set now the argument 'elements' to something other than None, you still get a matrix with same size as the original assembly but just DOF's related to other elements ignored in the assembly phase (rows and columns are zero).

As to what's ‘actually wanted’, I'm not exactly sure yet as the work on coupling subdomains #163 remains nebulous. (More on that shortly.) I mean, I don't think I've decided or worked out exactly how to set up a problem with multiple subdomains or domains.

I believe the changes in #166 are exactly what's wanted to do ex26 properly, but ex26 is very much a toy. I can certainly imagine domain-decomposition style approaches in which problems are to be solved (or at least assembled) in subdomains, and then coupled again #163, but whether ex26 is usefully indicative of them remains to be seen.

I plan to construct two more examples for #163 that should (begin to) elucidate this.

• Currently ex26 does not re-solve ex17, it solves a different problem. It might be useful or of interest to re-solve ex17 by setting up the problems on the two subdomains and then coupling them. One approach is the Poincaré–Steklov operator. I might start with this because I know how to do it and think it should be pretty easy. I think that this is equivalent to introducing a Lagrange multiplier for the interfacial flux and then solving its Schur complement.
  • This would differ from ex26 which is a Dirichlet problem in solving a Neumann problem in the wire. (With consequent issue of avoiding the null space, though that's hardly insurmountable.)
• Conjugate heat transfer, extending something like the forced convection example #162 by tacking a solid conducting slab below the floor. Contrary to earlier plans in #163, this is probably done more simply, more like ex17 and less like the previous, by:
  • Solve the Navier–Stokes equation in the fluid subdomain, like #160 ex26, but then
  • solve for the temperature monolithically, as in ex17, just with the advection only active in the fluid subdomain.

The second example is certainly justified as I'd use something based on it in real applications for sure. The first one less so as it's not clear what advantage it has over the monolithic solution in ex17. Maybe required for problems with ‘contact resistance’ in which there is a jump-discontinuity in the solution across the interface.

There are lots of other ways to do both problems. If one of those turns out to be better and then ex26 isn't demonstrating one of the substeps, then it can be jettisoned.

The changes in #166 are still correct and at the very least pretty harmless. They only very slightly increase the code-complexity and do eliminate the asymmetry with the facets in FacetBasis.