-
Notifications
You must be signed in to change notification settings - Fork 11
/
FemCase.m
170 lines (165 loc) · 6.35 KB
/
FemCase.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
classdef FemCase < handle
% FemCase Class
% Works a the collection of data structures that represent a FEM
% problem. It includes a mesh, the physical problem, it's boundary
% conditions and the results.
% Main Methods:
% - solve
% - eigen_values
properties
mesh
physics
loads
bc
initial_dis
dis
reactions
end
methods
function obj = FemCase(mesh_in,physics_in)
require(isa(mesh_in,'Mesh'), ...
'ArgumentError: mesh_in not type Mesh');
require(isa(physics_in,'Physics'), ...
'ArgumentError: mesh_in not type Physics');
obj.mesh = mesh_in;
obj.physics = physics_in;
% Create the rest of the Properties
obj.bc = obj.compound_function(false);
obj.loads = obj.compound_function(0);
obj.initial_dis = obj.compound_function(0);
obj.dis = obj.compound_function(0);
obj.reactions = obj.compound_function(0);
end
function [Z,D] = eigen_values(fem,mode_number)
% [V,D] = eigen_values(fem,mode_number)
% Solves the eigen value problem for K and M.
% V [n_dof x n_dof][Float]: Mode Shape Matrix
% D [n_dof x n_dof][Float]: Diagonal matrix with omega^2
% mode_number [Int]: Number of modes to be calculated
require(~isempty(fem.physics.m), ...
'Error: physics needs a mass matrix');
F = ~fem.bc.all_dofs;
S = fem.S;
M = fem.M;
[V1,D] = eigs(S(F,F),M(F,F),mode_number,'SM');
Z = fem.compound_function_array(0,mode_number);
dof_list = Z(1).all_dofs; % Starts empty
for n = 1:mode_number
dof_list(F) = V1(:,n); % Is filled
Z(n).dof_list_in(dof_list);
end
end
function solve(fem)
% Side Effects
% Solves the problem set by physics.k and the mesh.
% BC and Loads to DOF form
L = fem.loads.all_dofs;
F = ~fem.bc.all_dofs;
% Create Stiffness
S = fem.mesh.assembly_matrix(fem.physics.dofs_per_node, ...
fem.physics.dofs_per_ele, ...
fem.physics.k);
D = zeros(size(S,1),1);
I = fem.initial_dis.non_zero_index;
if any(I)
ID = fem.initial_dis.all_dofs;
NI = ~I;
FF = F(NI);
DD = D(NI);
SS = S(NI,NI);
LL = L(NI) - S(NI,I)*ID(I);
else
FF = F;
DD = D;
SS = S;
LL = L;
end
% For unit reasons (i.e. piezoelectric constant, some terms in
% the stiffness matrix (S) are several orders of magnitude
% larger than others. It is desirable to improve cond(S) before
% solving the system.
if (condest(SS(FF,FF)) > 1e8)
P = diag(sqrt(diag(SS)));
inv_P = inv(P);
SS = inv_P'*SS*inv_P;
LP = inv_P*LL;
DD(FF) = inv_P(FF,FF)*(SS(FF,FF) \ LP(FF));
% Change to true if output is wanted
if (false)
condest(SS(FF,FF))
size(SS(FF,FF))
sprank(SS(FF,FF))
end
else
DD(FF) = SS(FF,FF) \ LL(FF);
end
if any(I)
D(NI) = DD;
D(I) = ID(I);
else
D = DD;
end
fem.dis.dof_list_in(D);
fem.reactions.dof_list_in(S*D);
edge = fem.mesh.find_nodes(@(x,y,z) (abs(x)<1e-5));
fem.reactions.node_vals.vals(edge,3);
sum(abs(fem.reactions.node_vals.vals));
sum(fem.loads.node_vals.vals);
end
%% Wrappers
function L = compound_function_array(fem,filler,n_row)
% L = compound_function_array(fem,filler,n_row,n_col)
% Creates an array of compound functions that match the dofs in
% the FemCase
% n_row: [Int]: Number of rows
% n_col: [Int]: Number of columns
L = CompoundFunction.empty(n_row,0);
for i = 1:n_row
L(i) = fem.compound_function(filler);
end
end
function S = S(fem)
% S = S(fem)
% Wrapper for Stiffness Matrix
S = fem.mesh.assembly_matrix(fem.physics.dofs_per_node, ...
fem.physics.dofs_per_ele, ...
fem.physics.k);
end
function M = M(fem)
% M = M(fem)
% Wrapper for Mass Matrix
M = fem.mesh.assembly_matrix(fem.physics.dofs_per_node, ...
0, ... %% HARDCODED for mechanical vibration
fem.physics.m);
%% Adds point masses if needed
if ~isempty(fem.mesh.nodes_with_mass)
for i = 1:length(fem.mesh.nodes_with_mass)
n = fem.mesh.nodes_with_mass(i);
m = fem.mesh.mass_values(i);
n_dofs = fem.mesh.node_dofs(fem.physics.dofs_per_node,n);
n_dofs = n_dofs(1:3);
M(n_dofs,n_dofs) = M(n_dofs,n_dofs) + eye(3)*m;
end
end
end
function obj = compound_function(fem,filler)
% Wrapper for class CompoundFunction
obj = CompoundFunction(filler,fem.mesh.n_nodes, ...
fem.physics.dofs_per_node,fem.mesh.n_ele,fem.physics.dofs_per_ele);
end
function plot(fem,scale)
if ~isempty(fem.dis) && ~isempty(fem.reactions)
fem.mesh.plot_displacements(scale,fem.dis.node_vals.vals(:,1:3))
hold on
fem.mesh.plot_vector_field( ...
fem.reactions.node_vals.vals(:,1:3),'k');
else
hold on
fem.mesh.plot
end
fem.mesh.plot_vector_field(fem.bc.node_vals.vals(:,1:3),'b')
fem.mesh.plot_vector_field(fem.loads.node_vals.vals(:,1:3),'r')
hold off
end
end
end