From d3c69d2e9ee0c0529ae56c6b55c8197a967fc779 Mon Sep 17 00:00:00 2001
From: Ernesto Oquelis <ernesto.oquelis@gmail.com>
Date: Tue, 25 Sep 2018 21:37:55 +0200
Subject: [PATCH] it seems like it works now
---
TEST.py | 7 ++---
TEST_2.py | 22 ++++++-------
fem_2d.py | 93 ++++++++++++++++++++++++++++++++-----------------------
3 files changed, 68 insertions(+), 54 deletions(-)
diff --git a/TEST.py b/TEST.py
index 46d6810..060362c 100644
--- a/TEST.py
+++ b/TEST.py
@@ -8,7 +8,7 @@ def round_trip_connect(start, end):
return [(i, i+1) for i in range(start, end)] + [(end, start)]
def refine_grid(vertices, area):
- max_area = 0.3 ########################################################################
+ max_area = 1 ########################################################################
return bool(area>max_area)
def createMesh(vertices, edges) :
@@ -32,10 +32,7 @@ def plot(mesh_points, mesh_elements) :
plt.show()
-vertices = [(0,0),
- (1,0),
- (1,1),
- (0,1)]
+vertices = [(-1,-1),(1,-1),(1,1),(-1,1)]
edges = round_trip_connect(0, len(vertices)-1)
diff --git a/TEST_2.py b/TEST_2.py
index d065e10..d41bebb 100644
--- a/TEST_2.py
+++ b/TEST_2.py
@@ -6,12 +6,12 @@ mesh_elements = np.array([[1, 2, 4],
[4, 2, 3],
[0, 1, 4]], dtype='int32')
-mesh_points = np.array([[0. , 0. ],
- [1. , 0. ],
- [1. , 1. ],
- [0. , 1. ],
- [0.5, 0.5]], dtype='float64')
-'''
+mesh_points = np.array([[-1., -1.],
+ [ 1., -1.],
+ [ 1., 1.],
+ [-1., 1.],
+ [ 0., 0.]], dtype='float64')
+
STR = Structure2DPlaneStress(mesh_elements, mesh_points,
elastic_modulus = 10e6,
poisson_ratio = 0.3,
@@ -21,13 +21,13 @@ STR = Structure2DPlaneStress(mesh_elements, mesh_points,
STRUCT = list()
STRUCT.append(STR)
-STR.nodes_list[0].constraint_2d = Constraint2D([1,1])
-STR.nodes_list[1].constraint_2d = Constraint2D([1,1])
-STR.nodes_list[2].force_2d = Force2D([1,1])
-STR.nodes_list[3].constraint_2d = Constraint2D([1,1])
+STR.node_list[0].constraint_2d = Constraint2D([1,1])
+STR.node_list[1].constraint_2d = Constraint2D([1,1])
+STR.node_list[2].force_2d = Force2D([1,1])
+STR.node_list[3].constraint_2d = Constraint2D([1,1])
STR.solve()
STR.set_displacements()
-
+'''
print(STR.structure_displacement)
'''
\ No newline at end of file
diff --git a/fem_2d.py b/fem_2d.py
index dad79ed..6ca0359 100644
--- a/fem_2d.py
+++ b/fem_2d.py
@@ -88,6 +88,7 @@ class Element2D:
self.node_2 = node_2
self.node_3 = node_3
self.node_ids = node_ids
+ self.k_e = None
const_a = np.linalg.norm(self.node_1.position - self.node_2.position)
const_b = np.linalg.norm(self.node_2.position - self.node_3.position)
@@ -112,7 +113,9 @@ class Element2D:
[0.0, b_12, 0.0, b_22, 0.0, b_32],
[b_21, b_11, b_22, b_21, b_32, b_31]])
- return self.area * thickness * np.dot(np.dot(b_mat.T, c_mat), b_mat)
+ self.k_e = self.area * thickness * np.dot(np.dot(b_mat.T, c_mat), b_mat)
+
+ return self.k_e
def print_details(self):
"""
@@ -132,9 +135,10 @@ class Structure2DPlaneStress:
"""
def __init__(self, mesh_elements, mesh_points, elastic_modulus,
poisson_ratio, thickness):
- self.node_list = list()
- self.element_list = list()
+
self.k_es = None
+ self.k_global = None
+ self.r_global = None
self.thickness = thickness
self.structure_displacement = None
self.c_mat = (elastic_modulus / (1 - poisson_ratio ** 2) *
@@ -142,56 +146,69 @@ class Structure2DPlaneStress:
[poisson_ratio, 1, 0],
[0, 0, (1-poisson_ratio)/2]]))
- for el in mesh_elements:
- ke =
-
-'''
+ # Defines list of nodes while creating them
+ self.node_list = list()
+ for node in mesh_points:
+ self.node_list.append(Node2D(node))
+ # Defines list of elements while creating them
+ self.element_list = list()
+ for element_nodes in mesh_elements:
+ self.element_list.append(Element2D(self.node_list[element_nodes[0]],
+ self.node_list[element_nodes[1]],
+ self.node_list[element_nodes[2]],
+ element_nodes))
+
def solve(self):
"""
Solves linear system of equations.
"""
+
+ num_nodes = len(self.node_list)
# Defines empty k matrix and r vector
- k_global = np.zeros(((2 * self.mesh_points.shape[0]),
- (2 * self.mesh_points.shape[0])))
- r_global = np.zeros(((2 * self.mesh_points.shape[0]), 1))
+ k_global = np.zeros(((2 * num_nodes), (2 * num_nodes)))
+ r_global = np.zeros(((2 * num_nodes), 1))
# Assembles load vector
- for node_id in range((2 * self.mesh_points.shape[0])):
- if self.nodes_list[node_id].force_2d is not None:
- r_global[node_id * 2] = self.nodes_list[node_id].force_2d.value[0]
- r_global[node_id * 2 + 1] = self.nodes_list[node_id].force_2d.value[1]
-
- # Assembles stiffness matrix
- for row in self.mesh_elements:
- k_e = elementStiffnessMatrix(self.mesh_points, row, dummyMaterial)
- for i in range(len(row)) :
- for j in range(len(row)):
- K_global[2*row[i], 2*row[j]] += G[2*i,2*j]
- K_global[(2*row[i] + 1) ,(2*row[j] + 1)] += G[(2*i + 1) , (2*j + 1)]
-
+ for node_id in range(num_nodes):
+ 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:
+ k_e = element.get_k_e_mat(self.thickness, self.c_mat)
+ element_node_ids = element.node_ids
+
+ # This works because the node numbers directly correspond to the
+ # DOF numbers through:
+ # x --> 2*DOFno.
+ # y --> 2*DOFno. + 1
+ for i in range(3):
+ for j in range(3):
+ k_global[2 * element_node_ids[i],
+ 2 * element_node_ids[j]] += k_e[2 * i, 2 * j]
+ k_global[(2 * element_node_ids[i] + 1),
+ (2 * element_node_ids[j] + 1)] += k_e[(2 * i + 1), (2 * j + 1)]
+ k_global[(2 * element_node_ids[i] + 1),
+ (2 * element_node_ids[j])] += k_e[(2 * i + 1), (2 * j)]
+ k_global[(2 * element_node_ids[i]),
+ (2 * element_node_ids[j] + 1)] += k_e[(2 * i), (2 * j + 1)]
+
+ print(k_global)
+ self.k_global = k_global
self.structure_displacement = np.linalg.solve(k_global, r_global)
def set_displacements(self):
"""
Sets displacements on the nodes according to solution of linear system.
"""
- for node_id in range(self.num_nodes):
+ 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.nodes_list[node_id].displacement = disp
-'''
-'''
- def select_displacements(self):
- """
- Sets displacements on the nodes according to solution of linear system.
- """
- for node_id in range(self.get_number_nodes()):
- dof = np.array(self.get_node(node_id).get_dof_numbers()).reshape(3,)
- displacement = np.zeros((3,))
- for i in range(3):
- if dof[i] != -1:
- displacement[i] = self.structure_displacement[int(dof[i])]
- self.get_node(node_id).set_displacement(displacement)
+ self.node_list[node_id].displacement = disp
+'''
def print_details(self):
"""
Prints details of the structure object.
--
GitLab