inv free compress square

Author: liminchen

inv_free/MassSpringEnergy.py

@@ -0,0 +1,175 @@
+import utils
+import numpy as np
+import math
+
+def polar_svd(F):
+    [U, s, VT] = np.linalg.svd(F)
+    if np.linalg.det(U) < 0:
+        U[:, 1] = -U[:, 1]
+        s[1] = -s[1]
+    if np.linalg.det(VT) < 0:
+        VT[1, :] = -VT[1, :]
+        s[1] = -s[1]
+    return [U, s, VT]
+
+def dPsi_div_dsigma(s, mu, lam):
+    ln_sigma_prod = math.log(s[0] * s[1])
+    inv0 = 1.0 / s[0]
+    dPsi_dsigma_0 = mu * (s[0] - inv0) + lam * inv0 * ln_sigma_prod
+    inv1 = 1.0 / s[1]
+    dPsi_dsigma_1 = mu * (s[1] - inv1) + lam * inv1 * ln_sigma_prod
+    return [dPsi_dsigma_0, dPsi_dsigma_1]
+
+def d2Psi_div_dsigma2(s, mu, lam):
+    ln_sigma_prod = math.log(s[0] * s[1])
+    inv2_0 = 1 / (s[0] * s[0])
+    d2Psi_dsigma2_00 = mu * (1 + inv2_0) - lam * inv2_0 * (ln_sigma_prod - 1)
+    inv2_1 = 1 / (s[1] * s[1])
+    d2Psi_dsigma2_11 = mu * (1 + inv2_1) - lam * inv2_1 * (ln_sigma_prod - 1)
+    d2Psi_dsigma2_01 = lam / (s[0] * s[1])
+    return [[d2Psi_dsigma2_00, d2Psi_dsigma2_01], [d2Psi_dsigma2_01, d2Psi_dsigma2_11]]
+
+def B_left_coef(s, mu, lam):
+    sigma_prod = s[0] * s[1]
+    return (mu + (mu - lam * math.log(sigma_prod)) / sigma_prod) / 2
+
+def Psi(F, mu, lam):
+    J = np.linalg.det(F)
+    lnJ = math.log(J)
+    return mu / 2 * (np.trace(np.transpose(F).dot(F)) - 2) - mu * lnJ + lam / 2 * lnJ * lnJ
+
+def dPsi_div_dF(F, mu, lam):
+    FinvT = np.transpose(np.linalg.inv(F))
+    return mu * (F - FinvT) + lam * math.log(np.linalg.det(F)) * FinvT
+
+def d2Psi_div_dF2(F, mu, lam):
+    [U, sigma, VT] = polar_svd(F)
+
+    Psi_sigma_sigma = utils.make_PD(d2Psi_div_dsigma2(sigma, mu, lam))
+
+    B_left = B_left_coef(sigma, mu, lam)
+    Psi_sigma = dPsi_div_dsigma(sigma, mu, lam)
+    B_right = (Psi_sigma[0] + Psi_sigma[1]) / (2 * max(sigma[0] + sigma[1], 1e-6))
+    B = utils.make_PD([[B_left + B_right, B_left - B_right], [B_left - B_right, B_left + B_right]])
+
+    M = np.array([[0, 0, 0, 0]] * 4)
+    M[0, 0] = Psi_sigma_sigma[0, 0]
+    M[0, 3] = Psi_sigma_sigma[0, 1]
+    M[1, 1] = B[0, 0]
+    M[1, 2] = B[0, 1]
+    M[2, 1] = B[1, 0]
+    M[2, 2] = B[1, 1]
+    M[3, 0] = Psi_sigma_sigma[1, 0]
+    M[3, 3] = Psi_sigma_sigma[1, 1]
+
+    dP_div_dF = np.array([[0, 0, 0, 0]] * 4)
+    for j in range(0, 2):
+        for i in range(0, 2):
+            ij = j * 2 + i
+            for s in range(0, 2):
+                for r in range(0, 2):
+                    rs = s * 2 + r
+                    dP_div_dF[ij, rs] = M[0, 0] * U[i, 0] * VT[0, j] * U[r, 0] * VT[0, s] \
+                        + M[0, 3] * U[i, 0] * VT[0, j] * U[r, 1] * VT[1, s] \
+                        + M[1, 1] * U[i, 0] * VT[1, j] * U[r, 0] * VT[1, s] \
+                        + M[1, 2] * U[i, 0] * VT[1, j] * U[r, 1] * VT[0, s] \
+                        + M[2, 1] * U[i, 1] * VT[0, j] * U[r, 0] * VT[1, s] \
+                        + M[2, 2] * U[i, 1] * VT[0, j] * U[r, 1] * VT[0, s] \
+                        + M[3, 0] * U[i, 1] * VT[1, j] * U[r, 0] * VT[0, s] \
+                        + M[3, 3] * U[i, 1] * VT[1, j] * U[r, 1] * VT[1, s]
+    return dP_div_dF
+
+def deformation_grad(x, elemVInd, IB):
+    F = [x[elemVInd[1]] - x[elemVInd[0]], x[elemVInd[2]] - x[elemVInd[0]]]
+    return np.transpose(F).dot(IB)
+
+def dPsi_div_dx(P, IB):  # applying chain-rule, dPsi_div_dx = dPsi_div_dF * dF_div_dx
+    dPsi_dx_2 = P[0, 0] * IB[0, 0] + P[0, 1] * IB[0, 1]
+    dPsi_dx_3 = P[1, 0] * IB[0, 0] + P[1, 1] * IB[0, 1]
+    dPsi_dx_4 = P[0, 0] * IB[1, 0] + P[0, 1] * IB[1, 1]
+    dPsi_dx_5 = P[1, 0] * IB[1, 0] + P[1, 1] * IB[1, 1]
+    return [np.array([-dPsi_dx_2 - dPsi_dx_4, -dPsi_dx_3 - dPsi_dx_5]), np.array([dPsi_dx_2, dPsi_dx_3]), np.array([dPsi_dx_4, dPsi_dx_5])]
+
+def d2Psi_div_dx2(dP_div_dF, IB):  # applying chain-rule, d2Psi_div_dx2 = dF_div_dx^T * d2Psi_div_dF2 * dF_div_dx (note that d2F_div_dx2 = 0)
+    intermediate = np.array([[0.0, 0.0, 0.0, 0.0]] * 6)
+    for colI in range(0, 4):
+        _000 = dP_div_dF[0, colI] * IB[0, 0]
+        _010 = dP_div_dF[0, colI] * IB[1, 0]
+        _101 = dP_div_dF[2, colI] * IB[0, 1]
+        _111 = dP_div_dF[2, colI] * IB[1, 1]
+        _200 = dP_div_dF[1, colI] * IB[0, 0]
+        _210 = dP_div_dF[1, colI] * IB[1, 0]
+        _301 = dP_div_dF[3, colI] * IB[0, 1]
+        _311 = dP_div_dF[3, colI] * IB[1, 1]
+        intermediate[2, colI] = _000 + _101
+        intermediate[3, colI] = _200 + _301
+        intermediate[4, colI] = _010 + _111
+        intermediate[5, colI] = _210 + _311
+        intermediate[0, colI] = -intermediate[2, colI] - intermediate[4, colI]
+        intermediate[1, colI] = -intermediate[3, colI] - intermediate[5, colI]
+    result = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] * 6)
+    for colI in range(0, 6):
+        _000 = intermediate[colI, 0] * IB[0, 0]
+        _010 = intermediate[colI, 0] * IB[1, 0]
+        _101 = intermediate[colI, 2] * IB[0, 1]
+        _111 = intermediate[colI, 2] * IB[1, 1]
+        _200 = intermediate[colI, 1] * IB[0, 0]
+        _210 = intermediate[colI, 1] * IB[1, 0]
+        _301 = intermediate[colI, 3] * IB[0, 1]
+        _311 = intermediate[colI, 3] * IB[1, 1]
+        result[2, colI] = _000 + _101
+        result[3, colI] = _200 + _301
+        result[4, colI] = _010 + _111
+        result[5, colI] = _210 + _311
+        result[0, colI] = -_000 - _101 - _010 - _111
+        result[1, colI] = -_200 - _301 - _210 - _311
+    return result
+
+def val(x, e, vol, IB, mu, lam):
+    sum = 0.0
+    for i in range(0, len(e)):
+        F = deformation_grad(x, e[i], IB[i])
+        sum += vol[i] * Psi(F, mu[i], lam[i])
+    return sum
+
+def grad(x, e, vol, IB, mu, lam):
+    g = np.array([[0.0, 0.0]] * len(x))
+    for i in range(0, len(e)):
+        F = deformation_grad(x, e[i], IB[i])
+        P = vol[i] * dPsi_div_dF(F, mu[i], lam[i])
+        g_local = dPsi_div_dx(P, IB[i])
+        for j in range(0, 3):
+            g[e[i][j]] += g_local[j]
+    return g
+
+def hess(x, e, vol, IB, mu, lam):
+    IJV = [[0] * (len(e) * 36), [0] * (len(e) * 36), np.array([0.0] * (len(e) * 36))]
+    for i in range(0, len(e)):
+        F = deformation_grad(x, e[i], IB[i])
+        dP_div_dF = vol[i] * d2Psi_div_dF2(F, mu[i], lam[i])
+        local_hess = d2Psi_div_dx2(dP_div_dF, IB[i])
+        for xI in range(0, 3):
+            for xJ in range(0, 3):
+                for dI in range(0, 2):
+                    for dJ in range(0, 2):
+                        ind = i * 36 + (xI * 3 + xJ) * 4 + dI * 2 + dJ
+                        IJV[0][ind] = e[i][xI] * 2 + dI
+                        IJV[1][ind] = e[i][xJ] * 2 + dJ
+                        IJV[2][ind] = local_hess[xI * 2 + dI, xJ * 2 + dJ]
+    return IJV
+
+def init_step_size(x, e, p):
+    alpha = 1
+    for i in range(0, len(e)):
+        x21 = x[e[i][1]] - x[e[i][0]]
+        x31 = x[e[i][2]] - x[e[i][0]]
+        p21 = p[e[i][1]] - p[e[i][0]]
+        p31 = p[e[i][2]] - p[e[i][0]]
+        detT = np.linalg.det(np.transpose([x21, x31]))
+        a = np.linalg.det(np.transpose([p21, p31])) / detT
+        b = (np.linalg.det(np.transpose([x21, p31])) + np.linalg.det(np.transpose([p21, x31]))) / detT
+        c = 0.9  # solve for alpha that first brings the new volume to 0.1x the old volume for slackness
+        critical_alpha = utils.smallest_positive_real_root_quad(a, b, c)
+        if critical_alpha > 0:
+            alpha = min(alpha, critical_alpha)
+    return alpha

inv_free/NeoHookeanEnergy.py

@@ -10,12 +10,13 @@
 # simulation setup
 side_len = 1
 rho = 1000      # density of square
-k = 2e4         # spring stiffness
+E = 1e5         # Young's modulus
+nu = 0.4        # Poisson's ratio
 n_seg = 4       # num of segments per side of the square
 h = 0.01        # time step size in s
 DBC = [(n_seg + 1) * (n_seg + 1)]       # dirichlet node index
 DBC_v = [np.array([0.0, -0.5])]         # dirichlet node velocity
-DBC_limit = [np.array([0.0, -0.6])]     # dirichlet node limit position
+DBC_limit = [np.array([0.0, -0.7])]     # dirichlet node limit position
 ground_n = np.array([0.0, 1.0])         # normal of the slope
 ground_n /= np.linalg.norm(ground_n)    # normalize ground normal vector just in case
 ground_o = np.array([0.0, -1.0])        # a point on the slope  
@@ -26,12 +27,15 @@
 x = np.append(x, [[0.0, side_len * 0.6]], axis=0)   # ceil origin (with normal [0.0, -1.0])
 v = np.array([[0.0, 0.0]] * len(x))                 # velocity
 m = [rho * side_len * side_len / ((n_seg + 1) * (n_seg + 1))] * len(x)  # calculate node mass evenly
-# rest length squared
-l2 = []
+# rest shape basis, volume, and lame parameters
+vol = [0.0] * len(e)
+IB = [np.array([[0.0, 0.0]] * 2)] * len(e)
 for i in range(0, len(e)):
-    diff = x[e[i][0]] - x[e[i][1]]
-    l2.append(diff.dot(diff))
-k = [k] * len(e)    # spring stiffness
+    TB = [x[e[i][1]] - x[e[i][0]], x[e[i][2]] - x[e[i][0]]]
+    vol[i] = np.linalg.det(np.transpose(TB)) / 2
+    IB[i] = np.linalg.inv(np.transpose(TB))
+mu_lame = [0.5 * E / (1 + nu)] * len(e)
+lam = [E * nu / ((1 + nu) * (1 - 2 * nu))] * len(e)
 # identify whether a node is Dirichlet
 is_DBC = [False] * len(x)
 for i in DBC:
@@ -64,14 +68,16 @@ def screen_projection(x):
         screen_projection([x[-1][0] - 3.0, x[-1][1]]))   # ceil
     for eI in e:
         pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[0]]), screen_projection(x[eI[1]]))
+        pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[1]]), screen_projection(x[eI[2]]))
+        pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[2]]), screen_projection(x[eI[0]]))
     for xId in range(0, len(x) - 1):
         xI = x[xId]
         pygame.draw.circle(screen, (0, 0, 255), screen_projection(xI), 0.1 * side_len / n_seg * scale)
 
     pygame.display.flip()   # flip the display
 
     # step forward simulation and wait for screen refresh
-    [x, v] = time_integrator.step_forward(x, e, v, m, l2, k, ground_n, ground_o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, 1e-2)
+    [x, v] = time_integrator.step_forward(x, e, v, m, vol, IB, mu_lame, lam, ground_n, ground_o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, 1e-2)
     time_step += 1
     pygame.time.wait(int(h * 1000))
 

inv_free/simulator.py

@@ -8,20 +8,16 @@ def generate(side_length, n_seg):
         for j in range(0, n_seg + 1):
             x[i * (n_seg + 1) + j] = [-side_length / 2 + i * step, -side_length / 2 + j * step]
 
-    # connect the nodes with edges
+    # connect the nodes with triangle elements
     e = []
-    # horizontal edges
-    for i in range(0, n_seg):
-        for j in range(0, n_seg + 1):
-            e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j])
-    # vertical edges
-    for i in range(0, n_seg + 1):
-        for j in range(0, n_seg):
-            e.append([i * (n_seg + 1) + j, i * (n_seg + 1) + j + 1])
-    # diagonals
     for i in range(0, n_seg):
         for j in range(0, n_seg):
-            e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1])
-            e.append([(i + 1) * (n_seg + 1) + j, i * (n_seg + 1) + j + 1])
+            # triangulate each cell following a symmetric pattern:
+            if (i % 2)^(j % 2):
+                e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j, i * (n_seg + 1) + j + 1])
+                e.append([(i + 1) * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1, i * (n_seg + 1) + j + 1])
+            else:
+                e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1])
+                e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1, i * (n_seg + 1) + j + 1])
 
     return [x, e]

inv_free/square_mesh.py

@@ -7,13 +7,13 @@
 from scipy.sparse.linalg import spsolve
 
 import InertiaEnergy
-import MassSpringEnergy
+import NeoHookeanEnergy
 import GravityEnergy
 import BarrierEnergy
 import FrictionEnergy
 import SpringEnergy
 
-def step_forward(x, e, v, m, l2, k, n, o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, tol):
+def step_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, tol):
     x_tilde = x + v * h     # implicit Euler predictive position
     x_n = copy.deepcopy(x)
     mu_lambda = BarrierEnergy.compute_mu_lambda(x, n, o, contact_area, mu)  # compute mu * lambda for each node using x^n
@@ -23,54 +23,54 @@ def step_forward(x, e, v, m, l2, k, n, o, contact_area, mu, is_DBC, DBC, DBC_v,
             DBC_target.append(x_n[DBC[i]] + h * DBC_v[i])
         else:
             DBC_target.append(x_n[DBC[i]])
-    DBC_stiff = 10  # initialize stiffness for DBC springs
+    DBC_stiff = 1000  # initialize stiffness for DBC springs
 
     # Newton loop
     iter = 0
-    E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
-    [p, DBC_satisfied] = search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
+    E_last = IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
+    [p, DBC_satisfied] = search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
     while (LA.norm(p, inf) / h > tol) | (sum(DBC_satisfied) != len(DBC)):   # also check whether all DBCs are satisfied
         print('Iteration', iter, ':')
         print('residual =', LA.norm(p, inf) / h)
 
         if (LA.norm(p, inf) / h <= tol) & (sum(DBC_satisfied) != len(DBC)):
             # increase DBC stiffness and recompute energy value record
             DBC_stiff *= 2
-            E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
+            E_last = IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
 
         # filter line search
-        alpha = BarrierEnergy.init_step_size(x, n, o, p)  # avoid interpenetration and tunneling
-        while IP_val(x + alpha * p, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h) > E_last:
+        alpha = min(BarrierEnergy.init_step_size(x, n, o, p), NeoHookeanEnergy.init_step_size(x, e, p))  # avoid interpenetration, tunneling, and inversion
+        while IP_val(x + alpha * p, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h) > E_last:
             alpha /= 2
         print('step size =', alpha)
 
         x += alpha * p
-        E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
-        [p, DBC_satisfied] = search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
+        E_last = IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
+        [p, DBC_satisfied] = search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
         iter += 1
 
     v = (x - x_n) / h   # implicit Euler velocity update
     return [x, v]
 
-def IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
+def IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
     return InertiaEnergy.val(x, x_tilde, m) + h * h * (     # implicit Euler
-        MassSpringEnergy.val(x, e, l2, k) + 
+        NeoHookeanEnergy.val(x, e, vol, IB, mu_lame, lam) + 
         GravityEnergy.val(x, m) + 
         BarrierEnergy.val(x, n, o, contact_area) + 
         FrictionEnergy.val(v, mu_lambda, h, n)
     ) + SpringEnergy.val(x, m, DBC, DBC_target, DBC_stiff)
 
-def IP_grad(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
+def IP_grad(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
     return InertiaEnergy.grad(x, x_tilde, m) + h * h * (    # implicit Euler
-        MassSpringEnergy.grad(x, e, l2, k) + 
+        NeoHookeanEnergy.grad(x, e, vol, IB, mu_lame, lam) + 
         GravityEnergy.grad(x, m) + 
         BarrierEnergy.grad(x, n, o, contact_area) + 
         FrictionEnergy.grad(v, mu_lambda, h, n)
     ) + SpringEnergy.grad(x, m, DBC, DBC_target, DBC_stiff)
 
-def IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
+def IP_hess(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
     IJV_In = InertiaEnergy.hess(x, x_tilde, m)
-    IJV_MS = MassSpringEnergy.hess(x, e, l2, k)
+    IJV_MS = NeoHookeanEnergy.hess(x, e, vol, IB, mu_lame, lam)
     IJV_B = BarrierEnergy.hess(x, n, o, contact_area)
     IJV_F = FrictionEnergy.hess(v, mu_lambda, h, n)
     IJV_S = SpringEnergy.hess(x, m, DBC, DBC_target, DBC_stiff)
@@ -84,9 +84,9 @@ def IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_
     H = sparse.coo_matrix((IJV[2], (IJV[0], IJV[1])), shape=(len(x) * 2, len(x) * 2)).tocsr()
     return H
 
-def search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h):
-    projected_hess = IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h)
-    reshaped_grad = IP_grad(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h).reshape(len(x) * 2, 1)
+def search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, v, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h):
+    projected_hess = IP_hess(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h)
+    reshaped_grad = IP_grad(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h).reshape(len(x) * 2, 1)
     # check whether each DBC is satisfied
     DBC_satisfied = [False] * len(x)
     for i in range(0, len(DBC)):