diff --git a/fem_2d.py b/fem_2d.py
index 6ca0359c6cf51e10a68839d0c3288d81ecf3bc2a..4b8da5ca9ce73087ced23f9248f4a93ff17fcb12 100644
--- a/fem_2d.py
+++ b/fem_2d.py
@@ -158,6 +158,38 @@ class Structure2DPlaneStress:
                                                self.node_list[element_nodes[2]],
                                                element_nodes))
 
+    def reduce(self, k_global, r_global):
+        import numpy as np
+
+        to_delete = list()
+        to_keep = list()
+
+        num_nodes = len(self.node_list)
+        for node_id in range(num_nodes):
+            constr = self.node_list[node_id].constraint_2d
+            if constr is not None:
+                for ix, con in enumerate(constr.value):
+                    if con:
+                        to_delete.append(2 * node_id + ix)
+        to_keep = np.arange(2 * num_nodes)
+        to_keep = list(to_keep)
+        for element in to_delete:
+            to_keep.remove(element)
+        new_size = len(to_keep)
+
+        k_new = np.zeros((new_size, new_size))
+        r_new = np.zeros((new_size, 1))
+        print(r_new)
+
+        for ix, i in enumerate(to_keep):
+            r_new[ix,1] = r_global[i,1]
+            for jx, j in enumerate(to_keep):
+                k_new[ix,jx] = r_global[i,j]
+                
+
+        return k_new, r_new
+
+
     def solve(self):
         """
         Solves linear system of equations.
@@ -173,7 +205,6 @@ class Structure2DPlaneStress:
             if self.node_list[node_id].force_2d is not None:
                 r_global[node_id * 2] = self.node_list[node_id].force_2d.value[0]
                 r_global[node_id * 2 + 1] = self.node_list[node_id].force_2d.value[1]
-        self.r_global = r_global
 
         # Assembles stiffness matrixes
         for element in self.element_list:
@@ -195,8 +226,10 @@ class Structure2DPlaneStress:
                     k_global[(2 * element_node_ids[i]),
                              (2 * element_node_ids[j] + 1)] += k_e[(2 * i), (2 * j + 1)]
         
-        print(k_global)
+
+        k_global, r_global = self.reduce(k_global, r_global)
         self.k_global = k_global
+        self.r_global = r_global
         self.structure_displacement = np.linalg.solve(k_global, r_global)
 
     def set_displacements(self):