On 08/31/2016 04:45 PM, Ronghai Wu wrote:
>
>
> 在 2016年8月30日星期二 UTC+2下午9:05:43，Robert Cimrman写道：
>>
>> On 08/30/2016 03:21 PM, Ronghai Wu wrote:
>>> (2) In the interactive linear elasticity example [1], we can have access
>> to
>>> the displacement, stress and strain arrays by
>>>
>>>
>>> vec = pb.solve()
>>> strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
>>> stress = pb.evaluate('ev_cauchy_stress.2.Omega(u)', mode='el_avg')
>>>
>>> Can I get access to the whole (global) displacement, stress and strain
>> arrays in parallel example?
>>
>> Yes, after the call to pl.create_gather_to_zero() you have the whole
>> solution
>> vector in the task 0 (sol = psol_full[...].copy()....). Then you can
>> evaluate
>> global strain and stress as usual.
>>
>>
> As I do the following things with 2D mesh 50*50 in [2]
>
>         sol = psol_full[...].copy()
>         print 'sol is:',sol
>         print 'sol shape is:',sol.shape
>         u = FieldVariable('u', 'parameter', field1,
>               primary_var_name='(set-to-None)')
>         remap = gfds[0].id_map
>         ug = sol[remap]
>         print 'ug is:', ug
>         print 'ug shape is:', ug.shape
>         strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
>         print 'strain is:', strain
>         print 'strain shape is:', strain
>
> I got {sol shape is: (7500,)} and {ug shape is: (5000,)}. Is {ug} the
> displacement at element center? and {ug[0]} is the displacement in x-axis?

sol contains both the displacements (2 components) and the pressure (it is the
Biot problem example, right?), that is why it is bigger by 3/2. ug are the
usual displacements in nodes in the correct sfepy ordering.

> But {strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg')}
> does not work, with ValueError: argument u not found!
> This way works for interactive linear elasticity example, why it does not
> work here?

It does not work, because pb is a task-local Problem, with variables named u_i,
p_i, instead of u, p.

In order to evaluate a term globally, do something like this:

         # Set the actual global solution vector to u.
         u.set_data(ug)

         integral = Integral('i', order=2*(max(order_u, order_p)))
         term = Term.new('ev_cauchy_strain(u)', integral, u.field.region, u=u)
         term.setup()

         strain = term.evaluate(mode='el_avg')

         print 'strain is:', strain
         print 'strain shape is:', strain.shape

Does that work for you?

r.

>>>
>>> [1]
>> http://sfepy.org/doc/tutorial.html#interactive-example-linear-elasticity
>>> [2]
>>>
>> http://sfepy.org/doc-devel/examples/multi_physics/biot_parallel_interactive.html#multi-physics-biot-parallel-interactive
>>>
>>>
>>> Thanks
>>> Ronghai
>>>
>>
>>
>