Merge remote-tracking branch 'upstream/master'

app/HVSCMesh/area_adaptive/RB_adaptive.py

@@ -2,7 +2,7 @@
 import gmsh
 
 #from fealpy.mesh import TriangleMesh
-from fealpy.mesh import TriangleMesh
+from fealpy.old.mesh import TriangleMesh
 
 class Mesh:
     def __init__(self):
@@ -86,7 +86,7 @@ def __init__(self):
 #gmsh.option.setNumber("Mesh.MeshSizeMax", 20)
 gmsh.model.mesh.optimize('', True)
 gmsh.model.mesh.generate(2)
-#gmsh.fltk.run()
+gmsh.fltk.run()
 
 ntags, vxyz, _ = gmsh.model.mesh.getNodes()
 node = vxyz.reshape((-1,3))

app/HVSCMesh/test_case/mesh_quality_test.py

@@ -3,6 +3,7 @@
 
 from fealpy.mesh import TriangleMesh
 from fealpy.mesh.mesh_quality import RadiusRatioQuality
+from app.HVSCMesh.radius_ratio_objective import RadiusRatioSumObjective
 
 node = np.array([[0.0,0.0],[2.0,0.0],[1,np.sqrt(3)]],dtype=np.float64)
 cell = np.array([[0,1,2]],dtype=np.int_)
@@ -14,15 +15,27 @@
 node[cell[-1,1]] = node[cell[-1,1]]+[-0.1,0.15]
 node[cell[-1,2]] = node[cell[-1,2]]+[0,-0.15]
 
-localEdge = mesh.localEdge
-v = [node[cell[:,j],:] - node[cell[:,i],:] for i,j in localEdge]
-area = np.cross(v[2],-v[1])/2
-
+# mesh_quality test
 mesh_quality = RadiusRatioQuality(mesh) 
 
 quality = mesh_quality(node)
-print(quality)
+grad = mesh_quality.jac(node)
+A,B = mesh_quality.hess(node)
+
+# radius_ratio_objective test
+sum_quality = RadiusRatioSumObjective(mesh_quality)
+quality = sum_quality(node)
+grad = sum_quality.jac(node)
+A,B = sum_quality.hess(node)
+
+#print(quality)
+#print(grad)
+#print(A.toarray())
+#print(B.toarray())
+
 fig = plt.figure()
 axes = fig.gca()
 mesh.add_plot(axes,aspect=1)
 plt.show()
+
+

app/soptx/linear_elasticity/hexahedrom_mesh.py

@@ -108,23 +108,23 @@ def dirichlet(self, points: TensorLike) -> TensorLike:
                             in arbitrary order Lagrange finite element space on HexahedronMesh.")
 parser.add_argument('--backend',
                     choices=['numpy', 'pytorch'], 
-                    default='numpy', type=str,
+                    default='pytorch', type=str,
                     help='Specify the backend type for computation, default is "pytorch".')
 parser.add_argument('--degree', 
-                    default=1, type=int, 
+                    default=2, type=int, 
                     help='Degree of the Lagrange finite element space, default is 1.')
 parser.add_argument('--solver',
                     choices=['cg', 'spsolve'],
                     default='cg', type=str,
                     help='Specify the solver type for solving the linear system, default is "cg".')
 parser.add_argument('--nx', 
-                    default=2, type=int, 
+                    default=8, type=int, 
                     help='Initial number of grid cells in the x direction, default is 2.')
 parser.add_argument('--ny',
-                    default=2, type=int,
+                    default=8, type=int,
                     help='Initial number of grid cells in the y direction, default is 2.')
 parser.add_argument('--nz',
-                    default=2, type=int,
+                    default=8, type=int,
                     help='Initial number of grid cells in the z direction, default is 2.')
 args = parser.parse_args()
 
@@ -134,14 +134,14 @@ def dirichlet(self, points: TensorLike) -> TensorLike:
 bm.set_backend(args.backend)
 
 nx, ny, nz = args.nx, args.ny, args.nz
-mesh = HexahedronMesh.from_box(box=pde.domain(), nx=nx, ny=ny, nz=nz, device='cpu')
+mesh = HexahedronMesh.from_box(box=pde.domain(), nx=nx, ny=ny, nz=nz, device='cuda')
 
 p = args.degree
 
 tmr = timer("FEM Solver")
 next(tmr)
 
-maxit = 4
+maxit = 1
 errorType = ['$|| u  - u_h ||_{L2}$', '$|| u -  u_h||_{l2}$']
 errorMatrix = bm.zeros((len(errorType), maxit), dtype=bm.float64)
 NDof = bm.zeros(maxit, dtype=bm.int32)

example/experiment/fem/doublelaplace_fem_dirichlet.py

@@ -62,8 +62,6 @@
 u = (sp.sin(2*sp.pi*y)*sp.sin(2*sp.pi*x))**2
 pde = DoubleLaplacePDE(u) 
 ulist = get_flist(u)[:3]
-import ipdb
-ipdb.set_trace()
 mesh = TriangleMesh.from_box([0,1,0,1], n, n)
 NDof = bm.zeros(maxit, dtype=bm.float64)
 
@@ -94,6 +92,7 @@
     lform.add_integrator(ScalarSourceIntegrator(pde.source, q=p+4))
 
     A = bform.assembly()
+    #print(space.number_of_global_dofs())
     F = lform.assembly()
     tmr.send(f'第{i}次矩组装时间')
 
@@ -105,7 +104,10 @@
     A, F = bc1.apply(A, F)  
     tmr.send(f'第{i}次边界处理时间')
     A = A.to_scipy()
-    print(type(A))
+
+    from numpy.linalg import cond
+    print(gdof)
+    print(cond(A.toarray()))
     #A = coo_matrix(A)
     #A = csr_matrix((A.values(), A.indices()),A.shape)
     uh[:] = bm.tensor(spsolve(A, F))

example/experiment/fem/maxwell_eq_n1fem_dirichlet.py

@@ -64,7 +64,9 @@
 args = parser.parse_args()
 p = args.degree
 maxit = args.maxit
+maxit = 3
 dim = args.dim
+dim = 3
 backend = args.backend
 bm.set_backend(backend)
 if dim == 2:

example/experiment/fem/maxwell_eq_n2fem_dirichlet.py

@@ -65,9 +65,11 @@
 args = parser.parse_args()
 p = args.degree
 maxit = args.maxit
+maxit = 3
 dim = args.dim
+dim = 3
 backend = args.backend
-bm.set_backend("pytorch")
+bm.set_backend(backend)
 if dim == 2:
     pde = PDE2d()
 else:
@@ -79,8 +81,8 @@
              '$||\\nabla \\times u - \\nabla_h \\times u_h||_0$ with k=2',
              '$||\\nabla \\times u - \\nabla_h \\times u_h||_0$ with k=3',
              '$||\\nabla \\times u - \\nabla_h \\times u_h||_0$ with k=4']
-errorMatrix = bm.zeros((len(errorType), maxit), dtype=bm.float64)
-NDof = bm.zeros(maxit, dtype=bm.float64)
+errorMatrix = np.zeros((len(errorType), maxit), dtype=np.float64)
+NDof = np.zeros(maxit, dtype=np.float64)
 
 tmr = timer()
 next(tmr)

example/experiment/fem/surface_poisson_lfem_example.py

@@ -10,7 +10,7 @@
 from fealpy.backend import backend_manager as bm
 
 from fealpy.pde.surface_poisson_model import SurfaceLevelSetPDEData
-from fealpy.geometry import SphereSurface
+from fealpy.geometry.implicit_surface import SphereSurface
 from fealpy.mesh.triangle_mesh import TriangleMesh
 from fealpy.mesh.lagrange_triangle_mesh import LagrangeTriangleMesh
 from fealpy.functionspace.parametric_lagrange_fe_space import ParametricLagrangeFESpace
@@ -63,8 +63,8 @@
 surface = SphereSurface()
 tmesh = TriangleMesh.from_unit_sphere_surface()
 mesh = LagrangeTriangleMesh.from_triangle_mesh(tmesh, p, surface=surface)
-#fname = f"sphere_test.vtu"
-#mesh.to_vtk(fname=fname)
+fname = f"sphere_test.vtu"
+mesh.to_vtk(fname=fname)
 
 space = ParametricLagrangeFESpace(mesh, p=sdegree)
 #tmr.send(f'第{i}次空间时间')

example/opt/Levy_example.py

@@ -0,0 +1,52 @@
+from fealpy.backend import backend_manager as bm
+from fealpy.opt.Levy import levy
+device = 'cpu'
+
+# 定义后端
+bm.set_backend('pytorch')
+# 定义设备
+device = 'cuda'
+
+import time
+start_time = time.perf_counter()
+
+n = 500000
+num = 100
+iter_num = 1
+
+L = bm.zeros(1 + n * iter_num, 2, num, device=device)
+B = bm.zeros(1 + n * iter_num, 2, num, device=device)
+L[0, :, :] = L[0, :, :] + 10
+B[0, :, :] = B[0, :, :] + 10
+
+for i in range(iter_num):
+    bround_steps = bm.random.randn(n, 2, num, device=device)
+    levy_steps = levy(n, 2, 1.5, num, device)
+
+    levy_steps_tensor = bm.array(levy_steps, device=device, dtype=bm.float32)
+    bround_steps_tensor = bm.array(bround_steps, device=device, dtype=bm.float32)
+
+    L[n * i : (i + 1) * n, :, :] = bm.cumsum(levy_steps_tensor, axis=0)
+    B[n * i : (i + 1) * n, :, :] = bm.cumsum(bround_steps_tensor, axis=0)
+
+end_time = time.perf_counter()
+running_time = end_time - start_time
+print("Running time:", running_time)
+
+import matplotlib.pyplot as plt
+from matplotlib.markers import MarkerStyle
+L_cpu = L.cpu().numpy()
+B_cpu = B.cpu().numpy()
+
+for particle in range(0, 1):
+    L_x_vals = L_cpu[:, 0, particle]  
+    L_y_vals = L_cpu[:, 1, particle]  
+    B_x_vals = B_cpu[:, 0, particle]  
+    B_y_vals = B_cpu[:, 1, particle]  
+    plt.figure()
+    plt.plot(L_x_vals, L_y_vals, '-', color = 'blue')
+    plt.plot(B_x_vals, B_y_vals, '-', color = 'black')
+    plt.scatter(L_x_vals[0], L_y_vals[0], marker=MarkerStyle('*'), color='red', s=500) 
+plt.xticks([])
+plt.yticks([])
+plt.show()