dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

On 24/03/11 04:38, walker@xxxxxxxxxxxx wrote:
> 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.  

I remember!

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

Thanks. It's very helpful to hear your experience.

I'm wondering if it's all as simple as adding the geometric info to
ufl::cell. It makes sense that a cell can fully describe its own
geometry.  We could have:

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

which would be a general case in which a ufc::cell would be queried from
inside tabulate_tensor to provide the Jacobian, etc. For efficiency, we
could also have

  element = FiniteElement("affine triangle", 1, 1)

which would assume an affine map and therefore not query ufc::cell for
the transformation data.

Garth

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