diff --git a/app/HVSCMesh/area_adaptive/RB_adaptive.py b/app/HVSCMesh/area_adaptive/RB_adaptive.py index ed31115ff..6b223cafc 100644 --- a/app/HVSCMesh/area_adaptive/RB_adaptive.py +++ b/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)) diff --git a/app/HVSCMesh/test_case/mesh_quality_test.py b/app/HVSCMesh/test_case/mesh_quality_test.py index 66c5c7f88..f9f4ba828 100644 --- a/app/HVSCMesh/test_case/mesh_quality_test.py +++ b/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() + + diff --git a/app/soptx/linear_elasticity/hexahedrom_mesh.py b/app/soptx/linear_elasticity/hexahedrom_mesh.py index 3b1d77436..efbe03db1 100644 --- a/app/soptx/linear_elasticity/hexahedrom_mesh.py +++ b/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) diff --git a/example/experiment/fem/doublelaplace_fem_dirichlet.py b/example/experiment/fem/doublelaplace_fem_dirichlet.py index b62e6d1cf..958d20231 100644 --- a/example/experiment/fem/doublelaplace_fem_dirichlet.py +++ b/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)) diff --git a/example/experiment/fem/maxwell_eq_n1fem_dirichlet.py b/example/experiment/fem/maxwell_eq_n1fem_dirichlet.py index e1cd080d7..10df8b7d9 100644 --- a/example/experiment/fem/maxwell_eq_n1fem_dirichlet.py +++ b/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: diff --git a/example/experiment/fem/maxwell_eq_n2fem_dirichlet.py b/example/experiment/fem/maxwell_eq_n2fem_dirichlet.py index f6fffe86b..3719ac329 100644 --- a/example/experiment/fem/maxwell_eq_n2fem_dirichlet.py +++ b/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) diff --git a/example/experiment/fem/surface_poisson_lfem_example.py b/example/experiment/fem/surface_poisson_lfem_example.py index 37c78f901..c36e8814d 100644 --- a/example/experiment/fem/surface_poisson_lfem_example.py +++ b/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}次空间时间') diff --git a/example/iopt/Levy_example.py b/example/iopt/Levy_example.py deleted file mode 100644 index 5c92df173..000000000 --- a/example/iopt/Levy_example.py +++ /dev/null @@ -1,51 +0,0 @@ -from fealpy.backend import backend_manager as bm -from fealpy.opt.Levy import levy -import time -import matplotlib.pyplot as plt -from matplotlib.markers import MarkerStyle -# bm.set_backend('pytorch') - -start_time = time.perf_counter() - -# 定义设备 -device = 'cpu' -# device = 'cuda' - -n = 1000000 -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): - levy_steps = levy(n, 2, 1.5, num, device) - randn_steps = bm.random.randn(n, 2, num) - - levy_steps_tensor = bm.array(levy_steps, device=device, dtype=bm.float32) - randn_steps_tensor = bm.array(randn_steps, device=device, dtype=bm.float32) - - L[n * i : (i + 1) * n, :, :] = bm.cumsum(levy_steps_tensor, dim=0) + 10 - B[n * i : (i + 1) * n, :, :] = bm.cumsum(randn_steps_tensor, dim=0) + 10 - -end_time = time.perf_counter() -running_time = end_time - start_time -print("Running time:", running_time) - -# 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.show() \ No newline at end of file diff --git a/example/iopt/BFS_example.py b/example/opt/BFS_example.py similarity index 100% rename from example/iopt/BFS_example.py rename to example/opt/BFS_example.py diff --git a/example/iopt/DFS_example.py b/example/opt/DFS_example.py similarity index 100% rename from example/iopt/DFS_example.py rename to example/opt/DFS_example.py diff --git a/example/opt/Levy_example.py b/example/opt/Levy_example.py new file mode 100644 index 000000000..74e19c6fd --- /dev/null +++ b/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() \ No newline at end of file diff --git a/example/iopt/MultiAGV_example.py b/example/opt/MultiAGV_example.py similarity index 100% rename from example/iopt/MultiAGV_example.py rename to example/opt/MultiAGV_example.py diff --git a/example/iopt/README.md b/example/opt/README.md similarity index 100% rename from example/iopt/README.md rename to example/opt/README.md diff --git a/example/iopt/ant_example.py b/example/opt/ant_example.py similarity index 96% rename from example/iopt/ant_example.py rename to example/opt/ant_example.py index d80abf6d7..d1f42867b 100644 --- a/example/iopt/ant_example.py +++ b/example/opt/ant_example.py @@ -1,133 +1,133 @@ -import time -from fealpy.backend import backend_manager as bm -import matplotlib.pyplot as plt -from fealpy.iopt.ANT_TSP import calD, Ant_TSP -# bm.set_backend('pytorch') - -# 导入数据(34个城市) -citys = bm.array([ - [101.7400, 6.5600], - [112.5300, 37.8700], - [121.4700, 31.2300], - [119.3000, 26.0800], - [106.7100, 26.5700], - [103.7300, 36.0300], - [111.6500, 40.8200], - [120.1900, 30.2600], - [121.3000, 25.0300], - [106.5500, 29.5600], - [106.2700, 38.4700], - [116.4000, 39.9000], - [118.7800, 32.0400], - [114.1700, 22.3200], - [104.0600, 30.6700], - [108.9500, 34.2700], - [117.2000, 39.0800], - [117.2700, 31.8600], - [113.5400, 22.1900], - [102.7300, 25.0400], - [113.6500, 34.7600], - [123.3800, 41.8000], - [114.3100, 30.5200], - [113.2300, 23.1600], - [91.1100, 29.9700], - [117.0000, 36.6500], - [125.3500, 43.8800], - [113.0000, 28.2100], - [110.3500, 20.0200], - [87.6800, 43.7700], - [114.4800, 38.0300], - [126.6300, 45.7500], - [115.8900, 28.6800], - [108.3300, 22.8400] -]) - -# 距离矩阵 -D = calD(citys) - -# 初始化参数 -m = 10 # 蚂蚁数量 -alpha = 1 # 信息素重要程度因子 -beta = 5 # 启发函数重要程度因子 -rho = 0.5 # 信息素挥发因子 -Q = 1 # 常系数 -Eta = 1 / D # 启发函数 -iter_max = 100 # 最大迭代次数 - - -text1 = Ant_TSP(m, alpha, beta, rho, Q, Eta, D, iter_max) - -n = D.shape[0] -Tau = bm.ones((n, n), dtype=bm.float64) -Table = bm.zeros((m, n), dtype=int) # 路径记录表，每一行代表一个蚂蚁走过的路径 -Route_best = bm.zeros((iter_max, n), dtype=int) # 各代最佳路径 -Length_best = bm.zeros(iter_max) # 各代最佳路径的长度 - -# 设置循环次数 -l = 10 - -# 存储每次循环的结果 -shortest_lengths = [] -shortest_routes = [] - -# 循环迭代寻找最佳路径 -start_time = time.time() # 记录开始时间 - -for i in range(l): - # 迭代寻找最佳路径 - Length_best, Route_best = text1.cal(n, Tau, Table, Route_best, Length_best) - - # 结果显示 - Shortest_Length = bm.min(Length_best, keepdims = True) - index = bm.argmin(Length_best) - Shortest_Route = Route_best[index] - - # 在每次迭代结束后，将最短距离和最短路径添加到列表中 - shortest_lengths.append(Shortest_Length) - - SR_list = Shortest_Route.tolist() - result_list = SR_list + [SR_list[0]] - shortest_routes.append(bm.array(result_list)) - -end_time = time.time() # 记录结束时间 -elapsed_time = (end_time - start_time) / l # 计算运行时间 - -# 找到最短距离的索引 -min_distance_index = bm.argmin(bm.array(shortest_lengths)) - -# 根据索引找到最短路径 -Best_length = shortest_lengths[min_distance_index] -Best_path = shortest_routes[min_distance_index] - -# 找到最短距离迭代 -closest_index = bm.argmin(bm.abs(bm.array(shortest_lengths) - Best_length)) - -# 输出最短距离和对应路径 -print(f"运行时间: {elapsed_time}秒") -print('Shortest distance:', Best_length) -print('Shortest path:', Best_path) - -''' -# 绘图 -plt.figure(1) -x_values = bm.concatenate((citys[Best_path, 0], [citys[Best_path[0], 0]])) -y_values = bm.concatenate((citys[Best_path, 1], [citys[Best_path[0], 1]])) -plt.plot(x_values, y_values, 'o-') - -for i in range(citys.shape[0]): - plt.text(citys[i, 0], citys[i, 1], f' {i + 1}') -plt.scatter(citys[Best_path[0], 0], citys[Best_path[0], 1], color='red', s=100) -plt.xlabel('x of city location') -plt.ylabel('y of city location') -plt.title(f'Ant(shortest distance): {Best_length}\ntime: {elapsed_time}') # - -plt.figure(2) -plt.plot(range(iter_max), Length_best, 'b') -plt.legend(['shortest distance']) -plt.xlabel('number of iterations') -plt.ylabel('distance') -plt.title('changes of iterations') - -plt.show() - -''' +import time +from fealpy.backend import backend_manager as bm +import matplotlib.pyplot as plt +from fealpy.iopt.ANT_TSP import calD, Ant_TSP +# bm.set_backend('pytorch') + +# 导入数据(34个城市) +citys = bm.array([ + [101.7400, 6.5600], + [112.5300, 37.8700], + [121.4700, 31.2300], + [119.3000, 26.0800], + [106.7100, 26.5700], + [103.7300, 36.0300], + [111.6500, 40.8200], + [120.1900, 30.2600], + [121.3000, 25.0300], + [106.5500, 29.5600], + [106.2700, 38.4700], + [116.4000, 39.9000], + [118.7800, 32.0400], + [114.1700, 22.3200], + [104.0600, 30.6700], + [108.9500, 34.2700], + [117.2000, 39.0800], + [117.2700, 31.8600], + [113.5400, 22.1900], + [102.7300, 25.0400], + [113.6500, 34.7600], + [123.3800, 41.8000], + [114.3100, 30.5200], + [113.2300, 23.1600], + [91.1100, 29.9700], + [117.0000, 36.6500], + [125.3500, 43.8800], + [113.0000, 28.2100], + [110.3500, 20.0200], + [87.6800, 43.7700], + [114.4800, 38.0300], + [126.6300, 45.7500], + [115.8900, 28.6800], + [108.3300, 22.8400] +]) + +# 距离矩阵 +D = calD(citys) + +# 初始化参数 +m = 10 # 蚂蚁数量 +alpha = 1 # 信息素重要程度因子 +beta = 5 # 启发函数重要程度因子 +rho = 0.5 # 信息素挥发因子 +Q = 1 # 常系数 +Eta = 1 / D # 启发函数 +iter_max = 100 # 最大迭代次数 + + +text1 = Ant_TSP(m, alpha, beta, rho, Q, Eta, D, iter_max) + +n = D.shape[0] +Tau = bm.ones((n, n), dtype=bm.float64) +Table = bm.zeros((m, n), dtype=int) # 路径记录表，每一行代表一个蚂蚁走过的路径 +Route_best = bm.zeros((iter_max, n), dtype=int) # 各代最佳路径 +Length_best = bm.zeros(iter_max) # 各代最佳路径的长度 + +# 设置循环次数 +l = 10 + +# 存储每次循环的结果 +shortest_lengths = [] +shortest_routes = [] + +# 循环迭代寻找最佳路径 +start_time = time.time() # 记录开始时间 + +for i in range(l): + # 迭代寻找最佳路径 + Length_best, Route_best = text1.cal(n, Tau, Table, Route_best, Length_best) + + # 结果显示 + Shortest_Length = bm.min(Length_best, keepdims = True) + index = bm.argmin(Length_best) + Shortest_Route = Route_best[index] + + # 在每次迭代结束后，将最短距离和最短路径添加到列表中 + shortest_lengths.append(Shortest_Length) + + SR_list = Shortest_Route.tolist() + result_list = SR_list + [SR_list[0]] + shortest_routes.append(bm.array(result_list)) + +end_time = time.time() # 记录结束时间 +elapsed_time = (end_time - start_time) / l # 计算运行时间 + +# 找到最短距离的索引 +min_distance_index = bm.argmin(bm.array(shortest_lengths)) + +# 根据索引找到最短路径 +Best_length = shortest_lengths[min_distance_index] +Best_path = shortest_routes[min_distance_index] + +# 找到最短距离迭代 +closest_index = bm.argmin(bm.abs(bm.array(shortest_lengths) - Best_length)) + +# 输出最短距离和对应路径 +print(f"运行时间: {elapsed_time}秒") +print('Shortest distance:', Best_length) +print('Shortest path:', Best_path) + +''' +# 绘图 +plt.figure(1) +x_values = bm.concatenate((citys[Best_path, 0], [citys[Best_path[0], 0]])) +y_values = bm.concatenate((citys[Best_path, 1], [citys[Best_path[0], 1]])) +plt.plot(x_values, y_values, 'o-') + +for i in range(citys.shape[0]): + plt.text(citys[i, 0], citys[i, 1], f' {i + 1}') +plt.scatter(citys[Best_path[0], 0], citys[Best_path[0], 1], color='red', s=100) +plt.xlabel('x of city location') +plt.ylabel('y of city location') +plt.title(f'Ant(shortest distance): {Best_length}\ntime: {elapsed_time}') # + +plt.figure(2) +plt.plot(range(iter_max), Length_best, 'b') +plt.legend(['shortest distance']) +plt.xlabel('number of iterations') +plt.ylabel('distance') +plt.title('changes of iterations') + +plt.show() + +''' diff --git a/example/iopt/gridmap.py b/example/opt/gridmap.py similarity index 100% rename from example/iopt/gridmap.py rename to example/opt/gridmap.py diff --git a/example/iopt/pso_example.py b/example/opt/pso_example.py similarity index 100% rename from example/iopt/pso_example.py rename to example/opt/pso_example.py diff --git a/fealpy/experimental/tests/fem/test_scalar_mlaplace_source_integrator.py b/fealpy/experimental/tests/fem/test_scalar_mlaplace_source_integrator.py deleted file mode 100644 index 563e33770..000000000 --- a/fealpy/experimental/tests/fem/test_scalar_mlaplace_source_integrator.py +++ /dev/null @@ -1,42 +0,0 @@ -import ipdb -import numpy as np -import pytest -import sympy as sp -from fealpy.experimental.backend import backend_manager as bm -from fealpy.experimental.mesh.tetrahedron_mesh import TetrahedronMesh -from fealpy.experimental.functionspace.cm_conforming_fe_space3d import CmConformingFESpace3d -from fealpy.experimental.fem import LinearForm -from fealpy.experimental.fem.scalar_mlaplace_source_integrator import ScalarMLaplaceSourceIntegrator -from fealpy.experimental.tests.fem.mthlaplace_integrator_data import * -from fealpy.experimental.pde.biharmonic_triharmonic_3d import get_flist -class TestScalarMLaplaceSourceIntegrator: - @pytest.mark.parametrize("backend", ['numpy','pytorch']) - #@pytest.mark.parametrize("data", grad_m) - def test_scalar_mlaplace_source_integrator(self, backend): - bm.set_backend(backend) - def f(p): - x, y, z = p[..., 0], p[..., 1], p[..., 2] - return x**4*y**6+z**11 - x = sp.Symbol('x') - y = sp.Symbol('y') - z = sp.Symbol('z') - f = x**4*y**6+z**11 - flist = get_flist(f) - - mesh = TetrahedronMesh.from_box([0,1,0,1,0,1], nx=1, ny=1, nz=1) - node = mesh.entity('node') - isCornerNode = bm.zeros(len(node),dtype=bm.bool) - for n in bm.array([[0,0,0],[1,0,0],[0,1,0],[1,1,0],[1,1,1],[0,1,1],[0,0,1],[1,0,1]], dtype=bm.float64): - isCornerNode = isCornerNode | (bm.linalg.norm(node-n[None, :], axis=1)<1e-10) - space = CmConformingFESpace3d(mesh, p=11, m=1, isCornerNode=isCornerNode) - gdof = space.number_of_global_dofs() - lform = LinearForm(space) - integrator = ScalarMLaplaceSourceIntegrator(2, flist[2], q=14) - lform.add_integrator(integrator) - F = lform.assembly() -if __name__ == "__main__": - #pytest.main(['test_scalar_mlaplace_source_integrator.py', "-q", "-k","test_scalar_mlaplace_source_integrator", "-s"]) - t = TestScalarMLaplaceSourceIntegrator() - t.test_scalar_mlaplace_source_integrator("numpy") - - diff --git a/fealpy/fem/mlaplace_bernstein_integrator.py b/fealpy/fem/mlaplace_bernstein_integrator.py new file mode 100644 index 000000000..75ba9d774 --- /dev/null +++ b/fealpy/fem/mlaplace_bernstein_integrator.py @@ -0,0 +1,134 @@ +from typing import Optional +from scipy.special import factorial, comb +from fealpy.functionspace.functional import symmetry_index, symmetry_span_array + +from ..backend import backend_manager as bm +from ..typing import TensorLike, Index, _S + +from ..mesh import HomogeneousMesh +from ..functionspace.space import FunctionSpace as _FS +from ..utils import process_coef_func +from ..functional import bilinear_integral +from .integrator import ( + LinearInt, OpInt, CellInt, + enable_cache, + assemblymethod, + CoefLike +) + +class MLaplaceBernsteinIntegrator(LinearInt, OpInt, CellInt): + r"""The diffusion integrator for function spaces based on homogeneous meshes.""" + """ + m 求导次数 + """ + def __init__(self, m: int=2,coef: Optional[CoefLike]=None, q: + Optional[int]=None, *, + index: Index=_S, + batched: bool=False, + method: Optional[str]=None) -> None: + method = 'assembly' if (method is None) else method + super().__init__(method=method) + self.m = m + self.coef = coef + self.q = q + self.index = index + self.batched = batched + + @enable_cache + def to_global_dof(self, space: _FS) -> TensorLike: + return space.cell_to_dof()[self.index] + + @enable_cache + def fetch(self, space: _FS): + m = self.m + index = self.index + mesh = getattr(space, 'mesh', None) + + if not isinstance(mesh, HomogeneousMesh): + raise RuntimeError("The gradmIntegrator only support spaces on" + f"homogeneous meshes, but {type(mesh).__name__}is" + "not a subclass of HomoMesh.") + + cm = mesh.entity_measure('cell', index=index) + q = space.p+3 if self.q is None else self.q + qf = mesh.quadrature_formula(q, 'cell') + bcs, ws = qf.get_quadrature_points_and_weights() + gmphi = space.grad_m_basis(bcs, m=m) + return bcs, ws, gmphi, cm, index + + def assembly(self, space: _FS) -> TensorLike: + """ + @parm space: berstein space + """ + m = self.m + mesh = getattr(space, 'mesh', None) + p = space.p + + GD = mesh.geo_dimension() + NC = mesh.number_of_cells() + cellmeasure = mesh.entity_measure('cell') + multiIndex = mesh.multi_index_matrix(p-m, GD) + mmultiIndex = mesh.multi_index_matrix(m, GD) + pmultiIndex = mesh.multi_index_matrix(p, GD) + l = multiIndex.shape[0] # 矩阵的行数 + B = bm.zeros((l, l), dtype=space.ftype) + + # 积分公式 + def integrator(a): # a 为多重指标 + value = factorial(int(GD))*bm.prod(bm.array(factorial(bm.array(a))))/factorial(int(bm.sum(a)+GD)) + return value * cellmeasure[0] + + # 计算B^{\alpha-beta}矩阵 + for i in range(l): + for j in range(l): + alpha = multiIndex[i] + multiIndex[j] + c = (1/bm.prod(bm.array(factorial(multiIndex[i]))))*(1/bm.prod(bm.array(factorial(multiIndex[j])))) + B[i, j] = c * integrator(alpha) + + gmB = bm.zeros((NC, pmultiIndex.shape[0], pmultiIndex.shape[0]), dtype=space.ftype) + glambda = mesh.grad_lambda() + for beta in mmultiIndex: + idx1 = bm.where(bm.all(pmultiIndex-beta[None,:] >= 0, axis=1))[0] + multiidx1 = mesh.multi_index_matrix(bm.sum(beta), GD-1) + symidx, num = symmetry_index(GD, bm.sum(beta)) + + # \Lambda^{\beta} + glambda1 = symmetry_span_array(glambda, beta).reshape(NC, -1)[:, symidx] + idx11 = bm.broadcast_to(idx1[:, None], (l, l)) + c1 = factorial(int(bm.sum(beta)))/bm.prod(bm.array(factorial(beta))) + for theta in mmultiIndex: + idx2 = bm.where(bm.all(pmultiIndex-theta[None,:] >= 0, axis=1))[0] + glambda2 = symmetry_span_array(glambda, theta).reshape(NC, -1)[:, symidx] + c2 = factorial(int(bm.sum(theta)))/bm.prod(bm.array(factorial(theta))) + + # \Lambda^{\beta}: \Lambda^{\theta} + symglambda = bm.einsum('ci,ci,i->c', glambda1, glambda2, num) + idx22 = bm.broadcast_to(idx2[None, :], (l, l)) + B1 = B * c1 * c2 + BB = B1[None, :, :] * symglambda[:, None, None] + bm.add_at(gmB, (slice(None), idx11, idx22), BB) + gmB = gmB * factorial(int(p)) * factorial(int(p)) + return gmB + + + + + + + + + + + + + + + + + + + + + + + diff --git a/fealpy/fem/mthlaplace_integrator.py b/fealpy/fem/mthlaplace_integrator.py index eeae8eeba..3de360186 100644 --- a/fealpy/fem/mthlaplace_integrator.py +++ b/fealpy/fem/mthlaplace_integrator.py @@ -68,4 +68,19 @@ def assembly(self, space: _FS) -> TensorLike: #return bilinear_integral(gmphi1, gmphi, ws, cm, coef, # batched=self.batched) + @assemblymethod('without_numerical_integration') + def assembly_without_numerical_integration(self, space: _FS): + from .bilinear_form import BilinearForm + from .mlaplace_bernstein_integrator import MLaplaceBernsteinIntegrator + m = self.m + q = space.p+3 if self.q is None else self.q + cmcoeff = space.coeff + bgm = MLaplaceBernsteinIntegrator(m=m, q=q).assembly(space.bspace) + print(bgm.dtype) + print(cmcoeff.dtype) + M = bm.einsum('cil,clm,cpm->cip', cmcoeff, bgm, cmcoeff) + return M + + + diff --git a/fealpy/fem/utils.py b/fealpy/fem/utils.py index e7e241de1..baf062d5a 100644 --- a/fealpy/fem/utils.py +++ b/fealpy/fem/utils.py @@ -67,20 +67,4 @@ def shear_strain(gphi: TensorLike, indices: TensorLike, *, out: out = bm.set_at(out, (..., cursor, indices[:, j]), gphi[..., :, i]) cursor += 1 - # # TODO: Provide a unified implementation that is not backend-specific - # if bm.backend_name == 'numpy' or bm.backend_name == 'pytorch': - # for i in range(0, GD-1): - # for j in range(i+1, GD): - # out[..., cursor, indices[:, i]] = gphi[..., :, j] - # out[..., cursor, indices[:, j]] = gphi[..., :, i] - # cursor += 1 - # elif bm.backend_name == 'jax': - # for i in range(0, GD-1): - # for j in range(i+1, GD): - # out = out.at[..., cursor, indices[:, i]].set(gphi[..., :, j]) - # out = out.at[..., cursor, indices[:, j]].set(gphi[..., :, i]) - # cursor += 1 - # else: - # raise NotImplementedError("Backend is not yet implemented.") - return out diff --git a/fealpy/functionspace/cm_conforming_fe_space3d.py b/fealpy/functionspace/cm_conforming_fe_space3d.py index 8116f0934..d302eb30d 100644 --- a/fealpy/functionspace/cm_conforming_fe_space3d.py +++ b/fealpy/functionspace/cm_conforming_fe_space3d.py @@ -8,7 +8,7 @@ from ..backend import TensorLike from ..mesh.mesh_base import Mesh from .space import FunctionSpace -from .bernstein_fe_space import BernsteinFESpace +from fealpy.functionspace.bernstein_fe_space import BernsteinFESpace from .functional import symmetry_span_array, symmetry_index, span_array from scipy.special import factorial, comb from fealpy.decorator import barycentric @@ -260,6 +260,20 @@ def cell_to_dof(self): ## cell c2dof[:, ldof-cidof:] = c2id return c2dof + def is_boundary_dof(self, threshold=None): #TODO:threshold 未实现 + mesh = self.mesh + gdof = self.number_of_global_dofs() + isbdnode = mesh.boundary_node_flag() + isbdedge = mesh.boundary_edge_flag() + isbdFace = mesh.boundary_face_flag() + bd1 = self.node_to_dof()[isbdnode] + bd2 = self.edge_to_internal_dof()[isbdedge] + bd3 = self.face_to_internal_dof()[isbdFace] + bd = bm.concatenate([bd1.flatten(), bd2.flatten(), bd3.flatten()]) + isBdDof = bm.zeros(gdof, dtype=bm.bool) + isBdDof[bd] = True + return isBdDof + def get_corner(self): """ @brief 获取角点, 角边, 不太对啊 @@ -540,10 +554,10 @@ def coefficient_matrix(self): coeff[:, 4*ndof:4*ndof+6*edof] = bm.einsum('cji, cjk->cik',coeff2, coeff[:, 4*ndof:4*ndof+6*edof]) return coeff[:, : , dof2num] - def basis(self, bcs): + def basis(self, bcs, index=_S): coeff = self.coeff bphi = self.bspace.basis(bcs) - return bm.einsum('cil, cql -> cqi', coeff, bphi) + return bm.einsum('cil, cql -> cqi', coeff, bphi)[:,:,index] def grad_m_basis(self, bcs, m): coeff = self.coeff @@ -705,7 +719,58 @@ def interpolation(self, flist): S3 = self.dof_index["cell"] fI[c2id] = bcoeff[:, S3] return fI + + def boundary_interpolate(self, gD, uh, threshold=None): + mesh = self.mesh + m = self.m + p = self.p + isCornerNode = self.isCornerNode + isBdNode = mesh.boundary_node_flag() + isBdEdge = mesh.boundary_edge_flag() + isBdFace = mesh.boundary_face_flag() + + node = mesh.entity('node')[isBdNode] + edge = mesh.entity('edge')[isBdEdge] + face = mesh.entity('face')[isBdFace] + + NN = len(node) + NE = len(edge) + NF = len(face) + + nodeframe, edgeframe, faceframe = self.get_frame() + nodeframe = nodeframe[isBdNode] + edgeframe = edgeframe[isBdEdge] + faceframe = faceframe[isBdFace] + + n2id = self.node_to_dof()[isBdNode] + e2id = self.edge_to_internal_dof()[isBdEdge] + f2id = self.face_to_internal_dof()[isBdFace] + # 顶点 + uh[n2id[:, 0]] = gD[0](node) + k = 1 + for r in range(1, 4*m+1): + val = gD[r](node) + symidx, num = symmetry_index(3, r) + #bdnidxmap = + idx = bdnidxmap[coridx] + #print(idx) + return + + + + + + + + + + + + + + + diff --git a/fealpy/functionspace/first_nedelec_fe_space_3d.py b/fealpy/functionspace/first_nedelec_fe_space_3d.py index d0d2904f2..b1802915f 100644 --- a/fealpy/functionspace/first_nedelec_fe_space_3d.py +++ b/fealpy/functionspace/first_nedelec_fe_space_3d.py @@ -24,6 +24,7 @@ def __init__(self, mesh, p): self.multiindex = mesh.multi_index_matrix(p,3) self.ftype = mesh.ftype self.itype = mesh.itype + self.device = mesh.device def number_of_local_dofs(self, doftype='all'): p = self.p @@ -74,7 +75,7 @@ def face_to_dof(self): edge = self.mesh.entity('edge') face = self.mesh.entity('face') - f2d = bm.zeros((NF, fdof), dtype=self.itype) + f2d = bm.zeros((NF, fdof), device=self.device,dtype=self.itype) #f2d[:, :eldof*3] = e2dof[f2e].reshape(NF, eldof*3) f2d = bm.set_at(f2d,(slice(None),slice(eldof*3)),e2dof[f2e].reshape(NF, eldof*3)) s = [1, 0, 0] @@ -224,7 +225,7 @@ def basis(self, bcs, index=_S): c2esign = mesh.cell_to_edge_sign() #(NC, 3, 2) - l = bm.zeros((4, )+bcs[None,: ,0,None, None].shape, dtype=self.ftype) + l = bm.zeros((4, )+bcs[None,: ,0,None, None].shape,device=self.device, dtype=self.ftype) l = bm.set_at(l,(0),bcs[None, :,0,None, None]) l = bm.set_at(l,(1),bcs[None, :,1,None, None]) @@ -300,7 +301,7 @@ def face_internal_basis(self, bcs, index=_S): fdof = self.dof.number_of_local_dofs("face") glambda = mesh.grad_face_lambda() - val = bm.zeros((NF,) + bcs.shape[:-1]+(fdof, 3), dtype=self.ftype) + val = bm.zeros((NF,) + bcs.shape[:-1]+(fdof, 3),device=self.device, dtype=self.ftype) l = bm.zeros((3, )+bcs[None, :,0, None, None].shape, dtype=self.ftype) l = bm.set_at(l,(0),bcs[None, :,0, None, None]) @@ -338,7 +339,7 @@ def face_basis(self, bcs, index=_S): # edge basis phi = self.bspace.basis(bcs, p=p) multiIndex = self.mesh.multi_index_matrix(p, 2) - val = bm.zeros((NF,) + bcs.shape[:-1]+(ldof, 3), dtype=self.ftype) + val = bm.zeros((NF,) + bcs.shape[:-1]+(ldof, 3),device=self.device, dtype=self.ftype) for i in range(3): phie = phi[:, :, multiIndex[:, i]==0] c2esi = c2esign[:, i] @@ -384,7 +385,7 @@ def curl_basis(self, bcs): phi = self.bspace.basis(bcs, p=p) gphi = self.bspace.grad_basis(bcs, p=p) multiIndex = self.mesh.multi_index_matrix(p, 3) - val = bm.zeros((NC,)+bcs.shape[:-1]+(ldof, 3), dtype=self.ftype) + val = bm.zeros((NC,)+bcs.shape[:-1]+(ldof, 3),device=self.device, dtype=self.ftype) # edge basis locEdgeDual = bm.tensor([[2, 3], [1, 3], [1, 2], [0, 3], [0, 2], [0, 1]]) diff --git a/fealpy/functionspace/parametric_lagrange_fe_space.py b/fealpy/functionspace/parametric_lagrange_fe_space.py index f219cb96c..5968b5c5b 100644 --- a/fealpy/functionspace/parametric_lagrange_fe_space.py +++ b/fealpy/functionspace/parametric_lagrange_fe_space.py @@ -26,6 +26,7 @@ def __init__(self, mesh: _MT, p, q=None, spacetype='C'): self.dof = LinearMeshCFEDof(mesh, p) self.multi_index_matrix = mesh.multi_index_matrix + self.device = mesh.device self.GD = mesh.geo_dimension() self.TD = mesh.top_dimension() diff --git a/fealpy/functionspace/second_nedelec_fe_space_2d.py b/fealpy/functionspace/second_nedelec_fe_space_2d.py index c66edd5d6..8da8e07cf 100644 --- a/fealpy/functionspace/second_nedelec_fe_space_2d.py +++ b/fealpy/functionspace/second_nedelec_fe_space_2d.py @@ -96,7 +96,7 @@ def cell_to_dof(self): c2e = self.mesh.cell_to_edge() c2esign = self.mesh.cell_to_face_sign() - c2d = bm.zeros((NC, ldof), dtype=self.itype) + c2d = bm.zeros((NC, ldof), dtype=bm.int64) # c2d[:, e2ldof] : (NC, 3, p+1), e2dof[c2e] : (NC, 3, p+1) tmp = e2dof[c2e] # tmp[~c2esign] = tmp[~c2esign, ::-1] diff --git a/fealpy/functionspace/second_nedelec_fe_space_3d.py b/fealpy/functionspace/second_nedelec_fe_space_3d.py index c54402f6b..883ffea9e 100644 --- a/fealpy/functionspace/second_nedelec_fe_space_3d.py +++ b/fealpy/functionspace/second_nedelec_fe_space_3d.py @@ -26,7 +26,7 @@ def __init__(self, mesh, p): self.multiindex3 = mesh.multi_index_matrix(p,3) self.ftype = mesh.ftype self.itype = mesh.itype - + self.device=mesh.device def edge_to_local_face_dof(self): multiindex = self.multiindex2 @@ -53,7 +53,7 @@ def face_to_local_dof(self): ldof = self.number_of_local_dofs() fdof = self.number_of_local_dofs('faceall') - f2ld = bm.zeros((4, fdof), dtype=self.itype) + f2ld = bm.zeros((4, fdof),device=self.device, dtype=self.itype) eldof = self.edge_to_local_face_dof() nldof = bm.tensor([[eldof[(i+1)%3, 0], eldof[(i+2)%3, -1]] for i in range(3)]) @@ -91,10 +91,10 @@ def cell_to_dof(self): cdof = self.number_of_local_dofs('cell') fdof = self.number_of_local_dofs('faceall') - isndof = bm.zeros(ldof, dtype=bm.bool) + isndof = bm.zeros(ldof, device=self.device, dtype=bm.bool) isndof[f2ld] = True - c2d = bm.zeros((NC, ldof), dtype=bm.int64)# + c2d = bm.zeros((NC, ldof), device=self.device, dtype=bm.int64)# idx = bm.zeros((NC, 3), dtype=self.itype) fe = bm.tensor([[0, 1], [0, 2], [1, 2]], dtype=self.itype) #局部边 for i in range(4): @@ -235,7 +235,7 @@ def basis(self, bc,index=_S): c2v = self.basis_vector() #(NC, ldof, GD) shape = bc.shape[:-1] - val = bm.zeros((NC,)+ shape + (ldof, GD), dtype=self.ftype) + val = bm.zeros((NC,)+ shape + (ldof, GD),device=self.device, dtype=self.ftype) bval = self.lspace.basis(bc) #(NC, NQ, ldof//3) @@ -260,7 +260,7 @@ def face_basis(self, bc, index=_S): f2v = self.face_basis_vector(index=index)#(NF, ldof, GD) NF = len(f2v) shape = bc.shape[:-1] - val = bm.zeros((NF,) + shape+(ldof, GD), dtype=self.ftype) + val = bm.zeros((NF,) + shape+(ldof, GD), device=self.device, dtype=self.ftype) bval = self.lspace.basis(bc) #(NF, NQ, ldof//3) @@ -293,7 +293,7 @@ def face_basis_vector(self, index=_S): NF = len(f2e) - f2v = bm.zeros((NF, ldof, GD), dtype=self.ftype) + f2v = bm.zeros((NF, ldof, GD),device=self.device, dtype=self.ftype) # f2v[:, :ldof//2] = e2t[:, 0, None] #(NF, ldof//2, 3) # f2v[:, ldof//2:] = e2n[:, 0, None] f2v = bm.set_at(f2v,(slice(None),slice(None,ldof//2)),e2t[:, 0, None]) @@ -340,7 +340,7 @@ def face_dof_vector(self, index=_S): NF = len(f2e) - f2v = bm.zeros((NF, ldof, GD), dtype=self.ftype) + f2v = bm.zeros((NF, ldof, GD), device=self.device,dtype=self.ftype) # f2v[:, :ldof//2] = e2t[:, 0, None] #(NF, ldof//2, 3) # f2v[:, ldof//2:] = e2n[:, 0, None] f2v = bm.set_at(f2v,(slice(None),slice(None,ldof//2)),e2t[:, 0, None]) @@ -366,7 +366,7 @@ def curl_basis(self, bc): c2v = self.basis_vector()#(NC, ldof, GD) sgval = self.lspace.grad_basis(bc) #(NQ, NC, lldof, GD) - val = bm.zeros((NC,)+(bc.shape[0], )+(ldof, GD), dtype=self.ftype) + val = bm.zeros((NC,)+(bc.shape[0], )+(ldof, GD),device=self.device, dtype=self.ftype) # val[..., :ldof//3, :] = bm.cross(sgval, c2v[:,None, :ldof//3, :]) # val[..., ldof//3:2*(ldof//3), :] = bm.cross(sgval, c2v[:,None, ldof//3:2*(ldof//3), :]) @@ -544,7 +544,7 @@ def set_dirichlet_bc(self, gD, uh, threshold=None, q=None): uh[face2dof[:, ldof//2:]] = bm.sum(gval*f2v[:, ldof//2:], axis=-1) # uh[face2dof[:, :ldof//2]] =0 # uh[face2dof[:, ldof//2:]] =0 - isDDof = bm.zeros(gdof, dtype=bm.bool) + isDDof = bm.zeros(gdof, device=self.device,dtype=bm.bool) isDDof[face2dof] = True return uh,isDDof @@ -564,7 +564,7 @@ def set_neumann_bc(self,gD): points = mesh.bc_to_point(bcs)[isbdFace] n = mesh.face_unit_normal()[isbdFace] hval = gD(points,n) - vec = bm.zeros(gdof, dtype=self.ftype) + vec = bm.zeros(gdof, device=self.device,dtype=self.ftype) k = bm.einsum('fqg, fqlg,q,f->fl', hval, bphi,ws,fm) # (NF, ldof) bm.add.at(vec,face2dof,k) #bm.scatter_add(vec,face2dof,k) diff --git a/fealpy/geometry/functional.py b/fealpy/geometry/functional.py index b20b76549..e4f7ec270 100644 --- a/fealpy/geometry/functional.py +++ b/fealpy/geometry/functional.py @@ -118,3 +118,68 @@ def find_cut_point(phi: Callable, p0: TensorLike, p1: TensorLike) -> TensorLike: flag = flag[continue_signal] return pc + + +def project(imfun, p0:TensorLike, maxit=200, tol=1e-13, returngrad=False, returnd=False): + + eps = bm.finfo(float).eps + p = p0 + value = imfun(p) + s = bm.sign(value) + grad = imfun.gradient(p) + lg = bm.sum(grad**2, axis=-1, keepdims=True) + grad /= lg + grad *= value[..., None] + pp = p - grad + v = s[..., None]*(pp - p0) + d = bm.sqrt(bm.sum(v**2, axis=-1, keepdims=True)) + d *= s[..., None] + + g = imfun.gradient(pp) + g /= bm.sqrt(bm.sum(g**2, axis=-1, keepdims=True)) + g *= d + p = p0 - g + + k = 0 + while True: + value = imfun(p) + grad = imfun.gradient(p) + lg = bm.sqrt(bm.sum(grad**2, axis=-1, keepdims=True)) + grad /= lg + + v = s[..., None]*(p0 - p) + d = bm.sqrt(bm.sum(v**2, axis=-1)) + isOK = d < eps + d[isOK] = 0 + v[isOK] = grad[isOK] + v[~isOK] /= d[~isOK][..., None] + d *= s + + ev = grad - v + e = bm.max(bm.sqrt((value/lg.reshape(lg.shape[0:-1]))**2 + bm.sum(ev**2, axis=-1))) + if e < tol: + break + else: + k += 1 + if k > maxit: + break + grad /= lg + grad *= value[..., None] + pp = p - grad + v = s[..., None]*(pp - p0) + d = bm.sqrt(bm.sum(v**2, axis=-1, keepdims=True)) + d *= s[..., None] + + g = imfun.gradient(pp) + g /= bm.sqrt(bm.sum(g**2, axis=-1, keepdims=True)) + g *= d + p = p0 - g + + if (returnd is True) and (returngrad is True): + return p, d, grad + elif (returnd is False) and (returngrad is True): + return p, grad + elif (returnd is True) and (returngrad is False): + return p, d + else: + return p diff --git a/fealpy/geometry/implicit_surface.py b/fealpy/geometry/implicit_surface.py new file mode 100644 index 000000000..6014da4f9 --- /dev/null +++ b/fealpy/geometry/implicit_surface.py @@ -0,0 +1,151 @@ +from typing import Optional +from ..backend import backend_manager as bm +from ..backend import TensorLike +from .geometry_base import GeometryBase +from .functional import project + +class SphereSurface(): + def __init__(self, center=[0.0, 0.0, 0.0], radius=1.0): + self.center = center + self.radius = radius + r = radius + radius/10 + x = center[0] + y = center[1] + z = center[2] + self.box = [x-r, x+r, y-r, y+r, z-r, z+r] + + def __call__(self, *args): + if len(args) == 1: + p, = args + x = p[..., 0] + y = p[..., 1] + z = p[..., 2] + elif len(args) == 3: + x, y, z = args + else: + raise ValueError("the args must be a N*3 or (X, Y, Z)") + + cx, cy, cz = self.center + r = self.radius + return bm.sqrt((x - cx)**2 + (y - cy)**2 + (z - cz)**2) - r + + def gradient(self, p): + l = bm.sqrt(bm.sum((p - self.center)**2, axis=-1)) + n = (p - self.center)/l[..., None] + return n + + def unit_normal(self, p): + return self.gradient(p) + + def hessian(self, p): + x = p[..., 0] + y = p[..., 1] + z = p[..., 2] + shape = p.shape[:-1]+(3, 3) + H = bm.zeros(shape, dtype=bm.float64) + L = bm.sqrt(bm.sum(p*p, axis=-1)) + L3 = L**3 + H[..., 0, 0] = 1/L-x**2/L3 + H[..., 0, 1] = -x*y/L3 + H[..., 1, 0] = H[..., 0, 1] + H[..., 0, 2] = - x*z/L3 + H[..., 2, 0] = H[..., 0, 2] + H[..., 1, 1] = 1/L - y**2/L3 + H[..., 1, 2] = -y*z/L3 + H[..., 2, 1] = H[..., 1, 2] + H[..., 2, 2] = 1/L - z**2/L3 + return H + + def jacobi_matrix(self, p): + H = self.hessian(p) + n = self.unit_normal(p) + p[:], d = self.project(p) + + J = -(d[..., None, None]*H + bm.einsum('...ij, ...ik->...ijk', n, n)) + J[..., range(3), range(3)] += 1 + return J + + def tangent_operator(self, p): + pass + + def project(self, p, maxit=200, tol=1e-8): + d = self(p) + p = p - d[..., None]*self.unit_normal(p) + return p, d + + def init_mesh(self, meshtype='tri', returnnc=False, p=None): + if meshtype == 'tri': + t = (bm.sqrt(5) - 1)/2 + node = bm.array([ + [ 0, 1, t], + [ 0, 1,-t], + [ 1, t, 0], + [ 1,-t, 0], + [ 0,-1,-t], + [ 0,-1, t], + [ t, 0, 1], + [-t, 0, 1], + [ t, 0,-1], + [-t, 0,-1], + [-1, t, 0], + [-1,-t, 0]], dtype=bm.float64) + cell = bm.array([ + [6, 2, 0], + [3, 2, 6], + [5, 3, 6], + [5, 6, 7], + [6, 0, 7], + [3, 8, 2], + [2, 8, 1], + [2, 1, 0], + [0, 1,10], + [1, 9,10], + [8, 9, 1], + [4, 8, 3], + [4, 3, 5], + [4, 5,11], + [7,10,11], + [0,10, 7], + [4,11, 9], + [8, 4, 9], + [5, 7,11], + [10,9,11]], dtype=bm.int32) + node, d = self.project(node) + if returnnc: + return node, cell + else: + if p is None: + from fealpy.mesh.triangle_mesh import TriangleMesh + return TriangleMesh(node, cell) + else: + from fealpy.old.mesh.backup import LagrangeTriangleMesh + return LagrangeTriangleMesh(node, cell, p=p, surface=self) + + elif meshtype == 'quad': + node = bm.array([ + (-1, -1, -1), + (-1, -1, 1), + (-1, 1, -1), + (-1, 1, 1), + (1, -1, -1), + (1, -1, 1), + (1, 1, -1), + (1, 1, 1)], dtype=bm.float64) + cell = bm.array([ + (0, 1, 4, 5), + (6, 7, 2, 3), + (2, 3, 0, 1), + (4, 5, 6, 7), + (1, 3, 5, 7), + (2, 0, 6, 4)], dtype=bm.int32) + node, d = self.project(node) + if returnnc: + return node, cell + else: + if p is None: + from fealpy.mesh.quadrangle_mesh import QuadrangleMesh + return QuadrangleMesh(node, cell) + else: + from fealpy.old.mesh.backup import LagrangeQuadrangleMesh + return LagrangeQuadrangleMesh(node, cell, p=p, surface=self) + diff --git a/fealpy/mesh/lagrange_quadrangle_mesh.py b/fealpy/mesh/lagrange_quadrangle_mesh.py new file mode 100644 index 000000000..6f3da3207 --- /dev/null +++ b/fealpy/mesh/lagrange_quadrangle_mesh.py @@ -0,0 +1,76 @@ +from typing import Union, Optional, Sequence, Tuple, Any + +from ..backend import backend_manager as bm +from ..typing import TensorLike, Index, _S +from .. import logger + +from .utils import estr2dim +from .mesh_base import TensorMesh +from .quadrangle_mesh import QuadrangleMesh + +class LagrangeQuadrangleMesh(TensorMesh): + def __init__(self, node: TensorLike, cell: TensorLike, p=1, surface=None, + construct=False): + super().__init__(TD=2, itype=cell.dtype, ftype=node.dtype) + + kwargs = bm.context(cell) + self.p = p + self.node = node + self.cell = cell + self.surface = surface + + self.localEdge = self.generate_local_lagrange_edges(p) + self.localFace = self.localEdge + self.ccw = bm.tensor([0, 1, 2, 3], **kwargs) + self.localCell = None + + if construct: + self.construct() + + self.meshtype = 'lquad' + + + def generate_local_lagrange_edges(self, p: int) -> TensorLike: + """ + Generate the local edges for Lagrange elements of order p. + """ + TD = self.top_dimension() + multiIndex = bm.multi_index_matrix(p, TD) + + localEdge = bm.zeros((4, p+1), dtype=bm.int32) + localEdge[3, :], = bm.where(multiIndex[:, 2] == 0) + localEdge[2,:] = bm.flip(bm.where(multiIndex[:, 1] == 0)[1]) + localEdge[1,:] = bm.flip(bm.where(multiIndex[:, 1] == 0)[0]) + localEdge[0, :], = bm.where(multiIndex[:, 0] == 0) + + return localEdge + + def reference_cell_measure(self): + return 1 + + + @classmethod + def from_quadrangle_mesh(cls, mesh, p: int, surface=None): + node = mesh.interpolation_points(p) + cell = mesh.cell_to_ipoint(p) + if surface is not None: + node, _ = surface.project(node) + + lmesh = cls(node, cell, p=p, construct=True) + + lmesh.edge2cell = mesh.edge2cell # (NF, 4) + lmesh.cell2edge = mesh.cell_to_edge() + #lmesh.edge = mesh.edge_to_ipoint(p) + return lmesh + + def quadrature_formula(self, q, etype: Union[int, str] = 'cell'): + from ..quadrature import GaussLegendreQuadrature, TensorProductQuadrature + if isinstance(etype, str): + etype = estr2dim(self, etype) + qf = GaussLegendreQuadrature(q, dtype=self.ftype, device=self.device) + if etype == 2: + return TensorProductQuadrature((qf, qf)) + elif etype == 1: + return qf + else: + raise ValueError(f"entity type: {etype} is wrong!") diff --git a/fealpy/mesh/triangle_mesh.py b/fealpy/mesh/triangle_mesh.py index 09759b450..83090d503 100644 --- a/fealpy/mesh/triangle_mesh.py +++ b/fealpy/mesh/triangle_mesh.py @@ -68,8 +68,11 @@ def quadrature_formula(self, q: int, etype: Union[int, str]='cell', kwargs = {'dtype': self.ftype, 'device': self.device} if etype == 2: - from ..quadrature import TriangleQuadrature - quad = TriangleQuadrature(q, **kwargs) + from ..quadrature.stroud_quadrature import StroudQuadrature + if q > 9: + quad = StroudQuadrature(2, q) + else: + quad = TriangleQuadrature(q, **kwargs) elif etype == 1: from ..quadrature import GaussLegendreQuadrature quad = GaussLegendreQuadrature(q, **kwargs) diff --git a/fealpy/opt/Levy.py b/fealpy/opt/Levy.py index b83eb3859..1dd6cf9b2 100644 --- a/fealpy/opt/Levy.py +++ b/fealpy/opt/Levy.py @@ -5,7 +5,7 @@ def levy(n, m, beta, Num, device): num = gamma(1 + beta) * bm.sin(bm.array(bm.pi * beta / 2)) den = gamma((1 + beta) / 2) * beta * 2 ** ((beta - 1) / 2) sigma_u = (num / den) ** (1 / beta) - u = bm.random.randn([n, m, Num], device=device) * sigma_u - v = bm.random.randn([n, m, Num], device=device) + u = bm.random.randn(n, m, Num, device=device) * sigma_u + v = bm.random.randn(n, m, Num, device=device) z = u / (abs(v) ** (1 / beta)) return z \ No newline at end of file diff --git a/test/fem/test_mlaplace_berstesin_integrator.py b/test/fem/test_mlaplace_berstesin_integrator.py new file mode 100644 index 000000000..67c6d3810 --- /dev/null +++ b/test/fem/test_mlaplace_berstesin_integrator.py @@ -0,0 +1,34 @@ +import ipdb +import numpy as np +import pytest + +from fealpy.backend import backend_manager as bm +from fealpy.mesh import TriangleMesh +from fealpy.mesh.tetrahedron_mesh import TetrahedronMesh +from fealpy.functionspace.cm_conforming_fe_space3d import CmConformingFESpace3d +from fealpy.fem import BilinearForm +from fealpy.fem.mlaplace_bernstein_integrator import MLaplaceBernsteinIntegrator + +class TestMLaplaceBernsteinIntegrator: + @pytest.mark.parametrize("backend", ['numpy','pytorch']) + #@pytest.mark.parametrize("data", grad_m) + def test_mlaplace_berstein_integrator(self, backend, data): + bm.set_backend(backend) + mesh = TetrahedronMesh.from_box([0,1,0,1,0,1],1,1,1) + node = mesh.entity('node') + isCornerNode = np.zeros(len(node), dtype=np.bool) + for n in np.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],[1,1,0],[0,1,1],[1,0,1],[1,1,1]], dtype=np.float64): + isCornerNode = isCornerNode | (np.linalg.norm(node - n[None, :], axis=1) < 1e-10) + space = CmConformingFESpace3d(mesh, 11, 1, isCornerNode) + + bform = BilinearForm(space) + bform.add_integrator(MLaplaceBernsteinIntegrator(m=2, coef=1,q=14)) + A = bform.assembly() + + + + +if __name__ == "__main__": + t = TestMLaplaceBernsteinIntegrator() + #t.test_mlaplace_berstein_integrator('pytorch', None) + t.test_mlaplace_berstein_integrator('numpy', None) diff --git a/test/fem/test_mthlaplace_integrator.py b/test/fem/test_mthlaplace_integrator.py index e82512d6c..e7a52948c 100644 --- a/test/fem/test_mthlaplace_integrator.py +++ b/test/fem/test_mthlaplace_integrator.py @@ -10,6 +10,7 @@ from mthlaplace_integrator_data import * + class TestgradmIntegrator: @pytest.mark.parametrize("backend", ['numpy','pytorch']) @pytest.mark.parametrize("data", grad_m) @@ -34,5 +35,37 @@ def test_grad_m_integrator(self, backend, data): np.testing.assert_allclose(bm.to_numpy(FM), data["FM"], atol=1e-14) #np.testing.assert_allclose(bm.to_numpy(M.toarray()), data["M"], atol=1e-14) + @pytest.mark.parametrize("backend", ['numpy','pytorch']) + @pytest.mark.parametrize("data", grad_m) + def test_assembly_without_numerical_mlaplace_integrator(self, backend, data): + from fealpy.mesh.tetrahedron_mesh import TetrahedronMesh + from fealpy.functionspace.cm_conforming_fe_space3d import CmConformingFESpace3d + from fealpy.fem.mthlaplace_integrator import MthLaplaceIntegrator + + bm.set_backend(backend) + mesh = TetrahedronMesh.from_box([0,1,0,1,0,1], nx=1, ny=1, nz=1) + node = mesh.entity('node') + isCornerNode = bm.zeros(len(node),dtype=bm.bool) + for n in bm.array([[0,0,0],[1,0,0],[0,1,0],[1,1,0],[1,1,1],[0,1,1],[0,0,1],[1,0,1]], dtype=bm.float64): + isCornerNode = isCornerNode | (bm.linalg.norm(node-n[None, :], axis=1)<1e-10) + space = CmConformingFESpace3d(mesh, p=9, m=1, isCornerNode=isCornerNode) + + bform = BilinearForm(space) + integrator = MthLaplaceIntegrator(m=2, coef=1, q=14, method='without_numerical_integration') + bform.add_integrator(integrator) + M = bform.assembly() + + bform = BilinearForm(space) + integrator = MthLaplaceIntegrator(m=2, coef=1, q=14) + bform.add_integrator(integrator) + M1 = bform.assembly() + + + np.testing.assert_allclose(bm.to_numpy(M.toarray()), bm.to_numpy(M1.toarray()), atol=1e-14) + print('11111') + if __name__ == "__main__": - pytest.main(['test_grad_m_integrator.py', "-q", "-k","test_grad_m_integrator", "-s"]) + #pytest.main(['test_grad_m_integrator.py', "-q", "-k","test_grad_m_integrator", "-s"]) + t = TestgradmIntegrator() + t.test_assembly_without_numerical_mlaplace_integrator('numpy', grad_m[0]) + t.test_assembly_without_numerical_mlaplace_integrator('pytorch', grad_m[0]) diff --git a/test/fem/test_scalar_mlaplace_source_integrator.py b/test/fem/test_scalar_mlaplace_source_integrator.py index 16ba1516e..fe37b7098 100644 --- a/test/fem/test_scalar_mlaplace_source_integrator.py +++ b/test/fem/test_scalar_mlaplace_source_integrator.py @@ -7,7 +7,7 @@ from fealpy.functionspace.cm_conforming_fe_space3d import CmConformingFESpace3d from fealpy.fem import LinearForm from fealpy.fem.scalar_mlaplace_source_integrator import ScalarMLaplaceSourceIntegrator -#from fealpy.tests.fem.mthlaplace_integrator_data import * +from fealpy.tests.fem.mthlaplace_integrator_data import * from fealpy.pde.biharmonic_triharmonic_3d import get_flist class TestScalarMLaplaceSourceIntegrator: @pytest.mark.parametrize("backend", ['numpy','pytorch']) diff --git a/test/functionspace/test_cm_fe_space_3d.py b/test/functionspace/test_cm_fe_space_3d.py index 6c81ef4a4..9d3bc907e 100644 --- a/test/functionspace/test_cm_fe_space_3d.py +++ b/test/functionspace/test_cm_fe_space_3d.py @@ -5,7 +5,7 @@ from fealpy.mesh.tetrahedron_mesh import TetrahedronMesh from fealpy.quadrature.stroud_quadrature import StroudQuadrature from fealpy.functionspace.cm_conforming_fe_space3d import CmConformingFESpace3d -from fealpy.test.functionspace.cm_fe_space_3d_data import * +from cm_fe_space_3d_data import * from fealpy.pde.biharmonic_triharmonic_3d import get_flist import sympy as sp import ipdb @@ -86,7 +86,7 @@ def test_coefficient_matrix(self, data, backend): isCornerNode = np.zeros(len(node), dtype=np.bool) for n in np.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],[1,1,0],[0,1,1],[1,0,1],[1,1,1]], dtype=np.float64): isCornerNode = isCornerNode | (np.linalg.norm(node - n[None, :], axis=1) < 1e-10) - space = CmConformingFESpace3d(mesh, 9, 1, isCornerNode) + space = CmConformingFESpace3d(mesh, 11, 1, isCornerNode) coefficient_matrix = space.coefficient_matrix() #(6, 364, 364) np.testing.assert_allclose(coefficient_matrix[0, 180], data["cell0"], atol=1e-14) @@ -172,7 +172,7 @@ def test_source_vector(self, backend): for n in np.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],[1,1,0],[0,1,1],[1,0,1],[1,1,1]], dtype=np.float64): isCornerNode = isCornerNode | (np.linalg.norm(node - n[None, :], axis=1) < 1e-10) space = CmConformingFESpace3d(mesh, 9, 1, isCornerNode) - from fealpy.experimental.pde.biharmonic_triharmonic_3d import get_flist + from fealpy.pde.biharmonic_triharmonic_3d import get_flist x = sp.symbols('x') y = sp.symbols('y') z = sp.symbols('z') @@ -185,8 +185,8 @@ def f(p): y = p[..., 1] z = p[..., 2] return x**5*y**4+z**9 - from fealpy.experimental.fem import BilinearForm - from fealpy.experimental.fem import LinearForm, ScalarSourceIntegrator, ScalarMassIntegrator + from fealpy.fem import BilinearForm + from fealpy.fem import LinearForm, ScalarSourceIntegrator, ScalarMassIntegrator lform = LinearForm(space) lform.add_integrator(ScalarSourceIntegrator(f, q=14)) F = lform.assembly() @@ -210,17 +210,17 @@ def test_matrix(self, backend): for n in np.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],[1,1,0],[0,1,1],[1,0,1],[1,1,1]], dtype=np.float64): isCornerNode = isCornerNode | (np.linalg.norm(node - n[None, :], axis=1) < 1e-10) space = CmConformingFESpace3d(mesh, 11, 1, isCornerNode) - from fealpy.experimental.pde.biharmonic_triharmonic_3d import get_flist + from fealpy.pde.biharmonic_triharmonic_3d import get_flist x = sp.symbols('x') y = sp.symbols('y') z = sp.symbols('z') f = x**7*y**2+z**2 get_flist = get_flist(f) - from fealpy.experimental.fem import BilinearForm - from fealpy.experimental.fem import LinearForm - from fealpy.experimental.fem.scalar_mlaplace_source_integrator import ScalarMLaplaceSourceIntegrator - from fealpy.experimental.fem.mthlaplace_integrator import MthLaplaceIntegrator + from fealpy.fem import BilinearForm + from fealpy.fem import LinearForm + from fealpy.fem.scalar_mlaplace_source_integrator import ScalarMLaplaceSourceIntegrator + from fealpy.fem.mthlaplace_integrator import MthLaplaceIntegrator lform = LinearForm(space) lform.add_integrator(ScalarMLaplaceSourceIntegrator(2, get_flist[2], q=14)) gF = lform.assembly() @@ -236,6 +236,26 @@ def test_matrix(self, backend): #np.testing.assert_allclose(aa, gF, atol=1e-14,rtol=0) + @pytest.mark.parametrize("backend", ["numpy", "pytorch"]) + #@pytest.mark.parametrize("data", interpolation) + def test_boundary_interpolation(self, backend): + bm.set_backend(backend) + mesh = TetrahedronMesh.from_box([0,1,0,1,0,1],1,1,1) + node = mesh.entity('node') + isCornerNode = np.zeros(len(node), dtype=np.bool) + for n in np.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],[1,1,0],[0,1,1],[1,0,1],[1,1,1]], dtype=np.float64): + isCornerNode = isCornerNode | (np.linalg.norm(node - n[None, :], axis=1) < 1e-10) + space = CmConformingFESpace3d(mesh, 11, 1, isCornerNode) + uh = space.function() + from fealpy.pde.biharmonic_triharmonic_3d import get_flist + x = sp.symbols('x') + y = sp.symbols('y') + z = sp.symbols('z') + f = x**7*y**2+z**2 + gD = get_flist(f) + + boundary_interpolation = space.boundary_interpolate(gD, uh) + @@ -255,5 +275,6 @@ def test_matrix(self, backend): #t.test_basis(coefficient_matrix, "numpy") #t.test_interpolation(interpolation, "numpy") #t.test_interpolation(interpolation, "pytorch") - t.test_source_vector("pytorch") - t.test_matrix("pytorch") + #t.test_source_vector("pytorch") + #t.test_matrix("pytorch") + t.test_boundary_interpolation("numpy")