On 24/03/11 19:10, [email protected] 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, [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. >>
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. >> 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. >> >> -- >> 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. > 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 <[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
