[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r4800 - in /trunk/getfem: interface/tests/matlab/ src/
From: |
Yves . Renard |
Subject: |
[Getfem-commits] r4800 - in /trunk/getfem: interface/tests/matlab/ src/ |
Date: |
Wed, 05 Nov 2014 16:01:58 -0000 |
Author: renard
Date: Wed Nov 5 17:01:57 2014
New Revision: 4800
URL: http://svn.gna.org/viewcvs/getfem?rev=4800&view=rev
Log:
minor corrections
Modified:
trunk/getfem/interface/tests/matlab/demo_laplacian.m
trunk/getfem/interface/tests/matlab/demo_plasticity.m
trunk/getfem/src/getfem_plasticity.cc
Modified: trunk/getfem/interface/tests/matlab/demo_laplacian.m
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_laplacian.m?rev=4800&r1=4799&r2=4800&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_laplacian.m (original)
+++ trunk/getfem/interface/tests/matlab/demo_laplacian.m Wed Nov 5
17:01:57 2014
@@ -21,6 +21,8 @@
gamma0 = 0.001; % Nitsche's method parameter gamma0 (gamma = gamma0*h)
r = 1e8; % Penalization parameter
draw = true;
+quadrangles = true;
+K = 2; % Degree of the finite element method
asize = size(who('automatic_var654'));
if (asize(1)) draw = false; end;
@@ -28,17 +30,24 @@
% trace on;
gf_workspace('clear all');
NX = 20;
-m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
-%m=gf_mesh('import','structured','GT="GT_QK(2,1)";SIZES=[1,1];NOISED=1;NSUBDIV=[1,1];')
+if (quadrangles)
+ m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
+else
+
m=gf_mesh('import','structured',sprintf('GT="GT_PK(2,1)";SIZES=[1,1];NOISED=0;NSUBDIV=[%d,%d];',
NX, NX))
+end
% create a mesh_fem of for a field of dimension 1 (i.e. a scalar field)
mf = gf_mesh_fem(m,1);
-% assign the Q2 fem to all convexes of the mesh_fem,
-gf_mesh_fem_set(mf,'fem',gf_fem('FEM_QK(2,2)'));
+% assign the QK or PK fem to all convexes of the mesh_fem, and define an
+% integration method
+if (quadrangles)
+ gf_mesh_fem_set(mf,'fem',gf_fem(sprintf('FEM_QK(2,%d)', K)));
+ mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,6)'));
+else
+ gf_mesh_fem_set(mf,'fem',gf_fem(sprintf('FEM_PK(2,%d)', K)));
+ mim = gf_mesh_im(m, gf_integ('IM_TRIANGLE(6)'));
+end
-% Integration which will be used
-mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,4)'));
-%mim = gf_mesh_im(m,
gf_integ('IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,5),4)'));
% detect the border of the mesh
border = gf_mesh_get(m,'outer faces');
% mark it as boundary #1
@@ -49,9 +58,11 @@
end
% interpolate the exact solution
-Uexact = gf_mesh_fem_get(mf, 'eval', { 'y.*(y-1).*x.*(x-1)+x.^5' });
-% its second derivative
-F = gf_mesh_fem_get(mf, 'eval', { '-(2*(x.^2+y.^2)-2*x-2*y+20*x.^3)' });
+% Uexact = gf_mesh_fem_get(mf, 'eval', { '10*y.*(y-1).*x.*(x-1)+10*x.^5' });
+Uexact = gf_mesh_fem_get(mf, 'eval', { 'x.*sin(2*pi*x).*sin(2*pi*y)' });
+% Opposite of its laplacian
+% F = gf_mesh_fem_get(mf, 'eval', {
'-(20*(x.^2+y.^2)-20*x-20*y+200*x.^3)' });
+F = gf_mesh_fem_get(mf, 'eval', { '4*pi*(2*pi*x.*sin(2*pi*x) -
cos(2*pi*x)).*sin(2*pi*y)' });
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mf);
@@ -74,18 +85,18 @@
U = gf_model_get(md, 'variable', 'u');
if (draw)
- subplot(2,1,1); gf_plot(mf,U,'mesh','on','contour',.01:.01:.1);
+ subplot(2,1,1); gf_plot(mf,U,'mesh','off');
colorbar; title('computed solution');
- subplot(2,1,2); gf_plot(mf,U-Uexact,'mesh','on');
+ subplot(2,1,2); gf_plot(mf,Uexact-U,'mesh','on');
colorbar;title('difference with exact solution');
end
-err = gf_compute(mf, U-Uexact, 'H1 norm', mim);
+err = gf_compute(mf, Uexact-U, 'H1 norm', mim);
disp(sprintf('H1 norm of error: %g', err));
-if (err > 2E-5)
+if (err > 2E-4)
error('Laplacian test: error to big');
end
Modified: trunk/getfem/interface/tests/matlab/demo_plasticity.m
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_plasticity.m?rev=4800&r1=4799&r2=4800&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_plasticity.m (original)
+++ trunk/getfem/interface/tests/matlab/demo_plasticity.m Wed Nov 5
17:01:57 2014
@@ -28,14 +28,14 @@
%%
new_test = 0;
-test_tangent_matrix = 1;
+test_tangent_matrix = 0;
% Initialize used data
L = 100;
H = 20;
% f = [0 -330]';
-f = [10000 0]';
+f = [1 0]';
t = [0 0.9032 1 1.1 1.3 1.5 1.7 1.74 1.7 1.5 1.3 1.1 1 0.9032 0.7 0.5 0.3 0.1
0];
% Create the mesh
@@ -70,14 +70,20 @@
CVtop = find(CV > 250); % Retrieve index of convex located at the top
% Definition of Lame coeff
-lambda(CVbottom,1) = 121150; % Stell
-lambda(CVtop,1) = 84605; % Iron
+lambda(CVbottom,1) = 12.1150; % Steel
+lambda(CVtop,1) = 8.4605; % Iron
%lambda = 121150;
-mu(CVbottom,1) = 80769; %Stell
-mu(CVtop,1) = 77839; % Iron
+mu(CVbottom,1) = 8.0769; %Steel
+mu(CVtop,1) = 7.7839; % Iron
%mu = 80769;
-von_mises_threshold(CVbottom) = 7000;
-von_mises_threshold(CVtop) = 8000;
+von_mises_threshold(CVbottom) = 0.7000;
+von_mises_threshold(CVtop) = 0.8000;
+
+if (new_test == 1)
+ lambda(CVtop,1) = 12.1150;
+ mu(CVtop,1) = 8.0769;
+ von_mises_threshold(CVtop) = 0.7000;
+end;
% Assign boundary numbers
set(m,'boundary',1,fleft); % for Dirichlet condition
@@ -110,16 +116,18 @@
H = mu(1)/2;
set(md, 'add initialized data', 'H', [H]);
- coeff_long = '(2*lambda*H)/((2*mu+H)*(2*meshdim*lambda+4*mu+2*H))';
- B_inv =
sprintf('((2*mu/(2*mu+H))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim)
+ (%s)*(Id(meshdim)@Id(meshdim)))', coeff_long);
- B =
'((1+H/(2*mu))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim) +
(-lambda*H/(2*mu*(meshdim*lambda+2*mu)))*(Id(meshdim)@Id(meshdim)))';
- ApH =
'((2*mu+H)*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim) +
(lambda)*(Id(meshdim)@Id(meshdim)))';
+ coeff_long = '((lambda)*(H))/((2*(mu)+(H))*(meshdim*(lambda)+2*(mu)+(H)))';
+ B_inv =
sprintf('((2*(mu)/(2*(mu)+(H)))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim)
+ (%s)*(Id(meshdim)@Id(meshdim)))', coeff_long);
+ B =
'((1+(H)/(2*(mu)))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim)
+ (-(lambda)*(H)/(2*(mu)*(meshdim*lambda+2*(mu))))*(Id(meshdim)@Id(meshdim)))';
+ ApH =
'((2*(mu)+(H))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim) +
(lambda)*(Id(meshdim)@Id(meshdim)))';
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,
'))');
+ expr_sigma = strcat('(', B_inv, '*(Von_Mises_projection(((H)*', Enp1, ')+(',
ApH, '*(',Enp1,'+',En,')) + (', B, '*sigma), von_mises_threshold) + H*', Enp1,
'))');
+
+ % expr_sigma2 = 'Von_Mises_projection(Grad_u+Grad_u'', von_mises_threshold)';
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]');
else
% Declare that sigma is a data of the system on mf_sigma
@@ -148,9 +156,9 @@
end
if (test_tangent_matrix)
- gf_model_get(md, 'test tangent matrix', 1E-8, 10, 0.0001);
+ gf_model_get(md, 'test tangent matrix', 1E-8, 10, 0.00000001);
end;
-
+
% Solve the system
get(md, 'solve','lsolver', 'superlu', 'lsearch', 'simplest', 'alpha min',
0.8, 'very noisy', 'max_iter', 100, 'max_res', 1e-6);
Modified: trunk/getfem/src/getfem_plasticity.cc
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_plasticity.cc?rev=4800&r1=4799&r2=4800&view=diff
==============================================================================
--- trunk/getfem/src/getfem_plasticity.cc (original)
+++ trunk/getfem/src/getfem_plasticity.cc Wed Nov 5 17:01:57 2014
@@ -885,16 +885,14 @@
base_matrix tau(N, N), tau_D(N, N);
gmm::copy(args[0]->as_vector(), tau.as_vector());
scalar_type s = (*(args[1]))[0];
-
-
+
scalar_type tau_m = gmm::mat_trace(tau) / scalar_type(N);
gmm::copy(tau, tau_D);
for (size_type i = 0; i < N; ++i) tau_D(i,i) -= tau_m;
scalar_type norm_tau_D = gmm::mat_euclidean_norm(tau_D);
- if (norm_tau_D != scalar_type(0))
- gmm::scale(tau_D, std::min(norm_tau_D, s) / norm_tau_D);
+ if (norm_tau_D > s) gmm::scale(tau_D, s / norm_tau_D);
for (size_type i = 0; i < N; ++i) tau_D(i,i) += tau_m;
@@ -918,23 +916,21 @@
switch(nder) {
case 1:
- if (norm_tau_D < s) {
+ if (norm_tau_D <= s) {
gmm::clear(result.as_vector());
for (size_type i = 0; i < N; ++i)
for (size_type j = 0; j < N; ++j)
result(i,j,i,j) = scalar_type(1);
} else {
- if (norm_tau_D != scalar_type(0)) {
- for (size_type i = 0; i < N; ++i)
- for (size_type j = 0; j < N; ++j)
- for (size_type m = 0; m < N; ++m)
- for (size_type n = 0; n < N; ++n)
- result(i,j,m,n)
- = s * (-tau_D(i,j) * tau_D(m,n)
- + (i == m && j == n) ? scalar_type(1) :
scalar_type(0)
- - (i == j && m == n) ?
scalar_type(1)/scalar_type(N)
- : scalar_type(0)) /
norm_tau_D;
- } else gmm::clear(result.as_vector());
+ for (size_type i = 0; i < N; ++i)
+ for (size_type j = 0; j < N; ++j)
+ for (size_type m = 0; m < N; ++m)
+ for (size_type n = 0; n < N; ++n)
+ result(i,j,m,n)
+ = s * (-tau_D(i,j) * tau_D(m,n)
+ + ((i == m && j == n) ? scalar_type(1) :
scalar_type(0))
+ - ((i == j && m == n) ?
scalar_type(1)/scalar_type(N)
+ : scalar_type(0))) / norm_tau_D;
for (size_type i = 0; i < N; ++i)
for (size_type j = 0; j < N; ++j)
result(i,i,j,j) += scalar_type(1)/scalar_type(N);
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r4800 - in /trunk/getfem: interface/tests/matlab/ src/,
Yves . Renard <=