Galvanotaxis/elasticForces_modified.m at main · CoposGroup/Galvanotaxis · 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
% Calina Copos
% Elastic force computation on triangulated mesh

function F = elasticForces_modified(T,X0,X,kk,ell0,cond_elastic)

    N = length(X);
    Ntris = length(T);

    % zero vectors
    Fc = zeros(N,2);
    dA_T = zeros(Ntris,1);
    dA = zeros(N,1);

    % fetch edges
    TR = TriRep(T,X0(:,1),X0(:,2));
    E = edges(TR);
    dL = zeros(length(E),1);

    midpt = mean(X0(:,1));

    % area elements
    dA_T = 0.5*( (X(T(:,2),1)-X(T(:,1),1)).*(X(T(:,3),2)-X(T(:,1),2)) - ...
        (X(T(:,2),2)-X(T(:,1),2)).*(X(T(:,3),1)-X(T(:,1),1)) );
    for k=1:Ntris
        dA(T(k,1)) = dA(T(k,1)) + dA_T(k)/3;
        dA(T(k,2)) = dA(T(k,2)) + dA_T(k)/3;
        dA(T(k,3)) = dA(T(k,3)) + dA_T(k)/3;
    end

    % for each edge
    for k=1:length(E)
        id1 = E(k,1);
        id2 = E(k,2);
        dl = sqrt((X(id2,1)-X(id1,1))^2 + (X(id2,2)-X(id1,2))^2);
        dl0 = ell0(k);

        dk = (8*kk/(3*dl0))*(dA(id1)+dA(id2))/2;
        dL(k) = dl;

        ddA = (dA(id1)+dA(id2))/2;

        % nonzero resting length
        if cond_elastic == 0
            if dl<1.21*dl0
                Fc(id1,1) = Fc(id1,1) + (dk/ddA)*(dl/dl0 - 1.0)*(X(id2,1)-X(id1,1))/dl;
                Fc(id1,2) = Fc(id1,2) + (dk/ddA)*(dl/dl0 - 1.0)*(X(id2,2)-X(id1,2))/dl;
                Fc(id2,1) = Fc(id2,1) + (dk/ddA)*(dl/dl0 - 1.0)*(X(id1,1)-X(id2,1))/dl;
                Fc(id2,2) = Fc(id2,2) + (dk/ddA)*(dl/dl0 - 1.0)*(X(id1,2)-X(id2,2))/dl;
            end
        elseif cond_elastic ~= 0
            if (cond_elastic == 1 && (X0(id1,1)<midpt && X0(id2,1)<midpt) )
                Fc(id1,1) = Fc(id1,1) + (dk/ddA)*dl*(X(id2,1)-X(id1,1))/dl;
                Fc(id1,2) = Fc(id1,2) + (dk/ddA)*dl*(X(id2,2)-X(id1,2))/dl;
                Fc(id2,1) = Fc(id2,1) + (dk/ddA)*dl*(X(id1,1)-X(id2,1))/dl;
                Fc(id2,2) = Fc(id2,2) + (dk/ddA)*dl*(X(id1,2)-X(id2,2))/dl;
            elseif (cond_elastic == 2 && (X0(id1,1)>=midpt && X0(id2,1)>=midpt) )
                Fc(id1,1) = Fc(id1,1) + (dk/ddA)*dl*(X(id2,1)-X(id1,1))/dl;
                Fc(id1,2) = Fc(id1,2) + (dk/ddA)*dl*(X(id2,2)-X(id1,2))/dl;
                Fc(id2,1) = Fc(id2,1) + (dk/ddA)*dl*(X(id1,1)-X(id2,1))/dl;
                Fc(id2,2) = Fc(id2,2) + (dk/ddA)*dl*(X(id1,2)-X(id2,2))/dl;
            elseif (cond_elastic == 3)
                Fc(id1,1) = Fc(id1,1) + (dk/ddA)*dl*(X(id2,1)-X(id1,1))/dl;
                Fc(id1,2) = Fc(id1,2) + (dk/ddA)*dl*(X(id2,2)-X(id1,2))/dl;
                Fc(id2,1) = Fc(id2,1) + (dk/ddA)*dl*(X(id1,1)-X(id2,1))/dl;
                Fc(id2,2) = Fc(id2,2) + (dk/ddA)*dl*(X(id1,2)-X(id2,2))/dl;
            else
                Fc(id1,1) = Fc(id1,1) + (dk/ddA)*(dl/dl0 - 1.0)*(X(id2,1)-X(id1,1))/dl;
                Fc(id1,2) = Fc(id1,2) + (dk/ddA)*(dl/dl0 - 1.0)*(X(id2,2)-X(id1,2))/dl;
                Fc(id2,1) = Fc(id2,1) + (dk/ddA)*(dl/dl0 - 1.0)*(X(id1,1)-X(id2,1))/dl;
                Fc(id2,2) = Fc(id2,2) + (dk/ddA)*(dl/dl0 - 1.0)*(X(id1,2)-X(id2,2))/dl;
            end

        end

    F(:,1) = Fc(:,1);
    F(:,2) = Fc(:,2);
    end
end