mmmm maybe it works ok now

fem_2d.py

@@ -139,6 +139,7 @@ class Structure2DPlaneStress:
         self.k_es = None
         self.k_global = None
         self.r_global = None
+        self.kept = None
         self.thickness = thickness
         self.structure_displacement = None
         self.c_mat = (elastic_modulus / (1 - poisson_ratio ** 2) *
@@ -179,12 +180,13 @@ class Structure2DPlaneStress:
 
         k_new = np.zeros((new_size, new_size))
         r_new = np.zeros((new_size, 1))
-        print(r_new)
+
+        self.kept = to_keep
 
         for ix, i in enumerate(to_keep):
-            r_new[ix,1] = r_global[i,1]
+            r_new[ix] = r_global[i]
             for jx, j in enumerate(to_keep):
-                k_new[ix,jx] = r_global[i,j]
+                k_new[ix,jx] = k_new[ix,jx] + k_global[i,j]
                 
 
         return k_new, r_new
@@ -236,10 +238,14 @@ class Structure2DPlaneStress:
         """
         Sets displacements on the nodes according to solution of linear system.
         """
-        num_nodes = len(self.node_list)
-        for node_id in range(num_nodes):
-            disp = self.structure_displacement[(2 * node_id):(2 * node_id + 1)]
-            self.node_list[node_id].displacement = disp
+        for ix, for_node_id in enumerate(self.kept):
+            node_id = int(for_node_id/2)
+            even = bool(for_node_id%2)
+            disp = self.structure_displacement[ix]
+            if even:
+                self.node_list[node_id].displacement[1] = disp
+            else:
+                self.node_list[node_id].displacement[0] = disp
 
 '''
     def print_details(self):