In [1]:
import matplotlib.pyplot as plt
import meshio
import numpy as np
import pyvista
import torch
from scipy.optimize import bisect
from scipy.spatial import KDTree
from tqdm import tqdm
from torchfem import Assembly, ReferencePoint, Solid
from torchfem.data import get_data
from torchfem.materials import IsotropicElasticity3D
from torchfem.plot_utils import show_plotter
torch.set_default_dtype(torch.float64)
In [2]:
# Material model (Ti-6Al-4V) in imperial units
material = IsotropicElasticity3D(E=16500.0, nu=0.342)
# Qualitative colormap for domains
cmap = plt.get_cmap("tab10", 7)
Mesh¶
The domain is meshed with linear tetrahedrons in gmsh and subdivided in seven geometrical domains.
In [3]:
mesh = meshio.read(get_data("ge_bracket.vtu"))
elements = torch.tensor(mesh.cells[0].data)
nodes = torch.tensor(mesh.points)
domain = torch.tensor(mesh.cell_data["gmsh:geometrical"][0])
# Assembly of the solid model, a reference point in the bore and one per bolt
solid = Solid(nodes, elements, material)
pin = ReferencePoint([0.0, 0.0, 0.0])
mnts = [ReferencePoint(nodes[elements[domain == d]].mean((0, 1))) for d in [1, 3, 4, 5]]
bracket = Assembly([solid, pin, *mnts])
The reference points are coupled to the inner surfaces of the bore holes via kinematic coupling. Rotation around the pin or bold axes is free.
In [4]:
# Couple each mount to the inner surface of its hole.
R_hole = 6.0
for mnt in mnts:
hole = torch.linalg.norm(nodes[:, :2] - mnt.nodes[0, :2], dim=1) < R_hole
bracket.coupling(solid, hole, mnt)
mnt.constraints[0, :] = True
mnt.constraints[0, 5] = False
# Couple reference point to inner nodes at bore.
R_bore = 10.0
interface = torch.sqrt(nodes[:, 1] ** 2 + nodes[:, 2] ** 2) < R_bore
bracket.coupling(solid, interface, pin, dofs=[1, 2])
pin.constraints[0, 0] = True
bracket.plot(element_property={"Domain": domain}, cmap=cmap)
EmbeddableWidget(value='<iframe srcdoc="<!doctype html>\n<html lang="en">\n <head>\n <meta chars…
In [5]:
# The four load cases of the challenge.
load_cases = torch.zeros(4, 1, 6)
load_cases[0, 0, 2] = 8000.0
load_cases[1, 0, 1] = -8500.0
load_cases[2, 0, 1] = -9500 * np.sin(np.deg2rad(42))
load_cases[2, 0, 2] = 9500 * np.cos(np.deg2rad(42))
load_cases[3, 0, 5] = 5000.0
# Apply load case 1 for testing
pin.forces[:] = load_cases[0]
# Solve
u, f, σ, F, α = bracket.solve(rtol=0.01, verbose=True)
--- torch-fem | solve ------------------------------------------------------------------
model Assembly | 6 parts | 214,524 elem | 127,674 dof | float64
machine Apple M1 Pro | 8 threads | 16 GB RAM
solver cg | iterative | amg | scipy | cpu
newton rtol 1e-02 | atol 1e-06 | <=10 it
----------------------------------------------------------------------------------------
Increment Load factor Steps Iterations Residual Wall time
1 1 1 1 7.85e-07 4.78 s
----------------------------------------------------------------------------------------
converged | 1 increment | 1 iteration | 4.79 s
In [6]:
# Plot, with the pin and its coupling drawn on top
bracket.plot(u=u, node_property={"u": u})
EmbeddableWidget(value='<iframe srcdoc="<!doctype html>\n<html lang="en">\n <head>\n <meta chars…
Optimization parameters¶
We define the optimization parameters, i.e. the volume fraction, the penalization factor, the move limit, the filter radius, and the stiffness of void relative to solid material.
The stiffness is interpolated as
$$ \mathbf{C}(\rho) = \mathbf{C}_\text{min} + \rho^p (\mathbf{C}_\text{max} - \mathbf{C}_\text{min}) $$ with $$ \mathbf{C}_\text{min} = \text{soft} \cdot \mathbf{C}_\text{max} $$
Void keeps a fixed fraction of the solid stiffness, so the density is free to approach zero without the stiffness following it down.
In [7]:
volfrac = 0.15
p = 3
move = 0.2
R = 5.0
soft = 1e-3
In [8]:
# Design elements are only the ones in domain 6
design_elements = elements[domain == 6]
vols = solid.integrate_field()[domain == 6]
In [9]:
# Initial, minimum, and maximum values of design variables.
rho_0 = volfrac * torch.ones(len(design_elements))
rho_min = 1e-3 * torch.ones_like(rho_0)
rho_max = torch.ones_like(rho_0)
# Volume fraction
V_0 = volfrac * vols.sum()
# Analytical gradient of the stiffness matrix
k0 = solid.k0()[domain == 6].clone()
C0 = solid.material.C[domain == 6].clone()
C_min = soft * C0
# Precompute filter weights using scipy's KDTree
ecenters = torch.mean(nodes[design_elements], dim=1)
tree = KDTree(ecenters)
d = tree.sparse_distance_matrix(tree, R, output_type="coo_matrix")
# Convert to torch sparse tensor
H = torch.sparse_coo_tensor(
indices=torch.stack([torch.as_tensor(d.row), torch.as_tensor(d.col)]),
values=R - torch.as_tensor(d.data),
size=(len(design_elements), len(design_elements)),
)
H_sum = H.sum(dim=0).to_dense()
Optimization with optimality constraints.¶
This should take around 10 minutes to run...
In [10]:
rho = [rho_0]
history = []
# Iterate solutions
for k in tqdm(range(20)):
# Interpolate stiffness between void and solid
solid.material.C[domain == 6] = C_min + torch.einsum(
"n,nijkl->nijkl", rho[k] ** p, C0 - C_min
)
sensitivity = torch.zeros_like(rho[k])
compliance = 0.0
# Iterate over load cases
for lc in load_cases:
# Apply load case at the pin
pin.forces[:] = lc
# Compute solution
u_k, f_k, _, _, _ = bracket.solve(rtol=0.01)
# Evaluation of compliance
compliance += torch.inner(f_k[0].ravel(), u_k[0].ravel())
# Compute analytical sensitivities
u_j = u_k[0][design_elements].reshape(len(design_elements), -1)
w_k = torch.einsum("...i, ...ij, ...j", u_j, k0, u_j)
sensitivity += -p * rho[k] ** (p - 1.0) * w_k
# Filter sensitivities (if r provided)
sensitivity = H @ (rho[k] * sensitivity / vols) / H_sum / (rho[k] / vols)
# For a certain value of mu, apply the iteration scheme
def make_step(mu):
G_k = -sensitivity / mu
upper = torch.min(rho_max, (1 + move) * rho[k])
lower = torch.max(rho_min, (1 - move) * rho[k])
rho_trial = G_k**0.5 * rho[k]
return torch.maximum(torch.minimum(rho_trial, upper), lower)
# Constraint function
def g(mu):
rho_k = make_step(mu)
return torch.inner(rho_k, vols) - V_0
# Find the root of g(mu)
mu = bisect(g, 1e-10, 100.0)
# Update design history
rho.append(make_step(mu))
history.append(compliance)
100%|██████████| 20/20 [08:15<00:00, 24.77s/it]
In [11]:
plt.plot(history)
plt.xlabel("Iteration")
plt.ylabel("Summed Compliance")
plt.grid()
plt.xlim(0, len(history) - 1)
plt.ylim(bottom=0)
plt.show()
Result¶
In [12]:
# Non-design elements are solid.
rho_final = torch.ones(len(elements))
rho_final[domain == 6] = (H @ rho[-1]) / H_sum
plotter = pyvista.Plotter()
solid.plot(plotter=plotter, opacity=0.2, show_edges=False)
solid.plot(
plotter=plotter,
element_property={"rho": rho_final},
clip=("rho", 0.5),
show_edges=False,
color="lightblue",
)
show_plotter(plotter)
EmbeddableWidget(value='<iframe srcdoc="<!doctype html>\n<html lang="en">\n <head>\n <meta chars…