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, [email protected] 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 <[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. >> >>> >> >>> 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 <[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
