From: Roy Stogner <roystgnr@ic...>  20140407 03:50:31

On Sun, 6 Apr 2014, David Knezevic wrote: > Is it possible to impose a linear constraint on the variables in a > system using a DirichletBoundary object? Currently, no. > The case I'm thinking of is in elasticity, where you want set the > displacement in the normal (or tangential) direction on a surface. In > that case we want "n \dot (u,v,w) = normal_displacement" on the boundary. > > I've only used DirichletBoundary to impose independent conditions on > each variable, so I assume that a condition like "n \dot (u,v,w) = > normal_displacement" isn't possible yet? The lower level heterogeneous DofConstraints code will correctly impose such conditions once they've been consistently calculated for each DoF (this actually gets complicated at corners/edges...), but there's no way yet to tell the library to calculate such conditions. That's on my TODO list as well, but not high enough that I'd expect to get to it alone any time soon; I'd love to help you hash it out. > Roy: Do you have any comments on where to start, and how much work > it would require? Where to start is "what do we want the API to look like"; I'm not even sure if this can cleanly fit in DirichletBoundary or if we'd want to create some kind of new CoupledDirichletBoundary (or betternamed) addition. Most of the hard work in code can be done just by yetagain cloning (or for God's sake finally refactoring) our constraint loop code. The only catch I can see is when multiple constraints apply to the same dof. E.g. if two normal displacements act on a rotated 2D corner, we want that to correctly boil down to "the corner position is completely fixed" rather than "a constraint gets ignored" or "an error gets thrown"... and I'm not even sure if that's a new catch; I recall having to deal with the same problem when DirichletBoundary and PeriodicBoundary mix. But I don't know if there's an analog here to the solution I ended up going with there.  Roy 
From: David Knezevic <dknezevic@se...>  20140406 19:55:41

Is it possible to impose a linear constraint on the variables in a system using a DirichletBoundary object? The case I'm thinking of is in elasticity, where you want set the displacement in the normal (or tangential) direction on a surface. In that case we want "n \dot (u,v,w) = normal_displacement" on the boundary. I've only used DirichletBoundary to impose independent conditions on each variable, so I assume that a condition like "n \dot (u,v,w) = normal_displacement" isn't possible yet? If that's the case, I'd be happy to work on it. Roy: Do you have any comments on where to start, and how much work it would require? Thanks, David 
From: Roy Stogner <roystgnr@ic...>  20140407 03:50:31

On Sun, 6 Apr 2014, David Knezevic wrote: > Is it possible to impose a linear constraint on the variables in a > system using a DirichletBoundary object? Currently, no. > The case I'm thinking of is in elasticity, where you want set the > displacement in the normal (or tangential) direction on a surface. In > that case we want "n \dot (u,v,w) = normal_displacement" on the boundary. > > I've only used DirichletBoundary to impose independent conditions on > each variable, so I assume that a condition like "n \dot (u,v,w) = > normal_displacement" isn't possible yet? The lower level heterogeneous DofConstraints code will correctly impose such conditions once they've been consistently calculated for each DoF (this actually gets complicated at corners/edges...), but there's no way yet to tell the library to calculate such conditions. That's on my TODO list as well, but not high enough that I'd expect to get to it alone any time soon; I'd love to help you hash it out. > Roy: Do you have any comments on where to start, and how much work > it would require? Where to start is "what do we want the API to look like"; I'm not even sure if this can cleanly fit in DirichletBoundary or if we'd want to create some kind of new CoupledDirichletBoundary (or betternamed) addition. Most of the hard work in code can be done just by yetagain cloning (or for God's sake finally refactoring) our constraint loop code. The only catch I can see is when multiple constraints apply to the same dof. E.g. if two normal displacements act on a rotated 2D corner, we want that to correctly boil down to "the corner position is completely fixed" rather than "a constraint gets ignored" or "an error gets thrown"... and I'm not even sure if that's a new catch; I recall having to deal with the same problem when DirichletBoundary and PeriodicBoundary mix. But I don't know if there's an analog here to the solution I ended up going with there.  Roy 
From: David Knezevic <dknezevic@se...>  20140407 18:21:23

On 04/06/2014 11:50 PM, Roy Stogner wrote: > > On Sun, 6 Apr 2014, David Knezevic wrote: > >> Is it possible to impose a linear constraint on the variables in a >> system using a DirichletBoundary object? > > Currently, no. > >> The case I'm thinking of is in elasticity, where you want set the >> displacement in the normal (or tangential) direction on a surface. In >> that case we want "n \dot (u,v,w) = normal_displacement" on the >> boundary. >> >> I've only used DirichletBoundary to impose independent conditions on >> each variable, so I assume that a condition like "n \dot (u,v,w) = >> normal_displacement" isn't possible yet? > > The lower level heterogeneous DofConstraints code will correctly > impose such conditions once they've been consistently calculated for > each DoF (this actually gets complicated at corners/edges...), but > there's > no way yet to tell the library to calculate such conditions. > > That's on my TODO list as well, but not high enough that I'd expect to > get to it alone any time soon; I'd love to help you hash it out. OK, sounds good. >> Roy: Do you have any comments on where to start, and how much work >> it would require? > > Where to start is "what do we want the API to look like"; I'm not even > sure if this can cleanly fit in DirichletBoundary or if we'd want to > create some kind of new CoupledDirichletBoundary (or betternamed) > addition. I'd be happy with a CoupledDirichletBoundary class. Or perhaps LinearConstraintDirichletBoundary (if that isn't too verbose...), and if nonlinear couplings are needed later on, they could handled in a NonlinearConstraintDirichletBoundary class? Then, for each boundary constraint we want to apply we could attach functions f_i, i=1,...,n+1 to the the boundary condition object to give the constraint: \sum_{i=1}^{n} f_i(x,y,z) u_i(x,y,z) = f_{n+1}(x,y,z), where u_i is the i^th variable in our system. I guess this could be done with something like: dirichlet_bc.attach_constraint( std::map< unsigned int, FunctionBase* > lhs_functions, FunctionBase* rhs_function) The map keys would indicate which variable each lhs_function applies to, and any "unspecified" function would be assumed to be zero. How does that sound to you? > Most of the hard work in code can be done just by yetagain cloning > (or for God's sake finally refactoring) our constraint loop code. The > only catch I can see is when multiple constraints apply to the same > dof. E.g. if two normal displacements act on a rotated 2D corner, we > want > that to correctly boil down to "the corner position is completely > fixed" rather than "a constraint gets ignored" or "an error gets > thrown"... and I'm not even sure if that's a new catch; I recall > having to deal with the same problem when DirichletBoundary and > PeriodicBoundary mix. But I don't know if there's an analog here to > the solution I ended up going with there. OK, sounds good. I would hope that the solver can handle the case of multiple constraints on a dof. If the constraints are consistent, then the solver should do the right thing (i.e. completely fix the corner position in your example), and if the constraints are inconsistent, then the solver should complain about an illposed system. (The "double normal displacement on a corner node" would be a good test case for this.) David 
From: Roy Stogner <roystgnr@ic...>  20140407 20:08:36

On Mon, 7 Apr 2014, David Knezevic wrote: > On 04/06/2014 11:50 PM, Roy Stogner wrote: >> >> Where to start is "what do we want the API to look like"; I'm not even >> sure if this can cleanly fit in DirichletBoundary or if we'd want to >> create some kind of new CoupledDirichletBoundary (or betternamed) >> addition. > > > I'd be happy with a CoupledDirichletBoundary class. Or perhaps > LinearConstraintDirichletBoundary (if that isn't too verbose...), and if > nonlinear couplings are needed later on, they could handled in a > NonlinearConstraintDirichletBoundary class? > > Then, for each boundary constraint we want to apply we could attach > functions f_i, i=1,...,n+1 to the the boundary condition object to give > the constraint: > > \sum_{i=1}^{n} f_i(x,y,z) u_i(x,y,z) = f_{n+1}(x,y,z), > > where u_i is the i^th variable in our system. I guess this could be done > with something like: > > dirichlet_bc.attach_constraint( std::map< unsigned int, FunctionBase* > > lhs_functions, FunctionBase* rhs_function) > > The map keys would indicate which variable each lhs_function applies to, > and any "unspecified" function would be assumed to be zero. > > How does that sound to you? All reasonable. >> Most of the hard work in code can be done just by yetagain cloning >> (or for God's sake finally refactoring) our constraint loop code. The >> only catch I can see is when multiple constraints apply to the same >> dof. E.g. if two normal displacements act on a rotated 2D corner, we >> want >> that to correctly boil down to "the corner position is completely >> fixed" rather than "a constraint gets ignored" or "an error gets >> thrown"... and I'm not even sure if that's a new catch; I recall >> having to deal with the same problem when DirichletBoundary and >> PeriodicBoundary mix. But I don't know if there's an analog here to >> the solution I ended up going with there. > > OK, sounds good. > > I would hope that the solver can handle the case of multiple constraints > on a dof. Nope. The solver should easily handle multiple constraints *affecting* a dof, but since we write each constraint row into the matrix in place of the row corresponding to some dof's test function, we need to have a way of choosing which dof gets used for which constraint. > If the constraints are consistent, then the solver should do > the right thing (i.e. completely fix the corner position in your > example), and if the constraints are inconsistent, then the solver > should complain about an illposed system. (The "double normal > displacement on a corner node" would be a good test case for this.) The first nasty case I can think of is the combination of a simple constraint on a vertical or horizontal wall (deltax = 0 or deltay = 0, obviously gets applied to the deltax dof or the deltay dof respectively) with a previously applied compound constraint on a neighboring diagonal wall. Do we try to make sure that the compound constraint gets set on the dof which won't be subsequently constrained? Do we try to "move" compound constraints automatically when we're overwriting them?  Roy 
From: David Knezevic <dknezevic@se...>  20140407 21:05:35

>> >> I would hope that the solver can handle the case of multiple constraints >> on a dof. > > Nope. The solver should easily handle multiple constraints > *affecting* a dof, but since we write each constraint row into the > matrix in place of the row corresponding to some dof's test function, > we need to have a way of choosing which dof gets used for which > constraint. Ah, OK. Yes, I can see that makes things trickier since it doesn't seem like there's a general way to choose which dof to overwrite with the constraint row. >> If the constraints are consistent, then the solver should do >> the right thing (i.e. completely fix the corner position in your >> example), and if the constraints are inconsistent, then the solver >> should complain about an illposed system. (The "double normal >> displacement on a corner node" would be a good test case for this.) > > The first nasty case I can think of is the combination of a simple > constraint on a vertical or horizontal wall (deltax = 0 or deltay = 0, > obviously gets applied to the deltax dof or the deltay dof > respectively) with a previously applied compound constraint on a > neighboring diagonal wall. Do we try to make sure that the compound > constraint gets set on the dof which won't be subsequently > constrained? Do we try to "move" compound constraints automatically > when we're overwriting them? This does sound nasty... Do you see a good general solution here? This makes me wonder if I should use Lagrange multipliers (via SCALAR variables) to impose the type of constraint I have in mind instead. That way we add the extra constraint rows to the system rather than modifying existing rows, which seems like it would avoid some complications. Though the downside of Lagrange multipliers is that it adds dofs to the system and turns it into a saddlepoint problem, so certainly not a free lunch! David 
From: Roy Stogner <roystgnr@ic...>  20140407 21:30:38

On Mon, 7 Apr 2014, David Knezevic wrote: >> Nope. The solver should easily handle multiple constraints >> *affecting* a dof, but since we write each constraint row into the >> matrix in place of the row corresponding to some dof's test function, >> we need to have a way of choosing which dof gets used for which >> constraint. > > Ah, OK. Yes, I can see that makes things trickier since it doesn't seem like > there's a general way to choose which dof to overwrite with the constraint > row. Right. As long as there's just one constraint per dof things should work out fine, but once you get legitimate overlap you've got to have a way to resolve it. >> The first nasty case I can think of is the combination of a simple >> constraint on a vertical or horizontal wall (deltax = 0 or deltay = 0, >> obviously gets applied to the deltax dof or the deltay dof >> respectively) with a previously applied compound constraint on a >> neighboring diagonal wall. Do we try to make sure that the compound >> constraint gets set on the dof which won't be subsequently >> constrained? Do we try to "move" compound constraints automatically >> when we're overwriting them? > > This does sound nasty... Do you see a good general solution here? I think it should always be safe to "move" a constraint onto a notyetconstrained dof. The only remaining catch is the interaction with PeriodicBoundary constraints, where the constraints can involve more than one processor at once. > This makes me wonder if I should use Lagrange multipliers (via SCALAR > variables) to impose the type of constraint I have in mind instead. That way > we add the extra constraint rows to the system rather than modifying existing > rows, which seems like it would avoid some complications. Though the downside > of Lagrange multipliers is that it adds dofs to the system and turns it into > a saddlepoint problem, so certainly not a free lunch! I'd prefer to try to get the compound constraints working, but I'll understand if you want to back out. If you do want to futz with the Lagrange multiplier formulation, try playing around with the PETSc fieldsplit stuff I recently hooked us up to? I was never able to make it more efficient than tuned block Jacobi + ILU on incompressible flow problems, but it might do much better with a Schur block that only needs O((1/h)^{d1}) instead of O((1/h)^d) dofs.  Roy 
From: David Knezevic <dknezevic@se...>  20140407 21:44:04

>>> The first nasty case I can think of is the combination of a simple >>> constraint on a vertical or horizontal wall (deltax = 0 or deltay = 0, >>> obviously gets applied to the deltax dof or the deltay dof >>> respectively) with a previously applied compound constraint on a >>> neighboring diagonal wall. Do we try to make sure that the compound >>> constraint gets set on the dof which won't be subsequently >>> constrained? Do we try to "move" compound constraints automatically >>> when we're overwriting them? >> >> This does sound nasty... Do you see a good general solution here? > > I think it should always be safe to "move" a constraint onto a > notyetconstrained dof. The only remaining catch is the interaction > with PeriodicBoundary constraints, where the constraints can involve > more than one processor at once. OK, that approach doesn't sound too nasty... > >> This makes me wonder if I should use Lagrange multipliers (via SCALAR >> variables) to impose the type of constraint I have in mind instead. >> That way we add the extra constraint rows to the system rather than >> modifying existing rows, which seems like it would avoid some >> complications. Though the downside of Lagrange multipliers is that it >> adds dofs to the system and turns it into a saddlepoint problem, so >> certainly not a free lunch! > > I'd prefer to try to get the compound constraints working, but I'll > understand if you want to back out. OK. I do like the sound of compound constraints, so I'm happy to work on this too. Sounds like your "move the constraint" approach could work well. Do you have a first item on the todo list here? Maybe getting it to work in the case where we don't have any "conflicting constraints" first? As I said, I'm happy to work on this (though I'm not too familiar with the guts of the constraint code...) > If you do want to futz with the > Lagrange multiplier formulation, try playing around with the PETSc > fieldsplit stuff I recently hooked us up to? I was never able to make > it more efficient than tuned block Jacobi + ILU on incompressible flow > problems, but it might do much better with a Schur block that only > needs O((1/h)^{d1}) instead of O((1/h)^d) dofs. If I do use Lagrange multipliers for Dirichlet BCs, then I'll be sure to have a look at this. Thanks, David 
From: Roy Stogner <roystgnr@ic...>  20140408 02:41:31

On Mon, 7 Apr 2014, David Knezevic wrote: >>> This does sound nasty... Do you see a good general solution here? >> >> I think it should always be safe to "move" a constraint onto a >> notyetconstrained dof. The only remaining catch is the interaction >> with PeriodicBoundary constraints, where the constraints can involve >> more than one processor at once. > > OK, that approach doesn't sound too nasty... What if we get *really* ugly... If we define a constraint involving variables A and B that gets placed on B, then a constraint involving B and C that gets placed on C, then a constraint involving C... then we need to be able to move constraints recursively. But I think we can prove that recursive movement plus a tree search is enough to handle any possible (consistent) scenario... > Do you have a first item on the todo list here? Maybe getting it to work in > the case where we don't have any "conflicting constraints" first? Definitely! > As I said, I'm happy to work on this (though I'm not too familiar with the > guts of the constraint code...) Thanks.  Roy 