{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Incompressible elasticity"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Introduction\n",
"\n",
"Mixed elements can be used to overcome locking when the material becomes\n",
"incompressible. However, for an element to be stable, it needs to fulfill\n",
"the LBB condition.\n",
"In this example we will consider two different element formulations\n",
"- linear displacement with linear pressure approximation (does *not* fulfill LBB)\n",
"- quadratic displacement with linear pressure approximation (does fulfill LBB)\n",
"The quadratic/linear element is also known as the Taylor-Hood element.\n",
"We will consider Cook's Membrane with an applied traction on the right hand side."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Commented program\n",
"\n",
"What follows is a program spliced with comments."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Ferrite, Tensors\n",
"using BlockArrays, SparseArrays, LinearAlgebra"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"First we generate a simple grid, specifying the 4 corners of Cooks membrane."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function create_cook_grid(nx, ny)\n",
" corners = [\n",
" Vec{2}((0.0, 0.0)),\n",
" Vec{2}((48.0, 44.0)),\n",
" Vec{2}((48.0, 60.0)),\n",
" Vec{2}((0.0, 44.0)),\n",
" ]\n",
" grid = generate_grid(Triangle, (nx, ny), corners)\n",
" # facesets for boundary conditions\n",
" addfacetset!(grid, \"clamped\", x -> norm(x[1]) ≈ 0.0)\n",
" addfacetset!(grid, \"traction\", x -> norm(x[1]) ≈ 48.0)\n",
" return grid\n",
"end;"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"Next we define a function to set up our cell- and FacetValues."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function create_values(interpolation_u, interpolation_p)\n",
" # quadrature rules\n",
" qr = QuadratureRule{RefTriangle}(3)\n",
" facet_qr = FacetQuadratureRule{RefTriangle}(3)\n",
"\n",
" # cell and FacetValues for u\n",
" cellvalues_u = CellValues(qr, interpolation_u)\n",
" facetvalues_u = FacetValues(facet_qr, interpolation_u)\n",
"\n",
" # cellvalues for p\n",
" cellvalues_p = CellValues(qr, interpolation_p)\n",
"\n",
" return cellvalues_u, cellvalues_p, facetvalues_u\n",
"end;"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"We create a DofHandler, with two fields, `:u` and `:p`,\n",
"with possibly different interpolations"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function create_dofhandler(grid, ipu, ipp)\n",
" dh = DofHandler(grid)\n",
" add!(dh, :u, ipu) # displacement\n",
" add!(dh, :p, ipp) # pressure\n",
" close!(dh)\n",
" return dh\n",
"end;"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"We also need to add Dirichlet boundary conditions on the `\"clamped\"` facetset.\n",
"We specify a homogeneous Dirichlet bc on the displacement field, `:u`."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function create_bc(dh)\n",
" dbc = ConstraintHandler(dh)\n",
" add!(dbc, Dirichlet(:u, getfacetset(dh.grid, \"clamped\"), x -> zero(x), [1, 2]))\n",
" close!(dbc)\n",
" return dbc\n",
"end;"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"The material is linear elastic, which is here specified by the shear and bulk moduli"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"struct LinearElasticity{T}\n",
" G::T\n",
" K::T\n",
"end"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"Now to the assembling of the stiffness matrix. This mixed formulation leads to a blocked\n",
"element matrix. Since Ferrite does not force us to use any particular matrix type we will\n",
"use a `BlockedArray` from `BlockArrays.jl`."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function doassemble(\n",
" cellvalues_u::CellValues,\n",
" cellvalues_p::CellValues,\n",
" facetvalues_u::FacetValues,\n",
" K::SparseMatrixCSC, grid::Grid, dh::DofHandler, mp::LinearElasticity\n",
" )\n",
" f = zeros(ndofs(dh))\n",
" assembler = start_assemble(K, f)\n",
" nu = getnbasefunctions(cellvalues_u)\n",
" np = getnbasefunctions(cellvalues_p)\n",
"\n",
" fe = BlockedArray(zeros(nu + np), [nu, np]) # local force vector\n",
" ke = BlockedArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix\n",
"\n",
" # traction vector\n",
" t = Vec{2}((0.0, 1 / 16))\n",
" # cache ɛdev outside the element routine to avoid some unnecessary allocations\n",
" ɛdev = [zero(SymmetricTensor{2, 2}) for i in 1:getnbasefunctions(cellvalues_u)]\n",
"\n",
" for cell in CellIterator(dh)\n",
" fill!(ke, 0)\n",
" fill!(fe, 0)\n",
" assemble_up!(ke, fe, cell, cellvalues_u, cellvalues_p, facetvalues_u, grid, mp, ɛdev, t)\n",
" assemble!(assembler, celldofs(cell), ke, fe)\n",
" end\n",
"\n",
" return K, f\n",
"end;"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"The element routine integrates the local stiffness and force vector for all elements.\n",
"Since the problem results in a symmetric matrix we choose to only assemble the lower part,\n",
"and then symmetrize it after the loop over the quadrature points."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facetvalues_u, grid, mp, ɛdev, t)\n",
"\n",
" n_basefuncs_u = getnbasefunctions(cellvalues_u)\n",
" n_basefuncs_p = getnbasefunctions(cellvalues_p)\n",
" u▄, p▄ = 1, 2\n",
" reinit!(cellvalues_u, cell)\n",
" reinit!(cellvalues_p, cell)\n",
"\n",
" # We only assemble lower half triangle of the stiffness matrix and then symmetrize it.\n",
" for q_point in 1:getnquadpoints(cellvalues_u)\n",
" for i in 1:n_basefuncs_u\n",
" ɛdev[i] = dev(symmetric(shape_gradient(cellvalues_u, q_point, i)))\n",
" end\n",
" dΩ = getdetJdV(cellvalues_u, q_point)\n",
" for i in 1:n_basefuncs_u\n",
" divδu = shape_divergence(cellvalues_u, q_point, i)\n",
" δu = shape_value(cellvalues_u, q_point, i)\n",
" for j in 1:i\n",
" Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * mp.G * ɛdev[i] ⊡ ɛdev[j] * dΩ\n",
" end\n",
" end\n",
"\n",
" for i in 1:n_basefuncs_p\n",
" δp = shape_value(cellvalues_p, q_point, i)\n",
" for j in 1:n_basefuncs_u\n",
" divδu = shape_divergence(cellvalues_u, q_point, j)\n",
" Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ\n",
" end\n",
" for j in 1:i\n",
" p = shape_value(cellvalues_p, q_point, j)\n",
" Ke[BlockIndex((p▄, p▄), (i, j))] += - 1 / mp.K * δp * p * dΩ\n",
" end\n",
"\n",
" end\n",
" end\n",
"\n",
" symmetrize_lower!(Ke)\n",
"\n",
" # We integrate the Neumann boundary using the FacetValues.\n",
" # We loop over all the facets in the cell, then check if the facet\n",
" # is in our `\"traction\"` facetset.\n",
" for facet in 1:nfacets(cell)\n",
" if (cellid(cell), facet) ∈ getfacetset(grid, \"traction\")\n",
" reinit!(facetvalues_u, cell, facet)\n",
" for q_point in 1:getnquadpoints(facetvalues_u)\n",
" dΓ = getdetJdV(facetvalues_u, q_point)\n",
" for i in 1:n_basefuncs_u\n",
" δu = shape_value(facetvalues_u, q_point, i)\n",
" fe[i] += (δu ⋅ t) * dΓ\n",
" end\n",
" end\n",
" end\n",
" end\n",
" return\n",
"end\n",
"\n",
"function symmetrize_lower!(Ke)\n",
" for i in 1:size(Ke, 1)\n",
" for j in (i + 1):size(Ke, 1)\n",
" Ke[i, j] = Ke[j, i]\n",
" end\n",
" end\n",
" return\n",
"end;"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"To evaluate the stresses after solving the problem we once again loop over the cells in\n",
"the grid. Stresses are evaluated in the quadrature points, however, for\n",
"export/visualization you typically want values in the nodes of the mesh, or as single data\n",
"points per cell. For the former you can project the quadrature point data to a finite\n",
"element space (see the example with the `L2Projector` in Post processing and\n",
"visualization). In this example we choose to compute the mean\n",
"value of the stress within each cell, and thus end up with one data point per cell. The\n",
"mean value is computed as\n",
"$$\n",
"\\bar{\\boldsymbol{\\sigma}}_i = \\frac{1}{ |\\Omega_i|}\n",
"\\int_{\\Omega_i} \\boldsymbol{\\sigma}\\, \\mathrm{d}\\Omega, \\quad\n",
"|\\Omega_i| = \\int_{\\Omega_i} 1\\, \\mathrm{d}\\Omega\n",
"$$\n",
"where $\\Omega_i$ is the domain occupied by cell number $i$, and $|\\Omega_i|$ the volume\n",
"(area) of the cell. The integrals are evaluated using numerical quadrature with the help\n",
"of cellvalues for u and p, just like in the assembly procedure.\n",
"\n",
"Note that even though all strain components in the out-of-plane direction are zero (plane\n",
"strain) the stress components are not. Specifically, $\\sigma_{33}$ will be non-zero in\n",
"this formulation. Therefore we expand the strain to a 3D tensor, and then compute the (3D)\n",
"stress tensor."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function compute_stresses(\n",
" cellvalues_u::CellValues, cellvalues_p::CellValues,\n",
" dh::DofHandler, mp::LinearElasticity, a::Vector\n",
" )\n",
" ae = zeros(ndofs_per_cell(dh)) # local solution vector\n",
" u_range = dof_range(dh, :u) # local range of dofs corresponding to u\n",
" p_range = dof_range(dh, :p) # local range of dofs corresponding to p\n",
" # Allocate storage for the stresses\n",
" σ = zeros(SymmetricTensor{2, 3}, getncells(dh.grid))\n",
" # Loop over the cells and compute the cell-average stress\n",
" for cc in CellIterator(dh)\n",
" # Update cellvalues\n",
" reinit!(cellvalues_u, cc)\n",
" reinit!(cellvalues_p, cc)\n",
" # Extract the cell local part of the solution\n",
" for (i, I) in pairs(celldofs(cc))\n",
" ae[i] = a[I]\n",
" end\n",
" # Loop over the quadrature points\n",
" σΩi = zero(SymmetricTensor{2, 3}) # stress integrated over the cell\n",
" Ωi = 0.0 # cell volume (area)\n",
" for qp in 1:getnquadpoints(cellvalues_u)\n",
" dΩ = getdetJdV(cellvalues_u, qp)\n",
" # Evaluate the strain and the pressure\n",
" ε = function_symmetric_gradient(cellvalues_u, qp, ae, u_range)\n",
" p = function_value(cellvalues_p, qp, ae, p_range)\n",
" # Expand strain to 3D\n",
" ε3D = SymmetricTensor{2, 3}((i, j) -> i < 3 && j < 3 ? ε[i, j] : 0.0)\n",
" # Compute the stress in this quadrature point\n",
" σqp = 2 * mp.G * dev(ε3D) - one(ε3D) * p\n",
" σΩi += σqp * dΩ\n",
" Ωi += dΩ\n",
" end\n",
" # Store the value\n",
" σ[cellid(cc)] = σΩi / Ωi\n",
" end\n",
" return σ\n",
"end;"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"Now we have constructed all the necessary components, we just need a function\n",
"to put it all together."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "solve (generic function with 1 method)"
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"function solve(ν, interpolation_u, interpolation_p)\n",
" # material\n",
" Emod = 1.0\n",
" Gmod = Emod / 2(1 + ν)\n",
" Kmod = Emod * ν / ((1 + ν) * (1 - 2ν))\n",
" mp = LinearElasticity(Gmod, Kmod)\n",
"\n",
" # Grid, dofhandler, boundary condition\n",
" n = 50\n",
" grid = create_cook_grid(n, n)\n",
" dh = create_dofhandler(grid, interpolation_u, interpolation_p)\n",
" dbc = create_bc(dh)\n",
"\n",
" # CellValues\n",
" cellvalues_u, cellvalues_p, facetvalues_u = create_values(interpolation_u, interpolation_p)\n",
"\n",
" # Assembly and solve\n",
" K = allocate_matrix(dh)\n",
" K, f = doassemble(cellvalues_u, cellvalues_p, facetvalues_u, K, grid, dh, mp)\n",
" apply!(K, f, dbc)\n",
" u = K \\ f\n",
"\n",
" # Compute the stress\n",
" σ = compute_stresses(cellvalues_u, cellvalues_p, dh, mp, u)\n",
" σvM = map(x -> √(3 / 2 * dev(x) ⊡ dev(x)), σ) # von Mise effective stress\n",
"\n",
" # Export the solution and the stress\n",
" filename = \"cook_\" *\n",
" (interpolation_u == Lagrange{RefTriangle, 1}()^2 ? \"linear\" : \"quadratic\") *\n",
" \"_linear\"\n",
"\n",
" VTKGridFile(filename, grid) do vtk\n",
" write_solution(vtk, dh, u)\n",
" for i in 1:3, j in 1:3\n",
" σij = [x[i, j] for x in σ]\n",
" write_cell_data(vtk, σij, \"sigma_$(i)$(j)\")\n",
" end\n",
" write_cell_data(vtk, σvM, \"sigma von Mises\")\n",
" end\n",
" return u\n",
"end"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"We now define the interpolation for displacement and pressure. We use (scalar) Lagrange\n",
"interpolation as a basis for both, and for the displacement, which is a vector, we\n",
"vectorize it to 2 dimensions such that we obtain vector shape functions (and 2nd order\n",
"tensors for the gradients)."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Lagrange{RefTriangle, 2}()^2"
},
"metadata": {},
"execution_count": 11
}
],
"cell_type": "code",
"source": [
"linear_p = Lagrange{RefTriangle, 1}()\n",
"linear_u = Lagrange{RefTriangle, 1}()^2\n",
"quadratic_u = Lagrange{RefTriangle, 2}()^2"
],
"metadata": {},
"execution_count": 11
},
{
"cell_type": "markdown",
"source": [
"All that is left is to solve the problem. We choose a value of Poissons\n",
"ratio that results in incompressibility ($ν = 0.5$) and thus expect the\n",
"linear/linear approximation to return garbage, and the quadratic/linear\n",
"approximation to be stable."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"u1 = solve(0.5, linear_u, linear_p);\n",
"u2 = solve(0.5, quadratic_u, linear_p);"
],
"metadata": {},
"execution_count": 12
},
{
"cell_type": "markdown",
"source": [
"---\n",
"\n",
"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.3"
},
"kernelspec": {
"name": "julia-1.11",
"display_name": "Julia 1.11.3",
"language": "julia"
}
},
"nbformat": 4
}