在 2016年8月21日星期日 UTC+2下午10:37:37，Robert Cimrman写道：

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.