thanks for your pointers there

On 9/28/05, Michael Schindler <m-schindler@users.sourceforge.net > wrote:

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)=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 pseudo-DOF-number that collects the inhomoheneity in the
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>
        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
                    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;

But the easiest thing is to consider if you really need the the
inhomogeneous constraints. I had the simple case where I wanted a
pressure-driven flow, and the pressure had to fulfill similar
inhomogeneous constraints. But this could easily be avoided by adding
a constant gradient in the Navier-Stokes 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==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 dof-coupling 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
dof-coupling is built. This is exactly what happens to the

I take it you mean the dof_map, not the DofMap._dof_coupling.

hanging-node 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.

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ös.


Wout Ruijter