dolfin team mailing list archive

[walker@xxxxxxxxxxxx]

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.

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

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