Dear Fabien Amiot,
This is compeltely normal that you have a vanishing rhs since there is no load
on your model (volumic or surface forces). The incompressibility brick
introduces a new unknown (the pressure) and term to ensure that the
displacement is (linearly) incompressible. In fact, I did not understand well
what you intend to do.
Yves.
----- Original Message -----
From: "fabien amiot"<address@hidden>
To: address@hidden
Sent: Thursday, October 4, 2012 4:02:20 PM
Subject: [Getfem-users] eigenstrain
Dear all,
I'm rather new to Getfem (Matlab-interface), and I'm trying to simulate
a (linear elastic) bimorph :
%% DATA
% cantilever length
L=1;
% cantilever thickness
H=0.1;
% thickness ratio
r=0.1;
% Lame parameters
Young=60E9; nu=0.3;
lambda=Young.*nu./((1+nu).*(1-2.*nu));
mu=0.5.*Young./(1+nu);
lambda_layer=lambda;
mu_layer=mu;
%% Build the mesh
% cantilever
canti_m=gf_mesh('regular simplices',0:(L./10):L, (-H):(H./4):0);
% layer
layer_m=gf_mesh('regular simplices',0:(L./10):L, 0:(r.*H./4):(r.*H));
% retrieve the CVFIDs for the 2 regions
canti_cvfids = gf_mesh_get(canti_m, 'region',gf_mesh_get(canti_m,
'regions'));
layer_cvfids = gf_mesh_get(layer_m, 'region',gf_mesh_get(layer_m,
'regions'));
% merge the two meshes
gf_mesh_set(canti_m, 'merge', layer_m);
% set the region numbers
gf_mesh_set(canti_m, 'region', 1, canti_cvfids);
gf_mesh_set(canti_m, 'region', 2, layer_cvfids);
% integration method
mim=gf_mesh_im(canti_m);
gf_mesh_im_set(mim, 'integ',gf_integ('IM_TRIANGLE(6)'));
% fem object for displacement
mfu=gf_mesh_fem(canti_m,2);
gf_mesh_fem_set(mfu, 'fem',gf_fem('FEM_PK(2,2)'));
% fem object for pressure (used to impose the free strain as a
incompressibility condition)
mfp=gf_mesh_fem(canti_m);
gf_mesh_fem_set(mfp, 'fem',gf_fem('FEM_PK(2,1)'));
%% build elastic model
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mfu);
gf_model_set(md, 'add initialized data', 'lambda', [lambda]);
gf_model_set(md, 'add initialized data', 'mu', [mu]);
gf_model_set(md, 'add initialized data', 'lambda_layer', [lambda_layer]);
gf_model_set(md, 'add initialized data', 'mu_layer', [mu_layer]);
% set the elastic parameters for the cantilever
gf_model_set(md, 'add isotropic linearized elasticity brick', ...
mim, 'u', 'lambda', 'mu',1);
% set the elastic parameters for the layer
gf_model_set(md, 'add isotropic linearized elasticity brick', ...
mim, 'u', 'lambda_layer', 'mu_layer',2);
% set the boundary condition (clamping)
P=gf_mesh_get(canti_m,'pts');
fleft =gf_mesh_get(canti_m,'faces from pid',find(abs(P(1,:))<1e-6));
gf_mesh_set(canti_m,'boundary',3,fleft);
gf_model_set(md, 'add Dirichlet condition with multipliers', ...
mim, 'u', mfu, 3);
The idea is to describe the eigenstrain by exploiting the linear
incompressibility brick :
% set the loading
gf_model_set(md, 'add fem variable', 'p', mfp);
gf_model_set(md, 'add initialized data', 'coef', 1.);
gf_model_set(md, 'add linear incompressibility brick', mim, 'u', 'p', 2,
'coef');
But this produces a zero rhs (as seen from gf_model_get(md,'rhs') ). I
guess I miss either an important thing or an efficient alternative way
to achieve this...
Best regards,
Fabien