encapsulating Taylor-Hood

[gdmcbain]

I'm thinking about how to structure a refactoring of some of the common code for Taylor–Hood elements, used for Stokes and Navier–Stokes equations, e.g.

• ex18 ‘Stokes equations’
• ex24 ‘Backward-facing step’
• ex27 ‘Navier–Stokes equations’

as well as work in progress like #171 towards the three outstanding issues on these equations #169 #170 #240.

Each of these begins by repeating the same code to define the elements and interior bases.

scikit-fem/docs/examples/ex18.py

Lines 13 to 16 in f5cc86f

	element = {'u': ElementVectorH1(ElementTriP2()),
	'p': ElementTriP1()}
	basis = {variable: InteriorBasis(mesh, e, intorder=3)
	for variable, e in element.items()}

There are some other common operations too. I wonder to what extent these commonalities might be encapsulated.

They might need to be parameterized for intorder (as mentioned in #240) and also for the type of mesh; although all the examples thus far (excluding ex29 the one-dimensional Orr–Sommerfeld equation which is a bit special) all use MeshTri, I think they should also work for MeshQuad, MeshTet, and MeshHex.

Given a mesh, would it make sense for a function or class to automatically generate the appropriate pair of elements and interior bases for velocity and pressure? Something like if isinstance(mesh, MeshTri): … elif instance(mesh, MeshQuad): …?

Would this sit in skfem.models?

Or could this be getting too specialized and should it be left outside the scikit-fem package?

[kinnala]

I've actually experimented with related ideas in some early prototypes of the code (back before I was using Matlab). The idea was basically similar to ElementVectorH1: to have something like ElementComposite which takes as a parameter two elements e = ElementComposite(ElementVectorH1(ElementTriP2), ElementTriP1) and then initializing basis objects only for this composite object e. Then creating matrices could be done directly for e via asm(bilinf, InteriorBasis(e)) where the definition of bilinf includes unpacking the unknowns from the "vectorial unknown" as ux, uy, p = u. I don't have any reference code at hand but it is definitely doable.

[kinnala]

#332 adds support for composite elements.

Given a mesh, would it make sense for a function or class to automatically generate the appropriate pair of elements and interior bases for velocity and pressure? Something like if isinstance(mesh, MeshTri): … elif instance(mesh, MeshQuad): …?

I think it might have educational value to e.g. have predefined stable elements for mixed Poisson, Stokes, incompressible elasticity, etc. I suppose for now skfem.models would suffice.