On Thu, 3 Mar 2011, Boyce Griffith wrote:
> On 3/2/11 11:55 AM, Roy Stogner wrote:
>>
>> On Wed, 2 Mar 2011, Boyce Griffith wrote:
>>
>>> On 3/1/11 6:41 PM, Roy Stogner wrote:
>>>>
>>>> On Thu, 2 Dec 2010, Boyce Griffith wrote:
>>>>
>>>>> The patched code does the "right thing" only when
>>>>> asymmetric_constraint_rows=true when calling
>>>>> DofMap::constrain_element_{matrix,vector,matrix_and_vector}(). This is
>>>>> the default, but if asymmetric_constraint_rows=false, the code will
>>>>> silently do the "wrong thing," probably yielding bizarre results. One
>>>>> approach would be to keep track of whether a particular constraint row
>>>>> is provided by the user or not, and to require
>>>>> asymmetric_constraint_rows==true when such constraints are present.
>>>>
>>>> Thanks for the heads-up. I've got some FEMSystem-based codes that
>>>> could probably be adjusted pretty easily to test that case, which we'd
>>>> want to fix before committing.
>>>
>>> Let me know if you want me to think about what to do for this case, or
>>> if it isn't clear what goes wrong if asymmetric_constraint_rows=false.
>>
>> I haven't looked at the problem too closely yet, but it might be
>> better to think about the math before the code.
>>
>> For the inhomogeneous cases: we basically turn K*x = F, x=C*y into
>> C^T*K*C*y=C^T*F. That's still a symmetric matrix. It's a smaller
>> symmetric matrix in the case where length(y) < length(x), but we don't
>> actually condense the matrix, so we end up having indeterminate
>> hanging node entries in y. In the symmetric case we make the matrix a
>> diagonal for those entries, and then we fix them up later after the
>> linear solve with enforce_constraints_exactly. In the asymmetric case
>> we reuse those entries' lines with the constraint equations so that
>> the solver will enforce them (if only inexactly).
>>
>> So now we've got x=C*y+H
>>
>> Which means if we do things analogously in the symmetric case we want
>>
>> C^T*K*(C*y+H)=C^T*F
>> C^T*K*C*y = C^T*(F-K*H)
>>
>> Hmm... this isn't what your patch does, and when I think about it we
>> probably *don't* want to do things analogously (it's nice that in the
>> homogeneous case x is also just y with some entries omitted, but the
>> "x=C*y+H" model doesn't preserve that property). Is there any
>> similarly clean way to formulate what your patch is doing?
>
> I'll try to think about this next week. I will also need to review what the
> patch actually does. :-)
Me too, actually. Looking at the code it appeared that the symmetric
homogeneous case should still work just fine, with any potential
breakage limited to the symmetric inhomogeneous case. But actually
running ex26 with your patch, the convergence is lost.
---
Roy
------------------------------------------------------------------------------
Colocation vs. Managed Hosting
A question and answer guide to determining the best fit
for your organization - today and in the future.
http://p.sf.net/sfu/internap-sfd2d
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
