dolfin team mailing list archive

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

On 25/03/11 12:58, Anders Logg wrote:
> On Fri, Mar 25, 2011 at 10:28:20AM +0000, Garth N. Wells wrote:
>>
>>
>> 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.
> 
> We would then need a function space that can represent those kinds of
> fields (NURBS). 

Great. We all agree on this.

> If we add such functionality (to DOLFIN and FFC), we
> could also use the NURBS as basis functions ("isogeometric analysis")
> which seems to be popular these days.
>

There are different ways to go from some geometric representation to a
function that represents the mapping on each cell, and the choice will
often be *forced* on the user by the format of the geometry data that is
provided/available. This is why need maximum flexibility. It's not
feasible to support in FFC all possible maps. That's why I suggest
(which I think is along the lines of what Shawn suggests) that the map
be associated with a ufc::cell, with a member function that can return
the Jacobian. For simple cases (e.g. isoparmetric elements using FE
basis functions), FFC could generate the necessary ufc::cell code. The
input would be either a functional representation, or cells with
'higher-order' nodes, from which the map can be constructed. For
arbitrary maps, a user can overlaod the ufc::cell::jacobian function.
When the geometry features of DOLFIN are improved, DOLFIN may eventually
support some maps, e.g. for certain splines.

Associating a map with a cell provides a clean separation between the
geometry of a domain and a variational formulation.

Garth


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