Orientation optimization¶
In [1]:
import matplotlib.pyplot as plt
import torch
from tqdm import tqdm
from torchfem.rotations import planar_rotation
from torchfem.examples import get_example_file
from torchfem.io import import_mesh
from torchfem.materials import OrthotropicElasticityPlaneStress
torch.set_default_dtype(torch.float64)
In [2]:
mat = OrthotropicElasticityPlaneStress(
E_1=100000.0, E_2=10000.0, nu_12=0.1, G_12=5000.0
)
In [3]:
# Import mesh
plate = import_mesh(get_example_file("plate_hole.vtk"), mat)
plate.plot()
Boundary conditions¶
In [4]:
plate.constraints[plate.nodes[:, 0] == 0.0] = True
plate.forces[plate.nodes[:, 0] >= 0.199] = 10.0
plate.forces[:, 1] = 0.0
plate.plot()
Solve initial configuration¶
In [5]:
u, f, sigma, epsilon, state = plate.solve()
plate.plot(u=u, node_property=u[:, 0])
Optimization¶
Target function is the strain energy¶
In [6]:
def target_function(phi):
# Recompute stiffness due to changed orientations
R = planar_rotation(phi)
plate.material = mat.vectorize(plate.n_elem).rotate(R)
# Solve
u, f, _, _, _ = plate.solve()
# Return compliance
return torch.inner(u.ravel(), f.ravel())
The optimization¶
In [7]:
phi = torch.zeros((len(plate.elements)), requires_grad=True)
optimizer = torch.optim.Adam([phi], lr=0.1)
energies = []
for _ in tqdm(range(100)):
optimizer.zero_grad()
objective = target_function(phi)
energies.append(objective.detach().item())
objective.backward()
optimizer.step()
0%| | 0/100 [00:00<?, ?it/s]
3%|▎ | 3/100 [00:00<00:04, 20.17it/s]
6%|▌ | 6/100 [00:00<00:04, 20.41it/s]
9%|▉ | 9/100 [00:00<00:04, 20.01it/s]
12%|█▏ | 12/100 [00:00<00:04, 19.77it/s]
14%|█▍ | 14/100 [00:00<00:04, 19.75it/s]
17%|█▋ | 17/100 [00:00<00:04, 20.17it/s]
20%|██ | 20/100 [00:00<00:03, 20.59it/s]
23%|██▎ | 23/100 [00:01<00:03, 20.72it/s]
26%|██▌ | 26/100 [00:01<00:03, 21.04it/s]
29%|██▉ | 29/100 [00:01<00:03, 21.09it/s]
32%|███▏ | 32/100 [00:01<00:03, 21.11it/s]
35%|███▌ | 35/100 [00:01<00:03, 20.92it/s]
38%|███▊ | 38/100 [00:01<00:02, 20.90it/s]
41%|████ | 41/100 [00:01<00:02, 20.91it/s]
44%|████▍ | 44/100 [00:02<00:02, 20.94it/s]
47%|████▋ | 47/100 [00:02<00:02, 20.87it/s]
50%|█████ | 50/100 [00:02<00:02, 20.79it/s]
53%|█████▎ | 53/100 [00:02<00:02, 20.56it/s]
56%|█████▌ | 56/100 [00:02<00:02, 20.42it/s]
59%|█████▉ | 59/100 [00:02<00:01, 20.61it/s]
62%|██████▏ | 62/100 [00:03<00:01, 20.59it/s]
65%|██████▌ | 65/100 [00:03<00:01, 20.76it/s]
68%|██████▊ | 68/100 [00:03<00:01, 20.69it/s]
71%|███████ | 71/100 [00:03<00:01, 20.85it/s]
74%|███████▍ | 74/100 [00:03<00:01, 20.67it/s]
77%|███████▋ | 77/100 [00:03<00:01, 20.51it/s]
80%|████████ | 80/100 [00:03<00:00, 20.61it/s]
83%|████████▎ | 83/100 [00:04<00:00, 20.69it/s]
86%|████████▌ | 86/100 [00:04<00:00, 20.87it/s]
89%|████████▉ | 89/100 [00:04<00:00, 20.83it/s]
92%|█████████▏| 92/100 [00:04<00:00, 20.83it/s]
95%|█████████▌| 95/100 [00:04<00:00, 20.61it/s]
98%|█████████▊| 98/100 [00:04<00:00, 20.63it/s]
100%|██████████| 100/100 [00:04<00:00, 20.67it/s]
In [8]:
plt.plot(energies, ".-k")
plt.title("Optimization history")
plt.xlabel("Iteration")
plt.ylabel("Compliance")
plt.show()
In [9]:
# Compute optimized displacements
u, f, _, _, _ = plate.solve()
# Plot orientations
plate.plot(u=u, node_property=u[:, 0], linewidth=0.2, orientation=phi)
