Skip to content

diffuse_weights

diffuse_weights(Vv, Tv, phi, bI, dt=None, normalize=True)

Performs a diffusion on the tet mesh Vv, Tv at nodes bI for time dt.

Parameters:

Name Type Description Default
Vv (n, 3) float numpy array

Mesh vertices

 required
Tv (t, 4) int numpy array

Mesh tets

 required
phi (c, b) float numpy array

Quantity to diffuse

 required
bI (c, b) int numpy array

Indices at diffusion points

 required
dt float

Time to diffuse for

 None
normalize bool

Whether to normalize the weights

 True

Returns:

Name Type Description
W (n, b) float numpy array

Diffused quantities over entire mesh
Source code in src\fast_cody\diffuse_weights.py 
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
def diffuse_weights(Vv, Tv, phi, bI,  dt=None, normalize=True):
    """ Performs a diffusion on the tet mesh Vv, Tv at nodes bI for time dt.

    Parameters
    ----------
    Vv : (n, 3) float numpy array
        Mesh vertices
    Tv : (t, 4) int numpy array
        Mesh tets
    phi : (c, b) float numpy array
        Quantity to diffuse
    bI : (c, b) int numpy array
        Indices at diffusion points
    dt : float
        Time to diffuse for
    normalize : bool
        Whether to normalize the weights

    Returns
    -------
    W : (n, b) float numpy array
        Diffused quantities over entire mesh

    """

    if (dt is None):
        dt = np.mean(igl.edge_lengths(Vv, Tv)) ** 2

    L = laplacian(Vv, Tv)
    M = igl.massmatrix(Vv, Tv)

    Q = L * dt + M

    ii = np.setdiff1d(  np.arange(Q.shape[0]), bI)
    # selection matrix for indices bI
    Qii = Q[ii, :][:, ii]
    Qib = Q[ii, :][:, bI]

    Wii = fast_cody.umfpack_lu_solve(Qii, -Qib @ phi)
    W = np.zeros((L.shape[0], Wii.shape[1]))
    W[ii, :] = Wii
    W[bI, :] = phi
    # W = gpt.min_quad_with_fixed(L*dt + M, k=bI, y=phi)

    # normalize weights so that max is 1 and min is 0
    # W = (W - np.min(W, axis=0)[:, None]) / (np.max(W, axis=0)[:, None] - np.min(W, axis=0)[:, None])

    # ps.init()
    # # pc = ps.register_point_cloud("pc", Vs)
    # # pc_CI = ps.register_point_cloud("pc_cI", Vs[CI])
    # mesh = ps.register_volume_mesh("mesh", Vv, Tv)
    # mesh.add_scalar_quantity("w", W.flatten())
    # ps.show()


    # normalize between 0 and 1
    if W.ndim == 1:
        W = W[:, None]
    if normalize:
        W = (W - np.min(W, axis=0)) / (np.max(W, axis=0) - np.min(W, axis=0))
    # WeightsViewer(Vs, Ts, Ws, period=1)
    # WeightsViewer(Vv, Fv, W)

    # ps.init()
    # volms = ps.register_volume_mesh("vol", Vv, Tv)
    # volms.add_scalar_quantity("W", W[:, 0])
    #
    # ps.show()

    # ps.init()
    # mesh = ps.register_volume_mesh("mesh", Vv, Tv)
    # mesh.add_scalar_quantity("weights", W[:, 0], enabled=True, cmap='coolwarm')
    # ps.show()
    # WeightsViewer(Vv, Tv, W)
    return W