In [1]:
import meshio as mo
import pyvista as pv
import torch
from torchfem import SolidHeat
from torchfem.data import get_data
from torchfem.materials import IsotropicConductivity3D
from torchfem.plot_utils import embed_pyvista_animation
# Set default data type to double precision
torch.set_default_dtype(torch.float64)
3D unstructured mesh with linear hexahedrons read from Abaqus input file¶
In [2]:
kappa, rho, cp = 500.0, 7850.0, 500.0
In [3]:
inp_path = get_data("extruded_hex_mesh.inp")
# material
material = IsotropicConductivity3D(kappa=kappa, rho=rho * cp)
mesh = mo.read(inp_path)
# Create model
nodes = torch.tensor(mesh.points) * 1.0e-3 # SI
elements = torch.tensor(mesh.cells_dict["hexahedron"])
box_heat = SolidHeat(nodes, elements, material)
# Assign boundary conditions
box_heat.constraints[mesh.point_sets["temp_left"]] = True
box_heat.constraints[mesh.point_sets["temp_right"]] = True
# box_heat.constraints[mesh.point_sets["insulation"]] = True
box_heat.displacements[mesh.point_sets["temp_left"], 0] = 0.0
box_heat.displacements[mesh.point_sets["temp_right"], 0] = 20.0
# box_heat.forces[nodes[:, 1] == 1.0, 0] = 0.0
# box_heat.forces[mesh.point_sets["insulation"]] = 5.0
# solver settings
end_time = 5.0
delta_t = 0.2
# Results are returned at exactly these times
times = torch.arange(0.0, end_time + delta_t, delta_t)
# Solve
temp, rfl, hf, temp_grad, _ = box_heat.time_integration(times, delta_t, verbose=True)
─── torch-fem · time integration ───────────────────────────────────────────────────────
model SolidHeat · 360 elem · 584 dof · float64
machine AMD EPYC 7763 64-Core Processor · 2 threads · 16 GB RAM
solver spsolve · direct · scipy · cpu
newton rtol 1e-08 · atol 1e-06 · ≤100 it · Δt ≤ 0.2
────────────────────────────────────────────────────────────────────────────────────────
Time step Time Steps Iterations Residual Wall time
1 0.2 1 1 3.70e-15 0.01 s
2 0.4 1 1 4.44e-15 0.01 s
3 0.6 1 1 3.26e-15 0.01 s
4 0.8 1 1 2.83e-15 0.01 s
5 1 1 1 2.52e-15 0.01 s
6 1.2 1 1 2.89e-15 0.01 s
7 1.4 1 1 2.76e-15 0.01 s
8 1.6 1 1 2.63e-15 0.01 s
9 1.8 1 1 2.95e-15 0.01 s
10 2 1 1 3.09e-15 0.01 s
11 2.2 1 1 2.91e-15 0.01 s
12 2.4 1 1 3.06e-15 0.01 s
13 2.6 1 1 2.84e-15 0.01 s
14 2.8 1 1 3.22e-15 0.01 s
15 3 1 1 3.38e-15 0.01 s
16 3.2 1 1 2.81e-15 0.01 s
17 3.4 1 1 3.57e-15 0.01 s
18 3.6 1 1 2.94e-15 0.01 s
19 3.8 1 1 3.77e-15 0.01 s
20 4 1 1 3.71e-15 0.01 s
21 4.2 1 1 3.47e-15 0.01 s
22 4.4 1 1 3.47e-15 0.01 s
23 4.6 1 1 3.18e-15 0.01 s
24 4.8 1 1 3.45e-15 0.01 s
25 5 1 1 2.96e-15 0.01 s
────────────────────────────────────────────────────────────────────────────────────────
converged · 25 time steps · 25 iterations · 0.31 s
Plot the animated results.
In [4]:
pl = pv.Plotter()
def animate_temp(i):
pl.clear()
box_heat.plot(
node_property={"Temperature": temp[i]},
show_undeformed=True,
opacity=torch.linspace(0.2, 1.0, len(temp)),
colormap="magma",
clim=[temp.min(), temp.max()],
plotter=pl,
lighting=False,
)
frames = range(len(temp))
embed_pyvista_animation(pl, animate_temp, frames)
Out[4]:
