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