Skip to content
Snippets Groups Projects
Commit 6107ee67 authored by Ernesto Oquelis's avatar Ernesto Oquelis
Browse files

last update of today, life is good

parent d4613bea
No related branches found
No related tags found
1 merge request!4last update of today, life is good
......@@ -8,7 +8,6 @@ Written on: 15.Sep.18
Last edited on: 25.Sep.18
"""
import numpy as np
......@@ -164,7 +163,9 @@ class Structure2DPlaneStress:
element_nodes))
def reduce(self, k_global, r_global):
import numpy as np
"""
Missing docstring...
"""
to_delete = list()
to_keep = list()
......@@ -173,9 +174,9 @@ class Structure2DPlaneStress:
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):
for i_x, con in enumerate(constr.value):
if con:
to_delete.append(2 * node_id + ix)
to_delete.append(2 * node_id + i_x)
to_keep = np.arange(2 * num_nodes)
to_keep = list(to_keep)
for element in to_delete:
......@@ -187,11 +188,10 @@ class Structure2DPlaneStress:
self.kept = to_keep
for ix, i in enumerate(to_keep):
r_new[ix] = r_global[i]
for jx, j in enumerate(to_keep):
k_new[ix,jx] = k_new[ix,jx] + k_global[i,j]
for i_x, i in enumerate(to_keep):
r_new[i_x] = r_global[i]
for j_x, j in enumerate(to_keep):
k_new[i_x, j_x] = k_new[i_x, j_x] + k_global[i, j]
return k_new, r_new
......@@ -231,7 +231,6 @@ class Structure2DPlaneStress:
(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)]
k_global, r_global = self.reduce(k_global, r_global)
self.k_global = k_global
......@@ -242,10 +241,10 @@ class Structure2DPlaneStress:
"""
Sets displacements on the nodes according to solution of linear system.
"""
for ix, for_node_id in enumerate(self.kept):
for i_x, 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]
disp = self.structure_displacement[i_x]
if even:
self.node_list[node_id].displacement[1] = disp
else:
......@@ -291,8 +290,8 @@ class Structure2DPlaneStress:
for element_id in range(len(self.element_list)):
area = self.element_list[element_id].area
ids = self.element_list[element_id].node_ids
print('{:<10}{:>10.2E}{:>10}{:>10}{:>10}'.format(element_id, area,
ids[0], ids[1], ids[2]))
print('{:<10}{:>10.2E}{:>10}{:>10}{:>10}'.format(element_id, area,
ids[0], ids[1], ids[2]))
def print_results(self):
"""
......@@ -304,7 +303,7 @@ class Structure2DPlaneStress:
for node_id in range(len(self.node_list)):
displacement = self.node_list[node_id].displacement
print('{:<10}{:>15.6E}{:>15.6E}'.format(node_id, displacement[0],
displacement[1]))
displacement[1]))
'''
class Visualizer:
......@@ -460,29 +459,41 @@ class Visualizer:
self.viewer.add_object(poly_set)
'''
def round_trip_connect(start, end):
"""
Missing docstring...
"""
return [(i, i+1) for i in range(start, end)] + [(end, start)]
def refine_grid(vertices, area):
"""
Missing docstring...
"""
max_area = 1 ########################################################################
return bool(area>max_area)
return bool(area > max_area)
def createMesh(vertices, edges) :
def create_mesh(vertices, edges):
"""
Missing docstring...
"""
from meshpy import triangle
gridInfo = triangle.MeshInfo()
gridInfo.set_points(vertices)
gridInfo.set_facets(edges)
mesh = triangle.build(gridInfo, refinement_func=refine_grid)
mesh_points = np.array(mesh.points)
grid_info = triangle.MeshInfo()
grid_info.set_points(vertices)
grid_info.set_facets(edges)
mesh = triangle.build(grid_info, refinement_func=refine_grid)
mesh_points = np.array(mesh.points)
mesh_elements = np.array(mesh.elements)
return mesh_points, mesh_elements
def plot(mesh_points, mesh_elements):
def plot_2d(mesh_points, mesh_elements):
"""
Missing docstring...
"""
import matplotlib.pyplot as plt
plt.triplot(mesh_points[:, 0], mesh_points[:, 1], mesh_elements)
plt.gca().set_aspect('equal')
for x,y in enumerate(mesh_points) :
plt.text(y[0] - 0.01, y[1] + 0.01, x, ha='right')
for x_plot, y_plot in enumerate(mesh_points):
plt.text(y_plot[0] - 0.01, y_plot[1] + 0.01, x_plot, ha='right')
plt.show()
"""
This is a missing docstring...
"""
from fem_2d import createMesh, round_trip_connect
from fem_2d import create_mesh, round_trip_connect
def simple_square():
vertices = [(-1,-1),(1,-1),(1,1),(-1,1)]
edges = round_trip_connect(0, len(vertices)-1)
return createMesh(vertices,edges)
"""
Missing docstring.
"""
vertices = [(-1, -1), (1, -1), (1, 1), (-1, 1)]
edges = round_trip_connect(0, len(vertices) - 1)
return create_mesh(vertices, edges)
"""
Missing docstring...
"""
import numpy as np
import matplotlib.pyplot as plt
from fem_2d import Structure2DPlaneStress, Constraint2D, Force2D, plot
from fem_2d import Structure2DPlaneStress, Constraint2D, Force2D
from models_2d import simple_square
# This test the simple square model
SQ_mesh_points, SQ_mesh_elements = simple_square()
SQ_MESH_POINTS, SQ_MESH_ELEMENTS = simple_square()
STR = Structure2DPlaneStress(SQ_mesh_elements, SQ_mesh_points,
elastic_modulus = 10e6,
poisson_ratio = 0.3,
thickness = 0.1)
STR = Structure2DPlaneStress(SQ_MESH_ELEMENTS, SQ_MESH_POINTS,
elastic_modulus=10e6,
poisson_ratio=0.3,
thickness=0.1)
STRUCT = list()
STRUCT.append(STR)
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.node_list[4].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.node_list[4].constraint_2d = Constraint2D([1, 1])
STR.solve()
STR.set_displacements()
STR.print_results()
\ No newline at end of file
STR.print_results()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment