ffc team mailing list archive

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

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