← Back to team overview

ffc team mailing list archive

Re: slight modification of FFC and UFL

 



On Mon, 14 Sep 2009, Anders Logg wrote:

On Sun, Aug 23, 2009 at 10:08:57PM +0200, Kristian Oelgaard wrote:
Quoting Shawn Walker <walker@xxxxxxxxxxxxxxx>:


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.

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.

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

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

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?


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.

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
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.
:) )

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. :)

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.

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.

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.

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

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

Follow ups

References