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. 