From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-01-19 15:25:22
|
Excellent! Will this be the default behavior and then let specific finite elements (e.g. LAGRANGE) overload or re-implement the method? The L2 projection method is perfect since it addresses the problem for arbitrary element types. In fact, it can be extended to use the normal equations to handle p mismatch on an edge as well. In general (see ~benkirk/Documents/lab_presentations/amr/slides.pdf in the CFDLab) A u_c = B u_p u_c = (A^-1 B) u_p u_c = C u_p There is no assumption about corner DOFs in libMesh, the implementation in fe_lagrange.C just recognizes that for Lagrange-type elements A=I, so u_c = C u_p directly. Line 622 in fe_lagrange.C simply says keep constraints that look like this: u_1 = 0.5*u_2 + 0.5*u_3 and ignore ones that look like this: u_1 = 0.9999 * u_1 + 1.e-6 * u_2 + 1.e-6 * u_3 Which we can recognize as u_1 = u_1, which is not a constraint at all... So, for the general case the constraint matrix C = (A^1 B) will have some rows that are like the identity matrix (within floating-point error, of course) which should be ignored and other rows that are valid DOF constraints. Clear? -Ben -----Original Message----- From: lib...@li... [mailto:lib...@li...] On Behalf Of Roy Stogner Sent: Monday, January 17, 2005 9:30 PM To: lib...@li... Subject: [Libmesh-devel] compute_constraints I'm writing a "fallback" implementation of compute_constraints; without knowing the details of every single degree of freedom it should be possible to use local L2 projections (of face values for C0 elements like the Hierarchic, face values and fluxes for C1 elements like the Clough-Tocher) to generate correct constraint matrices. One concern I have: is libMesh assuming that corner degrees of freedom and edge/face degrees of freedom should be equal when one element's corner node is a less refined element's edge or face node? If not, then I must be misunderstanding how the Lagrange compute_constraints is working, since it throws away identity values. If so, then that's going to be a problem for non-Lagrange elements whose edge and corner degrees of freedom are essentially unrelated. --- Roy Stogner ------------------------------------------------------- The SF.Net email is sponsored by: Beat the post-holiday blues Get a FREE limited edition SourceForge.net t-shirt from ThinkGeek. It's fun and FREE -- well, almost....http://www.thinkgeek.com/sfshirt _______________________________________________ Libmesh-devel mailing list Lib...@li... https://lists.sourceforge.net/lists/listinfo/libmesh-devel |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-01-19 17:53:33
|
Right, the local numbering does not match, but the global numbering should (which is all that really matters). I just added an assertion to, uh, assert this. I would continue with the implementation as you have it planned... I would like to keep the implementation in terms of parent-child relationships wherever possible to protect locality. We can handle the complications with recursive constraints later. As for the FE template madness, I don't have a good answer... :-) Just something that seemed like a good idea at the time. Like anything else, though, it could be changed. There are very few (none?) things in the library I consider set in stone... -Ben -----Original Message----- From: Roy Stogner [mailto:roy...@ic...] Sent: Wednesday, January 19, 2005 11:11 AM To: KIRK, BENJAMIN (JSC-EG) (NASA) Cc: lib...@li... Subject: RE: [Libmesh-devel] compute_constraints On Wed, 19 Jan 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Excellent! Will this be the default behavior and then let specific > finite elements (e.g. LAGRANGE) overload or re-implement the method? At the moment I'm making this a Clough-specific behavior, but making it the default after it's tested is my plan - reimplementation will just require adding more cases to the switch statements in fe_interface.C too. Remind me again why the finite element types are using template specialization instead of inheritance? ;-) > The L2 projection method is perfect since it addresses the problem for > arbitrary element types. In fact, it can be extended to use the > normal equations to handle p mismatch on an edge as well. I'd certainly like to see continuous hp-refinable function spaces in libMesh, but I'm not sure how much of a help this will be. Shouldn't the hierarchic basis functions be able to handle p mismatch for C0 elements? If element 1 has a polynomial degree p and adjacent element 2 has higher degree q, then isn't the constraint just "element 2 shape functions with degree p+1 through q and support on edge 1-2 get dropped"? No, on second thought it's more complicated than that - basis functions can have higher polynomial degree perpendicular to a face than they do in the face, and the case of simultaneous h refinement gets even uglier. Should I be preparing for this case in the automatic code, then? Right now I'm just doing the projection between parent and child degrees of freedom; perhaps I should be doing neighbor and child instead. > directly. Line 622 in fe_lagrange.C simply says keep constraints that > look like this: > > u_1 = 0.5*u_2 + 0.5*u_3 > > and ignore ones that look like this: > > u_1 = 0.9999 * u_1 + 1.e-6 * u_2 + 1.e-6 * u_3 > > Which we can recognize as u_1 = u_1, which is not a constraint at > all... What about lines that look like this? u_1 = 0.99999 * u_3 That's not on the identity diagonal, yet it's the sort of situation I would expect to happen when a nodal degree of freedom on a Lagrangian quadratic child element corresponds to a midedge or midface DoF on the parent: the coefficient is still 1, but the local numbering doesn't match. --- Roy Stogner |
From: Roy S. <roy...@ic...> - 2005-01-19 18:19:38
|
On Wed, 19 Jan 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > > What about lines that look like this? > > > > u_1 = 0.99999 * u_3 > > > > That's not on the identity diagonal, yet it's the sort of > > situation I would expect to happen when a nodal degree of freedom > > on a Lagrangian quadratic child element corresponds to a midedge > > or midface DoF on the parent: the coefficient is still 1, but the > > local numbering doesn't match. > Right, the local numbering does not match, but the global numbering > should (which is all that really matters). I just added an > assertion to, uh, assert this. Should the global numbering match? That's what's got me worried with the non-Lagrangian elements. Take the Clough-Tocher triangles for example: they've got three DoFs at each node (function value and gradient components) and one DoF at each side (normal flux, which generally won't correspond exactly to any corner DoF). So what happens in the global numbering when a node is both a side node for a less refined element and a corner node for more refined elements? Does libMesh assume that the edge DoF and the first of the node DoFs are equal valued? Does it give me all four global DoFs there and leave it up to me to constrain away three? > I would continue with the implementation as you have it planned... > I would like to keep the implementation in terms of parent-child > relationships wherever possible to protect locality. We can handle > the complications with recursive constraints later. Will do. > As for the FE template madness, I don't have a good answer... :-) > Just something that seemed like a good idea at the time. Like > anything else, though, it could be changed. There are very few > (none?) things in the library I consider set in stone... I don't see any pressing reason to change it; I was just curious about the decision. The tradeoffs between templates and OOP are usually efficiency related, but it looked to me like the templated FE code was just trading virtual functions for switch statements, and I didn't think either construct was particularly faster than the other. --- Roy Stogner |
From: Roy S. <roy...@ic...> - 2005-01-19 17:57:29
|
On Wed, 19 Jan 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Excellent! Will this be the default behavior and then let specific finite > elements (e.g. LAGRANGE) overload or re-implement the method? At the moment I'm making this a Clough-specific behavior, but making it the default after it's tested is my plan - reimplementation will just require adding more cases to the switch statements in fe_interface.C too. Remind me again why the finite element types are using template specialization instead of inheritance? ;-) > The L2 projection method is perfect since it addresses the problem for > arbitrary element types. In fact, it can be extended to use the normal > equations to handle p mismatch on an edge as well. I'd certainly like to see continuous hp-refinable function spaces in libMesh, but I'm not sure how much of a help this will be. Shouldn't the hierarchic basis functions be able to handle p mismatch for C0 elements? If element 1 has a polynomial degree p and adjacent element 2 has higher degree q, then isn't the constraint just "element 2 shape functions with degree p+1 through q and support on edge 1-2 get dropped"? No, on second thought it's more complicated than that - basis functions can have higher polynomial degree perpendicular to a face than they do in the face, and the case of simultaneous h refinement gets even uglier. Should I be preparing for this case in the automatic code, then? Right now I'm just doing the projection between parent and child degrees of freedom; perhaps I should be doing neighbor and child instead. > directly. Line 622 in fe_lagrange.C simply says keep constraints that look > like this: > > u_1 = 0.5*u_2 + 0.5*u_3 > > and ignore ones that look like this: > > u_1 = 0.9999 * u_1 + 1.e-6 * u_2 + 1.e-6 * u_3 > > Which we can recognize as u_1 = u_1, which is not a constraint at all... What about lines that look like this? u_1 = 0.99999 * u_3 That's not on the identity diagonal, yet it's the sort of situation I would expect to happen when a nodal degree of freedom on a Lagrangian quadratic child element corresponds to a midedge or midface DoF on the parent: the coefficient is still 1, but the local numbering doesn't match. --- Roy Stogner |