dolfin team mailing list archive

[Peter Brune <prbrune@xxxxxxxxx>]

Reversing the affine transformation is a temporary measure; it's because I
haven't modified FFC to stop generating the affine jacobian.

It all depends on what level things are done at as to what gets complicated;
It seems that generating this code gets complicated either way, so here's
how it looks:

My present musings work entirely within UFL, and I can easily add geometric
terms at the form level like this, but will only work with
affinely-transforming elements -- a large class, but we're going for the
general case.  There are two ways we could go from there -- the one I was
considering with this question was keeping this in UFL and extracting the
transformation type from the basis functions (for non-affine basis
functions, at least).  This appears to be possible, but really quite a lot
of code, as done in the likes of the quadraturetransformer in FFC.

The other option is to generate code for the jacobian in FFC/Syfi/etc.  I
examined the generation code and determined that the second is probably the
best option for generality's sake.  However, it's still not pretty.  There's
a lot of code wrapped up in the affine transformation from FIAT onward, as
you are probably quite aware. :)

- Peter

On Tue, Jul 7, 2009 at 10:50 AM, Shawn Walker <walker@xxxxxxxxxxxxxxx>wrote:

>
> On Tue, 7 Jul 2009, Peter Brune wrote:
>
>  Ah, let me explain further: The reason I have the two vector spaces --
>> affine and iso, is so that I can
>> reverse the affine transformation before applying the isoparametric
>> transformation.  Today I found a small
>> bug in my isoparametric transformer, and now I'm sure stokes works; the
>> pretty example showing the
>> differences between refined affine and unrefined isoparametric is coming
>> along!  We want to get rid of having
>> to transform back to the reference using the affine coefficient space
>> before this is ready to be used.  I
>> have been looking through UFL and FFC for a good way to drop in a new
>> jacobian.  Passing the new jacobian to
>> FFC is really easy using metadata... using it is harder.
>>
>
> I'm a little unclear why you need to reverse the transformation.
>
>  So, anyone else have any tips on that?  Here's the MAJOR issue I can see
>> cropping up -- I can pepper the form
>> itself with the new Jacobian in order to transform function (derivatives)
>> affinely, but other transformations
>> (co/contravariant Piola anyone?) will not work like this.  I NEED to give
>> FFC something that can be
>> manipulated like the pullback/pushforward between geometries to use the
>> more interesting element types, or
>> even more interesting forms.   I need to be able to extract information
>> like how a given basis function (and
>> its derivatives) transform at the UFL form-transformation level, OR I need
>> to be able to change the Jacobian
>> (per-measure).
>>
>> Peter
>>
>
> It sounds like this is more complicated than it needs to be.  But maybe I
> am missing something.
>
> Currently, in the FFC generated form file, the jacobian inverse terms enter
> into the quadrature loop as constants.  All one has to do is replace those
> constants by the variable jacobian terms that are determined by your local
> map.   I did a manual modification of this and it seemed to work.
>
> My view of this is that FFC would need to be extended to generate this
> code, which shouldn't be too bad given that you can pass a
> Function/FunctionSpace to the form file that represents the higher order
> mesh geometry.
>
> The main point here is that you can pass a Function and FunctionSpace as a
> separate object to the form file, which we missed before.
>
> - Shawn
>
>
>  On Tue, Jul 7, 2009 at 8:25 AM, Shawn Walker <walker@xxxxxxxxxxxxxxx>
>> wrote:
>>      ok.  So, this is why you had those two spaces defined in the .ufl
>> file. That way the form knows
>>      how to handle it.
>>
>>      This sounds pretty good!  :)
>>
>>      - Shawn
>>
>>
>> On Tue, 7 Jul 2009, Peter Brune wrote:
>>
>>      It becomes a vector coefficient to the form.
>>
>>      - Peter
>>
>>      On Tue, Jul 7, 2009 at 12:04 AM, Shawn Walker <
>> walker@xxxxxxxxxxxxxxx> wrote:
>>           Actually, I just realized a problem with what I was proposing.
>>
>>           When you assemble the matrix, you loop through each triangle.
>>  But the higher order
>>      mesh function
>>           is defined on a different function space.  How does the
>> tabulate_tensor routine know
>>      which
>>           triangle to access in the higher order mesh function?
>>
>>           - Shawn
>>
>>
>>      On Mon, 6 Jul 2009, Anders Logg wrote:
>>
>>           On Mon, Jul 06, 2009 at 09:34:43AM -0400, Shawn Walker 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).
>>
>>                       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.
>>
>>
>>           I'm on vacation so I also won't have much to say about this
>> right now,
>>           but the form-based approach is very appealing since (1) it is a
>> simple
>>           layer on top of existing functionality (so we don't need to
>> modify
>>           the form compilers) and (2) it is more general since we can use
>> any
>>           Function to map the mesh from the reference geometry.
>>
>>           We still have the problem of reading/storing the geometry and we
>>           concluded a while back that we didn't want to use Function for
>> this,
>>           since it would create a circular dependency between the Mesh and
>>           Function classes. But if the mapping of the mesh can be stored
>>           separately from the mesh (as a coefficient in the form) then I
>> guess
>>           this is no longer a problem. (?)
>>
>>           --
>>           Anders
>>
>>
>>      _______________________________________________
>>      UFL-dev mailing list
>>      UFL-dev@xxxxxxxxxx
>>      http://fenics.org/mailman/listinfo/ufl-dev
>>
>>
>>
>>
>>
>>