Dear all, I finally found thanks to Yves. In my case, the Dirichlet condition matrix matrix needs to validate the condition : H^2=k*H If not, the sliding is blocked. H = [x1 x2; x3 x4] Moreover, [x1 x2] needs to be a vector orthogonal with the sliding direction of the boundary with the Dirichlet condition. Here is a way I set my H matrix, with 'Angle' the angle between the sliding direction and [1 0].
x2 = -1; x1 =tan(Angle); x4 = (x2./sqrt(x1))^2; x3 = x2; Thank you again for your help, Simon De : SIMON AMEYE - U510180 Envoyé : mercredi 31 octobre 2018 10:58 À : '[email protected]' <[email protected]> Objet : RE: [Getfem-users] Sliding angle for Dirichlet condition Dear Yves, Thank you for your quick answer. I tried H = [-2 1;0 0] but when the condition I found is not respected, the boundary is completely blocked. It only slides when this condition is respected : A^2=k*A. For example : A = 10.9899 -3.3151 -3.3151 1.0000 If this property is not true, the boundary is blocked. I came to this conclusion trying a lot of combinations. The problem is that I can't find the relation between the matrix definition and the sliding vector. I want to precise that I use GetFem without having installed the full version (basic installation). I attach my code below. Simon rot_mesh = gfMesh('from string', StringMesh); %% Initialisation gf_workspace('clear all') % Numerical parameters global E_rotor Nu_rotor E = E_rotor; Nu = Nu_rotor; lambda = E*Nu/((1+Nu)*(1-2*Nu)); mu = E/(2*(1+Nu)); %% Declaration of the mesh to GetFem rot_mesh = gfMesh('from string', StringMesh); %% Setting the var description mim=gfMeshIm(rot_mesh); set(mim, 'integ',gfInteg('IM_TRIANGLE(6)')); mfu=gfMeshFem(rot_mesh,2); set(mfu, 'fem',gfFem('FEM_PK(2,1)')); mfd=gfMeshFem(rot_mesh); set(mfd, 'fem',gfFem('FEM_PK(2,1)')); mf0=gfMeshFem(rot_mesh); set(mf0, 'fem',gfFem('FEM_PK(2,0)')); mfdu=gfMeshFem(rot_mesh); set(mfdu, 'fem',gfFem('FEM_PK_DISCONTINUOUS(2,1)')); %% Boundary calculation P=get(rot_mesh, 'pts'); pidlow=find(abs(P(2,:))<1e-6); %the bottom boundary can also be found using 'LowEdgeNodes' var flow =get(rot_mesh,'faces from pid',pidlow); ftop=get(rot_mesh,'faces from pid',(TopEdgeNodes'+1)); %% Assign boundary numbers LOW_BOUND = 1; TOP_BOUND = 2; rot_mesh.set_region(LOW_BOUND, flow); rot_mesh.set_region(TOP_BOUND, ftop); %% 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'); %% Centrifugal force source term % 2D Centrifugal Force Model : F = Rho * W^2 * Distance_with_center .*[x or y].*Projection_On_[x or y] FXX = get(mfd, 'eval', {[num2str(Rho_rotor),'.*',num2str(W),'.^2.*(x.^2+y.^2).^0.5 .*x./((x.^2+y.^2).^0.5)']}); FYY = get(mfd, 'eval', {[num2str(Rho_rotor),'.*',num2str(W),'.^2.*(x.^2+y.^2).^0.5 .*y./((x.^2+y.^2).^0.5)']}); gf_model_set(md, 'add initialized fem data', 'VolumicData', mfd, [FXX;FYY]); gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData'); % Referencing cource terms to the right number gf_model_set(md, 'add initialized data', ['VolumicData' num2str(Init_For_Mag_Term+i)], [Magn(i).PressureX,Magn(i).PressureY]); gf_model_set(md, 'add source term brick', mim, 'u', ['VolumicData' num2str(Init_For_Mag_Term+i)],Init_For_Mag_Term+i); %% Sliding dirichlet conditions % H matrix building : H_top is a projector (Hat) mathix : H^2=k*h and must be symetrical % The H matrix rules the sliding angle of the boundary x4 = 1;x1 =TopAngle;k = x4+x1;x2 = sqrt(x1)*(-sqrt(k-x1));x3 = x2; gf_model_set(md, 'add initialized data', 'H_LOW', [0 0;0 1]); gf_model_set(md, 'add initialized data', 'VECTOR_LOW', [0;0]); gf_model_set(md, 'add initialized data', 'H_TOP', [x1 x2;x3 x4]); gf_model_set(md, 'add initialized data', 'VECTOR_TOP', [0;0]); gf_model_set(md, 'add generalized Dirichlet condition with multipliers', mim, 'u', mfu, LOW_BOUND,'VECTOR_LOW', 'H_LOW'); gf_model_set(md, 'add generalized Dirichlet condition with multipliers', mim, 'u', mfu, TOP_BOUND,'VECTOR_TOP', 'H_TOP'); %% Solver 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); Max_VM = max(VM); %Pa Max_Displacement = max(U); %m -----Message d'origine----- De : Yves Renard <[email protected]<mailto:[email protected]>> Envoyé : mercredi 31 octobre 2018 09:43 À : SIMON AMEYE - U510180 <[email protected]<mailto:[email protected]>> Cc : getfem-users <[email protected]<mailto:[email protected]>> Objet : Re: [Getfem-users] Sliding angle for Dirichlet condition >>> Real sender address / Reelle adresse d expedition : >>> [email protected]<mailto:[email protected]> <<< ********************************************************************** Dear Simon, A priori, H is not necessarilly a hat matrix. In that case, you can juste take an orthognal vector, for instance [-2 1] and set H = [-2 1;0 0]], this should work, I think. Best regards, Yves ----- Original Message ----- From: "SIMON AMEYE" <[email protected]<mailto:[email protected]>> To: "getfem-users" <[email protected]<mailto:[email protected]>> Sent: Tuesday, October 30, 2018 1:52:54 PM Subject: [Getfem-users] Sliding angle for Dirichlet condition C1-Non sensitive ________________________________ Hi all, It has been a long time I have tried to understand the H matrix for Dirichlet conditions. I work with a 2D mesh, and I would like my boundary to slide freely according to a vector, let's say [1 2]. I think I understood that H needs to be a Hat matrix. How to do so ? I tried to construct a hat matrix this way, but the "Angle" is not respected : x4 = 1; x1 = x4/(pi/2/Angle); k = x4+x1; x2 = sqrt(x1)*(-sqrt(k-x1)); x3 = x2; gf_model_set(md, 'add initialized data', 'H_TOP', [x1 x2;x3 x4]); gf_model_set(md, 'add initialized data', 'VECTOR_TOP', [0;0]); gf_model_set(md, 'add generalized Dirichlet condition with multipliers', mim, 'u', mfu, TOP_BOUND,'VECTOR_TOP', 'H_TOP'); Simon
