Hi Ronghai,
On 11/06/2014 01:16 PM, Ronghai Wu wrote:
Dear sfepy users and developers,
I am using sfepy for solving an 2D inclusion elastic problem. So, I implemented based on "linear_elasticity/prestress_fibres" example in interactive way. However, the prestress in my problem is not constant, so I manually assign values. I have two questions:
- why when the integral order is 3, the array shape of prestress.datas[( 'Omega', 'i')][0]['val'] is (nxny,4,3,1), so when I assign by prestress_value, I have to make it the same shape. If change integral order to 1, the shape changes to (nxny,8,3,1), why?
This seems strange. You have quadrilateral elements, right?
./script/plot_quadratures.py -g 2_4 -n 3 sfepy: geometry: 2_4 order: 3 num. points: 4 true_order: 3
./script/plot_quadratures.py -g 2_4 -n 1 sfepy: geometry: 2_4 order: 1 num. points: 3 true_order: 2
So for the order one (or two) there should be 3 points, not 8.
- When I use mesh smaller than 7070, no problem. When 9090 very slow. When 120*120 "MemoryError" comes out. Thanks a lot. My codes and errors are as follow:
field = Field.from_args('displacement', np.float64, 'vector', omega, 1) u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u') prestress_value = np.zeros((nx*ny,3,1)) m = Material('m', lam=lambda_, mu=mu) prestress = Material('prestress', val=prestress_value)
This causes the memory error - you set the (nx*ny,3,1) as a constant value parameter to the material, so it is then repeated in pb.update_materials() below into all quadrature points in the whole mesh.
I do not recommend to work with the datas attribute directly - instead,
use a
function:
def eval_prestress(ts, coors, mode, **kwargs): """Evaluate/return the prestress here"""
m = Material('m', function=eval_prestress)
or, alternatively, use a placeholder constant above (for example prestress =
Material('prestress', val=0.0)) and then modify the datas by using
Material.set_all_data(). The former way (using a user-defined function) is
cleaner, though.
integral = Integral('i', order=3) t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u )', integral, omega, m=m , v=v, u=u) t2 = Term.new('dw_lin_prestress(prestress.val, v )',integral, omega, prestress=prestress, v=v) eq = Equation('balance', t1-t2 ) eqs = Equations([eq]) ls = ScipyDirect({}) nls_status = IndexedStruct() nls = Newton({'i_max' : 1,'eps_a' : 1e-6,'problem' : 'nonlinear'}, lin_solver=ls, status=nls_status) pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls) pb.update_materials() prestress_value = np.zeros((nx*ny,4,3,1))
###### #
this is how I assign values to prestress, eta is a known scalar field.
###### #
for i in range(nxny): for j in range(4): prestress_value[i][j] = np.array([[stress00_11], [stress00_22], [ stress00_12]]) eta.value[i] prestress.datas[('Omega', 'i')][0]['val'] = prestress_value pb.update_materials()
###### #
Let us know if the above helps. If not, send here the complete script, so that we can try it.
Cheers, r.