{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Example creating the stiffness for a linear elasticity element using tensors or matrices" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Ferrite\n", "using Tensors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stiffness" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Stiffness using normal matrices\n", "function ke_element_mat!{T, dim}(Ke, X::Vector{Vec{dim, T}}, fe_values::CellScalarValues{dim}, Ee, B, DB, BDB)\n", " n_basefuncs = getnbasefunctions(fe_values)\n", " @assert length(X) == n_basefuncs\n", " \n", " reinit!(fe_values, X)\n", " for q_point in 1:getnquadpoints(fe_values)\n", " for i in 1:n_basefuncs\n", " dNdx = shape_gradient(fe_values, q_point, i)[1]\n", " dNdy = shape_gradient(fe_values, q_point, i)[2]\n", " dNdz = shape_gradient(fe_values, q_point, i)[3]\n", "\n", " B[1, i * 3-2] = dNdx\n", " B[2, i * 3-1] = dNdy\n", " B[3, i * 3-0] = dNdz\n", " B[4, 3 * i-1] = dNdz\n", " B[4, 3 * i-0] = dNdy\n", " B[5, 3 * i-2] = dNdz\n", " B[5, 3 * i-0] = dNdx\n", " B[6, 3 * i-2] = dNdy\n", " B[6, 3 * i-1] = dNdx\n", " end\n", " \n", " A_mul_B!(DB, Ee, B)\n", " At_mul_B!(BDB, B, DB)\n", " scale!(BDB, getdetJdV(fe_values, q_point))\n", " for p in 1:size(Ke,1)\n", " for q in 1:size(Ke,2)\n", " Ke[p, q] += BDB[p, q]\n", " end\n", " end\n", " end\n", " \n", " return Ke\n", "end;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Stiffness using scalar values\n", "function ke_element!{T,dim}(Ke, X::Vector{Vec{dim, T}}, fe_values::CellScalarValues{dim}, C)\n", " n_basefuncs = getnbasefunctions(fe_values)\n", " @assert length(X) == n_basefuncs\n", " reinit!(fe_values, X)\n", " @inbounds for q_point in 1:getnquadpoints(fe_values)\n", " for a in 1:n_basefuncs\n", " for b in 1:n_basefuncs\n", " ∇ϕa = shape_gradient(fe_values, q_point, a)\n", " ∇ϕb = shape_gradient(fe_values, q_point, b)\n", " Ke_e = dotdot(∇ϕa, C, ∇ϕb) * getdetJdV(fe_values, q_point)\n", " for d1 in 1:dim, d2 in 1:dim\n", " Ke[dim*(a-1) + d1, dim*(b-1) + d2] += Ke_e[d1,d2]\n", " end\n", " end\n", " end\n", " end\n", " return Ke\n", "end;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Stiffness using vector values\n", "function ke_element2!{T,dim}(Ke, X::Vector{Vec{dim, T}}, fe_values::CellVectorValues{dim}, C)\n", " n_basefuncs = getnbasefunctions(fe_values)\n", " @assert length(X) * dim == n_basefuncs\n", " reinit!(fe_values, X)\n", " ɛ = [zero(SymmetricTensor{2, dim, T}) for i in 1:n_basefuncs]\n", " @inbounds for q_point in 1:getnquadpoints(fe_values)\n", " for i in 1:n_basefuncs\n", " ɛ[i] = symmetric(shape_gradient(fe_values, q_point, i)) \n", " end\n", " dΩ = getdetJdV(fe_values, q_point)\n", " for i in 1:n_basefuncs\n", " ɛC = ɛ[i] ⊡ C\n", " for j in 1:n_basefuncs\n", " Ke[i, j] += (ɛC ⊡ ɛ[j]) * dΩ\n", " end\n", " end\n", " end\n", " return Ke\n", "end;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "E = 200e9\n", "ν = 0.3\n", "λ = E*ν / ((1 + ν) * (1 - 2ν))\n", "μ = E / (2(1 + ν))\n", "δ(i,j) = i == j ? 1.0 : 0.0\n", "g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))\n", "\n", "C = SymmetricTensor{4, 3}(g)\n", "\n", "\n", "M = λ/ν * (1 - ν)\n", "\n", "Cmat = [ M λ λ 0.0 0.0 0.0;\n", " λ M λ 0.0 0.0 0.0;\n", " λ λ M 0.0 0.0 0.0;\n", " 0.0 0.0 0.0 μ 0.0 0.0;\n", " 0.0 0.0 0.0 0.0 μ 0.0;\n", " 0.0 0.0 0.0 0.0 0.0 μ]\n", "\n", "\n", "interpolation = Lagrange{3, RefCube, 1}()\n", "quad_rule = QuadratureRule{3, RefCube}(1)\n", "values = CellScalarValues(quad_rule, interpolation);\n", "vector_values = CellVectorValues(quad_rule, interpolation);\n", "\n", "# Generate some coordinates\n", "x = [-1.0 -1.0 -1.0;\n", " 1.0 -1.0 -1.0;\n", " 1.0 1.0 -1.0;\n", " -1.0 1.0 -1.0;\n", " -1.0 -1.0 1.0;\n", " 1.0 -1.0 1.0;\n", " 1.0 1.0 1.0;\n", " -1.0 1.0 1.0;]\n", "x = x .+ 0.05 * rand()\n", "x_vec = reinterpret(Vec{3, Float64}, x, (8,));\n", "\n", "n_basefunctions = getnbasefunctions(vector_values)\n", "Ke = zeros(n_basefunctions, n_basefunctions)\n", "Ke2 = copy(Ke)\n", "Ke3 = copy(Ke)\n", "\n", "B = zeros(6, n_basefunctions)\n", "DB = zeros(6, n_basefunctions)\n", "BDB = zeros(n_basefunctions, n_basefunctions);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fill!(Ke, 0)\n", "fill!(Ke2, 0)\n", "fill!(Ke3, 0)\n", "ke_element!(Ke2, x_vec, values, C)\n", "ke_element2!(Ke3, x_vec, vector_values, C);\n", "ke_element_mat!(Ke, x_vec, values, Cmat, B, DB, BDB);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[1m\u001b[32mTest Passed\n", "\u001b[39m\u001b[22m" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Test\n", "@test norm(Ke - Ke2) / norm(Ke) < 1e-14\n", "@test norm(Ke - Ke3) / norm(Ke) < 1e-14" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Stiffness successful\n" ] } ], "source": [ "println(\"Stiffness successful\")" ] } ], "metadata": { "kernelspec": { "name": "julia-0.6", "display_name": "Julia 0.6.0", "language": "julia" }, "language_info": { "mimetype": "application/julia", "file_extension": ".jl", "version": "0.6.0", "name": "julia" } }, "nbformat": 4, "nbformat_minor": 1 }