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.temperatures[mesh.point_sets["temp_left"], 0] = 0.0
box_heat.temperatures[mesh.point_sets["temp_right"], 0] = 20.0
# box_heat.heat_flux[nodes[:, 1] == 1.0, 0] = 0.0
# box_heat.heat_flux[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)
────────────────────────────────────────────────────────────────────────────────────────
model SolidHeat · 360 elem · 584 dof · float64
machine x86_64 · 2 threads · 16 GB RAM
solver direct · scipy · cpu
newton rtol 1e-08 · atol 1e-06 · <=100 it · dt <= 0.2
────────────────────────────────────────────────────────────────────────────────────────
Time step Time Steps Newton Krylov Residual Wall time
1 0.2 1 1 - 8.66e-15 0.01 s
2 0.4 1 1 - 5.74e-15 0.01 s
3 0.6 1 1 - 4.80e-15 0.01 s
4 0.8 1 1 - 3.46e-15 0.01 s
5 1 1 1 - 3.60e-15 0.01 s
6 1.2 1 1 - 3.41e-15 0.01 s
7 1.4 1 1 - 2.75e-15 0.01 s
8 1.6 1 1 - 3.56e-15 0.01 s
9 1.8 1 1 - 3.35e-15 0.01 s
10 2 1 1 - 3.28e-15 0.01 s
11 2.2 1 1 - 2.84e-15 0.01 s
12 2.4 1 1 - 3.03e-15 0.01 s
13 2.6 1 1 - 3.03e-15 0.01 s
14 2.8 1 1 - 2.94e-15 0.01 s
15 3 1 1 - 3.35e-15 0.01 s
16 3.2 1 1 - 3.14e-15 0.01 s
17 3.4 1 1 - 2.95e-15 0.01 s
18 3.6 1 1 - 3.22e-15 0.01 s
19 3.8 1 1 - 3.51e-15 0.01 s
20 4 1 1 - 3.23e-15 0.01 s
21 4.2 1 1 - 3.05e-15 0.01 s
22 4.4 1 1 - 3.34e-15 0.01 s
23 4.6 1 1 - 3.18e-15 0.01 s
24 4.8 1 1 - 3.44e-15 0.01 s
25 5 1 1 - 3.32e-15 0.01 s
────────────────────────────────────────────────────────────────────────────────────────
converged · 25 time steps · 25 Newton iterations · 0.16 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]:
