Investigate ex21 and add new convergence test

[kinnala]

Checked the eigenvalues in ex17 and they seem to be a bit high for my taste, so I'm not sure whether the equations are correct or not.

Found an analytical test case from https://cfwebprod.sandia.gov/cfdocs/CompResearch/docs/elast.pdf.

Todo: implement a convergence test.

[kinnala]

Either this was renamed to ex21 at some point or there was a typo.

[gdmcbain]

Picking up here from the duplicate #318.

I recall that I also noticed some anomalies in gdmcbain/fenics-tuto-in-skfem#3 in translating the FEniCS tutorial 06 elastostatic cantilever problem. There it was FEniCS that was inaccurate but the anomaly went away when its mesh was refined.

I'm a bit worried by this issue. I'll try to get to bottom of it. I assume it's not scipy.sparse.linalg.eigsh, although there are different ways to use this; I'll check that by trying a direct solver on a big machine, on a smaller problem if required. What else? Reduced quadrature? Basis functions or their gradients? Definition of linear elastic forms? Everything, I guess, and all of them could have further consequences! I'll look into this, as a matter of urgency.

[gdmcbain]

Ta for Day & Romera (2004). Rather than exact solutions to a continuous problem though, I'm particularly interested in exact or at least accurate eigensolutions of the discrete finite element problem on a definite mesh with definite elements.

There's some on a very simple geometry and mesh in:

• Ghodge, V., Bhattu, A. & B.Patil, S. (2018). Vibration analysis of beams. International Journal of Engineering Trends and Technology, 55, 81–86. doi:10.14445/22315381/IJETT-V55P216

It's a 550 × 50 × 5 mm cuboidal cantilever, fixed at one of the small ends, meshed with 44×4×1 hexahedra. I presume these are S2 #314. The eigenfrequencies quoted for this very coarse mesh are similar to the analytical values from one-dimensional Euler–Bernoulli beam theory. I've had these checked in ANSYS and it gives the same values for 'linear' elements which I presume are the same as ElementHex1. I'll attach the scikit-fem implementation here once I tidy it up; however, in brief, it's giving the same sort of misbehaviour as ex21. (Whereas, as noted in the duplicate #318, ex03 seems fine.)

[gdmcbain]

An even coarser example:

• Petyt, M. (1990). Introduction to Finite Element Vibration Analysis. Cambridge: Cambridge University Press, doi:10.1017/CBO9780511624292

Example 5.3 (§5.8) is a cube fixed on one side, reduced to a quarter (½×½×1) by the two planes of symmetry, and meshes with two cubical ElementHex1. Table 5.4 gives the first four modes (swaying, torsion, longitudinal, swaying) for the FEM and exact; the first three are within 10%.

[gdmcbain]

I'll also check the thermal eigenvalues of beams.msh and its refinements; that would be a good test, touching everything except the linear elastic bilinear form.

[gdmcbain]

…and the mass bilinear form, since the problem is scalar rather than vector. So also the element is just ElementTetP1 rather than ElementVectorH1(ElementTetP1()).

The thermal version, beams_thermal.py.txt, seems to behave fine: with -r 0, -r 1, and -r 2, the thermal eigenvalues are:

• [ 0.6078574 1.84443219 2.57644335 6.70555591 7.16710113 14.33344489]
• [ 0.59104062 1.7761232 2.54165659 6.58790639 6.85026825 13.8012704 ]
• [ 0.58206895 1.74027063 2.52427733 6.53693403 6.69840297 13.58157918]

Here's the fundamental thermal mode on the second refinement.

[gdmcbain]

I assume it's not scipy.sparse.linalg.eigsh, although there are different ways to use this; I'll check that by trying a direct solver on a big machine, on a smaller problem if required.

Back on the original elastic problem on the original mesh, I repeated the eigencalculation with

L = eigh(Kint.toarray(), Mint.toarray())[0][:len(L)]
print('Natural frequencies (eigh):', np.sqrt(L) / 2 / np.pi)

and obtained

Natural frequencies (eigh): [ 80.17036319 179.80324866 214.18727331 215.65390949 272.58643357
 297.40884306]

which is identical to the results of eigsh.

Similarly with .refine(1),

Natural frequencies (eigh): [ 61.02158325 135.15096073 149.1143932  182.71785978 210.27388124
 230.48736141]

identical to eigsh.

(Recall from #318 that the fundamental frequency is estimated as 80, 61, 47, … on the 0th, 1st, 2nd, … refined meshes.)

[gdmcbain]

As a test of skfem.model.elasticity.linear_elasticity, I varied the refinement on an elastostatic example, ex11, from 0 (not refining) to 4 and measured the extremal nodal displacement in the y or z direction (these being identical by symmetry) with

np.linalg.norm(u[ib.nodal_dofs], np.inf, axis=1)

The results:

refine	displacement
0	0
1	0.03461538
2	0.03078956
3	0.03066449
4	0.03055899

which seems reasonable.

N. B.:—This example uses ElementHex1 rather than ElementTetP1.

[gdmcbain]

Modifying ex11 to use ElementTetP1:

m = MeshTet.init_tensor(*[(0., 1.)]*3)
m.refine(4)
e1 = ElementTetP1()    

gives:

refine	v	w
0	0	0
1	0.02899299	0.02562738
2	0.02869832	0.02799793
3	0.02997764	0.02962972
4	0.03049255	0.03031104

which is not quite as good as the original hexahedra but still clearly convergent to the same limit.

[gdmcbain]

Here's a script to reproduce the example of Ghodge et al (2018): ghodge.py.txt but using ElementHex1 (instead of ElementHexS2 #314). Call with python ghodge.py to get the 44×4×1 MeshHex, and optionally pass --refine n to multiply the numbers of cells in the x and y directions by 2**n.

It's actually not so bad:

refine	f0	f1	f2
0	25.78368761	136.29237714	161.61759815
1	17.80901038	111.57197041	134.67460098
2	15.07335499	94.41605083	134.23872412
3	14.28790026	89.49109303	134.12060936

which does look as though it's tending towards the values [13.476, 84.458, 134.76] quoted in the paper.

The refinement here is only in the x and y directions, there's always only one cell through the z-thickness. This, I think, explains why the third mode seems to be converging faster than the first two: it's a 'yawing' mode, in the xy-plane.

Perhaps what's required for ex21 is MeshHex rather than MeshTet and then implementation of ElementHex2 or ElementHexS2. Is there anything else going on here?

[kinnala]

Some guesses:

• locking (refine more; some engineering codes [such as Abaqus] counter this via automatically applying various tricks)
• too large difference in the max entries of M and K (to counter: use e.g. different units)

I should try learning more about eigenvalue problems.

[gdmcbain]

Tricks?! How interesting. (Right now I wish that such examples weren't so interesting because I have some to solve, but, nevertheless, how interesting!)

What counts as tricks? Different quadratures? Mixed methods involving stress or strain instead of the current 'displacement' formulation? Exotic elements? Macroelements?

too large difference in the max entries of M and K (to counter: use e.g. different units)

I did try that, yes. I think scipy.sparse.linalg.eigsh is clever enough not to need that but I always think it's a good practice anyway to use appropriate units. Also the ratio E/ρ can be factored into the eigenvalues, even if one doesn't fully nondimensionalize.

I should try learning more about eigenvalue problems.

Yes, but eigsh &c. are pretty good. There is SLEPc too #236, but I think that the test of scipy.sparse.linalg.eigsh versus eigh means that the issue isn't in the numerical solution of the algebraic generalized eigenvalue problem but is already in the pair of matrices. I've also convinced myself that the forms and assembly &c. are correct, and so, yes, that it's something like ‘locking’ and in need of ‘tricks’. Surely someone has published some of these tricks and they aren't all locked up in closed-source codes; I'm continuing to look into the literature.

Thanks!

[gdmcbain]

O. K. I've found an easy fix: ElementTetP2.

The use of quadratic displacement formulated finite elements significantly improve the performance of the tetrahedral as well as the hexahedral elements.

• Benzley, S. E., Perry, E., Merkley, K., Clark, B. & Sjaardema, G. (1995). A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis. In T. Tautges (ed.), In Proceedings, 4th International Meshing Roundtable (pp. 179–191).

This works here too:

• Unrefined: f0 = 41.78
• refine(1): f0 = 41.36

There are fancier tricks (reduced quadrature, bubble functions, …) but this'll do for now; I'll modify ex21 accordingly.

[gdmcbain]

A further thought on this (I'll launch a new issue if it goes anywhere). The switch from ElementTetP1 to ElementTetP2 #322 enables a ‘reasonable’ solution on the received mesh, but the technique doesn't go much too much further. On my desktop, I can m.refine(1) (and redefine m.boundaries['fixed'] using .facets_satisfying) but no further without running out of memory.

I suspect what's happening is that eigsh, because sigma=0., is performing repeated direct solves using the stiffness matrix. If so, this eventually needs to be replaced by a preconditioned iteration. While without a shift eigsh does permit the mass matrix to be replaced by an operator which acts like its inverse, I don't think this works for the stiffness matrix. I'm wondering whether a simple expedient might be to swap the mass and stiffness matrices, set sigma back to its default None, replace the stiffness matrix with LinearOperator embodying a preconditioned Krylov solver, and then take the reciprocals of the eigenvalues. (The eigenvectors are unaffected by such a swap.)

In short, ‘I should try learning more about eigenvalue problems’ too.

[gdmcbain]

There's no need for skullduggery like swapping the stiffness and mass, eigsh accepts OPinv: LinearOperator which, in case sigma=0., acts like the inverse of the stiffness matrix.

[gdmcbain]

The 'fixed' boundary conditions are only applied to nodal degrees of freedom,

scikit-fem/docs/examples/ex21.py

Lines 28 to 32 in 0973571

	D = np.concatenate((
	dofs['fixed'].nodal['u^1'],
	dofs['fixed'].nodal['u^2'],
	dofs['fixed'].nodal['u^3'],
	))

but since the elements were upgraded to quadratic in #322, that's not all of them.