dolfin team mailing list archive

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

On 24/03/11 12:50, Kristian Ølgaard wrote:
> On 24 March 2011 13:09, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
>>
>>
>> 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.
> 
> I think there is already something in UFL that we could use for this.
> The default (affine) triangle is:
> 
> R2 = Space(2)
> triangle      = Cell("triangle", 1, R2)
> 
> Changing the degree to '2':
> 
> triangle      = Cell("triangle", 2, R2)
> 
> will imply a higher order geometry, of course we need to implement
> support in the form compilers to generate actual code.
> 

But what kind of mapping does it imply? We want to generalised beyond
using the FE basis functions for the map.

Garth

> Kristian
> 
>> 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
>>>>
>>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~dolfin
>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help   : https://help.launchpad.net/ListHelp
>>