3.0.2 version update

.gitignore

@@ -41,6 +41,7 @@ dist
 *.png
 *.fig
 *.vtk
+*.vts
 *.so
 *.swn
 *.swo

app/PathPlanning/Grid_maps.py

@@ -1,4 +1,3 @@
-from yaml import MappingStartEvent
 from fealpy.backend import backend_manager as bm
 from fealpy.opt.optimizer_base import opt_alg_options
 import scipy.sparse as sp
@@ -9,26 +8,26 @@
 
 # 二值栅格地图数据
 def grid_0():
-    Map=bm.array([[0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0],
-              [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0],
-              [1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0],
-              [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
-              [0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-              [0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0],
-              [0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
-              [1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1],
-              [0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0],
-              [0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
-              [0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
-              [0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0],
-              [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
-              [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-              [0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
-              [0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1],
-              [0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1],
-              [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
-              [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1],
-              [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
+    Map = bm.array([[0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0],
+                    [1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
+                    [0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0],
+                    [0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
+                    [1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1],
+                    [0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0],
+                    [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
+                    [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1],
+                    [0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1],
+                    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
+                    [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1],
+                    [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
     return Map
 
 

app/PathPlanning/main.py

@@ -108,13 +108,13 @@ def __init__(self, num, NP=20, lb=0, ub=1, MaxIters = 50):
             'PSO': ParticleSwarmOpt,
             'GWO': GreyWolfOptimizer,
             'ACO': AntColonyOptAlg,
-            # 'HO': HippopotamusOptAlg,
+            'HO': HippopotamusOptAlg,
             # 'CPO': CrestedPorcupineOpt,
-            'BKA':BlackwingedKiteAlg,
-            'BOA':ButterflyOptAlg,
-            'CS':CuckooSearchOpt,
-            'DE':DifferentialEvolution,
-            'ETO':ExponentialTrigonometricOptAlg,
+            'BKA': BlackwingedKiteAlg,
+            'BOA': ButterflyOptAlg,
+            'CS': CuckooSearchOpt,
+            'DE': DifferentialEvolution,
+            'ETO': ExponentialTrigonometricOptAlg,
         }
         self.results = {}
 
@@ -143,7 +143,7 @@ def optimize(self):
             gbest, gbest_f = optimizer.run()
             end_time = time.perf_counter()
             running_time = end_time - start_time
-
+            print(gbest)
             result = textMAP.calresult(gbest)
             result["path"] = [x for x, y in zip(result["path"], result["path"][1:] + [None]) if x != y]
 
@@ -154,7 +154,7 @@ def optimize(self):
             }
             print(f'The best optimal value by {algo_name} is: {float(gbest_f):.4f}  Time: {end_time - start_time:.4f} seconds')
 
-        textMAP.printMAP(self.results)
+        # textMAP.printMAP(self.results)
 
 
 
@@ -166,5 +166,5 @@ def optimize(self):
     tsp_optimizer.optimize()
 
     #Shortest path problem
-    short_optimizer = ShortestPathApp(num)
-    short_optimizer.optimize()
+    # short_optimizer = ShortestPathApp(num)
+    # short_optimizer.optimize()

app/fracturex/fracturex/cases/phase_field/model0_example.py

@@ -6,7 +6,7 @@
 from fealpy.mesh import TriangleMesh
 from fealpy.old.geometry.domain_2d import SquareWithCircleHoleDomain
 
-from app.fracturex.fracturex.phasefield.main_solver import MainSolve
+from app.fracturex.fracturex.phasefield.main_solve import MainSolve
 
 from fealpy.utils import timer
 
@@ -133,7 +133,7 @@ def is_dirchlet_boundary(self, p):
 mesh = TriangleMesh.from_domain_distmesh(domain, maxit=100)
 
 
-ms = MainSolve(mesh=mesh, material_params=model.params, p=p)
+ms = MainSolve(mesh=mesh, material_params=model.params)
 tmr.send('init')
 
 # 拉伸模型边界条件
@@ -146,11 +146,11 @@ def is_dirchlet_boundary(self, p):
 if bm.backend_name == 'pytorch':
     ms.auto_assembly_matrix()
 
-
+#ms.set_scipy_solver()
 ms.output_timer()
 
 ms.save_vtkfile(fname=vtkname)
-ms.solve(maxit=maxit)
+ms.solve(p=p, maxit=maxit)
 
 tmr.send('stop')
 tmr.send(None)

app/fracturex/fracturex/cases/phase_field/model3d.py

@@ -1,7 +1,9 @@
 from fealpy.backend import backend_manager as bm
+import torch
+import numpy as np
 from fealpy.mesh import TetrahedronMesh, HexahedronMesh
 
-from app.fracturex.fracturex.phasefield.main_solver import MainSolver
+from app.fracturex.fracturex.phasefield.main_solve import MainSolve
 from fealpy.utils import timer
 
 import time
@@ -64,9 +66,10 @@ def is_dirchlet_boundary(self, p):
         default='HybridModel', type=str,
         help='有限元方法, 默认为 HybridModel.')
 
+
 parser.add_argument('--mesh_type',
-        default='tet', type=str,
-        help='网格类型, 默认为 tet.')
+        default='tri', type=str,
+        help='网格类型, 默认为 tri.')
 
 parser.add_argument('--enable_adaptive',
         default=False, type=bool,
@@ -81,8 +84,8 @@ def is_dirchlet_boundary(self, p):
         help='网格加密方法, 默认为 bisect.')
 
 parser.add_argument('--n',
-        default=2, type=int,
-        help='初始网格加密次数, 默认为 2.')
+        default=6, type=int,
+        help='初始网格加密次数, 默认为 6.')
 
 parser.add_argument('--vtkname',
         default='test', type=str,
@@ -92,6 +95,18 @@ def is_dirchlet_boundary(self, p):
         default=True, type=bool,
         help='是否保存 vtk 文件, 默认为 False.')
 
+parser.add_argument('--force_type',
+        default='y', type=str,
+        help='Force type, default is y.')
+
+parser.add_argument('--gpu', 
+        default=False, type=bool,
+        help='是否使用 GPU, 默认为 False.')
+
+parser.add_argument('--cupy', 
+        default=False, type=bool,
+        help='是否使用cupy求解.')
+
 args = parser.parse_args()
 p= args.degree
 maxit = args.maxit
@@ -101,15 +116,18 @@ def is_dirchlet_boundary(self, p):
 marking_strategy = args.marking_strategy
 refine_method = args.refine_method
 n = args.n
-vtkname = args.vtkname
 save_vtkfile = args.save_vtkfile
 vtkname = args.vtkname +'_' + args.mesh_type + '_'
-
+force_type = args.force_type
+gpu = args.gpu
+cupy = args.cupy
 
 tmr = timer()
 next(tmr)
 start = time.time()
 bm.set_backend(backend)
+if gpu:
+    bm.set_default_device('cuda')
 model = square_with_circular_notch_3d()
 
 if args.mesh_type == 'hex':
@@ -121,33 +139,43 @@ def is_dirchlet_boundary(self, p):
 
 mesh.uniform_refine(n=n)
 
+
 fname = args.mesh_type + '_3d_square_with_a_notch_init.vtu'
 mesh.to_vtk(fname=fname)
 
-ms = MainSolver(mesh=mesh, material_params=model.params, p=p, model_type=model_type)
+ms = MainSolve(mesh=mesh, material_params=model.params, model_type=model_type)
 tmr.send('init')
 
-if enable_adaptive:
-    print('Enable adaptive refinement.')
-    ms.set_adaptive_refinement(marking_strategy=marking_strategy, refine_method=refine_method)
-
-# 拉伸模型边界条件
 ms.add_boundary_condition('force', 'Dirichlet', model.is_force_boundary, model.is_z_force(), 'z')
 
 
 # 固定位移边界条件
 ms.add_boundary_condition('displacement', 'Dirichlet', model.is_dirchlet_boundary, 0)
 
-ms.solve(maxit=maxit, vtkname=vtkname)
+
+if bm.backend_name == 'pytorch':
+    ms.auto_assembly_matrix()
+if cupy:
+    ms.set_cupy_solver()
+
+ms.output_timer()
+ms.save_vtkfile(fname=vtkname)
+ms.solve(p=p, maxit=maxit)
 
 tmr.send('stop')
+tmr.send(None)
 end = time.time()
 
-force = ms.Rforce
-disp = ms.force_value
-tname = args.mesh_type + '_p' + str(p) + '_' + 'model3d_disp.txt'
+force = ms.get_residual_force()
+disp = model.is_force()
+
+ftname = 'force_'+args.mesh_type + '_p' + str(p) + '_' + 'model3d_disp.pt'
+
+torch.save(force, ftname)
+#np.savetxt('force'+tname, bm.to_numpy(force))
+tname = 'params_'+args.mesh_type + '_p' + str(p) + '_' + 'model3d_disp.txt'
 with open(tname, 'w') as file:
-    file.write(f'force: {force},\n time: {end-start},\n degree:{p},\n, backend:{backend},\n, model_type:{model_type},\n, enable_adaptive:{enable_adaptive},\n, marking_strategy:{marking_strategy},\n, refine_method:{refine_method},\n, n:{n},\n, maxit:{maxit},\n, vtkname:{vtkname}\n')
+    file.write(f'\n time: {end-start},\n degree:{p},\n, backend:{backend},\n, model_type:{model_type},\n, enable_adaptive:{enable_adaptive},\n, marking_strategy:{marking_strategy},\n, refine_method:{refine_method},\n, n:{n},\n, maxit:{maxit},\n, vtkname:{vtkname}\n')
 fig, axs = plt.subplots()
 plt.plot(disp, force, label='Force')
 plt.xlabel('Displacement Increment')

app/fracturex/fracturex/cases/phase_field/square_domian_with_fracture.py

@@ -6,7 +6,7 @@
 from fealpy.mesh import TriangleMesh, QuadrangleMesh
 
 
-from app.fracturex.fracturex.phasefield.main_solver import MainSolve
+from app.fracturex.fracturex.phasefield.main_solve import MainSolve
 from fealpy.utils import timer
 import time
 import matplotlib.pyplot as plt
@@ -120,6 +120,10 @@ def adaptive_mesh(self, mesh, d0=0.49, d1=1.01, h=0.005):
         default=False, type=bool,
         help='是否使用 GPU, 默认为 False.')
 
+parser.add_argument('--cupy', 
+        default=False, type=bool,
+        help='是否使用cupy求解.')
+
 args = parser.parse_args()
 p= args.degree
 maxit = args.maxit
@@ -133,6 +137,7 @@ def adaptive_mesh(self, mesh, d0=0.49, d1=1.01, h=0.005):
 vtkname = args.vtkname +'_' + args.mesh_type + '_'
 force_type = args.force_type
 gpu = args.gpu
+cupy = args.cupy
 
 tmr = timer()
 next(tmr)
@@ -163,8 +168,7 @@ def adaptive_mesh(self, mesh, d0=0.49, d1=1.01, h=0.005):
 fname = args.mesh_type + '_square_with_a_notch_init.vtu'
 mesh.to_vtk(fname=fname)
 
-
-ms = MainSolve(mesh=mesh, material_params=model.params, p=p, model_type=model_type)
+ms = MainSolve(mesh=mesh, material_params=model.params, model_type=model_type)
 tmr.send('init')
 
 '''
@@ -189,17 +193,19 @@ def adaptive_mesh(self, mesh, d0=0.49, d1=1.01, h=0.005):
 
 if bm.backend_name == 'pytorch':
     ms.auto_assembly_matrix()
+if cupy:
+    ms.set_cupy_solver()
 
 ms.output_timer()
 ms.save_vtkfile(fname=vtkname)
-ms.solve(maxit=maxit)
+ms.solve(p=p, maxit=maxit)
 
 tmr.send('stop')
 tmr.send(None)
 end = time.time()
 
-force = ms.Rforce
-disp = ms.force_value
+force = ms.get_residual_force()
+
 
 ftname = 'force_'+args.mesh_type + '_p' + str(p) + '_' + 'model1_disp.pt'
 
@@ -208,6 +214,12 @@ def adaptive_mesh(self, mesh, d0=0.49, d1=1.01, h=0.005):
 tname = 'params_'+args.mesh_type + '_p' + str(p) + '_' + 'model1_disp.txt'
 with open(tname, 'w') as file:
     file.write(f'\n time: {end-start},\n degree:{p},\n, backend:{backend},\n, model_type:{model_type},\n, enable_adaptive:{enable_adaptive},\n, marking_strategy:{marking_strategy},\n, refine_method:{refine_method},\n, n:{n},\n, maxit:{maxit},\n, vtkname:{vtkname}\n')
+
+if force_type == 'y':
+    disp = model.is_y_force()
+elif force_type == 'x':
+    disp = model.is_x_force()
+
 fig, axs = plt.subplots()
 plt.plot(disp, force, label='Force')
 plt.xlabel('Displacement Increment')