Hello everyone,
I?m using the Matlab interface of GetFEM and I tried to build a very
simple three-dimensional beam. It is fixed on the left side. On the right
side I have two alternatives to bend it:
1. A force
2. A prescribed displacement
My question is: How can I compute the reaction force on the left side,
when I use alternative 2?
I´ve found something similar in an earlier post (
http://www.mail-archive.com/[email protected]/msg00273.html ), but I´m
not able to adapt it for my problem using the Matlab interface.
I added my matlab-code below. I?m pretty new to GetFEM and FEM at all,
hence I´m not sure if I set all the conditions correctly, especially
concerning the use of multipliers and penalization.
Thank you in advance,
Julia Guenther
%% ------------------------------- Mesh ----------------------------------
L = 10;
m = gfMesh('cartesian', 0:.5:L, 0:.5:1, 0:.5:1);
%% ------------------------------ FEM ------------------------------------
n = gf_mesh_get(m, 'dim');
k = 1; % Degree
% FEM u
f = gf_fem(sprintf('FEM_QK(%d,%d)',n,k));
mfu = gf_mesh_fem(m,3); gf_mesh_fem_set(mfu, 'fem', f);
% FEM Von Mises
f2 = gf_fem(sprintf('FEM_QK_DISCONTINUOUS(%d,1)', n));
mfvm = gf_mesh_fem(m,1); gf_mesh_fem_set(mfvm, 'fem', f2);
%% -------------------------------- IM -----------------------------------
k_i = 2*k;
INTEG_TYPE = sprintf('IM_GAUSS_PARALLELEPIPED(%d,%d)',n,k_i);
im = gf_integ(INTEG_TYPE); mim = gf_mesh_im(m,im);
%% ------------------------------ Regions --------------------------------
P = get(m,'pts');
boundary_left = gf_mesh_get(m, 'faces from pid', find(abs(P(1,:))<1e-6));
boundary_right = gf_mesh_get(m, 'faces from pid', find(abs(P(1,:) -
L)<1e-6));
gf_mesh_set(m, 'region', 1, boundary_left);
gf_mesh_set(m, 'region', 2, boundary_right);
%% ------------------------------- Data ---------------------------------
E = 200000;
nu = 0.3;
% Lamé
lambda = nu*E/((1+nu)*(1-2*nu));
mu = E/(2*(1+nu));
%% ------------------------------ 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 isotropic linearized elasticity brick', ...
mim, 'u', 'lambda', 'mu');
% homogenous Dirichlet on region 1
gf_model_set(md, 'add Dirichlet condition with penalization', ...
mim, 'u', 1e10, 1);
alt = 2; % 1 Force, 2 Displacement
% alternative 1: Force on region 2
if alt ==1
F = [0 -100 0];
gf_model_set(md, 'add initialized data', 'Force', F);
gf_model_set(md, 'add source term brick', mim, 'u', 'Force', 2);
else
% alternative 2: Displacement on region 2
diri = [0 -1.5 0];
gf_model_set(md, 'add initialized data', 'dirichletdata', diri);
gf_model_set(md, 'add Dirichlet condition with multipliers', ...
mim, 'u', mfu, 2, 'dirichletdata');
end
%% -------------------------- Solving -------------------------------
gf_model_get(md, 'solve');
u = gf_model_get(md, 'variable', 'u');
VM = gf_model_get(md, 'compute isotropic linearized Von Mises or Tresca',
'u', 'lambda', 'mu', mfvm);
f = figure;
gf_plot(mfvm,VM, 'deformation',u, 'deformation_mf',mfu, 'deformed_mesh',
'on', 'cvlst', get(m, 'outer faces'), 'deformation_scale', 1);
colorbar;
xlabel('x')
ylabel('y')
zlabel('z')
set(f, 'Units', 'normalized', 'Position', [0.2, 0.1, 0.7, 0.7]);_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
