ufl team mailing list archive

[Peter Brune <prbrune@xxxxxxxxx>]

On Mon, Jul 6, 2009 at 8:34 AM, Shawn Walker <walker@xxxxxxxxxxxxxxx> wrote:

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

For the moment my examples involve creating a trivial coordinate function
and transforming it to match the desired geometry.  This may involve
projection onto a surface, or using the coordinates as parameters of a new
coordinate function. (think: warping a relatively smooth bluff into a square
mesh)  This is then easily interpolated into a discrete function of any
order.  Because we can analytically construct the Jacobian from this, we can
use the transformation machinery in UFL to fold the geometry into the form.


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

Yes, of course!  Because I intend on defining it in terms of per-subdomain
measures, we can even vary the order of the geometry from mesh-piece to
mesh-piece (by having per-order parts in the form).  My present demos don't
distinguish between affine and curved (that I'll clean up and get out) do
things like wrap a pipe around the origin, considering the "reference"
coordinates to be r and theta.  I could just as easily add the affine part
to the form.  I'll try to make this part of my cleaned-up example.


>
>
>  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'm working on some good test examples right now, but the demos I have will
involve superparametric vs. a ludicrously refined mesh.  I'm able to run
both Poisson and Stokes examples (which I'm cleaning up and error
checking).  There could be a number of ways of reading these in, from
reading in a "higher-order" mesh format directly into a function space on
the "boundary-cell" subdomain and an affine space merely defined by the
geometry as usual, to storing the coordinates separately.  Our goal is to do
some trivial projection to an "exact" geometry.


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


- Peter