On 24/03/11 12:50, Kristian Ølgaard wrote: > On 24 March 2011 13:09, Garth N. Wells <[email protected]> wrote: >> >> >> On 24/03/11 04:38, [email protected] 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. > > I think there is already something in UFL that we could use for this. > The default (affine) triangle is: > > R2 = Space(2) > triangle = Cell("triangle", 1, R2) > > Changing the degree to '2': > > triangle = Cell("triangle", 2, R2) > > will imply a higher order geometry, of course we need to implement > support in the form compilers to generate actual code. >
But what kind of mapping does it imply? We want to generalised beyond using the FE basis functions for the map. Garth > Kristian > >> 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 <[email protected] >>>>>> <mailto:[email protected]>> 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 : [email protected] >>>>>> <mailto:[email protected]> >>>>>> > > 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 : [email protected] >>>>>> <mailto:[email protected]> >>>>>> > Unsubscribe : https://launchpad.net/~dolfin >>>>>> <https://launchpad.net/%7Edolfin> >>>>>> > More help : https://help.launchpad.net/ListHelp >>>>>> >>>>>> >>>> >>>> _______________________________________________ >>>> Mailing list: https://launchpad.net/~dolfin >>>> Post to : [email protected] >>>> Unsubscribe : https://launchpad.net/~dolfin >>>> More help : https://help.launchpad.net/ListHelp >>>> >>>> >>> >>> >>> >>> >>> >>> >> >> _______________________________________________ >> Mailing list: https://launchpad.net/~dolfin >> Post to : [email protected] >> Unsubscribe : https://launchpad.net/~dolfin >> More help : https://help.launchpad.net/ListHelp >> _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
