-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathEX662.m
218 lines (180 loc) · 9.25 KB
/
EX662.m
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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
%-----------------------------------------------------------------------------
% Example 6.6.2
% to solve the two-dimensional Laplace's equation given as
% u,xx + u,yy =0, 0 < x < 5, 0 < y < 10
% u(x,0) = 0, u(x,10) = 100sin(pi*x/10),
% u(0,y) = 0, u,x(5,y) = 0
% using isoparametric four-node quadrilateral elements
% (see Fig. 5.9.2 for the finite element mesh)
%
% Variable descriptions
% k = element matrix
% f = element vector
% kk = system matrix
% ff = system vector
% gcoord = coordinate values of each node
% nodes = nodal connectivity of each element
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofs in 'bcdof'
% point2 - integration (or sampling) points
% weight2 - weighting coefficients
% nglx - number of integration points along x-axis
% ngly - number of integration points along y-axis
% xcoord - x coordinate values of nodes
% ycoord - y coordinate values of nodes
% jacob2 - jacobian matrix
% shape - four-node quadrilateral shape functions
% dhdr - derivatives of shape functions w.r.t. natural coord. r
% dhds - derivatives of shape functions w.r.t. natural coord. s
% dhdx - derivatives of shape functions w.r.t. physical coord. x
% dhdy - derivatives of shape functions w.r.t. physical coord. y
%-----------------------------------------------------------------------------
clear
%------------------------------------
% input data for control parameters
%------------------------------------
nel=16; % number of elements
nnel=4; % number of nodes per element
ndof=1; % number of dofs per node
nnode=25; % total number of nodes in system
nglx=2; ngly=2; % use 2x2 integration rule
sdof=nnode*ndof; % total system dofs
edof=nnel*ndof; % dofs per element
%---------------------------------------------
% input data for nodal coordinate values
% gcoord(i,j) where i->node no. and j->x or y
%---------------------------------------------
gcoord(1,1)=0.0; gcoord(1,2)=0.0; gcoord(2,1)=1.25; gcoord(2,2)=0.0;
gcoord(3,1)=2.5; gcoord(3,2)=0.0; gcoord(4,1)=3.75; gcoord(4,2)=0.0;
gcoord(5,1)=5.0; gcoord(5,2)=0.0; gcoord(6,1)=0.0; gcoord(6,2)=2.5;
gcoord(7,1)=1.25; gcoord(7,2)=2.5; gcoord(8,1)=2.5; gcoord(8,2)=2.5;
gcoord(9,1)=3.75; gcoord(9,2)=2.5; gcoord(10,1)=5.0; gcoord(10,2)=2.5;
gcoord(11,1)=0.0; gcoord(11,2)=5.0; gcoord(12,1)=1.25; gcoord(12,2)=5.0;
gcoord(13,1)=2.5; gcoord(13,2)=5.0; gcoord(14,1)=3.75; gcoord(14,2)=5.0;
gcoord(15,1)=5.0; gcoord(15,2)=5.0; gcoord(16,1)=0.0; gcoord(16,2)=7.5;
gcoord(17,1)=1.25; gcoord(17,2)=7.5; gcoord(18,1)=2.5; gcoord(18,2)=7.5;
gcoord(19,1)=3.75; gcoord(19,2)=7.5; gcoord(20,1)=5.0; gcoord(20,2)=7.5;
gcoord(21,1)=0.0; gcoord(21,2)=10.; gcoord(22,1)=1.25; gcoord(22,2)=10.;
gcoord(23,1)=2.5; gcoord(23,2)=10.; gcoord(24,1)=3.75; gcoord(24,2)=10.;
gcoord(25,1)=5.0; gcoord(25,2)=10.;
%---------------------------------------------------------
% input data for nodal connectivity for each element
% nodes(i,j) where i-> element no. and j-> connected nodes
%---------------------------------------------------------
nodes(1,1)=1; nodes(1,2)=2; nodes(1,3)=7; nodes(1,4)=6;
nodes(2,1)=2; nodes(2,2)=3; nodes(2,3)=8; nodes(2,4)=7;
nodes(3,1)=3; nodes(3,2)=4; nodes(3,3)=9; nodes(3,4)=8;
nodes(4,1)=4; nodes(4,2)=5; nodes(4,3)=10; nodes(4,4)=9;
nodes(5,1)=6; nodes(5,2)=7; nodes(5,3)=12; nodes(5,4)=11;
nodes(6,1)=7; nodes(6,2)=8; nodes(6,3)=13; nodes(6,4)=12;
nodes(7,1)=8; nodes(7,2)=9; nodes(7,3)=14; nodes(7,4)=13;
nodes(8,1)=9; nodes(8,2)=10; nodes(8,3)=15; nodes(8,4)=14;
nodes(9,1)=11; nodes(9,2)=12; nodes(9,3)=17; nodes(9,4)=16;
nodes(10,1)=12; nodes(10,2)=13; nodes(10,3)=18; nodes(10,4)=17;
nodes(11,1)=13; nodes(11,2)=14; nodes(11,3)=19; nodes(11,4)=18;
nodes(12,1)=14; nodes(12,2)=15; nodes(12,3)=20; nodes(12,4)=19;
nodes(13,1)=16; nodes(13,2)=17; nodes(13,3)=22; nodes(13,4)=21;
nodes(14,1)=17; nodes(14,2)=18; nodes(14,3)=23; nodes(14,4)=22;
nodes(15,1)=18; nodes(15,2)=19; nodes(15,3)=24; nodes(15,4)=23;
nodes(16,1)=19; nodes(16,2)=20; nodes(16,3)=25; nodes(16,4)=24;
%-------------------------------------
% input data for boundary conditions
%-------------------------------------
bcdof(1)=1; % first node is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=2; % second node is constrained
bcval(2)=0; % whose described value is 0
bcdof(3)=3; % third node is constrained
bcval(3)=0; % whose described value is 0
bcdof(4)=4; % 4th node is constrained
bcval(4)=0; % whose described value is 0
bcdof(5)=5; % 5th node is constrained
bcval(5)=0; % whose described value is 0
bcdof(6)=6; % 6th node is constrained
bcval(6)=0; % whose described value is 0
bcdof(7)=11; % 11th node is constrained
bcval(7)=0; % whose described value is 0
bcdof(8)=16; % 16th node is constrained
bcval(8)=0; % whose described value is 0
bcdof(9)=21; % 21st node is constrained
bcval(9)=0; % whose described value is 0
bcdof(10)=22; % 22nd node is constrained
bcval(10)=38.2683; % whose described value is 38.2683
bcdof(11)=23; % 23rd node is constrained
bcval(11)=70.7107; % whose described value is 70.7107
bcdof(12)=24; % 24th node is constrained
bcval(12)=92.3880; % whose described value is 92.3880
bcdof(13)=25; % 25th node is constrained
bcval(13)=100; % whose described value is 100
%-----------------------------------------
% initialization of matrices and vectors
%-----------------------------------------
ff=zeros(sdof,1); % initialization of system force vector
kk=zeros(sdof,sdof); % initialization of system matrix
index=zeros(nnel*ndof,1); % initialization of index vector
%-----------------------------------------------------------
% loop for computation and assembly of element matrices
%-----------------------------------------------------------
[point2,weight2]=feglqd2(nglx,ngly); % sampling points & weights
for iel=1:nel % loop for the total number of elements
for i=1:nnel
nd(i)=nodes(iel,i); % extract connected node for (iel)-th element
xcoord(i)=gcoord(nd(i),1); % extract x value of the node
ycoord(i)=gcoord(nd(i),2); % extract y value of the node
end
k=zeros(edof,edof); % initialization of element matrix to zero
%--------------------------------
% numerical integration
%--------------------------------
for intx=1:nglx
x=point2(intx,1); % sampling point in x-axis
wtx=weight2(intx,1); % weight in x-axis
for inty=1:ngly
y=point2(inty,2); % sampling point in y-axis
wty=weight2(inty,2) ; % weight in y-axis
[shape,dhdr,dhds]=feisoq4(x,y); % compute shape functions and
% derivatives at sampling point
jacob2=fejacob2(nnel,dhdr,dhds,xcoord,ycoord); % compute Jacobian
detjacob=det(jacob2); % determinant of Jacobian
invjacob=inv(jacob2); % inverse of Jacobian matrix
[dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob); % derivatives w.r.t.
% physical coordinate
%------------------------------
% compute element matrix
%------------------------------
for i=1:edof
for j=1:edof
k(i,j)=k(i,j)+(dhdx(i)*dhdx(j)+dhdy(i)*dhdy(j))*wtx*wty*detjacob;
end
end
end
end % end of numerical integration loop
index=feeldof(nd,nnel,ndof);% extract system dofs associated with element
%----------------------------------
% assemble element matrices
%----------------------------------
kk=feasmbl1(kk,k,index);
end % end of element loops
%-----------------------------
% apply boundary conditions
%-----------------------------
[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);
%----------------------------
% solve the matrix equation
%----------------------------
fsol=kk\ff;
%---------------------
% analytical solution
%---------------------
for i=1:nnode
x=gcoord(i,1); y=gcoord(i,2);
esol(i)=100*sinh(0.31415927*y)*sin(0.31415927*x)/sinh(3.1415927);
end
%------------------------------------
% print both exact and fem solutions
%------------------------------------
num=1:1:sdof;
store=[num' fsol esol']
%---------------------------------------------------------------