Examples

Elastic Hex20

import meshio

import pyfebio as feb

# import a mesh from gmsh format
from_gmsh = meshio.gmsh.read("../../assets/gmsh/hex20.msh")
# use translation function to convert meshio object to an febio mesh
mesh = feb.mesh.translate_meshio(from_gmsh)

# creat a base febio model
my_model = feb.model.Model(mesh=mesh)

# loop over the discovered Elements (parts) and assign materials
# and solid domains
for i, part in enumerate(my_model.mesh.elements):
    mat = feb.material.NeoHookean(
        name=part.name,
        id=i + 1,
        E=feb.material.MaterialParameter(text=1.0 * (i + 1)),
        v=feb.material.MaterialParameter(text=0.3),
    )
    my_model.material.add_material(mat)
    my_model.mesh_domains.add_solid_domain(
        feb.meshdomains.SolidDomain(name=part.name, mat=part.name)
    )

# fix the bottom nodes in space
fix_bottom = my_model.boundary.add_bc(
    feb.boundary.BCZeroDisplacement(node_set="bottom", x_dof=1, y_dof=1, z_dof=1)
)

# let's twist the top nodes by pi radians about their central z-axis
twist_top = my_model.boundary.add_bc(
    feb.boundary.BCRigidDeformation(
        node_set="top", pos="0.5,0.5,0.0", rot=feb.boundary.Value(lc=1, text="0.0,0.0,3.14")
    )
)

# the load curve for our twist
my_model.load_data.add_load_curve(
    feb.loaddata.LoadCurve(id=1, points=feb.loaddata.CurvePoints(points=["0.0,0.0", "1.0,1.0"]))
)

# save the model
my_model.save("elastic_hex20.feb")

# run the model
feb.model.run_model("elastic_hex20.feb")

Biphasic Hex20

from meshio.gmsh import read

import pyfebio as feb

# import mesh from gmsh format
from_gmsh = read("../../assets/gmsh/hex20.msh")
# translate meshio mesh to febio mesh
mesh = feb.mesh.translate_meshio(from_gmsh)

# initialize a biphasic febio modelk
my_model = feb.model.BiphasicModel(mesh=mesh)

# loop over Elements (parts) and assign biphasic materials
# and also assign solid domains
for i, part in enumerate(my_model.mesh.elements):
    mat = feb.material.BiphasicMaterial(
        name=part.name,
        id=i + 1,
        solid=feb.material.NeoHookean(id=i + 1),
        permeability=feb.material.ConstantIsoPerm(
            perm=feb.material.MaterialParameter(text=1e-3 * (i + 1))
        ),
    )
    my_model.material.add_material(mat)
    my_model.mesh_domains.add_solid_domain(
        feb.meshdomains.SolidDomain(name=part.name, mat=part.name)
    )

# fix the bottom nodes in space
fix_bottom = my_model.boundary.add_bc(
    feb.boundary.BCZeroDisplacement(node_set="bottom", x_dof=1, y_dof=1, z_dof=1)
)

# fix the bottom nodes in fluid pressure
fix_bottom_fluid = my_model.boundary.add_bc(feb.boundary.BCZeroFluidPressure(node_set="bottom"))

# fix the bottom nodes in fluid flux
fix_bottom_flux = my_model.boundary.add_bc(feb.boundary.BCZeroFluidPressure(node_set="bottom"))

# set zero fluid pressure bc on top nodes to allow for free-draining
drain_top = my_model.boundary.add_bc(feb.boundary.BCZeroFluidPressure(node_set="top"))

# displace the top nodes by -0.5 mm in z
move_top = my_model.boundary.add_bc(
    feb.boundary.BCPrescribedDisplacement(
        node_set="top", dof="z", value=feb.boundary.Value(lc=1, text=-0.5)
    )
)

# load curve to apply displacement
my_model.load_data.add_load_curve(
    feb.loaddata.LoadCurve(
        id=1, points=feb.loaddata.CurvePoints(points=["0.0,0.0", "0.1,1.0", "1.0,1.0"])
    )
)

# add some additonial variables to plotfile
my_model.output.add_plotfile(
    feb.output.OutputPlotfile(
        all_vars=[
            feb.output.Var(type="displacement"),
            feb.output.Var(type="effective fluid pressure"),
            feb.output.Var(type="nodal fluid flux"),
        ]
    )
)

# save the model to disk
my_model.save("biphasic_hex20.feb")
# run the model
feb.model.run_model("biphasic_hex20.feb")

Sliding Contact

import meshio

import pyfebio as feb

# read a 27 node hex mesh in gmsh format
from_gmsh = meshio.gmsh.read("../../assets/gmsh/hex27_contact.msh")
# translate gmsh meshio object to an febio mesh
mesh = feb.mesh.translate_meshio(from_gmsh)

# initialize and febio model with default settings
my_model = feb.model.Model(mesh=mesh)

# loop over the mesh Elements (parts) and assign materials
# and solid domains
for i, part in enumerate(my_model.mesh.elements):
    mat = feb.material.NeoHookean(
        name=part.name,
        id=i + 1,
        E=feb.material.MaterialParameter(text=1.0),
        v=feb.material.MaterialParameter(text=0.3),
    )
    my_model.material.add_material(mat)
    my_model.mesh_domains.add_solid_domain(
        feb.meshdomains.SolidDomain(name=part.name, mat=part.name)
    )

# add a surface pair for the contact definition
my_model.mesh.add_surface_pair(
    feb.mesh.SurfacePair(name="contact", primary="bottom-box-top", secondary="top-box-bottom")
)

# fix the bottom nodes of the bottom box in space
my_model.boundary.add_bc(
    feb.boundary.BCZeroDisplacement(node_set="bottom-box-bottom", x_dof=1, y_dof=1, z_dof=1)
)

# fix the top nodes of the top box in the y dimension
my_model.boundary.add_bc(
    feb.boundary.BCZeroDisplacement(node_set="top-box-top", x_dof=0, y_dof=1, z_dof=0)
)

# move the top nodes of the top box in the the dimension
my_model.boundary.add_bc(
    feb.boundary.BCPrescribedDisplacement(
        node_set="top-box-top", dof="z", value=feb.boundary.Value(lc=1, text=-0.15)
    )
)

# move the top nodes of the top box in the x dimension
my_model.boundary.add_bc(
    feb.boundary.BCPrescribedDisplacement(
        node_set="top-box-top", dof="x", value=feb.boundary.Value(lc=1, text=0.3)
    )
)
# add a sliding contact definition. enforce it with augmented lagrangian multiplier
my_model.contact.add_contact(
    feb.contact.SlidingElastic(
        name="sliding_contact", surface_pair="contact", auto_penalty=1, laugon="AUGLAG", two_pass=1
    )
)

# load curve controlling the top nodes displacement
my_model.load_data.add_load_curve(
    feb.loaddata.LoadCurve(id=1, points=feb.loaddata.CurvePoints(points=["0,0", "1,1"]))
)

# add variables to the plotfile output
my_model.output.add_plotfile(
    feb.output.OutputPlotfile(
        all_vars=[
            feb.output.Var(type="displacement"),
            feb.output.Var(type="contact pressure"),
            feb.output.Var(type="contact gap"),
        ]
    )
)
# save the model to disk
my_model.save("contact.feb")
# run the model
feb.model.run_model("contact.feb")

Adaptive Remeshing

import meshio

import pyfebio as feb

# import an 8 node hex mesh in gmsh format
from_gmsh = meshio.gmsh.read("../../assets/gmsh/hex8.msh")

# convert the meshio object to febio mesh
mesh = feb.mesh.translate_meshio(from_gmsh)

# crete a base febio model but override the control section to have constant time step
my_model = feb.model.Model(
    mesh=mesh, control=feb.control.Control(time_steps=10, step_size=0.1, time_stepper=None)
)

# loop over all Elements in the mesh and assign materials and solid domains
for i, part in enumerate(my_model.mesh.elements):
    mat = feb.material.NeoHookean(
        name=part.name,
        id=i + 1,
        E=feb.material.MaterialParameter(text=10.0 * (i + 1)),
        v=feb.material.MaterialParameter(text=0.3),
    )
    my_model.material.add_material(mat)
    my_model.mesh_domains.add_solid_domain(
        feb.meshdomains.SolidDomain(name=part.name, mat=part.name)
    )

# fix the bottom nodes
fix_bottom = my_model.boundary.add_bc(
    feb.boundary.BCZeroDisplacement(node_set="bottom", x_dof=1, y_dof=1, z_dof=1)
)

# displace the top nodes in z by -0.5 mm
move_top = my_model.boundary.add_bc(
    feb.boundary.BCPrescribedDisplacement(
        node_set="top", dof="z", value=feb.boundary.Value(lc=1, text=0.5)
    )
)

# load curve controlling the top nodes displacement
my_model.load_data.add_load_curve(
    feb.loaddata.LoadCurve(id=1, points=feb.loaddata.CurvePoints(points=["0.0,0.0", "1.0,1.0"]))
)

# creat are mesh adaptor to refine the mesh based
# on the stress error criterion
adaptor = feb.meshadaptor.HexRefineAdaptor(
    elem_set="bottom-layer",
    max_iters=2,
    criterion=feb.meshadaptor.RelativeErrorCriterion(
        error=0.01, data=feb.meshadaptor.StressCriterion()
    ),
)
my_model.meshadaptor.add_adaptor(adaptor)

# we have to add the stress error output to the plotfile so our mesh
# adaptor works
my_model.output.add_plotfile(
    feb.output.OutputPlotfile(
        all_vars=[
            feb.output.Var(type="displacement"),
            feb.output.Var(type="stress error"),
        ]
    )
)

# save our madel to disk
my_model.save("mesh_adapt.feb")
# run the model
feb.model.run_model("mesh_adapt.feb")