ffc team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxxx>]

Quoting Peter Brune <prbrune@xxxxxxxxx>:

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

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.

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

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
> hard-working gräsänkling for just two more weeks. :)
>
> - Peter
>
> On Sun, Aug 23, 2009 at 3:08 PM, Kristian Oelgaard
> <k.b.oelgaard@xxxxxxxxxx>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. :)
> >
> > 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
> > > >
> > > >
> > > >
> > > >
> >
> >
> >
>