dolfin team mailing list archive

[Peter Brune <prbrune@xxxxxxxxx>]

It's cool!  I'm making baby steps towards being able to do this.  In the
mean time it seems that I can compute a lot of things with what I have now
(I'll send a patch once my demo is nice enough)  You pretty much pointed me
towards where to look at ENUMATH.  Have a good vacation.

- Peter

On Tue, Jul 7, 2009 at 3:20 PM, Anders Logg <logg@xxxxxxxxx> wrote:

> 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
>
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkpTrfYACgkQTuwUCDsYZdGOvQCdEjt/J8LsaGmtrgQyRzQvH4Wq
> fMQAn29qPJEi4c0Ns7rouOBwAZ6CFnoX
> =iCb3
> -----END PGP SIGNATURE-----
>
>