[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r4822 - /trunk/getfem/interface/tests/matlab/demo_dynam
From: |
farshid . dabaghi |
Subject: |
[Getfem-commits] r4822 - /trunk/getfem/interface/tests/matlab/demo_dynamic_plasticity.m |
Date: |
Wed, 26 Nov 2014 11:53:27 -0000 |
Author: fdabaghi
Date: Wed Nov 26 12:53:26 2014
New Revision: 4822
URL: http://svn.gna.org/viewcvs/getfem?rev=4822&view=rev
Log:
add general alpha method
Modified:
trunk/getfem/interface/tests/matlab/demo_dynamic_plasticity.m
Modified: trunk/getfem/interface/tests/matlab/demo_dynamic_plasticity.m
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_dynamic_plasticity.m?rev=4822&r1=4821&r2=4822&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_dynamic_plasticity.m
(original)
+++ trunk/getfem/interface/tests/matlab/demo_dynamic_plasticity.m Wed Nov
26 12:53:26 2014
@@ -41,6 +41,13 @@
NX = 50;
NY = 20;
+
+% alpha is parametr of the generalized integration algorithms The
+% The choice alpha = 1/2 yields the mid point method and alpha = 1 leads to
+% backward Euler integration
+alpha = 1;
+
+
% f = [0 -600]';
%t = [0 0.5 0.6 0.7 0.8 0.9 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0];
f = [15000 0]';
@@ -139,7 +146,12 @@
% gf_model_set(md, 'add fem data', 'Previous_u', mf_u);
mim_data = gf_mesh_im_data(mim, -1, [N, N]);
gf_model_set(md, 'add im data', 'sigma', mim_data);
-
+
+
+ % Declare that alpha is a data of the system
+
+set(md, 'add initialized data', 'alpha', [alpha]);
+
set(md, 'add initialized data', 'H', [H]);
Is = 'Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim)';
@@ -150,8 +162,15 @@
ApH = sprintf('((2*(mu)+(H))*(%s) + (lambda)*(%s))', Is, IxI);
Enp1 = '((Grad_u+Grad_u'')/2)';
En = '((Grad_Previous_u+Grad_Previous_u'')/2)';
-expr_sigma = strcat('(', B_inv, '*(Von_Mises_projection((-(H)*', Enp1, ')+(',
ApH, '*(',Enp1,'-',En,')) + (', B, '*sigma), von_mises_threshold) + H*', Enp1,
'))');
-
+ %expression de sigma for Implicit Euler method
+ %expr_sigma = strcat('(', B_inv, '*(Von_Mises_projection((-(H)*', Enp1,
')+(', ApH, '*(',Enp1,'-',En,')) + (', B, '*sigma), von_mises_threshold) + H*',
Enp1, '))');
+
+ %expression de sigma for generalized alpha algorithms
+ expr_sigma = strcat('(', B_inv,
'*(Von_Mises_projection((',B,'*(1-alpha)*sigma)+(-(H)*(((1-alpha)*',En,')+(alpha*',
Enp1, ')))+(alpha*', ApH, '*(',Enp1,'-',En,')) + (alpha*', ...
+ B, '*sigma), von_mises_threshold) + (H)*(((1-alpha)*',En,')+(alpha*',
Enp1, '))))');
+
+
+
gf_model_set(md, 'add nonlinear generic assembly brick', mim,
strcat(expr_sigma, ':Grad_Test_u'));
% gf_model_set(md, 'add finite strain elasticity brick', mim, 'u',
'SaintVenant Kirchhoff', '[lambda; mu]');
@@ -164,6 +183,9 @@
set(md, 'add source term brick', mim, 'u', 'VolumicData', 2);
+
+
+
% interpolate the initial data
U0 = get(md, 'variable', 'u');
V0 = 0*U0;
@@ -176,6 +198,23 @@
% Initialisation of the acceleration 'Previous_Dot2_u'
gf_model_set(md, 'perform init time derivative', dt/20.);
gf_model_get(md, 'solve');
+
+
+
+if(alpha==0.5)
+ nbdofu = gf_mesh_fem_get(mf_u, 'nbdof');
+ M = gf_asm('mass matrix', mim, mf_u);
+ gf_model_set(md, 'add explicit matrix', 'u', 'u', M/(dt*dt*alpha*alpha));
+ ind_rhs = gf_model_set(md, 'add explicit rhs', 'u', zeros(nbdofu,1));
+ MV0=M*V0;
+ F=get(mf_data, 'eval',{f(1,1)*sin(0);f(2,1)*sin(0)});
+ FF = gf_asm('volumic source', mim, mf_u, mf_data, F);
+ K = gf_asm('elastoplasticity', mim, mf_u, mf_data,
ones(nbdofd,1)*clambda, ones(nbdofd,1)*cmu); %%%%%% ???????????
+ MA0 = FF-K*U0;
+end
+
+
+
VM=zeros(1,get(mf_vm, 'nbdof'));
@@ -184,6 +223,20 @@
for t = 0:dt:T
coeff = sin(16*t);
+
+ if(alpha==0.5)
+ F = get(mf_data, 'eval',{f(1,1)*coeff;f(2,1)*coeff});
+ FF = gf_asm('volumic source', mim, mf_u, mf_data, F);
+ D = (M*U0 + dt*MV0)/(dt*dt*alpha*alpha) + (1-alpha)*MA0/alpha;
+ gf_model_set(md, 'set private rhs', ind_rhs, D);
+ gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', niter);
+ U = gf_model_get(md, 'variable', 'u');
+ MV = ((M*U1 - M*U0)/dt -(1-alpha)*MV0)/alpha;
+ MA = ((MV1-MV0)/dt - (1-alpha)*MA0)/alpha;
+ else
+
+
+
disp(sprintf('step %d, coeff = %g', step , coeff));
@@ -194,12 +247,12 @@
U = gf_model_get(md, 'variable', 'u');
V = gf_model_get(md, 'variable', 'Dot_u');
-
+ end
if (test_tangent_matrix)
gf_model_get(md, 'test tangent matrix', 1E-8, 10, 0.000001);
- end;
+ end;
@@ -207,15 +260,15 @@
gf_model_set(md, 'variable', 'sigma', sigma);
gf_model_set(md, 'variable', 'Previous_u', U);
- M = gf_asm('mass matrix', mim, mf_vm);
+ M_VM = gf_asm('mass matrix', mim, mf_vm);
L = gf_asm('generic', mim, 1, 'sqrt(3/2)*Norm(Deviator(sigma))*Test_vm',
-1, 'sigma', 0, mim_data, sigma, 'vm', 1, mf_vm, zeros(gf_mesh_fem_get(mf_vm,
'nbdof'),1));
- VM = (M \ L)';
+ VM = (M_VM \ L)';
coeff1='-lambda/(2*mu*(meshdim*lambda+2*mu))';
coeff2='1/(2*mu)';
Ainv=sprintf('(%s)*(%s) + (%s)*(%s)', coeff1, IxI, coeff2, Is);
Ep = sprintf('(Grad_u+Grad_u'')/2 - (%s)*sigma', Ainv);
L = gf_asm('generic', mim, 1, sprintf('Norm(%s)*Test_vm', Ep), -1,
'sigma', 0, mim_data, sigma, 'u', 0, mf_u, U, 'vm', 1, mf_vm,
zeros(gf_mesh_fem_get(mf_vm, 'nbdof'),1), 'mu', 0, mf_data, mu, 'lambda', 0,
mf_data, lambda);
- plast = (M \ L)';
+ plast = (M_VM \ L)';
Epsilon_u = gf_model_get(md, 'interpolation', '((Grad_u+Grad_u'')/2)',
mim_data);
nb_gauss_pt_per_element = size(sigma, 2) / (N*N*gf_mesh_get(m, 'nbcvs'));
@@ -230,6 +283,9 @@
end
sigma_fig(1,step)=sigma(N*N*ind_gauss_pt + 1);
Epsilon_u_fig(1,step)=Epsilon_u(N*N*ind_gauss_pt + 1);
+
+
+
@@ -255,15 +311,25 @@
plot(Epsilon_u_fig, sigma_fig,'r','LineWidth',2)
xlabel('Strain');
ylabel('Stress')
- axis([-0.25 0.25 -16500 16500 ]);
+ axis([-0.3 0.3 -16500 16500 ]);
title(sprintf('step %d / %d, coeff = %g', step,size([0:dt:T],2) ,
coeff));
pause(0.1);
end
-
step= step+ 1;
- gf_model_set(md, 'shift variables for time integration');
+
+ if(alpha==0.5)
+
+ U0 = U;
+ MV0 = MV;
+ MA0 = MA;
+ else
+
+
+ gf_model_set(md, 'shift variables for time integration');
+
+ end
end;
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r4822 - /trunk/getfem/interface/tests/matlab/demo_dynamic_plasticity.m,
farshid . dabaghi <=