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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
|
% Copyright (C) 2017-2020 Yves Renard.
%
% 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.
% -*- matlab -*- (enables emacs matlab mode)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters for nonlinear membrane
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% pde parameters : %%%%%
% geometry
LX =50;%cm % size in X.
LY =25;% % size in Y.
NX =10;% % space step.
NY =6;% % space step.
% material properties
P0 =1500; % young modulus, X' axis (N/cm)
P1 =0.01; % poison XY coeff ( v YX will be computed as vYX=Ey*vxy/Ex)
P2 =1500; % young modulus, Y' axis (N/cm)
P3 =0; % shear modulus (N/cm), if 0, computed as G=Ex/(2* (1 + v yx))
% Bdy Limit (0 to set dirichlet on the face, 1 to set neumann BL on the face, -1 not to set any BL on the face )
bdyUp=0;
bdyLow=0;
bdyLeft=0;
bdyRight=0;
%specify type of boundary condition to impose on specific elements (set 0 for dirichlet, 1 for neumann, -1 not to set any BL on elements)
bdy_type=-1;
% specify specific elements and faces on which to impose BL ( set face nbr to -1 if not used)
bdy_element1= 9;%431 ;
bdy_face1= 1;%0;
bdy_element2= 29;%432 ;
bdy_face2= 1;%1;
bdy_element3=235;% 435 ;
bdy_face3=-1;%0 ;
bdy_element4=218;% 436 ;
bdy_face4= -1;%1 ;
bdy_element5=219;%439 ;
bdy_face5= -1;%0;
bdy_element6= -1;
bdy_face6= -1;
bdy_element7= -1 ;
bdy_face7= -1 ;
bdy_element8= -1 ;
bdy_face8= -1 ;
%imposed displacement on Dirichlet BL
dx=0;
dy=0;
dz=0;
%set to 1 to reverse imposed displacements on opposite boundaries
%(This will inverse the sign of imposed disp, dx => -dx if x < half max x, dy => -dy if y < half max y,
%to allow to impose an expansion of the domain (one border is moved towards positive values of the corresponding axis, the opposite border if moved
%the same distance towards negative values), set this param to 1
opposite_bdy_reversed=0;
%set initial random vertical displ on free dofs
%set to 1 to impose random initial displacement to avoid tg matrix singularities along Z' axis
%(needed if no pretension)
INITIAL_DISP=1;
initialDispAmplitude=1e-1; %amplitude of initial displacement
% print convex description
PRINT_CONVEXES=0; %set to 1 to print convex description
PRINT_STRESSES=0; %set to 1 to print stresses
%
% applied forces
%
%pretension (could also help avoid convergence problems due to wrinkling)
PRETENSION_X=0;% N/cm pretension on X' axis
PRETENSION_Y=0;%1e-6;% N/cm pretension on Y' axis
%source forces is applied to neumann limit region if src_type=1,otherwize volumic src term on all convexes
src_type=0;
% source forces in N/cm2 if src_type=0 (volumic forces) or in N/cm if src_type=1 (bdy forces)
FORCEX = 0; % Amplitude of the external force
FORCEY = 0;
FORCEZ = 10;
% punctual forces
PUNCTUAL_DOF1=0;% dof nbr to set punctual force
PUNCTUAL_FORCE1=0; % N
PUNCTUAL_DOF2=0;% %dof nbr to set punctual force
PUNCTUAL_FORCE2=0;
%%%%% discretisation parameters : %%%%%
MESH_TYPE = 'GT_PK(2,1)';
%%MESH_TYPE = 'GT_QK(2,1)'; % linear rectangles
%MESH_TYPE = 'GT_PRISM(3,1)'; % 3D prisms
MESH_NOISED = 0; % Set to one if you want to "shake" the mesh
FEM_TYPE = 'FEM_PK(2,1)';
%FEM_TYPE = 'FEM_PK(3,2)'; % P1 for triangles
%FEM_TYPE = 'FEM_PK_WITH_CUBIC_BUBBLE(3, 1)';
%FEM_TYPE = 'FEM_PK_WITH_CUBIC_BUBBLE(3, 2)';
%%FEM_TYPE = 'FEM_QK(2,1)'; % Q1 fem for quadrangles
%FEM_TYPE = 'FEM_PRODUCT(FEM_PK(2,1),FEM_PK(1,1))'; % tensorial product of FEM for prisms
%FEM_TYPE = 'FEM_PK_HIERARCHICAL(2,2)'; % Hierarchical PK on simplexes
%FEM_TYPE = 'FEM_PK_HIERARCHICAL_COMPOSITE(2,1,2)'; % Hierarchical PK with s divisions
FEM_TYPE_P = 'FEM_PK(2,1)'; % P1 for triangles
%FEM_TYPE_P = 'FEM_PK_DISCONTINUOUS(3,1)'; % Discontinuous P1 for triangles
% DATA_FEM_TYPE must be defined if your main FEM is not Lagrangian
%DATA_FEM_TYPE = 'FEM_PK(2,1)';
VONMISES_FEM_TYPE = 'FEM_PK_DISCONTINUOUS(2,1)';%must be discontinuous!, ! degree 1 gives vtk reading error!!
%INTEGRATION_CT='IM_TRIANGLE(6)';
%INTEGRATION = 'IM_TETRAHEDRON(6)'
INTEGRATION = 'IM_TRIANGLE(3)'; % quadrature rule for polynomials up PK fem
% to degree 6 on triangles
%INTEGRATION = 'IM_EXACT_SIMPLEX(2)'; % exact integration on triangles
%INTEGRATION = 'IM_NC(2,6)'; % newton-cotes of degree 6 on triangles
%INTEGRATION = 'IM_NC_PARALLELEPIPED(2,6)'; % newton-cotes, degree 6, QK fem
% quadrangles
%INTEGRATION ='IM_QUAD(7)';
%INTEGRATION = 'IM_NC_PRISM(3,12)'; % newton-cotes, degree 12, prims
%INTEGRATION = 'IM_NC_PRISM(2,8)'; % test jyh
%INTEGRATION = 'IM_GAUSS1D(10)'; % Gauss-Legendre integration on the
% segment of order 10
%INTEGRATION = 'IM_GAUSSLOBATTO1D(10)'; % Gauss-Lobatto-Legendre
% integration on the segment
% of order 10
%INTEGRATION = 'IM_GAUSS_PARALLELEPIPED(3,5)'; % Product of two
% IM_GAUSS1D(10) (for
% quadrangles)
%INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_GAUSS1D(5), 3)';
%INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(7), 3)';
RESIDUAL = 1E-8; % residu for iterative solvers.
MAXITER = 100;
DIRICHLET_VERSION=0;
NBSTEP =10;
%%%%% saving parameters %%%%%
ROOTFILENAME = 'nonlinear_membrane'; % Root of data files.
VTK_EXPORT = 1 % export solution to a ascii .vtk file ?
NOISY=1;
|
