Hi Robert:

I’ve been able to get the stiffness matrix under your help, but I still get several puzzles about the linear system components [1] (8):

1. Is the linear system in the matrix formula: Kx=r ? I exported both K and r to matlab to used the right division to get x: x = r \ K. But this gave me a different result from variables = pb.solve() and x = variables.vec. Analysis showed that pb.solve() gave the appropriate result, what’s the difference between the two approaches?

2. In my project, the measurements (knowns) are Φ on surface (∂Ω)  and the unknown is the source distribution Q in attached ProblemDescription.docx. So I need to separate the vectors Φ and Q (nodal values of Φ(x) and Q(x)) explicitly in the matrix form, which means I need to get the matrix F: (FQ is the discretized form of in (2), where Q(x) is the source distribution, s(x) is the test function). How can I get this matrix?

3. I found that the solver in SfePy was much faster than the right division in matlab, so I’m curious that which function does SfePy adopts to solve the linear system. Besides, my project requires the inversion of a sub-matrix of the stiffness matrix K, does SfePy offer such methods or which function do you recommend to invert such a banded sparse matrix?