ufl team mailing list archive

[Peter Brune <prbrune@xxxxxxxxx>]

I'm working on some problems with sub/super/isoparametric elements and have
very quickly implemented it entirely using UFL and a transformation of the
form to include a geometric coefficient.  This is done using the
transformation framework to append the geometric information to the parts of
the form that require transformation.

I had previously written calculations to include the coefficient space for
geometry by hand, but I talked with Anders at ENUMATH and he told me more
about the Transformer machinery.  Right now, a Poisson form with this looks
like:

cell = triangle

iso_element = VectorElement("Lagrange", cell, 4)
aff_element = VectorElement("Lagrange", cell, 1)
element = FiniteElement("Lagrange", cell, 4)

iso_func = Function(iso_element)
affine_func = Function(aff_element)

u = TrialFunction(element)
v = TestFunction(element)
f = Function(element)

b = inner(grad(u), grad(v))
K = v*f

J = dot(inv(grad(iso_func)), grad(affine_func))
detJ = det(J)

a = apply_geometry(b, J)
L = apply_geometry(K, J)
a = a*dx
L = L*dx

Where apply_geometry applies a Transformer to the form in order to include
the geometric coefficients.

The benefits of this approach are:
1. General function spaces for the geometry -- no reliance on the somewhat
contradictory concept of "extra vertices." on a simplex.
2. Uses all the already existing mechanisms for compilation and optimization

Right now I'm stuck transforming the affinely-transformed components back to
the reference and applying the map.  The ideal would be appending something
to the measure, which is then appended to the form.  This might look like:

J = grad(iso_func)
a = inner(grad(u), grad(v))*dx(0) + inner(grad(u), grad(v))*dx(1, jacobian =
J)
L = v*f*dx(0) + v*f*dx(1, jacobian = J)

With the form compiler then omitting the generation of the affine Jacobian.
Like this we can easily have the higher-order geometry only defined on, say,
boundary cells where we have a higher-order geometry defined.  Otherwise the
affine form can be used.

Thoughts on how I should go about this?  I'm still generalizing my
transformer, but have run a demos of simple forms (Poisson, Stokes) with no
real problems; now I'm moving onto what I actually want to do like this, but
improving the interface would be nice eventually.

- Peter Brune