Hello Michael,

On 9/29/05, Michael Schindler <m-schindler@users.sourceforge.net> wrote:
Hello Wout,
On 29.09.05, Wout Ruijter wrote:
> On 9/28/05, Michael Schindler <m-schindler@users.sourceforge.net> 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 pseudo-DOF-number 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,


in dof_map_constraints.C. It needs to know of the right-hand side.
I had also to enlarge

in order to also provide the rhs vector:

and so on ...
I need this for slip boundary conditions, where I have no

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

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
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 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 non-homogeneous 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 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
> > hanging-node 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?


  Michael Schindler

"A mathematician is a device for turning coffee into theorems"
  Paul Erdös.

This SF.Net email is sponsored by:
Power Architecture Resource Center: Free content, downloads, discussions,
and more. http://solutions.newsforge.com/ibmarch.tmpl
Libmesh-users mailing list