From: Vijay S. M. <vi...@gm...> - 2009-02-24 19:10:54
|
Hi, Are there any higher than second order Finite Elements in Libmesh ? I obviously do not see them in the class docs but just wondering if there was some trick to apply say a 5th order Lagrange basis to a QUAD4 elem. Then would Libmesh automatically create the extra dofs needed to make this unisolvent ?! If this is gibberish and libmesh does not have higher than second-order elements, do let me know. If you have suggestions to implement higher order elements, I would be glad to hear them also. Thanks, Vijay |
From: Roy S. <roy...@ic...> - 2009-02-24 19:31:09
|
On Tue, 24 Feb 2009, Vijay S. Mahadevan wrote: > Are there any higher than second order Finite Elements in Libmesh ? I > obviously do not see them in the class docs but just wondering if > there was some trick to apply say a 5th order Lagrange basis to a > QUAD4 elem. Then would Libmesh automatically create the extra dofs > needed to make this unisolvent ?! There are higher order elements, just not higher order Lagrange elements; try the HIERARCHIC basis instead if you want C0 higher order polynomials. You can even do adaptive p refinement with them. (probably hp, too, but I wouldn't guarantee that's working) There is a limitation on geometric element compatibility, though: you need to have a node on every face or edge which corresponds to the center of support of a basis function. So to use your example, a QUAD4 won't work for most higher order elements (cubic HERMITEs being the only exception) because there's no nodes to store edge degrees of freedom. There's also a limitation on output: our .xda/.xdr formats will save all your higher order data to machine precision, but at least the visualization formats I use are first-order-only; I'm not sure we can even plot quadratics except by turning them into 2^d linears first. > If this is gibberish and libmesh does not have higher than > second-order elements, do let me know. In 2D, nth order HIERARCHICs should work on every second order geometric element. In 3D, they work on hexes only at the moment. > If you have suggestions to implement higher order elements, I would > be glad to hear them also. Adding support for HIERARCHIC tets and pyramids would first require adding the geometric elements to support them: a PYRAMID5 doesn't have the edge or face nodes necessary to support even second order elements, and a TET10 doesn't have the face nodes. We need to round things out with a PYRAMID14, PYRAMID19, and TET14 one of these days. If the impression you're getting is that we've each written exactly enough higher order stuff for our own apps to work, you'd be right. None of the missing bits I've mentioned would be hard to add to libMesh if you need them; they're just tedious enough that nobody's done so yet. --- Roy |
From: Vijay S. M. <vi...@gm...> - 2009-02-24 20:22:10
|
> There are higher order elements, just not higher order Lagrange > elements; try the HIERARCHIC basis instead if you want C0 higher order > polynomials. You can even do adaptive p refinement with them. > (probably hp, too, but I wouldn't guarantee that's working) Thanks for the quick reply. I did think about Hierarchic basis but I'm currently working on a few diffusion problems and Lagrange basis seems to work quite well. I will try to switch my implementation to use Hierarchics instead and see if there is an increase/decrease in performance/accuracy. Thanks for the suggestion. I am currently looking for adaptive p or h but hp would be future work. So I do not care about hp adaptivity for now. And hopefully, I will be able to figure out the visualization bugs for higher order solution with VTK since that is the format I use most frequently. > In 2D, nth order HIERARCHICs should work on every second order > geometric element. In 3D, they work on hexes only at the moment. Well, that should be enough for me now to get started. If I add something, I'll send a patch. > If the impression you're getting is that we've each written exactly > enough higher order stuff for our own apps to work, you'd be right. I understand that all the developers are doing your PhDs and it is more important to finish your research before adding all the features you can imagine, into a library, that selected few people use everyday. I am using LibMesh for my research too and just wanted to get an idea on what I would need to add/change, if I need higher than second order elements. Thanks for all the suggestions thus far Roy. I will email you if I have more questions regarding this topic. Cheers, Vijay On Tue, Feb 24, 2009 at 1:31 PM, Roy Stogner <roy...@ic...> wrote: > > On Tue, 24 Feb 2009, Vijay S. Mahadevan wrote: > >> Are there any higher than second order Finite Elements in Libmesh ? I >> obviously do not see them in the class docs but just wondering if >> there was some trick to apply say a 5th order Lagrange basis to a >> QUAD4 elem. Then would Libmesh automatically create the extra dofs >> needed to make this unisolvent ?! > > There are higher order elements, just not higher order Lagrange > elements; try the HIERARCHIC basis instead if you want C0 higher order > polynomials. You can even do adaptive p refinement with them. > (probably hp, too, but I wouldn't guarantee that's working) > > There is a limitation on geometric element compatibility, though: you > need to have a node on every face or edge which corresponds to the > center of support of a basis function. So to use your example, a > QUAD4 won't work for most higher order elements (cubic HERMITEs being > the only exception) because there's no nodes to store edge degrees of > freedom. > > There's also a limitation on output: our .xda/.xdr formats will save > all your higher order data to machine precision, but at least the > visualization formats I use are first-order-only; I'm not sure we can > even plot quadratics except by turning them into 2^d linears first. > >> If this is gibberish and libmesh does not have higher than >> second-order elements, do let me know. > > In 2D, nth order HIERARCHICs should work on every second order > geometric element. In 3D, they work on hexes only at the moment. > >> If you have suggestions to implement higher order elements, I would >> be glad to hear them also. > > Adding support for HIERARCHIC tets and pyramids would first require > adding the geometric elements to support them: a PYRAMID5 doesn't have > the edge or face nodes necessary to support even second order > elements, and a TET10 doesn't have the face nodes. We need to round > things out with a PYRAMID14, PYRAMID19, and TET14 one of these days. > > If the impression you're getting is that we've each written exactly > enough higher order stuff for our own apps to work, you'd be right. > None of the missing bits I've mentioned would be hard to add to > libMesh if you need them; they're just tedious enough that nobody's > done so yet. > --- > Roy > |
From: Roy S. <roy...@ic...> - 2009-02-24 20:26:52
|
On Tue, 24 Feb 2009, Vijay S. Mahadevan wrote: > I understand that all the developers are doing your PhDs It's gotten worse than that now that the primary developers are all graduated. E.g. the adaptive p refinement was something I did on a lark as a break from my dissertation research. It's harder to find time for things like that now that most of my work isn't libMesh related and the libMesh stuff needs to be either reasonably related to my job or on my own time. --- Roy |
From: Ondrej C. <on...@ce...> - 2009-03-01 11:47:39
|
Hi Roy! On Tue, Feb 24, 2009 at 12:26 PM, Roy Stogner <roy...@ic...> wrote: > > On Tue, 24 Feb 2009, Vijay S. Mahadevan wrote: > >> I understand that all the developers are doing your PhDs > > It's gotten worse than that now that the primary developers are all > graduated. E.g. the adaptive p refinement was something I did on a > lark as a break from my dissertation research. It's harder to find > time for things like that now that most of my work isn't libMesh > related and the libMesh stuff needs to be either reasonably related to > my job or on my own time. Just curious --- what do you work on now? Congratulations to your graduation! Ondrej |
From: Kirk, B. (JSC-EG311) <ben...@na...> - 2009-03-01 20:30:46
|
>>> I understand that all the developers are doing your PhDs >> >> It's gotten worse than that now that the primary developers are all >> graduated. E.g. the adaptive p refinement was something I did on a >> lark as a break from my dissertation research. It's harder to find >> time for things like that now that most of my work isn't libMesh >> related and the libMesh stuff needs to be either reasonably related to >> my job or on my own time. It's not all doom and gloom, though, a number of us are using libMesh and derived application codes to solve problems at work... So at least I am confident there is a place for the library well into the future. The issue is that these tend to be more "industrial-type" applications, where higher-order elements are often not used for various reasons (non-smooth solutions, sharp complex geometry, etc...) So what we really need is someone who wants to become a primary developer and then jumps to a pure applied math position. (Does Varis read this thread? ;-) ) -Ben |
From: Roy S. <roy...@ic...> - 2009-03-01 16:31:50
|
On Sun, 1 Mar 2009, Ondrej Certik wrote: > On Tue, Feb 24, 2009 at 12:26 PM, Roy Stogner <roy...@ic...> wrote: >> >> On Tue, 24 Feb 2009, Vijay S. Mahadevan wrote: >> >>> I understand that all the developers are doing your PhDs >> >> It's gotten worse than that now that the primary developers are all >> graduated. E.g. the adaptive p refinement was something I did on a >> lark as a break from my dissertation research. It's harder to find >> time for things like that now that most of my work isn't libMesh >> related and the libMesh stuff needs to be either reasonably related to >> my job or on my own time. > > Just curious --- what do you work on now? A postdoctoral fellowship here: http://www.ices.utexas.edu/centers/pecos/ I'm working on hypersonic flow code - we're stuck using an old NASA Fortran code to meet some early deliverables, but it's suboptimal for our coupling and UQ needs, and we'll probably be using (an expanded version of) some libMesh hypersonics code of Ben's for the long term work. --- Roy |
From: Ondrej C. <on...@ce...> - 2009-03-02 15:24:43
|
On Sun, Mar 1, 2009 at 11:31 AM, Roy Stogner <roy...@ic...> wrote: > > On Sun, 1 Mar 2009, Ondrej Certik wrote: > >> On Tue, Feb 24, 2009 at 12:26 PM, Roy Stogner <roy...@ic...> >> wrote: >>> >>> On Tue, 24 Feb 2009, Vijay S. Mahadevan wrote: >>> >>>> I understand that all the developers are doing your PhDs >>> >>> It's gotten worse than that now that the primary developers are all >>> graduated. E.g. the adaptive p refinement was something I did on a >>> lark as a break from my dissertation research. It's harder to find >>> time for things like that now that most of my work isn't libMesh >>> related and the libMesh stuff needs to be either reasonably related to >>> my job or on my own time. >> >> Just curious --- what do you work on now? > > A postdoctoral fellowship here: > > http://www.ices.utexas.edu/centers/pecos/ > > I'm working on hypersonic flow code - we're stuck using an old NASA > Fortran code to meet some early deliverables, but it's suboptimal for > our coupling and UQ needs, and we'll probably be using (an expanded > version of) some libMesh hypersonics code of Ben's for the long term > work. I see, thanks for the update. I myself moved to Nevada/Reno to this group: http://hpfem.math.unr.edu/ On Sun, Mar 1, 2009 at 3:30 PM, Kirk, Benjamin (JSC-EG311) <ben...@na...> wrote: >>>> I understand that all the developers are doing your PhDs >>> >>> It's gotten worse than that now that the primary developers are all >>> graduated. E.g. the adaptive p refinement was something I did on a >>> lark as a break from my dissertation research. It's harder to find >>> time for things like that now that most of my work isn't libMesh >>> related and the libMesh stuff needs to be either reasonably related to >>> my job or on my own time. > > It's not all doom and gloom, though, a number of us are using libMesh and > derived application codes to solve problems at work... So at least I am > confident there is a place for the library well into the future. The issue > is that these tend to be more "industrial-type" applications, where > higher-order elements are often not used for various reasons (non-smooth > solutions, sharp complex geometry, etc...) Well, we do higher order FEM and we *are* interested in industrial type applications, e.g. non-smooth solutions, sharp complex geometry etc. But we just put our codes to the web and so far only the 2D code is production ready, but we are working very hard on the 3D code too and we are close. Our codes are here: http://hpfem.org/ If I find time, I'd like to create Python interface to libmesh, because I used it in the past and I would like to compare libmesh h-adaptivity with our hp-adaptivity. But I would like to build libmesh without any external solvers, as I would like to use the same solver for both our code and libmesh, called from Python. I did some work on that in the past, I'll have a look at the newest libmesh, if it's possible to build it like that. Ondrej |
From: Roy S. <roy...@ic...> - 2009-03-02 17:16:08
|
On Mon, 2 Mar 2009, Ondrej Certik wrote: > If I find time, I'd like to create Python interface to libmesh, > because I used it in the past and I would like to compare libmesh > h-adaptivity with our hp-adaptivity. If you're already doing that, you might want to go one step farther and compare libmesh hp-adaptivity with yours. The former is theoretically working but I've rarely gotten the convergence rates I expected out of it, so I suspect there's some subtle bugs left. --- Roy |
From: Ondrej C. <on...@ce...> - 2009-03-02 22:09:37
|
On Mon, Mar 2, 2009 at 12:16 PM, Roy Stogner <roy...@ic...> wrote: > > > On Mon, 2 Mar 2009, Ondrej Certik wrote: > >> If I find time, I'd like to create Python interface to libmesh, >> because I used it in the past and I would like to compare libmesh >> h-adaptivity with our hp-adaptivity. > > If you're already doing that, you might want to go one step farther > and compare libmesh hp-adaptivity with yours. The former is > theoretically working but I've rarely gotten the convergence rates I > expected out of it, so I suspect there's some subtle bugs left. Ah, yes, absolutely! I didn't know that libmesh also has hp-adaptivity. Then I will definitely give it a shot. Ondrej |
From: Roy S. <roy...@ic...> - 2009-03-02 22:36:07
|
On Mon, 2 Mar 2009, Ondrej Certik wrote: > On Mon, Mar 2, 2009 at 12:16 PM, Roy Stogner <roy...@ic...> wrote: >> >> >> On Mon, 2 Mar 2009, Ondrej Certik wrote: >> >>> If I find time, I'd like to create Python interface to libmesh, >>> because I used it in the past and I would like to compare libmesh >>> h-adaptivity with our hp-adaptivity. >> >> If you're already doing that, you might want to go one step farther >> and compare libmesh hp-adaptivity with yours. The former is >> theoretically working but I've rarely gotten the convergence rates I >> expected out of it, so I suspect there's some subtle bugs left. > > Ah, yes, absolutely! I didn't know that libmesh also has > hp-adaptivity. Well, that might be overstating the case. ;-) libMesh has the ability to adaptively refine p-hierarchic element types in h and/or p, which technically might qualify as "has (isotropic) hp-adaptivity"... but libMesh lacks a good heuristic for deciding which type of adaptation to do. IIRC the only strategy classes I wrote were "h adapt if you touch a prespecified singularity, p adapt otherwise" and "try to see whether an h-coarsening or a p-coarsening would be worse, then h-adapt or p-adapt respectively instead". The former seemed to give exponential but probably non-optimal convergence in a couple benchmark problems; the latter didn't even go exponential. > Then I will definitely give it a shot. Thanks! I would have liked to spend more time on this myself, but eventually dissertation work took precedence. For that matter, I'd love to have an excuse to play with it now, but I still haven't figured out how to make Ben's hypersonics formulations work perfectly with hierarchic elements. --- Roy |
From: Ondrej C. <on...@ce...> - 2009-03-02 22:50:18
|
Hi, I CCed to our hpfem list, so that people in our group can have a look at this too. On Mon, Mar 2, 2009 at 5:35 PM, Roy Stogner <roy...@ic...> wrote: > > > On Mon, 2 Mar 2009, Ondrej Certik wrote: > >> On Mon, Mar 2, 2009 at 12:16 PM, Roy Stogner <roy...@ic...> >> wrote: >>> >>> >>> On Mon, 2 Mar 2009, Ondrej Certik wrote: >>> >>>> If I find time, I'd like to create Python interface to libmesh, >>>> because I used it in the past and I would like to compare libmesh >>>> h-adaptivity with our hp-adaptivity. >>> >>> If you're already doing that, you might want to go one step farther >>> and compare libmesh hp-adaptivity with yours. The former is >>> theoretically working but I've rarely gotten the convergence rates I >>> expected out of it, so I suspect there's some subtle bugs left. >> >> Ah, yes, absolutely! I didn't know that libmesh also has >> hp-adaptivity. > > Well, that might be overstating the case. ;-) libMesh has the > ability to adaptively refine p-hierarchic element types in h and/or p, > which technically might qualify as "has (isotropic) hp-adaptivity"... > but libMesh lacks a good heuristic for deciding which type of > adaptation to do. IIRC the only strategy classes I wrote were "h > adapt if you touch a prespecified singularity, p adapt otherwise" and > "try to see whether an h-coarsening or a p-coarsening would be worse, > then h-adapt or p-adapt respectively instead". The former seemed to > give exponential but probably non-optimal convergence in a couple > benchmark problems; the latter didn't even go exponential. > >> Then I will definitely give it a shot. > > Thanks! I would have liked to spend more time on this myself, but > eventually dissertation work took precedence. > > For that matter, I'd love to have an excuse to play with it now, but I > still haven't figured out how to make Ben's hypersonics formulations > work perfectly with hierarchic elements. That's awesome that you are interested in this. At this moment, I am quite busy with other things, but I am very interested in this comparison. Does the isotropic hp-adaptivity work in 3d? The difficult part is in choosing the right candidate for hp-adaptivity, and I know that the guys in my group spent months in getting it right, so I am curious in comparing it with your version to see if it's different. I want to see the graph, then I will believe. :) On Sun, Mar 1, 2009 at 3:30 PM, Kirk, Benjamin (JSC-EG311) <ben...@na...> wrote: > is that these tend to be more "industrial-type" applications, where > higher-order elements are often not used for various reasons (non-smooth > solutions, sharp complex geometry, etc...) In fact, hp-fem performs the best exactly with solutions that are both non-smooth and sharp somewhere (it uses a low polynomial order there) and very smooth somewhere else in the domain (it uses a high polynomial order there). But as I said, it's tough to get everything into a production ready state, but it's exciting when we get there finally both in 2D and in 3D. Ondrej |
From: Roy S. <roy...@ic...> - 2009-03-02 23:28:32
|
On Mon, 2 Mar 2009, Ondrej Certik wrote: > That's awesome that you are interested in this. At this moment, I am > quite busy with other things, but I am very interested in this > comparison. Does the isotropic hp-adaptivity work in 3d? It works for any elements we've got hierarchic basis functions available on, but right now that means just hexes in 3D. > <ben...@na...> wrote: >> is that these tend to be more "industrial-type" applications, where >> higher-order elements are often not used for various reasons (non-smooth >> solutions, sharp complex geometry, etc...) > > In fact, hp-fem performs the best exactly with solutions that are both > non-smooth and sharp somewhere (it uses a low polynomial order there) > and very smooth somewhere else in the domain (it uses a high > polynomial order there). Yeah; locally non-smooth solutions should be fine for hp elements. Complex geometry is a little trickier, unless you've got mechanisms for curved element sides (libMesh does, but only quadratics) and proper refinement into exactly specified geometry (libMesh doesn't unless user code adds it). But I think Ben left out the biggest reasons why industrial applications don't get solved with hp: formulation fragility and code complexity. Some numerical formulations work fine at low order and require significant analysis and redesign before achieving optimal convergence (or even converging at all) at high order. Most industrial software is complex, and rearchitecting to get faster and more accurate solves is deemed (rightly or wrongly) less productive than making a GUI nicer. --- Roy |
From: Kirk, B. (JSC-EG311) <ben...@na...> - 2009-03-03 02:39:14
|
> <ben...@na...> wrote: >> is that these tend to be more "industrial-type" applications, where >> higher-order elements are often not used for various reasons (non-smooth >> solutions, sharp complex geometry, etc...) > > In fact, hp-fem performs the best exactly with solutions that are both > non-smooth and sharp somewhere (it uses a low polynomial order there) > and very smooth somewhere else in the domain (it uses a high > polynomial order there). > > But as I said, it's tough to get everything into a production ready > state, but it's exciting when we get there finally both in 2D and in > 3D. my myopic view is tainted because my applications are nearly hyperbolic and everything of interest is downstream of a shockwave. so even for flows over blunt bodies where the shock layer is smooth, at some point upstream pollution becomes dominant. it is of course theoretically possible to use a low order element and allow h-adaptivity to resolve the shock, but you need ridiculously small elements to mask the pollution. of course, I'd love a demonstrated counter example... -Ben |