Truss Implementation #167
Labels
enhancement
New feature or request
Milestone
Comments
|
This is a very first proposal: # -*- coding: utf-8 -*-
"""
Created on Tue Dec 21 22:40:07 2021
@author: adutz
"""
import numpy as np
import felupe as fe
from scipy.sparse import csr_matrix as sparsematrix
mesh = fe.Cube(n=2)
mesh.cells = np.array([[0,7],[0,3],[1,4],[6,2]])
mesh.update(mesh.cells)
class RegionTruss(fe.Region):
def __init__(self, mesh):
super().__init__(mesh, fe.Line(), fe.GaussLegendre(1, 1), grad=False)
self.dX = (mesh.points[mesh.cells[:, 1]] - mesh.points[mesh.cells[:, 0]]).T
self.L = np.sqrt(np.einsum("i...,i...->...", self.dX, self.dX))
self.I = fe.math.identity(self.dX.reshape(*self.dX.shape, 1)).T[0].T
class FieldTruss(fe.Field):
def __init__(self, region, values=0):
super().__init__(region, dim=region.mesh.dim, values=values)
def dx(self):
du = (self.values[mesh.cells[:, 1]] - self.values[mesh.cells[:, 0]]).T
return self.region.dX + du
def l(self):
return np.sqrt(np.einsum("i...,i...->...", self.dx(), self.dx()))
def n(self):
return self.dx() / np.linalg.norm(self.dx(), axis=0)
def m(self):
return fe.math.dya(self.n(), self.n(), mode=1)
def eye(self):
m = self.m()
return fe.math.identity(m.reshape(*m.shape, 1)).T[0].T
def extract(self, grad=True, sym=False, add_identity=True):
stretch = self.l() / self.region.L
if add_identity:
return stretch
else:
return stretch - 1
class LinearElastic1D:
"1D linear elasticity."
def __init__(self, E=None, A=None):
self.E = E
self.A = A
def stress(self, stretch, E=None):
E = E if E is not None else self.E
strain = stretch - 1
return E * strain
def elasticity(self, stretch=None, E=None):
E = E if E is not None else self.E
return E
class IntegralFormTruss:
def __init__(self, fun, v, A, u=None):
self.fun = fun
self.v = v
self.A = A
self.u = u
def integrate(self):
if self.u is None:
n = self.v.n()
r_aie = self.fun * self.A * np.stack((-n, n))
return r_aie
else:
l = self.v.l()
I = self.v.region.I
L = self.v.region.L
m = self.v.m()
k = self.A / L * self.fun * m - self.A / l * (I - m)
k = k.reshape(*m.shape, 1, 1)
K_aibke = np.einsum("ikeab->aibke", np.block([[k, -k],[-k, k]]))
return K_aibke
def assemble(self, values=None):
if values is None:
values = self.integrate()
if self.u is None:
indices = self.v.indices.ai
else:
eai = self.v.indices.eai
ebk = self.u.indices.eai
eaibk0 = np.repeat(eai, ebk.shape[1] * self.u.dim)
eaibk1 = np.tile(ebk, (1, eai.shape[1] * self.v.dim, 1)).ravel()
indices = (eaibk0, eaibk1)
return sparsematrix((values.ravel(), indices))
region = RegionTruss(mesh)
u = FieldTruss(region, values=np.random.rand(*mesh.points.shape) / 10)
stretch = u.extract()
umat = LinearElastic1D(E=1.0, A=1.0)
r = IntegralFormTruss(fun=umat.stress(stretch), v=u, A=1).assemble()
K = IntegralFormTruss(fun=umat.elasticity(stretch), v=u, A=1, u=u).assemble() |
|
rename the attributes for better readability... |
|
2nd version of the proposal: # -*- coding: utf-8 -*-
"""
Created on Tue Dec 21 22:40:07 2021
@author: adutz
"""
import numpy as np
import felupe as fe
from scipy.sparse import csr_matrix as sparsematrix
mesh = fe.Cube(n=2)
mesh.cells = np.array([[0,7],[0,3],[1,4],[6,2]])
mesh.update(mesh.cells)
class RegionTruss(fe.Region):
def __init__(self, mesh):
super().__init__(mesh, fe.Line(), fe.GaussLegendre(1, 1), grad=False)
self.truss = (mesh.points[mesh.cells[:, 1]] - mesh.points[mesh.cells[:, 0]]).T
self.length = np.sqrt(np.einsum("i...,i...->...", self.truss, self.truss))
self.identity = fe.math.identity(self.truss.reshape(*self.truss.shape, 1)).T[0].T
class FieldTruss(fe.Field):
def __init__(self, region, values=0):
super().__init__(region, dim=region.mesh.dim, values=values)
def truss(self):
du = (self.values[mesh.cells[:, 1]] - self.values[mesh.cells[:, 0]]).T
return self.region.truss + du
def length(self):
return np.sqrt(np.einsum("i...,i...->...", self.truss(), self.truss()))
def normal(self):
return self.truss() / np.linalg.norm(self.truss(), axis=0)
def extract(self, grad=True, sym=False, add_identity=True):
stretch = self.length() / self.region.length
if add_identity:
return stretch
else:
return stretch - 1
def grad(self, sym=False):
return self.extract(add_identity=False)
class LinearElastic1D:
"1D linear elasticity."
def __init__(self, E=None, A=None):
self.E = E
self.A = A
def stress(self, stretch, E=None):
E = E if E is not None else self.E
strain = stretch - 1
return E * strain
def elasticity(self, stretch=None, E=None):
E = E if E is not None else self.E
return E
class IntegralFormTruss:
def __init__(self, fun, v, A, u=None):
self.fun = fun
self.v = v
self.A = A
self.u = u
def integrate(self):
n = self.v.normal()
if self.u is None:
r_aie = self.fun * self.A * np.stack((-n, n))
return r_aie
else:
l = self.v.length()
I = self.v.region.identity
L = self.v.region.length
m = fe.math.dya(n, n, mode=1)
k = self.A / L * self.fun * m - self.A / l * (I - m)
k = k.reshape(*m.shape, 1, 1)
K_aibke = np.einsum("ikeab->aibke", np.block([[k, -k],[-k, k]]))
return K_aibke
def assemble(self, values=None, parallel=False):
if values is None:
values = self.integrate()
if self.u is None:
indices = self.v.indices.ai
else:
eai = self.v.indices.eai
ebk = self.u.indices.eai
eaibk0 = np.repeat(eai, ebk.shape[1] * self.u.dim)
eaibk1 = np.tile(ebk, (1, eai.shape[1] * self.v.dim, 1)).ravel()
indices = (eaibk0, eaibk1)
return sparsematrix((values.ravel(), indices))
region = RegionTruss(mesh)
u = FieldTruss(region, values=np.random.rand(*mesh.points.shape) / 10)
stretch = u.extract()
umat = LinearElastic1D(E=1.0, A=1.0)
r = IntegralFormTruss(fun=umat.stress(stretch), v=u, A=1).assemble()
K = IntegralFormTruss(fun=umat.elasticity(stretch), v=u, A=1, u=u).assemble() |
|
final proposal, ready for implementation: # -*- coding: utf-8 -*-
"""
Created on Tue Dec 21 22:40:07 2021
@author: adutz
"""
import numpy as np
import felupe as fe
from scipy.sparse import csr_matrix as sparsematrix
class RegionTruss(fe.Region):
def __init__(self, mesh):
super().__init__(mesh, fe.Line(), fe.GaussLegendre(1, 1), grad=False)
self.truss = (mesh.points[mesh.cells[:, 1]] - mesh.points[mesh.cells[:, 0]]).T
self.length = np.sqrt(np.einsum("i...,i...->...", self.truss, self.truss))
self.identity = fe.math.identity(self.truss.reshape(*self.truss.shape, 1)).T[0].T
class FieldTruss(fe.Field):
def __init__(self, region, values=0):
super().__init__(region, dim=region.mesh.dim, values=values)
def truss(self):
change = (self.values[mesh.cells[:, 1]] - self.values[mesh.cells[:, 0]]).T
return self.region.truss + change
def length(self):
return np.sqrt(np.einsum("i...,i...->...", self.truss(), self.truss()))
def normal(self):
return self.truss() / np.linalg.norm(self.truss(), axis=0)
def extract(self, grad=True, sym=False, add_identity=True):
stretch = self.length() / self.region.length
if add_identity:
return stretch
else:
return stretch - 1
def grad(self, sym=False):
return self.extract(add_identity=False)
class LinearElastic1D:
"1D linear elasticity."
def __init__(self, E=None, A=None):
self.E = E
self.A = A
def stress(self, stretch, E=None):
E = E if E is not None else self.E
return E * (stretch - 1)
def elasticity(self, stretch=None, E=None):
E = E if E is not None else self.E
return E
class IntegralFormTruss:
def __init__(self, fun, v, A, u=None):
self.fun = fun
self.v = v
self.A = A
self.u = u
def integrate(self):
n = self.v.normal()
if self.u is None:
r_aie = self.fun * self.A * np.stack((-n, n))
return r_aie
else:
l = self.v.length()
I = self.v.region.identity
L = self.v.region.length
m = fe.math.dya(n, n, mode=1)
k = self.A / L * self.fun * m - self.A / l * (I - m)
k = k.reshape(*m.shape, 1, 1)
K_aibke = np.einsum("ikeab->aibke", np.block([[k, -k],[-k, k]]))
return K_aibke
def assemble(self, values=None, parallel=False):
if values is None:
values = self.integrate()
if self.u is None:
indices = self.v.indices.ai
else:
eai = self.v.indices.eai
ebk = self.u.indices.eai
eaibk0 = np.repeat(eai, ebk.shape[1] * self.u.dim)
eaibk1 = np.tile(ebk, (1, eai.shape[1] * self.v.dim, 1)).ravel()
indices = (eaibk0, eaibk1)
return sparsematrix((values.ravel(), indices))
mesh = fe.Cube(n=2)
mesh.cells = np.array([[0,7],[0,3],[1,4],[6,2]])
mesh.update(mesh.cells)
mesh.cell_type = "line"
region = RegionTruss(mesh)
u = FieldTruss(region)
umat = LinearElastic1D(E=1.0, A=1.0)
r = IntegralFormTruss(fun=umat.stress(u.extract()), v=u, A=umat.A).assemble()
K = IntegralFormTruss(fun=umat.elasticity(u.extract()), v=u, A=umat.A, u=u).assemble()
bounds = {
"fixed": fe.Boundary(u, fx=lambda x: x==0),
"right": fe.Boundary(u, fx=lambda x: x==1, skip=(1, 0, 0)),
"move": fe.Boundary(u, fx=lambda x: x==1, skip=(0, 1, 1), value=0.1),
}
dof0, dof1 = fe.dof.partition(u, bounds)
ext0 = fe.dof.apply(u, bounds, dof0)
system = fe.solve.partition(u, K, dof1, dof0, r)
u += fe.solve.solve(*system, ext0)
r = IntegralFormTruss(fun=umat.stress(u.extract()), v=u, A=umat.A).assemble()
np.linalg.norm(r.toarray()[dof1])
fe.save(
region,
u,
cell_data={"force": [umat.stress(u.extract()) * umat.A]}
) |
|
Although this implementation just works as expected I think it focusses too much on structural mechanics. Will try to re-formulate this in the future to be more general (for 1d and 2d problems, eventually with unit normal and edge vectors, etc.) |
|
I think this will be a great addition to the mechanics-namespace as SolidBodyTruss. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add a truss element (see TrussPy)
(Truss element from TrussPy)
Todo:
RegionTrussThe text was updated successfully, but these errors were encountered: