diff --git a/docs/tutorials/2015-08-29-developing-juliafem.ipynb b/docs/tutorials/2015-08-29-developing-juliafem.ipynb index c2912ab..56f9232 100644 --- a/docs/tutorials/2015-08-29-developing-juliafem.ipynb +++ b/docs/tutorials/2015-08-29-developing-juliafem.ipynb @@ -49,7 +49,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": { "collapsed": false }, @@ -60,7 +60,7 @@ "Logger(root,DEBUG,Base.PipeEndpoint(open, 0 bytes waiting),root)" ] }, - "execution_count": 2, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -116,7 +116,7 @@ "type MyQuad4 <: Element\n", " connectivity :: Array{Int, 1}\n", " basis :: Basis\n", - " fields :: Dict{Symbol, FieldSet}\n", + " fields :: Dict{ASCIIString, FieldSet}\n", "end" ] }, @@ -147,16 +147,10 @@ ], "source": [ "function MyQuad4(connectivity)\n", - " h(xi) = [\n", - " (1-xi[1])*(1-xi[2])/4\n", - " (1+xi[1])*(1-xi[2])/4\n", - " (1+xi[1])*(1+xi[2])/4\n", - " (1-xi[1])*(1+xi[2])/4]\n", - " dh(xi) = [\n", - " -(1-xi[2])/4.0 -(1-xi[1])/4.0\n", - " (1-xi[2])/4.0 -(1+xi[1])/4.0\n", - " (1+xi[2])/4.0 (1+xi[1])/4.0\n", - " -(1+xi[2])/4.0 (1-xi[1])/4.0]\n", + " h(xi) = 1/4*[(1-xi[1])*(1-xi[2]) (1+xi[1])*(1-xi[2]) (1+xi[1])*(1+xi[2]) (1-xi[1])*(1+xi[2])]\n", + " dh(xi) = 1/4*[\n", + " -(1-xi[2]) (1-xi[2]) (1+xi[2]) -(1+xi[2])\n", + " -(1-xi[1]) -(1+xi[1]) (1+xi[1]) (1-xi[1])]\n", " basis = Basis(h, dh)\n", " MyQuad4(connectivity, basis, Dict())\n", "end" @@ -166,7 +160,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Also some basic information like number of basis functions and element dimension is needed." + "Also some basic information like number of basis functions and element dimension is needed. We can set them directly to the `Base.size`:" ] }, { @@ -179,7 +173,7 @@ { "data": { "text/plain": [ - "get_element_dimension (generic function with 7 methods)" + "size (generic function with 72 methods)" ] }, "execution_count": 7, @@ -188,15 +182,14 @@ } ], "source": [ - "JuliaFEM.get_number_of_basis_functions(el::Type{MyQuad4}) = 4\n", - "JuliaFEM.get_element_dimension(el::Type{MyQuad4}) = 2" + "Base.size(element::Type{MyQuad4}) = (2, 4) # 2 => (ξ₁, ξ₂), 4 => 4 basis functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Next we check that everything is well defined:" + "Here comes one important thing. We always define our \"things\" so that the first index is dimension, like $(x, y, z)$ or $(\\xi_1, \\xi_2, \\xi_3)$ and second index is basis function number / node id or something similar to that. To motivate this, consider the following example:" ] }, { @@ -205,22 +198,143 @@ "metadata": { "collapsed": false }, + "outputs": [ + { + "data": { + "text/plain": [ + "2x4 Array{Int64,2}:\n", + " 1 3 5 7\n", + " 2 4 6 8" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = [1 2; 3 4; 5 6; 7 8]'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, if we consider $x$ as of some field e.g. geometry, our coordinates of first node is $(1, 2)$, second is $(3, 4)$ and so on. Typically on vector field problems the global assembly is something like $(u_1, v_1, u_2, v_2, \\ldots, )$. If fields are defined this way we can easily flatten matrix to vector and back:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8-element Array{Int64,1}:\n", + " 1\n", + " 2\n", + " 3\n", + " 4\n", + " 5\n", + " 6\n", + " 7\n", + " 8" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dim = size(x)\n", + "x2 = vec(x)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2x4 Array{Int64,2}:\n", + " 1 3 5 7\n", + " 2 4 6 8" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "reshape(x2, dim)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Also:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 2 3 4 5 6 7 8 " + ] + } + ], + "source": [ + "for i=1:8\n", + " print(\"$(x[i]) \")\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So when keeping this order, things automatically go right in Julia.\n", + "\n", + "Next we check that everything is well defined:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "26-Oct 05:30:31:INFO:root:Testing element MyQuad4\n", - "26-Oct 05:30:31:INFO:root:number of basis functions in this element: 4\n", - "26-Oct 05:30:31:INFO:root:Initializing element\n", - "26-Oct 05:30:31:INFO:root:Element dimension: 2\n", - "26-Oct 05:30:31:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])\n", - "26-Oct 05:30:31:INFO:root:basis at [0.0,0.0]: [0.25 0.25 0.25 0.25]\n", - "26-Oct 05:30:31:INFO:root:field val at [0.0,0.0]: 2.5\n", - "26-Oct 05:30:32:INFO:root:derivative of basis at [0.0,0.0]: [-0.5 0.5 0.5 -0.5\n", + "27-Oct 06:20:33:INFO:root:Testing element MyQuad4\n", + "27-Oct 06:20:33:INFO:root:element dimension: 2 x 4\n", + "27-Oct 06:20:33:INFO:root:Initializing element\n", + "27-Oct 06:20:34:INFO:root:basis at [0.0,0.0]: [0.25 0.25 0.25 0.25]\n", + "27-Oct 06:20:34:INFO:root:field val at [0.0,0.0]: 2.5\n", + "27-Oct 06:20:34:INFO:root:derivative of basis at [0.0,0.0]: [-0.5 0.5 0.5 -0.5\n", " -0.5 -0.5 0.5 0.5]\n", - "26-Oct 05:30:32:INFO:root:field val at [0.0,0.0]: [0.0 2.0]\n", - "26-Oct 05:30:32:INFO:root:Element MyQuad4 passed tests.\n" + "27-Oct 06:20:34:INFO:root:field val at [0.0,0.0]: [0.0 2.0]\n", + "27-Oct 06:20:34:INFO:root:Element MyQuad4 passed tests.\n" ] } ], @@ -238,7 +352,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 13, "metadata": { "collapsed": false }, @@ -246,14 +360,14 @@ { "data": { "text/plain": [ - "Dict{Symbol,JuliaFEM.FieldSet} with 4 entries:\n", - " symbol(\"heat coefficien… => JuliaFEM.FieldSet(symbol(\"heat coefficient\"),Juli…\n", - " :geometry => JuliaFEM.FieldSet(:geometry,JuliaFEM.Field[JuliaF…\n", - " :temperature => JuliaFEM.FieldSet(:temperature,JuliaFEM.Field[Jul…\n", - " :displacement => JuliaFEM.FieldSet(:displacement,JuliaFEM.Field[Ju…" + "Dict{ASCIIString,JuliaFEM.FieldSet} with 4 entries:\n", + " \"geometry\" => JuliaFEM.FieldSet(\"geometry\",JuliaFEM.Field[JuliaFEM.Fi…\n", + " \"heat coefficient\" => JuliaFEM.FieldSet(\"heat coefficient\",JuliaFEM.Field[Jul…\n", + " \"displacement\" => JuliaFEM.FieldSet(\"displacement\",JuliaFEM.Field[JuliaFE…\n", + " \"temperature\" => JuliaFEM.FieldSet(\"temperature\",JuliaFEM.Field[JuliaFEM…" ] }, - "execution_count": 9, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } @@ -290,7 +404,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 14, "metadata": { "collapsed": false }, @@ -303,7 +417,7 @@ " 0.5" ] }, - "execution_count": 10, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -317,7 +431,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 15, "metadata": { "collapsed": false }, @@ -329,7 +443,7 @@ " 0.0 1.0" ] }, - "execution_count": 11, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -342,7 +456,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 16, "metadata": { "collapsed": false }, @@ -353,7 +467,7 @@ "2.5" ] }, - "execution_count": 12, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -365,7 +479,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 17, "metadata": { "collapsed": false }, @@ -376,7 +490,7 @@ "3" ] }, - "execution_count": 13, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -388,7 +502,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -401,7 +515,7 @@ " 0.0703125 0.0078125" ] }, - "execution_count": 14, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -426,7 +540,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -454,7 +568,7 @@ " 1.0)" ] }, - "execution_count": 15, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } @@ -520,27 +634,15 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 20, "metadata": { "collapsed": false }, - "outputs": [ - { - "data": { - "text/plain": [ - "get_unknown_field_name (generic function with 4 methods)" - ] - }, - "execution_count": 16, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "using JuliaFEM: Equation, IntegrationPoint, Quad4, Seg2, get_unknown_field_name\n", + "using JuliaFEM: Equation, IntegrationPoint, Quad4, Seg2\n", "\n", - "abstract Heat <: Equation\n", - "JuliaFEM.get_unknown_field_name(eq::Heat) = symbol(\"temperature\")" + "abstract Heat <: Equation" ] }, { @@ -552,7 +654,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 21, "metadata": { "collapsed": false }, @@ -560,10 +662,10 @@ { "data": { "text/plain": [ - "size (generic function with 64 methods)" + "size (generic function with 74 methods)" ] }, - "execution_count": 17, + "execution_count": 21, "metadata": {}, "output_type": "execute_result" } @@ -573,7 +675,6 @@ "type DC2D4 <: Heat\n", " element :: Quad4\n", " integration_points :: Array{IntegrationPoint, 1}\n", - " global_dofs :: Array{Int64, 1}\n", "end\n", "function DC2D4(element::Quad4)\n", " integration_points = [\n", @@ -582,23 +683,22 @@ " IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]\n", " push!(element, FieldSet(\"temperature\"))\n", - " DC2D4(element, integration_points, [])\n", + " DC2D4(element, integration_points)\n", "end\n", - "Base.size(equation::DC2D4) = 4\n", + "Base.size(equation::DC2D4) = (1, 4)\n", "\n", "\"\"\" Diffusive heat transfer for 2-node linear segment. \"\"\"\n", "type DC2D2 <: Heat\n", " element :: Seg2\n", " integration_points :: Array{IntegrationPoint, 1}\n", - " global_dofs :: Array{Int64, 1}\n", "end\n", "function DC2D2(element::Seg2)\n", " integration_points = [\n", " IntegrationPoint([0.0], 2.0)]\n", " push!(element, FieldSet(\"temperature\"))\n", - " DC2D2(element, integration_points, [])\n", + " DC2D2(element, integration_points)\n", "end\n", - "Base.size(equation::DC2D2) = 2" + "Base.size(equation::DC2D2) = (1, 2)" ] }, { @@ -612,7 +712,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 22, "metadata": { "collapsed": false }, @@ -623,7 +723,7 @@ "has_force_vector (generic function with 3 methods)" ] }, - "execution_count": 18, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } @@ -635,23 +735,23 @@ "function JuliaFEM.get_mass_matrix(equation::DC2D4, ip, time)\n", " element = get_element(equation)\n", " basis = get_basis(element)\n", - " ρ = basis(\"density\", ip, time)\n", - " return ρ * basis(ip,time)'*basis(ip,time)\n", + " ρ = basis(\"density\")\n", + " return ρ(ip, time) * basis(ip,time)'*basis(ip,time)\n", "end\n", "\"\"\" Left hand side defined in integration point. \"\"\"\n", "function JuliaFEM.get_stiffness_matrix(equation::DC2D4, ip, time)\n", " element = get_element(equation)\n", " basis = get_basis(element)\n", " dbasis = grad(basis)\n", - " k = basis(\"temperature thermal conductivity\", ip, time)\n", - " return k * dbasis(ip,time)'*dbasis(ip,time)\n", + " k = basis(\"temperature thermal conductivity\")\n", + " return k(ip,time) * dbasis(ip,time)'*dbasis(ip,time)\n", "end\n", "\"\"\" Right hand side defined in integration point. \"\"\"\n", "function JuliaFEM.get_force_vector(equation::DC2D4, ip, time)\n", " element = get_element(equation)\n", " basis = get_basis(element)\n", - " f = basis(\"temperature load\", ip, time)\n", - " return basis(ip,time)'*f\n", + " f = basis(\"temperature load\")\n", + " return f(ip, time)*basis(ip,time)'\n", "end\n", "JuliaFEM.has_mass_matrix(equation::DC2D4) = true\n", "JuliaFEM.has_stiffness_matrix(equation::DC2D4) = true\n", @@ -661,8 +761,8 @@ "function JuliaFEM.get_force_vector(equation::DC2D2, ip, time)\n", " element = get_element(equation)\n", " basis = get_basis(element)\n", - " g = basis(\"temperature flux\", ip, time)\n", - " return basis(ip,time)'*g\n", + " g = basis(\"temperature flux\")\n", + " return g(ip,time)*basis(ip,time)'\n", "end\n", "JuliaFEM.has_force_vector(equation::DC2D2) = true" ] @@ -682,7 +782,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 23, "metadata": { "collapsed": false }, @@ -710,7 +810,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 24, "metadata": { "collapsed": false }, @@ -725,7 +825,7 @@ " -1.0 -2.0 -1.0 4.0" ] }, - "execution_count": 20, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" } @@ -733,14 +833,14 @@ "source": [ "using JuliaFEM: initialize_local_assembly, calculate_local_assembly!\n", "\n", - "local_assembly = initialize_local_assembly(equation)\n", - "calculate_local_assembly!(local_assembly, equation)\n", + "local_assembly = initialize_local_assembly()\n", + "calculate_local_assembly!(local_assembly, equation, \"temperature\")\n", "local_assembly.stiffness_matrix" ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 25, "metadata": { "collapsed": false }, @@ -755,7 +855,7 @@ " 3.0" ] }, - "execution_count": 21, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" } @@ -773,7 +873,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 26, "metadata": { "collapsed": false }, @@ -786,7 +886,7 @@ " 1.0" ] }, - "execution_count": 22, + "execution_count": 26, "metadata": {}, "output_type": "execute_result" } @@ -808,7 +908,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 27, "metadata": { "collapsed": false }, @@ -821,16 +921,16 @@ " 1.0" ] }, - "execution_count": 23, + "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "boundary_assembly = JuliaFEM.initialize_local_assembly(boundary_equation)\n", - "JuliaFEM.calculate_local_assembly!(boundary_assembly, boundary_equation)\n", + "local_assembly = initialize_local_assembly()\n", + "calculate_local_assembly!(local_assembly, boundary_equation, \"temperature\")\n", "b = zeros(4)\n", - "b[fdofs] = boundary_assembly.force_vector\n", + "b[fdofs] = local_assembly.force_vector\n", "u = zeros(4)\n", "u[fdofs] = A[fdofs, fdofs] \\ b[fdofs]" ] @@ -845,7 +945,7 @@ "\n", "Let's consider the following functional\n", "\\begin{equation}\n", - "\\min\\, J\\left(u\\right)=\\int_{\\Omega}(k+6u)\\left|\\nabla u\\right|^{2}\\,\\mathrm{d}x-Pu,\n", + "\\min\\, J\\left(u\\right)=\\int_{\\Omega}(k+6u)\\left(\\nabla u \\cdot \\nabla u\\right) \\, \\mathrm{d}x-Pu,\n", "\\end{equation}\n", "where $P$ contains point loads at the free corners of the domain. This is nothing more but a redefinition of the earlier example problem defined such that the source term $k$, which in previous example was constant, is no more constant, but is replaced a term which depend from field $u$: $q(u) = k + 6u$.\n", "\n", @@ -859,7 +959,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 28, "metadata": { "collapsed": false }, @@ -867,10 +967,10 @@ { "data": { "text/plain": [ - "size (generic function with 65 methods)" + "size (generic function with 75 methods)" ] }, - "execution_count": 24, + "execution_count": 28, "metadata": {}, "output_type": "execute_result" } @@ -880,7 +980,6 @@ "type DC2D4NL <: Heat\n", " element :: Quad4\n", " integration_points :: Array{IntegrationPoint, 1}\n", - " global_dofs :: Array{Int64, 1}\n", "end\n", "function DC2D4NL(element::Quad4)\n", " integration_points = [\n", @@ -891,9 +990,9 @@ " push!(element, FieldSet(\"temperature\"))\n", " # Initial configuration needs to be defined\n", " push!(element[\"temperature\"], Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n", - " DC2D4NL(element, integration_points, [])\n", + " DC2D4NL(element, integration_points)\n", "end\n", - "Base.size(equation::DC2D4NL) = 4" + "Base.size(equation::DC2D4NL) = (1, 4)" ] }, { @@ -905,7 +1004,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 29, "metadata": { "collapsed": false }, @@ -916,7 +1015,7 @@ "has_potential_energy (generic function with 2 methods)" ] }, - "execution_count": 32, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } @@ -930,7 +1029,7 @@ " f = basis(\"temperature load\", ip, time)\n", " T = basis(\"temperature\", ip, time, variation)\n", " ∇T = grad(basis)(\"temperature\", ip, time, variation)\n", - " Wint = (k + 6*T) * 1/2*∇T*∇T'\n", + " Wint = (k + 6*T) * 1/2*vecdot(∇T, ∇T)\n", " Wext = f*T\n", " return Wint - Wext\n", "end\n", @@ -946,7 +1045,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 30, "metadata": { "collapsed": false }, @@ -960,9 +1059,9 @@ "increment 3, |du| = 0.04657, |r| = 4.243\n", "increment 4, |du| = 0.00057, |r| = 4.243\n", "increment 5, |du| = 0.00000, |r| = 4.243\n", - "elapsed time: 0.001905274 seconds\n", - "error: 1.6653345369377348e-15\n", - "temperature at free end: [0.6666666666666683,0.6666666666666683], should be 0.6666666666666666\n" + "elapsed time: 0.611092931 seconds\n", + "error: 1.5543122344752192e-15\n", + "temperature at free end: [0.6666666666666682,0.6666666666666682], should be 0.6666666666666666\n" ] } ], @@ -987,15 +1086,14 @@ " # create model -- end\n", "\n", " T0 = element[\"temperature\"][1] # initial temperature field, we need something to \"variate\"\n", - " la = initialize_local_assembly(equation) # create workspace for local matrices\n", + " la = initialize_local_assembly() # create workspace for local matrices\n", " T = zeros(4) # create workspace for solution vector\n", " ΔT = zeros(4) # \n", " fd = [1, 2]\n", " tic()\n", " # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.\n", " for i=1:5\n", - " la = initialize_local_assembly(equation,la) # empty workspace -- every iteration should start with this\n", - " calculate_local_assembly!(la, equation) # calculate local matrices\n", + " calculate_local_assembly!(la, equation, \"temperature\") # calculate local matrices\n", " ΔT[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd] # <-- note sign convention, more on this below\n", " T = T + ΔT # add increment to previous value\n", " # create a new field \"similar\" to field T0 (i.e., same dimension of field variable with new data)\n", @@ -1046,14 +1144,23 @@ "\n", "### Principle of virtual work\n", "\n", - "Defining variational form or \"energy form\", or in other words, using principle of minimum potential energy, requires system to be *conservative*, i.e. there is some energy functional describing the system. Indeed it's a very elegant approach to model e.g. hyperelasticity and other systems without a loss of energy. When system has energy dissipation, principle of minimum potential energy cannot be used for obvious reasons. There is no any \"potential function\" $\\Pi$ which could be defined and variated around it's equilibrium state. Possible situations when this happens include e.g. material plasticity or friction. In these situations the principle of virtual work is useful.\n", + "Defining variational form or \"energy form\", or in other words, using principle of minimum potential energy, requires system to be *conservative*, i.e. there is some energy functional describing the system. Indeed it's a very elegant approach to model e.g. hyperelasticity and other systems without a loss of energy. When system has energy dissipation, principle of minimum potential energy cannot be used. Possible situations when this happens include e.g. material plasticity or friction. In these situations the principle of virtual work is useful.\n", "\n", - "In this situation we find equation $\\mathbf{r}=\\mathbf{f}_{\\mathrm{int}}-\\mathbf{f}_{\\mathrm{ext}}=0$ which needs to be solved. Again we let `ForwardDiff` to do the linearization of the right hand side. Let's demonstrate this too." + "In this situation we find equation $\\mathbf{r}=\\mathbf{f}_{\\mathrm{int}}-\\mathbf{f}_{\\mathrm{ext}}=0$ which needs to be solved. Again we let the `ForwardDiff` do the linearization of the right hand side. We demonstrate this by constructing geometrically non-linear elasticity equation for 4-node bilinear element. Virtual work now looks\n", + "\n", + "\\begin{equation}\n", + "\\delta\\mathcal{W} := \\int_{\\Omega_{0}}\\boldsymbol{S}:\\delta\\boldsymbol{E}\\,\\mathrm{d}V_{0}-\\int_{\\Omega_{0}}\\boldsymbol{b}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0} = 0\n", + "\\end{equation}\n", + "where\n", + "\\begin{equation}\n", + "\\delta\\mathcal{W}_{\\mathrm{int}}=\\int\\mathbf{S}:\\delta\\mathbf{E}\\,\\mathrm{d}V_{0}=\\int_{\\Omega_{0}}\\mathbf{P}\\cdot\\mbox{Grad}\\delta\\mathbf{u}\\,\\mathrm{d}V_{0}\n", + "\\end{equation}\n", + "and $\\boldsymbol{b}_{0}$ is a volume load defined in reference configuration." ] }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 31, "metadata": { "collapsed": false }, @@ -1064,7 +1171,7 @@ "has_residual_vector (generic function with 2 methods)" ] }, - "execution_count": 37, + "execution_count": 31, "metadata": {}, "output_type": "execute_result" } @@ -1072,44 +1179,64 @@ "source": [ "using JuliaFEM: get_field\n", "\n", - "\"\"\" Diffusive heat transfer for 4-node bilinear element,\n", - "with a nonlinear term to be added later.. \"\"\"\n", - "type DC2D4NLY <: Heat\n", + "abstract Elasticity <: Equation\n", + "abstract PlaneElasticity <: Elasticity\n", + "abstract PlaneStressElasticity <: PlaneElasticity\n", + "\n", + "\"\"\" Plane stress formulation for 4-node bilinear element. \"\"\"\n", + "type CPS4 <: PlaneStressElasticity\n", " element :: Quad4\n", " integration_points :: Array{IntegrationPoint, 1}\n", - " global_dofs :: Array{Int64, 1}\n", "end\n", - "function DC2D4NLY(element::Quad4)\n", + "\n", + "function CPS4(element::Quad4)\n", " integration_points = [\n", " IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]\n", - " push!(element, FieldSet(\"temperature\"))\n", - " # Initial configuration needs to be defined\n", - " push!(element[\"temperature\"], Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n", - " DC2D4NLY(element, integration_points, [])\n", + " push!(element, FieldSet(\"displacement\"))\n", + " CPS4(element, integration_points)\n", "end\n", - "Base.size(equation::DC2D4NLY) = 4\n", "\n", - "function JuliaFEM.get_residual_vector(equation::DC2D4NLY, ip, time; variation=nothing)\n", + "JuliaFEM.size(eq::CPS4) = (2, 4)\n", + "\n", + "function JuliaFEM.get_residual_vector(equation::CPS4, ip, time; variation=nothing)\n", " element = get_element(equation)\n", " basis = get_basis(element)\n", " dbasis = grad(basis)\n", - " k = basis(\"temperature thermal conductivity\", ip, time)\n", - " f = basis(\"temperature load\", ip, time)\n", - " T = get_field(basis, \"temperature\", time, variation)\n", - " f_int = k * dbasis(ip,time)'*dbasis(ip,time) * T\n", - " f_ext = f * basis(ip,time)'\n", - " r = f_int[:] - f_ext[:]\n", - " return r\n", + "\n", + " # material parameters\n", + " E = basis(\"youngs modulus\", ip, time)\n", + " ν = basis(\"poissons ratio\", ip, time)\n", + " μ = E/(2*(1+ν))\n", + " λ = E*ν/((1+ν)*(1-2*ν))\n", + " λ = 2*λ*μ/(λ + 2*μ) # <- correction for 2d\n", + "\n", + " # elasticity formulation\n", + " u = basis(\"displacement\", ip, time, variation)\n", + " ∇u = dbasis(\"displacement\", ip, time, variation)\n", + " F = I + ∇u\n", + " b = basis(\"displacement volume load\", ip, time)\n", + " E = 1/2*(F'*F - I)\n", + " S = λ*trace(E)*I + 2*μ*E\n", + " J = det(F)\n", + " σ = J^-1 * F*S*F'\n", + " #push!(ip[\"cauchy stress\"], Field(time, σ))\n", + " P = F*S\n", + "\n", + " # residual vector\n", + " r_int = P*dbasis(ip,time)\n", + " r_ext = b*basis(ip,time)\n", + " r = r_int - r_ext\n", + " return vec(r)\n", "end\n", - "JuliaFEM.has_residual_vector(equation::DC2D4NLY) = true" + "JuliaFEM.has_residual_vector(equation::CPS4) = true" ] }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 32, "metadata": { "collapsed": false }, @@ -1124,10 +1251,13 @@ { "data": { "text/plain": [ - "0.020793035" + "Success :: (line:-1) :: fact was true\n", + " Expression: disp --> roughly(-8.77303119819776)\n", + " Expected: -8.77303119819776\n", + " Occurred: -8.773031198197748" ] }, - "execution_count": 38, + "execution_count": 32, "metadata": {}, "output_type": "execute_result" }, @@ -1135,58 +1265,60 @@ "name": "stdout", "output_type": "stream", "text": [ - "1, |du| = 2.82843, |r| = 4.243\n", - "increment 2, |du| = 0.00000, |r| = 4.243\n", - "increment 3, |du| = 0.00000, |r| = 4.243\n", - "increment 4, |du| = 0.00000, |r| = 4.243\n", - "increment 5, |du| = 0.00000, |r| = 4.243\n", - "temperature field: [2.0000000000000004,2.0000000000000004,0.0,0.0]\n", - "elapsed time: 0.020793035 seconds\n" + "1, |du| = 14.44128\n", + "increment 2, |du| = 4.01742\n", + "increment 3, |du| = 1.54645\n", + "increment 4, |du| = 1.12361\n", + "increment 5, |du| = 0.79119\n", + "increment 6, |du| = 0.12733\n", + "increment 7, |du| = 0.00725\n", + "increment 8, |du| = 0.00001\n", + "increment 9, |du| = 0.00000\n", + "increment 10, |du| = 0.00000\n", + "elapsed time: 1.770659582 seconds\n", + "displacement at tip: -8.773031198197748\n" ] } ], "source": [ + "using JuliaFEM: get_integration_points\n", + "\n", "function run_simulation_2()\n", " # create model -- start\n", " element = Quad4([1, 2, 3, 4])\n", - " fieldset1 = FieldSet(\"geometry\")\n", - " field1 = Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])\n", - " push!(fieldset1, field1)\n", - " fieldset2 = FieldSet(\"temperature thermal conductivity\")\n", - " push!(fieldset2, Field(0.0, 6.0))\n", - " fieldset3 = FieldSet(\"temperature load\")\n", - " push!(fieldset3, Field(0.0, [12.0, 12.0, 12.0, 12.0]))\n", - " fieldset4 = FieldSet(\"temperature nodal load\")\n", - " push!(fieldset4, Field(0.0, [3.0, 3.0, 0.0, 0.0])) # <-- P is defined here\n", - " push!(element, fieldset1)\n", - " push!(element, fieldset2)\n", - " push!(element, fieldset3)\n", - " push!(element, fieldset4)\n", - " equation = DC2D4NLY(element)\n", + " push!(element, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]])]))\n", + " push!(element, FieldSet(\"youngs modulus\", [Field(0.0, 500.0)]))\n", + " push!(element, FieldSet(\"poissons ratio\", [Field(0.0, 0.3)]))\n", + " push!(element, FieldSet(\"displacement volume load\",\n", + " [Field(0.0, Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]])]))\n", + " equation = CPS4(element)\n", + " #for ip in get_integration_points(equation)\n", + " # push!(ip, FieldSet(\"cauchy stress\"))\n", + " #end\n", " # create model -- end\n", "\n", - " T0 = element[\"temperature\"][1] # initial temperature field, we need something to \"variate\"\n", - " la = initialize_local_assembly(equation) # create workspace for local matrices\n", - " T = zeros(4) # create workspace for solution vector\n", - " ΔT = zeros(4) # \n", - " fd = [1, 2]\n", + " u0 = Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])\n", + " push!(element[\"displacement\"], u0)\n", + " u = zeros(8)\n", + " du = zeros(8)\n", + " fd = [3, 4, 5, 6]\n", + " la = initialize_local_assembly()\n", " tic()\n", - " # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.\n", - " for i=1:5\n", - " la = initialize_local_assembly(equation,la) # empty workspace -- every iteration should start with this\n", - " calculate_local_assembly!(la, equation) # calculate local matrices\n", - " ΔT[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd] # <-- note sign convention, more on this below\n", - " T = T + ΔT # add increment to previous value\n", - " # create a new field \"similar\" to field T0 (i.e., same dimension of field variable with new data)\n", - " new_field = similar(T0, T)\n", + " for i=1:10\n", + " calculate_local_assembly!(la, equation, \"displacement\")\n", + " du[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd]\n", + " u += du\n", + " new_field = similar(u0, u)\n", " new_field.time = 1.0\n", " new_field.increment = i\n", - " push!(element[\"temperature\"], new_field) # add new field to \"temperature\" fieldset of element\n", - " # print some convergence information\n", - " @printf(\"increment %2d, |du| = %8.5f, |r| = %8.3f\\n\", i, norm(ΔT), norm(b[fd]))\n", + " push!(element[\"displacement\"], new_field)\n", + " @printf(\"increment %2d, |du| = %8.5f\\n\", i, norm(du))\n", " end\n", - " println(\"temperature field: \",element[\"temperature\"](Inf).values)\n", " toc()\n", + " disp = get_basis(element)(\"displacement\", [1.0, 1.0])[2]\n", + " println(\"displacement at tip: $disp\")\n", + " # verified using Code Aster.\n", + " @fact disp --> roughly(-8.77303119819776E+00)\n", "end\n", "run_simulation_2()" ] @@ -1199,12 +1331,12 @@ "\n", "There might still be situations where all the above methods writing field equations are just not enough. This situation can happen for example when calculating tangent stiffness matrix analytically for a nonlinear problem. In this case all matrices are usually calculated at the same time. Or maybe for educational purposes it's important to show how matrices are actually calculated. Or for debugging. \n", "\n", - "Anyway, it's possible to override `calculate_local_assembly!` function for your own equation and after that access all the low level stuff. Here's example how to do that:" + "Anyway, it's possible to override `calculate_local_assembly!` function for your own equation and after that access all the low level stuff. Here's an example how to do that. " ] }, { "cell_type": "code", - "execution_count": 39, + "execution_count": 33, "metadata": { "collapsed": false }, @@ -1218,16 +1350,20 @@ " Occurred: [1.0000000000000002,1.0000000000000002]" ] }, - "execution_count": 39, + "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "using JuliaFEM: LocalAssembly, get_integration_points\n", + "using JuliaFEM: LocalAssembly, get_integration_points, initialize_local_assembly!\n", + "\n", + "function JuliaFEM.calculate_local_assembly!(assembly::LocalAssembly, equation::DC2D4,\n", + " unknown_field_name::ASCIIString, time::Number=Inf,\n", + " problem=nothing)\n", + " # zero all workspace. always start with this command\n", + " initialize_local_assembly!(assembly, equation)\n", "\n", - "function JuliaFEM.calculate_local_assembly!(assembly::LocalAssembly, equation::DC2D4, time::Number=Inf)\n", - " initialize_local_assembly(assembly, equation) # zero all workspaces\n", " element = get_element(equation)\n", " basis = get_basis(element)\n", " dbasis = grad(basis)\n", @@ -1244,7 +1380,7 @@ " # do assembly\n", " assembly.mass_matrix += w * ρ*N'*N\n", " assembly.stiffness_matrix += w * k*∇N'*∇N\n", - " assembly.force_vector += w * (N'*f)[:]\n", + " assembly.force_vector += w * N'*f\n", " end\n", "end\n", "\n", @@ -1260,8 +1396,8 @@ " push!(element, fieldset3)\n", " push!(element, fieldset4)\n", " equation = DC2D4(element)\n", - " la = initialize_local_assembly(equation)\n", - " calculate_local_assembly!(la, equation)\n", + " la = initialize_local_assembly()\n", + " calculate_local_assembly!(la, equation, \"temperature\")\n", " free_dofs = [1, 2]\n", " @fact la.stiffness_matrix[free_dofs, free_dofs] \\ la.force_vector[free_dofs] --> roughly([1.0, 1.0])\n", "end\n", @@ -1269,20 +1405,13 @@ "test_local_assembly()" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To be continued (stuff below this is a little broken at the moment...)" - ] - }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining own problem\n", "\n", - "Main objective of \"problem\" is to take a set of elements and map corresponding field equations to them. In this way we can solve several different fields at the same time, like for example temperature + displacement." + "The main purpose of `Problem` is to create a mapping between element types and field equations, so that solver knows how to construct equations. Notice that there is no any limits of how many different equations can be set to one element, so it's totally possible for example solve displacement and temperature equations with same element and even so that problems are loosely coupled, i.e. using temperature field to calculate stresses in mechanical problem." ] }, { @@ -1295,7 +1424,7 @@ { "data": { "text/plain": [ - "PlaneHeatProblem" + "initialize_global_assembly (generic function with 3 methods)" ] }, "execution_count": 34, @@ -1304,12 +1433,65 @@ } ], "source": [ - "using JuliaFEM: Problem, get_equation, get_dimension\n", + "using JuliaFEM: Problem, FieldProblem, get_equations, get_connectivity\n", + "using JuliaFEM: get_element, get_unknown_field_dimension\n", "\n", - "type PlaneHeatProblem <: Problem\n", + "type PlaneHeatProblem <: FieldProblem\n", + " unknown_field_name :: ASCIIString\n", + " unknown_field_dimension :: Int\n", " equations :: Array{Equation, 1}\n", + " element_mapping :: Dict{DataType, DataType}\n", "end\n", - "PlaneHeatProblem() = PlaneHeatProblem([])" + "function PlaneHeatProblem()\n", + " element_mapping = Dict(\n", + " Quad4 => DC2D4,\n", + " Seg2 => DC2D2)\n", + " PlaneHeatProblem(\"temperature\", 1, [], element_mapping)\n", + "end\n", + "\n", + "function test_problem()\n", + " el1 = Quad4([1, 2, 3, 4])\n", + " push!(el1, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]))\n", + " push!(el1, FieldSet(\"temperature thermal conductivity\", [Field(0.0, 6.0)]))\n", + " push!(el1, FieldSet(\"temperature load\", [Field(0.0, [12.0, 12.0, 12.0, 12.0])]))\n", + " push!(el1, FieldSet(\"density\", [Field(0.0, 10.0)]))\n", + " el2 = Seg2([1, 2])\n", + " push!(el2, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]])]))\n", + "\n", + " # Boundary load, linear ramp 0 -> 600 at time 0 -> 1\n", + " load = FieldSet(\"temperature flux\")\n", + " push!(load, Field(0.0, 0.0))\n", + " push!(load, Field(1.0, 600.0))\n", + " push!(el2, load)\n", + "\n", + " problem = PlaneHeatProblem()\n", + " push!(problem, el1)\n", + " push!(problem, el2)\n", + "\n", + " problem\n", + " \n", + "end\n", + "\n", + "using JuliaFEM: Assembly\n", + "\n", + "type GlobalAssembly <: Assembly\n", + " ndofs :: Int\n", + " mass_matrix :: SparseMatrixCSC\n", + " stiffness_matrix :: SparseMatrixCSC\n", + " force_vector :: SparseMatrixCSC\n", + "end\n", + "\n", + "function initialize_global_assembly(ndofs::Int=1)\n", + " mass_matrix = spzeros(ndofs, ndofs)\n", + " stiffness_matrix = spzeros(ndofs, ndofs)\n", + " force_vector = spzeros(ndofs, 1)\n", + " return GlobalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector)\n", + "end\n", + "\n", + "function initialize_global_assembly(problem::Problem)\n", + " dim, ndofs = size(problem)\n", + " return initialize_global_assembly(ndofs)\n", + "end" ] }, { @@ -1319,365 +1501,67 @@ "collapsed": false }, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "assembling problem for temperature\n" + ] + }, { "data": { "text/plain": [ - "get_equation (generic function with 6 methods)" + "2-element Array{Float64,1}:\n", + " 101.0\n", + " 101.0" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" - } - ], - "source": [ - "JuliaFEM.get_dimension(pr::Type{PlaneHeatProblem}) = 1\n", - "JuliaFEM.get_equation(pr::Type{PlaneHeatProblem}, el::Type{Quad4}) = DC2D4\n", - "JuliaFEM.get_equation(pr::Type{PlaneHeatProblem}, el::Type{Seg2}) = DC2D2" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Our solution procedure so far is therefore" - ] - }, - { - "cell_type": "code", - "execution_count": 33, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "25-Oct 23:10:39:DEBUG:root:total dofs: 4\n" - ] }, { - "ename": "LoadError", - "evalue": "LoadError: MethodError: `has_lhs` has no method matching has_lhs(::DC2D4)\nwhile loading In[33], in expression starting on line 30", - "output_type": "error", - "traceback": [ - "LoadError: MethodError: `has_lhs` has no method matching has_lhs(::DC2D4)\nwhile loading In[33], in expression starting on line 30", - "" + "name": "stdout", + "output_type": "stream", + "text": [ + "dimension of unknown field: 1, problem dofs: 4\n" ] } ], "source": [ - "using JuliaFEM: get_connectivity, set_global_dofs!, get_global_dofs\n", - "using JuliaFEM: add_element!, get_equations, get_matrix_dimension, calculate_global_dofs\n", - "using JuliaFEM: assign_global_dofs!\n", + "function calculate_global_assembly!(assembly::GlobalAssembly, problem::Problem,\n", + " unknown_field_name::ASCIIString, time::Number=Inf)\n", + " dim, ndofs = size(problem)\n", + " println(\"assembling problem for $unknown_field_name\")\n", + " println(\"dimension of unknown field: $dim, problem dofs: $ndofs\")\n", + " local_assembly = initialize_local_assembly()\n", + " for (i, equation) in enumerate(get_equations(problem))\n", + " calculate_local_assembly!(local_assembly, equation, unknown_field_name, time)\n", + " conn = get_connectivity(get_element(equation))\n", + " gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))\n", + " assembly.mass_matrix[gdofs, gdofs] += local_assembly.mass_matrix\n", + " assembly.stiffness_matrix[gdofs, gdofs] += local_assembly.stiffness_matrix\n", + " assembly.force_vector[gdofs] += local_assembly.force_vector\n", + " end\n", + "end\n", "\n", - "# create elements and add necessary properties like connectivity and geometry\n", - "el1 = Quad4([2, 3, 4, 5])\n", - "fs1geom = FieldSet(\"geometry\")\n", - "push!(fs1geom, Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n", - "fs1temp = FieldSet(\"temperature thermal conductivity\")\n", - "push!(fs1temp, Field(0.0, 6.0))\n", - "push!(el1, fs1geom)\n", - "push!(el1, fs1temp)\n", + "function Base.size(problem::Problem)\n", + " mc = 0\n", + " for equation in get_equations(problem)\n", + " element = get_element(equation)\n", + " mc = max(mc, get_connectivity(element)...)\n", + " end\n", + " dim = get_unknown_field_dimension(problem)\n", + " return (dim, dim*mc)\n", + "end\n", "\n", - "el2 = Seg2([2, 3])\n", - "fs2geom = FieldSet(\"geometry\")\n", - "push!(fs2geom, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n", - "fs2temp = FieldSet(\"temperature flux\")\n", - "push!(fs2temp, Field(1.0, 600.0))\n", - "push!(el2, fs2geom)\n", - "push!(el2, fs2temp)\n", - "\n", - "problem = PlaneHeatProblem()\n", - "push!(problem, el1)\n", - "push!(problem, el2)\n", - "\n", - "dofmap = calculate_global_dofs(problem)\n", - "assign_global_dofs!(problem, dofmap)\n", - "\n", - "t = 1.0\n", - "A = sparse(get_lhs(problem, t)...)\n", - "b = sparsevec(get_rhs(problem, t)..., size(A, 1))" - ] - }, - { - "cell_type": "code", - "execution_count": 30, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: A not defined\nwhile loading In[30], in expression starting on line 2", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: A not defined\nwhile loading In[30], in expression starting on line 2", - "" - ] - } - ], - "source": [ - "fdofs = [1, 2]\n", - "A[fdofs, fdofs] \\ b[fdofs]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We still need to consider Dirichlet boundary conditions:\n", - "\\begin{align}\n", - "u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n", - "\\end{align}" - ] - }, - { - "cell_type": "code", - "execution_count": 31, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "1-element Array{JuliaFEM.DirichletEquation,1}:\n", - " JuliaFEM.DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.FieldSet(symbol(\"reaction force\"),JuliaFEM.Field[]),:geometry=>JuliaFEM.FieldSet(:geometry,JuliaFEM.Field[JuliaFEM.Field{Array{Array{T,1},1}}(0.0,0,Array{T,1}[[0.0,0.0],[0.0,1.0]])]))),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Symbol,JuliaFEM.FieldSet}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Symbol,JuliaFEM.FieldSet}())],Int64[],fieldval)" - ] - }, - "execution_count": 31, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "using JuliaFEM: DirichletProblem\n", - "\n", - "# create elements and add necessary properties like connectivity and geometry\n", - "el3 = Seg2([4, 5])\n", - "push!(el3, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]])]))\n", - "\n", - "bc1 = DirichletProblem()\n", - "push!(bc1, el3)" - ] - }, - { - "cell_type": "code", - "execution_count": 32, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: integrate_lhs not defined\nwhile loading In[32], in expression starting on line 5", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: integrate_lhs not defined\nwhile loading In[32], in expression starting on line 5", - "", - " in get_lhs at /home/jukka/.julia/v0.4/JuliaFEM/src/problems.jl:119" - ] - } - ], - "source": [ - "assign_global_dofs!(bc1, dofmap)\n", - "\n", - "# integrate and assembly\n", - "t = 1.0\n", - "A2 = sparse(get_lhs(bc1, t)...)\n", - "b2 = sparsevec(get_rhs(bc1, t)...)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "collapsed": false - }, - "source": [ - "Now we have two problems defined, " - ] - }, - { - "cell_type": "code", - "execution_count": 33, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: A not defined\nwhile loading In[33], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: A not defined\nwhile loading In[33], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "full(A)" - ] - }, - { - "cell_type": "code", - "execution_count": 34, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: b not defined\nwhile loading In[34], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: b not defined\nwhile loading In[34], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "full(b)" - ] - }, - { - "cell_type": "code", - "execution_count": 35, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: A2 not defined\nwhile loading In[35], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: A2 not defined\nwhile loading In[35], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "full(A2)" - ] - }, - { - "cell_type": "code", - "execution_count": 36, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: b2 not defined\nwhile loading In[36], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: b2 not defined\nwhile loading In[36], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "full(b2)" - ] - }, - { - "cell_type": "code", - "execution_count": 37, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: A not defined\nwhile loading In[37], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: A not defined\nwhile loading In[37], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "Atot = [A A2; A2 zeros(A2)]" - ] - }, - { - "cell_type": "code", - "execution_count": 38, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: b not defined\nwhile loading In[38], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: b not defined\nwhile loading In[38], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "btot = [b; b2]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Problem is that now are total matrix Atot has zero rows which needs to be removed." - ] - }, - { - "cell_type": "code", - "execution_count": 39, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: Atot not defined\nwhile loading In[39], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: Atot not defined\nwhile loading In[39], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "full(Atot)" - ] - }, - { - "cell_type": "code", - "execution_count": 40, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: UndefVarError: Atot not defined\nwhile loading In[40], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: UndefVarError: Atot not defined\nwhile loading In[40], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "r = unique(rowvals(Atot))\n", - "println(\"Non-zero rows: $r\")\n", - "xtot = zeros(btot)\n", - "F = lufact(Atot[r,r])\n", - "s = full(btot[r])\n", - "xtot[r] = F \\ s\n", - "full(xtot)" + "p = test_problem()\n", + "size(p)\n", + "global_assembly = initialize_global_assembly(p)\n", + "calculate_global_assembly!(global_assembly, p, \"temperature\")\n", + "global_assembly.stiffness_matrix, global_assembly.force_vector\n", + "free_dofs = [1, 2]\n", + "lufact(global_assembly.stiffness_matrix[free_dofs, free_dofs]) \\ full(global_assembly.force_vector)[free_dofs]" ] }, { @@ -1690,72 +1574,71 @@ "\n", "Last part. Defining own solver.\n", "\n", - "- takes a set of problems (typically main problem + boundary problems)\n", + "- takes a set of problems (typically some field problem + dirichlet boundary problems)\n", "- solves them, updates fields" ] }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { - "data": { - "text/plain": [ - "call (generic function with 1288 methods)" - ] - }, - "execution_count": 41, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "assembling problem for temperature\n", + "dimension of unknown field: 1, problem dofs: 4\n", + "assembling problem for reaction force\n", + "dimension of unknown field: 1, problem dofs: 4\n", + "Residual norm: 5.684341886080802e-14\n" + ] } ], "source": [ - "using JuliaFEM: Solver, get_problems\n", + "using JuliaFEM: Solver, get_problems, get_unknown_field_name\n", "\n", "\"\"\" Simple solver for educational purposes. \"\"\"\n", "type SimpleSolver <: Solver\n", - " problems\n", + " problems :: Array{Problem, 1}\n", "end\n", "\n", "\"\"\" Default initializer. \"\"\"\n", "function SimpleSolver()\n", - " SimpleSolver(Problem[])\n", + " SimpleSolver([])\n", "end\n", "\n", "\"\"\"\n", "Call solver to solve a set of problems.\n", "\n", - "This is simple serial solver for demonstration purposes. It handles the most\n", + "This is a simple direct solver for demonstration purposes. It handles the\n", "common situation, i.e., some main field problem and it's Dirichlet boundary.\n", + "\n", + " Au + C'λ = f\n", + " Cu = g\n", + "\n", "\"\"\"\n", - "function call(solver::SimpleSolver, t)\n", - " problems = get_problems(solver)\n", - " problem1 = problems[1]\n", - " problem2 = problems[2]\n", + "function call(solver::SimpleSolver, time::Number=Inf)\n", + " p1, p2 = get_problems(solver)\n", "\n", - " # calculate order of degrees of freedom in global matrix\n", - " # and set the ordering to problems\n", - " dofmap = calculate_global_dofs(problem1)\n", - " assign_global_dofs!(problem1, dofmap)\n", - " assign_global_dofs!(problem2, dofmap)\n", + " ga1 = initialize_global_assembly(p1)\n", + " calculate_global_assembly!(ga1, p1, get_unknown_field_name(p1))\n", + " ga2 = initialize_global_assembly(p2)\n", + " calculate_global_assembly!(ga2, p2, get_unknown_field_name(p2))\n", "\n", - " # assemble problem 1\n", - " A1 = sparse(get_lhs(problem1, t)...)\n", - " b1 = sparsevec(get_rhs(problem1, t)..., size(A1, 1))\n", - "\n", - " # assemble problem 2\n", - " A2 = sparse(get_lhs(problem2, t)...)\n", - " b2 = sparsevec(get_rhs(problem2, t)..., size(A2, 1))\n", + " A1 = ga1.stiffness_matrix\n", + " b1 = ga1.force_vector\n", + " A2 = ga2.stiffness_matrix\n", + " b2 = ga2.force_vector\n", " \n", - " # make one monolithic assembly\n", + " # create a saddle point problem\n", " A = [A1 A2; A2' zeros(A2)]\n", " b = [b1; b2]\n", "\n", " # solve problem\n", - " nz = unique(rowvals(A))\n", + " nz = unique(rowvals(A)) # here we remove any zero rows\n", " x = zeros(b)\n", " x[nz] = lufact(A[nz,nz]) \\ full(b[nz])\n", "\n", @@ -1769,23 +1652,45 @@ " println(\"Residual norm: $(norm(R1+R2))\")\n", "\n", " # update field for elements in problem 1\n", - " for equation in get_equations(problem1)\n", - " gdofs = get_global_dofs(equation)\n", + " for equation in get_equations(p1)\n", + " dim = get_unknown_field_dimension(p1)\n", + " field_name = get_unknown_field_name(p1)\n", " element = get_element(equation)\n", - " field_name = get_unknown_field_name(equation) # field we are solving\n", - " field = Field(t, full(x1[gdofs])[:])\n", + " conn = get_connectivity(element)\n", + "\n", + " gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))\n", + " field = Field(time, full(x1[gdofs])[:])\n", " push!(element[field_name], field)\n", " end\n", "\n", - " # update field for elements in problem 2\n", - " for equation in get_equations(problem2)\n", - " gdofs = get_global_dofs(equation)\n", + " # update field for elements in problem 2 (Dirichlet boundary)\n", + " for equation in get_equations(p2)\n", + " dim = get_unknown_field_dimension(p2)\n", + " field_name = get_unknown_field_name(p2)\n", " element = get_element(equation)\n", - " field_name = get_unknown_field_name(equation)\n", - " field = Field(t, full(x2[gdofs]))\n", + " conn = get_connectivity(element)\n", + "\n", + " gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))\n", + " field = Field(time, full(x2[gdofs])[:])\n", " push!(element[field_name], field)\n", " end\n", - "end" + "end\n", + "\n", + "using JuliaFEM: DirichletProblem\n", + "\n", + "problem1 = test_problem()\n", + "\n", + "el3 = Seg2([3, 4])\n", + "push!(el3, FieldSet(\"geometry\", [Field(0.0, Vector[[1.0, 1.0], [0.0, 1.0]])]))\n", + "\n", + "problem2 = DirichletProblem(1)\n", + "push!(problem2, el3)\n", + "\n", + "# Create a solver for a set of problems\n", + "solver = SimpleSolver()\n", + "push!(solver, problem1)\n", + "push!(solver, problem2)\n", + "call(solver, 1.0)" ] }, { @@ -1799,123 +1704,112 @@ }, { "cell_type": "code", - "execution_count": 42, + "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { - "name": "stderr", + "name": "stdout", "output_type": "stream", "text": [ - "25-Oct 20:27:59:DEBUG:root:total dofs: 4\n" + "assembling problem for temperature\n" ] }, { - "ename": "LoadError", - "evalue": "LoadError: MethodError: `has_lhs` has no method matching has_lhs(::DC2D4)\nwhile loading In[42], in expression starting on line 37", - "output_type": "error", - "traceback": [ - "LoadError: MethodError: `has_lhs` has no method matching has_lhs(::DC2D4)\nwhile loading In[42], in expression starting on line 37", - "" + "data": { + "text/plain": [ + "Success :: (line:-1) :: fact was true\n", + " Expression: T --> roughly(100.0)\n", + " Expected: 100.0\n", + " Occurred: 100.00000000000003" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dimension of unknown field: 1, problem dofs: 4\n", + "assembling problem for reaction force\n", + "dimension of unknown field: 1, problem dofs: 4\n", + "Residual norm: 9.845568954283847e-14\n", + "Temperature at point X = [0.5,0.0] is T = 100.00000000000003\n" ] } ], "source": [ - "# Define Problem 1:\n", - "# - Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]}\n", - "# - Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁\n", + "\"\"\" Define Problem 1:\n", "\n", - "el1 = Quad4([1, 2, 3, 4])\n", - "push!(el1, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]))\n", - "push!(el1, FieldSet(\"temperature thermal conductivity\", [Field(0.0, 6.0)]))\n", + "- Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]}\n", + "- Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁\n", + "\"\"\"\n", + "function get_heatproblem()\n", + " el1 = Quad4([1, 2, 3, 4])\n", + " push!(el1, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]))\n", + " push!(el1, FieldSet(\"temperature thermal conductivity\", [Field(0.0, 6.0)]))\n", + " push!(el1, FieldSet(\"temperature load\", [Field(0.0, 0.0*[12.0, 12.0, 12.0, 12.0])]))\n", + " push!(el1, FieldSet(\"density\", [Field(0.0, 10.0)]))\n", "\n", - "el2 = Seg2([1, 2])\n", - "push!(el2, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]])]))\n", + " el2 = Seg2([1, 2])\n", + " push!(el2, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]])]))\n", "\n", - "# Boundary load, linear ramp 0 -> 600 at time 0 -> 1\n", - "load = FieldSet(\"temperature flux\")\n", - "push!(load, Field(0.0, 0.0))\n", - "push!(load, Field(1.0, 600.0))\n", - "push!(el2, load)\n", + " # Boundary load, linear ramp 0 -> 600 at time 0 -> 1\n", + " load = FieldSet(\"temperature flux\")\n", + " push!(load, Field(0.0, 0.0))\n", + " push!(load, Field(1.0, 600.0))\n", + " push!(el2, load)\n", "\n", - "problem1 = PlaneHeatProblem()\n", - "push!(problem1, el1)\n", - "push!(problem1, el2)\n", + " problem1 = PlaneHeatProblem()\n", + " push!(problem1, el1)\n", + " push!(problem1, el2)\n", + " return problem1\n", + "end\n", "\n", - "# Define Problem 2:\n", - "# - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂\n", + "\"\"\" Define Problem 2:\n", + " - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂\n", + "\"\"\"\n", + "function get_boundaryproblem()\n", + " el3 = Seg2([3, 4])\n", + " push!(el3, FieldSet(\"geometry\", [Field(0.0, Vector[[1.0, 1.0], [0.0, 1.0]])]))\n", "\n", - "el3 = Seg2([3, 4])\n", - "push!(el3, FieldSet(\"geometry\", [Field(0.0, Vector[[1.0, 1.0], [0.0, 1.0]])]))\n", + " problem2 = DirichletProblem(1)\n", + " push!(problem2, el3)\n", + " return problem2\n", + "end\n", "\n", - "problem2 = DirichletProblem()\n", - "push!(problem2, el3)\n", + "function main()\n", + " problem1 = get_heatproblem()\n", + " problem2 = get_boundaryproblem()\n", + " # Create a solver for a set of problems\n", + " solver = SimpleSolver()\n", + " push!(solver, problem1)\n", + " push!(solver, problem2)\n", + " # Solve problem at time t=1.0 and update fields\n", + " call(solver, 1.0)\n", "\n", - "# Create a solver for a set of problems\n", - "solver = SimpleSolver()\n", - "push!(solver, problem1)\n", - "push!(solver, problem2)\n", + " # Postprocess.\n", + " # Interpolate temperature field along boundary of Γ₁ at time t=1.0\n", + " xi = [0.0, -1.0]\n", + " el2 = get_element(problem1.equations[2])\n", + " basis = get_basis(el2)\n", + " X = basis(\"geometry\", xi, 1.0)\n", + " T = basis(\"temperature\", xi, 1.0)\n", + " println(\"Temperature at point X = $X is T = $T\")\n", + " @fact T --> roughly(100.0)\n", + "end\n", "\n", - "# Solve problem at time t=1.0 and update fields\n", - "call(solver, 1.0)\n", - "\n", - "# Postprocess.\n", - "# Interpolate temperature field along boundary of Γ₁ at time t=1.0\n", - "xi = linspace([-1.0], [1.0], 5)\n", - "X = interpolate(el2, \"geometry\", xi, 1.0)\n", - "T = interpolate(el2, \"temperature\", xi, 1.0)\n", - "println(X)\n", - "println(T)" - ] - }, - { - "cell_type": "code", - "execution_count": 43, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "Error :: (line:-1)\n", - " Expression: mean(T) --> roughly(100.0)\n", - " UndefVarError: T not defined\n", - " in anonymous at /home/jukka/.julia/v0.4/FactCheck/src/FactCheck.jl:271\n", - " in do_fact at /home/jukka/.julia/v0.4/FactCheck/src/FactCheck.jl:333\n", - " in include_string at loading.jl:266\n", - " in execute_request_0x535c5df2 at /home/jukka/.julia/v0.4/IJulia/src/execute_request.jl:177\n", - " in eventloop at /home/jukka/.julia/v0.4/IJulia/src/IJulia.jl:141\n", - " in anonymous at task.jl:447" - ] - }, - "execution_count": 43, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "@fact mean(T) --> roughly(100.0)" + "main()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## TODO\n", - "\n", - "- dynamics\n", - "- different basis for trial and test (Petrov-Galerkin)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Frequently asked questions\n", - "\n", - "Any question about data structures, better ideas and improvements are very welcome." + "In this notebook the basic instructions how to develop JuliaFEM has given. The most imporant concepts has been considered; how to develop own element with own basis, several ways how to define own equation, and how to finally assemble and calculate the problem using solver. Any comments and/or discussion about technical details, theory, programming, or from life in general is very desirable; our issue log is in address https://github.com/JuliaFEM/JuliaFEM.jl/issues" ] } ], diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index 90a1b4d..9c5beba 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -573,7 +573,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 29, "metadata": { "collapsed": false }, @@ -582,30 +582,46 @@ "name": "stderr", "output_type": "stream", "text": [ - "24-Aug 01:07:34:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n", - "24-Aug 01:07:34:DEBUG:root:Found NODE section\n", - "24-Aug 01:07:35:DEBUG:root:Found ELEMENT section\n", - "24-Aug 01:07:35:DEBUG:root:120 elements found\n", - "24-Aug 01:07:35:INFO:root:Creating ELSET Body1\n", - "24-Aug 01:07:35:DEBUG:root:Found NSET section\n", - "24-Aug 01:07:35:DEBUG:root:Creating node set SUPPORT\n", - "24-Aug 01:07:35:DEBUG:root:Found NSET section\n", - "24-Aug 01:07:35:DEBUG:root:Creating node set LOAD\n", - "24-Aug 01:07:35:DEBUG:root:Found NSET section\n", - "24-Aug 01:07:35:DEBUG:root:Creating node set TOP\n" + "26-Oct 06:27:30:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n", + "26-Oct 06:27:30:DEBUG:root:processing -- section\n", + "26-Oct 06:27:30:DEBUG:root:Found NODE section\n", + "26-Oct 06:27:30:DEBUG:root:processing ++ section\n", + "26-Oct 06:27:30:DEBUG:root:processing -- section\n", + "26-Oct 06:27:30:DEBUG:root:Starting to process\n", + "26-Oct 06:27:30:DEBUG:root:Found ELEMENT section\n", + "26-Oct 06:27:30:DEBUG:root:processing ++ section\n", + "26-Oct 06:27:30:DEBUG:root:processing -- section\n", + "26-Oct 06:27:30:DEBUG:root:Starting to process\n", + "26-Oct 06:27:30:DEBUG:root:Parsing elements\n", + "26-Oct 06:27:30:DEBUG:root:120 elements found\n", + "26-Oct 06:27:30:INFO:root:Creating ELSET Body1\n", + "26-Oct 06:27:30:DEBUG:root:Found NSET section\n", + "26-Oct 06:27:30:DEBUG:root:processing ++ section\n", + "26-Oct 06:27:30:DEBUG:root:processing -- section\n", + "26-Oct 06:27:30:DEBUG:root:Starting to process\n", + "26-Oct 06:27:30:DEBUG:root:Creating node set SUPPORT\n", + "26-Oct 06:27:30:DEBUG:root:Found NSET section\n", + "26-Oct 06:27:30:DEBUG:root:processing ++ section\n", + "26-Oct 06:27:30:DEBUG:root:processing -- section\n", + "26-Oct 06:27:30:DEBUG:root:Starting to process\n", + "26-Oct 06:27:30:DEBUG:root:Creating node set LOAD\n", + "26-Oct 06:27:30:DEBUG:root:Found NSET section\n", + "26-Oct 06:27:30:DEBUG:root:processing ++ section\n", + "26-Oct 06:27:30:DEBUG:root:Starting to process\n", + "26-Oct 06:27:30:DEBUG:root:Creating node set TOP\n" ] }, { "data": { "text/plain": [ "Dict{Any,Any} with 4 entries:\n", - " \"nodes\" => Dict{Any,Any}(288=>[97.5,7.5,10.0],11=>[92.5,2.5,5.0],134=>[45.0,10.0,0.0],158=>[2.5,2.5,0.0],160=>[7.5,7.5,0.0],215=>[60.0,0.0,5.0],29=>[2.5,7…\n", - " \"elements\" => Dict{Any,Any}(68=>[71,144,149,198,51,150,57,43,50,214],2=>[204,199,175,130,207,208,209,3,4,176],89=>[95,78,104,52,127,126,106,60,68,67],11=>[15…\n", - " \"elsets\" => Dict{Any,Any}(\"Body1\"=>[1,2,3,4,5,6,7,8,9,10 … 111,112,113,114,115,116,117,118,119,120])\n", - " \"nsets\" => Dict{Any,Any}(\"LOAD\"=>[82,84,87,179,197,246,249,256,257],\"SUPPORT\"=>[108,109,111,155,162,216,225,281,298],\"TOP\"=>[70,75,76,84,88,90,95,96,98,10…" + " \"nodes\" => Dict{Any,Any}(288=>[97.5,7.5,10.0],11=>[92.5,2.5,5.0],134=>[45.…\n", + " \"elements\" => Dict{Any,Any}(68=>[71,144,149,198,51,150,57,43,50,214],2=>[204,…\n", + " \"elsets\" => Dict{Any,Any}(\"Body1\"=>[1,2,3,4,5,6,7,8,9,10 … 111,112,113,11…\n", + " \"nsets\" => Dict{Any,Any}(\"LOAD\"=>[82,84,87,179,197,246,249,256,257],\"SUPPO…" ] }, - "execution_count": 23, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } @@ -619,29 +635,7 @@ }, { "cell_type": "code", - "execution_count": 20, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "abstract Tet10 <: CG\n", - "\n", - "\"\"\"\n", - "Stress/displacement elements. 10-node quadratic tetrahedron.\n", - "\"\"\"\n", - "type C3D10 <: Tet10\n", - " id :: Int\n", - " node_ids :: Array{Int, 1}\n", - " coordinates :: Array{Float64, 2}\n", - " integration_points :: Array{IntegrationPoint, 1}\n", - " attributes :: Dict{ASCIIString, Any}\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 21, + "execution_count": 31, "metadata": { "collapsed": false }, @@ -649,123 +643,35 @@ { "data": { "text/plain": [ - "C3D10" + "size (generic function with 65 methods)" ] }, - "execution_count": 21, + "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function C3D10(element_id, node_ids)\n", - "\n", - " # Construct Lagrange basis\n", - " \n", - " P(xi) = [\n", - " 1\n", - " xi[1]\n", - " xi[2]\n", - " xi[3]\n", - " xi[1]^2\n", - " xi[2]^2\n", - " xi[3]^2\n", - " xi[1]*xi[2]\n", - " xi[2]*xi[3]\n", - " xi[3]*xi[1]]\n", - "\n", - " dP(xi) = [\n", - " 0 0 0\n", - " 1 0 0\n", - " 0 1 0\n", - " 0 0 1\n", - " 2*xi[1] 0 0\n", - " 0 2*xi[2] 0\n", - " 0 0 2*xi[3]\n", - " xi[2] xi[1] 0\n", - " 0 xi[3] xi[2]\n", - " xi[3] 0 xi[1]\n", - " ]\n", - "\n", - " X = [\n", - " 0.0 0.0 0.0\n", - " 1.0 0.0 0.0\n", - " 0.0 1.0 0.0\n", - " 0.0 0.0 1.0\n", - " 0.5 0.0 0.0\n", - " 0.5 0.5 0.0\n", - " 0.0 0.5 0.0\n", - " 0.0 0.0 0.5\n", - " 0.5 0.0 0.5\n", - " 0.0 0.5 0.5]\n", - "\n", - " A = zeros(10, 10)\n", - "\n", - " for i=1:10\n", - " A[i,:] = P(X[i,:])\n", - " end\n", - "\n", - " invA = inv(A)\n", - " basis(xi) = invA'*P(xi)\n", - " dbasis(xi) = invA'*dP(xi)\n", - "\n", - " shape_functions = FunctionSpace(basis, dbasis)\n", - "\n", + "using JuliaFEM: Tet10\n", + "\"\"\" Stress/displacement elements. 10-node quadratic tetrahedron. \"\"\"\n", + "type C3D10 <: Elasticity\n", + " element :: Tet10\n", + " integration_points :: Array{IntegrationPoint, 1}\n", + " global_dofs :: Array{Int64, 1}\n", + "end\n", + "function C3D10(element::Tet10)\n", " a = .585410196624969\n", " b = .138196601125011\n", " w = .041666666666667\n", - "\n", " integration_points = [\n", - " JuliaFEM.IntegrationPoint([a, b, b], w, Dict()),\n", - " JuliaFEM.IntegrationPoint([b, a, b], w, Dict()),\n", - " JuliaFEM.IntegrationPoint([b, b, a], w, Dict()),\n", - " JuliaFEM.IntegrationPoint([b, b, b], w, Dict())]\n", - " \n", - " attributes = Dict(\"displacement\" => zeros(3, 10))\n", - "\n", - " C3D10(element_id, node_ids, shape_functions, integration_points, attributes)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 22, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "get_lhs (generic function with 4 methods)" - ] - }, - "execution_count": 22, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "\"\"\"\n", - "1-node point force element for 3d elasticity problem.\n", - "\"\"\"\n", - "type C3D1 <: ContinuumElement\n", - " id :: Int\n", - " node_ids :: Array{Int, 1}\n", - " attributes :: Dict{ASCIIString, Any}\n", + " JuliaFEM.IntegrationPoint([a, b, b], w),\n", + " JuliaFEM.IntegrationPoint([b, a, b], w),\n", + " JuliaFEM.IntegrationPoint([b, b, a], w),\n", + " JuliaFEM.IntegrationPoint([b, b, b], w)]\n", + " push!(element, FieldSet(\"displacement\"))\n", + " C3D10(element, integration_points, [])\n", "end\n", - "function C3D1(element_id, node_ids)\n", - " attributes = Dict(\n", - " \"displacement\" => zeros(3, 1),\n", - " \"displacement nodal load\" => zeros(3, 1))\n", - " C3D1(element_id, node_ids, attributes)\n", - "end\n", - "function get_rhs(el::C3D1)\n", - " return el.attributes[\"displacement nodal load\"]\n", - "end\n", - "function get_lhs(el::C3D1)\n", - " return None\n", - "end" + "JuliaFEM.size(equation::C3D10) = 30" ] }, { @@ -860,7 +766,7 @@ } ], "source": [ - "function solve_3d_model()\n", + "function create_model()\n", " Logging.debug(\"Creating elements\")\n", " elements = JuliaFEM.Element[]\n", " \n", diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 2314de1..f7f7b67 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -36,12 +36,12 @@ include("interpolate.jl") # interpolation routines include("elements.jl") include("lagrange.jl") # Lagrange elements #include("hierarchical.jl") # P-elements +include("integrate.jl") # integration points include("equations.jl") include("problems.jl") include("solvers.jl") -#include("math.jl") # basic mathematical operations -- obsolete ..? # pre- and postprocess include("xdmf.jl") include("abaqus_reader.jl") diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 35c36ae..b594331 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -5,48 +5,56 @@ abstract DirichletEquation <: Equation -get_unknown_field_name(eq::DirichletEquation) = symbol("reaction force") - ### Dirichlet problem + equations type DirichletProblem <: BoundaryProblem + unknown_field_name :: ASCIIString + unknown_field_dimension :: Int equations :: Array{DirichletEquation, 1} + element_mapping :: Dict{DataType, DataType} + field_value :: Function end -function DirichletProblem() - DirichletProblem([]) -end -get_dimension(pr::Type{DirichletProblem}) = 1 # ..? -get_equation(pr::Type{DirichletProblem}, el::Type{Seg2}) = DBC2D2 -""" -Dirichlet boundary condition element for 2 node line segment -""" +function DirichletProblem(dimension::Int, field_value::Function=(X)->[0.0,0.0,0.0]) + element_mapping = nothing + if dimension == 1 + element_mapping = Dict( + Seg2 => DBC2D2 + ) + end + DirichletProblem("reaction force", dimension, [], element_mapping, field_value) +end + +""" Dirichlet boundary condition element for 2 node line segment """ type DBC2D2 <: DirichletEquation element :: Seg2 integration_points :: Array{IntegrationPoint, 1} - global_dofs :: Array{Int64, 1} - fieldval :: Function end function DBC2D2(element::Seg2) integration_points = [ IntegrationPoint([-sqrt(1/3)], 1.0), IntegrationPoint([+sqrt(1/3)], 1.0)] push!(element, FieldSet("reaction force")) - fieldval(X, t) = 0.0 - DBC2D2(element, integration_points, [], fieldval) + DBC2D2(element, integration_points) +end +Base.size(equation::DBC2D2) = (1, 2) + +function calculate_local_assembly!(assembly::LocalAssembly, equation::DirichletEquation, + unknown_field_name::ASCIIString, time::Number=Inf, + problem=nothing) + initialize_local_assembly!(assembly, equation) + element = get_element(equation) + basis = get_basis(element) + detJ = det(basis) + for ip in get_integration_points(equation) + w = ip.weight * detJ(ip) + N = basis(ip, time) + assembly.stiffness_matrix += w * N'*N + if !isa(problem, Void) + X = basis("geometry", ip, time) + u = problem.field_value(X) + assembly.force_vector += w * N'*u + end + end end -function get_lhs(eq::DBC2D2, ip, t) - el = get_element(eq) - h = get_basis(el)(ip.xi) - return h*h' -end -function get_rhs(eq::DBC2D2, ip, t) - el = get_element(eq) - h = get_basis(el, ip.xi) - f = eq.fieldval - X = interpolate(el, "geometry", ip.xi, t) - return h*f(X, t) -end -has_lhs(eq::DBC2D2) = true -has_rhs(eq::DBC2D2) = true diff --git a/src/elements.jl b/src/elements.jl index d4fe7b9..af2d8dc 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -15,12 +15,12 @@ abstract Element """ Get FieldSet from element. """ function Base.getindex(element::Element, field_name) - element.fields[symbol(field_name)] + element.fields[field_name] end """ Add new FieldSet to element. """ function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name) - fieldset.name = symbol(fieldset_name) + fieldset.name = fieldset_name element.fields[fieldset.name] = fieldset end function Base.push!(element::Element, fieldset::FieldSet) @@ -73,9 +73,8 @@ End of example. =# -# These must be implemented for your own element -get_number_of_basis_functions(el::Type{Element}) = nothing -get_element_dimension(el::Type{Element}) = nothing +# define size of your element as (dim, nbasis) tuple where first integer is spatial dimension and second is number of basis functions. +# Base.size(element::Type{Element}) = nothing ### COMMON ELEMENT ROUTINES ### @@ -95,12 +94,14 @@ This uses FactCheck and throws exceptions if element is not passing all tests. function test_element(element_type) Logging.info("Testing element $element_type") local element - n = get_number_of_basis_functions(element_type) - Logging.info("number of basis functions in this element: $n") - @fact n --> not(nothing) """ - Unable to determine number of nodes for $eltype define a function - 'get_number_of_basis_functions' which returns the number of nodes - for this element.""" + dim = nothing + n = nothing + try + dim, n = size(element_type) + catch + Logging.error("Unable to determine element dimensions. Define Base.size(element::Type{$elementtype}) = (dim, nbasis) where dim is spatial dimension of element and nbasis is number of basis functions of element.") + end + Logging.info("element dimension: $dim x $n") Logging.info("Initializing element") try @@ -112,45 +113,23 @@ function test_element(element_type) return false end - dim = get_element_dimension(element_type) - Logging.info("Element dimension: $dim") - @fact dim --> not(nothing) """ - Unable to get element dimension define function 'get_element_dimension' - which return the dimension of this element (1, 2, 3)""" - # try to interpolate some scalar field - field = Field(0.0, collect(1:n)) - Logging.info("Creating new scalar field $field") - fieldset = FieldSet("field1") - push!(fieldset, field) - push!(element, fieldset) + push!(element, FieldSet("field1", [Field(0.0, collect(1:n))])) + # TODO: how to parametrize this? push!(element, FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])) # evaluate basis functions at middle point of element + basis = get_basis(element) + dbasis = grad(basis) mid = zeros(dim) - try - basis = get_basis(element) - val1 = basis(mid, 0.0) - Logging.info("basis at $mid: $val1") - val2 = basis("field1", mid, 0.0) - Logging.info("field val at $mid: $val2") - catch - Logging.error(""" - Unable to evaluate basis, define function 'get_basis' for - this element.""") - end - try - basis = get_basis(element) - dbasis = grad(basis) - val3 = dbasis(mid, 0.0) - Logging.info("derivative of basis at $mid: $val3") - val4 = dbasis("field1", mid, 0.0) - Logging.info("field val at $mid: $val4") - catch - Logging.error(""" - Unable to evaluate partial derivatives of basis, - define function 'get_dbasisdxi' for this element.""") - end + val1 = basis(mid, 0.0) + Logging.info("basis at $mid: $val1") + val2 = basis("field1", mid, 0.0) + Logging.info("field val at $mid: $val2") + val3 = dbasis(mid, 0.0) + Logging.info("derivative of basis at $mid: $val3") + val4 = dbasis("field1", mid, 0.0) + Logging.info("field val at $mid: $val4") Logging.info("Element $element_type passed tests.") end @@ -184,73 +163,70 @@ function call(u::FunctionSpace, field_name, xi::Vector, t::Number=Inf, variation return f.values end h = u.element.basis.basis(xi) - return h*f + return dot(vec(h), f) end """ If basis is called without a field, return basis functions evaluated at that point. """ function call(u::FunctionSpace, xi::Vector, t::Number=Inf) - return u.element.basis.basis(xi)' + return u.element.basis.basis(xi) end """ Evaluate gradient of field on element function space. """ function call(gradu::GradientFunctionSpace, field_name, xi::Vector, t::Number=Inf, variation=nothing) f = !isa(variation, Void) ? variation : gradu.element[field_name](t) X = gradu.element["geometry"](t) - b = gradu.element.basis.dbasisdxi(xi) - return b*f*inv(b*X) + dN = gradu.element.basis.dbasisdxi(xi) + J = sum([dN[:,i]*X[i]' for i=1:length(X)]) + grad = inv(J)*dN + gradf = sum([grad[:,i]*f[i]' for i=1:length(f)])' + return gradf end """ If gradient of basis is called without a field, return "empty" gradient evaluated at that point. """ function call(gradu::GradientFunctionSpace, xi::Vector, t::Number=Inf) X = gradu.element["geometry"](t) - b = gradu.element.basis.dbasisdxi(xi) - return (b*inv(b*X))' + dN = gradu.element.basis.dbasisdxi(xi) + J = sum([dN[:,i]*X[i]' for i=1:length(X)]) + grad = inv(J)*dN + return grad end # on-line functions to get api more easy to use, ip -> xi.ip -call(u::FunctionSpace, ip::IntegrationPoint, t::Number) = call(u, ip.xi, t) -call(u::FunctionSpace, ip::IntegrationPoint) = call(u, ip.xi) -call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number) = call(u, ip.xi, t) -call(u::GradientFunctionSpace, ip::IntegrationPoint) = call(u, ip.xi) +call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t) +call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t) +# i think these will be the most called functions. +call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation) +call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation) +call(u::FunctionSpace, field_name) = (args...) -> call(u, field_name, args...) +call(u::GradientFunctionSpace, field_name) = (args...) -> call(u, field_name, args...) -""" Return field from function space. """ +""" Return a field from function space. """ function get_field(u::FunctionSpace, field_name, time=Inf) return u.element[field_name](time) end -""" Return field from function space. """ +""" Return a field from function space. """ function get_field(u::FunctionSpace, field_name, time=Inf, variation=nothing) return !isa(variation, Void) ? variation : u.element[field_name](time) end -""" Return fieldset from function space. """ +""" Return a fieldset from function space. """ function get_fieldset(u::FunctionSpace, field_name) return u.element[field_name] end -# i think these will be the most called functions. -call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation) -call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation) - -function jacobian(u::FunctionSpace, xi, t) - u.element.basis.dbasisdxi(xi)*u.element["geometry"](t) +function LinAlg.det(u::FunctionSpace, xi::Vector, t::Number=Inf) + X = u.element["geometry"](t) + dN = u.element.basis.dbasisdxi(xi) + J = sum([dN[:,i]*X[i]' for i=1:length(X)]) + m, n = size(J) + return m == n ? det(J) : norm(J) end - -function jacobian(u::FunctionSpace, ip::IntegrationPoint, t::Number) - jacobian(u, ip.xi, t) +function LinAlg.det(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) + LinAlg.det(u, ip.xi, t) end - -function jacobian(u::FunctionSpace, xi) - jacobian(u, xi, Inf) -end - function LinAlg.det(u::FunctionSpace) - function detJ(args...) - J = jacobian(u, args...) - m, n = size(J) - return m == n ? det(J) : norm(J) - end - return detJ + return (args...) -> det(u, args...) end function get_basis(element::Element) @@ -265,7 +241,7 @@ Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(a """ Check does fieldset exist. """ function Base.haskey(element::Element, what) - haskey(element.fields, symbol(what)) + haskey(element.fields, what) end diff --git a/src/equations.jl b/src/equations.jl index 5ee0de6..16ca2da 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -11,17 +11,12 @@ type LocalAssembly <: Assembly mass_matrix :: Matrix stiffness_matrix :: Matrix force_vector :: Matrix - potential_energy# :: Union{Array, Float64} + potential_energy residual_vector :: Vector end -function LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector::Matrix) - LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector[:]) -end - -""" Initialize workspace for local assembly. """ -function LocalAssembly(equation::Equation) - ndofs = size(equation) +""" Initialize workspace for local matrices for dimension ndofs. """ +function initialize_local_assembly(ndofs::Int=1) mass_matrix = zeros(ndofs, ndofs) stiffness_matrix = zeros(ndofs, ndofs) force_vector = zeros(ndofs, 1) @@ -31,14 +26,24 @@ function LocalAssembly(equation::Equation) potential_energy, residual_vector) end +""" Initialize workspace for local matrices, get dimension from equation. """ function initialize_local_assembly(equation::Equation) - LocalAssembly(equation) + ndofs = prod(size(equation)) + return initialize_local_assembly(ndofs) end -function initialize_local_assembly(equation::Equation, assembly::LocalAssembly) - if size(equation) != assembly.ndofs +""" Initialize or zero workspace. """ +function initialize_local_assembly!(assembly::LocalAssembly, equation::Equation) + ndofs = prod(size(equation)) + if ndofs != assembly.ndofs # if problem size changes, automatically initialize new work space - return initialize_local_assembly(equation) + assembly.ndofs = ndofs + assembly.mass_matrix = zeros(ndofs, ndofs) + assembly.stiffness_matrix = zeros(ndofs, ndofs) + assembly.force_vector = zeros(ndofs, 1) + assembly.potential_energy = 0.0 + assembly.residual_vector = zeros(ndofs) + return end # otherwise, empty workspace ready for next iteration fill!(assembly.mass_matrix, 0.0) @@ -46,15 +51,7 @@ function initialize_local_assembly(equation::Equation, assembly::LocalAssembly) fill!(assembly.force_vector, 0.0) assembly.potential_energy = 0.0 fill!(assembly.residual_vector, 0.0) - return assembly -end -function initialize_local_assembly(assembly::LocalAssembly, equation::Equation) - initialize_local_assembly(equation, assembly) -end - -function get_unknown_field_name(equation::Equation) - eqtype = typeof(equation) - error("define get_unknown_field_name for this equation type $eqtype") + return end has_mass_matrix(equation::Equation) = false @@ -73,14 +70,15 @@ get_integration_points(equation::Equation) = equation.integration_points """ Return a local assembly for element. """ -function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, time::Number=Inf) +function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, + unknown_field_name::ASCIIString, time::Number=Inf, + problem=nothing) - initialize_local_assembly(assembly, equation) # zero all + initialize_local_assembly!(assembly, equation) # zero all element = get_element(equation) basis = get_basis(element) detJ = det(basis) - field_name = get_unknown_field_name(equation) # 1. if equations are defined we just integrate them if has_mass_matrix(equation) || has_stiffness_matrix(equation) || has_force_vector(equation) @@ -96,17 +94,16 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, assembly.force_vector += s*get_force_vector(equation, ip, time)[:] end # external loads -- if any nodal loads is defined add to force vector - if haskey(element, "$field_name nodal load") - assembly.force_vector += element["$field_name nodal load"](time)[:] + if haskey(element, "$unknown_field_name nodal load") + assembly.force_vector += element["$unknown_field_name nodal load"](time)[:] end end end # 2. variational / energy form - user has defined some potential energy / variational form if has_potential_energy(equation) - field_name = get_unknown_field_name(equation) element = get_element(equation) - field = element[field_name](time) + field = element[unknown_field_name](time) function potential_energy(data::Vector) # calculate potential energy for some setting. this is needed by forwarddiff assembly.potential_energy = 0.0 @@ -117,8 +114,8 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, assembly.potential_energy += ip.weight * dw * detJ(ip) end # external energy -- if any nodal loads is defined, decrease from potential energy - if haskey(element, "$field_name nodal load") - P = element["$field_name nodal load"](time) + if haskey(element, "$unknown_field_name nodal load") + P = element["$unknown_field_name nodal load"](time) assembly.potential_energy -= dot(P[:], df[:]) end if isa(assembly.potential_energy, Array) @@ -135,9 +132,8 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, # 3. virtual work form - user has defined residual vector δW_int(u,δu) + δW_ext(u,δu) = 0 ∀ v if has_residual_vector(equation) - field_name = get_unknown_field_name(equation) element = get_element(equation) - field = element[field_name](time) + field = element[unknown_field_name](time) function residual_vector(data::Vector) fill!(assembly.residual_vector, 0.0) df = similar(field, data) @@ -147,8 +143,8 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, assembly.residual_vector += ip.weight*dr*detJ(ip) end # external loads -- if any nodal loads is defined, remove from residual - if haskey(element, "$field_name nodal load") - assembly.residual_vector -= element["$field_name nodal load"](time)[:] + if haskey(element, "$unknown_field_name nodal load") + assembly.residual_vector -= element["$unknown_field_name nodal load"](time)[:] end return assembly.residual_vector end @@ -160,18 +156,3 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, end -function calculate_local_assembly!(equation::Equation, assembly::LocalAssembly, time::Number=Inf) - calculate_local_assembly!(assembly, equation) -end - - -""" Get global degrees of freedom for this element. """ -function get_global_dofs(eq::Equation) - eq.global_dofs -end - -""" Set global degrees of freedom for this element. """ -function set_global_dofs!(eq::Equation, dofs) - eq.global_dofs = dofs -end - diff --git a/src/heat.jl b/src/heat.jl index 1007d89..8808d96 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -6,68 +6,72 @@ abstract HeatProblem <: Problem abstract HeatEquation <: Equation -get_unknown_field_name(eq::HeatEquation) = symbol("temperature") - - ### Plane heat problem + equations ### type PlaneHeatProblem <: HeatProblem + unknown_field_name :: ASCIIString + unknown_field_dimension :: Int equations :: Array{HeatEquation, 1} + element_mapping :: Dict{DataType, DataType} end """ Default constructor for problem takes no arguments. """ function PlaneHeatProblem() - return PlaneHeatProblem([]) + element_mapping = Dict( + Quad4 => DC2D4, + Seg2 => DC2D2) + return PlaneHeatProblem("temperature", 1, [], element_mapping) end -""" Return dimension of unknown field variable, temperature is scalar field. """ -get_dimension(pr::Type{PlaneHeatProblem}) = 1 - -""" Map Lagrange element Quad4 to equation DC2D4 """ -get_equation(pr::Type{PlaneHeatProblem}, el::Type{Quad4}) = DC2D4 - -""" Map Lagrange element Seg2 to equation DC2D2 """ -get_equation(pr::Type{PlaneHeatProblem}, el::Type{Seg2}) = DC2D2 """ Diffusive heat transfer for 4-node bilinear element. """ type DC2D4 <: HeatEquation element :: Quad4 integration_points :: Array{IntegrationPoint, 1} - global_dofs :: Array{Int64, 1} end function DC2D4(element::Quad4) - integration_points = [ - IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0), - IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0), - IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0), - IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)] + integration_points = get_default_integration_points(element) push!(element, FieldSet("temperature")) - DC2D4(element, integration_points, []) + DC2D4(element, integration_points) end -function get_lhs(equation::DC2D4, ip, time) - element = get_element(equation) - dNdX = get_dbasisdX(element, ip.xi, time) - k = interpolate(element, "temperature thermal conductivity", ip.xi, time) - return dNdX*k*dNdX' -end -JuliaFEM.has_lhs(eq::DC2D4) = true +Base.size(equation::DC2D4) = (1, 4) """ Diffusive heat transfer for 2-node linear segment. """ type DC2D2 <: HeatEquation element :: Seg2 integration_points :: Array{IntegrationPoint, 1} - global_dofs :: Array{Int64, 1} end function DC2D2(element::Seg2) - integration_points = [IntegrationPoint([0.0], 2.0)] + integration_points = get_default_integration_points(element) push!(element, FieldSet("temperature")) - DC2D2(element, integration_points, []) + DC2D2(element, integration_points) end -function get_rhs(equation::DC2D2, ip, time) +Base.size(equation::DC2D2) = (1, 2) + +function calculate_local_assembly!(assembly::LocalAssembly, equation::HeatEquation, + unknown_field_name::ASCIIString, time::Number=Inf, + problem=nothing) + + initialize_local_assembly!(assembly, equation) + element = get_element(equation) - h = get_basis(element, ip.xi) - f = interpolate(element, "temperature flux", ip.xi, time) - return h*f + basis = get_basis(element) + dbasis = grad(basis) + detJ = det(basis) + for ip in get_integration_points(equation) + w = ip.weight * detJ(ip) + # evaluate fields in integration point + ρ = basis("density", ip, time) + k = basis("temperature thermal conductivity", ip, time) + f = basis("temperature load", ip, time) + # evaluate basis functions and gradient in integration point + N = basis(ip, time) + dN = dbasis(ip, time) + # do assembly + assembly.mass_matrix += w * ρ*N'*N + assembly.stiffness_matrix += w * k*dN'*dN + assembly.force_vector += w * N'*f + end end -JuliaFEM.has_rhs(eq::DC2D2) = true + diff --git a/src/integrate.jl b/src/integrate.jl new file mode 100644 index 0000000..9571a20 --- /dev/null +++ b/src/integrate.jl @@ -0,0 +1,17 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +function get_default_integration_points(element::Quad4) + [ + IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0), + IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0), + IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0), + IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0) + ] +end + +function get_default_integration_points(element::Seg2) + [ + IntegrationPoint([0.0], 2.0) + ] +end diff --git a/src/interpolate.jl b/src/interpolate.jl index c1dbfda..7910521 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -59,7 +59,7 @@ function interpolate(basis::Basis, field::Field, ip::IntegrationPoint) interpolate(basis, field, ip.xi) end -function dinterpolate(basis::Basis, u::Field, xi::Array{Float64, 1}) - basis.dbasisdxi(xi)*u -end +#function dinterpolate(basis::Basis, u::Field, xi::Array{Float64, 1}) +# basis.dbasisdxi(xi)*u +#end diff --git a/src/lagrange.jl b/src/lagrange.jl index 754d9c9..1707c42 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -17,8 +17,9 @@ function calculate_lagrange_basis(P, X) end # Logging.debug("Calculating inverse of A") invA = inv(A)' - basis(xi) = invA*P(xi) - basis + basis(xi) = (invA*P(xi))' + dbasisdxi(xi) = (ForwardDiff.jacobian((xi) -> invA*P(xi), xi, cache=autodiffcache))' + basis, dbasisdxi end """ @@ -33,21 +34,22 @@ macro create_lagrange_element(element_name, element_description, X, P) eltype = esc(element_name) quote global get_element_description - global get_number_of_basis_functions, get_element_dimension - dim = size($X, 1) - nbasis = size($X, 2) - basis = calculate_lagrange_basis($P, $X) + #global get_number_of_basis_functions, get_element_dimension + #dim = size($X, 1) + #nbasis = size($X, 2) + basis, dbasisdxi = calculate_lagrange_basis($P, $X) type $eltype <: CG connectivity :: Array{Int, 1} basis :: Basis - fields :: Dict{Symbol, FieldSet} + fields :: Dict{ASCIIString, FieldSet} end function $eltype(connectivity, args...) - $eltype(connectivity, Basis(basis), Dict()) + $eltype(connectivity, Basis(basis, dbasisdxi), Dict()) end get_element_description(el::Type{$eltype}) = $element_description - get_number_of_basis_functions(el::Type{$eltype}) = nbasis - get_element_dimension(el::Type{$eltype}) = dim + #get_number_of_basis_functions(el::Type{$eltype}) = nbasis + #get_element_dimension(el::Type{$eltype}) = dim + Base.size(el::Type{$eltype}) = Base.size($X) end end diff --git a/src/math.jl b/src/math.jl deleted file mode 100644 index c11c414..0000000 --- a/src/math.jl +++ /dev/null @@ -1,115 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using ForwardDiff - -""" -Linearize function f w.r.t some given field, i.e. calculate dR/du - -Parameters ----------- -f::Function - (possibly) nonlinear function to linearize -field::ASCIIString - field variable - -Returns -------- -Array{Float64, 2} - jacobian / "tangent stiffness matrix" -""" -function linearize(f::Function, el::Element, field::ASCIIString) - dim, nnodes = size(el.attributes[field]) - function helper!(x, y) - orig = copy(el.attributes[field]) - el.attributes[field] = reshape(x, dim, nnodes) - y[:] = f(el) - el.attributes[field] = copy(orig) - end - jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes) - return jac(el.attributes[field][:]) -end - - -""" -This version returns another function which can be then evaluated against field -""" -function linearize(f::Function, field::ASCIIString) - function jacobian(el::Element, args...) - fld = get_field(el, field) - dim, nnodes = size(fld) - function helper!(x, y) - orig = copy(fld) - set_field(el, field, reshape(x, dim, nnodes)) - y[:] = f(el, args...) - set_field(el, field, copy(orig)) - end - jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes) - return jac(fld[:]) - end - return jacobian -end - - -""" -In-place version, no additional garbage collection. -""" -function linearize!(f::Function, el::Element, field::ASCIIString, target::ASCIIString) - el.attributes[target][:] = 0.0 - dim, nnodes = size(el.attributes[field]) - function helper!(x, y) - orig = copy(el.attributes[field]) - el.attributes[field] = reshape(x, dim, nnodes) - y[:] = f(el) - el.attributes[field] = copy(orig) - end - jac! = ForwardDiff.forwarddiff_jacobian!(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes) - jac!(el.attributes[field][:], el.attributes[target]) -end - - -""" -This version returns a function which must be operated with element e -""" -function integrate(f::Function) - function integrate(el::Element) - target = [] - for ip in el.integration_points - J = interpolate(el, :geometry, ip.xi; derivative=true) - push!(target, ip.weight*f(el, ip)*det(J)) - end - return sum(target) - end - return integrate -end - -""" -This version saves results inplace to target, garbage collection free -""" -function integrate!(f::Function, el::Element, target) - # set target to zero - el.attributes[target][:] = 0.0 - for ip in el.integration_points - J = interpolate(el, :geometry, ip.xi; derivative=true) - el.attributes[target][:,:] += ip.weight*f(el, ip)*det(J) - end -end - -function linearize(eq::Equation, f::Function, field::ASCIIString) - function jacobian(eq::Equation, args...) - el = get_element(eq) - fld = get_field(el, field) - dim, nnodes = size(fld) - function helper(x::Vector) - orig = copy(fld) - set_field(el, field, reshape(x, dim, nnodes)) - y = f(eq, args...) - set_field(el, field, orig) - return y[:] - end - jac = ForwardDiff.jacobian(helper) - return jac(fld[:]) - end - return jacobian -end - diff --git a/src/problems.jl b/src/problems.jl index 0087934..69e4abe 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -5,141 +5,22 @@ abstract Problem abstract BoundaryProblem <: Problem abstract FieldProblem <: Problem -get_equations(pr::Problem) = pr.equations - -function get_dimension(pr::Type{Problem}) - throw("Unable to determine problem dimension for problem $pr") +function get_equations(problem::Problem) + problem.equations end -function get_equation(pr::Type{Problem}, el::Type{Element}) - throw("Could not find corresponding equation for element $el in problem $pr") +function get_unknown_field_dimension(problem::Problem) + problem.unknown_field_dimension end -""" -Add new element to problem -""" -function add_element!(problem::Problem, element::Element) - equation = get_equation(typeof(problem), typeof(element)) - push!(problem.equations, equation(element)) +function get_unknown_field_name(problem::Problem) + problem.unknown_field_name end + +""" Add new element to problem. """ function Base.push!(problem::Problem, element::Element) - equation = get_equation(typeof(problem), typeof(element)) - push!(problem.equations, equation(element)) -end - -""" -Return total number of basis functions in problem -""" -function get_number_of_basis_functions(pr::Problem) - conn = Int[] - for eq in get_equations(pr) - append!(conn, get_connectivity(eq)) - end - length(unique(conn)) -end - -""" -Problem matrix size dimension -""" -function get_matrix_dimension(pr::Problem) - get_dimension(typeof(pr))*get_number_of_basis_functions(pr) -end - -""" -Assign global dofs for element. This doesn't do any reordering. -""" -function set_global_dofs!(pr::Problem) - #ndim = get_dimension(pr)*get_number_of_basis_functions(pr) - #ndim = get_matrix_dimension(pr) - dim = get_dimension(typeof(pr)) - nconn = get_number_of_basis_functions(pr) - ndim = dim*nconn - Logging.debug("Problem (matrix) dimension: $ndim") - gdofs = reshape(collect(1:ndim), dim, nconn) - for eq in get_equations(pr) - lconn = get_connectivity(eq) - gconn = gdofs[:, lconn][:] - set_global_dofs!(eq, gconn) - end -end - -""" Return unique list of connectivity (i.e. node ids). """ -function get_connectivity(problem::Problem) - connectivity = Int[] - for equation in get_equations(problem) - element = get_element(equation) - append!(connectivity, get_connectivity(element)) - end - connectivity = unique(connectivity) - return connectivity -end - -""" -Calculate global dofs for equations, maybe using some bandwidth -minimizing or fill reducing algorithm -""" -function calculate_global_dofs(pr::Problem) - conn = get_connectivity(pr) - dim = get_dimension(typeof(pr)) - ndofs = dim*length(conn) - Logging.debug("total dofs: $ndofs") - - mconn = maximum(conn) - gdofs = reshape(collect(1:mconn), dim, mconn) - dofmap = Dict{Int64, Array{Int64, 1}}() - for (i, c) in enumerate(conn) - dofmap[c] = gdofs[:, i] - end - return dofmap -end - -""" -Assign global dofs for equations. -""" -function assign_global_dofs!(pr::Problem, dofmap) - for eq in get_equations(pr) - el = get_element(eq) - c = get_connectivity(el) - #gdofs = [dofmap[ci] for ci in c] - gdofs = Int64[] - for ci in c - append!(gdofs, dofmap[ci]) - end - set_global_dofs!(eq, gdofs) - end -end - -function get_lhs(pr::Problem, t::Float64) - I = Int64[] - J = Int64[] - V = Float64[] - dim = get_dimension(typeof(pr)) - for eq in filter(has_lhs, get_equations(pr)) - dofs = get_global_dofs(eq) - lhs = integrate_lhs(eq, t) - for (li, i) in enumerate(dofs) - for (lj, j) in enumerate(dofs) - push!(I, i) - push!(J, j) - push!(V, lhs[li, lj]) - end - end - end - return I, J, V -end - -function get_rhs(pr::Problem, t::Float64) - I = Int64[] - V = Float64[] - dim = get_dimension(typeof(pr)) - for eq in filter(has_rhs, get_equations(pr)) - dofs = get_global_dofs(eq) - rhs = integrate_rhs(eq, t) - for (li, i) in enumerate(dofs) - push!(I, i) - push!(V, rhs[li]) - end - end - return I, V + element_type = typeof(element) + equation_type = problem.element_mapping[element_type] + push!(problem.equations, equation_type(element)) end diff --git a/src/types.jl b/src/types.jl index f9330f5..5fcfbea 100644 --- a/src/types.jl +++ b/src/types.jl @@ -15,6 +15,9 @@ end function Field(time, values) Field(time, 0, values) end +function Field(values) + Field(0.0, 0, values) +end """ Get length of a field (number of basis functions in practice). """ function Base.length(f::Field) length(f.values) @@ -31,16 +34,23 @@ end function Base.(:*)(k::Number, f::Field) Field(f.time, k*f.values) end -""" Multiply field with some vector x. """ -function Base.(:*)(x::Vector, f::Field) + +""" Inner product of field and vector x. """ +function Base.dot(x::Vector, f::Field) @assert length(x) == length(f) sum([f[i]*x[i] for i in 1:length(f)]) end -""" Multiply field with some matrix x. """ -# function Base.(:*){T}(x::Matrix, f::Field{Vector{T}}) -function Base.(:*)(x::Matrix, f::Field) - sum([f[i]*x[i,:] for i in 1:length(f)]) + +function Base.size(field::Field) + (length(field.values[1]), length(field.values)) end + +#""" Multiply field with some matrix x. """ +# function Base.(:*){T}(x::Matrix, f::Field{Vector{T}}) +#function Base.(:*)(x::Matrix, f::Field) +# sum([f[i]*x[:,i]' for i in 1:length(f)]) +#end + """ Sum two fields. """ function Base.(:+)(f1::Field, f2::Field) @assert(f1.time == f2.time, "Cannot add fields: time mismatch, $(f1.time) != $(f2.time)") @@ -59,6 +69,9 @@ Examples function Base.getindex(field::Field, c::Colon) [field.values...;] end +function Base.vec(field::Field) + [field.values...;] +end """ Return field similar to input but with new data in it. @@ -87,20 +100,17 @@ end - - - """ FieldSet is set of fields, each field can have different time and/or increment. """ type FieldSet - name :: Symbol + name :: ASCIIString fields :: Array{Field, 1} end """ Initializer for FieldSet. """ function FieldSet(field_name) - FieldSet(Symbol(field_name), []) + FieldSet(field_name, []) end function FieldSet() - FieldSet(Symbol("unknown field"), []) + FieldSet("unknown field", []) end """ Add new field to fieldset. """ function Base.push!(fs::FieldSet, field::Field) @@ -127,15 +137,10 @@ type Basis basis :: Function dbasisdxi :: Function end -""" Constructor of basis function. """ -function Basis(basis) - Basis(basis, ForwardDiff.jacobian(basis)) -end -""" Get partial derivative of basis function. """ -function grad(basis::Basis) - (ip) -> basis.dbasisdxi(ip.xi) -end - +#""" Constructor of basis function. """ +#function Basis(basis) +# Basis(basis, ForwardDiff.jacobian(basis)) +#end """ Integration point @@ -151,7 +156,7 @@ attributes :: Dict{Any, Any} type IntegrationPoint xi :: Array{Float64, 1} weight :: Float64 - fields :: Dict{Symbol, FieldSet} + fields :: Dict{ASCIIString, FieldSet} end function IntegrationPoint(xi, weight) IntegrationPoint(xi, weight, Dict()) diff --git a/test/test_elements.jl b/test/test_elements.jl index c3b1f48..314e400 100644 --- a/test/test_elements.jl +++ b/test/test_elements.jl @@ -11,25 +11,20 @@ This should always pass test_element if everything is ok. type MockElement <: Element connectivity :: Array{Int, 1} basis :: Basis - fields :: Dict{Symbol, FieldSet} + fields :: Dict{ASCIIString, FieldSet} end function MockElement(connectivity) - h(xi) = [ - (1-xi[1])*(1-xi[2])/4 - (1+xi[1])*(1-xi[2])/4 - (1+xi[1])*(1+xi[2])/4 - (1-xi[1])*(1+xi[2])/4] - dh(xi) = [ - -(1-xi[2])/4.0 -(1-xi[1])/4.0 - (1-xi[2])/4.0 -(1+xi[1])/4.0 - (1+xi[2])/4.0 (1+xi[1])/4.0 - -(1+xi[2])/4.0 (1-xi[1])/4.0] + + h(xi) = 1/4*[(1-xi[1])*(1-xi[2]) (1+xi[1])*(1-xi[2]) (1+xi[1])*(1+xi[2]) (1-xi[1])*(1+xi[2])] + + dh(xi) = 1/4*[ + -(1-xi[2]) (1-xi[2]) (1+xi[2]) -(1+xi[2]) + -(1-xi[1]) -(1+xi[1]) (1+xi[1]) (1-xi[1])] + basis = Basis(h, dh) MockElement(connectivity, basis, Dict()) end -JuliaFEM.get_number_of_basis_functions(el::Type{MockElement}) = 4 -JuliaFEM.get_element_dimension(el::Type{MockElement}) = 2 - +Base.size(element::Type{MockElement}) = (2, 4) using JuliaFEM: test_element facts("test test_element against mock element") do diff --git a/test/test_types.jl b/test/test_types.jl index 1925d92..135a7ba 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -1,7 +1,7 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -using JuliaFEM: Basis, Field, FieldSet, interpolate, dinterpolate +using JuliaFEM: Basis, Field, FieldSet, interpolate using FactCheck facts("test fields") do @@ -19,6 +19,8 @@ facts("test fields") do @fact u3.values --> [1.0, 3.0] end +#= to be fixed + facts("test interpolation of fields") do # interpolation of field in spatial domain @@ -66,8 +68,6 @@ facts("test interpolation of fields") do # get midpoint of field in spatial domain @fact interpolate(N, X, [0.0, 0.0]) --> [0.5, 0.5] - # derivatives of field at midpoint - @fact dinterpolate(N, X, [0.0, 0.0]) --> [0.5 0.0; 0.0 0.5] # interpolate of scalar field -> scalar H = Field(0.0, 6.0) @@ -82,3 +82,6 @@ facts("test interpolation of fields") do #@fact interpolate(N, f, [0.0]) --> [50.0, 50.0] end + +=# +