```from __future__ import division

import numpy as np

import pdb

import scipy.sparse
from scipy.sparse import coo_matrix
import cvxopt
import cvxopt.cholmod

from utils import deleterowcol

class Problem(object):

@staticmethod
def lk(E=1.):
"""element stiffness matrix"""
nu = 0.3
k = np.array([0.5 - nu / 6., 0.125 + nu / 8., -0.25 - nu / 12.,
-0.125 + 0.375 * nu, -0.25 + nu / 12., -0.125 - nu / 8., nu / 6.,
0.125 - 0.375 * nu])
KE = E / (1 - nu**2) * np.array([
[k[0], k[1], k[2], k[3], k[4], k[5], k[6], k[7]],
[k[1], k[0], k[7], k[6], k[5], k[4], k[3], k[2]],
[k[2], k[7], k[0], k[5], k[6], k[3], k[4], k[1]],
[k[3], k[6], k[5], k[0], k[7], k[2], k[1], k[4]],
[k[4], k[5], k[6], k[7], k[0], k[1], k[2], k[3]],
[k[5], k[4], k[3], k[2], k[1], k[0], k[7], k[6]],
[k[6], k[3], k[4], k[1], k[2], k[7], k[0], k[5]],
[k[7], k[2], k[1], k[4], k[3], k[6], k[5], k[0]]])
return KE

def build_indices(self, nelx, nely):
""" FE: Build the index vectors for the for coo matrix format. """
self.KE = self.lk()
self.edofMat = np.zeros((nelx * nely, 8), dtype=int)
for elx in range(nelx):
for ely in range(nely):
el = ely + elx * nely
n1 = (nely + 1) * elx + ely
n2 = (nely + 1) * (elx + 1) + ely
self.edofMat[el, :] = np.array([2 * n1 + 2, 2 * n1 + 3,
2 * n2 + 2, 2 * n2 + 3, 2 * n2, 2 * n2 + 1, 2 * n1,
2 * n1 + 1])
# Construct the index pointers for the coo format
self.iK = np.kron(self.edofMat, np.ones((8, 1))).flatten()
self.jK = np.kron(self.edofMat, np.ones((1, 8))).flatten()

def __init__(self, nelx, nely, penal, bc):
# Problem size
self.nelx = nelx
self.nely = nely

# Max and min stiffness
self.Emin = 1e-9
self.Emax = 1.0

# SIMP penalty
self.penal = penal

# dofs:
self.ndof = 2 * (nelx + 1) * (nely + 1)

# FE: Build the index vectors for the for coo matrix format.
self.build_indices(nelx, nely)

# BC's and support (half MBB-beam)
dofs = np.arange(2 * (nelx + 1) * (nely + 1))
self.fixed = bc.get_fixed_nodes()
self.free = np.setdiff1d(dofs, self.fixed)

# Solution and RHS vectors
self.f = bc.get_forces()
self.u = np.zeros(self.f.shape)

# Per element compliance
self.ce = np.zeros(nely * nelx)

def compute_displacements(self, xPhys):
# Setup and solve FE problem
sK = ((self.KE.flatten()[np.newaxis]).T * (
self.Emin + (xPhys)**self.penal *
(self.Emax - self.Emin))).flatten(order='F')
K = scipy.sparse.coo_matrix((sK, (self.iK, self.jK)),
shape=(self.ndof, self.ndof)).tocsc()
# Remove constrained dofs from matrix and convert to coo
K = deleterowcol(K, self.fixed, self.fixed).tocoo()
# Solve system
K1 = cvxopt.spmatrix(K.data, K.row.astype(np.int), K.col.astype(np.int))
B = cvxopt.matrix(self.f[self.free, :])
cvxopt.cholmod.linsolve(K1, B)
self.u[self.free, :] = np.array(B)[:, :]

def compute_compliance(self, xPhys, dc):
# Compute compliance and its gradient
u = self.u[:, 0][self.edofMat].reshape(-1, 8)
self.ce[:] = (u.dot(self.KE) * u).sum(1)
obj = ((self.Emin + xPhys**self.penal *
(self.Emax - self.Emin)) * self.ce).sum()
dc[:] = (-self.penal * xPhys**(self.penal - 1.0) *
(self.Emax - self.Emin)) * self.ce
return obj
```