|
From: | SIMON AMEYE |
Subject: | [Getfem-users] Forces recovery |
Date: | Wed, 24 Jan 2018 09:32:49 +0000 |
Hi all, As I was not able to use the contact brick between two meshes, I now try to use another way: I want to compute the forces at the boundary (dirichlet condition) with my first part. And then, I want to use those forces on a second finite element analysis. Everything is ok except I am not able to associate the forces I get with the correct nodes. I use this code : % End of the fem computation gf_model_set(md,
'add generalized Dirichlet condition with multipliers', mim,
'u', mfu, 1,'VECTOR',
'H'); 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', mfdu); % Forces recovery (from previous getfem++ thread) % get the tangent matrix
tangent_matrix = gf_model_get(md,
'tangent_matrix');
% get the multipliers
mult = gf_model_get(md,
'variable',
'mult_on_u'); % mult2 = gf_model_get(md, 'variable', 'mult_on_u');
% get the number of multipliers and DOFs
nb_mult = size(mult,2);
nb_dof = gf_mesh_fem_get(mfu,'nbdof');
% part of the tangent matrix concerning the multipliers
mult_matrix = zeros(nb_dof,nb_mult);
for
i = 1:nb_dof for
j = 1:nb_mult mult_matrix(i,j) = tangent_matrix(nb_mult+i,j);
end
end
% computing the nodal forces by multiplying the multipliers
% with the right part of the tangent matrix
nodalforce = -(mult_matrix*transpose(mult));
% extract the x-, y- components
for
i = 0:(size(nodalforce,1)/2)-1 nodalforce_x(i+1) = nodalforce(2*i+1);
nodalforce_y(i+1) = nodalforce(2*i+2);
end % Associate the forces with the nodes [nodalforce_x,nodalforce_y] = GetTheForce(md,mfu,'mult_on_u'); F = [nodalforce_x' nodalforce_y']; ID = gf_mesh_get(magn_mesh,
'pid'); Coor = gf_mesh_get(magn_mesh,
'pts',ID)'; Coor = Coor(1:end,:); % Plot the forces with segments on nodes gf_plot(mfdu,VM,'deformed_mesh','on',
'deformation',U,'deformation_mf',mfu,'refine',
4, 'deformation_scale',DEFO_SCALE*100*0*mod(i,2));hold
on;
plot([Coor2(:,1)';Coor2(:,1)'+F(:,1)'*1e-4],[Coor2(:,2)';Coor2(:,2)'+F(:,2)'*1e-4]) Result : As you can see with this result, the forces are not well located as there are forces outside the boundary. (The boundary used for the dirichlet condition is the outside of the mesh) Do anyone have a solution ? Thank you again, Simon
|
[Prev in Thread] | Current Thread | [Next in Thread] |