{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Hyperelasticity\n", "Author: Jørgen S. Dokken and Garth N. Wells\n", "\n", "This section shows how to solve the hyperelasticity problem for deformation of a beam.\n", "\n", "We will also show how to create a constant boundary condition for a vector function space.\n", "\n", "We start by importing DOLFINx and some additional dependencies.\n", "Then, we create a slender cantilever consisting of hexahedral elements and create the function space `V` for our unknown." ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": { "tags": [] }, "outputs": [], "source": [ "from dolfinx import log, default_scalar_type\n", "from dolfinx.fem.petsc import NonlinearProblem\n", "from dolfinx.nls.petsc import NewtonSolver\n", "import pyvista\n", "import numpy as np\n", "import ufl\n", "\n", "from mpi4py import MPI\n", "from dolfinx import fem, mesh, plot\n", "L = 20.0\n", "domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron)\n", "V = fem.functionspace(domain, (\"Lagrange\", 2, (domain.geometry.dim, )))" ] }, { "cell_type": "markdown", "id": "2", "metadata": {}, "source": [ "We create two python functions for determining the facets to apply boundary conditions to" ] }, { "cell_type": "code", "execution_count": null, "id": "3", "metadata": {}, "outputs": [], "source": [ "def left(x):\n", " return np.isclose(x[0], 0)\n", "\n", "\n", "def right(x):\n", " return np.isclose(x[0], L)\n", "\n", "\n", "fdim = domain.topology.dim - 1\n", "left_facets = mesh.locate_entities_boundary(domain, fdim, left)\n", "right_facets = mesh.locate_entities_boundary(domain, fdim, right)" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "Next, we create a marker based on these two functions" ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [], "source": [ "# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two\n", "marked_facets = np.hstack([left_facets, right_facets])\n", "marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])\n", "sorted_facets = np.argsort(marked_facets)\n", "facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "We then create a function for supplying the boundary condition on the left side, which is fixed." ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [], "source": [ "u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)" ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "To apply the boundary condition, we identity the dofs located on the facets marked by the `MeshTag`." ] }, { "cell_type": "code", "execution_count": null, "id": "9", "metadata": {}, "outputs": [], "source": [ "left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))\n", "bcs = [fem.dirichletbc(u_bc, left_dofs, V)]" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "Next, we define the body force on the reference configuration (`B`), and nominal (first Piola-Kirchhoff) traction (`T`)." ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [], "source": [ "B = fem.Constant(domain, default_scalar_type((0, 0, 0)))\n", "T = fem.Constant(domain, default_scalar_type((0, 0, 0)))" ] }, { "cell_type": "markdown", "id": "12", "metadata": {}, "source": [ "Define the test and solution functions on the space $V$" ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "v = ufl.TestFunction(V)\n", "u = fem.Function(V)" ] }, { "cell_type": "markdown", "id": "14", "metadata": {}, "source": [ "Define kinematic quantities used in the problem" ] }, { "cell_type": "code", "execution_count": null, "id": "15", "metadata": {}, "outputs": [], "source": [ "# Spatial dimension\n", "d = len(u)\n", "\n", "# Identity tensor\n", "I = ufl.variable(ufl.Identity(d))\n", "\n", "# Deformation gradient\n", "F = ufl.variable(I + ufl.grad(u))\n", "\n", "# Right Cauchy-Green tensor\n", "C = ufl.variable(F.T * F)\n", "\n", "# Invariants of deformation tensors\n", "Ic = ufl.variable(ufl.tr(C))\n", "J = ufl.variable(ufl.det(F))" ] }, { "cell_type": "markdown", "id": "16", "metadata": {}, "source": [ "Define the elasticity model via a stored strain energy density function $\\psi$, and create the expression for the first Piola-Kirchhoff stress:" ] }, { "cell_type": "code", "execution_count": null, "id": "17", "metadata": {}, "outputs": [], "source": [ "# Elasticity parameters\n", "E = default_scalar_type(1.0e4)\n", "nu = default_scalar_type(0.3)\n", "mu = fem.Constant(domain, E / (2 * (1 + nu)))\n", "lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))\n", "# Stored strain energy density (compressible neo-Hookean model)\n", "psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2\n", "# Stress\n", "# Hyper-elasticity\n", "P = ufl.diff(psi, F)" ] }, { "cell_type": "markdown", "id": "18", "metadata": {}, "source": [ "```{admonition} Comparison to linear elasticity\n", "To illustrate the difference between linear and hyperelasticity, the following lines can be uncommented to solve the linear elasticity problem.\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "19", "metadata": {}, "outputs": [], "source": [ "# P = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I" ] }, { "cell_type": "markdown", "id": "20", "metadata": {}, "source": [ "Define the variational form with traction integral over all facets with value 2. We set the quadrature degree for the integrals to 4." ] }, { "cell_type": "code", "execution_count": null, "id": "21", "metadata": {}, "outputs": [], "source": [ "metadata = {\"quadrature_degree\": 4}\n", "ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)\n", "dx = ufl.Measure(\"dx\", domain=domain, metadata=metadata)\n", "# Define form F (we want to find u such that F(u) = 0)\n", "F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)" ] }, { "cell_type": "markdown", "id": "22", "metadata": {}, "source": [ "As the varitional form is non-linear and written on residual form, we use the non-linear problem class from DOLFINx to set up required structures to use a Newton solver." ] }, { "cell_type": "code", "execution_count": null, "id": "23", "metadata": {}, "outputs": [], "source": [ "problem = NonlinearProblem(F, u, bcs)" ] }, { "cell_type": "markdown", "id": "24", "metadata": {}, "source": [ "and then create and customize the Newton solver" ] }, { "cell_type": "code", "execution_count": null, "id": "25", "metadata": {}, "outputs": [], "source": [ "solver = NewtonSolver(domain.comm, problem)\n", "\n", "# Set Newton solver options\n", "solver.atol = 1e-8\n", "solver.rtol = 1e-8\n", "solver.convergence_criterion = \"incremental\"\n" ] }, { "cell_type": "markdown", "id": "26", "metadata": {}, "source": [ "We create a function to plot the solution at each time step." ] }, { "cell_type": "code", "execution_count": null, "id": "27", "metadata": {}, "outputs": [], "source": [ "pyvista.start_xvfb()\n", "plotter = pyvista.Plotter()\n", "plotter.open_gif(\"deformation.gif\", fps=3)\n", "\n", "topology, cells, geometry = plot.vtk_mesh(u.function_space)\n", "function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)\n", "\n", "values = np.zeros((geometry.shape[0], 3))\n", "values[:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))\n", "function_grid[\"u\"] = values\n", "function_grid.set_active_vectors(\"u\")\n", "\n", "# Warp mesh by deformation\n", "warped = function_grid.warp_by_vector(\"u\", factor=1)\n", "warped.set_active_vectors(\"u\")\n", "\n", "# Add mesh to plotter and visualize\n", "actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10])\n", "\n", "# Compute magnitude of displacement to visualize in GIF\n", "Vs = fem.functionspace(domain, (\"Lagrange\", 2))\n", "magnitude = fem.Function(Vs)\n", "us = fem.Expression(ufl.sqrt(sum([u[i]**2 for i in range(len(u))])), Vs.element.interpolation_points)\n", "magnitude.interpolate(us)\n", "warped[\"mag\"] = magnitude.x.array" ] }, { "cell_type": "markdown", "id": "28", "metadata": {}, "source": [ "Finally, we solve the problem over several time steps, updating the z-component of the traction" ] }, { "cell_type": "code", "execution_count": null, "id": "29", "metadata": {}, "outputs": [], "source": [ "log.set_log_level(log.LogLevel.INFO)\n", "tval0 = -1.5\n", "for n in range(1, 10):\n", " T.value[2] = n * tval0\n", " num_its, converged = solver.solve(u)\n", " assert (converged)\n", " u.x.scatter_forward()\n", " print(f\"Time step {n}, Number of iterations {num_its}, Load {T.value}\")\n", " function_grid[\"u\"][:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))\n", " magnitude.interpolate(us)\n", " warped.set_active_scalars(\"mag\")\n", " warped_n = function_grid.warp_by_vector(factor=1)\n", " warped.points[:, :] = warped_n.points\n", " warped.point_data[\"mag\"][:] = magnitude.x.array\n", " plotter.update_scalar_bar_range([0, 10])\n", " plotter.write_frame()\n", "plotter.close()" ] }, { "cell_type": "markdown", "id": "30", "metadata": {}, "source": [ "\"gif\"" ] } ], "metadata": { "jupytext": { "formats": "ipynb,py:light" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }