dolfin team mailing list archive

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

On 24/03/11 19:10, walker@xxxxxxxxxxxx wrote:
> On Thu, March 24, 2011 1:32 pm, Anders Logg wrote:
>> On Thu, Mar 24, 2011 at 05:19:03PM +0000, Garth N. Wells wrote:
>>>
>>>
>>> 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.
>>
>> We've had this argument many times before... The thing that I've
>> always found attractive with Peter's approach is that one can
>> represent the geometry (the embedding of the mesh) as a field on top of
>> a standard mesh (with straight edges). So it doesn't require adding
>> tons of new data to the mesh class.
>>

How would this work for very complicated meshes, e.g. a car, engine,
etc? In these cases one doesn't have a straightforward functional
description of the geometry, but just points that lie of the boundary or
a complex collection of spline.

>> I've been attracted to this approach ever since Matt first presented
>> the design for the Sieve where one of the main points is to representl
>> the coordinates of the mesh as a field, or in Sieve-speak a section of
>> a vector bundle... ;-)
>>
>> But this argument may be worth less these days as we can no longer
>> store fields (Functions) to file.
>>
>> --
>> Anders
> 
> Is this really different from having a finite element space for the
> mapping?  You can define a finite element function in that space to store
> the higher order mesh data.  Or maybe have a separate higher order
> geometry class that can do that OR some other representation like NURBS.
> 

This is what I would like - the flexibility to use a functional
representation (including functions other than FE functions, such as
NURBS) or an interpolated approach with higher-order vertices.

Garth

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