Re: [Getfem-users] compute a reaction force

[Yves Renard]

Dear Julia Guenther,

There is at leat two means to obtain the force density on a Dirichlet boundary. The first one is to prescribe the Dirichlet condition with a multiplier. The value of the multiplier directly gives the force which has been necessary to prescribe the Dirichlet condition (but of course this adds some unknowns). The second mean is to compute the residual of the problem without this dirichlet condition. this means in that case to disable the corresponding brick and ask for the assembly and then compute the residual. But in this latter case, this will give the result on the same finite element method as the one used for the displacement. So you have to take the value on the boundary of the result. It could be simpler in some situations but may be not in your case.

Yves.



Le 17/12/2013 15:42, Julia Guenther a écrit :

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]);

_______________________________________________
Getfem-users mailing list
address@hidden
https://mail.gna.org/listinfo/getfem-users

-- 

  Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------