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:', strainI 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....
Thanks Ronghai