Towards less abuse of asm

[gatling-nrl]

Ok, so yeah, I like to inflict syntactic mischief with Functional and asm. Abuse, some would say. Have said, in fact.

But its just so useful for poking around at things...

Is there any appetite for an actual iterator that traverses the data structures similar to how I've been using asm?

This is a wide open question... obviously you can get quite far with what exists already. But I admit it's side effects, and it's ugly... And I'm not too clear at all on syntax or anything... its just an idea.

The first use case I was thinking of is all this plotting and poking around I've been doing. But then you can look at my train wreck in #851 and see how maybe it could be quite useful when trying to extend skfem.

[gdmcbain]

Abuse, some would say. Have said, in fact.

Ha, yes. Just in case this should be taken wrongly, ‘abuse of notation’ is a common locution in metamathematics and not at all derogatory:

In mathematics, abuse of notation occurs when an author uses a mathematical notation in a way that is not entirely formally correct, but which might help simplify the exposition or suggest the correct intuition (while possibly minimizing errors and confusion at the same time). ) (Wikipedia)

and I agree that your (ab)uses have helped simplify the exposition and are helping suggest the correct intuition. So let's continue following those suggestions towards the correct intuition.

Now, as to the idea of an iterator: interesting, not bad, maybe a bit un-Pythonic in that usually in NumPy we try to avoid explicit iteration over ndarrays in favour of built-in vectorized operations, but perhaps there might be a case to be made for its use in development as an explicit way which can be asserted to do the same thing as a less obvious einsum? Yes, I can see that.

I've been away on holidays and see that I've still quite a bit to catch up on, but before that I did begin on working to create somewhat similar figures on triangular and quadrilateral meshes, indicating quadrature points, normals to facets, &c., without recourse to side-effects of asm. This was in ‘render ElementTriRT0 field’ #824, but I think you already saw that towards the end of last year.

Rather than iterators, that accessed ndarrays of FacetBasis.global_coordinates() and FacetBasis.normals. Do we need iterators rather than arrays? If one really wants to iterate, one can always iterate over an ndarray with nditer. Are there properties or attributes that we would like to iterate over that aren't yet as easily available as .global_coordinates() and .normals?

I'll take a look at #851 and after that I hope to have more concrete suggestions. It's something that I'm very interested in anyway. We really want to be able to calculate the jump and flux across an interface between subdomains as in ex26 and that means being able to consistently reorient the constitutive facets #821 or at least know what their orientation is so that it can be accounted for.

[gdmcbain]

Imagine that ex26 were extended to solve not just in the core subdomain but also in the annulus. if a Dirichlet condition on temperature were used on the interface in both problems, the jump would be zero, but the fluxes wouldn't match. Unless we somehow managed to specify the temperature at just the right values.

Conversely, if we had specified the flux on the interface as a Neumann condition for the separate subdomains, the temperatures wouldn't match. The difference in temperature at the interface from thw two sides is called the jump. (Note here that i mean the local flux density at each point in the interface, not the integarted flux, or "current" in the electrical analogue.)

Sometimes we don't want the jumo to vanish; e.g. for thermal contact resistance, in which there is a jump proportional to minus the flux.

A specified mismatch in fluxes is also possible. This could corrspond to an interfacial source, generalizing a point source.

See my 2012 slides from the Fourth FreeFem workshop : Domain decomposition techniques for interfacial discontinuities. This is all really easy in FreeFEM because it's possible to define (n-1)-dimensional bits of the mesh (basically "boundaries" but not necessarily actually on the boundary, really just subsets of facets) that are consistently oriented.

I hadn't conceived of these interfaces as bounding subdomains, although they do in special cases like your isolated island subdomains. I thought of them mostly as separating two subdomains; not always the same thing, e.g. if three subdomains meet at a point (in the plane, or at a curve in space). I have also even had cases (in FreeFEM) in which there was an internal curvilinear (well, polygonal anyway) segment inside a two-dimensional domain. This represented a thin electrically heated sheet embedded in a thermally conducting solid. This kind of topology is easy to specify in Gmsh (I was using the old GEO syntax then, before the new Python API) and easy to represent in Gmsh's MSH formats using lower-dimensional elements. It's then easy to read in to scikit-fem as something a bit like a boundary (even though it doesn't bound anything) but the infrastructure isn't set up for this. The obvious thing is that the facets aren't consistently oriented.

It's not clear to me that #851 is well set up to handle such cases. It might not get in the way, but I would like the orientation part factored out so that I can make sense of the orientation of any continuous collection of facets. This needn't hold up your current effort, but I'm flagging that either I might have to use something else later or refactor it if I return to these interfacial discontinuity problems.

[kinnala]

Thanks @gdmcbain for an interesting comment. I think there are many ways to hack InteriorFacetBasis to use a consistent normal, however, they won't be very general. Not every collection of facets is continuous so trying to orient the normals by default is not maybe the correct approach?

[gdmcbain]

I'm starting to think that that's right, yes. We do see any example of a discontinuous set of facets in the application of SubdomainFacetBasis to a disconnected subdomain.

I'm starting to wonder whether what's really wanted for the completion of ex26 is something very like mortar elements but supported on facets of the n-domensional mesh rather than a separate (n-1)-dinensional mesh on an independent hypersurface. Thus this "mortar" would be connecting compatibly meshed subdomains; connecting but possibly also dividing with inhomigeneous jump and flux-mismatch conditions.

I'm wondering whether it would be better to proceed more cautiously, solving examples with the library as is putting all the extra code (identification and orientation of facets, calculations of local jumps and flux-mismatches) in the top level script. Then if it appears that there is something there that is going to be of general interest and that can be reasonably factored into a function, then that can be considered for inclusion in the library. This is how I was proceeding in ex26, e.g. at

scikit-fem/docs/examples/ex26.py

Lines 34 to 37 in b24b078

	annulus = np.unique(basis.element_dofs[:, mesh.subdomains['annulus']])
	temperature = basis.zeros()
	core = basis.complement_dofs(annulus)
	core_basis = Basis(mesh, basis.elem, elements=mesh.subdomains['core'])

A possible next step in that example would be to calculate the local flux out of the core to use as a Neumann condition for the annulus, then compute the jump across the interface and iterate, recrossing the Dirichlet condition for the core on the interface until the jump vanishes. Such Dirichlet–Neumann conditions don't always converge, depending on the ratio of thermal conductivities. Better is to write a function that maps core interfacial temperature to interfacial jump and then pass that to a solver—same structure as the Uzawa-CG sklver for preseure–incompressibility in

scikit-fem/docs/examples/ex30.py

Lines 96 to 101 in b24b078

	def dilatation(pressure: np.ndarray) -> np.ndarray:
	"""compute the dilatation corresponding to a guessed pressure"""
	return -B @ flow(pressure)
	dilatation0 = dilatation(basis['p'].zeros())

…or, going further, including everything in the system matrix and solving that, perhaps with preconditioned iterations, as in

scikit-fem/docs/examples/ex32.py

Lines 139 to 149 in b24b078

	def precondition(uvp):
	uv, p = np.split(uvp, [Aint.shape[0]])
	return np.concatenate([Apc @ uv, p / diagQ])
	M = LinearOperator(Kint.shape, precondition, dtype=Q.dtype)
	velocity, pressure = np.split(
	solve(Kint, fint, u, I,
	solver=solver_iter_krylov(minres, verbose=True, M=M)),
	[basis['u'].N])

But first step is orientation. #821 I think we've got a couple of techniques for that now #824 #851, but let's factor these out and clean them up and not abuse asm, which was not necessary in #824 though it didn't yet go on to calculate the one-sided fluxes.

One question to think about is what function space the one-sided fluxes and flux-mismatch should be in. It's like a Lagrange multiplier space. For techniques as described above it's most simply the same finite dimension as the jump so that the solving makes sense.

[gdmcbain]

A special case of what I'm describing here could restrict the space of core-interfacial temperatures to ℝ; I.e. to be a single undetermined value. This would be associated with the integrated flux. This would be a rational resolution of #751 and maybe the ultimate goal of #851 though I'm not exactly sure what that is.

[kinnala]

Technically the normal vector is a finite element function defined by ElementVector(ElementTriSkeletonP0()). However, I'm not happy about how these skeleton elements are implemented and I mean to do something about it.