1
+ function [] = FEM_Class_1D();
2
+
3
+ % uses \ to get inverse
4
+ % arranges as NQ * NQ matrix
5
+
6
+ % ONE degree of freedom per node.
7
+ close all
8
+ clear all
9
+ clc
10
+
11
+ InputData ;
12
+ Stiffness ;
13
+ ModifyForBCon ;
14
+ StressCalc ;
15
+ ReactionCal ;
16
+ Output ;
17
+ % ------------------------ function InputData ---------------------------
18
+ function [] = InputData();
19
+ global NN NE NM NDIM NEN NDN
20
+ global ND NL NCH NPR NMPC NBW
21
+ global X NOC F AREA MAT DT S
22
+ global PM NU U MPC BT STRESS REACT
23
+ global CNST
24
+ global TITLE FILE1 FILE2
25
+ global LINP LOUT
26
+ global NQ EL
27
+
28
+ FILE1 = input(' Input Data File Name ' ,' s' );
29
+ LINP = fopen(FILE1 ,' r' );
30
+ FILE2 = input(' Output Data File Name ' ,' s' );
31
+ LOUT = fopen(FILE2 ,' w' );
32
+
33
+ TITLE = fgets(LINP )
34
+ DUMMY = fgets(LINP )
35
+ TMP = str2num(fgets(LINP ))
36
+ [NN , NE , NM , NDIM , NEN , NDN ] = deal(TMP(1 ),TMP(2 ),TMP(3 ),TMP(4 ),TMP(5 ),TMP(6 ));
37
+
38
+ NQ = NDN * NN ;
39
+
40
+ DUMMY = fgets(LINP );
41
+ TMP = str2num(fgets(LINP ));
42
+ [ND , NL , NPR , EL ]= deal(TMP(1 ),TMP(2 ), TMP(3 ), TMP(4 ));
43
+
44
+ % ----- Coordinates -----
45
+ DUMMY = fgets(LINP );
46
+ for I= 1 : NN
47
+ TMP = str2num(fgets(LINP ));
48
+ [N , X(N ,: )]=deal(TMP(1 ),TMP(2 : 1 + NDIM ));
49
+ end
50
+ % ----- Connectivity -----
51
+ DUMMY = fgets(LINP );
52
+ for I= 1 : NE
53
+ TMP = str2num(fgets(LINP ));
54
+ [N ,NOC(N ,: ), MAT(N ,: ), AREA(N ,: )] = ...
55
+ deal(TMP(1 ),TMP(2 : 1 + NEN ), TMP(2 + NEN ), TMP(3 + NEN ));
56
+ end
57
+
58
+ % ----- Specified Displacements -----
59
+ DUMMY = fgets(LINP );
60
+ for I= 1 : ND
61
+ TMP = str2num(fgets(LINP ));
62
+ [NU(I ,: ),U(I ,: )] = deal(TMP(1 ), TMP(2 ));
63
+ end
64
+ % ----- Component Loads -----
65
+ DUMMY = fgets(LINP );
66
+ F = zeros(NQ ,1 );
67
+ for I= 1 : NL
68
+ TMP = str2num(fgets(LINP ));
69
+ [N ,F(N )]=deal(TMP(1 ),TMP(2 ));
70
+ end
71
+
72
+ % ----- Material Properties -----
73
+ DUMMY = fgets(LINP );
74
+ for I= 1 : NM
75
+ TMP = str2num(fgets(LINP ));
76
+ [N , PM(N ,: )] = deal(TMP(1 ), TMP(2 : NPR + 1 ));
77
+ end
78
+
79
+ % ------------------------ function Stiffness ---------------------------
80
+ function []=Stiffness();
81
+ global NN NE NM NDIM NEN NDN
82
+ global ND NL NCH NPR NMPC NBW
83
+ global X NOC F AREA MAT DT S
84
+ global PM NU U MPC BT STRESS REACT
85
+ global CNST
86
+ global TITLE FILE1 FILE2
87
+ global LINP LOUT
88
+ global NQ EL
89
+
90
+ % ----- First initialize stiffness matrix
91
+ S = zeros(NQ ,NQ );
92
+
93
+ % ----- Stiffness Matrix modified -----
94
+ for N = 1 : NE
95
+ p = NOC(N , 1 );
96
+ s = NOC(N , 2 );
97
+ I3 = MAT(N );
98
+
99
+
100
+ EAL = PM(I3 , 1 ) * AREA(N ) / EL ;
101
+
102
+ % ----------- Element Stiffness Matrix SE() -----------
103
+
104
+ SE= [1 - 1 ; - 1 1 ];
105
+
106
+
107
+ disp(sprintf(' ..... Adding %d th Element Stiffness to Global Locations' ,N ));
108
+
109
+
110
+ row= 0 ;
111
+ for I= [p ,s ]
112
+ col= 0 ;
113
+ row= row + 1 ;
114
+ for J= [p ,s ];
115
+ col= col + 1 ;
116
+ S(I ,J )=S(I ,J )+ (SE(row ,col )*EAL );
117
+ end
118
+ end
119
+ end
120
+ Stiffness_matrix= S
121
+
122
+
123
+ % ------------------------ function ModifyForBCon ---------------------------
124
+ function []=ModifyForBCon();
125
+ global NN NE NM NDIM NEN NDN
126
+ global ND NL NCH NPR NMPC NBW
127
+ global X NOC F AREA MAT DT S
128
+ global PM NU U MPC BT STRESS REACT
129
+ global CNST J
130
+ global TITLE FILE1 FILE2 FF
131
+ global LINP LOUT
132
+ global NQ
133
+
134
+
135
+
136
+ % ----- Modify for Boundary Conditions -----
137
+ % --- Displacement BC ---
138
+ BCC= [];
139
+ SS= S ;
140
+ FF= F ;
141
+ FO= F
142
+ v= [];
143
+ BCC= zeros(NQ ,1 )
144
+ for I = 1 : ND
145
+ NI = NU(I );
146
+ v= [v ,NI ];
147
+ BC= S(: ,NI );
148
+ BCC= BCC +(BC * U(I ));
149
+ end
150
+
151
+ u= [];
152
+ for I= 1 : NQ
153
+ if I ~= v
154
+ u= [u I ];
155
+ end
156
+ end
157
+
158
+ FF= F(u );
159
+ FF= FF - BCC(u );
160
+ F(v )=U ;
161
+ SS(v ,: )=[];
162
+ SS(: ,v )=[];
163
+ [FF ]=SS \ FF
164
+ F(u )=FF
165
+
166
+ % ------------------------ function StressCalc ---------------------------
167
+ function []=StressCalc();
168
+ global NN NE NM NDIM NEN NDN
169
+ global ND NL NCH NPR NMPC NBW
170
+ global X NOC F AREA MAT DT S
171
+ global PM NU U MPC BT STRESS REACT
172
+ global CNST
173
+ global TITLE FILE1 FILE2
174
+ global LINP LOUT
175
+ global NQ EL
176
+
177
+
178
+ % ----- Stress Calculation modified for one degree of freedom system -----
179
+
180
+ for I = 1 : NE
181
+ I1 = NOC(I , 1 );
182
+ I2 = NOC(I , 2 );
183
+ I3 = MAT(I );
184
+
185
+ J1 = I1 ;
186
+ K1 = I2 ;
187
+ DLT = F(K1 ) - F(J1 );
188
+ STRESS(I ) = PM(I3 , 1 ) * (DLT / EL );
189
+ end
190
+
191
+
192
+
193
+ % ------------------------ function ReactionCal ---------------------------
194
+ function []=ReactionCal();
195
+ global NN NE NM NDIM NEN NDN
196
+ global ND NL NCH NPR NMPC NBW
197
+ global X NOC F AREA MAT DT S
198
+ global PM NU U MPC BT STRESS REACT
199
+ global CNST
200
+ global TITLE FILE1 FILE2
201
+ global LINP LOUT
202
+ global NQ
203
+
204
+ % ----- Reaction Calculation -----
205
+ for I = 1 : ND
206
+ NI = NU(I );
207
+ BC= S(: ,NI )
208
+ REACT(I )=BC ' *F ;
209
+ end
210
+
211
+ % ------------------------ function Output ---------------------------
212
+ function []=Output();
213
+ global NN NE NM NDIM NEN NDN
214
+ global ND NL NMPC NBW
215
+ global X NOC F AREA MAT DT S
216
+ global PM NU U MPC BT STRESS REACT
217
+ global CNST
218
+ global TITLE FILE1 FILE2
219
+ global LINP LOUT
220
+
221
+ disp(TITLE );
222
+ fprintf(LOUT ,' %s\n ' ,TITLE );
223
+ disp(' NODE# DISPLACEMENT' );
224
+ fprintf(LOUT ,' NODE# DISPLACEMENT\n ' );
225
+ for I= 1 : NN
226
+ disp(sprintf(' %d %15.5G ' , I , F(I )));
227
+ fprintf(LOUT ,' %d %15.5G\n ' , I , F(I ));
228
+ end
229
+ % ----- Stress -----
230
+ disp(' ELEM# STRESS' );
231
+ fprintf(LOUT ,' ELEM# STRESS\n ' );
232
+ for N= 1 : NE
233
+ disp(sprintf(' %d %15.5G ' ,N , STRESS(N )));
234
+ fprintf(LOUT ,' %d %15.5G\n ' ,N , STRESS(N ));
235
+ end
236
+ % ----- Reaction -----
237
+ disp(' NODE# REACTION' );
238
+ fprintf(LOUT ,' NODE# REACTION\n ' );
239
+ for I= 1 : ND
240
+ N = NU(I );
241
+ disp(sprintf(' %d %15.5G ' ,N , REACT(I )));
242
+ fprintf(LOUT ,' %d %15.5G\n ' , N , REACT(I ));
243
+ end
244
+
245
+ fclose(LOUT );
246
+ disp(sprintf(' RESULTS ARE IN FILE : %s ' , FILE2 ));
0 commit comments