Updated Gradient

Homework 1/Problem1.ipynb

@@ -1 +1,80 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from IPython.display import clear_output # Only for Python
 
+%Gradient of Elastic Energies
+
+def gradEb(xkm1, ykm1, xk, yk, xkp1, ykp1, curvature0, l_k, EI):
+    """
+    Returns the derivative of bending energy E_k^b with respect to
+    x_{k-1}, y_{k-1}, x_k, y_k, x_{k+1}, and y_{k+1}.
+
+    Parameters:
+    xkm1, ykm1 : float
+        Coordinates of the previous node (x_{k-1}, y_{k-1}).
+    xk, yk : float
+        Coordinates of the current node (x_k, y_k).
+    xkp1, ykp1 : float
+        Coordinates of the next node (x_{k+1}, y_{k+1}).
+    curvature0 : float
+        Discrete natural curvature at node (xk, yk).
+    l_k : float
+        Voronoi length of node (xk, yk).
+    EI : float
+        Bending stiffness.
+
+    Returns:
+    dF : np.ndarray
+        Derivative of bending energy.
+    """
+
+    # Nodes in 3D
+    node0 = np.array([xkm1, ykm1, 0.0])
+    node1 = np.array([xk, yk, 0])
+    node2 = np.array([xkp1, ykp1, 0])
+
+    # Unit vectors along z-axis
+    m2e = np.array([0, 0, 1])
+    m2f = np.array([0, 0, 1])
+
+    kappaBar = curvature0
+
+    # Initialize gradient of curvature
+    gradKappa = np.zeros(6)
+
+    # Edge vectors
+    ee = node1 - node0
+    ef = node2 - node1
+
+    # Norms of edge vectors
+    norm_e = np.linalg.norm(ee)
+    norm_f = np.linalg.norm(ef)
+
+    # Unit tangents
+    te = ee / norm_e
+    tf = ef / norm_f
+
+    # Curvature binormal
+    kb = 2.0 * np.cross(te, tf) / (1.0 + np.dot(te, tf))
+
+    chi = 1.0 + np.dot(te, tf)
+    tilde_t = (te + tf) / chi
+    tilde_d2 = (m2e + m2f) / chi
+
+    # Curvature
+    kappa1 = kb[2]
+
+    # Gradient of kappa1 with respect to edge vectors
+    Dkappa1De = 1.0 / norm_e * (-kappa1 * tilde_t + np.cross(tf, tilde_d2))
+    Dkappa1Df = 1.0 / norm_f * (-kappa1 * tilde_t - np.cross(te, tilde_d2))
+
+    # Populate the gradient of kappa
+    gradKappa[0:2] = -Dkappa1De[0:2]
+    gradKappa[2:4] = Dkappa1De[0:2] - Dkappa1Df[0:2]
+    gradKappa[4:6] = Dkappa1Df[0:2]
+
+    # Gradient of bending energy
+    dkappa = kappa1 - kappaBar
+    dF = gradKappa * EI * dkappa / l_k
+
+    return dF