ffc team mailing list archive

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

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

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 <walker@xxxxxxxxxxxxxxx> wrote:
>
> On Sun, 23 Aug 2009, Kristian Oelgaard wrote:
>
> > Quoting Shawn Walker <walker@xxxxxxxxxxxxxxx>:
> >
> >> 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.
> >
> > 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
>
> 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
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/ffc-dev