[Getfem-users] compute a reaction force

[Julia Guenther]

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/address@hidden ), 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]);