def METHOD_NAME(self, q=None):
    # 积分公式
    qf = self.integrator if q is None else self.mesh1.integrator(q, etype='cell')
    # bcs : (NQ, n)
    # ws : (NQ, )
    bcs, ws = qf.get_quadrature_points_and_weights() 
    rm = self.mesh1.reference_cell_measure() # 参考单元测度
    d = self.mesh1.first_fundamental_form(bcs)
    d = np.sqrt(np.linalg.det(d))
    d *= self.model.stress(bcs, self.ce, self.rho, self.gamma, self.mesh0, self.mesh1)
    gphi = self.cspace.grad_basis(bcs) # (NQ, NC, ldof, GD)
    phi = self.dspace.basis(bcs)
    M0 = np.einsum('q, qci, qcj, qc->cij', ws*rm, gphi[..., 0], phi, d) # (NC, ldof0, ldof1)
    M1 = np.einsum('q, qci, qcj, qc->cij', ws*rm, gphi[..., 1], phi, d) # (NC, ldof0, ldof1)
    c2d0 = self.cspace.cell_to_dof()
    c2d1 = self.dspace.cell_to_dof()
    I = np.broadcast_to(c2d0[:, :, None], shape=M0.shape)
    J = np.broadcast_to(c2d1[:, None, :], shape=M0.shape)
    gdof0 = self.cspace.number_of_global_dofs()
    gdof1 = self.dspace.number_of_global_dofs()
    M0 = csr_matrix(
            (M0.flat, (I.flat, J.flat)), 
            shape=(gdof0, gdof1)
            )
    M1 = csr_matrix(
            (M1.flat, (I.flat, J.flat)), 
            shape=(gdof0, gdof1)
            )
    return M0, M1