To debug a GetFEM code, I am doing “hand-made” computations on a very simple example. Before moving to the contact algorithm, I need to make sure that things work well for the following toy problem without contact.
The unit square (0,1)\times(0,1) is occupied by a linear elastic isotropic material. I apply a vertical surface load on the upper side y=1, say, \mathbf g=(0,g) (let’s forget about units), and the left side x=0 is clamped.
I discretize the problem with only two triangular elements (\mathbb P_1 elements), meaning that I have the following global node numbering:
(0,0) : 1
(1,0) : 2
(1,1) : 3
(0,1) : 4
so that the mesh is made up by the two right triangles 1-2-4 and 2-3-4.
Local indexing of the nodes is always counterclockwise, of course (starting with the right angle).
I want to compute by hand the value of the displacement field at nodes 2 and 3, which are the only free nodes (and actually the surface load only “works” on node 3).
The local stiffness matrix (independent of the element) I obtain is the following 6x6 matrix:
K_e = {1 \over 2} \begin{pmatrix} \lambda + 3\mu & \lambda + \mu & -\lambda - 2\mu & -\mu & -\mu & -\lambda \\ \lambda + \mu & \lambda + 3\mu & -\lambda & -\mu & -\mu & -\lambda - 2\mu \\ -\lambda - 2\mu & -\lambda & \lambda + 2\mu & 0 & 0 & \lambda \\ -\mu & -\mu & 0 & \mu & \mu & 0 \\ -\mu & -\mu & 0 & \mu & \mu & 0 \\ -\lambda & -\lambda - 2\mu & \lambda & 0 & 0 & \lambda + 2\mu \end{pmatrix} ,
where \lambda and \mu are the usual Lamé moduli.
Let’s denote the displacement at a node i as (u_i, v_i). Assembling the global stiffness matrix, the global right-hand side, and taking account of the homogeneous Dirichlet boundary conditions on nodes 1 and 4 (so that u_1=v_1=u_4=v_4=0), I obtain the following reduced linear system (4 scalar unknowns instead of 8):
K_{glob} U = F_{glob},
where the unknown is U = (u_2, v_2, u_3, v_3)^\top, F_{glob} = (0,0,0,g/2)^\top, and the following (reduced) global stiffness matrix:
K_{glob} = {1\over 2} \begin{pmatrix} \lambda + 3\mu & 0 & -\mu & -\mu \\ 0 & \lambda + 3\mu & -\lambda & -\lambda - 2\mu \\ -\mu & -\lambda & \lambda + 3\mu & \lambda + \mu \\ -\mu & -\lambda - 2\mu & \lambda + \mu & \lambda + 3\mu \end{pmatrix} .
I am assuming plane strains, and for the values E=20 and \nu = 0.2 of Young’s modulus and Poisson’s ratio, and g=-5 for the vertical (downwards) surface load, I get the following solution:
U = (u_2,v_2,u_3,v_3)^\top = (-0.0983193, -0.337815, 0.138655, -0.49916)^\top.
Nevertheless, using a GetFEM test script (the mesh and global node numbering are slightly different in this case), I get the following values, after properly rearranging values for comparison:
U = (u_2,v_2,u_3,v_3)^\top = (-0.0605042, -0.337815, 0.07563025, -0.39831933)^\top.
What are these differences due to? Am I missing something in the assembly procedure?
I am certain that the values obtained with GetFEM are the correct ones, but I cannot figure out what is wrong in the procedure I carried out.