Examples package#
Examples index
Submodules#
Examples.Example48 module#
Examples.Example49 module#
Examples.Example50 module#
Examples.Example50-2 module#
Examples.Punto3 module#
Examples.Punto6 module#
Examples.example1 module#
Creation of 2D elements#
Coords and gdls (degrees of freedom) are given by a Numpy ndarray matrix.
In the coordinate matrix, each row represents a node of the element and each column a dimension. For example, a 2D triangular element of 3 nodes must have a 3x2 matrix.
In the gdls matrix, each row represents a variable of the node and each column a node. For example, a 2D triangular element of 3 nodes and 2 variables per node (plane stress) must have a 2x3 gdls matrix.
In this example several characteristics are tested:
The element creation and transformations are tested over 2 triangular (2D), 2 quadrilateral (2D) and 2 lineal (1D) elements.
The isBetwwen methods gives the oportunity to check if a given set of points is are inside a element.
The inverseMapping method allows to give a set of global coordinates and convert them to natural coordinates.
The jacobian graph allows to verigfy the numerical estability of the element
Coordinate trasformation and shape functions#
Serendipity (8 nodes 2D) element coordinate transformation and shape functions.#
Quadrilateral (4 nodes 2D) element coordinate transformation and shape functions.#
Triangular (6 nodes 2D) element coordinate transformation and shape functions.#
Triangular (3 nodes 2D) element coordinate transformation and shape functions.#
Line (3 nodes 1D) element coordinate transformation and shape functions.#
Line (2 nodes 2D) element coordinate transformation and shape functions.#
Point inside element (works in all dimensions)#
Test if given points are inside the element.#
Jacobian graphs#
Serendipity (8 nodes 2D) element jacobian graph.#
Quadrilateral (4 nodes 2D) element jacobian graph.#
Triangular (6 nodes 2D) element jacobian graph.#
Triangular (3 nodes 2D) element jacobian graph.#
Code#
from FEM.Elements.E2D.Serendipity import Serendipity
from FEM.Elements.E2D.Quadrilateral import Quadrilateral
from FEM.Elements.E2D.LTriangular import LTriangular
from FEM.Elements.E2D.QTriangular import QTriangular
from FEM.Elements.E1D.LinealElement import LinealElement
from FEM.Elements.E1D.QuadraticElement import QuadraticElement
elements = []
coords = [[1, 1], [3, 2], [3.5, 3], [0.5, 4], [2, 1.5],
[3.25, 2.5], [(3.5+.5)/2, 3.5], [(0.5+1)/2, 5/2]]
gdl = np.array([[1, 2, 3, 4, 5, 6, 7, 8]])
eRS = Serendipity(coords, gdl)
eRS2 = Serendipity(coords, gdl)
elements.append(eRS)
coords = [[1, 1], [3, 2], [3.5, 3], [0.5, 4]]
gdl = np.array([[1, 2, 3, 4]])
eRC = Quadrilateral(coords, gdl)
elements.append(eRC)
# [[-1.5, -1], [1, 0], [0, 5], [-0.25, -0.5], [0.5, 2.5], [-0.75, 2.0]]
coords = [[-1.5, -1], [1, 0], [0, 5]]
for i in range(len(coords)-1):
coords += [[(coords[i][0]+coords[i+1][0])/2,
(coords[i][1]+coords[i+1][1])/2]]
coords += [[(coords[i+1][0]+coords[0][0])/2,
(coords[i+1][1]+coords[0][1])/2]]
gdl = np.array([[1, 2, 3, 4, 5, 6]])
eTC = QTriangular(coords, gdl)
elements.append(eTC)
coords = [[1.4, 1.5], [3.1, 2.4], [2, 3]]
gdl = np.array([[1, 2, 3]])
eTL = LTriangular(coords, gdl)
elements.append(eTL)
coords = [[3], [4], [5]]
gdl = np.array([[3, 4, 5]])
e = QuadraticElement(coords, gdl, 3)
elements.append(e)
coords = [[3], [5]]
gdl = np.array([[3, 5]])
e = LinealElement(coords, gdl, 3)
elements.append(e)
# Coordinate transformation
for i, e in enumerate(elements):
e.draw()
plt.show()
# Point inside element
p_test = np.array([1.25, 5.32, 3.1, 3.5, 5.0])
result = e.isInside(p_test)
e.draw()
plt.gca().plot(p_test[result], [0.0] *
len(p_test[result]), 'o', c='g', label='Inside')
plt.gca().plot(p_test[np.invert(result)], [0.0] *
len(p_test[np.invert(result)]), 'o', c='r', label='Not Inside')
plt.legend()
plt.legend()
plt.show()
# Inverse Mapping
z = eTL.inverseMapping(np.array([[1.3, 2.5, 3.5], [1.5, 2.6, 8.5]]))
# z = eTL.inverseMapping(np.array([[1.3],[1.5]]))
print(z)
# print(eTL.isInside(np.array([3.5,2.5])))
# Jacobian Graphs
for i, e in enumerate(elements):
e.jacobianGraph()
plt.show()
Examples.example10 module#
Examples.example11 module#
Examples.example13 module#
Examples.example14 module#
Examples.example15 module#
Examples.example16 module#
Examples.example17 module#
Examples.example18 module#
Examples.example19 module#
Examples.example2 module#
Geometry creation using triangles. Torsion2D#
An I shape section when a unitary rotation is applied.
The input mesh is created using the Delaunay class using second order elements.
Geomery#
Input geometry created using Delaunay triangulations.#
Result#
Analysis results.#
Code#
import matplotlib.pyplot as plt
from FEM.Torsion2D import Torsion2D
from FEM.Geometry import Delaunay, Geometry2D
a = 0.3
b = 0.3
tw = 0.05
tf = 0.05
E = 200000
v = 0.27
G = E / (2 * (1 + v))
phi = 1
# Geometry perimeter
vertices = [
[0, 0],
[a, 0],
[a, tf],
[a / 2 + tw / 2, tf],
[a / 2 + tw / 2, tf + b],
[a, tf + b],
[a, 2 * tf + b],
[0, 2 * tf + b],
[0, tf + b],
[a / 2 - tw / 2, tf + b],
[a / 2 - tw / 2, tf],
[0, tf],
]
# Input string to the delaunay triangulation.
# o=2 second order elements.
# a=0.00003 Maximum element area.
params = Delaunay._strdelaunay(constrained=True, delaunay=True,
a=0.00003, o=2)
# Mesh reation.
geometria = Delaunay(vertices, params)
# geometria.exportJSON('Examples/Mesh_tests/I_test.json') # Save the mesh to a file.
# geometria = Geometry2D.importJSON('Examples/Mesh_tests/I_test.json') # Load mesh from file.
# Show the geometry.
geometria.show()
plt.show()
# Create the Torsion2D analysis.
O = Torsion2D(geometria, G, phi)
O.solve() # All finite element steps.
plt.show()
Examples.example20 module#
Examples.example21 module#
Examples.example22 module#
Examples.example23 module#
Examples.example24 module#
Examples.example25 module#
Examples.example26 module#
Examples.example27 module#
Examples.example28 module#
Examples.example29 module#
Examples.example3 module#
Geometry creation using triangles. Torsion2D#
An square shape section when a unitary rotation is applied.
The input mesh is created using the Delaunay class using second order elements.
Geomery#
Input geometry created using Delaunay triangulations.#
Result#
Analysis results.#
Code#
import matplotlib.pyplot as plt
from FEM.Torsion2D import Torsion2D
from FEM.Geometry import Geometry2D, Delaunay
a = 1
b = 1
E = 200000
v = 0.27
G = E/(2*(1+v))*0+1
phi = 1
# Use these lines to create the input geometry and file
# coords = np.array([[0, 0], [1, 0], [1, 1], [0, 1.0]])
# params = Delaunay._strdelaunay(a=0.0001, o=2)
# geo = Delaunay(coords, params)
# geo.exportJSON('Examples/Mesh_tests/Square_torsion.json')
geometria = Geometry2D.importJSON(
'Examples/Mesh_tests/Square_torsion.json')
geometria.show()
plt.show()
O = Torsion2D(geometria, G, phi, verbose=True)
O.solve()
plt.show()
Examples.example30 module#
Examples.example31 module#
Examples.example32 module#
Examples.example33 module#
Examples.example34 module#
Examples.example35 module#
Examples.example36 module#
Examples.example37 module#
Examples.example39 module#
Examples.example4 module#
Ordinary diferential equation 1D#
The current file ust the general formulation of the EDO1D Class using custom function coeficients.
Geomery#
Result#
Analysis results.#
Code#
import matplotlib.pyplot as plt
from FEM.EDO1D import EDO1D
from FEM.Geometry import Lineal
# Define functions of the diferential equation
def a(x): return (x[0]**2-2)
def c(x): return (x[0]-3)
def f(x): return (x[0]**2-2)*6+(x[0]-3)*(3*x[0]**2)
# Define border conditions. List of border conditions. In the first node, value=0.0, in the last node, value = 3.0
cbe = [[0, 0.0], [-1, 3.0]]
lenght = 1
n = 500
o = 2
geometria = Lineal(lenght, n, o)
O = EDO1D(geometria, a, c, f)
O.cbe = cbe
O.solve()
plt.show()
Examples.example40 module#
Examples.example41 module#
Examples.example42 module#
Examples.example43 module#
Examples.example44 module#
Examples.example45 module#
Examples.example46 module#
Examples.example47 module#
Examples.example5 module#
Examples.example51 module#
Examples.example52 module#
Examples.example53 module#
Examples.example54 module#
Examples.example6 module#
Test of point inside volumen (3D elements)#
This example creates a 3D Brick element and test if a given set of points are inside the element.
Coords and gdls (degrees of freedom) are given by a Numpy ndarray matrix.
In the coordinate matrix, each row represents a node of the element and each column a dimension. For example, a 2D triangular element of 3 nodes must have a 3x2 matrix.
In the gdls matrix, each row represents a variable of the node and each column a node. For example, a 2D triangular element of 3 nodes and 2 variables per node (plane stress) must have a 2x3 gdls matrix.
Code#
from FEM.Elements.E3D.Brick import Brick
import numpy as np
coords = np.array([[0, 0, 0.0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [
0, 0, 1.0], [1, 0, 1], [1.5, 1.5, 1.5], [0, 1, 1]])
coords2 = coords + 1.5
gdl = np.array([[0, 0, 0, 0, 0, 0, 0, 0]])
e = Brick(coords=coords, gdl=gdl)
e2 = Brick(coords=coords2, gdl=gdl)
domain = e.T(e.domain.T)[0]
domain2 = e2.T(e.domain.T)[0]
points = np.array(
np.array([[-1, 0, 0], [0.5, 0.5, 0.5]]).tolist()+domain.tolist())
r1 = e.isInside(points)
r2 = e.isInside(domain2)
print(r1, r2)
a = 0
Examples.example7 module#
Creation of 3D elements#
3D Elasticity beam problem using 8 node brick elements.
Code#
import numpy as np
import matplotlib.pyplot as plt
from FEM.Elasticity3D import Elasticity
from FEM.Geometry.Geometry import Geometry
E = 21000000.0
v = 0.2
h = 0.6
b = 0.3
L = 2.5
gamma = 23.54
_a = L
_b = h
_c = b
nx = 50
ny = 8
nz = 7
dx = _a/nx
dy = _b/ny
dz = _c/nz
coords = []
for i in range(nx+1):
x = i*dx
for j in range(ny+1):
y = j*dy
for k in range(nz+1):
z = k*dz
coords += [[x, y, z]]
dicc = []
def node(i, j, k): return i*(ny+1)*(nz+1)+j*(nz+1)+k
for i in range(nx):
for j in range(ny):
for k in range(nz):
node1 = node(i, j, k)
node2 = node(i+1, j, k)
node3 = node(i+1, j+1, k)
node4 = node(i, j+1, k)
node5 = node(i, j, k+1)
node6 = node(i+1, j, k+1)
node7 = node(i+1, j+1, k+1)
node8 = node(i, j+1, k+1)
dicc += [[node1, node2, node3, node4,
node5, node6, node7, node8]]
def fy(x): return -gamma
geometria = Geometry(dicc, coords, ["B1V"]*len(dicc), nvn=3, fast=True)
cbe = []
for i in range(len(coords)):
if 0.0 == coords[i][0]:
cbe += [[i*3, 0.0]]
cbe += [[i*3+1, 0.0]]
cbe += [[i*3+2, 0.0]]
geometria.cbe = cbe
O = Elasticity(geometria, E, v, gamma, fy=fy, verbose=True)
O.solve()
O.exportJSON('weird_beam.json')