From: Wout Ruijter <woutruijter@gm...>  20050927 20:06:00
Attachments:
Message as HTML

Dear sirs, Three questions regarding periodic boundary conditions: 1  Michael already showed how to implement periodic bc's of type u(l)=3Du(= r), how do you add u(l)u(r)=3Dx? Is that a matter of modifying the righthandsi= de? I mean, is the rhs taken into account for constrained dof? 2  Is there any way in which of two elements with the same error estimate only one will be refined? If this is not the case I would be able to ensure matching meshes by modifying the error estimates, or should I modify the refinement flags directly? 3 Am I right to presume that in the case that no interpolated constraints are used the _dof_coupling can be used for x=3D=3D0? I guess this could speed up the whole matter in cases with many couplings. Thanks in advance for your opinion Wout 
From: Roy Stogner <roystgnr@ic...>  20050927 21:22:18

On Tue, 27 Sep 2005, Wout Ruijter wrote: > Three questions regarding periodic boundary conditions: > 1  Michael already showed how to implement periodic bc's of type u(l)=u(r), > how do you add u(l)u(r)=x? Is that a matter of modifying the righthandside? > I mean, is the rhs taken into account for constrained dof? You'll have to do so explicitly. I don't recall whether Michael finally ended up using a penalty method or building libMesh constraints to do his boundary conditions, but only the former method will work for your case. I don't think the DofConstraintRow objects can include right hand side terms. > 2  Is there any way in which of two elements with the same error estimate > only one will be refined? If you have a mesh on which 400 elements have error 2 and 600 have error 1, then AFAIK telling libMesh to refine the worst 30% of elements will cause a random 300 elements from the first group to be flagged for refinement. Normally this is barely noticeable: any elements "left behind" by one refinement step are likely to be the first elements in line for refinement on the next step. If you need two meshes to match exactly or one mesh to be perfectly symmetric, though, this will be a problem. > If this is not the case I would be able to ensure matching meshes by > modifying the error estimates, or should I modify the refinement > flags directly? Either way is good; I'm not sure which way would be quicker to code. To be safe make sure to modify the flags conservatively: only set new refinement flags and remove coarsening flags, never set coarsening or remove refinement. > 3 Am I right to presume that in the case that no interpolated constraints > are used the _dof_coupling can be used for x==0? > I guess this could speed up the whole matter in cases with many couplings. I don't think I understand this question.  Roy Stogner 
From: Michael Povolotskyi <povolotskyi@in...>  20050928 09:05:45

Roy Stogner wrote: > On Tue, 27 Sep 2005, Wout Ruijter wrote: > >> Three questions regarding periodic boundary conditions: >> 1  Michael already showed how to implement periodic bc's of type >> u(l)=u(r), >> how do you add u(l)u(r)=x? Is that a matter of modifying the >> righthandside? >> I mean, is the rhs taken into account for constrained dof? > > > You'll have to do so explicitly. I don't recall whether Michael > finally ended up using a penalty method or building libMesh > constraints to do his boundary conditions, but only the former method > will work for your case. I don't think the DofConstraintRow objects > can include right hand side terms. > Yes, I implemented only libMesh constraints without the right hand side terms. Michael. 
From: Michael Schindler <mschindler@us...>  20050928 11:38:01

Hello, although I am probably not the Michael that Wout refers to ;) I have also tried to implement the constraints for several purposes On 27.09.05, Roy Stogner wrote: > On Tue, 27 Sep 2005, Wout Ruijter wrote: > > >Three questions regarding periodic boundary conditions: > >1  Michael already showed how to implement periodic bc's of type > >u(l)=u(r), > >how do you add u(l)u(r)=x? Is that a matter of modifying the > >righthandside? > >I mean, is the rhs taken into account for constrained dof? Exactly. At the moment only homogeneous constraints are supported. If you really need inhomogeneous constraints you need two things: 1. A marker pseudoDOFnumber that collects the inhomoheneity in the ConstraintsRow 2. patch the source code of libmesh. The function that really does the constraint resolving. But the easiest thing is to consider if you really need the the inhomogeneous constraints. I had the simple case where I wanted a pressuredriven flow, and the pressure had to fulfill similar inhomogeneous constraints. But this could easily be avoided by adding a constant gradient in the NavierStokes equation. What kind of system do you need the constraints for? > >2  Is there any way in which of two elements with the same error estimate > >only one will be refined? > > If you have a mesh on which 400 elements have error 2 and 600 have > error 1, then AFAIK telling libMesh to refine the worst 30% of > elements will cause a random 300 elements from the first group to be > flagged for refinement. Normally this is barely noticeable: any > elements "left behind" by one refinement step are likely to be the > first elements in line for refinement on the next step. If you need > two meshes to match exactly or one mesh to be perfectly symmetric, > though, this will be a problem. > > >If this is not the case I would be able to ensure matching meshes by > >modifying the error estimates, or should I modify the refinement > >flags directly? > > Either way is good; I'm not sure which way would be quicker to code. > To be safe make sure to modify the flags conservatively: only set new > refinement flags and remove coarsening flags, never set coarsening or > remove refinement. > > >3 Am I right to presume that in the case that no interpolated constraints > >are used the _dof_coupling can be used for x==0? > >I guess this could speed up the whole matter in cases with many couplings. This is a serious problem. In most cases the _dof_coupling is already built when you impose the constraints. Then, the Petsc solver will have to enlarge its dofcoupling every time you fill the matrix at an unexpected position. This makes the story damn slow. Therefore, you will need to create the constraints _before_ the dofcoupling is built. This is exactly what happens to the hangingnode constraints. I took me a long time to recognize this issue, and maybe there is another solution, but I ended up by drilling a hole into the complete library. I now have something like "attach_constrain_function" which works similar as "attach_assemble_function" and which is called from inside the library right after the hanging nodes have been built. Michael Schindler  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 
From: Wout Ruijter <woutruijter@gm...>  20050929 10:27:13
Attachments:
Message as HTML

Micheal, thanks for your pointers there On 9/28/05, Michael Schindler <mschindler@...> wrote: > > Hello, > > although I am probably not the Michael that Wout refers to ;) I have > also tried to implement the constraints for several purposes It seems that being a Micheal is periodically coupled to an abnormal interest in periodic boundary conditions ;) On 27.09.05, Roy Stogner wrote: > > On Tue, 27 Sep 2005, Wout Ruijter wrote: > > > > >Three questions regarding periodic boundary conditions: > > >1  Michael already showed how to implement periodic bc's of type > > >u(l)=3Du(r), > > >how do you add u(l)u(r)=3Dx? Is that a matter of modifying the > > >righthandside? > > >I mean, is the rhs taken into account for constrained dof? > > Exactly. At the moment only homogeneous constraints are supported. If > you really need inhomogeneous constraints you need two things: > 1. A marker pseudoDOFnumber that collects the inhomoheneity in the > ConstraintsRow > 2. patch the source code of libmesh. The function that really does the > constraint resolving. I'm trying, in dof_map.h >>>>>>> //typedef std::map<unsigned int, Real> DofConstraintRow; class DofConstraintRow : public std::map<unsigned int,Real> { public: Real rhs; DofConstraintRow(){rhs =3D0.0;}; }; <<<<<<< and in dof_map_constraints.C >>>>>>> for (unsigned int i=3D0; i<elem_dofs.size(); i++) if (this>is_constrained_dof(elem_dofs[i])) { // If the DOF is constrained DofConstraints::const_iterator pos =3D _dof_constraints.find(elem_dofs[i]); // find the constraint row assert (pos !=3D _dof_constraints.end()); const DofConstraintRow& constraint_row =3D pos>second; // assert (!constraint_row.empty()); // check that it is not empty // this should not be used because a // constraint of type u(n)=3Drhs // should be legal rhs(i) =3D constraint_row.rhs; // rhs(i) =3D 0.0; // std::cout<< "rhs =3D " <<rhs(i)<<std::endl; } <<<<<<<< But the easiest thing is to consider if you really need the the > inhomogeneous constraints. I had the simple case where I wanted a > pressuredriven flow, and the pressure had to fulfill similar > inhomogeneous constraints. But this could easily be avoided by adding > a constant gradient in the NavierStokes equation. > What kind of system do you need the constraints for? I need them for elastostatics, so nodes on opposite sides of a unit cell moving such that the unit cell remains repeatable. > >2  Is there any way in which of two elements with the same error > estimate > > >only one will be refined? > > > > If you have a mesh on which 400 elements have error 2 and 600 have > > error 1, then AFAIK telling libMesh to refine the worst 30% of > > elements will cause a random 300 elements from the first group to be > > flagged for refinement. Normally this is barely noticeable: any > > elements "left behind" by one refinement step are likely to be the > > first elements in line for refinement on the next step. If you need > > two meshes to match exactly or one mesh to be perfectly symmetric, > > though, this will be a problem. > > > > >If this is not the case I would be able to ensure matching meshes by > > >modifying the error estimates, or should I modify the refinement > > >flags directly? > > > > Either way is good; I'm not sure which way would be quicker to code. > > To be safe make sure to modify the flags conservatively: only set new > > refinement flags and remove coarsening flags, never set coarsening or > > remove refinement. > > > > >3 Am I right to presume that in the case that no interpolated > constraints > > >are used the _dof_coupling can be used for x=3D=3D0? > > >I guess this could speed up the whole matter in cases with many > couplings. > > This is a serious problem. In most cases the _dof_coupling is already > built when you impose the constraints. Then, the Petsc solver will > have to enlarge its dofcoupling every time you fill the matrix at an > unexpected position. This makes the story damn slow. > Therefore, you will need to create the constraints _before_ the > dofcoupling is built. This is exactly what happens to the I take it you mean the dof_map, not the DofMap._dof_coupling. hangingnode constraints. I took me a long time to recognize this issue, and maybe there is > another solution, but I ended up by drilling a hole into the complete > library. I now have something like "attach_constrain_function" which work= s > similar as "attach_assemble_function" and which is called from inside > the library right after the hanging nodes have been built. Thanks for pointing this out, this seems a difficult one, but I think your solution would be the only way to have the desirable behaviour that: 1. The periodic or any new constraints override the hanging node couplings 2. They are taken into account when building the dof map Michael Schindler > >  > "A mathematician is a device for turning coffee into theorems" > Paul Erd=F6s. Greetings Wout Ruijter 
From: Wout Ruijter <woutruijter@gm...>  20050928 11:16:54
Attachments:
Message as HTML

Roy, Thanks for the comments On 9/27/05, Roy Stogner <roystgnr@...> wrote: > > On Tue, 27 Sep 2005, Wout Ruijter wrote: > > > Three questions regarding periodic boundary conditions: > > 1  Michael already showed how to implement periodic bc's of type > u(l)=3Du(r), > > how do you add u(l)u(r)=3Dx? Is that a matter of modifying the > righthandside? > > I mean, is the rhs taken into account for constrained dof? > > You'll have to do so explicitly. I don't recall whether Michael > finally ended up using a penalty method or building libMesh > constraints to do his boundary conditions, but only the former method > will work for your case. I don't think the DofConstraintRow objects > can include right hand side terms. I guess one could create loose nodes, or dofs, for that matter, which could be forced using a penalty method and which provide the constant that is needed. I see now that there has been a discussion on that before, I'll read some o= n it. > 2  Is there any way in which of two elements with the same error estimat= e > > only one will be refined? > > If you have a mesh on which 400 elements have error 2 and 600 have > error 1, then AFAIK telling libMesh to refine the worst 30% of > elements will cause a random 300 elements from the first group to be > flagged for refinement. Normally this is barely noticeable: any > elements "left behind" by one refinement step are likely to be the > first elements in line for refinement on the next step. If you need > two meshes to match exactly or one mesh to be perfectly symmetric, > though, this will be a problem. > > > If this is not the case I would be able to ensure matching meshes by > > modifying the error estimates, or should I modify the refinement > > flags directly? > > Either way is good; I'm not sure which way would be quicker to code. > To be safe make sure to modify the flags conservatively: only set new > refinement flags and remove coarsening flags, never set coarsening or > remove refinement. Of course. > 3 Am I right to presume that in the case that no interpolated constraint= s > > are used the _dof_coupling can be used for x=3D=3D0? > > I guess this could speed up the whole matter in cases with many > couplings. > > I don't think I understand this question. >  > Roy Stogner > Sorry, the last one is fairly vague, I mean, if a node is linked to 1 other node with u(l)=3Du(r), one could use the _dof_coupling instead of a constraint, whereas if it is linked to a couple of nodes (when it is a hanging node), or when u(l)u(r)=3Dx, x!=3D0,= one would need a constraint row to handle it. However, this may be fairly obsolete, because I don't think using _dof_coupling is a good idea, as it seems to allocate N^2 memory. Sorry for being thick, but I don't understand how the penalty method you suggested in an earlier mail would lead to a coupling of DOF (not a constraint) without having to do a lookup of the corresponding dof on the other side of the mesh and/or meddling with the sparsity structure. If it does that, I would happily implement and use it. Thanks a lot for your time, Wout 
From: Roy Stogner <roystgnr@ic...>  20050928 12:58:07

On Wed, 28 Sep 2005, Wout Ruijter wrote: > Sorry for being thick, but I don't understand how the penalty method you > suggested in an earlier mail would lead to a coupling of DOF (not a > constraint) without having to do a lookup of the corresponding dof on the > other side of the mesh and/or meddling with the sparsity structure. If it > does that, I would happily implement and use it. The penalty method would require a lookup of the corresponding element on the other side of the mesh, although not necessarily of every corresponding DoF. Such a lookup table would be fast to update after refinement, though  you'd just need to loop through the previous mesh's lookup table, delete any coarsened entries and check which children of refined elements correspond to which children of the parent's corresponding element.  Roy 
From: Michael Schindler <mschindler@us...>  20050929 11:00:48

Hello Wout, On 29.09.05, Wout Ruijter wrote: > On 9/28/05, Michael Schindler <mschindler@...> wrote: > On 27.09.05, Roy Stogner wrote: > > > On Tue, 27 Sep 2005, Wout Ruijter wrote: > > > > > > >Three questions regarding periodic boundary conditions: > > > >1  Michael already showed how to implement periodic bc's of type > > > >u(l)=u(r), > > > >how do you add u(l)u(r)=x? Is that a matter of modifying the > > > >righthandside? > > > >I mean, is the rhs taken into account for constrained dof? > > > > Exactly. At the moment only homogeneous constraints are supported. If > > you really need inhomogeneous constraints you need two things: > > 1. A marker pseudoDOFnumber that collects the inhomoheneity in the > > ConstraintsRow > > 2. patch the source code of libmesh. The function that really does the > > constraint resolving. > I'm trying, in dof_map.h > >>>>>>> > //typedef std::map<unsigned int, Real> DofConstraintRow; > > class DofConstraintRow : public std::map<unsigned int,Real> > { > public: > Real rhs; > DofConstraintRow(){rhs =0.0;}; > }; > <<<<<<< > and in dof_map_constraints.C > >>>>>>> > for (unsigned int i=0; i<elem_dofs.size(); i++) > if (this>is_constrained_dof(elem_dofs[i])) > { > // If the DOF is constrained > DofConstraints::const_iterator > pos = _dof_constraints.find(elem_dofs[i]); // find the constraint row > > assert (pos != _dof_constraints.end()); > > const DofConstraintRow& constraint_row = pos>second; > > // assert (!constraint_row.empty()); // check that it is not empty > // this should not be used because a > // constraint of type u(n)=rhs > // should be legal > > rhs(i) = constraint_row.rhs; > // rhs(i) = 0.0; > > // std::cout<< "rhs = " <<rhs(i)<<std::endl; > } > <<<<<<<< and of course, you will need to change the routine that really does the job, DofMap::constrain_element_matrix_and_vector in dof_map_constraints.C. It needs to know of the righthand side. I had also to enlarge build_constraint_matrix in order to also provide the rhs vector: build_constraint_matrix_and_vector and so on ... I need this for slip boundary conditions, where I have no alternative. > But the easiest thing is to consider if you really need the the > > inhomogeneous constraints. I had the simple case where I wanted a > > pressuredriven flow, and the pressure had to fulfill similar > > inhomogeneous constraints. But this could easily be avoided by adding > > a constant gradient in the NavierStokes equation. > > What kind of system do you need the constraints for? > > > I need them for elastostatics, so nodes on opposite sides of a unit cell > moving such that the unit cell remains repeatable. You have a position variable x(xi) for the nodes, right? x is the cartesian component, xi a reference coordinate, in your case also the cartesian one. To make the whole thing repeatable your constraints for nodes i and j (on a left/right boundary resp.) look like x_i = x_j + X where X is the periodicity vector. Why not replace the variable x by something that takes X already into account? When the left reference boundary is xi=0 and the right one xi=1, then you can define a new position variable p(xi) = x(xi) + xi * X grad(p) = grad(x) + X and you will end up with homogeneous constraints p_i = p_j This is exactly what I have done with a constant pressure gradient. And you can do that in arbitrary dimensions as long as you fill space with a Bravais grid. > > This is a serious problem. In most cases the _dof_coupling is already > > built when you impose the constraints. Then, the Petsc solver will > > have to enlarge its dofcoupling every time you fill the matrix at an > > unexpected position. This makes the story damn slow. > > Therefore, you will need to create the constraints _before_ the > > dofcoupling is built. This is exactly what happens to the > > hangingnode constraints. > > I take it you mean the dof_map, not the DofMap._dof_coupling. I am not quite sure at the moment. What I really mean is the sparsity pattern in the Petsc matrix. You do not really want to fill the matrix with values at unexpected positions. This makes Petsc slow. > I took me a long time to recognize this issue, and maybe there is > > another solution, but I ended up by drilling a hole into the complete > > library. I now have something like "attach_constrain_function" which works > > similar as "attach_assemble_function" and which is called from inside > > the library right after the hanging nodes have been built. > > Thanks for pointing this out, this seems a difficult one, but I think your > solution would be the only way to have the desirable behaviour that: > 1. The periodic or any new constraints override the hanging node couplings > 2. They are taken into account when building the dof map I have not found another way to do that. Greetings, Michael Schindler  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 
From: Wout Ruijter <woutruijter@gm...>  20050929 14:35:52
Attachments:
Message as HTML

Hello Michael, On 9/29/05, Michael Schindler <mschindler@...> wrote: > > Hello Wout, > On 29.09.05, Wout Ruijter wrote: > > On 9/28/05, Michael Schindler <mschindler@...> wrote= : > > On 27.09.05, Roy Stogner wrote: > > > > On Tue, 27 Sep 2005, Wout Ruijter wrote: > > > > > > > > >Three questions regarding periodic boundary conditions: > > > > >1  Michael already showed how to implement periodic bc's of type > > > > >u(l)=3Du(r), > > > > >how do you add u(l)u(r)=3Dx? Is that a matter of modifying the > > > > >righthandside? > > > > >I mean, is the rhs taken into account for constrained dof? > > > > > > Exactly. At the moment only homogeneous constraints are supported. If > > > you really need inhomogeneous constraints you need two things: > > > 1. A marker pseudoDOFnumber that collects the inhomoheneity in the > > > ConstraintsRow > > > 2. patch the source code of libmesh. The function that really does th= e > > > constraint resolving. > > > I'm trying, in dof_map.h > > >>>>>>> > > //typedef std::map<unsigned int, Real> DofConstraintRow; > > > > class DofConstraintRow : public std::map<unsigned int,Real> > > { > > public: > > Real rhs; > > DofConstraintRow(){rhs =3D0.0;}; > > }; > > <<<<<<< > > > and in dof_map_constraints.C > > >>>>>>> > > for (unsigned int i=3D0; i<elem_dofs.size(); i++) > > if (this>is_constrained_dof(elem_dofs[i])) > > { > > // If the DOF is constrained > > DofConstraints::const_iterator > > pos =3D _dof_constraints.find(elem_dofs[i]); // find the constraint row > > > > assert (pos !=3D _dof_constraints.end()); > > > > const DofConstraintRow& constraint_row =3D pos>second; > > > > // assert (!constraint_row.empty()); // check that it is not empty > > // this should not be used because a > > // constraint of type u(n)=3Drhs > > // should be legal > > > > rhs(i) =3D constraint_row.rhs; > > // rhs(i) =3D 0.0; > > > > // std::cout<< "rhs =3D " <<rhs(i)<<std::endl; > > } > > <<<<<<<< > > and of course, you will need to change the routine that really does > the job, > > DofMap::constrain_element_matrix_and_vector > > in dof_map_constraints.C. It needs to know of the righthand side. > I had also to enlarge > > build_constraint_matrix > in order to also provide the rhs vector: > build_constraint_matrix_and_vector > > and so on ... > I need this for slip boundary conditions, where I have no > alternative. > > > But the easiest thing is to consider if you really need the the > > > inhomogeneous constraints. I had the simple case where I wanted a > > > pressuredriven flow, and the pressure had to fulfill similar > > > inhomogeneous constraints. But this could easily be avoided by adding > > > a constant gradient in the NavierStokes equation. > > > What kind of system do you need the constraints for? > > > > > > I need them for elastostatics, so nodes on opposite sides of a unit cel= l > > moving such that the unit cell remains repeatable. > > You have a position variable x(xi) for the nodes, right? > x is the cartesian component, xi a reference coordinate, in your case > also the cartesian one. > To make the whole thing repeatable your constraints for nodes i and j > (on a left/right boundary resp.) look like > > x_i =3D x_j + X > > where X is the periodicity vector. > Why not replace the variable x by something that takes X already into > account? > When the left reference boundary is xi=3D0 and the right one xi=3D1, then > you can define a new position variable > > p(xi) =3D x(xi) + xi * X > grad(p) =3D grad(x) + X > > and you will end up with homogeneous constraints > > p_i =3D p_j > > This is exactly what I have done with a constant pressure gradient. > And you can do that in arbitrary dimensions as long as you fill space > with a Bravais grid. This seems a nice way to deal with it, however, this means I have to reformulate my assembly and postprocessing functions in terms of p, right? I'll have a look at it, but I will first finalize the nonhomogeneous constraints. > > This is a serious problem. In most cases the _dof_coupling is already > > > built when you impose the constraints. Then, the Petsc solver will > > > have to enlarge its dofcoupling every time you fill the matrix at an > > > unexpected position. This makes the story damn slow. > > > Therefore, you will need to create the constraints _before_ the > > > dofcoupling is built. This is exactly what happens to the > > > hangingnode constraints. > > > > I take it you mean the dof_map, not the DofMap._dof_coupling. > > I am not quite sure at the moment. What I really mean is the sparsity > pattern in the Petsc matrix. You do not really want to fill the matrix > with values at unexpected positions. This makes Petsc slow. Ok, I think I know where to look. > I took me a long time to recognize this issue, and maybe there is > > > another solution, but I ended up by drilling a hole into the complete > > > library. I now have something like "attach_constrain_function" which > works > > > similar as "attach_assemble_function" and which is called from inside > > > the library right after the hanging nodes have been built. > > > > Thanks for pointing this out, this seems a difficult one, but I think > your > > solution would be the only way to have the desirable behaviour that: > > 1. The periodic or any new constraints override the hanging node > couplings > > 2. They are taken into account when building the dof map > > I have not found another way to do that. I must say I find it quite a nice solution, apart from the fact that does not come standard. So why search for an alternative? Greetings Wout Greetings, > Michael Schindler > >  > "A mathematician is a device for turning coffee into theorems" > Paul Erd=F6s. > > >  > This SF.Net email is sponsored by: > Power Architecture Resource Center: Free content, downloads, discussions, > and more. http://solutions.newsforge.com/ibmarch.tmpl > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > 