Create animations

[adtzlr]

import felupe as fem
import pypardiso

radius = 60
height = 20
thickness = 4
layers = 3
details = 6
compression = 15

m = 1 + int(1 + height / (radius / (2 * (details - 1))))

circle = fem.Circle(radius=radius, n=details)
cylinder = circle.expand(n=m, z=height).add_runouts(
    normalize=True,
    axis=2,
    centerpoint=[0, 0, height / 2],
    exponent=3,
    values=[0.075, 0.075],
)
plate = circle.expand(n=2, z=thickness)
plates = [
    plate.translate(layer * height + thickness * (layer - 2), axis=2)
    for layer in range(1 + layers)
]
cylinders = [
    cylinder.translate(layer * height + thickness * (layer - 1), axis=2)
    for layer in range(layers)
]

plates = fem.MeshContainer(plates, merge=True).stack()
cylinders = fem.MeshContainer(cylinders, merge=True).stack()

container = fem.MeshContainer([plates, cylinders], merge=True)
container.plot(opacity=1.0, colors=["grey", "green"]).show()

field = fem.Field.from_mesh_container(container).as_container()
field = fem.FieldContainer([fem.Field(fem.RegionHexahedron(container.stack()), dim=3)])

regions = [
    fem.RegionHexahedron(container.meshes[0]),
    fem.RegionHexahedron(container.meshes[1]),
]
fields = [
    fem.FieldContainer([fem.Field(regions[0], dim=3)]),
    fem.FieldContainer([fem.Field(regions[1], dim=3)]),
]

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True, axis=2, sym=False)
umats = [
    fem.LinearElasticLargeStrain(E=2.1e5, nu=0.3),
    fem.NeoHooke(mu=1),
]
solids = [
    fem.SolidBody(umat=umats[0], field=fields[0]),
    fem.SolidBodyNearlyIncompressible(umat=umats[1], field=fields[1], bulk=5000),
]

move = fem.math.linsteps([0, -compression, 0], num=25)
ramp = {boundaries["move"]: move}
step = fem.Step(items=solids, ramp=ramp, boundaries=boundaries)

plotter = solids[0].plot("Principal Values of Cauchy Stress", off_screen=True)
plotter.open_gif("result.gif", fps=25)


def save(i, j, substep):
    global plotter
    plotter.clear()
    plotter = solids[0].plot(
        show_undeformed=False,
        show_edges=True,
        color="grey",
        opacity=1.0,
        plotter=plotter,
    )
    plotter = solids[1].plot(
        "Principal Values of Cauchy Stress",
        plotter=plotter,
        show_undeformed=False,
        clim=[-1, 1],
        below_color="black",
        above_color="purple",
    )
    plotter.add_text(
        f"FElupe v{fem.__version__}\n\nStep {i+1}, Substep {j+1}", font_size=8
    )
    plotter.write_frame()


job = fem.Job(steps=[step], callback=save)
job.evaluate(x0=field, tol=1e-1, solver=pypardiso.spsolve)

plotter.close()

[adtzlr]

It's also possible to evaluate only a few number of substeps and interpolate the displacements / strains on any signal.

result.mp4

import felupe as fem
import pypardiso
import numpy as np

radius = 60
height = 20
thickness = 4
layers = 2
details = 11
compression = 5

m = 1 + int(1 + height / (radius / (2 * (details - 1))))

circle = fem.Circle(radius=radius, n=details)
cylinder = circle.expand(n=m, z=height).add_runouts(
    normalize=True,
    axis=2,
    centerpoint=[0, 0, height / 2],
    exponent=3,
    values=[0.075, 0.075],
)
plate = circle.expand(n=2, z=thickness)
plates = [
    plate.translate(layer * height + thickness * (layer - 2), axis=2)
    for layer in range(1 + layers)
]
cylinders = [
    cylinder.translate(layer * height + thickness * (layer - 1), axis=2)
    for layer in range(layers)
]

plates = fem.MeshContainer(plates, merge=True).stack()
cylinders = fem.MeshContainer(cylinders, merge=True).stack()

container = fem.MeshContainer([plates, cylinders], merge=True)
# container.plot(opacity=1.0, colors=["grey", "green"]).show()

field = fem.Field.from_mesh_container(container).as_container()
field = fem.FieldContainer([fem.Field(fem.RegionHexahedron(container.stack()), dim=3)])

regions = [
    fem.RegionHexahedron(container.meshes[0]),
    fem.RegionHexahedron(container.meshes[1]),
]
fields = [
    fem.FieldContainer([fem.Field(regions[0], dim=3)]),
    fem.FieldContainer([fem.Field(regions[1], dim=3)]),
]

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True, axis=2, sym=False)
umats = [
    fem.LinearElasticLargeStrain(E=2.1e5, nu=0.3),
    fem.NeoHooke(mu=1),
]
solids = [
    fem.SolidBody(umat=umats[0], field=fields[0]),
    fem.SolidBodyNearlyIncompressible(umat=umats[1], field=fields[1], bulk=500),
]

move = fem.math.linsteps([0, -compression], num=3)
ramp = {boundaries["move"]: move}
step = fem.Step(items=solids, ramp=ramp, boundaries=boundaries)

values = []


def save(i, j, substep):
    values.append(substep.x[0].values)


job = fem.Job(steps=[step], callback=save)
job.evaluate(x0=field, tol=1e-1, solver=pypardiso.spsolve, parallel=True)
values = np.array(values)
clim = [0, field.evaluate.log_strain(tensor=False).mean(-2)[-1].max()]

from scipy.interpolate import griddata
from tqdm import tqdm

plotter = field.plot(
    "Principal Values of Logarithmic Strain",
    clim=clim,
    off_screen=True,
    show_undeformed=False,
)
plotter.open_movie("result.mp4", framerate=24)

time = fem.math.linsteps([0, 1, 2, 3, 4], num=24)
signal = -compression / 2 + np.sin(2 * np.pi * time) * compression / 2

for xi in tqdm(signal):
    field[0].values = griddata(move, values, xi, method="cubic")
    data = field.evaluate.log_strain(tensor=False).mean(-2)[-1]

    plotter.mesh.points[:] = field.region.mesh.points + field[0].values
    plotter.mesh[plotter.mesh.active_scalars_info.name] = data

    plotter.write_frame()

plotter.close()