dolfin team mailing list archive

[Peter Brune <prbrune@xxxxxxxxx>]

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.

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

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