diffpy/Node3.py at master · nolanlad/diffpy · 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
104
105
106
107
108
109
110
import numpy as np
import matplotlib.pyplot as plt
from utilities import *


'''
def is_independent(A,v):
    v = np.array(v)
    A2 = np.concatenate((A,np.transpose(v.reshape(1,-1))),axis=1)
    conj = np.transpose(A2)@A2
    if np.abs(np.linalg.det(conj)) != 0:
        if A2.shape[0] == A2.shape[1]:
            if np.abs(np.linalg.det(A2)) == 0:
                return A,False
        return A2,True
    return A,False
'''


def is_independent(A,v):
    ''' Checks the linear dependence between the vectos in A and v'''
    v = np.array(v)
    svdsol = svdsolve(A,v)
    diff = np.abs(A@svdsol - v)
    if np.max(diff) < 1e-2:
        return A,False
    else:
        return np.concatenate((A,np.transpose(v.reshape(1,-1))),axis=1),True

def get_node_soln(i,nodes,D,h,dim):
    '''Chosses linearly independent nodes to calculate the discrete derivative
    it is an attempt to make a very robust node choosing scheme that does not break
    under any circumstances'''
    n = nodes.get_node(i)
    diffx = n.get_dim('x') - nodes.get_dim('x')
    diffy = n.get_dim('y') - nodes.get_dim('y')
    dist = np.sqrt(diffx**2 + diffy**2)
    arr = np.argsort(dist)
    A = np.array(make_combos([diffx[arr[0]],diffy[arr[0]] ],D+h))
    A = np.transpose(A.reshape(1,-1))
    inds = [arr[0]]
    #print(A.shape)
    for i in arr[1:]:
        v = np.array(make_combos([diffx[arr[i]],diffy[arr[i]] ],D+h))
        A,stat = is_independent(A,v)
        if stat:
            inds.append(arr[i])
        if A.shape[0] == A.shape[1]:
            break
    thelist = make_combos2('xy',D+h)
    b = np.zeros(len(thelist))
    b[thelist.index(dim*D)] = factorial(D)
    soln = svdsolve(A,b)
    #soln = np.linalg.solve(A,b)
    return inds,soln

def test(n,nodes,dim,D,h):
    diffx = n.get_dim('x') - nodes.get_dim('x')
    diffy = n.get_dim('y') - nodes.get_dim('y')
    dist = np.sqrt(diffx**2 + diffy**2)
    arr = np.argsort(dist)
    A = np.array(make_combos([diffx[arr[0]],diffy[arr[0]] ],D+h))
    A = np.transpose(A.reshape(1,-1))
    inds = [arr[0]]
    #print(A.shape)
    #print(n.get_dim('x'),n.get_dim('y'))
    for i in arr[1:]:
        v = np.array(make_combos([diffx[arr[i]],diffy[arr[i]] ],D+h))
        A,stat = is_independent(A,v)
        if stat:
            inds.append(arr[i])
            #print(v)
            #print(diffx[arr[i]],diffy[arr[i]])
            #print(nodes.get_dim('x')[arr[i]],nodes.get_dim('y')[arr[i]])
        if A.shape[0] == A.shape[1]:
            break
    thelist = make_combos2('xy',D+h)
    b = np.zeros(len(thelist))
    b[thelist.index(dim*D)] = factorial(D)
    soln = svdsolve(A,b)
    return soln,inds,A,b

def make_stiffness_general(nodes,D,h,dim,labelname = None):
    n = nodes
    M = n.lennodes
    A = np.zeros((M,M))
    if D == 0:
        if labelname:
            ind = n.labnames[labelname]
            w = np.where(n.labels == ind)[0]
            for i in w:
                A[i][i] = 1.
            return A
        else:
            return np.eye(M)
    if not labelname:
        for i in range(0,n.lennodes):
            #print(i)
            #ids,soln = get_node_soln(n.get_node(i),n,dim,D,h)
            ids,soln = get_node_soln(i,n,D,h,dim)
            A[i][ids] = soln
    else:
        ind = n.labnames[labelname]
        w = np.where(n.labels == ind)[0]
        for i in w:
            ids,soln = get_node_soln(i,n,D,h,dim)
            #ids,soln = get_node_soln(n.get_node(i),n,dim,D,h)
            A[i][ids] = soln

    return A