More elements

[kinnala]

More elements would be nice, e.g.,

• [x] Simplical elements up to p=3
• [x] Second-order hexahedron
• [x] Curved elements through MappingIsoparametric
• [x] MINI

[gdmcbain]

I've got some Python code for one-dimensional Hermite elements of arbitrary order which are handy for fourth order equations like Euler–Bernoulli beam and Orr–Sommerfeld

• Mamou, M., & M. Khalid (2004). Finite element solution of the Orr–Sommerfeld equation using high precision Hermite elements: plane Poiseuille flow. International Journal for Numerical Methods in Fluids 44(7):721–735

but also very accurate for smooth second order systems. I haven't used them since discovering scikit-fem, but my plan is to rewrite them to fit in if possible, certainly before using them again. I saw that class Element has a class variable nodal_dofs which was encouraging.

I was unable to do that in UFL at the time because it explicitly excluded one-dimensional Hermite elements. Looking that up now though, I see that that was reversed earlier this month! ‘Hermite makes sense in one dimension’

[kinnala]

I can later give a short explanation how GlobalBasis uses the Element objects.

Hermite is not H1 conforming so the reference mappings in ElementH1 do not work directly. I think the derivative DOFs need a scaling of 1/h if a reference element is used.

Another approach is to assemble the element globally. Both approaches are viable

[gdmcbain]

I've forgotten what H1 conforming means (but I'll look it up) and yes there was a row vector

np.tile(h ** np.arange(self.dofs_per_node), 2)

that arose in a number of places.

It hadn't occurred to me to ‘assemble the element globally’; I'll have to think about that too.

[kinnala]

I'm sorry, don't worry about it, it's just math jargon. H1 elements typically have only continuity whereas H2 elements have also continuity of the derivative (which is the idea in Hermite, I suppose?).

The reference element approach is totally fine also. I'll come back to this in the evening or tomorrow. It's 8.30 AM and I need a cup of coffee.

[gdmcbain]

Yes, that's right; there are 2*nodal_dofs shape functions on each element and each has one derivative equal to one at one end (and the rest zero) so continuity of the derivatives is achieved up to order nodal_dofs-1.

nodal_dofs=2 is commonly used for the beam equation in elementary treatments but it's not a lot of extra work to encode the general case. Putting nodal_dofs=1 recovers ElementLineP1.

[kinnala]

I made some modifications to the GlobalBasis-Element interface, they are still in a local branch. They enable the Element implementation to decide whether (u, du) or, possibly, (u, du, ddu) can be used as parameters in the bilinear form.

The only difference Element implementations have to take into account is how Element.gbasis is called in e.g. InteriorBasis.init. Element.gbasis can now return an arbitrary number of parameters, e.g. (u), (u, du), (u, du, ddu), or whatever.

I call this the "order" of the element and it's defined as a ntuple, Element.order. Each entry in the tuple is an integer and defines the tensorial order of the i'th return parameter of Element.gbasis.

I'm still not happy with the resulting @bilinear_form decorator as I'd like it to support mixing element with different len(Element.order)'s.

I hope to commit this as soon as I get the decorators figured out. The takeaway message is that Hermite-type elements can now have ddu-symbols when defining the bilinear form.

[gdmcbain]

The way the old one-dimensional Hermite code handled this was to assume that the basis functions are polynomials and represent them with numpy.polynomial.Polynomial. The derivative of order m, another polynomial, is then obtained with the method .deriv(m). That worked O. K. in one dimension, but I'm not sure how well it would generalize to more. The basic idea though was that the basis functions were objects which could be asked for whatever order derivative rather than an explicit finite sequence of the first few derivatives.

The code:

    @property
    def basis(self) -> List[np.polynomial.Polynomial]:
        '''return basis functions on canonical element

        with len == d * 2, degree == 2 * d - 1, where d =
        self.dofs_per_node

        The i-th needs to be scaled by h**(i % self.dofs_per_node),
        where h is the length of the element.

        '''

        d = self.dofs_per_node
        coef = np.linalg.inv(
            np.vstack([
                np.hstack([np.diagflat(factorial(np.arange(d))),
                           np.zeros((d,)*2)]),
                np.cumprod(np.vstack([np.ones(2 * d),
                                      toeplitz(np.zeros(d - 1),
                                               np.arange(2 * d))]),
                           0)])).T
        return list(map(partial(np.polynomial.Polynomial,
                                domain=self.NODES, window=self.NODES),
                        coef))

I think it should be easy enough to get from that to the fixed-length version anyway. I don't know how often third or higher order derivatives arise in applications.

[kinnala]

What I described earlier (support for ddu's in @bilinear_form) was added in 52bee1aa3ba72970154abd6ee1adb228e7cd99bc. I'm happy that almost everything in this commit made things simpler and many parts ended up with less LOC. Only exception are the decorators: I'm slightly unhappy about how those turned out.

[kinnala]

To my surprise, assembling curved elements worked almost without modifications. I only had to add a quadratic nodal element (now ElementTriL2, where L stands for Lagrange). Check out this example for some wonky curved element action.

from skfem import *
import numpy as np
from copy import deepcopy

p = np.array([[0.0, 0.0],     #0
              [1.0, 0.0],     #1
              [0.0, 1.0],     #2
              [-1.0, 0.0],    #3
              [0.0, -1.0]]).T #4
t = np.array([[0, 1, 2],
              [0, 1, 4],
              [0, 3, 4],
              [0, 2, 3]]).T

m = MeshTri(p, t)
m.draw()

M = deepcopy(m)

def tri_2nd_order(m):
    import numpy as np
    from copy import deepcopy
    M = deepcopy(m)
    offset = np.max(m.t) + 1
    M.t = np.vstack((M.t, M.t2f + offset))
    px = np.sum(0.5*m.p[0, m.facets], axis=0)
    py = np.sum(0.5*m.p[1, m.facets], axis=0)
    M.p = np.hstack((M.p, np.array([px, py])))
    return M

M = tri_2nd_order(m)

M.p[:, 5:] = M.p[:, 5:]*np.sqrt(2)

e = ElementTriL2()
mapping = MappingIsoparametric(M, e)
basis = InteriorBasis(m, e, mapping, 3)

@bilinear_form
def laplace(u, du, v, dv, w):
    return du[0]*dv[0] + du[1]*dv[1]

@linear_form
def load(v, dv, w):
    return 1.0*v

A = asm(laplace, basis)
b = asm(load, basis)

D1 = m.boundary_nodes()
D2 = m.boundary_facets() + m.p.shape[1]
I = basis.complement_dofs(D1, D2)

x = 0*b
x[I] = solve(*condense(A, b, I=I))

if __name__ == "__main__":
    mv, xv = basis.refinterp(x, 3)
    mv.plot(xv, smooth=True)
    mv.show()

Maybe there should be some additional features to deal with the second-order meshes? The way I do it in this script is kinda inconvenient.

Note: For some mesh/refinement combinations the visualisation fails for some unknown reason, gives warning on zero Jacobian determinant. Maybe more Newton iterations are required in the inverse isoparametric mapping?

[gdmcbain]

Very nice.

I believe that it's possible to generate six-noded curved-triangular meshes like this in Gmsh and read them with meshio. Might have to check the convention for the ordering of the six nodes. I'll try this later.

On nomenclature: so are P2 elements restricted to those in which the edge-nodes are at the midpoints? I'm not sure that I'd realized that distinction previously. I hadn't heard of L2 elements, but this name might risk confusion with the L2 space, norm, projection, etc. I think that that L is for Lebesgue. I know that there are some kinds of non-Lagrangian elements (e.g. Hermite, Raviart–Thomas) but thought that P1, P2, …, Q1, Q2, … were all Lagrangian in the sense that the nodal degrees of freedom corresponded to ordinates and the shape functions were the corresponding polynomial cardinal basis.

An upcoming use of six-noded elements might be as the velocity part of the Taylor–Hood pair for Stokes flow―I saw fluid mechanics mentioned for new examples in #31. I think that they're usually P2. Does the isoparametric mapping change that?

I note that there's some diversity in how other codes read in meshes for higher order elements. FreeFem++ always reads in three-noded meshes (in 2-D, four for 3-D) and I believe inserts the midpoints if required for P2. Escript-Finley, on the other hand, insists on the six nodes already being present, e.g. in a Gmsh MSH file. It might be nice to support both options if possible.

[gdmcbain]

Higher-order meshes are generated in Gmsh with -order int on the command line, with int up to 5.

I wasn't able to get Gmsh reproduce exactly the mesh above from a GEO file, but here's a hand-modified one that I think is valid: kite2.tar.gz.

Of course this won't go into MeshTri.load at present since MeshTri.meshio_type == 'triangle' but cells.keys() == ['triangle6']. The contents of M.p & points[:, :2].T and M.t & cells['triangle6'].T look sort of similar though, with a permutation of p.

[kinnala]

I agree that the naming is confusing, L2 is not any standard nomenclature afaik. In fact, I'll probably replace the current ElementTriP2 implementation with the current implementation of ElementTriL2 so that ElementTriL2 is not needed.

Note that quadratic meshes are not required by the library in order to use quadratic elements. They are just useful for describing higher-order geometry.

[gdmcbain]

Easy one: Q0 #266 .

[kinnala]

Continued in #327