Bernoulli Console with Random Geometry

seed = 0  # integer or None
summary = {}  # we will this up during execution

import numpy as np

rs = np.random.RandomState()
if isinstance(seed, int):
    rs.seed(seed)


def brownian1d(T=1, N=100, mu=0.1, sigma=0.01, S0=20):
    dt = float(T)/N
    t = np.linspace(0, T, N)
    W = rs.standard_normal(size=N)
    W = np.cumsum(W)*np.sqrt(dt)  # standard brownian motion
    X = (mu-0.5*sigma**2)*t + sigma*W
    S = S0*np.exp(X)  # geometric brownian motion
    return S


params = dict(T=0.1, N=100, mu=0.5, sigma=0.1)

x = brownian1d(**params)
y = brownian1d(**params)
z = brownian1d(**params)

from polymesh import PolyData, CartesianFrame, LineData
import numpy as np

frame = CartesianFrame(dim=3)

coords = np.stack([x, y, z], axis=1)
topo = np.zeros((coords.shape[0]-1, 2)).astype(int)
topo[:, 0] = np.arange(topo.shape[0])
topo[:, 1] = topo[:, 0] + 1

mesh = PolyData(coords=coords, frame=frame)
mesh['lines'] = LineData(topo=topo, frame=frame)

mesh.plot(notebook=False)

Index and coordinate of the last node. Later, when the 2-noded mesh is transformed into a mesh of 3-noded elements, this information is going to be necessary to identify the free end.

i_last = len(coords) - 1
i_last, coords[i_last]

(99, array([21.41614825, 21.56705787, 20.62494753]))

The parameters of the following cell are up for manipulation. They define a simple rectangular cross section and an isotropic Hooke model governed by Young’s modulus and Poisson’s ratio.

# INPUT DATA
w, h = 10., 10.  # width and height of the rectangular cross section
F = -100.  # value of the vertical load at the free end
E = 210000.0  # Young's modulus
nu = 0.3  # Poisson's ratio

The model stiffness matrix is calculated automatically here:

import numpy as np
from numpy import pi as PI

# cross section
A = w * h  # area
Iy = w * h**3 / 12  # second moment of inertia around the y axis
Iz = w * h**3 / 12  # second moment of inertia around the z axis
Ix = Iy + Iz  # torsional inertia

# model stiffness matrix
G = E / (2 * (1 + nu))
Hooke = np.array([
    [E*A, 0, 0, 0],
    [0, G*Ix, 0, 0],
    [0, 0, E*Iy, 0],
    [0, 0, 0, E*Iz]
])

The B2 element

Run a linear analysis and do some postprocessing using 2-Noded Bernoulli members.

from sigmaepsilon.solid import Structure, LineMesh, PointData
from neumann.linalg import linspace, Vector
from polymesh.space import StandardFrame, PointCloud, frames_of_lines
from sigmaepsilon.solid.fem.cells import B2 as Beam

# space
GlobalFrame = StandardFrame(dim=3)

# support at the first, load at the last node
loads = np.zeros((coords.shape[0], 6))
fixity = np.zeros((coords.shape[0], 6)).astype(bool)
global_load_vector = Vector([0., 0, F], frame=GlobalFrame).show()
loads[-1, :3] = global_load_vector
fixity[0, :] = True

# pointdata
pd = PointData(coords=coords, frame=GlobalFrame,
               loads=loads, fixity=fixity)

# celldata
frames = frames_of_lines(coords, topo)
cd = Beam(topo=topo, frames=frames)

# set up mesh and structure
mesh = LineMesh(pd, cd, model=Hooke, frame=GlobalFrame)
structure = Structure(mesh=mesh)
structure.linsolve()

# postproc
# 1) displace the mesh
dofsol = structure.nodal_dof_solution(store='dofsol')
forces = structure.nodal_forces()
reactions = structure.reaction_forces()

Vertical displacement of the free end, external load at the free end and Vertical reaction force at the support:

dofsol[-1, 2], forces[-1, 2], reactions[0, 2]

(-1.985415419091222e-05, -100.0, 99.99999993859637)

Work of the external forces at the free end:

dofsol[-1] @ forces[-1]

0.0019854154190912217

internal_forces = structure.internal_forces()
internal_forces[0, 0, 2]

-99.81153149633388

The B3 element

We transform the mesh using the appropriate transformation function from PolyMesh.

from polymesh.topo.tr import L2_to_L3
coords, topo = L2_to_L3(coords, topo)

Rebuild and calculate the mesh using 3-Node elements.

from sigmaepsilon.solid.fem.cells import B3 as Beam

# support at the first, load at the last node
loads = np.zeros((coords.shape[0], 6))
fixity = np.zeros((coords.shape[0], 6)).astype(bool)
global_load_vector = Vector([0., 0, F], frame=GlobalFrame).show()
loads[i_last, :3] = global_load_vector
fixity[0, :] = True

# pointdata
pd = PointData(coords=coords, frame=GlobalFrame,
               loads=loads, fixity=fixity)

# celldata
frames = frames_of_lines(coords, topo)
cd = Beam(topo=topo, frames=frames)

# set up mesh and structure
mesh = LineMesh(pd, cd, model=Hooke, frame=GlobalFrame)
structure = Structure(mesh=mesh)
structure.linsolve()

# postproc
# 1) displace the mesh
dofsol = structure.nodal_dof_solution(store='dofsol')
forces = structure.nodal_forces()
reactions = structure.reaction_forces()

dofsol[i_last, 2], forces[i_last, 2], reactions[0, 2]

(-1.9854154472816174e-05, -100.0, 100.00000173423086)

dofsol[i_last] @ forces[i_last]

0.0019854154472816175

internal_forces = structure.internal_forces()
internal_forces[0, 0, 2]

-99.81153322028695