FEMBarComputer/StiffnessMatrixComputer.m at master · jose-antonio-torres/FEMBarComputer · GitHub

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
classdef StiffnessMatrixComputer < handle

    properties (Access = private)
        Kel
        R
        nElem
        x
        Tnod
        Td
        mat1
        mat2
        Kprima
    end

    properties (Access = public)
        K
    end

    methods (Access = public)

        function obj = StiffnessMatrixComputer(cParams)
            obj.init(cParams);
        end

        function compute(obj)
            obj.computeElementalStiffnessMatrices();
            obj.globalStiffnessMatrixAssembly();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.Kel   = zeros(4,4,12);
            obj.K     = zeros(16,16);
            obj.nElem = cParams.nElem;
            obj.x     = cParams.x;
            obj.Tnod  = cParams.Tnod;
            obj.Td    = cParams.Td;
            obj.mat1  = cParams.mat1;
            obj.mat2  = cParams.mat2;
        end

        function computeEuclidianMatrix(obj,s)
            Robj  = EuclidianMatrixComputer(s);
            Robj.compute();
            obj.R = Robj.R;
        end

        function computeElementalStiffnessMatrices(obj)
            Ke = zeros(4,4,obj.nElem);
            for iel=1:obj.nElem
                s.x        = obj.x;
                s.Tnod     = obj.Tnod;
                s.iel      = iel;
                s          = obj.computeBarLength(s);
                obj.computeEuclidianMatrix(s);
                obj.computeLocalElementalStiffnessMatrix(s.l);
                Ke(:,:,iel) = obj.computeGlobalElementalStiffnessMatrix();
            end
            obj.Kel = Ke;
        end

        function Kprima = computeLocalElementalStiffnessMatrix(obj,l)
            obj.Kprima = obj.mat1*obj.mat2/l*[1 -1; -1 1];
        end

        function Kelem = computeGlobalElementalStiffnessMatrix(obj)
            Kelem = obj.R'*obj.Kprima*obj.R;
        end

        function globalStiffnessMatrixAssembly(obj)
            KG = obj.K;
            for iel = 1:obj.nElem
                for i = 1:4
                    I = obj.Td(iel,i);
                    for j = 1:4
                        J = obj.Td(iel,j);
                        KG(I,J) = KG(I,J)+obj.Kel(i,j,iel);
                    end
                end
            end
            obj.K = KG;
        end

    end

    methods (Access = private, Static)

        function s = computeBarLength(s)
            BarElem = BarElemComputer(s);
            BarElem.compute();
            s.x1      = BarElem.x1;
            s.x2      = BarElem.x2;
            s.y1      = BarElem.y1;
            s.y2      = BarElem.y2;
            s.l       = BarElem.l;
        end

    end

end