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
