dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jul 07, 2009 at 02:12:51PM -0400, Shawn Walker wrote:
>
> On Tue, 7 Jul 2009, Peter Brune wrote:
>
>> 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
>
> Actually, I am not that aware of the innards of FFC.  I am a humble man!  
> :)    But I thought that since the application of the jacobian is already 
> structurally there, having a variable jacobian is not such a big deal.  
> We need an FFC person to come to our aid!  I was also hoping Anders would 
> comment on this.

I'm on vacation so I only have time to reply when something is either
(1) extremely urgent or (2) simple to respond to. This is neither so
it will have to wait... :-)

-- 
Anders


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

> _______________________________________________
> UFL-dev mailing list
> UFL-dev@xxxxxxxxxx
> http://fenics.org/mailman/listinfo/ufl-dev

Attachment: signature.asc
Description: Digital signature