dolfin team mailing list archive

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

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.

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
>