RheoKelvinVoigt1D.m 5.32 KB
Newer Older
Andreas Zilian's avatar
Andreas Zilian committed
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
classdef RheoKelvinVoigt1D < mafe.Element
    % Kelvin-Voigt rheological element, uni-directional (axial in 1 or 2, or rotational)

    properties
        node = mafe.Node.empty(0,2); % the two nodes of the rheo device
        sect = mafe.Section0D(); % the section descriptor
        tdof = mafe.DofType.Disp1; % the active degree of freedom (direction)
    end

    methods
        %% constructor
        function obj = RheoKelvinVoigt1D(nodes, section, tdof)
            if nargin >= 2
                obj.node = nodes;
                obj.sect = section;
            end
            if nargin >= 3
                obj.tdof = tdof;
            end
        end
        %% activation of dofs at the nodes of the element
        function activateDofTypes(self)
            ineed = [ self.tdof ];
            for i = 1:length(self.node)
                self.node(i).dof = [self.node(i).dof ineed];
            end
        end
        %% provide element index vector of dof
        function idx = getDofSysIndices(self)
Andreas Zilian's avatar
Andreas Zilian committed
30
			      dtype = find( [self.node.dof] == self.tdof );
Andreas Zilian's avatar
Andreas Zilian committed
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
            ndidx = [ self.node.idx ];
            idx = [ ndidx(dtype(1)), ndidx(dtype(2)) ];
        end
        %% calculate transformation matrix local->global
        function [L, T] = transform(self)
            L = 0.0;
            T = eye(2); % not needed in fact
        end
        %% calculation of element stiffness matrix
        function ele_mat = stiffness(self)
            % get parameters
            k = self.sect.k;
            % element matrix
            ele_mat = k*[ +1,  -1; ...
                          -1,  +1 ];
        end
        %% calculation of element mass matrix
        function ele_mat = mass(self)
            % local element matrix
            ele_mat = [ +0, +0; +0, +0 ];
        end
        %% calculation of element damping matrix
        function ele_mat = damping(self)
            % get parameters
            d = self.sect.d;
            % local element matrix
            ele_mat = d*[ +1, -1; ...
                          -1, +1 ];
Andreas Zilian's avatar
Andreas Zilian committed
59
60
61
62
63
64
            % provide damping matrix based on Rayleigh assumption
            % get parameters
            d_alpha = self.sect.d_alpha;
            d_beta  = self.sect.d_beta;
            % D = alpha * M + beta * K
            ele_mat = ele_mat + d_alpha * self.mass() + d_beta * self.stiffness();
Andreas Zilian's avatar
Andreas Zilian committed
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
        end
        %% calculation of element force vector
        function ele_vec = force(self)
            % local element vector
            ele_vec = [ +0; +0 ];
        end
        %------------------------------------------------------------------
        %% graphical output: plot element
        function plot(self, config, scale)
            % my nodes
            n1 = self.node(1);
            n2 = self.node(2);
            % coordinates and displacements
            x1 = [ n1.ref(1), n2.ref(1) ];
            if length(n1.ref) > 1 & length(n2.ref) > 1
                x2 = [ n1.ref(2), n2.ref(2) ];
            else
                x2 = [ 0, 0 ];
            end
            if self.tdof == mafe.DofType.Disp1
                u1 = [ n1.lhs( find(n1.dof == mafe.DofType.Disp1) ), n2.lhs( find(n2.dof == mafe.DofType.Disp1) ) ];
            else
                u1 = [ 0, 0 ];
            end
            if self.tdof == mafe.DofType.Disp2
                u2 = [ n1.lhs( find(n1.dof == mafe.DofType.Disp2) ), n2.lhs( find(n2.dof == mafe.DofType.Disp2) ) ];
            else
                u2 = [ 0, 0 ];
            end
            %
			% defaults
			lcolor = [0.6 0.6 0.6];
			mcolor = [0.2 0.2 0.2];
			lwidth = 0.5;
			mwidth = 10.;
			ecolor = [0.4 0.4 0.4];
			marker = '.';
			%
            switch config
                case 'deformed'
                    XX1 = x1+u1*scale;
                    XX2 = x2+u2*scale;
                    lcolor = [43.9, 50.6,  0.0]/100;
                    mcolor = [43.9, 50.6,  0.0]/100;
					falpha = 0.6;
                    lwidth = 4.0;
                    mwidth = 20.;
                case 'intensity'
                    XX1 = x1+u1*scale;
                    XX2 = x2+u2*scale;
                    normalforce = 0.5 * ( -self.internal(1) + self.internal(2) );
                    if normalforce < 0.0
                        lcolor = [15.3, 31.0, 54.5]/100;
                    elseif normalforce > 0.0
                        lcolor = [82.7, 58.4, 16.9]/100;
                    else
                        lcolor = [58.8, 65.5,  7.8]/100;
                    end
                    mcolor = [0.0 0.0 0.0 0.0];
					falpha = 0.6;
                    lwidth = 2.0;
                    mwidth = 0.1;
                otherwise
                    XX1 = x1;
                    XX2 = x2;
                    lcolor = [0.6 0.6 0.6];
                    mcolor = [0.6 0.6 0.6];
					falpha = 0.6;
                    lwidth = 2.0;
                    mwidth = 12.;
            end

			mafe.patchline(XX1, XX2, ...
                     'EdgeColor', lcolor, ...
                     'LineWidth', lwidth, ...
                     'EdgeAlpha', falpha);

            plot(XX1([1 end]), XX2([1 end]), ...
                     marker, ...
					 'color', mcolor, ...
					 'markersize', mwidth);
        end
        %% print element data
        function str = print(self)
            str = sprintf('%+4.3e ', [ -self.internal(1), ... % left F
                                        self.internal(2)  ]); % right F
        end
    end

end