From: Kirk, B. (JSC-EG) <Ben...@na...> - 2008-10-31 20:55:27
|
Ok, I'm with you and agree that function is probably the way to go. I'd think a function call like extract_... is what you need. The full matrix method embeds (or used to) the constraint in the row (at the cost of introducing asymmetry), so even without enforce_constraints_exactly() you should recover the constrained value to some accuracy commensurate with your linear solver tolerance. The constrain_element_residual() should do the trick if it replaces whatever residual value you have with Re(i) = u(i) - u_constrained For either dirichlet or hanging-node type constraints... Right? ----- Original Message ----- From: Derek Gaston <der...@in...> To: Kirk, Benjamin (JSC-EG) Cc: lib...@li... <lib...@li...> Sent: Fri Oct 31 15:36:47 2008 Subject: Re: [Libmesh-devel] Contraints for JFNK hanging nodes and Dirchlet BC's On Oct 31, 2008, at 2:00 PM, Benjamin Kirk wrote: > I use it to extract element-based Dofs, and they get properly > constrained. > Would the residual value then be the current value - the constraned > value, > as you mention? I didn't know that method existed. It looks like it's pretty much achieving the same goal that enforce_constraints_exactly() does.... it at least means that you are calculating your next residual using the properly constrained value. Unfortunately, this doesn't help with calculating a constrained value. For instance, let's say that on solve #23 your solver decides it's converged... and libMesh writes out the solution. Do you have any idea what the value of your solution at the hanging nodes is? Unless you called enforce_constraints_exactly() before you wrote out the solution.... you don't. One idea would be for me to overwrite my own residual before I add it into the global residual.... I guess I could call is_constrained() on every node on the element and if it is then I could overwrite what residual is there with u - constrained_value before I add it in..... this is essentially the work I was thinking constrain_element_residual() would do. Derek |
From: Kirk, B. (JSC-EG) <Ben...@na...> - 2008-10-31 21:11:46
|
Right now the form is U_constrained is some linear combination of U_unconstrained... What about generalizing it to '+constant', where the constant (which could be a function of time) is your dirichlet value? We can still provide some user-friendly way to set the constant for dirichlet bcs... (?) -Ben |
From: Derek G. <der...@in...> - 2008-10-31 21:45:39
|
On Oct 31, 2008, at 3:11 PM, Kirk, Benjamin (JSC-EG) wrote: > Right now the form is > > U_constrained is some linear combination of U_unconstrained... > > What about generalizing it to '+constant', where the constant (which > could be a function of time) is your dirichlet value? > I think you'll actually want to make it more like U_constrained = constant... but yeah... that's the idea. > We can still provide some user-friendly way to set the constant for > dirichlet bcs... > Now the questions shifts to how to provide that constant. My recommendation for now is to call a function that takes (x, y, z, t, boundary_id). The problem is that nodes (and DOFs) don't have a boundary id. We could add one... but then the question shifts to how to apply those boundary_id's. My vote is that during equation_system_init() / reinit() we read the BoundaryInfo object on the mesh and apply the boundary_id of any face to all of it's nodes. That way boundary id's are automatically propagated from "sidesets" after adaptivity. Trying to "prolong" nodesets is a bad bad bad... I can get into why if anyone is interested. Let me just say that I've been there... and it's not good. What do you think? Can we pay the price of putting an unsigned int on every node? Derek |
From: Derek G. <fri...@gm...> - 2008-10-31 22:03:15
|
Another option is that when you provide the constraint... you provide the constant (and it's stored in the constraint system). This isn't a bad idea... and will probably use less storage. Derek On Oct 31, 2008, at 3:44 PM, Derek Gaston wrote: > On Oct 31, 2008, at 3:11 PM, Kirk, Benjamin (JSC-EG) wrote: >> Right now the form is >> >> U_constrained is some linear combination of U_unconstrained... >> >> What about generalizing it to '+constant', where the constant >> (which could be a function of time) is your dirichlet value? >> > > I think you'll actually want to make it more like U_constrained = > constant... but yeah... that's the idea. > >> We can still provide some user-friendly way to set the constant for >> dirichlet bcs... >> > > Now the questions shifts to how to provide that constant. My > recommendation for now is to call a function that takes (x, y, z, t, > boundary_id). The problem is that nodes (and DOFs) don't have a > boundary id. We could add one... but then the question shifts to > how to apply those boundary_id's. My vote is that during > equation_system_init() / reinit() we read the BoundaryInfo object on > the mesh and apply the boundary_id of any face to all of it's > nodes. That way boundary id's are automatically propagated from > "sidesets" after adaptivity. Trying to "prolong" nodesets is a bad > bad bad... I can get into why if anyone is interested. Let me just > say that I've been there... and it's not good. > > What do you think? > > Can we pay the price of putting an unsigned int on every node? > > Derek > > ------------------------------------------------------------------------- > 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: Benjamin K. <ben...@na...> - 2008-10-31 22:08:30
|
> Another option is that when you provide the constraint... you provide the > constant (and it's stored in the constraint system). This isn't a bad idea... > and will probably use less storage. That's what I was envisioning... As for the node bc ids, > Now the questions shifts to how to provide that constant. My recommendation > for now is to call a function that takes (x, y, z, t, boundary_id). The > problem is that nodes (and DOFs) don't have a boundary id. We could add > one... but then the question shifts to how to apply those boundary_id's. My > vote is that during equation_system_init() / reinit() we read the BoundaryInfo > object on the mesh and apply the boundary_id of any face to all of it's nodes. > That way boundary id's are automatically propagated from "sidesets" after > adaptivity. Trying to "prolong" nodesets is a bad bad bad... I can get into > why if anyone is interested. Let me just say that I've been there... and it's > not good. > > What do you think? > > Can we pay the price of putting an unsigned int on every node? What if instead after each refinement/coarsening you (i) clear a node list, and then loop over active elements & their faces to build up a derived node bc list? You regenerate it instead of trying to prolong nodes... -Ben |
From: Derek G. <fri...@gm...> - 2008-11-01 22:20:58
|
On Fri, Oct 31, 2008 at 4:08 PM, Benjamin Kirk <ben...@na...>wrote: > > What if instead after each refinement/coarsening you (i) clear a node list, > and then loop over active elements & their faces to build up a derived node > bc list? > > You regenerate it instead of trying to prolong nodes... > Right... my point was that you should regenerate that list using _sidesets_ not nodesets... I think we're on the same page. So... should we have users attach a function that regenerates this node bc list? That way it can automatically be called when the mesh changes? There should also be a way to manually trigger the list regeneration... for instance, for Dirichlet BC's that change with time... user code will probably want to signal that the bc list should be rebuilt. Derek |
From: Derek G. <der...@in...> - 2008-10-31 21:04:25
|
On Oct 31, 2008, at 2:55 PM, Kirk, Benjamin (JSC-EG) wrote: > Ok, I'm with you and agree that function is probably the way to go. > I'd think a function call like extract_... is what you need. > > The full matrix method embeds (or used to) the constraint in the row > (at the cost of introducing asymmetry), so even without > enforce_constraints_exactly() you should recover the constrained > value to some accuracy commensurate with your linear solver tolerance. > > The constrain_element_residual() should do the trick if it replaces > whatever residual value you have with > > Re(i) = u(i) - u_constrained > > For either dirichlet or hanging-node type constraints... Right? > Well... Dirichlet constraints are a tad more difficult because you have to supply some other value that the library doesn't (yet) know about. If we could come up with a good API for telling the constraint system about the values.... then we'd be in business. Derek |