dolfin team mailing list archive

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

On 24/03/11 14:28, walker@xxxxxxxxxxxx wrote:
> On Thu, March 24, 2011 8:53 am, Peter Brune wrote:
>> Ah, this argument, again.  I continue to disagree.  While you could
>> certainly home-roll a FEM code that does this under assumptions like
>> constant quadrature order and small meshes,  this is going to fall down
>> hard
>> as soon as you have varying quadrature degree or type, which we always do.
>> This would also keep around way more state than the current assembly
>> paradigm allows, not to mention the per-quadrature-point storage
>> requirements of a full-blown tensor.
> 
> I don't see why constant quadrature order is so terrible.  Is it really
> necessary to have varying quadrature degree?  What do you mean be more
> state info?  This is just another locally defined function.  I am a little
> unclear about what you are saying.
> 

I don't see what the 'argument' is either. What is it exactly?

>> The other problem with this approach is that you severely restrict the
>> spaces your coordinates can live in by having some oddly-defined notion of
>> "higher-order vertices.  If you have some coefficient that is the
>> coordinates instead it's automatic.  I think you definitely should
>> recreate
>> the jacobian in tabulate_tensor; it's what we do now; it's what we should
>> do
>> if we don't want to store per-quadrature-point tensors on the DOLFIN side
>> when the DOLFIN side doesn't even know about quadrature points.
> 
> I am not proposing having higher order vertices.  I am saying treat the
> geometric map like **it has its own finite element space**!  If you want
> to have a special map for each bilinear form (so you can have varying quad
> degree) then fine.  But keep in mind that once you have a non-linear map
> for the geometry, you will (in general) not be able to compute the tensors
> exactly.  So you will need to fix the quad degree somehow.
> 

I agree with this. Forgetting forms and PDEs, we should be able able to
have a pure geometric representation. Say, for example, that I have a
mesh that uses some mapping and I want to know (possibly approximately)
its volume. The mesh and cells need to know all the geometry. Also, say
I want to know in which cell a point (in the real coordiantes) lies.
This is not related to forms.

Garth


>> Arguments about the interface can come later.  I'm sanity checking the
>> code
>> I updated to the latest FFC with a couple test cases and should have it in
>> a
>> repo somewhere public shortly.
>>
>> - Peter
> 
> I think I will let other people figure out how to put this into FEniCS. 
> Obviously, my experience is limited.  However, my code is available on my
> web page: https://www.math.lsu.edu/~walker/FELICITY.html
> 
> I am NOT advertising my code.  I am just offering it as an example of what
> to do or what NOT to do.
> 
> If you are curious about how the generated code looks, please have a look.
>  All you need to have is MATLAB and a C++ compiler.  So the instructions
> are:
> 
> 1. Open MATLAB.
> 2. type "mex -setup" to pick the C++ compiler (this is easy).
> 3. unzip my code to a directory.
> 4. within MATLAB, change to that directory.
> 5. run this script: "FELICITY_paths.m".  This will add all of the
> appropriate paths and unit tests.
> 6. then run this script: "test_FELICITY.m".  This runs all the unit tests
> I have so far.
> 7. It you want to see the code that is generated, then just go to the
> Unit_Test subdir under the Matrix_Assembly dir.
> 
> Steps 1-6 should not take more 5 mins (maximum).  I have tested it on
> windows and LINUX.  If it doesn't run, please let me know.
> 
> - Shawn
> 
>> On Thu, Mar 24, 2011 at 7:09 AM, 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.
>>>
>>> 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
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>
>>
> 
>