Merge pull request #1390 from tiantian0347/master

FractureX 程序测试

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

@@ -68,8 +68,8 @@ def is_dirchlet_boundary(self, p):
 
 
 parser.add_argument('--mesh_type',
-        default='tri', type=str,
-        help='网格类型, 默认为 tri.')
+        default='tet', type=str,
+        help='网格类型, 默认为 tet.')
 
 parser.add_argument('--enable_adaptive',
         default=False, type=bool,
@@ -167,7 +167,7 @@ def is_dirchlet_boundary(self, p):
 end = time.time()
 
 force = ms.get_residual_force()
-disp = model.is_force()
+disp = model.is_z_force()
 
 ftname = 'force_'+args.mesh_type + '_p' + str(p) + '_' + 'model3d_disp.pt'
 
@@ -177,6 +177,7 @@ def is_dirchlet_boundary(self, p):
 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')
 fig, axs = plt.subplots()
+disp = model.is_z_force()
 plt.plot(disp, force, label='Force')
 plt.xlabel('Displacement Increment')
 plt.ylabel('Residual Force')

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

@@ -19,7 +19,7 @@ def __init__(self):
         E = 210
         nu = 0.3
         Gc = 2.7e-3
-        l0 = 0.02
+        l0 = 0.015
         self.params = {'E': E, 'nu': nu, 'Gc': Gc, 'l0': l0}
 
     def is_y_force(self):
@@ -211,6 +211,7 @@ def adaptive_mesh(self, mesh, d0=0.49, d1=1.01, h=0.005):
 
 torch.save(force, ftname)
 #np.savetxt('force'+tname, bm.to_numpy(force))
+
 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')
@@ -219,6 +220,8 @@ def adaptive_mesh(self, mesh, d0=0.49, d1=1.01, h=0.005):
     disp = model.is_y_force()
 elif force_type == 'x':
     disp = model.is_x_force()
+else:
+    raise ValueError('Invalid force type.')
 
 fig, axs = plt.subplots()
 plt.plot(disp, force, label='Force')

app/fracturex/fracturex/phasefield/main_solve.py

@@ -75,7 +75,7 @@ def __init__(self, mesh, material_params: Dict,
         self.tmr = timer()
         next(self.tmr)
 
-    def initialization_settings(self, p: int = 1, q: int = None, ):
+    def initialize_settings(self, p: int = 1, q: int = None, ):
         """
         Initialize the settings for the problem.
         """
@@ -121,11 +121,11 @@ def solve(self, method: str = 'lfem', p: int = 1, q: int = None, maxit: int = 50
             VTK output file name, by default None.
         """
         self._method = method
-        self.initialization_settings(p=p, q=q)
+        self.initialize_settings(p=p, q=q)
         self._initialize_force_boundary()
         self._Rforce = bm.zeros_like(self._force_value)
 
-#        for i in range(1):
+        #for i in range(2):
         for i in range(len(self._force_value)-1):
             print('i', i)
             self._currt_force_value = self._force_value[i+1]
@@ -243,6 +243,7 @@ def solve_displacement(self) -> float:
         du = self.solver(A, R, atol=1e-20)
         uh += du[:]
         self.uh = uh
+        print('uh', uh)
 
         self.pfcm.update_disp(uh)
         tmr.send('disp_solver')
@@ -258,32 +259,42 @@ def solve_phase_field(self) -> float:
             The norm of the residual.
         """
         Gc, l0, d = self.Gc, self.l0, self.d
+
+        @barycentric
+        def diff_coef(bc, index):
+            gg_hd, c_d = self.CSDFunc.grad_grad_density_function(d(bc))
+            return Gc * l0 * 2 / c_d
+
+        @barycentric
+        def mass_coef1(bc, index):
+            gg_hd, c_d = self.CSDFunc.grad_grad_density_function(d(bc))
+            return gg_hd * Gc / (l0 * c_d)
+
 
         @barycentric
-        def coef(bc, index):
-            gg_gd = self.EDFunc.grad_grad_degradation_function(d)
+        def mass_coef2(bc, index):
+            gg_gd = self.EDFunc.grad_grad_degradation_function(d(bc))
             return gg_gd * self.pfcm.maximum_historical_field(bc)
 
         @barycentric
         def source_coef(bc, index):
-            gg_gd = self.EDFunc.grad_degradation_function_constant_coef()
-            return -1 * gg_gd * self.pfcm.maximum_historical_field(bc)
+            gc_gd = self.EDFunc.grad_degradation_function_constant_coef()
+            return -1 * gc_gd * self.pfcm.maximum_historical_field(bc)
 
         tmr = self.tmr
         tmr.send('phase_start')
 
-        gg_hd, c_d = self.CSDFunc.grad_grad_density_function(d)
-
         dbform = BilinearForm(self.space)
-        dbform.add_integrator(ScalarDiffusionIntegrator(coef=Gc * l0 * 2 / c_d, q=self.q))
-        dbform.add_integrator(ScalarMassIntegrator(coef=gg_hd * Gc / (l0 * c_d), q=self.q))
-        dbform.add_integrator(ScalarMassIntegrator(coef=source_coef, q=self.q))
+        dbform.add_integrator(ScalarDiffusionIntegrator(coef=diff_coef, q=self.q))
+        dbform.add_integrator(ScalarMassIntegrator(coef=mass_coef1, q=self.q))
+        dbform.add_integrator(ScalarMassIntegrator(coef=mass_coef2, q=self.q))
         A = dbform.assembly()
         tmr.send('phase_matrix_assemble')
 
         dlform = LinearForm(self.space)
-        dlform.add_integrator(ScalarSourceIntegrator(source=coef, q=self.q))
+        dlform.add_integrator(ScalarSourceIntegrator(source=source_coef, q=self.q))
         R = dlform.assembly()
+        print('R', R)
         R -= A @ d[:] 
 
 
@@ -293,7 +304,9 @@ def source_coef(bc, index):
         tmr.send('phase_apply_bc')
 
         dd = self.solver(A, R, atol=1e-20)
+        print('R', R)
         d += dd[:]
+        print('d', d)
 
 
         self.d = d
@@ -409,11 +422,6 @@ def mass_kernel_func2(u):
         def mass_grad_kernel_func2(u):
             return self.EDFunc.grad_grad_degradation_function(u)
 
-        @barycentric
-        def source_coef(bc, index):
-            gg_gd = self.EDFunc.grad_degradation_function_constant_coef()
-            return -1 * gg_gd * self.pfcm.maximum_historical_field(bc)
-
         diffusion_coef.kernel_func = diffusion_kernel_func
         mass_coef1.kernel_func = mass_kernel_func1
         mass_coef2.kernel_func = mass_kernel_func2
@@ -667,7 +675,7 @@ def set_energy_degradation(self, degradation_type='quadratic', EDfunc=None, **kw
             Additional parameters passed to the energy degradation function.
         """
         if EDfunc is not None:
-            self.EDFunc = EDfunc(degradation_type=degradation_type, **kwargs)
+            self.EDFunc = EDfunc(degradation_type='user_defined', **kwargs)
         else:
             self.EDFunc = EDFunc(degradation_type=degradation_type)
 
@@ -685,7 +693,7 @@ def set_crack_surface_density(self, density_type='AT2', CSDfunc=None, **kwargs):
             Additional parameters passed to the crack surface density function.
         """
         if CSDfunc is not None:
-            self.CSDFunc = CSDfunc(density_type=density_type, **kwargs)
+            self.CSDFunc = CSDfunc(density_type=='user_defined', **kwargs)
         else:
             self.CSDFunc = CSDFunc(density_type=density_type)
 

app/fracturex/fracturex/phasefield/phase_fracture_material.py

@@ -334,6 +334,16 @@ def heaviside(self, x):
         val[x < -1e-13] = 0
         return val
 
+    def linear_strain_value(self, bc):
+        """
+        Compute the linear strain tensor.
+        """
+        bc = bm.array([[1/3, 1/3, 1/3]], dtype=bm.float64)
+        uh = self.uh
+        guh = uh.grad_value(bc)
+        strain = 0.5 * (guh + bm.swapaxes(guh, -2, -1))
+        return strain
+
     @ barycentric
     def maximum_historical_field(self, bc):