On 11/15/2013 10:13 AM, Anastasia Konovalova wrote:
> Thanks for the answer, Robert! And excuse me for the long response.
>
> I'm not sure, that I understand you clearly.
> I have added changes from Git, and the computation of determinant in my
> programm is correct now. And your "black magic" code don't produce error
> any more. So you repaired this old bag, right?) And strain can be used for
> further computations now? Thank you very much!
Yes, it should work ok now.
> As for the terms and issues, I compute determinant and other values for the
> energy calculation. If it is possible to make terms for the energy density
> function for hyperelastic material, it would be great and I'll add the
> issue as soon as I can:)
Yes, that would be possible. Add the issue, please, so that it is not forgotten.
Cheers,
r.
> четверг, 7 ноября 2013 г., 16:34:59 UTC+6 пользователь Robert Cimrman
> написал:
>>
>> Hi Anastasia,
>>
>> great job, you have just found a very old bug!
>>
>> The following code (to be added in stress_strain() after i3 computation)
>> should
>> not raise errors but it does:
>>
>> # black magic here...
>> cache =
>> state.variables['u'].evaluate_cache['tl_common'][0][('__v_o1', 0,
>> 0, 'volume', None)]
>>
>> from sfepy.linalg import dot_sequences as dot
>> from sfepy.mechanics.tensors import get_full_indices
>>
>> c = dot(cache.mtx_f, cache.mtx_f, 'ATB')
>> e = 0.5 * (c - nm.eye(3))
>>
>> ii = get_full_indices(3)
>> c2 = cache.sym_c[..., ii, 0]
>> e2 = cache.green_strain[..., ii, 0]
>>
>> assert(nm.allclose(nm.sqrt(i3), cache.det_f, rtol=1e-15, atol=1e-15))
>> assert(nm.allclose(c, c2, rtol=1e-15, atol=1e-15))
>> assert(nm.allclose(e, e2, rtol=1e-15, atol=1e-15))
>>
>> The last line assertion fails, so the Green strain e2 as computed by the
>> dw_tl_he_neohook term call differs from e computed by definition (0.5 *
>> (F^T F
>> - I)).
>>
>> It has escaped me for so long as we have used the green strain only in
>> postprocessing to make "colorful images"... The Green strain has not been
>> used
>> in any further computations. The stresses are ok - they are computed from
>> the C
>> invariants directly.
>>
>> The problem (doubled out-of diagonal entries of E) is now fixed in the
>> repo and
>> your example gives positive i3 in all time steps.
>>
>> Thanks again for noticing the problem!
>>
>> Cheers,
>> r.
>>
>> PS: If you would like to have a term/terms for all the attributes of the
>> evaluate cache above (do 'print cache' to see them), add an issue, please.
>>
>> On 11/06/2013 08:45 AM, Anastasia Konovalova wrote:
>>> Hello Robert!
>>> I got a question about the energy calculation again :)
>>>
>>> If you remember, in my programm there is function get_Wh that is called
>> in
>>> postproc. This function compute the determinant of the right
>> Cauchy-Green
>>> deformation tensor. When I was trying to do a large deformation, a
>>> situation arose, where in some cells the determinant is very small or
>> even
>>> negative.
>>>
>>> I can analyze the results on these elements in Paraview, so I see that
>> the
>>> field of displacements of these cells are normal. I can deform my sample
>> by
>>> displacement, so I find out that the volume of element is normal too,
>>> determinant haven't been far from 1. But If I calculate the determinant
>>> from the components of strain tensor by hands, I'll get this negative
>>> value, so that the computation of the determinant is right too.
>>>
>>> It seems to me that this is the problem of the Green strain tensor
>>> computation.
>>>
>>> I attached the screenshot from Paraview, where you can see the
>> determinant
>>> value, the volume of element, and components of the Green strain tensor,
>>> the file of programm and a simple mesh, and file with the results (for
>>> paraview or some other similar programm) of one of the first steps, when
>>> the negative determinant appears.
>>>
>>> At the function get_Wh I get the the Green strain tensor, compute the
>> right
>>> Cauchy-Green deformation tensor (C = 2E+I), and get i3 =
>> nm.linalg.det(C).
>>> I need the sqrt of i3 to get a J = dV/dV_0, and the error appears there.
>>>
>>> Perhaps, it's too difficult to understand all this, so I will be glad to
>>> any advices. Thank you anyway.
>>>
>>
>>
>