-
Notifications
You must be signed in to change notification settings - Fork 13
/
Problem.m
128 lines (122 loc) · 7.21 KB
/
Problem.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
classdef Problem
% Class representing a FEM problem.
% Given the plate's dimensions, a finite element and a 'problemId', it
% defines the boundary conditions in its properties: 'fixeddof',
% 'freedof' and 'F'.
properties
nelx % number of elements along the x axis
nely % number of elements along the y axis
problemId % string identifing certain boundary conditions
DOFindex % global dofs' index
fixeddof % idexes of the constrained dofs
freedof % idexes of the free dofs
F % global vector of loads
end
methods
function obj = Problem(nelx, nely, element, problemId)
obj.nelx = nelx;
obj.nely = nely;
obj.problemId = problemId;
n = element.nodes;
ndof = element.ndof;
% loads and constraints
obj.F = sparse(ndof*(nely+1)*(nelx+1), 1);
alldof = 1:ndof*(nelx+1)*(nely+1);
switch obj.problemId
% static cases
case 'a'
obj.F(ndof*(nely+1)*nelx/2 + 1) = 1; % concentrated load
left = 1:ndof:ndof*nely+1; % supported left edge
bottom = ndof*nely+1:ndof*(nely+1):ndof*((nely+1)*(nelx+1)-1)+1; % supported bottom edge
right = ndof*(nely+1)*nelx+1:ndof:ndof*((nely+1)*(nelx+1)-1)+1; % supported right edge
obj.fixeddof = union(union(left, bottom), right);
case 'b'
% same as case 'a' with clamped edge
obj.F(ndof*(nely+1)*nelx/2 + 1) = 1; % concentrated load
left = 1:ndof*(nely+1); % clamped left edge
bottom = [];
for i = 1:ndof
bottom_new = ndof*nely+i:ndof*(nely+1):ndof*((nely+1)*(nelx+1)-1)+i; % clamped bottom edge
bottom = union(bottom_new, bottom);
end
right = ndof*(nely+1)*nelx+1:ndof*(nely+1)*(nelx+1); % clamped right edge
obj.fixeddof = union(union(left, bottom), right);
case 'c'
obj.F(1) = 1; % concentrated load on the left corner
obj.F((ndof*(nely+1)*nelx/2) + 1) = 1; % concentrated load in the middle
obj.F((ndof*(nely+1)*nelx) + 1) = 1; % concentrated load on right corner
bottom = [];
for i = 1:ndof
bottom_new = ndof*nely+i:ndof*(nely+1):ndof*((nely+1)*(nelx+1)-1)+i; % clamped bottom edge
bottom = union(bottom_new, bottom);
end
obj.fixeddof = bottom;
case 'd'
obj.F(ndof*(nely+1)*nelx + 1) = 1; % concentrated load right corner
obj.fixeddof = 1:ndof*(nely+1); % left edge clamped
case 'e'
obj.F(ndof*((nely+1)*nelx/2 + nely/2) + 1) = 1; % centered concentrated load
left = 1:ndof*(nely+1); % clamped left edge
bottom = [];
top = [];
for i = 1:ndof
bottom_new = ndof*nely+i:ndof*(nely+1):ndof*((nely+1)*(nelx+1)-1)+i; % clamped bottom edge
top_new = ndof*(nely+1)+i:ndof*(nely+1):ndof*(nely+1)*nelx+i;
bottom = union(bottom_new, bottom);
top = union(top_new, top);
end
right = ndof*(nely+1)*nelx+1:ndof*(nely+1)*(nelx+1); % clamped right edge
obj.fixeddof = union(union(union(left, right), bottom), top);
% dinamic cases
case 'd1'
left = 1:ndof:ndof*nely+1; % supported left edge
bottom = ndof*nely+1:ndof*(nely+1):ndof*((nely+1)*(nelx+1)-1)+1; % supported bottom edge
right = ndof*(nely+1)*nelx+1:ndof:ndof*((nely+1)*(nelx+1)-1)+1; % supported right edge
obj.fixeddof = union(union(left, bottom), right);
case 'd2'
left = 1:ndof:ndof*nely+1; % supported left edge
bottom = ndof*nely+1:ndof*(nely+1):ndof*((nely+1)*(nelx+1)-1)+1; % supported bottom edge
right = ndof*(nely+1)*nelx+1:ndof:ndof*((nely+1)*(nelx+1)-1)+1; % supported right edge
top = 1:ndof*(nely+1):ndof*(nely+1)*nelx+1; % supported top edge
obj.fixeddof = union(union(union(left, right), bottom), top);
case 'd3'
obj.fixeddof = 1:ndof*(nely+1); % left edge clamped
% test cases
case 'test1'
obj.F(ndof*(nely+1)*nelx/2 + 1) = 1; % concentrated load in the middle
left = 1:ndof:ndof*nely+1; % supported left edge
right = ndof*(nely+1)*nelx+1:ndof:ndof*((nely+1)*(nelx+1)-1)+1; % supported right edge
obj.fixeddof = union(left, right);
case 'test2'
obj.F(1:ndof:ndof*(nely+1)*(nelx+1)) = 0.001; % ~distributed load
obj.fixeddof = 1:ndof*(nely+1); % left edge clamped
case 'test3'
obj.F(ndof*((nely+1)*nelx/2 + nely/2) + 1) = 1; % centered concentrated load
left = 1:ndof:ndof*nely+1; % supported left edge
right = ndof*(nely+1)*nelx+1:ndof:ndof*((nely+1)*(nelx+1)-1)+1; % supported right edge
obj.fixeddof = union(left, right);
case 'test4'
obj.F(ndof*((nely+1)*nelx/2 + nely/2) + 1) = 1; % centered concentrated load
left = 1:ndof:ndof*nely+1; % supported left edge
bottom = ndof*nely+1:ndof*(nely+1):ndof*((nely+1)*(nelx+1)-1)+1; % supported bottom edge
right = ndof*(nely+1)*nelx+1:ndof:ndof*((nely+1)*(nelx+1)-1)+1; % supported right edge
top = 1:ndof*(nely+1):ndof*(nely+1)*nelx+1; % supported top edge
obj.fixeddof = union(union(union(left, right), bottom), top);
end
obj.freedof = setdiff(alldof, obj.fixeddof);
% create global dofs' index
obj.DOFindex = [];
for elx = 1:nelx
for ely = 1:nely
nodesnum = [(elx-1)*(nely+1) + ely + 1
(elx)*(nely+1) + ely + 1
(elx)*(nely+1) + ely
(elx-1)*(nely+1) + ely]; % global nodes numbers of the current element
for i = 1:n
obj.DOFindex = cat(2, obj.DOFindex, (nodesnum(i)-1)*ndof+1:nodesnum(i)*ndof);
end
end
end
end
end
end