{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Linear elasticity\n",
"\n",
"\n",
"\n",
"*Figure 1*: Linear elastically deformed 1mm $\\times$ 1mm Ferrite logo."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Introduction\n",
"\n",
"The classical first finite element problem to solve in solid mechanics is a linear balance\n",
"of momentum problem. We will use this to introduce a vector valued field, the displacements\n",
"$\\boldsymbol{u}(\\boldsymbol{x})$. In addition, some features of the\n",
"[`Tensors.jl`](https://github.com/Ferrite-FEM/Tensors.jl) toolbox are demonstrated.\n",
"\n",
"### Strong form\n",
"The strong form of the balance of momentum for quasi-static loading is given by\n",
"$$\n",
"\\begin{alignat*}{2}\n",
" \\mathrm{div}(\\boldsymbol{\\sigma}) + \\boldsymbol{b} &= 0 \\quad &&\\boldsymbol{x} \\in \\Omega \\\\\n",
" \\boldsymbol{u} &= \\boldsymbol{u}_\\mathrm{D} \\quad &&\\boldsymbol{x} \\in \\Gamma_\\mathrm{D} \\\\\n",
" \\boldsymbol{n} \\cdot \\boldsymbol{\\sigma} &= \\boldsymbol{t}_\\mathrm{N} \\quad &&\\boldsymbol{x} \\in \\Gamma_\\mathrm{N}\n",
"\\end{alignat*}\n",
"$$\n",
"where $\\boldsymbol{\\sigma}$ is the (Cauchy) stress tensor and $\\boldsymbol{b}$ the body force.\n",
"The domain, $\\Omega$, has the boundary $\\Gamma$, consisting of a Dirichlet part, $\\Gamma_\\mathrm{D}$, and\n",
"a Neumann part, $\\Gamma_\\mathrm{N}$, with outward pointing normal vector $\\boldsymbol{n}$.\n",
"$\\boldsymbol{u}_\\mathrm{D}$ denotes prescribed displacements on $\\Gamma_\\mathrm{D}$,\n",
"while $\\boldsymbol{t}_\\mathrm{N}$ the known tractions on $\\Gamma_\\mathrm{N}$.\n",
"\n",
"In this tutorial, we use linear elasticity, such that\n",
"$$\n",
"\\boldsymbol{\\sigma} = \\mathsf{C} : \\boldsymbol{\\varepsilon}, \\quad\n",
"\\boldsymbol{\\varepsilon} = \\left[\\mathrm{grad}(\\boldsymbol{u})\\right]^\\mathrm{sym}\n",
"$$\n",
"where $\\mathsf{C}$ is the 4th order elastic stiffness tensor and\n",
"$\\boldsymbol{\\varepsilon}$ the small strain tensor.\n",
"The colon, $:$, represents the double contraction,\n",
"$\\sigma_{ij} = \\mathsf{C}_{ijkl} \\varepsilon_{kl}$, and the superscript $\\mathrm{sym}$\n",
"denotes the symmetric part.\n",
"\n",
"### Weak form\n",
"The resulting weak form is given as follows: Find $\\boldsymbol{u} \\in \\mathbb{U}$ such that\n",
"$$\n",
"\\int_\\Omega\n",
" \\mathrm{grad}(\\delta \\boldsymbol{u}) : \\boldsymbol{\\sigma}\n",
"\\, \\mathrm{d}\\Omega\n",
"=\n",
"\\int_{\\Gamma}\n",
" \\delta \\boldsymbol{u} \\cdot \\boldsymbol{t}\n",
"\\, \\mathrm{d}\\Gamma\n",
"+\n",
"\\int_\\Omega\n",
" \\delta \\boldsymbol{u} \\cdot \\boldsymbol{b}\n",
"\\, \\mathrm{d}\\Omega\n",
"\\quad \\forall \\, \\delta \\boldsymbol{u} \\in \\mathbb{T},\n",
"$$\n",
"where $\\mathbb{U}$ and $\\mathbb{T}$ denote suitable trial and test function spaces.\n",
"$\\delta \\boldsymbol{u}$ is a vector valued test function and\n",
"$\\boldsymbol{t} = \\boldsymbol{\\sigma}\\cdot\\boldsymbol{n}$ is the traction vector on\n",
"the boundary. In this tutorial, we will neglect body forces (i.e. $\\boldsymbol{b} = \\boldsymbol{0}$) and the weak form reduces to\n",
"$$\n",
"\\int_\\Omega\n",
" \\mathrm{grad}(\\delta \\boldsymbol{u}) : \\boldsymbol{\\sigma}\n",
"\\, \\mathrm{d}\\Omega\n",
"=\n",
"\\int_{\\Gamma}\n",
" \\delta \\boldsymbol{u} \\cdot \\boldsymbol{t}\n",
"\\, \\mathrm{d}\\Gamma \\,.\n",
"$$\n",
"\n",
"### Finite element form\n",
"Finally, the finite element form is obtained by introducing the finite element shape functions.\n",
"Since the displacement field, $\\boldsymbol{u}$, is vector valued, we use vector valued shape functions\n",
"$\\delta\\boldsymbol{N}_i$ and $\\boldsymbol{N}_i$ to approximate the test and trial functions:\n",
"$$\n",
"\\boldsymbol{u} \\approx \\sum_{i=1}^N \\boldsymbol{N}_i (\\boldsymbol{x}) \\, \\hat{u}_i\n",
"\\qquad\n",
"\\delta \\boldsymbol{u} \\approx \\sum_{i=1}^N \\delta\\boldsymbol{N}_i (\\boldsymbol{x}) \\, \\delta \\hat{u}_i\n",
"$$\n",
"Here $N$ is the number of nodal variables, with $\\hat{u}_i$ and $\\delta\\hat{u}_i$ representing the $i$-th nodal value.\n",
"Using the Einstein summation convention, we can write this in short form as\n",
"$\\boldsymbol{u} \\approx \\boldsymbol{N}_i \\, \\hat{u}_i$ and $\\delta\\boldsymbol{u} \\approx \\delta\\boldsymbol{N}_i \\, \\delta\\hat{u}_i$.\n",
"\n",
"Inserting the these into the weak form, and noting that that the equation should hold for all $\\delta \\hat{u}_i$, we get\n",
"$$\n",
"\\underbrace{\\int_\\Omega \\mathrm{grad}(\\delta \\boldsymbol{N}_i) : \\boldsymbol{\\sigma}\\ \\mathrm{d}\\Omega}_{f_i^\\mathrm{int}} = \\underbrace{\\int_\\Gamma \\delta \\boldsymbol{N}_i \\cdot \\boldsymbol{t}\\ \\mathrm{d}\\Gamma}_{f_i^\\mathrm{ext}}\n",
"$$\n",
"Inserting the linear constitutive relationship, $\\boldsymbol{\\sigma} = \\mathsf{C}:\\boldsymbol{\\varepsilon}$,\n",
"in the internal force vector, $f_i^\\mathrm{int}$, yields the linear equation\n",
"$$\n",
"\\underbrace{\\left[\\int_\\Omega \\mathrm{grad}(\\delta \\boldsymbol{N}_i) : \\mathsf{C} : \\left[\\mathrm{grad}(\\boldsymbol{N}_j)\\right]^\\mathrm{sym}\\ \\mathrm{d}\\Omega\\right]}_{K_{ij}}\\ \\hat{u}_j = f_i^\\mathrm{ext}\n",
"$$\n",
"\n",
"## Implementation\n",
"First we load Ferrite, and some other packages we need."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Ferrite, FerriteGmsh, SparseArrays"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"As in the heat equation tutorial, we will use a unit square - but here we'll load the grid of the Ferrite logo!\n",
"This is done by downloading [`logo.geo`](logo.geo) and loading it using [FerriteGmsh.jl](https://github.com/Ferrite-FEM/FerriteGmsh.jl),"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Downloads: download\n",
"logo_mesh = \"logo.geo\"\n",
"asset_url = \"https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/\"\n",
"isfile(logo_mesh) || download(string(asset_url, logo_mesh), logo_mesh)\n",
"\n",
"FerriteGmsh.Gmsh.initialize() #hide\n",
"FerriteGmsh.Gmsh.gmsh.option.set_number(\"General.Verbosity\", 2) #hide\n",
"grid = togrid(logo_mesh);\n",
"FerriteGmsh.Gmsh.finalize(); #hide"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"The generated grid lacks the facetsets for the boundaries, so we add them by using Ferrite's\n",
"`addfacetset!`. It allows us to add facetsets to the grid based on coordinates.\n",
"Note that approximate comparison to 0.0 doesn't work well, so we use a tolerance instead."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"addfacetset!(grid, \"top\", x -> x[2] ≈ 1.0) # facets for which x[2] ≈ 1.0 for all nodes\n",
"addfacetset!(grid, \"left\", x -> abs(x[1]) < 1.0e-6)\n",
"addfacetset!(grid, \"bottom\", x -> abs(x[2]) < 1.0e-6);"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"### Trial and test functions\n",
"In this tutorial, we use the same linear Lagrange shape functions to approximate both the\n",
"test and trial spaces, i.e. $\\delta\\boldsymbol{N}_i = \\boldsymbol{N}_i$.\n",
"As our grid is composed of triangular elements, we need the Lagrange functions defined\n",
"on a `RefTriangle`. All currently available interpolations can be found under\n",
"`Interpolation`.\n",
"\n",
"Here we use linear triangular elements (also called constant strain triangles).\n",
"The vector valued shape functions are constructed by raising the interpolation\n",
"to the power `dim` (the dimension) since the displacement field has one component in each\n",
"spatial dimension."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"dim = 2\n",
"order = 1 # linear interpolation\n",
"ip = Lagrange{RefTriangle, order}()^dim; # vector valued interpolation"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"In order to evaluate the integrals, we need to specify the quadrature rules to use. Due to the\n",
"linear interpolation, a single quadrature point suffices, both inside the cell and on the facet.\n",
"In 2d, a facet is the edge of the element."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"qr = QuadratureRule{RefTriangle}(1) # 1 quadrature point\n",
"qr_face = FacetQuadratureRule{RefTriangle}(1);"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"Finally, we collect the interpolations and quadrature rules into the `CellValues` and `FacetValues`\n",
"buffers, which we will later use to evaluate the integrals over the cells and facets."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"cellvalues = CellValues(qr, ip)\n",
"facetvalues = FacetValues(qr_face, ip);"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"### Degrees of freedom\n",
"For distributing degrees of freedom, we define a `DofHandler`. The `DofHandler` knows that\n",
"`u` has two degrees of freedom per node because we vectorized the interpolation above."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"dh = DofHandler(grid)\n",
"add!(dh, :u, ip)\n",
"close!(dh);"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"> **Numbering of degrees of freedom**\n",
">\n",
"> A common assumption is that the numbering of degrees of freedom follows the global\n",
"> numbering of the nodes in the grid. This is *NOT* the case in Ferrite. For more\n",
"> details, see the Ferrite numbering rules."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Boundary conditions\n",
"We set Dirichlet boundary conditions by fixing the motion normal to the bottom and left\n",
"boundaries. The last argument to `Dirichlet` determines which components of the field should be\n",
"constrained. If no argument is given, all components are constrained by default."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ch = ConstraintHandler(dh)\n",
"add!(ch, Dirichlet(:u, getfacetset(grid, \"bottom\"), (x, t) -> 0.0, 2))\n",
"add!(ch, Dirichlet(:u, getfacetset(grid, \"left\"), (x, t) -> 0.0, 1))\n",
"close!(ch);"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"In addition, we will use Neumann boundary conditions on the top surface, where\n",
"we add a traction vector of the form\n",
"$$\n",
"\\boldsymbol{t}_\\mathrm{N}(\\boldsymbol{x}) = (20e3) x_1 \\boldsymbol{e}_2\\ \\mathrm{N}/\\mathrm{mm}^2\n",
"$$"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"traction(x) = Vec(0.0, 20.0e3 * x[1]);"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"On the right boundary, we don't do anything, resulting in a zero traction Neumann boundary.\n",
"In order to assemble the external forces, $f_i^\\mathrm{ext}$, we need to iterate over all\n",
"facets in the relevant facetset. We do this by using the `FacetIterator`."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "assemble_external_forces! (generic function with 1 method)"
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"function assemble_external_forces!(f_ext, dh, facetset, facetvalues, prescribed_traction)\n",
" # Create a temporary array for the facet's local contributions to the external force vector\n",
" fe_ext = zeros(getnbasefunctions(facetvalues))\n",
" for facet in FacetIterator(dh, facetset)\n",
" # Update the facetvalues to the correct facet number\n",
" reinit!(facetvalues, facet)\n",
" # Reset the temporary array for the next facet\n",
" fill!(fe_ext, 0.0)\n",
" # Access the cell's coordinates\n",
" cell_coordinates = getcoordinates(facet)\n",
" for qp in 1:getnquadpoints(facetvalues)\n",
" # Calculate the global coordinate of the quadrature point.\n",
" x = spatial_coordinate(facetvalues, qp, cell_coordinates)\n",
" tₚ = prescribed_traction(x)\n",
" # Get the integration weight for the current quadrature point.\n",
" dΓ = getdetJdV(facetvalues, qp)\n",
" for i in 1:getnbasefunctions(facetvalues)\n",
" Nᵢ = shape_value(facetvalues, qp, i)\n",
" fe_ext[i] += tₚ ⋅ Nᵢ * dΓ\n",
" end\n",
" end\n",
" # Add the local contributions to the correct indices in the global external force vector\n",
" assemble!(f_ext, celldofs(facet), fe_ext)\n",
" end\n",
" return f_ext\n",
"end"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"### Material behavior\n",
"Next, we need to define the material behavior, specifically the elastic stiffness tensor, $\\mathsf{C}$.\n",
"In this tutorial, we use plane strain linear isotropic elasticity, with Hooke's law as\n",
"$$\n",
"\\boldsymbol{\\sigma} = 2G \\boldsymbol{\\varepsilon}^\\mathrm{dev} + 3K \\boldsymbol{\\varepsilon}^\\mathrm{vol}\n",
"$$\n",
"where $G$ is the shear modulus and $K$ the bulk modulus. This expression can be written as\n",
"$\\boldsymbol{\\sigma} = \\mathsf{C}:\\boldsymbol{\\varepsilon}$, with\n",
"$$\n",
" \\mathsf{C} := \\frac{\\partial \\boldsymbol{\\sigma}}{\\partial \\boldsymbol{\\varepsilon}}\n",
"$$\n",
"The volumetric, $\\boldsymbol{\\varepsilon}^\\mathrm{vol}$,\n",
"and deviatoric, $\\boldsymbol{\\varepsilon}^\\mathrm{dev}$ strains, are defined as\n",
"$$\n",
"\\begin{align*}\n",
"\\boldsymbol{\\varepsilon}^\\mathrm{vol} &= \\frac{\\mathrm{tr}(\\boldsymbol{\\varepsilon})}{3}\\boldsymbol{I}, \\quad\n",
"\\boldsymbol{\\varepsilon}^\\mathrm{dev} &= \\boldsymbol{\\varepsilon} - \\boldsymbol{\\varepsilon}^\\mathrm{vol}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"Starting from Young's modulus, $E$, and Poisson's ratio, $\\nu$, the shear and bulk modulus are\n",
"$$\n",
"G = \\frac{E}{2(1 + \\nu)}, \\quad K = \\frac{E}{3(1 - 2\\nu)}\n",
"$$"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "166666.66666666663"
},
"metadata": {},
"execution_count": 11
}
],
"cell_type": "code",
"source": [
"Emod = 200.0e3 # Young's modulus [MPa]\n",
"ν = 0.3 # Poisson's ratio [-]\n",
"\n",
"Gmod = Emod / (2(1 + ν)) # Shear modulus\n",
"Kmod = Emod / (3(1 - 2ν)) # Bulk modulus"
],
"metadata": {},
"execution_count": 11
},
{
"cell_type": "markdown",
"source": [
"Finally, we demonstrate `Tensors.jl`'s automatic differentiation capabilities when\n",
"calculating the elastic stiffness tensor"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"C = gradient(ϵ -> 2 * Gmod * dev(ϵ) + 3 * Kmod * vol(ϵ), zero(SymmetricTensor{2, 2}));"
],
"metadata": {},
"execution_count": 12
},
{
"cell_type": "markdown",
"source": [
"### Element routine\n",
"To calculate the global stiffness matrix, $K_{ij}$, the element routine computes the\n",
"local stiffness matrix `ke` for a single element and assembles it into the global matrix.\n",
"`ke` is pre-allocated and reused for all elements.\n",
"\n",
"Note that the elastic stiffness tensor $\\mathsf{C}$ is constant.\n",
"Thus is needs to be computed and once and can then be used for all integration points."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "assemble_cell! (generic function with 1 method)"
},
"metadata": {},
"execution_count": 13
}
],
"cell_type": "code",
"source": [
"function assemble_cell!(ke, cellvalues, C)\n",
" for q_point in 1:getnquadpoints(cellvalues)\n",
" # Get the integration weight for the quadrature point\n",
" dΩ = getdetJdV(cellvalues, q_point)\n",
" for i in 1:getnbasefunctions(cellvalues)\n",
" # Gradient of the test function\n",
" ∇Nᵢ = shape_gradient(cellvalues, q_point, i)\n",
" for j in 1:getnbasefunctions(cellvalues)\n",
" # Symmetric gradient of the trial function\n",
" ∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q_point, j)\n",
" ke[i, j] += (∇Nᵢ ⊡ C ⊡ ∇ˢʸᵐNⱼ) * dΩ\n",
" end\n",
" end\n",
" end\n",
" return ke\n",
"end"
],
"metadata": {},
"execution_count": 13
},
{
"cell_type": "markdown",
"source": [
"### Global assembly\n",
"We define the function `assemble_global` to loop over the elements and do the global\n",
"assembly. The function takes the preallocated sparse matrix `K`, our DofHandler `dh`, our\n",
"`cellvalues` and the elastic stiffness tensor `C` as input arguments and computes the\n",
"global stiffness matrix `K`."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "assemble_global! (generic function with 1 method)"
},
"metadata": {},
"execution_count": 14
}
],
"cell_type": "code",
"source": [
"function assemble_global!(K, dh, cellvalues, C)\n",
" # Allocate the element stiffness matrix\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
" ke = zeros(n_basefuncs, n_basefuncs)\n",
" # Create an assembler\n",
" assembler = start_assemble(K)\n",
" # Loop over all cells\n",
" for cell in CellIterator(dh)\n",
" # Update the shape function gradients based on the cell coordinates\n",
" reinit!(cellvalues, cell)\n",
" # Reset the element stiffness matrix\n",
" fill!(ke, 0.0)\n",
" # Compute element contribution\n",
" assemble_cell!(ke, cellvalues, C)\n",
" # Assemble ke into K\n",
" assemble!(assembler, celldofs(cell), ke)\n",
" end\n",
" return K\n",
"end"
],
"metadata": {},
"execution_count": 14
},
{
"cell_type": "markdown",
"source": [
"### Solution of the system\n",
"The last step is to solve the system. First we allocate the global stiffness matrix `K`\n",
"and assemble it."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"K = allocate_matrix(dh)\n",
"assemble_global!(K, dh, cellvalues, C);"
],
"metadata": {},
"execution_count": 15
},
{
"cell_type": "markdown",
"source": [
"Then we allocate and assemble the external force vector."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"f_ext = zeros(ndofs(dh))\n",
"assemble_external_forces!(f_ext, dh, getfacetset(grid, \"top\"), facetvalues, traction);"
],
"metadata": {},
"execution_count": 16
},
{
"cell_type": "markdown",
"source": [
"To account for the Dirichlet boundary conditions we use the `apply!` function.\n",
"This modifies elements in `K` and `f`, such that we can get the\n",
"correct solution vector `u` by using solving the linear equation system $K_{ij} \\hat{u}_j = f^\\mathrm{ext}_i$,"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"apply!(K, f_ext, ch)\n",
"u = K \\ f_ext;"
],
"metadata": {},
"execution_count": 17
},
{
"cell_type": "markdown",
"source": [
"> **Numbering of degrees of freedom**\n",
">\n",
"> Once again, recall that numbering of degrees of freedom does *NOT* follow the global\n",
"> numbering of the nodes in the grid. Specifically, `u[2 * i - 1]` and `u[2 * i]` are\n",
"> *NOT* the displacements at node `i`."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Postprocessing\n",
"In this case, we want to analyze the displacements, as well as the stress field.\n",
"We calculate the stress in each quadrature point, and then export it in two different\n",
"ways:\n",
"1) Constant in each cell (matching the approximation of constant strains in each element).\n",
" Note that a current limitation is that cell data for second order tensors must be exported\n",
" component-wise (see issue #768)\n",
"2) Interpolated using the linear lagrange ansatz functions via the `L2Projector`."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function calculate_stresses(grid, dh, cv, u, C)\n",
" qp_stresses = [\n",
" [zero(SymmetricTensor{2, 2}) for _ in 1:getnquadpoints(cv)]\n",
" for _ in 1:getncells(grid)\n",
" ]\n",
" avg_cell_stresses = tuple((zeros(getncells(grid)) for _ in 1:3)...)\n",
" for cell in CellIterator(dh)\n",
" reinit!(cv, cell)\n",
" cell_stresses = qp_stresses[cellid(cell)]\n",
" for q_point in 1:getnquadpoints(cv)\n",
" ε = function_symmetric_gradient(cv, q_point, u, celldofs(cell))\n",
" cell_stresses[q_point] = C ⊡ ε\n",
" end\n",
" σ_avg = sum(cell_stresses) / getnquadpoints(cv)\n",
" avg_cell_stresses[1][cellid(cell)] = σ_avg[1, 1]\n",
" avg_cell_stresses[2][cellid(cell)] = σ_avg[2, 2]\n",
" avg_cell_stresses[3][cellid(cell)] = σ_avg[1, 2]\n",
" end\n",
" return qp_stresses, avg_cell_stresses\n",
"end\n",
"\n",
"qp_stresses, avg_cell_stresses = calculate_stresses(grid, dh, cellvalues, u, C);"
],
"metadata": {},
"execution_count": 18
},
{
"cell_type": "markdown",
"source": [
"We now use the the L2Projector to project the stress-field onto the piecewise linear\n",
"finite element space that we used to solve the problem."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"proj = L2Projector(Lagrange{RefTriangle, 1}(), grid)\n",
"stress_field = project(proj, qp_stresses, qr);\n",
"\n",
"color_data = zeros(Int, getncells(grid)) #hide\n",
"colors = [ #hide\n",
" \"1\" => 1, \"5\" => 1, # purple #hide\n",
" \"2\" => 2, \"3\" => 2, # red #hide\n",
" \"4\" => 3, # blue #hide\n",
" \"6\" => 4, # green #hide\n",
"] #hide\n",
"for (key, color) in colors #hide\n",
" for i in getcellset(grid, key) #hide\n",
" color_data[i] = color #hide\n",
" end #hide\n",
"end #hide"
],
"metadata": {},
"execution_count": 19
},
{
"cell_type": "markdown",
"source": [
"To visualize the result we export to a VTK-file. Specifically, an unstructured\n",
"grid file, `.vtu`, is created, which can be viewed in e.g.\n",
"[ParaView](https://www.paraview.org/)."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "VTKGridFile for the closed file \"linear_elasticity.vtu\"."
},
"metadata": {},
"execution_count": 20
}
],
"cell_type": "code",
"source": [
"VTKGridFile(\"linear_elasticity\", dh) do vtk\n",
" write_solution(vtk, dh, u)\n",
" for (i, key) in enumerate((\"11\", \"22\", \"12\"))\n",
" write_cell_data(vtk, avg_cell_stresses[i], \"sigma_\" * key)\n",
" end\n",
" write_projection(vtk, proj, stress_field, \"stress field\")\n",
" Ferrite.write_cellset(vtk, grid)\n",
" write_cell_data(vtk, color_data, \"colors\") #hide\n",
"end"
],
"metadata": {},
"execution_count": 20
},
{
"cell_type": "markdown",
"source": [
"We used the displacement field to visualize the deformed logo in *Figure 1*,\n",
"and in *Figure 2*, we demonstrate the difference between the interpolated stress\n",
"field and the constant stress in each cell."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"\n",
"\n",
"*Figure 2*: Vertical normal stresses (MPa) exported using the `L2Projector` (left)\n",
"and constant stress in each cell (right)."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Test #hide\n",
"linux_result = 0.31742879147646924 #hide\n",
"@test abs(norm(u) - linux_result) < 0.01 #hide\n",
"Sys.islinux() && @test norm(u) ≈ linux_result #hide\n",
"nothing #hide"
],
"metadata": {},
"execution_count": 21
},
{
"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
}