def METHOD_NAME(self, x):
    GD = self.mesh.geo_dimension()
    node0 = self.mesh.entity('node')
    node = np.full_like(node0, 0.0)
    node[self.isBdNode, :] = node0[self.isBdNode, :]
    NI = np.sum(~self.isBdNode)
    node[~self.isBdNode,0] = x[:NI]
    node[~self.isBdNode,1] = x[NI:]
    cell = self.mesh.entity('cell')
    NC = self.mesh.number_of_cells() 
    localEdge = self.mesh.ds.local_edge()
    v = [node[cell[:,j],:] - node[cell[:,i],:] for i,j in localEdge]
    l2 = np.zeros((NC, 3))
    for i in range(3):
        l2[:, i] = np.sum(v[i]**2, axis=1)
    l = np.sqrt(l2)
    p = l.sum(axis=1)
    q = l.prod(axis=1)
    area = np.cross(v[2], -v[1])/2
    METHOD_NAME = p*q/(16*area**2)
    A, B = self.grad_matrix(node=node)
    gradp = np.full_like(x, 0.0).reshape(GD, -1)
    gradp[0, :] = (A@node[:, 0] + B@node[:, 1])[~self.isBdNode]
    gradp[1, :] = (B.T@node[:, 0] + A@node[:, 1])[~self.isBdNode]
    return np.mean(METHOD_NAME), gradp.flat