dolfin team mailing list archive

[walker@xxxxxxxxxxxx]

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

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