Re: [Libmesh-users] Non-symmetric constrained element matrices in case of recursive constraints

[Roy S. <roy...@ic...>]

On Tue, 16 Jan 2018, David Knezevic wrote:

> I'm working with Jonas on this topic, and I just wanted to add a few
> - The original issue that Jonas ran into was indeed due to a "constraint
> cycle".

Sorry my email crossed with yours!

> - Our preferred fix is to detect constraint cycles and fix them ourselves,
> or at least to detect them and throw an error. What is the situation with
> error checking for constraint cycles at the moment?

I'm pretty sure there is none.

> This seems dangerous to me since libMesh can currently silently give
> the wrong answer.

Agreed.

> - If there is no error checking for constraint cycles, then I'd propose to
> add some optional error checking to DofMap

Certainly!  Except that I'd make it mandatory in dbg+devel mode.

> (and make sure this can be active in opt as well as dbg mode).

If that can be done sensibly, I suppose.  However, a dbg+devel-mode-only
check would be *easy* - dig into the function where we expand
recursive constraints, and add a quick libmesh_assert() that each of
the constrainer DoF indices doesn't equal the constrained DoF index at
each step at the expansion.  Doing that check in opt mode would be too
expensive, even optionally, IMHO, since you'd be testing your
do-we-do-testing flag a zillion times.

If you make a separate function to test constraints for recursion
before expanding them out, and only call that function after testing
an optional flag, that would be fine in opt mode, but more repetitive
to code.

> I guess this is as simple as a function
> that is called after process_constraints that checks that all entries in
> constrain rows are unconstrained.

Hmm... because *if* the recursion terminated properly, then we'd
expect that to be true after the expansion code was done?  That might
work.  I wouldn't trust it, though.  You're no longer proposing
testing whether the constraints are expected to break after
processing, you're proposing testing whether the processing looks
broken in a specific way.  It'll work for now, but then it may just
fail to trigger after the next time we refactor process_constraints().

I suppose either way you'll want to add a unit test or two that catch
cyclic constraints, and if we've got test coverage then I'm less
worried about breaking things.  So okay, this would be reasonable too.
---
Roy