an unsteady one-dimensional example

[gdmcbain]

Is it worth including an unsteady one-dimensional example? If not, this can be closed and can just be referred to here.

Is ex19 (‘Heat equation’) the only unsteady example so far? This is a reduction of it to one dimension. The finite element and analytical aspects are straightforward, but the plotting is different; one dimension is more different to two and three than the latter are to each other—we don't normally colour along a line for temperature in one dimension.

[gdmcbain]

[gdmcbain]

The main trick is that the plotting of the snapshot is changed from

scikit-fem/docs/examples/ex19.py

Line 117 in d08af2a

field.set_array(u)

to

scikit-fem/docs/examples/ex39.py

Lines 108 to 113 in 5b788e7

	# skfem.visuals.matplotlib.plot_meshline plots a segment for
	# each cell, separating them with a "None".
	field.set_ydata(
	np.vstack([u[mesh.t], np.full(mesh.t.shape[1:], None)]).flatten("F")
	)

This relies on the internal details of

scikit-fem/skfem/visuals/matplotlib.py

Lines 140 to 153 in d08af2a

	xs = []
	ys = []
	color = kwargs["color"] if "color" in kwargs else 'ko-'
	for y1, y2, s, t in zip(z[m.t[0]],
	z[m.t[1]],
	m.p[0, m.t[0]],
	m.p[0, m.t[1]]):
	xs.append(s)
	xs.append(t)
	xs.append(None)
	ys.append(y1)
	ys.append(y2)
	ys.append(None)

which it might be clearer to abstract, but I haven't attempted that here. (I have though replaced the for-loop with a NumPy equivalent in 66144ed.)

[gdmcbain]

Need to check generality of this. The previous uses only first two rows of m.t so should it be m.t[:2]? Need to check whether the first two points are the termini.

[gdmcbain]

Or is this only for two-point line-elements?

[gdmcbain]

So it's invoked for MeshLine

scikit-fem/skfem/visuals/matplotlib.py

Lines 130 to 131 in d08af2a

	@plot.register(MeshLine)
	def plot_meshline(m: MeshLine, z: ndarray, **kwargs):

which is MeshLine1

scikit-fem/skfem/mesh/__init__.py

Line 38 in d08af2a

MeshLine = MeshLine1

and that does only cover two-point line-segments defined by their termini.

However, it remains to consider the array of ordinates z, which might be defined on a basis other than ElementLineP1.

[gdmcbain]

Right, no, it fails for

element = ElementLineP2()

raising ValueError: dimension mismatch.

[gdmcbain]

…but that was only because there was a bunch of other stuff that was only going to work with P1; e.g.

scikit-fem/docs/examples/ex39.py

Line 56 in 5b788e7

u_init = np.cos(np.pi * mesh.p[0, :] / 2 / halfwidth)

scikit-fem/docs/examples/ex39.py

Line 82 in 5b788e7

ax = plot(mesh, u_init)

I've fixed all that now; the element can be changed from ElementLineP1 or ElementLineP2 without changing anything else.

[kinnala]

I don't know if this is needed though? Isn't ax.plot(m.p[0], z) enough? This could be just some copypasta from draw_mesh2d where appending None gives faster performance.

[kinnala]

Oh I see, this is because of the arbitrary order of m.p[0]. I think np.argsort(m.p[0]) is needed before plotting.

[gdmcbain]

Oh, when i read the old for-loop implementation, I thought it was deliberately drawing one line segment per cell, possibly so that this would work with

• unsorted .p or
• discontinuous Galerkin?

[gdmcbain]

I opened #676 to discuss whether plot_meshline needs revision.

[gdmcbain]

I don't know if this is needed though? Isn't ax.plot(m.p[0], z) enough?

It is indeed, or even using basis.doflocs and plotting all the ordinates, which is richer if using ElementLineP2. I've gone with that in bdc010b and reverted the changes to plot(MeshLine1).

That might be revisited in #676, or perhaps the resolution to that will end up being again, ‘isn't ax.plot(m.p[0], z) enough?’ and the function deprecated.

[kinnala]

Unsorted p must be supported but dG is done through plot_basis.

[gdmcbain]

I see that we don't have ElementLineDG yet. #675

[gdmcbain]

Figure:— Evolution of the fundamental cooling mode of a uniform slab with cold walls, rendered directly in matplotlib (rather than the wrapper in skfem.visuals)

[gdmcbain]

I'm thinking of following this with (besides the method of characteristics #600 for pure advection or advection–diffusion) something on a second-order IVP #531, e.g. a wave equation or possibly transverse vibration of an Euler–Bernoulli beam (ex34, 2a297e0). I don't think the second time derivative will pose any difficulties (it didn't in my in-house code years ago and scikit-fem is only better). A partial motivation for this (besides these two being of intrinsic interest) is the idea of treating strong advection with a wave-like reformulation.

That's not being held up by the present PR but will probably use a similar idiom for time-stepping and animation.

• Gresho, P. M., Sani, R. L. (1998). Incompressible Flow and the Finite Element Method. Volume One: Advection–Diffusion. Wiley. §2.7.7f ‘A wave equation method’
• Glowinski, R. & Pan, T.-W. (2019). Two decades of wave-like equation for the numerical simulation of incompressible viscous flow: A review. In B. N. Chetverushkin, W. Fitzgibbon, Y. Kuznetsov, P. Neittaanmäki, J. Periaux & O. Pironneau (eds.), Contributions to Partial Differential Equations and Applications (pp. 221--250). Springer. doi: 10.1007/978-3-319-78325-3_13

[kinnala]

Can you figure out how this could be run in the test suite? In particular, tests/test_examples.py should ideally run all examples so that they won't so easily break in the future. Maybe add some call to evolve to tests/test_examples.py?

[gdmcbain]

I've launched #678 for this.

[gdmcbain]

Subsumed by #679.