Cubes with different element types

All unit cubes are clamped at one end and subjected to a uniaxial strain.

import meshzoo
import torch

from torchfem import Solid
from torchfem.mesh import cube_hexa
from torchfem.elements import linear_to_quadratic
from torchfem.materials import IsotropicElasticity3D

# Set default data type to double precision
torch.set_default_dtype(torch.float64)

# Elastic material model
material = IsotropicElasticity3D(E=1000.0, nu=0.3)

3D Cube with linear hexahedrons

# Generate cube
nodes, elements = cube_hexa(5, 5, 5)

# Create model
box = Solid(nodes, elements, material)

# Assign boundary conditions
box.constraints[nodes[:, 0] == 0.0, :] = True
box.constraints[nodes[:, 0] == 1.0, 0] = True
box.displacements[nodes[:, 0] == 1.0, 0] = 0.1


# Solve
u, f, σ, F, α = box.solve()

# Plot
box.plot(u=u, node_property={"Disp": u})

3D Cube with quadratic hexahedrons

# Upgrade elements to quadratic
nodes, elements = linear_to_quadratic(nodes, elements)

# Create model
box = Solid(nodes, elements, material)

# Assign boundary conditions
box.constraints[nodes[:, 0] == 0.0, :] = True
box.constraints[nodes[:, 0] == 1.0, 0] = True
box.displacements[nodes[:, 0] == 1.0, 0] = 0.1

# Solve
u, f, σ, F, α = box.solve()

# Plot
box.plot(u=u, node_property={"Disp": u})

3D Cube with linear tetrahedrons

# Generate cube
points, cells = meshzoo.cube_tetra(
    torch.linspace(0.0, 1.0, 3),
    torch.linspace(0.0, 1.0, 3),
    torch.linspace(0.0, 1.0, 3),
)
nodes = torch.tensor(points)
elements = torch.tensor(cells.tolist())

box = Solid(nodes, elements, material)

# Assign boundary conditions
box.constraints[nodes[:, 0] == 0.0, :] = True
box.constraints[nodes[:, 0] == 1.0, 0] = True
box.displacements[nodes[:, 0] == 1.0, 0] = 0.1

u, f, σ, F, α = box.solve()

box.plot(u=u, node_property={"Disp": u})

3D Cube with quadratic tetrahedrons

nodes, elements = linear_to_quadratic(nodes, elements)

box = Solid(nodes, elements, material)

# Assign boundary conditions
box.constraints[nodes[:, 0] == 0.0, :] = True
box.constraints[nodes[:, 0] == 1.0, 0] = True
box.displacements[nodes[:, 0] == 1.0, 0] = 0.1

u, f, σ, F, α = box.solve()

box.plot(u=u, node_property={"Disp": u})