dolfin team mailing list archive

[Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx>]

On 24 March 2011 18:12, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
>
>
> 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.

So we could change the parameter 'degree' to 'mapping' (or equivalent)
instead and let it be an object of some sort which contains all the
information that we need in the form compilers.

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