From: John P. <jwp...@gm...> - 2008-11-06 16:33:19
|
Hi all, I just fixed up the quadrature rules for pyramids, which have been wrong since I first coded them up years ago!! If you've ever done anything with pyramids in the library, it was probably at the very least inaccurate and (most likely) totally wrong :) Anyway, I actually wanted to start a thread about (the lack of) pyramid refinement in the library. I've been working with some viz folks here at TACC recently, and after telling them "yeah we can do pyramids" I actually went back and checked... Anyway, there appears to be at least a couple issues to discuss: 1.) A possible refinement pattern for the pyramid involves 10 children: 4 sub-pyramids on the base, one in the apex, and one pyramid upside-down below the apex pyramid. Then, there are 4 tets filling in the gaps. (This is a little hard to visualize, but I can send you an image if you want to see what I mean.) One bad thing about this refinement pattern is that under repeated refinements, you end up with a little "ribbon" of pyramids running through the original element. So, we probably need to discuss other possibilities, but this one would at least get us started. 2.) *Any* refinement pattern for the pyramid will likely involve a splitting where some children have a different type than the parent. Since this would be the only element with such an inhomogeneous splitting, I'd like to discuss the best way to do this in the present context of the embedding matrix. Since embedding_matrix() is a virtual function, it can be redefined in Pyramid sub-classes to do whatever we want, e.g. we could place -1 entries in the actual tables and just ensure, through embedding_matrix(), that those entries are never accessed. The real issue seems to be determining where in the library, if anywhere, we've assumed that all the children have the same geometric type. I'd be glad to hear your comments there. 3.) The pyramid basis functions are a little weird (rational), and, for example, their evaluation at the apex is a little touchy. I need to do more testing on these, including checking the accuracy of integration for e.g. a mass matrix, and making sure the element Jacobian is well-defined everywhere. Anyone know anything about the accuracy of quadrature for functions which are ratios of polynomials? I asked Roy about some sort of macroelement splitting to get away from the rational basis functions, but I don't think that idea can really go anywhere. If you have read anything about other types of pyramid finite elements besides the one we've got, please let me know. Finally, to be more useful in incompressible Navier-Stokes calculations, were are gonna need a quadratic pyramid at some point. I have a paper on how to define arbitrary-order Lagrange pyramids with rational basis functions, so it's on my long term development plan. Your comments on any of the above are welcomed. -- John |
From: Benjamin K. <ben...@na...> - 2008-11-06 17:17:07
|
> 2.) *Any* refinement pattern for the pyramid will likely involve a > splitting where some children have a different type than the parent. > Since this would be the only element with such an inhomogeneous > splitting, I'd like to discuss the best way to do this in the present > context of the embedding matrix. Since embedding_matrix() is a > virtual function, it can be redefined in Pyramid sub-classes to do > whatever we want, e.g. we could place -1 entries in the actual tables > and just ensure, through embedding_matrix(), that those entries are > never accessed. The real issue seems to be determining where in the > library, if anywhere, we've assumed that all the children have the > same geometric type. I'd be glad to hear your comments there. I'll have to look way back in the logs. Before the embedding_matrix we actually did pyramids with the first split you proposed. I'd rather go to the second (only sub-pyramids on the base) for a reimplementation. One issue which has to be addressed is the recursive sideset stuff. That is why the pyramid has a 'nonstandard' side numbering: 0,1,2,3 are the triangular sides, and 4 is the square base. It has to be that way so the child tets (which only have 4 sides) can inherit bcs from the same side in the parent. So that is not a problem, but an issue that must be addressed. The only other thing that comes to mind is places like constraint calculation where we compute on a parent and child. We just need to be careful there is no assumption that they are the same type... I don't think this is an issue at all for Lagrange elements, because constraints are only calculated using the sides, which are of the same type. But perhaps some of Roy's more exotic FE spaces? -Ben |
From: Roy S. <roy...@ic...> - 2008-11-06 17:26:56
|
Benjamin Kirk wrote: > But perhaps some of > Roy's more exotic FE spaces? For now, you can forget about seeing very exotic FE spaces on our pyramids and tets: topological DoF storage only works if you've got a DofObject for every necessary topological connection, and we'd need elements with a midnode on every triangular face for that. So even after John's added his desired Pyramid14, we'll still need a Pyramid18 and a Tet14 before we can think about high p elements on them. --- Roy |
From: John P. <jwp...@gm...> - 2008-11-06 17:18:08
|
On Thu, Nov 6, 2008 at 11:03 AM, Roy Stogner <roy...@ic...> wrote: > John Peterson wrote: > >> 1.) A possible refinement pattern for the pyramid involves 10 >> children: 4 sub-pyramids on the base, one in the apex, and one pyramid >> upside-down below the apex pyramid. Then, there are 4 tets filling in >> the gaps. (This is a little hard to visualize, but I can send you an >> image if you want to see what I mean.) > > Sadly, this is probably the *easiest* alternative to visualize. > >> One bad thing about this >> refinement pattern is that under repeated refinements, you end up with >> a little "ribbon" of pyramids running through the original element. >> So, we probably need to discuss other possibilities, but this one >> would at least get us started. > > Replace those two central pyramids with four tets? Should be just as easy > and arguably just as good, especially if we do the same sort of > geometrically adaptive axis selection that we do in tet refinement now. Hey, that's a pretty good idea! I will move that to the top of the list of candidates. >> Since embedding_matrix() is a >> virtual function, it can be redefined in Pyramid sub-classes to do >> whatever we want, e.g. we could place -1 entries in the actual tables >> and just ensure, through embedding_matrix(), that those entries are >> never accessed. The real issue seems to be determining where in the >> library, if anywhere, we've assumed that all the children have the >> same geometric type. I'd be glad to hear your comments there. > > There is an assumption in FEBase::coarsened_dof_values that elements and > their children are oriented the same way - i.e. if child c has a side on > side s of its parent, then the corresponding side number on c is also s. > That's actually an assumption that we don't *have* to break with either > potential pyramid splitting, but in either case we'd have to be careful with > the tet node ordering. OK, interesting. I wasn't aware of that particular function in FEBase so I'm glad you mentioned it. For non-homogeneous refinement patterns, this assumption might not even make sense. A pyramid has five sides but its tet children will only have 4, for example. It probably wouldn't be too hard to remove that assumption from FEBase but it would almost certainly be a performance hit to do so... >> Anyone know anything about the >> accuracy of quadrature for functions which are ratios of polynomials? > > We can derive custom quadrature rules which would integrate a mass matrix > exactly... but would they then also integrate, say, a Laplacian matrix > exactly? The answer is an obvious "yes" for polynomial bases but I'd expect > a "no" for pyramids. That could be bad. > > What are we doing for them now? The current quadrature rules have accuracies like you would expect for 1D elements, since they are conical products of Gauss-like rules. So, for example, a 2x2x2 rule will integrate exactly all monomials of the form x^a y^b z^c, a+b+c <= 3. I have no idea what will happen when we try to integrate the rational basis functions... -- John |
From: Roy S. <roy...@ic...> - 2008-11-06 17:29:01
|
John Peterson wrote: > 1.) A possible refinement pattern for the pyramid involves 10 > children: 4 sub-pyramids on the base, one in the apex, and one pyramid > upside-down below the apex pyramid. Then, there are 4 tets filling in > the gaps. (This is a little hard to visualize, but I can send you an > image if you want to see what I mean.) Sadly, this is probably the *easiest* alternative to visualize. > One bad thing about this > refinement pattern is that under repeated refinements, you end up with > a little "ribbon" of pyramids running through the original element. > So, we probably need to discuss other possibilities, but this one > would at least get us started. Replace those two central pyramids with four tets? Should be just as easy and arguably just as good, especially if we do the same sort of geometrically adaptive axis selection that we do in tet refinement now. > Since embedding_matrix() is a > virtual function, it can be redefined in Pyramid sub-classes to do > whatever we want, e.g. we could place -1 entries in the actual tables > and just ensure, through embedding_matrix(), that those entries are > never accessed. The real issue seems to be determining where in the > library, if anywhere, we've assumed that all the children have the > same geometric type. I'd be glad to hear your comments there. Surprisingly, I don't think this will be as pervasive a problem as you might assume. We already reinit separate FE objects for parents and children when necessary, and we abstract the FE/Elem interaction away from code using them sufficiently well that there was never any need to make constant elem type assumptions. There is an assumption in FEBase::coarsened_dof_values that elements and their children are oriented the same way - i.e. if child c has a side on side s of its parent, then the corresponding side number on c is also s. That's actually an assumption that we don't *have* to break with either potential pyramid splitting, but in either case we'd have to be careful with the tet node ordering. > Anyone know anything about the > accuracy of quadrature for functions which are ratios of polynomials? We can derive custom quadrature rules which would integrate a mass matrix exactly... but would they then also integrate, say, a Laplacian matrix exactly? The answer is an obvious "yes" for polynomial bases but I'd expect a "no" for pyramids. That could be bad. What are we doing for them now? > Finally, to be more useful in incompressible Navier-Stokes > calculations, were are gonna need a quadratic pyramid at some point. > I have a paper on how to define arbitrary-order Lagrange pyramids with > rational basis functions, so it's on my long term development plan. That would be excellent. --- Roy |
From: Roy S. <roy...@ic...> - 2008-11-06 17:35:21
|
John Peterson wrote: > On Thu, Nov 6, 2008 at 11:03 AM, Roy Stogner <roy...@ic...> wrote: >> Replace those two central pyramids with four tets? Should be just as easy >> and arguably just as good, especially if we do the same sort of >> geometrically adaptive axis selection that we do in tet refinement now. > > Hey, that's a pretty good idea! I will move that to the top of the > list of candidates. The adaptive axis selection can probably go on the back burner, though. It was critical for Tets since the effect of bad axis selection there is cumulative. For Pyramids, a bad selection would only happen once and wouldn't be exacerbated by the subsequent Tet refinement. > For non-homogeneous refinement > patterns, this assumption might not even make sense. A pyramid has > five sides but its tet children will only have 4, for example. Ben appears to have already figured out this (in another context): the final pyramid side is the bottom, which Tets will never overlap. --- Roy |
From: John P. <jwp...@gm...> - 2008-11-07 17:47:57
|
There's some more strange things going on with the pyramid... In the definition of the geometric element (src/geom/cell_pyramid5.C) the numbering of the base is *clockwise* instead of following the usual counter-clockwise numbering used by e.g. Hexes. The build_side() code in the class is self-consistent with this numbering, but the numbering of the Lagrange basis functions is not, nor is the element defined in reference_elements/3D/one_pyramid.xda (this reference element should be inverted with respect to the definition in the Pyramid5 class, as far as I can tell.) Some questions... Remind me again what was the reason for numbering it this way? I'm sure there was one, I just can't remember what it was. I'll change the basis function numbering and the reference element definition. Is there any other obvious part of the library to fix up? Thanks, -- John |
From: Benjamin K. <ben...@na...> - 2008-11-07 18:05:56
|
> Remind me again what was the reason for numbering it this way? I'm > sure there was one, I just can't remember what it was. I'm not sure there was a reason, or at least I can't find anything that jogs my memory looking back through the logs. One possible culprit, though - the ordering of our pyramid matches exactly the ordering used in vtk. Would we be better off reorienting the base? In that case the side creation would have to change, as well as any read functions which take in a pyramid. The only sucky part here is keeping backward compatiblity in the .xda/.xdr mesh format. -Ben |
From: John P. <jwp...@gm...> - 2008-11-07 18:10:16
|
On Fri, Nov 7, 2008 at 12:05 PM, Benjamin Kirk <ben...@na...> wrote: >> Remind me again what was the reason for numbering it this way? I'm >> sure there was one, I just can't remember what it was. > > I'm not sure there was a reason, or at least I can't find anything that jogs > my memory looking back through the logs. One possible culprit, though - the > ordering of our pyramid matches exactly the ordering used in vtk. > > Would we be better off reorienting the base? In that case the side creation > would have to change, as well as any read functions which take in a pyramid. > > The only sucky part here is keeping backward compatiblity in the .xda/.xdr > mesh format. Yeah, we could either renumber the shape functions or change the numbering of the reference element. I figured the former offered the path of least resistance, but I could be convinced to do either. -- John |
From: Benjamin K. <ben...@na...> - 2008-11-07 18:10:31
|
>> Remind me again what was the reason for numbering it this way? I'm >> sure there was one, I just can't remember what it was. > > I'm not sure there was a reason, or at least I can't find anything that jogs > my memory looking back through the logs. One possible culprit, though - the > ordering of our pyramid matches exactly the ordering used in vtk. ..looks like it matches the ExodusII ordering too. Perhaps we just used it because someone else did? |
From: John P. <jwp...@gm...> - 2008-11-13 17:30:44
|
On Fri, Nov 7, 2008 at 11:37 AM, John Peterson <jwp...@gm...> wrote: > There's some more strange things going on with the pyramid... > > In the definition of the geometric element (src/geom/cell_pyramid5.C) > the numbering of the base is *clockwise* instead of following the > OK, I was smoking crack when I first wrote this. There is no clockwise node numbering in the pyramid!! Momentarily, I will check in changes that allow build_cube() to take PYRAMID5 as the element type. Similarly to the tet case, we initially build a hex mesh and then split each hex into 6 pyramids. -- John |
From: John P. <jwp...@gm...> - 2008-11-12 22:25:12
|
On Thu, Nov 6, 2008 at 11:18 AM, John Peterson <jwp...@gm...> wrote: > On Thu, Nov 6, 2008 at 11:03 AM, Roy Stogner <roy...@ic...> wrote: >> John Peterson wrote: >>> Anyone know anything about the >>> accuracy of quadrature for functions which are ratios of polynomials? >> >> We can derive custom quadrature rules which would integrate a mass matrix >> exactly... but would they then also integrate, say, a Laplacian matrix >> exactly? The answer is an obvious "yes" for polynomial bases but I'd expect >> a "no" for pyramids. That could be bad. >> >> What are we doing for them now? > > The current quadrature rules have accuracies like you would expect for > 1D elements, since they are conical products of Gauss-like rules. So, > for example, a 2x2x2 rule will integrate exactly all monomials of the > form x^a y^b z^c, a+b+c <= 3. I have no idea what will happen when we > try to integrate the rational basis functions... Just a quick update on the quadrature over pyramids stuff. After checking it with Maple, it appears that the "standard" 2nd/3rd-order quadrature rule can exactly integrate the Pyramid5 mass matrix. The laplace matrix, however, is a different story. I needed to go up to 6/7th-order quadrature before I could get 9-10 digits of precision from LibMesh. At first, this seems a little paradoxical since the Laplace matrix is usually the easier of the two, but with rational basis functions, the more derivatives you take the more poles you get in the denominator, and the harder it is to integrate the functions. Since the default quadrature rule is currently selected by the FEType without regard to the geometric element type, it's not immediately obvious how we should ensure the user gets accurate quadrature on pyramids. A couple options... 1.) Just remember that higher-than-normal order quadrature on pyramids is required and your answer may be inaccurate. AKA "do nothing" :-) 2.) Redefine, within the pyramid quadrature rules, the meaning of order. I.e. return a rule several orders higher than what the user requests. 3.) Research quadrature rules for rational functions. I have a few papers on this but haven't looked into it too much yet. Phillipe Devloo may be doing something special in his library, so I will check there as well... -- John |
From: Derek G. <fri...@gm...> - 2008-11-12 22:33:39
|
For now... I would vote for "do nothing"... maybe print a warning in debug mode. Having the library try to interpret what you "really" want might be trouble. Derek On Wed, Nov 12, 2008 at 3:24 PM, John Peterson <jwp...@gm...> wrote: > On Thu, Nov 6, 2008 at 11:18 AM, John Peterson <jwp...@gm...> > wrote: > > On Thu, Nov 6, 2008 at 11:03 AM, Roy Stogner <roy...@ic...> > wrote: > >> John Peterson wrote: > >>> Anyone know anything about the > >>> accuracy of quadrature for functions which are ratios of polynomials? > >> > >> We can derive custom quadrature rules which would integrate a mass > matrix > >> exactly... but would they then also integrate, say, a Laplacian matrix > >> exactly? The answer is an obvious "yes" for polynomial bases but I'd > expect > >> a "no" for pyramids. That could be bad. > >> > >> What are we doing for them now? > > > > The current quadrature rules have accuracies like you would expect for > > 1D elements, since they are conical products of Gauss-like rules. So, > > for example, a 2x2x2 rule will integrate exactly all monomials of the > > form x^a y^b z^c, a+b+c <= 3. I have no idea what will happen when we > > try to integrate the rational basis functions... > > Just a quick update on the quadrature over pyramids stuff. > > After checking it with Maple, it appears that the "standard" > 2nd/3rd-order quadrature rule can exactly integrate the Pyramid5 mass > matrix. The laplace matrix, however, is a different story. I needed > to go up to 6/7th-order quadrature before I could get 9-10 digits of > precision from LibMesh. > > At first, this seems a little paradoxical since the Laplace matrix is > usually the easier of the two, but with rational basis functions, the > more derivatives you take the more poles you get in the denominator, > and the harder it is to integrate the functions. Since the default > quadrature rule is currently selected by the FEType without regard to > the geometric element type, it's not immediately obvious how we should > ensure the user gets accurate quadrature on pyramids. A couple > options... > > 1.) Just remember that higher-than-normal order quadrature on pyramids > is required and your answer may be inaccurate. AKA "do nothing" :-) > 2.) Redefine, within the pyramid quadrature rules, the meaning of > order. I.e. return a rule several orders higher than what the user > requests. > 3.) Research quadrature rules for rational functions. I have a few > papers on this but haven't looked into it too much yet. Phillipe > Devloo may be doing something special in his library, so I will check > there as well... > > -- > John > > ------------------------------------------------------------------------- > This SF.Net email is sponsored by the Moblin Your Move Developer's > challenge > Build the coolest Linux based applications with Moblin SDK & win great > prizes > Grand prize is a trip for two to an Open Source event anywhere in the world > http://moblin-contest.org/redirect.php?banner_id=100&url=/ > _______________________________________________ > Libmesh-devel mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-devel > |
From: Roy S. <roy...@ic...> - 2008-11-12 22:39:46
|
On Wed, 12 Nov 2008, John Peterson wrote: > 1.) Just remember that higher-than-normal order quadrature on pyramids > is required and your answer may be inaccurate. AKA "do nothing" :-) Not a good answer, but probably a fair one. > 2.) Redefine, within the pyramid quadrature rules, the meaning of > order. I.e. return a rule several orders higher than what the user > requests. Maybe an "extra_pyramid_quadrature_order()" option? So we can give it a safe default but still let users override it for efficiency if they know what they're doing. > 3.) Research quadrature rules for rational functions. I have a few > papers on this but haven't looked into it too much yet. Phillipe > Devloo may be doing something special in his library, so I will check > there as well... This could be ideal. --- Roy |
From: John P. <jwp...@gm...> - 2008-11-12 22:45:40
|
On Wed, Nov 12, 2008 at 4:38 PM, Roy Stogner <roy...@ic...> wrote: > > On Wed, 12 Nov 2008, John Peterson wrote: > >> 1.) Just remember that higher-than-normal order quadrature on pyramids >> is required and your answer may be inaccurate. AKA "do nothing" :-) > > Not a good answer, but probably a fair one. > >> 2.) Redefine, within the pyramid quadrature rules, the meaning of >> order. I.e. return a rule several orders higher than what the user >> requests. > > Maybe an "extra_pyramid_quadrature_order()" option? So we can give it > a safe default but still let users override it for efficiency if they > know what they're doing. Yeah, this sounds like a good idea, at least until something useful comes out of option 3... Maybe a flag in QBase, similar to the way we let knowledgable users disallow the use of rules with negative weights? >> 3.) Research quadrature rules for rational functions. I have a few >> papers on this but haven't looked into it too much yet. Phillipe >> Devloo may be doing something special in his library, so I will check >> there as well... > > This could be ideal. -- John |