import matplotlib.pyplot as plt
import torch
from tqdm import tqdm

from torchfem import Shell
from torchfem.materials import OrthotropicElasticityPlaneStress
from torchfem.mesh import rect_tri
from torchfem.rotations import planar_rotation

torch.set_default_dtype(torch.float64)

G_t = 10000.0 / (2.0 * (1 + 0.3))
mat = OrthotropicElasticityPlaneStress(
    E_1=100000.0, E_2=10000.0, nu_12=0.1, G_12=5000.0, G_13=G_t, G_23=G_t
)

# Parameters
L = 100.0
F = -10.0
Zm = 100.0

# Define plate
nodes, elements = rect_tri(15, 15, 0.5 * L, 0.5 * L, variant="zigzag")
nodes = torch.hstack([nodes, torch.zeros((len(nodes), 1))])

# Create Shell model
plate = Shell(nodes, elements, mat)

# Boundaries
bottom = nodes[:, 1] < 0.01
left = nodes[:, 0] < 0.01

# Boundary conditions
plate.forces[0, 2] = F
plate.constraints[-1] = True
plate.constraints[left, 0] = True
plate.constraints[left, 4] = True
plate.constraints[left, 5] = True
plate.constraints[bottom, 1] = True
plate.constraints[bottom, 3] = True
plate.constraints[bottom, 5] = True

# Solve
u, f, σ, _, _ = plate.solve()


plate.plot(
    u[:, 0:3],
    node_property={"u": torch.linalg.norm(u[:, 0:3], dim=1)},
    mirror=(True, True, False),
    bcs=True,
)

EmbeddableWidget(value='<iframe srcdoc="<!doctype html>\n<html lang=&quot;en&quot;>\n  <head>\n    <meta chars…

def target_function(phi):
    # Recompute stiffnesses due to changed orientations
    R = planar_rotation(phi)
    plate.material = mat.vectorize(plate.n_elem).rotate(R)
    # Solve
    u, f, _, _, _ = plate.solve(differentiable_parameters=phi)
    # Return compliance
    return torch.inner(u.ravel(), f.ravel())

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()

100%|██████████| 100/100 [00:09<00:00, 10.22it/s]

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

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

# Compute direction
loc_dir = torch.stack([torch.cos(phi), -torch.sin(phi), torch.zeros_like(phi)]).T
dir = torch.einsum("...ij, ...j", plate.t.transpose(1, 2), loc_dir)

# Show optimized displacements and orientations
plate.plot(
    u[:, 0:3],
    orientations=dir.unsqueeze(1),
    node_property={"u": torch.linalg.norm(u[:, 0:3], dim=1)},
)

EmbeddableWidget(value='<iframe srcdoc="<!doctype html>\n<html lang=&quot;en&quot;>\n  <head>\n    <meta chars…

Anisotropic orientation optimization¶

Optimization¶

Target function is the strain energy¶

The optimization¶