dolfin team mailing list archive

[walker@xxxxxxxxxxxx]

Hello Anders and Garth.

You may not remember me, but I had some interest way back in seeing this
higher order geometry stuff put in.  However, it seemed like P. Brune had
a way so I backed off.   In the meantime, I decided to stick with my own
FEM toolbox setup and continued developing it.  I have higher order
elements implemented (in principle for any finite element that is defined)
and I can do this for 1-D curves in 3-D and 2-D surface meshes in 3-D
(which is another thing I would like to see in FEniCS), in addition to 2-D
and 3-D curved meshes.  I can even compute things like the total curvature
vector (i.e. 2nd fundamental form).  I even implemented some code
generation in my setup because I think you have a good idea there.  But my
goal is not to duplicate FEniCS (there is no way I could do that).  I just
needed these basic components for my own work.

Anyway, I would like to describe how I did it, b/c the lessons I learned
may be useful to you.

1. Garth is correct.  All local geometric transformations should be
computed separately from the tabulate tensor routine.  It has been awhile
since I perused the FEniCS code, but I assume that is still the same.  I
repeat, you should NOT compute the jacobian (or derivatives of the
jacobian, or the normal vector, etc...) in the tabulate tensor.  This has
the following advantages: (a) the code is more modular; you separate the
geometry transformation from the element calculation; (b) you can *reuse*
the jacobian calculation for each bilinear form that has to get assembled,
instead of recomputing it every time you enter tabulate_tensor.

I personally found this to work great for me.  And I would like to
emphasize that I do not see a 10-fold slow-down in assembly time (for a
fixed quad rule); I don't remember the timings, but it seemed to be
negligible b/c this calculation is done in cache.  Again, the local
mapping transformation for each element should have *nothing to do* with
the bilinear forms, etc.  It is a separate calculation (of course, you
need to take into account Piola transform when appropriate, but this is
minor).

2. When the user wants to specify the type of "higher order" geometry, I
would suggest the following.  You should have the user define a
"Geometric_Element" just like when they define regular finite elements,
i.e.

Geometric_Element = FiniteElement("triangle",1,1)

Or something like that (I don't remember your exact syntax).  In fact, you
could even make "Geometric_Element" a keyword.  And if the user does not
define it, then you default to first order lagrange continuous element
(which is what you do now).

I think you have everything already to do this.  The Geometric_Element is
just a special case of a regular finite element, except the element is
*always* the reference element.  You should be able to reuse all of your
FFC code generation.  Of course, computing higher order derivatives of
basis functions will be more complicated (i.e. the chain rule will give
more terms b/c of the nonlinear map), but I assume you are doing most of
that symbolically within Python.

3. One last thing I learned in doing this was I also separated the
transformations of the basis functions.  So the order of operations in the
assembly loop is:

(a) get the index for the current element.
(b) compute geometric transformations.
(c) compute basis function transformations.
(d) compute local element tensors.
(e) insert into global matrix.

At least for my setup, this gave a nice separation.  In fact, the tabulate
tensor routine is almost trivial since all the transformations are done.

I hope this info is helpful.  I am not trying to imply that my setup is
better; but maybe this gives a way for how to proceed.

p.s. is it possible yet to have finite element spaces defined on a
subdomain of codimension 1 simultaneously with other codimension 0 spaces
and be able to compute bilinear forms involving mixed elements?  I noticed
you have some support for integrating on subdomains, but does that include
this???

- Shawn

On Wed, March 23, 2011 7:27 pm, Anders Logg wrote:
> On Wed, Mar 23, 2011 at 11:24:37PM +0000, Garth N. Wells wrote:
>> We've heard about this work, but never seen it ;).
>> Can you post it somewhere? It could provide a discussion basis on how
to
>> move this along.
>> I think that we should be careful in asking FFC to do too much. It may
better to compute J, det(J), etc, outside tabulate_tensor(...)
>> functions.
>
> Aren't those computed using codesnippets just pasted in by FFC into the
code anyway? Perhaps those codesnippets can be expanded with higher order
codesnippets?
>
> --
> Anders
>
>
>> Garth
>> On 23/03/11 21:05, Peter Brune wrote:
>> > I saw this thread and had already started running the tests I had for
this with the latest FFC to see if anything had changed recently.  I never
made isoparametry through UFL public because I could never get
>> the
>> > efficiency to be reasonable.  The situation appears to have improved
a
>> > little bit from a year ago, but not much.  With the latest FFC, it's
about 10x slower (in 2D) to assemble Poisson with P1 for the
>> test/trial
>> > space and P2 for the coordinates than the same problem with affine
coordinates.  It's even worse if you include things like
>> > Piola-transformed surface normals into the mix.  When I was working
on
>> > this I only assembled a very small fraction of the elements (around a
cylinder in a flow) with the parametric Jacobian, so it worked OK for
small problems.
>> >
>> > I believe that it would do much, much better if the optimizing
quadrature compiler in FFC supported fractions.  This is necessary because
you need to apply the inverse isoparametric Jacobian, which includes 1 /
|J|, to any basis function derivatives in the form.
>> >
>> > For reference, here are the assembly times for the iso(super, I
>> suppose
>> > is more accurate)-parametric, optimized, and non-optimized Poisson
problems in 2D on a 256x256 square (with the parametric coordinates
>> just
>> > being the regular coordinates passed through).
>> >
>> > isoparametric assembly      |        3.2017      3.2017     1
optimized assembly          |       0.34515     0.34515     1 regular
assembly            |       0.36524     0.36524     1
>> >
>> > I got really deep in FFC trying to make this work with no success,
but
>> > this was before the rewrite.  I'd be willing to declare victory on
>> this
>> > one and submit my code if someone else were willing to make it fast
enough to use.  There's also the issue of how exactly to extend the
interface to support this in an elegant fashion.  Right now I just
>> call
>> > a function that takes a form and a parametric Jacobian, runs a
Transformer on the form, and spits out a new form.
>> >
>> > - Peter
>> >
>> >
>> > On Wed, Mar 23, 2011 at 3:30 PM, Anders Logg <logg@xxxxxxxxx
>> > <mailto:logg@xxxxxxxxx>> wrote:
>> >
>> >     Peter Brune claims to have solved this by a small addition to the
>> form
>> >     language that automatically expresses the curved elements as a
>> mapping
>> >     and expands appropriately (and invisible to the user) those
>> mappings
>> >     to yield a form that may then be assembled. The higher order
>> geometry
>> >     is then expressed as a vector-field on the mesh.
>> >
>> >     Perhaps Peter can be pushed to polish up on his code and submit
>> it.
>> >
>> >
>> >
>> >     On Wed, Mar 23, 2011 at 07:46:57PM +0000, Garth N. Wells wrote:
>> >     > We haven't really looked at this. It was discussed a while
back,
>> >     but no
>> >     > one has committed much time to it. We struggled to settle on an
appropriate abstraction to push on with.
>> >     >
>> >     > Garth
>> >     >
>> >     > On 23/03/11 18:40, Douglas Arnold wrote:
>> >     > > What is the status of curved (e.g., isoparametric) elements
in
>> >     dolfin?
>> >     > > I gather they are not implemented in the main branch.   Has
>> anyone
>> >     > > done anything with this can be used?  Is there any example
>> code?
>> >     > > (For example, if you want to
>> >     > > solve the Poisson problem in a disc and get better than 2nd
>> order
>> >     > > convergence, you need to do better than polygonal
>> approximation of
>> >     > > the disc.)
>> >     > >
>> >     > >  -- Doug
>> >     > >
>> >     > > _______________________________________________
>> >     > > Mailing list: https://launchpad.net/~dolfin
>> >     <https://launchpad.net/%7Edolfin>
>> >     > > Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> >     <mailto:dolfin@xxxxxxxxxxxxxxxxxxx>
>> >     > > Unsubscribe : https://launchpad.net/~dolfin
>> >     <https://launchpad.net/%7Edolfin>
>> >     > > More help   : https://help.launchpad.net/ListHelp
>> >     >
>> >     > _______________________________________________
>> >     > Mailing list: https://launchpad.net/~dolfin
>> >     <https://launchpad.net/%7Edolfin>
>> >     > Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> >     <mailto:dolfin@xxxxxxxxxxxxxxxxxxx>
>> >     > Unsubscribe : https://launchpad.net/~dolfin
>> >     <https://launchpad.net/%7Edolfin>
>> >     > More help   : https://help.launchpad.net/ListHelp
>> >
>> >
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>
>