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

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.

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,

DofMap::constrain_element_matrix_and_vector

in dof_map_constraints.C. It needs to know of the right-hand 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

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

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

Greetings

Wout

Greetings

Wout

Greetings,

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

Libmesh-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/libmesh-users