Hi Robert,
I am using sfepy for static linear elasticity problem in 2D. The mesh size is 200*200 (mesh file is attached) and equations are as following:
def D_function(self, ts, coors, mode, term=None, **kwargs): if mode != 'qp': return _, weights = term.integral.get_qp('2_4') self.n_qp = weights.shape[0] self.val_D = np.repeat(self.stiffness_tensor, self.n_qp, axis=0) return {'function' : self.val_D}
def prestress_function(self,ts, coors, mode, term=None, **kwargs): if mode != 'qp': return _, weights = term.integral.get_qp('2_4') self.n_qp = weights.shape[0] self.val_pre = np.repeat(self.prestress_value, self.n_qp, axis=0) return {'function' : self.val_pre}
self.m = Material('m', function=self.D_function) self.prestress_value = np.copy(self.stress_0) (self.prestress_value).shape = (-1, 3, 1)
self.prestress = Material('prestress', function=self.prestress_function) integral = Integral('i', order=2) self.t1 = Term.new('dw_lin_elastic(m.function, v, u )', integral, self.omega, m=self.m, v=self.v, u=self.u) self.t2 = Term.new('dw_lin_prestress(prestress.function, v )',integral, self.omega, prestress=self.prestress, v=self.v) self.t3 = Term.new('dw_surface_ltr(traction_lef.val, v)', integral, self.lef, traction_lef=self.traction_lef, v=self.v) self.t4 = Term.new('dw_surface_ltr(traction_rig.val, v)', integral, self.rig, traction_rig=self.traction_rig, v=self.v) self.t5 = Term.new('dw_surface_ltr(traction_top.val, v)', integral, self.top, traction_top=self.traction_top, v=self.v) self.t6 = Term.new('dw_surface_ltr(traction_bot.val, v)', integral, self.bot, traction_bot=self.traction_bot, v=self.v) self.fem_eq = Equation('balance', self.t1 - self.t2 - self.t3 - self.t4 - self.t5 - self.t6) self.eqs = Equations([self.fem_eq]) self.ls = ScipyDirect({}) self.nls_status = IndexedStruct() self.nls = Newton({'i_max' : 1,'eps_a' : 1e-8,'problem' : 'nonlinear'}, lin_solver=self.ls, status=self.nls_status) self.pb = Problem('elasticity', equations=self.eqs, nls=self.nls, ls=self.ls)
self.pb.time_update(ebcs=Conditions([fix_u_lef_corn]))
self.pb.update_materials() self.vec = self.pb.solve()
When I run the simulation, the information below are printed. I have several questions about it. I tried to find answers from source codes but failed. Could you please give me some help, thanks. (1) I do understand why "matrix shape: (79999, 79999)" and "matrix structural nonzeros: 1439965". what exactly this "matrix" is? what is its relation with mesh size (200*200)? (2) I do not understand why "updating materials..." did twice. (3) "solve: 6.18 [s]". I guess it should not take such a long time for a mesh size (200*200) problem. Is it because I generate mesh file in the wrong way or because of other reasons? (4) In which source code file the following information are printed?
sfepy: updating variables... sfepy: ...done sfepy: setting up dof connectivities... sfepy: ...done in 0.00 s sfepy: matrix shape: (79999, 79999) sfepy: assembling matrix graph... sfepy: ...done in 0.05 s sfepy: matrix structural nonzeros: 1439965 (2.25e-04% fill) sfepy: updating materials... sfepy: m sfepy: traction_bot sfepy: traction_rig sfepy: prestress sfepy: traction_top sfepy: traction_lef sfepy: ...done in 0.05 s sfepy: updating materials... sfepy: m sfepy: traction_bot sfepy: traction_rig sfepy: prestress sfepy: traction_top sfepy: traction_lef sfepy: ...done in 0.05 s sfepy: nls: iter: 0, residual: 4.741758e+01 (rel: 1.000000e+00) sfepy: rezidual: 0.04 [s] sfepy: solve: 6.18 [s] sfepy: matrix: 0.10 [s] sfepy: nls: iter: 1, residual: 5.936719e-14 (rel: 1.252008e-15) sfepy: equation "tmp": sfepy: ev_cauchy_strain.2.Omega(u) sfepy: updating materials... sfepy: ...done in 0.00 s
Regards Ronghai