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/[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 > ifalt ==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 -- Yves Renard ([email protected]) 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 ---------
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
