1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
% Copyright (C) 2015-2016 FABRE Mathieu, SECK Mamadou, DALLERIT Valentin,
%
% This file is a part of GetFEM++
%
% GetFEM++ is free software; you can redistribute it and/or modify it
% under the terms of the GNU Lesser General Public License as published
% by the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version along with the GCC Runtime Library
% Exception either version 3.1 or (at your option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
% License and GCC Runtime Library Exception for more details.
% You should have received a copy of the GNU Lesser General Public License
% along with this program; if not, write to the Free Software Foundation,
% Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
% Simple supported Mindlin-Reissner plate
Emodulus = 1; % Young Modulus
nu = 0.5; % Poisson Coefficient
epsilon = 0.001; % Plate thickness
kappa = 5/6; % Shear correction factor
f = -5*epsilon^3; % Prescribed force on the top of the plate
variant = 2; % 0 : not reduced, 1 : with reduced integration,
% 2 : MITC reduction
quadrangles = true; % Locking free only on quadrangle for the moment
K = 1; % Degree of the finite element method
with_Mindlin_brick = true; % Uses the Reissner-Mindlin predefined brick or not
dirichlet_version = 1; % 0 = simplification, 1 = with multipliers
% 2 = penalization
r = 1E8; % Penalization parameter.
plot_mesh = false;
draw_solution = true;
% trace on;
gf_workspace('clear all');
NX = 80;
if (quadrangles)
m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
else
m=gf_mesh('import','structured',sprintf('GT="GT_PK(2,1)";SIZES=[1,1];NOISED=0;NSUBDIV=[%d,%d];', NX, NX));
end
% Create a mesh_fem for a 2D vector field
mftheta = gf_mesh_fem(m,2);
mfu = gf_mesh_fem(m,1);
if (quadrangles)
gf_mesh_fem_set(mftheta,'fem',gf_fem(sprintf('FEM_QK(2,%d)', K)));
gf_mesh_fem_set(mfu,'fem',gf_fem(sprintf('FEM_QK(2,%d)', K)));
mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,6)'));
mim_reduced = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,1)'));
else
gf_mesh_fem_set(mftheta,'fem',gf_fem(sprintf('FEM_PK(2,%d)', K)));
gf_mesh_fem_set(mfu,'fem',gf_fem(sprintf('FEM_PK(2,%d)', K)));
mim = gf_mesh_im(m, gf_integ('IM_TRIANGLE(6)'));
mim_reduced = gf_mesh_im(m, gf_integ('IM_TRIANGLE(1)'));
end
% Detect the border of the mesh and assign it the boundary number 1
border = gf_mesh_get(m,'outer faces');
gf_mesh_set(m, 'boundary', 1, border);
if (plot_mesh)
gf_plot_mesh(m, 'regions', [1]); % the boundary edges appears in red
pause(1);
end
% Build the model
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mfu);
gf_model_set(md, 'add fem variable', 'theta', mftheta);
gf_model_set(md, 'add initialized data', 'E', Emodulus);
gf_model_set(md, 'add initialized data', 'nu', nu);
gf_model_set(md, 'add initialized data', 'epsilon', epsilon);
gf_model_set(md, 'add initialized data', 'kappa', kappa);
if (with_Mindlin_brick)
gf_model_set(md, 'add Mindlin Reissner plate brick', mim, mim_reduced, 'u', 'theta', 'E', 'nu', 'epsilon', 'kappa', variant);
else
gf_model_set(md, 'add elementary rotated RT0 projection', 'RT0_projection');
gf_model_set(md, 'add linear generic assembly brick', mim, '(E*epsilon*epsilon*epsilon*(1-nu)/(48 * (1 - nu*nu))) * ((Grad_theta+Grad_theta''):(Grad_Test_theta+Grad_Test_theta''))');
gf_model_set(md, 'add linear generic assembly brick', mim, '(E*epsilon*epsilon*epsilon*nu/(12 * (1 - nu*nu))) * (Div_theta*Div_Test_theta)');
if (variant == 0)
gf_model_set(md, 'add linear generic assembly brick', mim, '(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Grad_Test_u) + (E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Test_theta)');
elseif (variant == 1)
gf_model_set(md, 'add linear generic assembly brick', mim_reduced, '(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Grad_Test_u) + (E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Test_theta)');
else
gf_model_set(md, 'add linear generic assembly brick', mim, '(E*kappa*epsilon/(1 + nu)) * ((Grad_u + Elementary_transformation(theta,RT0_projection)).Grad_Test_u) + (E*kappa*epsilon/(1 + nu)) * ((Grad_u + Elementary_transformation(theta, RT0_projection)).(Elementary_transformation(Test_theta, RT0_projection)))');
end
end
gf_model_set(md, 'add initialized data', 'VolumicData', f);
gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');
gf_model_set(md, 'add initialized data', 'DirichletData', 0);
switch (dirichlet_version)
case 0,
gf_model_set(md, 'add Dirichlet condition with simplification', 'u', 1, 'DirichletData');
case 1,
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', mfu, 1, 'DirichletData');
case 2,
gf_model_set(md, 'add Dirichlet condition with penalization', mim, 'u', r, 1, 'DirichletData');
end
gf_model_get(md, 'solve');
U = gf_model_get(md, 'variable', 'u');
if (draw_solution)
gf_plot(mfu,U,'mesh','off', 'zplot', 'on');
colorbar; title('computed solution');
end
|
