[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r4869 - /trunk/getfem/interface/tests/matlab/demo_nitsc
From: |
Yves . Renard |
Subject: |
[Getfem-commits] r4869 - /trunk/getfem/interface/tests/matlab/demo_nitsche_contact_by_hand.m |
Date: |
Mon, 02 Mar 2015 14:43:52 -0000 |
Author: renard
Date: Mon Mar 2 15:43:52 2015
New Revision: 4869
URL: http://svn.gna.org/viewcvs/getfem?rev=4869&view=rev
Log:
small example
Added:
trunk/getfem/interface/tests/matlab/demo_nitsche_contact_by_hand.m
Added: trunk/getfem/interface/tests/matlab/demo_nitsche_contact_by_hand.m
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_nitsche_contact_by_hand.m?rev=4869&view=auto
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_nitsche_contact_by_hand.m (added)
+++ trunk/getfem/interface/tests/matlab/demo_nitsche_contact_by_hand.m Mon Mar
2 15:43:52 2015
@@ -0,0 +1,306 @@
+
+
+asize = size(who('automatic_var654777'));
+if (asize(1)) is_automatic = true; else is_automatic = false; end
+
+gf_workspace('clear all');
+
+approximation_type = 2; % 0 : Augmentend Lagrangian
+ % 1 : Nitsche (biased)
+ % 2 : Unbiased Nitsche method
+
+draw_mesh = true;
+ref_sol = 0
+max_res = 1E-8;
+ max_iter = 50; %N = 2;% 2 ou 3 dimensions
+R=0.25;
+dirichlet_val = 0;
+ f_coeff=0;
+ clambda = 1; % Lame coefficient
+cmu = 1; % Lame coefficient
+vertical_force = -0.3;
+ penalty_parameter = 0.00001;
+if (ref_sol == 0)
+ method = [-1]; % theta
+ gamma00 = [0.1]; %1/200
+ nxy =[30]; % 2D ->400 and 3D -> 30
+ ls_degree = 2;
+
+
+end
+if (ref_sol == 1)
+ method = [0];
+ gamma00 = [1/200];
+ if (N==2)
+ nxy=[2 5 10 15 27 37 50 60 75 80];
+ else
+ nxy = [5 10 12 15 20 25];
+ end
+ ls_degree = 2;
+
+
+end
+if (ref_sol == 2)
+ method = [0 1 -1];
+ gamma00 = [400 200 100 50 25 10 1 1/10 1/25 1/50 1/100 1/200 1/400];
+ nxy = [21];
+ ls_degree = 2;
+
+
+end
+for xx = 1:1:size(method,2)
+for yy = 1:1:size(gamma00,2)
+for zz = 1:1:size(nxy,2)
+theta = method(xx);
+gamma0 = gamma00(yy);
+NX = 4*floor(nxy(zz)/4);
+
+mo = gf_mesher_object('ball',[0 0],0.25)
+mesh1 = gf_mesh('generate', mo, 1/NX ,2) ;
+%mesh1 = gf_mesh('load',
'../../../tests/meshes/sphere_with_quadratic_tetra_400_elts.mesh');
+mesh2 = gf_mesh('regular simplices', -.5:(1/NX):.5, -0.5:(1/NX):-0.25);
+
+
+N = gf_mesh_get(mesh1, 'dim');
+
+mfu1 = gf_mesh_fem(mesh1, N); gf_mesh_fem_set(mfu1, 'classical fem', 2);
+mflambda1 = gf_mesh_fem(mesh1, 1); gf_mesh_fem_set(mflambda1, 'classical fem',
1);
+
+mfvm1 = gf_mesh_fem(mesh1); gf_mesh_fem_set(mfvm1, 'classical discontinuous
fem', 1);
+
+mfu2 = gf_mesh_fem(mesh2, N); gf_mesh_fem_set(mfu2, 'classical fem', 2);
+
+mfvm2 = gf_mesh_fem(mesh2); gf_mesh_fem_set(mfvm2, 'classical discontinuous
fem', 1);
+
+CONTACT_BOUNDARY1 = 1;
+
+fb1 = gf_mesh_get(mesh1, 'outer faces');
+normals1 = gf_mesh_get(mesh1, 'normal of faces', fb1);
+contact_boundary1=fb1(:,find(normals1(N,:) <= 2/NX));
+gf_mesh_set(mesh1,'region', CONTACT_BOUNDARY1,contact_boundary1);
+
+mim1 = gf_mesh_im(mesh1, 4);
+mim1_contact = gf_mesh_im(mesh1, 6);
+
+CONTACT_BOUNDARY2 = 2;
+
+border = gf_mesh_get(mesh2,'outer faces');
+ normals2 = gf_mesh_get(mesh2, 'normal of faces', border);
+
+ P = gf_mesh_get(mesh2, 'pts');
+ P=P(:,find(P(1,:)>=-0.25));
+ P=P(:,find(P(1,:)<=0.25)) ;
+ P=P(:,find(P(2,:)>=-0.25)) ;
+ Id=gf_mesh_get(mesh2, 'pid from coords', P);
+ contact_boundary2=gf_mesh_get(mesh2,'faces from pid',Id);
+
+ gf_mesh_set(mesh2, 'region', CONTACT_BOUNDARY2, contact_boundary2);
+
+
+dirichlet_boundary=border(:, find(normals2(N, :) < -0.01));
+DIRICHLET_BOUNDARY2 = 3;
+gf_mesh_set(mesh2, 'region', DIRICHLET_BOUNDARY2, dirichlet_boundary);
+
+mim2 = gf_mesh_im(mesh2, 4);
+mim2_contact = gf_mesh_im(mesh2, 6);
+
+ if (draw_mesh)
+ gf_plot_mesh(mesh2);
+ hold on
+ gf_plot_mesh(mesh1, 'regions', CONTACT_BOUNDARY1);
+ hold off
+ % pause;
+ end
+
+ %Elastic model
+md=gf_model('real');
+gf_model_set(md,'add fem variable', 'u1', mfu1);
+gf_model_set(md,'add fem variable', 'u2', mfu2);
+
+gf_model_set(md,'add initialized data', 'gamma0', gamma0);
+gf_model_set(md, 'add initialized data', 'theta', theta);
+gf_model_set(md, 'add initialized data', 'friction_coeff', [f_coeff]);
+
+gf_model_set(md, 'add initialized data', 'cmu', [cmu]);
+gf_model_set(md, 'add initialized data', 'clambda', [clambda]);
+gf_model_set(md, 'add isotropic linearized elasticity brick', mim1,
'u1','clambda', 'cmu');
+gf_model_set(md, 'add isotropic linearized elasticity brick', mim2,
'u2','clambda', 'cmu');
+
+gf_model_set(md, 'add initialized data', 'Fdata', [0 vertical_force]);
+gf_model_set(md, 'add source term brick', mim1, 'u1', 'Fdata');
+
+Ddata = zeros(1, N);
+gf_model_set(md, 'add initialized data', 'Ddata2', Ddata);
+gf_model_set(md, 'add Dirichlet condition with multipliers', mim2, 'u2', 1,
DIRICHLET_BOUNDARY2, 'Ddata2');
+
+if (N==2)
+ cpoints = [0,0 0,0.1]; % constrained points for 2d
+ cunitv = [1,0 1,0]; % corresponding constrained directions for 2d,
mieux avec [0, 0.1]
+else
+ cpoints = [0, 0, 0, 0, 0, 0, 0, 0, 0.1]; % constrained points for 3d
+ cunitv = [1, 0, 0, 0, 1, 0, 0, 1, 0]; % corresponding constrained
directions for 3d, mieux avec [0, 0.1]
+end
+
+gf_model_set(md, 'add initialized data', 'cpoints', cpoints);
+gf_model_set(md, 'add initialized data', 'cunitv', cunitv);
+gf_model_set(md, 'add pointwise constraints with multipliers', 'u1',
'cpoints', 'cunitv');
+
+gf_model_set(md, 'add initialized data', 'penalty_param1',
[penalty_parameter]);
+indmass = gf_model_set(md, 'add mass brick', mim1, 'u1', 'penalty_param1');
+
+
+
+if (approximation_type == 0)
+ gf_model_set(md, 'add initialized data', 'N1', [0;-1]);
+ gf_model_set(md, 'add_interpolate_transformation_from_expression', 'Proj1',
mesh1, mesh2, '[X(1);-0.25]');
+
+ gf_model_set(md, 'add filtered fem variable', 'lambda1', mflambda1,
CONTACT_BOUNDARY1);
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact,
'-lambda1*(Test_u1.N1)+lambda1*(Interpolate(Test_u2,Proj1).N1)',
CONTACT_BOUNDARY1);
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact,
'-(gamma0*element_size)*(lambda1 +
neg_part(lambda1-(1/(gamma0*element_size))*((u1-Interpolate(u2,Proj1)+X-Interpolate(X,Proj1)).N1)))*Test_lambda1',
CONTACT_BOUNDARY1);
+
+else
+
+
+ gf_model_set(md, 'add initialized data', 'N1', [0;-1]);
+ gf_model_set(md, 'add_interpolate_transformation_from_expression', 'Proj1',
mesh1, mesh2, '[X(1);-0.25]');
+ gamma='(gamma0*element_size)';
+ sigma_vh1_n ='(((clambda*Trace(Grad_Test_u1)*Id(qdim(u1)) +
cmu*(Grad_Test_u1 + Grad_Test_u1''))*Normal).N1)';
+ sigma_uh1_n ='(((clambda*Trace(Grad_u1)*Id(qdim(u1)) + cmu*(Grad_u1 +
Grad_u1''))*Normal).N1)';
+ Pn_gamma_u1=strcat('((u1-Interpolate(u2,Proj1)+X-Interpolate(X,Proj1)).N1) -
',gamma,'*(',sigma_uh1_n,')');
+
+ if (approximation_type == 2) coeff = '0.5*'; else coeff = ''; end
+
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact,
strcat('-',coeff,'theta*',gamma,'*(',sigma_uh1_n,')*(',sigma_vh1_n,')'),
CONTACT_BOUNDARY1);
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact,
strcat('-',coeff,'theta*',sigma_vh1_n,'*pos_part(',Pn_gamma_u1,')'),
CONTACT_BOUNDARY1);
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact,
strcat(coeff,'(1/',gamma,')*(Test_u1.N1)*pos_part(',Pn_gamma_u1,')'),
CONTACT_BOUNDARY1);
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact,
strcat('-',coeff,'(1/',gamma,')*(Interpolate(Test_u2,Proj1).N1)*pos_part(',Pn_gamma_u1,')'),CONTACT_BOUNDARY1);
+
+ if (approximation_type == 2)
+ gf_model_set(md, 'add initialized data', 'N2', [0;1]);
+ gf_model_set(md,'add_interpolate_transformation_from_expression', 'Proj2',
mesh2,mesh1, '[ X(1) ; 0.0001-sqrt(0.0625-sqr(X(1)))]');
+ sigma_vh2_n ='(((clambda*Trace(Grad_Test_u2)*Id(qdim(u2)) +
cmu*(Grad_Test_u2 + Grad_Test_u2'')).Normal).N2)';
+ sigma_uh2_n ='(((clambda*Trace(Grad_u2)*Id(qdim(u2)) + cmu*(Grad_u2 +
Grad_u2'')).Normal).N2)';
+ Pn_gamma_u2=strcat('((u2-Interpolate(u1,Proj2)).N2) -
',gamma,'*(',sigma_uh2_n,')-(Interpolate(X,Proj2)-X).N2');
+
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim2_contact,
strcat('-0.5*theta*',gamma,'*(',sigma_uh2_n,')*(',sigma_vh2_n,')'),
CONTACT_BOUNDARY2);
+
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim2_contact,
strcat('-0.5*theta*',sigma_vh2_n,'*pos_part(',Pn_gamma_u2,')'),
CONTACT_BOUNDARY2);
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim2_contact,
strcat('0.5*(1/',gamma,')*(Test_u2.N2)*pos_part(',Pn_gamma_u2,')'),
CONTACT_BOUNDARY2);
+ gf_model_set(md, 'add nonlinear generic assembly brick', mim2_contact,
strcat('-0.5*(1/',gamma,')*(Interpolate(Test_u1,Proj2).N2)*pos_part(',Pn_gamma_u2,')'),
CONTACT_BOUNDARY2);
+ end
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+solve=true;
+% gf_model_get(md, 'test tangent matrix term', 'u1', 'u2', 1e-6, niter, 10.0);
+% gf_model_get(md, 'test tangent matrix', 1e-6, 10, 10);
+
+niter= 100;
+gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', niter, 'noisy',
'lsearch', 'simplest', 'alpha min', 0.8);
+
+U1 = gf_model_get(md, 'variable', 'u1');
+UU1 = gf_model_get(md, 'variable', 'u1');
+VM1 = gf_model_get(md, 'compute_isotropic_linearized_Von_Mises_or_Tresca', ...
+ 'u1', 'clambda', 'cmu', mfvm1);
+U2 = gf_model_get(md, 'variable', 'u2');
+UU2 = gf_model_get(md, 'variable', 'u2');
+VM2 = gf_model_get(md, 'compute_isotropic_linearized_Von_Mises_or_Tresca', ...
+ 'u2', 'clambda', 'cmu', mfvm2);
+
+% plot figure
+ if (ref_sol==1)
+ sl1=gf_slice({'none'}, mesh1, 5);
+ P1=gf_slice_get(sl1,'pts'); dP1=gf_compute(mfu1,U1,'interpolate
on',sl1);
+ gf_slice_set(sl1, 'pts', P1 + dP1);
+ VMsl1=gf_compute(mfvm1,VM1,'interpolate on',sl1);
+ sl2=gf_slice({'none'}, mesh2, 5);
+ P2=gf_slice_get(sl2,'pts'); dP2=gf_compute(mfu2,U2,'interpolate
on',sl2);
+ gf_slice_set(sl2, 'pts', P2+dP2);
+ VMsl2=gf_compute(mfvm2,VM2,'interpolate on',sl2);
+ gf_plot_slice(sl1,'mesh','off','mesh_slice_edges','off','data',VMsl1);
+ hold on;
+ gf_plot_slice(sl2,'mesh','off','mesh_slice_edges','off','data',VMsl2);
+ hold off;
+ %else
+
sl1=gf_slice({'boundary',{'intersection',{'ball',-1,[0;0;0],R},{'planar',1,[0;0;0],[1;0;0]}}},
m, 5);
+
sl2=gf_slice({'boundary',{'intersection',{'planar',-1,[0;0;zc],[0;0;1]},{'planar',1,[0;0;0],[1;0;0]}}},
m, 5);
+ P1=gf_slice_get(sl1,'pts');P2=gf_slice_get(sl2,'pts');
+ dP1=gf_compute(mfu1,U1,'interpolate on',sl1);
+ dP2=gf_compute(mfu2,U2,'interpolate on',sl2);
+ gf_slice_set(sl1, 'pts', P1+dP1);
+ gf_slice_set(sl2, 'pts', P2+dP2);
+ VMsl1=gf_compute(mfvm,VM1,'interpolate on',sl1);
+ VMsl2=gf_compute(mfvm,VM2,'interpolate on',sl2);
+ set(gcf,'renderer','zbuffer');
+
h=gf_plot_slice(sl1,'mesh','off','mesh_faces','on','mesh_slice_edges','on','data',VMsl1);
+ hold on;
+
h=gf_plot_slice(sl2,'mesh','off','mesh_slice_edges','off','data',VMsl2);
+ hold off;
+ view(-55,10); axis on; camlight('headlight');
+ xlabel('x'); ylabel('y'); zlabel('z');
+ %title('3D');
+
+ end
+
+gf_plot(mfvm1,VM1,'mesh', 'off', 'deformed_mesh','on',
'deformation',U1,'deformation_mf',mfu1,'deformation_scale', 1, 'refine', 8);
colorbar;
+hold on
+gf_plot(mfvm2,VM2,'mesh', 'off', 'deformed_mesh','on',
'deformation',U2,'deformation_mf',mfu2,'deformation_scale', 1, 'refine', 8);
colorbar;
+
+return
+
+ if (ref_sol == 0)
+ map(1,:)=[0 200 255]; % bleu fonc??
+ map(2,:)=[0 255 255];
+ map(3,:)=[128 255 128];
+ map(4,:)=[255 255 0];
+ map(5,:)=[255 128 0];
+ map(6,:)=[255 0 0];
+ colormap(map./255)
+ gf_mesh_fem_get(mfu1, 'save', 'sol_ref_mesh_fem1','with_mesh');
+ gf_mesh_fem_get(mfu2, 'save', 'sol_ref_mesh_fem2','with_mesh');
+ save sol_de_reference1 UU1;
+ save sol_de_reference2 UU2;
+ else
+ meshref = gf_mesh('load', 'sol_ref_mesh_fem');
+ end
+end
+end
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r4869 - /trunk/getfem/interface/tests/matlab/demo_nitsche_contact_by_hand.m,
Yves . Renard <=