ufl team mailing list archive

[Shawn Walker <walker@xxxxxxxxxxxxxxx>]

On Mon, 6 Jul 2009, Peter Brune wrote:

> 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

How exactly do you represent the geometry as a function space?  We had 
thought about this before and there was a problem (because you need a mesh 
to create a function space in the furst place).

> 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.

This is great.  There is some support for reading in higher order meshes 
in DOLFIN now.  There is even a boolean parameter for saying which 
triangles are curved and which are affine.

Please keep in mind that you may want to have triangles that are NOT on 
the boundary to also be curved.  This is necessary if the mesh is highly 
anisotropic (picture a wing with a very anisotropic curved mesh at the 
boundary).

> 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

So, you have done this in UFL?  Does that mean just the notation is setup? 
I had done a hand modification of a poisson demo that reads in a higher 
order mesh and computes the stiffness matrix on two triangles (only one is 
curved) and compare it to a stiffness matrix computed by other means.  I 
don't know if this will be useful for you; what I did was a little hacky. 
Have you thought about how the higher order mesh would be stored in an XML 
file?  I can resend the example I made (that shows this), but it should be 
in the archive on DOLFIN.

I am traveling right now, so I won't be able to say much on this, 
unfortunately.

- Shawn