From: Nasser M. A. <nas...@mc...> - 2008-09-27 17:20:46
|
Hi: For some post-processing tasks I need to get access to gradients of variables not at the interior quadrature points, but on the nodes themselves. It seems like a simple enough thing to do, but I cannot find any existing functionality in Libmesh to do it. I trying to look through the source/documentation, but its been several hours and I am actually more confused now than when I started :) -- Nasser Mohieddin Abukhdeir Graduate Student (Materials Modeling Research Group) McGill University - Department of Chemical Engineering http://webpages.mcgill.ca/students/nabukh/web/ http://mmrg.chemeng.mcgill.ca/ |
From: David K. <dav...@gm...> - 2008-09-27 17:26:57
|
Well, the problem I think is that the gradients are not well-defined at node points, since finite element solutions are piecewise polynomials. One way to get an answer (John suggested this to me once) is to compute the gradients at quadrature points and then do an L2 projection of that solution, and then just sample the projected solution at the nodes. - Dave Nasser Mohieddin Abukhdeir wrote: > Hi: > For some post-processing tasks I need to get access to gradients of > variables not at the interior quadrature points, but on the nodes > themselves. It seems like a simple enough thing to do, but I cannot > find any existing functionality in Libmesh to do it. I trying to look > through the source/documentation, but its been several hours and I am > actually more confused now than when I started :) > |
From: Derek G. <fri...@gm...> - 2008-09-27 17:35:24
|
David's right.... On Sep 27, 2008, at 11:26 AM, David Knezevic wrote: > Well, the problem I think is that the gradients are not well-defined > at > node points, since finite element solutions are piecewise polynomials. Yep.. for your normal Lagrange elements the gradient is undefined on the element boundaries (including the nodes). Now, for C1 continuous elements (such as Clough-Toucher's, Hermite's, etc.) you should be able to get the value of the gradient at the nodes pretty easily: it should be in your solution vector. Obviously, I've never used these elements or I would know the answer to that... maybe Roy could fill us in. > One way to get an answer (John suggested this to me once) is to > compute > the gradients at quadrature points and then do an L2 projection of > that > solution, and then just sample the projected solution at the nodes. Yep... this is what's calle "Gradient Recovery". There are several methods for doing this... Derek |
From: Nasser M. A. <nas...@mc...> - 2008-09-27 18:24:58
|
I found a few references on gradient recovery, seems a bit expensive to do concurrently in my simulations (both directly or through some Galerkin method), so I'll just do it a posteriori as a separate FEM run. I found a decent paper describing direct and Galerkin methods. Switching to Hermite or any higher order elements is simply too expensive for my problems, at least with these gradient recovery methods I can delay the work until after the main computation is complete. BTW, here is a paper that I found useful in understanding these methods: http://www3.interscience.wiley.com/journal/110544657/abstract Nasser Mohieddin Abukhdeir Graduate Student (Materials Modeling Research Group) McGill University - Department of Chemical Engineering http://webpages.mcgill.ca/students/nabukh/web/ http://mmrg.chemeng.mcgill.ca/ Derek Gaston wrote: > David's right.... > > On Sep 27, 2008, at 11:26 AM, David Knezevic wrote: > > >> Well, the problem I think is that the gradients are not well-defined >> at >> node points, since finite element solutions are piecewise polynomials. >> > > Yep.. for your normal Lagrange elements the gradient is undefined on > the element boundaries (including the nodes). Now, for C1 continuous > elements (such as Clough-Toucher's, Hermite's, etc.) you should be > able to get the value of the gradient at the nodes pretty easily: it > should be in your solution vector. Obviously, I've never used these > elements or I would know the answer to that... maybe Roy could fill us > in. > > >> One way to get an answer (John suggested this to me once) is to >> compute >> the gradients at quadrature points and then do an L2 projection of >> that >> solution, and then just sample the projected solution at the nodes. >> > > Yep... this is what's calle "Gradient Recovery". There are several > methods for doing this... > > Derek > |
From: Roy S. <roy...@ic...> - 2008-09-29 00:15:48
|
On Sat, 27 Sep 2008, Derek Gaston wrote: > Now, for C1 continuous elements (such as Clough-Toucher's, > Hermite's, etc.) you should be able to get the value of the gradient > at the nodes pretty easily: it should be in your solution vector. > Obviously, I've never used these elements or I would know the answer > to that... maybe Roy could fill us in. Currently the C1 elements have degrees of freedom corresponding to x, y partial derivatives at vertex nodes, but they don't have any degrees of freedom which can give you the whole gradient at a second-order node. In general, the right way to get the well-defined gradient on a C1 node is the same way you'd get the cheapest O(h^p) gradient approximation on a C0 node: start from an adjoining element, and do a dof-weighted sum of its shape functions at that point. >> One way to get an answer (John suggested this to me once) is to >> compute the gradients at quadrature points and then do an L2 >> projection of that solution, and then just sample the projected >> solution at the nodes. > > Yep... this is what's calle "Gradient Recovery". There are several > methods for doing this... Yup. For instance: if you do just a sum (weighted by element size? angle? I forget) of the local gradients approaching a node on all elements touching that node, you do get an O(h^{p+1}) approximation that way. That's the first step to the popular Zienkiewicz-Zhu error estimator, and it may be good enough (and cheap enough) for whatever Nasser needs too. --- Roy |
From: Nasser M. A. <nas...@mc...> - 2008-09-29 00:42:06
|
I was thinking about one of those approaches, but it seems like the path of least resistance is to run postprocessing code afterwards that just uses FEM to solve N stationary problems, where N is the number of timesteps from the main simulation. I am already tight on resources (during the main simulation) and am a bit intimidated by writing code that deviates from the DiffSystem framework (at this point in my LibMesh career :). This approach has the benefit of: 1) Resulting in flexible post-processing code where I don't have to skimp on the number of variables because of memory constraints (that would exist if I was doing this during my main simulation). 2) I can change my mind as many times as I want and rerun this postprocessor in case I need different quantities. 3) This is an element independent implementation. Are there any drawbacks to doing this versus one of the Gradient Recovery methods? Nasser Mohieddin Abukhdeir Graduate Student (Materials Modeling Research Group) McGill University - Department of Chemical Engineering http://webpages.mcgill.ca/students/nabukh/web/ http://mmrg.chemeng.mcgill.ca/ Roy Stogner wrote: > On Sat, 27 Sep 2008, Derek Gaston wrote: > > >> Now, for C1 continuous elements (such as Clough-Toucher's, >> Hermite's, etc.) you should be able to get the value of the gradient >> at the nodes pretty easily: it should be in your solution vector. >> Obviously, I've never used these elements or I would know the answer >> to that... maybe Roy could fill us in. >> > > Currently the C1 elements have degrees of freedom corresponding to x, > y partial derivatives at vertex nodes, but they don't have any degrees > of freedom which can give you the whole gradient at a second-order > node. > > In general, the right way to get the well-defined gradient on a C1 > node is the same way you'd get the cheapest O(h^p) gradient > approximation on a C0 node: start from an adjoining element, and do a > dof-weighted sum of its shape functions at that point. > > >>> One way to get an answer (John suggested this to me once) is to >>> compute the gradients at quadrature points and then do an L2 >>> projection of that solution, and then just sample the projected >>> solution at the nodes. >>> >> Yep... this is what's calle "Gradient Recovery". There are several >> methods for doing this... >> > > Yup. For instance: if you do just a sum (weighted by element size? > angle? I forget) of the local gradients approaching a node on all > elements touching that node, you do get an O(h^{p+1}) approximation > that way. That's the first step to the popular Zienkiewicz-Zhu error > estimator, and it may be good enough (and cheap enough) for whatever > Nasser needs too. > --- > Roy > |
From: Roy S. <roy...@ic...> - 2008-09-29 00:50:44
|
On Sun, 28 Sep 2008, Nasser Mohieddin Abukhdeir wrote: > I was thinking about one of those approaches, but it seems like the path of > least resistance is to run postprocessing code afterwards that just uses FEM > to solve N stationary problems, where N is the number of timesteps from the > main simulation. That makes sense. That's the workflow I use for visualization - just save the results in a "lossless" bzipped xdr/xda format, then turn them into tecplot or gmv or whatever when I need to later. > I am already tight on resources (during the main simulation) and am > a bit intimidated by writing code that deviates from the DiffSystem > framework (at this point in my LibMesh career :). It's good to know that DiffSystem is making things easier for more new users. On the other hand, that's going to make me feel like a jerk when I break the API. Would it mess with your code too badly if I added an argument like "const FEMContext &context" to many of the DiffSystem methods (including the user-overridden ones), and then you had to access variables as "context.elem" instead of just "elem"? --- Roy |
From: Roy S. <roy...@ic...> - 2008-10-02 14:41:28
|
On Thu, 2 Oct 2008, Nasser Mohieddin Abukhdeir wrote: > I ran a simulation with the newly reformulated class structure, everything > is peachy. I'll get started on this PCAdaptiveTimeSolver over the weekend, > should I send an email to the devel list detailing this change and the game > plan for PCAdaptiveTimeSolver? Post the details and future plans on devel, but I just realized we should mention this on users too, since there's an API change that can affect people who follow the SVN head: Anyone who was previously instantiating an AdaptiveTimeSolver, we've renamed that due to refactoring. You'll want to build a TwostepTimeSolver instead - same functionality, just a different leaf class since we'll have more adaptive time solver implementations coming in the future. You'll also need to include "twostep_time_solver.h" for the new class. --- Roy |