Orientation optimization

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)

mat = OrthotropicElasticityPlaneStress(
    E_1=100000.0, E_2=10000.0, nu_12=0.1, G_12=5000.0
)

# Import mesh
plate = import_mesh(get_example_file("plate_hole.vtk"), mat)
plate.plot()

Boundary conditions

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

u, f, sigma, F, state = plate.solve()
plate.plot(u=u, node_property=u[:, 0])

Optimization

Target function is the strain energy

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

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]

  5%|▌         | 5/100 [00:00<00:02, 44.43it/s]

 10%|█         | 10/100 [00:00<00:02, 44.70it/s]

 15%|█▌        | 15/100 [00:00<00:01, 42.90it/s]

 20%|██        | 20/100 [00:00<00:01, 43.55it/s]

 25%|██▌       | 25/100 [00:00<00:01, 44.01it/s]

 30%|███       | 30/100 [00:00<00:01, 44.24it/s]

 35%|███▌      | 35/100 [00:00<00:01, 44.43it/s]

 40%|████      | 40/100 [00:00<00:01, 44.48it/s]

 45%|████▌     | 45/100 [00:01<00:01, 44.33it/s]

 50%|█████     | 50/100 [00:01<00:01, 44.42it/s]

 55%|█████▌    | 55/100 [00:01<00:01, 44.42it/s]

 60%|██████    | 60/100 [00:01<00:00, 44.53it/s]

 65%|██████▌   | 65/100 [00:01<00:00, 44.44it/s]

 70%|███████   | 70/100 [00:01<00:00, 44.39it/s]

 75%|███████▌  | 75/100 [00:01<00:00, 44.37it/s]

 80%|████████  | 80/100 [00:01<00:00, 44.60it/s]

 85%|████████▌ | 85/100 [00:01<00:00, 44.66it/s]

 90%|█████████ | 90/100 [00:02<00:00, 44.88it/s]

 95%|█████████▌| 95/100 [00:02<00:00, 44.89it/s]

100%|██████████| 100/100 [00:02<00:00, 45.10it/s]

100%|██████████| 100/100 [00:02<00:00, 44.47it/s]

plt.plot(energies, ".-k")
plt.title("Optimization history")
plt.xlabel("Iteration")
plt.ylabel("Compliance")
plt.show()

# Compute optimized displacements
u, f, _, _, _ = plate.solve()

# Plot orientations
plate.plot(u=u, node_property=u[:, 0], linewidth=0.2, orientation=phi)