Hi Ronghai,
On 08/19/2016 07:23 PM, Ronghai Wu wrote:
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)?
You have 201 * 201 * 2 = 80802 degrees of freedom (DOFs) - two displacements in each vertex, and some DOFs are probably fixed by fix_u_lef_corn boundary condition (what are your regions?), so the shape is (79999, 79999). Concerning the number of non-zeros, it corresponds to the mesh connectivity (how many nodes are connected to a node), and indicates how much memory the matrix storage would take.
(2) I do not understand why "updating materials..." did twice.
Because you called:
self.pb.update_materials()
Problem.update_materials() is called automatically in Problem.solve().
(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?
This depends on the linear solver used. You use ScipyDirect which is either superLU (rather slow), or umfpack (faster), provided scikit-umfpack and the umfpack libraries are installed.
You can try another solver (e.g. pyamg, or an iterative solver from petsc) to speed things up.
(4) In which source code file the following information are printed?
You can use 'git grep' to find the patterns if you have the git version of sfepy.
r.