All,
first: Thanks y'all for working together on this!
> I've been thinking about this procedure, too. Boundary values should
> usually be the weakest part (periodic b.c. or hanging nodes should
> dominate), and that would work fine for most user codes. The only thing
> we have to make very clear then is that hanging nodes should be inserted
> before boundary constraints. And of course document that dofs that
> already are constrained are not touched in the
> interpolate_boundary_values function. I can make this change if we all
> agree on that this is the best way to do things.
I think this is the procedure with the least surprise. For sure, there are
always ways to get the "wrong" result, depending on how you define "wrong":
- With Martin's implementation, the result does not satisfy boundary
conditions any more for general boundary functions
- If it had been done the other way around, the solution would not satisfy
hanging node constraints at this point and consequently would not
satisfy the regularity properties of the element chosen (e.g. would not
be continuous despite using a Q1 element).
I will note that the situation becomes completely hopeless if you consider
curved boundaries since the edge midpoint (i.e. the hanging node) does in
general not lie on the mother edge. Consequently, the solution will not be H1
confirming anyway, regardless of the order of operations. So after Martin's
change, the solution will be neither conforming, nor have the right boundary
values. However, I still think that his implementation does the right thing
because it should be easier to debug a case where the solution has the wrong
boundary value than where it violates constraints. In practice, I can't see
why anyone would care which operation wins.
The problem Luca originally reported would not have happened if we had had a
more useful interface for ConstraintMatrix. The way we build up constraints
piecemeal through add_line/add_entry/set_inhomogeneity was clearly
ill-conceived from the beginning. We should have a call
ConstraintMatrix::add_constraint (const unsigned int constrained_dof,
const std::vector<std::pair<unsigned int, double> > &entries,
const double inhomogeneity = 0);
or similar. This would simply throw an exception if we tried to set a
constraint on a DoF that's already constrained (or we would just overwrite
it, the semantics are open to discussion). I've left myself a TODO in
constraint_matrix.h -- maybe I get to it, if someone were to beat me to it
that would be ok as well.
Best & thanks again to all
W.
-------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.tamu.edu/~bangerth/
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
