ffc team mailing list archive

[Shawn Walker <walker@xxxxxxxxxxxxxxx>]

On Sun, 23 Aug 2009, 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.



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


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