From e01e3fb21627c772c11300a3d09a5d6291c6fc3a Mon Sep 17 00:00:00 2001
From: Ernesto Oquelis <ernesto.oquelis@gmail.com>
Date: Tue, 25 Sep 2018 20:25:49 +0200
Subject: [PATCH] update
---
TEST_2.py | 9 +++--
fem_2d.py | 107 +++++++++++-------------------------------------------
2 files changed, 27 insertions(+), 89 deletions(-)
diff --git a/TEST_2.py b/TEST_2.py
index 639c6c0..d065e10 100644
--- a/TEST_2.py
+++ b/TEST_2.py
@@ -11,11 +11,11 @@ mesh_points = np.array([[0. , 0. ],
[1. , 1. ],
[0. , 1. ],
[0.5, 0.5]], dtype='float64')
-
+'''
STR = Structure2DPlaneStress(mesh_elements, mesh_points,
- elastic_modulus = 10e6,
- poisson_ratio = 0.3,
- thickness = 0.1)
+ elastic_modulus = 10e6,
+ poisson_ratio = 0.3,
+ thickness = 0.1)
STRUCT = list()
@@ -30,3 +30,4 @@ 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 574aa29..dad79ed 100644
--- a/fem_2d.py
+++ b/fem_2d.py
@@ -62,29 +62,6 @@ class Node2D:
self.displacement = np.zeros((2,))
self.constraint_2d = None
self.force_2d = None
- self.dof_numbers = np.zeros((2, 1))
-
- def enumerate_dofs(self, start):
- """
- Enumerates dofs. All restricted movements are saved as -1. The dofs are
- numbered instead.
- """
- current_dof = start
-
- if self.constraint_2d is None:
- self.dof_numbers = np.array([current_dof, current_dof + 1],
- dtype='int32').reshape(2, 1)
- current_dof = current_dof + 2
-
- else:
- for i in range(2):
- if self.constraint[i]:
- self.dof_numbers[i] = -1
- else:
- self.dof_numbers[i] = current_dof
- current_dof = current_dof + 1
-
- return current_dof
def print_details(self):
"""
@@ -111,7 +88,6 @@ class Element2D:
self.node_2 = node_2
self.node_3 = node_3
self.node_ids = node_ids
- self.dof_numbers = 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)
@@ -121,7 +97,7 @@ class Element2D:
(const_s - const_b) *
(const_s - const_c))
- def get_k_mat(self, thickness, c_mat):
+ def get_k_e_mat(self, thickness, c_mat):
"""
Computes the local stiffness matrix of the element.
"""
@@ -138,24 +114,16 @@ class Element2D:
return self.area * thickness * np.dot(np.dot(b_mat.T, c_mat), b_mat)
- def enumerate_dofs(self):
- """
- Enumerates the dofs and concatenates the results on a single array.
- """
- self.dof_numbers = np.concatenate([self.node_1.dof_numbers,
- self.node_2.dof_numbers,
- self.node_3.dof_numbers])
-
def print_details(self):
"""
Pritns details of the element object.
"""
print('\nDetails of node 1:')
- self.node_1().print_node()
+ self.node_1().print_details()
print('Details of node 2:')
- self.node_2().print_node()
+ self.node_2().print_details()
print('Details of node 3:')
- self.node_3().print_node()
+ self.node_3().print_details()
class Structure2DPlaneStress:
@@ -164,26 +132,9 @@ class Structure2DPlaneStress:
"""
def __init__(self, mesh_elements, mesh_points, elastic_modulus,
poisson_ratio, thickness):
- self.nodes_list = list()
- self.elements_list = list()
- self.number_equations = int(0)
-
- nodes_array = mesh_points
- self.num_nodes = nodes_array.shape[0]
- self.num_eq = 2 * self.num_nodes
- for node_id in range(self.num_nodes):
- self.nodes_list.append(Node2D(nodes_array[node_id, :]))
-
- elements_array = mesh_elements
- self.num_elements = mesh_elements.shape[0]
- for element_id in range(self.num_elements):
- element_nodes = elements_array[element_id, :]
- self.elements_list.append(Element2D(self.nodes_list[element_nodes[0]],
- self.nodes_list[element_nodes[1]],
- self.nodes_list[element_nodes[2]],
- element_nodes))
- self.elements_list[-1].node_ids = element_nodes
-
+ self.node_list = list()
+ self.element_list = list()
+ self.k_es = None
self.thickness = thickness
self.structure_displacement = None
self.c_mat = (elastic_modulus / (1 - poisson_ratio ** 2) *
@@ -191,46 +142,32 @@ class Structure2DPlaneStress:
[poisson_ratio, 1, 0],
[0, 0, (1-poisson_ratio)/2]]))
- def enumerate_dofs(self):
- """
- Enumerates the dofs of all nodes on the structure.
- """
- current_dof = int(0)
- for i in range(self.num_nodes):
- current_dof = self.nodes_list[i].enumerate_dofs(current_dof)
- for j in range(self.num_elements):
- self.get_element(j).enumerate_dofs()
- return current_dof
-
+ for el in mesh_elements:
+ ke =
+
+'''
def solve(self):
"""
Solves linear system of equations.
"""
# Defines empty k matrix and r vector
- k_global = np.zeros((self.num_eq, self.num_eq))
- r_global = np.zeros((self.num_eq, 1))
+ 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))
# Assembles load vector
- for node_id in range(self.num_nodes):
+ 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 element_id in range(self.num_elements):
- k_g_ix = list()
- node_ids = self.elements_list[element_id].node_ids
- k_e = self.elements_list[element_id].get_k_mat(self.thickness,
- self.c_mat)
- for i_x in node_ids:
- k_g_ix.append(i_x * 2)
- k_g_ix.append(i_x * 2 + 1)
-
- for i in range(6):
- for j in range(6):
- k_global[k_g_ix[i], k_g_ix[j]] = ((k_global[k_g_ix[i],
- k_g_ix[j]]) +
- k_e[i, j])
+ 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)]
self.structure_displacement = np.linalg.solve(k_global, r_global)
@@ -242,7 +179,7 @@ class Structure2DPlaneStress:
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.
--
GitLab