From f12ae5adfcec898337c7f3eacd160437b2e0339d Mon Sep 17 00:00:00 2001
From: Ernesto Oquelis <ernesto.oquelis@gmail.com>
Date: Tue, 25 Sep 2018 22:37:01 +0200
Subject: [PATCH] almost working!!!
---
fem_2d.py | 37 +++++++++++++++++++++++++++++++++++--
1 file changed, 35 insertions(+), 2 deletions(-)
diff --git a/fem_2d.py b/fem_2d.py
index 6ca0359..4b8da5c 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):
--
GitLab