fix bug

test/opt/test_matrix_vector_product_gradient_optimizer.py

@@ -1,9 +1,13 @@
 import numpy as np
+import pytest
+import matplotlib.pyplot as plt
 from scipy.sparse.linalg import LinearOperator, cg
 from scipy.sparse import csr_matrix
-import pytest
+
 from fealpy.mesh import TriangleMesh
-from fealpy.opt import Problem,MatrixVectorProductGradientOptimizer
+from fealpy.opt import Problem, MatrixVectorProductGradientOptimizer
+
+import ipdb
 
 class TriMeshProblem(Problem):
     def __init__(self,mesh:TriangleMesh):
@@ -20,7 +24,6 @@ def quality(self, x):
         node = np.full_like(node0, 0.0)
 
         node[self.isBdNode, :] = node0[self.isBdNode, :]
-        #node[~self.isBdNode, :].T.flat[:] = x
         NI = np.sum(~self.isBdNode)
         node[~self.isBdNode,0] = x[:NI]
         node[~self.isBdNode,1] = x[NI:]
@@ -60,7 +63,6 @@ def grad_matrix(self, node=None):
         v0 = node[idxk] - node[idxj]
         v1 = node[idxi] - node[idxk]
         v2 = node[idxj] - node[idxi]
-
         area = 0.5*(-v2[:, [0]]*v1[:, [1]] + v2[:, [1]]*v1[:, [0]])
         l2 = np.zeros((NC, 3), dtype=np.float64)
         l2[:, 0] = np.sum(v0**2, axis=1)
@@ -115,19 +117,36 @@ def block_jacobi_preconditioner(self, x):
         b = B*node[:, 0] - A*node1[:, 1]
         node1[isFreeNode, 1], info = cg(A[np.ix_(isFreeNode, isFreeNode)], b[isFreeNode], x0=node[isFreeNode, 1], tol=1e-6)
 
-        return (node1[isFreeNode, :] - node[isFreeNode, :]).T.flat
+        return np.array((node1[isFreeNode, :] - node[isFreeNode, :]).T.flat)
 
 def test_triangle_mesh_opt():
-    #mesh = TriangleMesh.from_unit_circle_gmsh(h=0.1)
-    mesh = TriangleMesh.from_one_triangle('equ')
-    mesh.uniform_refine(n=2)
+    mesh = TriangleMesh.from_unit_circle_gmsh(h=0.1)
+    #mesh = TriangleMesh.from_one_triangle('equ')
+    #mesh.uniform_refine(n=3)
     area = mesh.entity_measure('cell')
     problem = TriMeshProblem(mesh)
+
     NDof = len(problem.x0)
-    problem.Preconditioner = LinearOperator((NDof, NDof), problem.block_jacobi_preconditioner)
-    problem.StepLength = 0.3
+    problem.Preconditioner = LinearOperator((NDof, NDof), problem.block_jacobi_preconditioner, dtype=np.float64)
+    problem.StepLength = 1.0
     opt = MatrixVectorProductGradientOptimizer(problem)
-    opt.run()
+    x, f, g = opt.run()
+
+    fig = plt.figure()
+    axes = fig.gca()
+    mesh.add_plot(axes)
+
+    node = mesh.entity('node')
+    isFreeNode = ~mesh.ds.boundary_node_flag()
+    n = len(x)//2 
+    node[isFreeNode,0] = x[:n]
+    node[isFreeNode,1] = x[n:]
+
+    fig = plt.figure()
+    axes = fig.gca()
+    mesh.add_plot(axes)
+
+    plt.show()
 
 if __name__ == "__main__":
     test_triangle_mesh_opt()