Subject: [FFC-dev] slight modification of FFC and UFL

From: Shawn Walker <[email-masked]>
Date: Sat, 22 Aug 2009 21:47:43 -0400 (EDT)

Hello.  I would like to know how hard it would be to modify FFC (and UFL) 
to do the following.

Here is the sample .ufl file:

Poisson.ufl
------------------------------------------
element = FiniteElement("Lagrange", "triangle", 2)
vector  = VectorElement("Lagrange", "triangle", 2)

v = TestFunction(element)
u = TrialFunction(element)
f = Function(element)
G = Function(vector)

a = inner(grad(v), grad(u))*dx
L = v*f*dx
------------------------------------------

And I would like to do the following in C++:

// in C++
a.G = Some_Function;

Basically, I just want to attach to `a' some external function that does 
NOT appear in the bilinear form.  How hard would that be?  I'm not sure 
what the UFL syntax should be.  This is for the higher order mesh stuff.

- Shawn

---

From: Kristian Oelgaard <[email-masked]>
Date: Sun, 23 Aug 2009 11:28:46 +0200

> Here is the sample .ufl file:

I'm sure it's possible, but does it make sense that a has a function attached on
which it does not depend?
The code that will be generated from this form, will not use G in any way, so
how do you plan on using it? For higher order mesh stuff, should one not just
use:
triangle = Cell("triangle", 2) # Cell() is defined in geometry.py in UFL

and then define the finite element as:

element = FiniteElement("CG", triangle, 1)

maybe we should/could attach the function G to the cell if you intend to use
this to compute the geometry?

Kristian

> - Shawn
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

From: Shawn Walker <[email-masked]>
Date: Sun, 23 Aug 2009 10:55:25 -0400 (EDT)

> this to compute the geometry?

Yes, that is exactly what I want to do!  I'm not sure completely about the 
Syntax though.  Would this make sense:

--------------------------------------
# this is just a normal lagrange polynomial over a straight triangle
G  = VectorElement("Lagrange", "triangle", 2)

# now define what the higher order element is
my_triangle = Cell("triangle", G)

# then define the element
element = FiniteElement("CG", my_triangle, 1)
--------------------------------------

That way it is clear that G influences what the shape of the triangle is, 
and it splits up the implementation in a nice way.  How about this?

- Shawn

---

From: Peter Brune <[email-masked]>
Date: Sun, 23 Aug 2009 10:20:44 -0500

I disagree!  I think that this in some sense returns us to the previous
objection with this method of defining the coordinates, in that it is a
circular definition.  What I recommend for standardizing this (and the
completion of what I've done so far with it) is attaching the coordinate
function not to the cell, but to the measure on the integral.  It could be
done using the existing metadata functionality.  This would look something
like:

coords = Function(iso_vector)

a = dot(grad(u), grad(v)*dx(0, {"coordinate_function" : coords})

or:

coords = Function(iso_vector)

a_t = dot(grad(u), grad(v))

a = a_t*dx(0) + a_t*dx(1, {"coordinate_function" : coords})

Which would allow for both an isoparametric and affine parts of the domain
as separate integrals, but would allow the same test, trial, and
coefficients to exist for both.  This would allow one to treat curved
boundaries and non-curved interiors, and is totally general -- if you want
part of your form to be radial or parabolic, you could do it.  This is
basically what I'm doing right now, only I've had to devise ways around the
present system by first inverting the affine jacobian and then applying the
isoparametric jacobian.  I apply this with the form transformation
machinery.  Also, the facet integral would have to be rewritten to
accomodate a full transformation rather than what is done now.  I've gotten
around this as well, but it took some work.  No offense, but it seems to
need a rewrite anyways, especially for the optimization. (I had it eat all
my memory and die when I tried to turn on optimization on a facet integral
that involved a tensor.)

I'm totally willing to sift through FFC to do this, but I sort of need to
know where to start.  Kristian, any hints?  (oh, and I still am looking at
the quadrature stuff, but it's been a busy summer for me here at KTH. :) )

- Peter

---

From: Shawn Walker <[email-masked]>
Date: Sun, 23 Aug 2009 12:11:49 -0400 (EDT)

On Sun, 23 Aug 2009, Peter Brune wrote:

> I disagree!  I think that this in some sense returns us to the previous objection with this method of
> defining the coordinates, in that it is a circular definition.  What I recommend for standardizing this (and
> the completion of what I've done so far with it) is attaching the coordinate function not to the cell, but to
> the measure on the integral.  It could be done using the existing metadata functionality.  This would look
> something like:

I may appear circular, but it would not be so in the implementation (i.e. 
the FFC code).  If you want, you could think of this as overloading the 
notation.

> a = a_t*dx(0) + a_t*dx(1, {"coordinate_function" : coords})

I am a little unclear on what this means.

> Which would allow for both an isoparametric and affine parts of the domain as separate integrals, but would
> allow the same test, trial, and coefficients to exist for both.  This would allow one to treat curved
> boundaries and non-curved interiors, and is totally general -- if you want part of your form to be radial or
> parabolic, you could do it.  This is basically what I'm doing right now, only I've had to devise ways around
> the present system by first inverting the affine jacobian and then applying the isoparametric jacobian.  I
> apply this with the form transformation machinery.  Also, the facet integral would have to be rewritten to
> accomodate a full transformation rather than what is done now.  I've gotten around this as well, but it took
> some work.  No offense, but it seems to need a rewrite anyways, especially for the optimization. (I had it
> eat all my memory and die when I tried to turn on optimization on a facet integral that involved a tensor.)

I think you could have affine and non-affine with the other approach too, 
but I guess having affine and non-affine interior facets would be a 
problem.  But I'm not sure how the notation you listed above actually 
deals with having affine and non-affine elements.  What exactly is the 
iso_vector?

Obviously, I am biased, but I think that having non-affine everywhere 
would be easier to implement in FFC.  That is basically what I did in that 
demo code I sent out.  And the syntax that I propose would do that.  I 
admit, it would be nice to have a variable quadrature rule over the 
different entities based on whether or not they are affine.  You might be 
able to do this with what I wrote below, but you would need some boolean 
functions telling you what was affine and what was not, and I admit this 
is wonky.  having non-affine everywhere is more expensive, but is not the 
end of the world.

- Shawn

> I'm totally willing to sift through FFC to do this, but I sort of need to know where to start.  Kristian, any
> hints?  (oh, and I still am looking at the quadrature stuff, but it's been a busy summer for me here at KTH.
> :) )
> 
> - Peter
> 
> On Sun, Aug 23, 2009 at 9:55 AM, Shawn Walker <[email-masked]> wrote:
>

---

From: Kristian Oelgaard <[email-masked]>
Date: Sun, 23 Aug 2009 22:08:57 +0200

Quoting Shawn Walker <[email-masked]>:

> objection with this method of
> > the completion of what I've done so far with it) is attaching the
> coordinate function not to the cell, but to

I think Peter is right, although I did not follow the discussion closely the
last time we had it. The implementation might not be circular, but the notation
would look like it in the ufl file, which makes the metadata approach better.
We just need to grab this information somewhere in FFC.

> as separate integrals, but would
> part of your form to be radial or
> around this as well, but it took

A rewrite of the facet integral itself or the optimisations? Did you use the
latest FFC version? If so can you send the form (or the part thereof) which
kills your memory?

> would be easier to implement in FFC.

Using affine mapping in one sub-domain and non-affine in another should be quite
easy using the syntax Peter does above. That also take care of the 'variable
quadrature rule' that you hint at below.

That is basically what I did in that
> know where to start.  Kristian, any

Well, if we go for the metadata approach, we should probably look for the
'coordinate_function' and set some flag that tells the transformers to generate
transformations using this function. Perhaps we should give Anders a chance to
throw in his two cents before doing anything drastically? That will give you a
last chance to go out and enjoy the Swedish summer before it ends. :)

Kristian

---

From: Shawn Walker <[email-masked]>
Date: Sun, 23 Aug 2009 16:23:59 -0400 (EDT)

> We just need to grab this information somewhere in FFC.

I'm fine with whatever implementation gets put in.  I just want to be able 
to do the following thing.  If I have a fluid domain that is deforming 
with some velocity field, then I want to be able to update the higher 
order mesh function with the velocity field (maybe I have to interpolate 
to compute the iso_vector function, but that is ok).  The only other 
thing that should be done, is the new higher order mesh function should 
be passed to the form, i.e. updating the mesh function should not require 
the bilinear/linear form to be redefined/recreated.

I'm not clear on the notation that Peter has.  Peter, is what I say above 
possible?

- Shawn

---

From: Peter Brune <[email-masked]>
Date: Sun, 23 Aug 2009 15:32:30 -0500

The form that was causing the problems was G2 navier stokes with
higher-order geometry AND penalty-enforced slip in it.  Oddly enough, it was
the slip (or to be fair, the penetration) that killed it.  And... I retried
it and it seems ok.. well then...

Anyways, yeah, Shawn, it should be possible -- you've still got the Function
defining your coordinates, which you can treat however you like, such as
solve->project->solve->whatever, and the operator should remain the same
barring you totally messing up your mesh and having to recreate it.

Kristian, as to your recommendation; unfortunately, I've already got the
metadata in the compiler, and am prying apart the quadraturetransformer and
quadraturegenerator looking where to generate this code, as well as writing
a prototype where I generate the higher-order jacobian using the
quadraturetransformer outside as a test.  I'm also looking into just
generating null affine jacobians for now, and imposing the geometry at the
form level.  However, this will NOT work with say, H_div or H_curl elements
where you need to be aware of how your basis functions transform to preform
the transformations. Oh well, it's dark outside already and I'm a
hard-working gräsänkling for just two more weeks. :)

- Peter

On Sun, Aug 23, 2009 at 3:08 PM, Kristian Oelgaard
<[email-masked]>wrote:

> > > the measure on the integral.  It could be done using the existing
> metadata
> > functionality.  This would look
> > > something like:
> >
> > > allow the same test, trial, and coefficients to exist for both.  This
> would
> > allow one to treat curved
> > > parabolic, you could do it.  This is basically what I'm doing right
> now,
> > only I've had to devise ways around
> > > the present system by first inverting the affine jacobian and then
> applying
> > the isoparametric jacobian.  I
> > > some work.  No offense, but it seems to need a rewrite anyways,
> especially
> > for the optimization. (I had it
> > > eat all my memory and die when I tried to turn on optimization on a
> facet
> > integral that involved a tensor.)
>
> > > hints?  (oh, and I still am looking at the quadrature stuff, but it's
> been
> > a busy summer for me here at KTH.
> > > :) )
>
> > > > how do you plan on using it? For higher order mesh stuff, should one
> not
> > just
> > > > use:
> > > > maybe we should/could attach the function G to the cell if you intend
> to
> > use

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 24 Aug 2009 11:45:43 +0200

> it and it seems ok.. well then...

OK, I know the quadrature optimisation can be memory intensive for complex forms
and it is on my todo list (although way down) to reduce the memory usage.
I just re-wrote the optimisations and it has greatly improved w.r.t. speed and
also a little with w.r.t. memory consumption, but there's still room for
improvement in this regard.

> the transformations.

Sounds like you're already where I would be if I was to give it a go, but if you
have specific questions I'll be happy to help. I think the higher order
jacobians should be computed where we apply the transformations to the
basisfunctions or functions (in the create_....(), this is probably where you
are already?) function and not at form level.

Kristian

Oh well, it's dark outside already and I'm a

---

From: Peter Brune <[email-masked]>
Date: Tue, 1 Sep 2009 10:46:32 -0500

So, I took a short break from this, but I'm back at it today and have a
couple questions:

1. Why do we not put the metadata for the integral as part of its string
description?  This is a one-line change, but I was wondering what the reason
for reconstructing the integrals without metadata was.
2. Eventually, this all will require other things like the
evaluate_basis_derivative functions to either know about the coordinates, or
to have some sort of a second implementation when associated with a form
(potentially with multiple types of coordinates)
3. By calculating the jacobian in this create_*, ... I'm not quite getting
what you're thinking.  Do you have an organizational idea for this that I'm
not seeing?  Right now I'm appending a coordinate coefficient to forms with
coordinate metadata, and am working on getting it to generate the jacobian
in such a way I could use it.
4.  Two options: either generate the Jacobian like we do now as a
codesnippet and apply it, OR create a preprocessing transformation-aware
transformer that puts the gradient of the coefficient in ALL the right
places (because we're doing this in FFC and know what the basis functions
transform like now).  I'm leaning towards the second because I've already
sort of done it once, and I don't see how one would be slower than the
other.   This might make the above problem with the other evaluation
functions becoming harder.

Thoughts?

- Peter

---

From: Shawn Walker <[email-masked]>
Date: Tue, 1 Sep 2009 13:05:21 -0400 (EDT)

It might be easiest to have FFC generate code for a specific type of 
geometry function.  The function can be anything, but FFC would make the 
code for a specific type (say, lagrange 3rd order).  And then the user 
would `attach' this function somehow to the form.  I would also suggest 
just making all the cells have higher order geometry (for now) than try to 
worry about having an if-then in the code to decide on how to compute the 
jacobian.  That can be cleaned up later.

I would like to point out that evaluating the basis functions and 
derivatives is not completely straightforward.  Even determining what cell 
a point belongs to may be difficult, although maybe this isn't such a big 
deal.  One could just find the cell assuming everything is straight and 
just accept the possible error.

But even if you find the right cell, you need to map that point to the 
reference cell first, because you really only know the basis functions on 
the reference cell (when dealing with iso-parametric stuff).  But mapping 
that point will require a non-linear solve (unless I am missing something 
obvious).

- Shawn

On Tue, 1 Sep 2009, Peter Brune wrote:

> So, I took a short break from this, but I'm back at it today and have a couple questions:
> 
> 1. Why do we not put the metadata for the integral as part of its string description?  This is a one-line
> change, but I was wondering what the reason for reconstructing the integrals without metadata was.
> 2. Eventually, this all will require other things like the evaluate_basis_derivative functions to either know
> about the coordinates, or to have some sort of a second implementation when associated with a form
> 3. By calculating the jacobian in this create_*, ... I'm not quite getting what you're thinking.  Do you have
> an organizational idea for this that I'm not seeing?  Right now I'm appending a coordinate coefficient to
> forms with coordinate metadata, and am working on getting it to generate the jacobian in such a way I could
> use it. 
> 4.  Two options: either generate the Jacobian like we do now as a codesnippet and apply it, OR create a
> preprocessing transformation-aware transformer that puts the gradient of the coefficient in ALL the right
> places (because we're doing this in FFC and know what the basis functions transform like now).  I'm leaning
> towards the second because I've already sort of done it once, and I don't see how one would be slower than
> the other.   This might make the above problem with the other evaluation functions becoming harder.
> 
> Thoughts?
> 
> - Peter
> 
>

---

From: Anders Logg <[email-masked]>
Date: Mon, 14 Sep 2009 09:44:38 +0200

> last chance to go out and enjoy the Swedish summer before it ends. :)

First, sorry for not weighing in on this earlier.

Second, I really like the form-based approach. The more we can
leverage UFL to do things the better.

So, if Peter thinks he knows how to handle this and if Shawn thinks it
will be enough for his purposes, then just go ahead and bring on the
patches. I'll try to be more responsive from now on.

Another thing to think about when we rethink the transformations/code
snippets is the possibility of extending to assembly over manifolds
(like a 2D surface embedded in R^3). We now make assumptions on the
domain being R^2.

Just one small thing. We're trying to put together a new release with
focus on parallelization so it might be good to hold off until we've
made that release so as not to delay it further.

--
Anders

---

From: Shawn Walker <[email-masked]>
Date: Mon, 14 Sep 2009 08:10:32 -0400 (EDT)

> patches. I'll try to be more responsive from now on.

It should probably be fine.  Once this is completely working, I will 
remove that code that I put in for reading in higher order mesh data. I 
would also like to test what Peter is doing, but I will wait until his 
contributions are merged.  I have an example stiffness matrix for P2 
geometry computed by unrelated code, so it is a good independent check of 
Peter's implementation.

> (like a 2D surface embedded in R^3). We now make assumptions on the

Adding to this, I think it would be good to ensure this works for a 2D 
manifold embedded in a 3D mesh.  I assume above you just meant a 
topological 2D mesh in 3D geometry.  For problems of surface tension 
driven flow, where there are two phases, it would be nice to handle a 2D 
interface between two bulk phases.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Mon, 14 Sep 2009 14:14:14 +0200

On Mon, Sep 14, 2009 at 08:10:32AM -0400, Shawn Walker wrote:
> remove that code that I put in for reading in higher order mesh
> data. I would also like to test what Peter is doing, but I will wait
> until his contributions are merged.  I have an example stiffness
> matrix for P2 geometry computed by unrelated code, so it is a good
> independent check of Peter's implementation.

ok, sounds good.

> Adding to this, I think it would be good to ensure this works for a
> 2D manifold embedded in a 3D mesh.  I assume above you just meant a
> driven flow, where there are two phases, it would be nice to handle
> a 2D interface between two bulk phases.

Yes, topologically 2D. It's already supported by the mesh library
(that's why we can extract a boundary mesh) but there are assumptions
in a few places on matching topological and geometric dimension. Since
we can in fact assemble over facets, this already works to some degree
but we need to think through where the limitations are (I haven't).

--
Anders