Reimplementation of skfem.mesh

[kinnala]

• Added MeshHex2
• Added MeshTet2
• Added MeshHex.to_meshtet
• Added MeshHex.element_finder so that Basis.interpolator works for hexahedral meshes. Fixes Implement MeshHex.element_finder #503
• Removed Mesh.{refine,scale,translate}
• Removed Mesh.define_boundary
• Added Mesh.with_boundaries and Mesh.with_subdomains
• Changed Mesh{Tri,Quad}.refined to clear subdomains and boundaries, and not update them based on the new indexing (this will be hopefully reintroduced when resolving load auxiliary data with meshes #271)
• Fixes Steps towards a "functional Mesh" #540
• Fixes Use Mesh.refined in docs, examples and tests #541
• Fixes Tasks before releasing 3.0.0 #546

[kinnala]

Let me add that in theory it supports arbitrarily high order meshes.

[kinnala]

One thing that must be verified is that the performance cost is not too large because it's using MappingIsoparametric always in place of MappingAffine.

[gdmcbain]

assert (a-A).sum() < 1e-10

Possibly not the best test since

a.sum() == A.sum() == 0

and almost numerically also for a2, A2, …, by definition the laplacian of a constant vanishing.

Maybe

assert (a - A).nnz == (a3 - A3).nnz == (a4 - A4).nnz == 0

but (a2 - A2).nnz == 452 and has data ranging from -504.5 to 251.6.

[gdmcbain]

It's not immediately obvious to me how this distinguishes quadrilaterals from tetrahedra, but I'll step through the (obviously working) example to see.

[gdmcbain]

Oh, is it just that quadrilaterals aren't included yet and will be handled by something like

scikit-fem/skfem/mapping/mapping_isoparametric.py

Line 202 in 93c8bc1

elif self.dim == 2 and self.mesh.t2f.shape[0] == 4:

[gdmcbain]

[np.array_equal(m.facets, M.facets) for m, M in [(m, M), (m2, M2), (m3, M3), (m4, M4)]]

[True, True, False, False]

[kinnala]

Thanks for pointing these out. It seems that the quadratic Grid still has some issues. It's a bit complicated to wrap my head around it, however, at least we have the working reference MeshTri2 to compare against.

[kinnala]

...

[kinnala]

Although it's very nice approach, currently the performance difference is so large that I won't replace Mesh unless I'm able to optimize it significantly:

In [1]: %timeit basis3 = InteriorBasis(m3, ElementTetP1())                      
142 ms ± 773 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [2]: %timeit Basis3 = InteriorBasis(M3, ElementTetP1())                      
775 ms ± 5.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

There are some bottlenecks that I immediately recognize. I'll start looking into it tomorrow.

[kinnala]

TODO hexahedron save reordering