diff --git a/docs/tutorials/2015-08-29-developing-juliafem.ipynb b/docs/tutorials/2015-08-29-developing-juliafem.ipynb index 56f9232..ab9fa46 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": 3, + "execution_count": 2, "metadata": { "collapsed": false }, @@ -60,7 +60,7 @@ "Logger(root,DEBUG,Base.PipeEndpoint(open, 0 bytes waiting),root)" ] }, - "execution_count": 3, + "execution_count": 2, "metadata": {}, "output_type": "execute_result" } @@ -89,13 +89,13 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ - "using JuliaFEM: Element, Field, FieldSet, Basis" + "using JuliaFEM: Element, Field, FieldSet, Basis, Quad4" ] }, { @@ -107,7 +107,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": { "collapsed": false }, @@ -129,7 +129,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -140,7 +140,7 @@ "MyQuad4" ] }, - "execution_count": 6, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } @@ -165,7 +165,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "metadata": { "collapsed": false }, @@ -173,10 +173,10 @@ { "data": { "text/plain": [ - "size (generic function with 72 methods)" + "size (generic function with 74 methods)" ] }, - "execution_count": 7, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } @@ -194,7 +194,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -207,7 +207,7 @@ " 2 4 6 8" ] }, - "execution_count": 8, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -225,7 +225,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -244,7 +244,7 @@ " 8" ] }, - "execution_count": 9, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } @@ -256,7 +256,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "metadata": { "collapsed": false }, @@ -269,7 +269,7 @@ " 2 4 6 8" ] }, - "execution_count": 10, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } @@ -287,7 +287,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": { "collapsed": false }, @@ -317,7 +317,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "metadata": { "collapsed": false }, @@ -326,15 +326,15 @@ "name": "stderr", "output_type": "stream", "text": [ - "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", + "28-Oct 04:26:43:INFO:root:Testing element MyQuad4\n", + "28-Oct 04:26:43:INFO:root:element dimension: 2 x 4\n", + "28-Oct 04:26:43:INFO:root:Initializing element\n", + "28-Oct 04:26:43:INFO:root:basis at [0.0,0.0]: [0.25 0.25 0.25 0.25]\n", + "28-Oct 04:26:43:INFO:root:field val at [0.0,0.0]: 2.5\n", + "28-Oct 04:26:43: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", - "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" + "28-Oct 04:26:43:INFO:root:field val at [0.0,0.0]: [0.0 2.0]\n", + "28-Oct 04:26:43:INFO:root:Element MyQuad4 passed tests.\n" ] } ], @@ -352,7 +352,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": { "collapsed": false }, @@ -367,7 +367,7 @@ " \"temperature\" => JuliaFEM.FieldSet(\"temperature\",JuliaFEM.Field[JuliaFEM…" ] }, - "execution_count": 13, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } @@ -383,28 +383,28 @@ "push!(geometry_field, [ 0.0, 1.0])\n", "geometry_fieldset = FieldSet(\"geometry\") # create fieldset \"geometry\"\n", "push!(geometry_fieldset, geometry_field) # add field to fieldset\n", - "push!(element, geometry_fieldset) # add fieldset to element\n", + "element[\"geometry\"] = geometry_fieldset # add fieldset to element\n", "\n", "temperature_fieldset = FieldSet(\"temperature\")\n", "push!(temperature_fieldset, Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n", "push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n", - "push!(element, temperature_fieldset)\n", + "element[\"temperature\"] = temperature_fieldset\n", "\n", "displacement_fieldset = FieldSet(\"displacement\")\n", "push!(displacement_fieldset, Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]))\n", "push!(displacement_fieldset, Field(1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]))\n", - "push!(element, displacement_fieldset)\n", + "element[\"displacement\"] = displacement_fieldset\n", "\n", "heat_coefficient_fieldset = FieldSet(\"heat coefficient\")\n", "push!(heat_coefficient_fieldset, Field(0.0, 2))\n", "push!(heat_coefficient_fieldset, Field(1.0, 3))\n", - "push!(element, heat_coefficient_fieldset)\n", + "element[\"heat coefficient\"] = heat_coefficient_fieldset\n", "element.fields" ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 13, "metadata": { "collapsed": false }, @@ -417,7 +417,7 @@ " 0.5" ] }, - "execution_count": 14, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } @@ -431,7 +431,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 14, "metadata": { "collapsed": false }, @@ -443,7 +443,7 @@ " 0.0 1.0" ] }, - "execution_count": 15, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -456,7 +456,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 15, "metadata": { "collapsed": false }, @@ -467,7 +467,7 @@ "2.5" ] }, - "execution_count": 16, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -477,35 +477,35 @@ "basis(\"heat coefficient\", [0.0, 0.0], 0.5)" ] }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "3" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# set time to Inf to get very last field value and -Inf to get first one.\n", + "basis(\"heat coefficient\", [0.0, 0.0], Inf)" + ] + }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, - "outputs": [ - { - "data": { - "text/plain": [ - "3" - ] - }, - "execution_count": 17, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# set time to Inf to get very last field value and -Inf to get first one.\n", - "basis(\"heat coefficient\", [0.0, 0.0], Inf)" - ] - }, - { - "cell_type": "code", - "execution_count": 18, - "metadata": { - "collapsed": false - }, "outputs": [ { "data": { @@ -515,7 +515,7 @@ " 0.0703125 0.0078125" ] }, - "execution_count": 18, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -540,7 +540,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -568,7 +568,7 @@ " 1.0)" ] }, - "execution_count": 19, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -599,7 +599,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "but we have an higher level interface for that also.\n", + "but we have something better than that. Keep on reading!\n", "\n", "### Summary of developing own elements\n", "\n", @@ -634,7 +634,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -654,7 +654,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 20, "metadata": { "collapsed": false }, @@ -662,10 +662,10 @@ { "data": { "text/plain": [ - "size (generic function with 74 methods)" + "size (generic function with 76 methods)" ] }, - "execution_count": 21, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -682,7 +682,9 @@ " 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", + " if !haskey(element, \"temperature\")\n", + " element[\"temperature\"] = FieldSet()\n", + " end\n", " DC2D4(element, integration_points)\n", "end\n", "Base.size(equation::DC2D4) = (1, 4)\n", @@ -695,7 +697,9 @@ "function DC2D2(element::Seg2)\n", " integration_points = [\n", " IntegrationPoint([0.0], 2.0)]\n", - " push!(element, FieldSet(\"temperature\"))\n", + " if !haskey(element, \"temperature\")\n", + " element[\"temperature\"] = FieldSet()\n", + " end\n", " DC2D2(element, integration_points)\n", "end\n", "Base.size(equation::DC2D2) = (1, 2)" @@ -712,7 +716,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 21, "metadata": { "collapsed": false }, @@ -723,7 +727,7 @@ "has_force_vector (generic function with 3 methods)" ] }, - "execution_count": 22, + "execution_count": 21, "metadata": {}, "output_type": "execute_result" } @@ -736,7 +740,7 @@ " element = get_element(equation)\n", " basis = get_basis(element)\n", " ρ = basis(\"density\")\n", - " return ρ(ip, time) * basis(ip,time)'*basis(ip,time)\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", @@ -744,14 +748,14 @@ " basis = get_basis(element)\n", " dbasis = grad(basis)\n", " k = basis(\"temperature thermal conductivity\")\n", - " return k(ip,time) * dbasis(ip,time)'*dbasis(ip,time)\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\")\n", - " return f(ip, time)*basis(ip,time)'\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", @@ -762,7 +766,7 @@ " element = get_element(equation)\n", " basis = get_basis(element)\n", " g = basis(\"temperature flux\")\n", - " return g(ip,time)*basis(ip,time)'\n", + " return g(ip, time)*basis(ip, time)'\n", "end\n", "JuliaFEM.has_force_vector(equation::DC2D2) = true" ] @@ -775,37 +779,68 @@ "\n", "First of all, equations are written in integration or gauss quadrature points for a certain reason. The reason is that for some incremental constitutive models internal variables are needed, which are saved to FieldSets inside integration points. We can always access these variables through `ip.fields` in the same way we can access all other field variables defined inside element via `element.fields`.\n", "\n", - "Secondly, starting point for our equation is that everything can be time-dependent. FieldSet is, like name suggests, a set of fields with different times and/or increments. By transferring time parameter inside equation makes it possible to access not only current field values but earlier values too. In this example field variables thermal conductivity $k$, load $f$ and flux $g$ are constant in time but this way they can depend easily from spatial or temporal dimension. Actually let's change the parameter $g$ such a way that it is linear ramp from $g(t) = 6t$.\n", + "Secondly, starting point for our equation is that everything can be time-dependent. FieldSet is, like name suggests, a set of fields with different times and/or increments. By transferring time parameter inside equation makes it possible to access not only current field values but earlier values too. In this example field variables thermal conductivity $k$ and load $f$ are constant in time but we are planning to set flux $g$ such a way that it is linear ramp from $g(t) = 6t$.\n", "\n", "If we want to try out this formulation, we must create element and assign this equation for it:" ] }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# if no time is given when defining fields it defaults to 0.0.\n", + "# create volume element\n", + "element = Quad4([1, 2, 3, 4])\n", + "fieldset1 = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n", + "fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(6.0)])\n", + "fieldset3 = FieldSet(\"temperature load\", [Field([12.0, 12.0, 12.0, 12.0])])\n", + "fieldset4 = FieldSet(\"density\", [Field(36.0)])\n", + "element[\"geometry\"] = fieldset1\n", + "element[\"temperature thermal conductivity\"] = fieldset2\n", + "element[\"temperature load\"] = fieldset3\n", + "element[\"density\"] = fieldset4\n", + "equation = DC2D4(element)\n", + "\n", + "# create boundary element with \n", + "boundary_element = Seg2([1, 2])\n", + "boundary_element[\"geometry\"] = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0]])])\n", + "# linear ramp from 1 to 6 in time 0 to 1\n", + "boundary_element[\"temperature flux\"] = FieldSet(\"temperature flux\", [Field(0.0, 0.0), Field(1.0, 6.0)])\n", + "boundary_equation = DC2D2(boundary_element);" + ] + }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "4x4 Array{Float64,2}:\n", + " 4.0 2.0 1.0 2.0\n", + " 2.0 4.0 2.0 1.0\n", + " 1.0 2.0 4.0 2.0\n", + " 2.0 1.0 2.0 4.0" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "# create volume element\n", - "element = Quad4([1, 2, 3, 4])\n", - "fieldset1 = FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n", - "fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(0.0, 6.0)])\n", - "fieldset3 = FieldSet(\"temperature load\", [Field(0.0, [12.0, 12.0, 12.0, 12.0])])\n", - "fieldset4 = FieldSet(\"density\", [Field(0.0, 1.0)])\n", - "push!(element, fieldset1)\n", - "push!(element, fieldset2)\n", - "push!(element, fieldset3)\n", - "push!(element, fieldset4)\n", - "equation = DC2D4(element)\n", + "using JuliaFEM: initialize_local_assembly, calculate_local_assembly!\n", "\n", - "# create boundary element with \n", - "boundary_element = Seg2([1, 2])\n", - "push!(boundary_element, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]])]))\n", - "# linear ramp from 1 to 6 in time 0 to 1\n", - "push!(boundary_element, FieldSet(\"temperature flux\", [Field(0.0, 0.0), Field(1.0, 6.0)]))\n", - "boundary_equation = DC2D2(boundary_element);" + "local_assembly = initialize_local_assembly()\n", + "calculate_local_assembly!(local_assembly, equation, \"temperature\", 0.0)\n", + "local_assembly.mass_matrix" ] }, { @@ -831,10 +866,6 @@ } ], "source": [ - "using JuliaFEM: initialize_local_assembly, calculate_local_assembly!\n", - "\n", - "local_assembly = initialize_local_assembly()\n", - "calculate_local_assembly!(local_assembly, equation, \"temperature\")\n", "local_assembly.stiffness_matrix" ] }, @@ -903,7 +934,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Notice that we created a separate boundary element which integrates heat flux. Defining boundary elements goes identically with other elements. With flux $g=6$ we should have constant $T=1$ on free boundary. This should equal to $u(x,y) = x$, i.e. 1d bar stretched from other side while other one is fixed." + "Notice that we created a separate boundary element which integrates heat flux. Defining boundary elements goes identically with other elements. With flux $g=6$ we should have constant $T=1$ on free boundary. This should equal to $u(x,y) = x$, i.e. 1d bar stretched from other end while the other one is fixed." ] }, { @@ -935,6 +966,42 @@ "u[fdofs] = A[fdofs, fdofs] \\ b[fdofs]" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Some quick notes. Passing time all around is not strictly needed and defaults to $t=+\\inf$, when the most recent field is picked from `FieldSet`. Also when defining `Field` without time, it defaults to $t=0$. It's however good practice to interpolate fields inside functions using time parameter which makes the equations more general. If we want to solve the last equation in time $t=0.5$ we should get $T=0.5$:" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2-element Array{Float64,1}:\n", + " 0.5\n", + " 0.5" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "local_assembly = initialize_local_assembly()\n", + "calculate_local_assembly!(local_assembly, boundary_equation, \"temperature\", 0.5) # <-- time added\n", + "b = zeros(4)\n", + "b[fdofs] = local_assembly.force_vector\n", + "u = zeros(4)\n", + "u[fdofs] = A[fdofs, fdofs] \\ b[fdofs]" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -959,7 +1026,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 29, "metadata": { "collapsed": false }, @@ -967,10 +1034,10 @@ { "data": { "text/plain": [ - "size (generic function with 75 methods)" + "size (generic function with 77 methods)" ] }, - "execution_count": 28, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } @@ -987,9 +1054,11 @@ " 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", + " if !haskey(element, \"temperature\")\n", + " element[\"temperature\"] = FieldSet()\n", + " # Initial configuration needs to be defined if using autodiff\n", + " push!(element[\"temperature\"], Field([0.0, 0.0, 0.0, 0.0]))\n", + " end\n", " DC2D4NL(element, integration_points)\n", "end\n", "Base.size(equation::DC2D4NL) = (1, 4)" @@ -1004,7 +1073,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 30, "metadata": { "collapsed": false }, @@ -1015,7 +1084,7 @@ "has_potential_energy (generic function with 2 methods)" ] }, - "execution_count": 29, + "execution_count": 30, "metadata": {}, "output_type": "execute_result" } @@ -1045,7 +1114,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 31, "metadata": { "collapsed": false }, @@ -1054,12 +1123,12 @@ "name": "stdout", "output_type": "stream", "text": [ - "increment 1, |du| = 1.41421, |r| = 4.243\n", - "increment 2, |du| = 0.42426, |r| = 4.243\n", - "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.611092931 seconds\n", + "increment 1, |du| = 1.41421, |r| = 2.121\n", + "increment 2, |du| = 0.42426, |r| = 2.121\n", + "increment 3, |du| = 0.04657, |r| = 2.121\n", + "increment 4, |du| = 0.00057, |r| = 2.121\n", + "increment 5, |du| = 0.00000, |r| = 2.121\n", + "elapsed time: 0.617104273 seconds\n", "error: 1.5543122344752192e-15\n", "temperature at free end: [0.6666666666666682,0.6666666666666682], should be 0.6666666666666666\n" ] @@ -1078,10 +1147,10 @@ " push!(fieldset3, Field(0.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", + " element[\"geometry\"] = fieldset1\n", + " element[\"temperature thermal conductivity\"] = fieldset2\n", + " element[\"temperature load\"] = fieldset3\n", + " element[\"temperature nodal load\"] = fieldset4\n", " equation = DC2D4NL(element)\n", " # create model -- end\n", "\n", @@ -1160,7 +1229,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 32, "metadata": { "collapsed": false }, @@ -1168,10 +1237,10 @@ { "data": { "text/plain": [ - "has_residual_vector (generic function with 2 methods)" + "has_residual_vector (generic function with 3 methods)" ] }, - "execution_count": 31, + "execution_count": 32, "metadata": {}, "output_type": "execute_result" } @@ -1195,7 +1264,12 @@ " 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(\"displacement\"))\n", + " if !haskey(element, \"displacement\")\n", + " element[\"displacement\"] = FieldSet()\n", + " # initial field must be defined if using autodiff\n", + " u0 = Field(Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])\n", + " push!(element[\"displacement\"], u0)\n", + " end\n", " CPS4(element, integration_points)\n", "end\n", "\n", @@ -1234,9 +1308,16 @@ "JuliaFEM.has_residual_vector(equation::CPS4) = true" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There's a tiny Newton solver which can be used to test single element models to verify formulations. It makes the life a bit easier." + ] + }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 33, "metadata": { "collapsed": false }, @@ -1245,7 +1326,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "increment " + "displacement at tip: -8.773031198197748\n" ] }, { @@ -1257,64 +1338,29 @@ " Occurred: -8.773031198197748" ] }, - "execution_count": 32, + "execution_count": 33, "metadata": {}, "output_type": "execute_result" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "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", + "using JuliaFEM: get_integration_points, solve!\n", "\n", "function run_simulation_2()\n", " # create model -- start\n", + " # fields can be created less verbose way. just keep in mind that if time is not defined\n", + " # it defaults to 0.0. fields can be set \"directly\" to the element but it creates a fieldset\n", + " # anyway\n", " element = Quad4([1, 2, 3, 4])\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", + " element[\"geometry\"] = Field(Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]])\n", + " element[\"youngs modulus\"] = Field(500.0)\n", + " element[\"poissons ratio\"] = Field(0.3)\n", + " element[\"displacement volume load\"] = Field(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", - " 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", - " 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[\"displacement\"], new_field)\n", - " @printf(\"increment %2d, |du| = %8.5f\\n\", i, norm(du))\n", - " end\n", - " toc()\n", + " free_dofs = [3, 4, 5, 6]\n", + " solve!(equation, \"displacement\", free_dofs) # launch a newton solver for single element\n", " disp = get_basis(element)(\"displacement\", [1.0, 1.0])[2]\n", " println(\"displacement at tip: $disp\")\n", " # verified using Code Aster.\n", @@ -1336,7 +1382,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 34, "metadata": { "collapsed": false }, @@ -1350,7 +1396,7 @@ " Occurred: [1.0000000000000002,1.0000000000000002]" ] }, - "execution_count": 33, + "execution_count": 34, "metadata": {}, "output_type": "execute_result" } @@ -1360,8 +1406,10 @@ "\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", + " problem=nothing) # <= notice problem\n", + "\n", + " # zero all workspace. always start with this command, otherwise\n", + " # local matrices are added to the old ones.\n", " initialize_local_assembly!(assembly, equation)\n", "\n", " element = get_element(equation)\n", @@ -1376,10 +1424,10 @@ " f = basis(\"temperature load\", ip, time)\n", " # evaluate basis functions and gradient in integration point\n", " N = basis(ip, time)\n", - " ∇N = dbasis(ip, time)\n", + " dN = dbasis(ip, time)\n", " # do assembly\n", " assembly.mass_matrix += w * ρ*N'*N\n", - " assembly.stiffness_matrix += w * k*∇N'*∇N\n", + " assembly.stiffness_matrix += w * k*dN'*dN\n", " assembly.force_vector += w * N'*f\n", " end\n", "end\n", @@ -1387,14 +1435,11 @@ "function test_local_assembly()\n", " # create volume element\n", " element = Quad4([1, 2, 3, 4])\n", - " fieldset1 = FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n", - " fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(0.0, 6.0)])\n", - " fieldset3 = FieldSet(\"temperature load\", [Field(0.0, [12.0, 12.0, 12.0, 12.0])])\n", - " fieldset4 = FieldSet(\"density\", [Field(0.0, 1.0)])\n", - " push!(element, fieldset1)\n", - " push!(element, fieldset2)\n", - " push!(element, fieldset3)\n", - " push!(element, fieldset4)\n", + " # fields can be created even more directly:\n", + " element[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n", + " element[\"temperature thermal conductivity\"] = 6.0\n", + " element[\"temperature load\"] = [12.0, 12.0, 12.0, 12.0]\n", + " element[\"density\"] = 36.0\n", " equation = DC2D4(element)\n", " la = initialize_local_assembly()\n", " calculate_local_assembly!(la, equation, \"temperature\")\n", @@ -1409,89 +1454,18 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Defining own problem\n", - "\n", - "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." + "Note that function takes also problem as input, which allows to access to rest of the equations and elements. From that in the next section." ] }, { - "cell_type": "code", - "execution_count": 34, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "initialize_global_assembly (generic function with 3 methods)" - ] - }, - "execution_count": 34, - "metadata": {}, - "output_type": "execute_result" - } - ], + "cell_type": "markdown", + "metadata": {}, "source": [ - "using JuliaFEM: Problem, FieldProblem, get_equations, get_connectivity\n", - "using JuliaFEM: get_element, get_unknown_field_dimension\n", + "## Defining own problem\n", "\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", - "function PlaneHeatProblem()\n", - " element_mapping = Dict(\n", - " Quad4 => DC2D4,\n", - " Seg2 => DC2D2)\n", - " PlaneHeatProblem(\"temperature\", 1, [], element_mapping)\n", - "end\n", + "The main purpose of `Problem` is to create a mapping between element types and field equations, so that solver knows how to construct equations. Loosely speaking `Problem` is a set of equations and instructions how they can be created from elements. Every time new element is pushed to problem, element mapping creates corresponding equations ready for integration.\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" + "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." ] }, { @@ -1501,67 +1475,33 @@ "collapsed": false }, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "assembling problem for temperature\n" - ] - }, { "data": { "text/plain": [ - "2-element Array{Float64,1}:\n", - " 101.0\n", - " 101.0" + "PlaneHeatProblem" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "dimension of unknown field: 1, problem dofs: 4\n" - ] } ], "source": [ - "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", + "using JuliaFEM: Problem, get_equations, get_connectivity\n", + "using JuliaFEM: get_element, get_unknown_field_dimension\n", "\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", + "type PlaneHeatProblem <: Problem\n", + " unknown_field_name :: ASCIIString\n", + " unknown_field_dimension :: Int\n", + " equations :: Array{Equation, 1}\n", + " element_mapping :: Dict{DataType, DataType}\n", "end\n", - "\n", - "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]" + "function PlaneHeatProblem(equations=[])\n", + " element_mapping = Dict(\n", + " Quad4 => DC2D4,\n", + " Seg2 => DC2D2)\n", + " PlaneHeatProblem(\"temperature\", 1, equations, element_mapping)\n", + "end" ] }, { @@ -1572,10 +1512,7 @@ "source": [ "## Developing own solver\n", "\n", - "Last part. Defining own solver.\n", - "\n", - "- takes a set of problems (typically some field problem + dirichlet boundary problems)\n", - "- solves them, updates fields" + "Solver is the part where all numerical crunching happens. Typically solver takes a set of problems. Each problem has all necessary information to get assembled. Solver then constructs the final problem to be solved. Typically this problem is some field problem + dirichlet boundaries but it can also be e.g. mixed problem. One essential task of solver is after solution to update corresponding fields to elements. Here we give a simple (yet parallel solution inside one machine) demo solver which solves a typical problem of one field and one dirichlet boundary. It isn't optimized any way but is appropriate as simple test solver." ] }, { @@ -1586,19 +1523,19 @@ }, "outputs": [ { - "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" - ] + "data": { + "text/plain": [ + "call (generic function with 1319 methods)" + ] + }, + "execution_count": 36, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "using JuliaFEM: Solver, get_problems, get_unknown_field_name\n", + "using JuliaFEM: Solver, get_problems, get_unknown_field_name, get_gdofs\n", + "using JuliaFEM: initialize_global_assembly, calculate_global_assembly!\n", "\n", "\"\"\" Simple solver for educational purposes. \"\"\"\n", "type SimpleSolver <: Solver\n", @@ -1624,9 +1561,9 @@ " p1, p2 = get_problems(solver)\n", "\n", " ga1 = initialize_global_assembly(p1)\n", - " calculate_global_assembly!(ga1, p1, get_unknown_field_name(p1))\n", + " calculate_global_assembly!(ga1, p1)\n", " ga2 = initialize_global_assembly(p2)\n", - " calculate_global_assembly!(ga2, p2, get_unknown_field_name(p2))\n", + " calculate_global_assembly!(ga2, p2)\n", "\n", " A1 = ga1.stiffness_matrix\n", " b1 = ga1.force_vector\n", @@ -1646,74 +1583,53 @@ " x1 = x[1:length(b1)]\n", " x2 = x[length(b1)+1:end]\n", "\n", - " # check residual\n", - " R1 = A1*x1 - b1\n", - " R2 = A2*x2 - b2\n", - " println(\"Residual norm: $(norm(R1+R2))\")\n", - "\n", " # update field for elements in problem 1\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", - " 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", + " field_name = get_unknown_field_name(p1)\n", + " gdofs = get_gdofs(p1, equation)\n", + " element_solution = full(x1[gdofs])\n", + " field = Field(time, element_solution)\n", " push!(element[field_name], field)\n", " end\n", "\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", - " 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", + " field_name = get_unknown_field_name(p2)\n", + " gdofs = get_gdofs(p2, equation)\n", + " element_solution = full(x2[gdofs])\n", + " field = Field(time, element_solution)\n", " push!(element[field_name], field)\n", " end\n", - "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)" + "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Summary\n", + "It can be seen that we could already use parallel processing to assemble each problem concurrently. Moreover, LU factorization is already parallelized inside one machine, thanks to BLAS.\n", "\n", "Next we collect all things together and solve Poisson equation in domain $\\Omega = [0,1]\\times[0,1]$ with some Neumann boundary on $\\Gamma_1$ and Dirichlet boundary on $\\Gamma_2$." ] }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { - "name": "stdout", + "name": "stderr", "output_type": "stream", "text": [ - "assembling problem for temperature\n" + "28-Oct 04:27:34:INFO:root:assembling problem for temperature\n", + "28-Oct 04:27:34:INFO:root:dimension of unknown field: 1, problem dofs: 4\n", + "28-Oct 04:27:35:INFO:root:assembling problem for reaction force\n", + "28-Oct 04:27:35:INFO:root:dimension of unknown field: 1, problem dofs: 4\n", + "28-Oct 04:27:35:INFO:root:Temperature at point X = [0.5,0.0] is T = 100.00000000000003\n" ] }, { @@ -1725,20 +1641,9 @@ " Occurred: 100.00000000000003" ] }, - "execution_count": 38, + "execution_count": 37, "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": [ @@ -1749,19 +1654,26 @@ "\"\"\"\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", + " # these might look like normal values but believe me, they\n", + " # are fields with temporal and spatial dimension\n", + " el1[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n", + " el1[\"temperature thermal conductivity\"] = 6.0\n", + " el1[\"temperature load\"] = zeros(4)\n", + " el1[\"density\"] = 36.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", - "\n", + " el2[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0]]\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", + " # yet another simplification, if field is given as a tuple,\n", + " # multiple fields are created. there is 1 second time step between\n", + " # each field. So the following is basically same as\n", + " # fieldset = FieldSet(\"temperature flux\")\n", + " # field1 = Field(0.0, 0.0)\n", + " # field2 = Field(1.0, 600.0)\n", + " # push!(fieldset, field1)\n", + " # push!(fieldset, field2)\n", + " # element[\"temperature flux\"] = fieldset\n", + " el2[\"temperature flux\"] = (0.0, 600.0)\n", "\n", " problem1 = PlaneHeatProblem()\n", " push!(problem1, el1)\n", @@ -1769,13 +1681,14 @@ " return problem1\n", "end\n", "\n", + "using JuliaFEM: DirichletProblem\n", + "\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[\"geometry\"] = Vector[[1.0, 1.0], [0.0, 1.0]]\n", " problem2 = DirichletProblem(1)\n", " push!(problem2, el3)\n", " return problem2\n", @@ -1798,7 +1711,7 @@ " 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", + " Logging.info(\"Temperature at point X = $X is T = $T\")\n", " @fact T --> roughly(100.0)\n", "end\n", "\n", @@ -1809,7 +1722,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "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" + "In this notebook the basic instructions how to develop JuliaFEM has been 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/src/JuliaFEM.jl b/src/JuliaFEM.jl index f7f7b67..be86887 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -6,10 +6,10 @@ This is JuliaFEM -- Finite Element Package """ module JuliaFEM -using Lexicon using Logging @Logging.configure(level=DEBUG) +using Lexicon using ForwardDiff autodiffcache = ForwardDiffCache() @@ -28,7 +28,6 @@ function Base.linspace(X1, X2, n) [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)] end - include("types.jl") # type definitions include("interpolate.jl") # interpolation routines @@ -36,21 +35,26 @@ include("interpolate.jl") # interpolation routines include("elements.jl") include("lagrange.jl") # Lagrange elements #include("hierarchical.jl") # P-elements -include("integrate.jl") # integration points +### EQUATIONS ### +include("integrate.jl") # default integration points for elements include("equations.jl") include("problems.jl") -include("solvers.jl") - -# pre- and postprocess -include("xdmf.jl") -include("abaqus_reader.jl") -#include("interfaces.jl") +### FORMULATIION ### include("dirichlet.jl") include("heat.jl") -#include("elasticity_solver.jl") +include("elasticity.jl") + +### ASSEMBLY + SOLVE ### +include("assembly.jl") +include("solvers.jl") + +# PRE AND POSTPROCESS +include("xdmf.jl") +include("abaqus_reader.jl") end # module FEM = JuliaFEM + diff --git a/src/assembly.jl b/src/assembly.jl new file mode 100644 index 0000000..019c603 --- /dev/null +++ b/src/assembly.jl @@ -0,0 +1,63 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# Functions to handle global assembly of problem + +""" Global assembly. """ +type GlobalAssembly <: Assembly + ndofs :: Int + mass_matrix :: SparseMatrixCSC + stiffness_matrix :: SparseMatrixCSC + force_vector :: SparseMatrixCSC +end + +""" Initialize global assembly of size ndofs. """ +function initialize_global_assembly(ndofs::Int=1) + mass_matrix = spzeros(ndofs, ndofs) + stiffness_matrix = spzeros(ndofs, ndofs) + force_vector = spzeros(ndofs, 1) + return GlobalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector) +end + +""" Initialize global assembly, get dimension from problem. """ +function initialize_global_assembly(problem::Problem) + dim, ndofs = size(problem) + return initialize_global_assembly(ndofs) +end + +""" Initialize or empty workspace for global assembly. """ +function initialize_global_assembly!(assembly::GlobalAssembly, problem::Problem) + ndofs = prod(size(problem)) + if ndofs != assembly.ndofs + # if problem size changes, automatically initialize new work space + assembly.ndofs = ndofs + assembly.mass_matrix = spzeros(ndofs, ndofs) + assembly.stiffness_matrix = spzeros(ndofs, ndofs) + assembly.force_vector = spzeros(ndofs, 1) + return + end + # otherwise, empty workspace ready for next iteration + fill!(assembly.mass_matrix, 0.0) + fill!(assembly.stiffness_matrix, 0.0) + fill!(assembly.force_vector, 0.0) + return +end + +""" Calculate global assembly for a problem. """ +function calculate_global_assembly!(assembly::GlobalAssembly, problem::Problem, time::Number=Inf) + + unknown_field_name = get_unknown_field_name(problem) + initialize_global_assembly!(assembly, problem) # zero all + dim, ndofs = size(problem) + Logging.info("assembling problem for $unknown_field_name") + Logging.info("dimension of unknown field: $dim, problem dofs: $ndofs") + local_assembly = initialize_local_assembly() + for (i, equation) in enumerate(get_equations(problem)) + calculate_local_assembly!(local_assembly, equation, unknown_field_name, time) + conn = get_connectivity(get_element(equation)) + gdofs = get_gdofs(problem, equation) + assembly.mass_matrix[gdofs, gdofs] += local_assembly.mass_matrix + assembly.stiffness_matrix[gdofs, gdofs] += local_assembly.stiffness_matrix + assembly.force_vector[gdofs] += local_assembly.force_vector + end +end diff --git a/src/dirichlet.jl b/src/dirichlet.jl index b594331..c63360d 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -34,7 +34,9 @@ function DBC2D2(element::Seg2) integration_points = [ IntegrationPoint([-sqrt(1/3)], 1.0), IntegrationPoint([+sqrt(1/3)], 1.0)] - push!(element, FieldSet("reaction force")) + if !haskey(element, "reaction force") + element["reaction force"] = FieldSet() + end DBC2D2(element, integration_points) end Base.size(equation::DBC2D2) = (1, 2) diff --git a/src/elasticity.jl b/src/elasticity.jl new file mode 100644 index 0000000..f517b7c --- /dev/null +++ b/src/elasticity.jl @@ -0,0 +1,130 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# Elasticity problems + +abstract ElasticityProblem <: Problem +abstract ElasticityEquation <: Equation + +### Formulation ### + +""" Calculate internal energy for elasticity equation. + +Override this to define your own material model. By default we use +Saint Venant-Kirchhoff material model, which is simply + + S(E) = λtr(E) + 2μE +""" +function get_internal_energy(equation::Equation, ip::IntegrationPoint, + time::Number, F::Matrix) + element = get_element(equation) + basis = get_basis(element) + dbasis = grad(basis) + + # material parameters + young = basis("youngs modulus", ip, time) + poisson = basis("poissons ratio", ip, time) + mu = young/(2*(1+poisson)) + lambda = young*poisson/((1+poisson)*(1-2*poisson)) + if isa(equation, PlaneStressElasticityEquation) + lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d + end + + # material model + E = 1/2*(F'*F - I) # strain + S = lambda*trace(E)*I + 2*mu*E + P = F*S + + return P*dbasis(ip, time) +end + +""" Elasticity equations. + +Formulation +----------- + +Field equation is: +∂u/∂t = ∇⋅f - b + +Weak form is: find u∈U such that ∀v in V + + δW := ∫ρ₀∂²u/∂t²⋅δu dV₀ + ∫S:δE dV₀ - ∫b₀⋅δu dV₀ - ∫t₀⋅δu dA₀ = 0 + +where + + ρ₀ = density + b₀ = displacement load + t₀ = displacement traction + +References +---------- + +https://en.wikipedia.org/wiki/Linear_elasticity +https://en.wikipedia.org/wiki/Finite_strain_theory +https://en.wikipedia.org/wiki/Stress_measures +https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid +https://en.wikipedia.org/wiki/Strain_energy_density_function +https://en.wikipedia.org/wiki/Plane_stress +https://en.wikipedia.org/wiki/Hooke's_law + +""" +function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint, + time::Number; variation=nothing) + + element = get_element(equation) + basis = get_basis(element) + dbasis = grad(basis) + + u = basis("displacement", ip, time, variation) + gradu = dbasis("displacement", ip, time, variation) + F = I + gradu # deformation gradient + + # residual vector - internal energy + r = get_internal_energy(equation, ip, time, F) + + # external forces - volume load + if haskey(element, "displacement load") + b = basis("displacement load", ip, time) + r -= b*basis(ip, time) + end + + return vec(r) +end +has_residual_vector(equation::ElasticityEquation) = true + +### Problem 1 - plane elasticity ### + +abstract PlaneElasticityProblem <: ElasticityProblem + +type PlaneStressElasticityProblem <: PlaneElasticityProblem + unknown_field_name :: ASCIIString + unknown_field_dimension :: Int + equations :: Array{ElasticityEquation, 1} + element_mapping :: Dict{DataType, DataType} +end + +function PlaneStressElasticityProblem(equations=[]) + element_mapping = Dict( + Quad4 => CPS4) + return PlaneStressElasticityProblem("displacement", 2, equations, element_mapping) +end + +### Equations ### + +abstract PlaneElasticityEquation <: ElasticityEquation +abstract PlaneStressElasticityEquation <: PlaneElasticityEquation + +type CPS4 <: PlaneStressElasticityEquation + element :: Quad4 + integration_points :: Array{IntegrationPoint, 1} +end +function CPS4(element::Quad4) + integration_points = get_default_integration_points(element) + if !haskey(element, "displacement") + element["displacement"] = FieldSet() + push!(element["displacement"], Field(Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])) + end + CPS4(element, integration_points) +end +Base.size(equation::CPS4) = (2, 4) + diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl deleted file mode 100644 index 946ee1c..0000000 --- a/src/elasticity_solver.jl +++ /dev/null @@ -1,333 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -module elasticity_solver - -using ForwardDiff - -using Logging -@Logging.configure(level=INFO) - -# Below this line is internal functions related to solver. They can be used -# directly if needed or using general interface combining data model and -# solver. - -""" -This is dummy function. Testing doctests and documentation. - -Parameters ----------- -x : Array{Float64, 1} - -Returns -------- -Array{float64, 1} - x + 1 - -Notes ------ -This is dummy function - -Raises ------- -Exception - if things are not going right - -Examples --------- ->>> a = [1.0, 2.0, 3.0] ->>> dummy(a) -[2.0, 3.0, 4.0] -""" -function dummy(a) - # not doing anything useful. - return a+1 -end - - -""" -Calculate local tangent stiffness matrix and residual force vector -R = T - F for elasticity problem. - -Parameters ----------- -X : Element coordinates -u : Displacement field -R : Residual force vector -K : Tangent stiffness matrix -basis : Basis functions -dbasis : Derivative of basis functions -lambda : Material parameter -mu : Material parameter -ipoints : integration points -iweights : integration weights - -Returns -------- -None - -Notes ------ -If material parameters are given in list, they are interpolated to gauss -points using shape functions. -""" -function calc_local_matrices!(X, u, R, K, basis, dbasis, lambda_, mu_, ipoints, iweights) - dim, nnodes = size(X) - I = eye(dim) - R[:,:] = 0.0 - - #dF = zeros(dim, dim) - - function calc_R!(u, R) - for m = 1:length(iweights) - w = iweights[m] - xi = ipoints[m, :] - # calculate material parameters - lambda = typeof(lambda_) == Float64 ? lambda_ : dot(lambda_, basis(xi)) - mu = typeof(mu_) == Float64 ? mu_ : dot(mu_, basis(xi)) - Jt = X*dbasis(xi) - detJ = det(Jt) - dbasisdX = dbasis(xi)*inv(Jt) - - gradu = u*dbasisdX - F = I + gradu # Deformation gradient - E = 1/2*(gradu' + gradu + gradu'*gradu) # Green-Lagrange strain tensor - S = lambda*trace(E)*I + 2*mu*E # PK2 stress tensor - P = F*S # PK1 stress tensor - - R[:,:] += w*P*dbasisdX'*detJ - end - end - - # herlper for tangent stiffness matrix - function R!(u, R) - R[:] = 0 - calc_R!(reshape(u, dim, nnodes), reshape(R, dim, nnodes)) - #calc_Wext!(reshape(u, 2, 4), reshape(R, 2, 4)) - end - Jacobian = ForwardDiff.forwarddiff_jacobian(R!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes) - - K[:, :] = Jacobian(reshape(u, dim*nnodes)) - R!(reshape(u, dim*nnodes), reshape(R, dim*nnodes)) - -end - - -""" -Assemble global stiffness matrix to I,J,V ready for sparse format - -Parameters ----------- -ke : local matrix -eldofs_ : Array - degrees of freedom -I,J,V : Arrays for sparse matrix - -Notes ------ -eldofs can also be node ids for convenience. In that case dimension -is calculated and eldofs are "extended" to problem dimension. -""" -function assemble!(ke, eldofs_, I, J, V) - n, m = size(ke) - dim = round(Int, n/length(eldofs_)) - @debug("problem dim = ", dim) - if dim == 1 - eldofs = eldofs_ - else - eldofs = Int64[] - for i in eldofs_ - for d in 1:dim - push!(eldofs, dim*(i-1)+d) - end - end - @debug("old eldofs", eldofs_) - @debug("new eldofs", eldofs) - end - for i in 1:n - for j in 1:m - push!(I, eldofs[i]) - push!(J, eldofs[j]) - push!(V, ke[i,j]) - end - end -end - - -""" -Assemble global RHS to I,V ready for sparse format - -Parameters ----------- -fe : local vector -eldofs_ : Array - degrees of freedom -I,V : Arrays for sparse matrix - -Notes ------ -eldofs can also be node ids for convenience. In that case dimension -is calculated and eldofs are "extended" to problem dimension. -""" -function assemble!(fe, eldofs_, I, V) - n = length(fe) - dim = round(Int, n/length(eldofs_)) - if dim == 1 - eldofs = eldofs_ - else - eldofs = Int64[] - for i in eldofs_ - for d in 1:dim - push!(eldofs, dim*(i-1)+d) - end - end - end - - for i in 1:n - push!(I, eldofs[i]) - push!(V, fe[i]) - end -end - - -""" -Eliminate Dirichlet boundary conditions from matrix - -Parameters ----------- -dirichletbc : array [dim x nnodes] -I, J, V : sparse matrix arrays - -Returns -------- -I, J, V : boundary conditions removed - -Notes ------ -pros: -- matrix assembly remains positive definite -cons: -- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.) -- if u != 0 in dirichlet boundary requires extra care - -Raises ------- -Exception, if displacement boundary conditions given, i.e. -DX=2 for some node, for example. - -""" -function eliminate_boundary_conditions(dirichletbc, I, J, V) - if any(dirichletbc .> 0) - throw("displacement boundary condition not supported") - end - # dofs to remove - free_dofs = find(isnan(dirichletbc)) - remove_dofs = find(!isnan(dirichletbc)) - @debug("Removing dofs: ", remove_dofs) - # this can be done more clever by removing corresponging indexes from I, J, and V - A = sparse(I, J, V) - A = A[free_dofs, free_dofs] - return findnz(A) -end - -""" -Eliminate Dirichlet boundary conditions from vector - -Parameters ----------- -dirichletbc : array [dim x nnodes] -I, V : sparse vector arrays - -Returns -------- -I, V : boundary conditions removed - -Notes ------ -pros: -- matrix assembly remains positive definite -cons: -- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.) -- if u != 0 in dirichlet boundary requires extra care - -Raises ------- -Exception, if displacement boundary conditions given, i.e. -DX=2 for some node, for example. -""" -function eliminate_boundary_conditions(dirichletbc, I, V) - if any(dirichletbc .> 0) - throw("displacement boundary condition not supported") - end - # dofs to remove - free_dofs = find(isnan(dirichletbc)) - remove_dofs = find(!isnan(dirichletbc)) - @debug("Removing dofs: ", remove_dofs) - # this can be done more clever by removing corresponging indexes from I, J, and V - A = sparsevec(I, V) - A = A[free_dofs] - @debug("new vector: ", A) - i, j, v = findnz(A) - return i, v -end - - - -""" -Solve one increment of elasticity problem -""" -function solve_elasticity_increment!(X, u, du, elmap, nodalloads, - dirichletbc, lambda, mu, N, dNdchi, ipoints, - iweights) - if length(size(elmap)) == 1 - # quick hack for just one element - elmap = elmap'' - end - nelnodes, nelements = size(elmap) - dim, nnodes = size(u) - dofs = dim*nelnodes - - Imat = Int64[] - Jmat = Int64[] - Vmat = Float64[] - Ivec = Int64[] - Vvec = Float64[] - - # FIXME: different number of nodes/element - R = zeros(dim, nelnodes) - Kt = zeros(dofs, dofs) - - # this can be parallelized - for i in 1:nelements - eldofs = elmap[:,i] - calc_local_matrices!(X[:, eldofs], u[:, eldofs], R, Kt, N, dNdchi, - lambda[eldofs], mu[eldofs], ipoints, iweights) - assemble!(Kt, eldofs, Imat, Jmat, Vmat) - assemble!(R, eldofs, Ivec, Vvec) - end - - # add additional neumann boundary conditions to force vector - for (i, nodal_load) in enumerate(nodalloads) - if nodal_load == 0 - continue - end - push!(Ivec, i) - push!(Vvec, -nodal_load) - end - - # Create sparse matrix and vector - A = sparse(Imat, Jmat, Vmat) - b = sparsevec(Ivec, Vvec) - - # Remove dirichlet boundary conditions - free_dofs = find(isnan(dirichletbc)) - #Imat, Jmat, Vmat = eliminate_boundary_conditions(dirichletbc, Imat, Jmat, Vmat) - b = b[free_dofs] - A = A[free_dofs, free_dofs] - - # solution - du[free_dofs] = lufact(A) \ -full(b) -end - - -end diff --git a/src/elements.jl b/src/elements.jl index af2d8dc..681dd6a 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -1,83 +1,11 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -#= -Related notebooks ------------------ - -2015-08-29-developing-juliafem.ipynb -=# - using FactCheck using ForwardDiff abstract Element -""" Get FieldSet from element. """ -function Base.getindex(element::Element, field_name) - element.fields[field_name] -end - -""" Add new FieldSet to element. """ -function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name) - fieldset.name = fieldset_name - element.fields[fieldset.name] = fieldset -end -function Base.push!(element::Element, fieldset::FieldSet) - element[fieldset.name] = fieldset -end - - -#= ELEMENT DEFINITIONS - -Example -------- - -This is example how to create new element. This is commented because I use code -generation for simple elements like Lagrage elements. Feel free to use -code generation but elements can be of course created manually too! - -abstract CG <: Element # create new element family "Continous Galerkin" - -type Quad4 <: CG - connectivity :: Array{Int, 1} - fields :: Dict{Any, Any} -end - -""" Default contructor. """ -Quad4(connectivity) = Quad4(connectivity, Dict{Any, Any}()) - -""" Return number of basis functions of this element. """ -get_number_of_basis_functions(el::Type{Quad4}) = 4 - -""" Return element dimension (length of xi vector). """ -get_element_dimension(el::Type{Quad4}) = 2 - -""" Return basis functions for this element (xi dim = 2, functions = 4). """ -function get_basis(el::Quad4, 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] -end - -""" Return partial derivatives of basis functions. """ -function get_dbasisdxi(el::Quad4, 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] -end - -End of example. - -=# - -# 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 ### - """ Test routine for element. If this passes, element interface is properly defined. @@ -114,9 +42,9 @@ function test_element(element_type) end # try to interpolate some scalar field - push!(element, FieldSet("field1", [Field(0.0, collect(1:n))])) + element["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]])])) + element["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) @@ -134,26 +62,94 @@ function test_element(element_type) Logging.info("Element $element_type passed tests.") end +""" Get FieldSet from element. """ +function Base.getindex(element::Element, field_name) + element.fields[field_name] +end + +"""Add new FieldSet to element. + +Examples +-------- +>>> field = Field(0.0, [1, 2, 3, 4]) +>>> fieldset = FieldSet("geometry", Field[field]) +>>> element["geometry"] = fieldset +JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])]))) +""" +function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name) + fieldset.name = fieldset_name + element.fields[fieldset.name] = fieldset +end + +"""Add new FieldSet to element. + +Examples +-------- +>>> field = Field(0.0, [1, 2, 3, 4]) +>>> element["geometry"] = field +JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])]))) +""" +function Base.setindex!(element::Element, field::Field, fieldset_name) + element[fieldset_name] = FieldSet(field) +end + +"""Add new FieldSet to element. + +Examples +-------- +>>> element["geometry"] = [1, 2, 3, 4] +JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])]))) +""" +function Base.setindex!(element::Element, field_data::Union{Number, Array}, fieldset_name) + element[fieldset_name] = Field(field_data) +end + +"""Add new FieldSet to element. + +Notes +----- +This last version takes tuple and each cell in tuple is converted to new field. +Time in field is 0.0, 1.0, ..., n + +Examples +-------- +>>> element["load"] = (1, 2) +JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("load"=>JuliaFEM.FieldSet("load",JuliaFEM.Field[JuliaFEM.Field{Int64}(0.0,0,1),JuliaFEM.Field{Int64}(1.0,0,2)]))) +""" +function Base.setindex!(element::Element, field_data::Tuple, fieldset_name) + fields = Field[Field(Float64(i-1), field) for (i,field) in enumerate(field_data)] + element.fields[fieldset_name] = FieldSet(fieldset_name, fields) +end function get_connectivity(el::Element) el.connectivity end -type MixedFunctionSpace +abstract AbstractFunctionSpace + +type FunctionSpace <: AbstractFunctionSpace + element :: Element +end + +type GradientFunctionSpace <: AbstractFunctionSpace + element :: Element +end + +type MixedFunctionSpace <: AbstractFunctionSpace element1 :: Element element2 :: Element end -type FunctionSpace - element :: Element +function get_basis(element::Element) + return FunctionSpace(element) end -type GradientFunctionSpace - element :: Element +function get_dbasis(element::Element) + return GradientFunctionSpace(element) end function grad(u::FunctionSpace) - GradientFunctionSpace(u.element) + return GradientFunctionSpace(u.element) end """ Evaluate field on element function space. """ @@ -195,8 +191,8 @@ end 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, ip::IntegrationPoint, t::Number=Inf, variation=nothing) = call(u, field_name, ip.xi, t, variation) +call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number=Inf, 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...) @@ -215,6 +211,7 @@ function get_fieldset(u::FunctionSpace, field_name) return u.element[field_name] end +""" Get a determinant of element in point ξ. """ function LinAlg.det(u::FunctionSpace, xi::Vector, t::Number=Inf) X = u.element["geometry"](t) dN = u.element.basis.dbasisdxi(xi) @@ -229,156 +226,13 @@ function LinAlg.det(u::FunctionSpace) return (args...) -> det(u, args...) end -function get_basis(element::Element) - return FunctionSpace(element) -end - Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...) Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...) Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...) Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) - v(args...) - """ Check does fieldset exist. """ function Base.haskey(element::Element, what) haskey(element.fields, what) end - -# FIXME: These two needs integration -- maybe not in elements.jl ..? -""" -Fit field s.t. || ∫ (Nᵢ(ξ)αᵢ - f(el, ξ)) dS || -> min! - -Parameters ----------- -f::Function - Needs to take (el::Element, xi::Vector) as argument -fixed_coeffs::Int[] - These coefficients are not changed during fitting -> constrained optimizatio -""" -function fit_field!(el::Element, field, f, fixed_coeffs=Int[]) - w = [ - 128/225, - (332+13*sqrt(70))/900, - (332+13*sqrt(70))/900, - (332-13*sqrt(70))/900, - (332-13*sqrt(70))/900] - xi = Vector[ - [0.0], - [ 1/3*sqrt(5 - 2*sqrt(10/7))], - [-1/3*sqrt(5 - 2*sqrt(10/7))], - [ 1/3*sqrt(5 + 2*sqrt(10/7))], - [-1/3*sqrt(5 + 2*sqrt(10/7))]] - n = get_number_of_basis_functions(el) - fld = get_field(el, field) - nfld = length(fld[1]) - #Logging.debug("dim of field $field: $nfld") - - M = zeros(n, n) - b = zeros(n, nfld) - for i=1:length(w) - detJ = get_detJ(el, xi[i]) - N = get_basis(el, xi[i]) - M += w[i]*N*N'*detJ - fi = f(el, xi[i]) - for j=1:nfld - b[:, j] += w[i]*N*fi[j]*detJ - end - end - - coeffs = zeros(n) - for j=1:nfld - for k=1:n - coeffs[k] = fld[k][j] - end - if length(fixed_coeffs) != 0 - # constrained problem, some coefficients are fixed - N = Int[] # rest of coeffs - S = Int[] # fixed coeffs - for i = 1:n - if i in fixed_coeffs - push!(S, i) - else - push!(N, i) - end - end - lhs = M[N,N] - rhs = b[N,j] - M[N,S]*coeffs[S] - coeffs[N] = lhs \ rhs - else - coeffs[:] = M \ b[:,j] - end - for k=1:n - fld[k][j] = coeffs[k] - end - end - set_field(el, field, fld) - return -end - - -""" -Fit field s.t. || ∫ ∂/∂ξ(∑Nᵢ(ξ)αᵢ)f(el, ξ) dS || -> min! -""" -function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[]) - w = [ - 128/225, - (332+13*sqrt(70))/900, - (332+13*sqrt(70))/900, - (332-13*sqrt(70))/900, - (332-13*sqrt(70))/900] - xi = Vector[ - [0.0], - [ 1/3*sqrt(5 - 2*sqrt(10/7))], - [-1/3*sqrt(5 - 2*sqrt(10/7))], - [ 1/3*sqrt(5 + 2*sqrt(10/7))], - [-1/3*sqrt(5 + 2*sqrt(10/7))]] - n = get_number_of_basis_functions(el) - fld = get_field(el, field) - nfld = length(fld[1]) - #Logging.debug("dim of field $field: $nfld") - - M = zeros(n, n) - b = zeros(n, nfld) - for i=1:length(w) - detJ = get_detJ(el, xi[i]) - dNdxi = get_dbasisdxi(el, xi[i]) - dNdX = dNdxi / detJ - M += w[i]*dNdX*dNdX'*detJ - fi = f(el, xi[i]) - for j=1:nfld - b[:, j] += w[i]*dNdX*fi[j]*detJ - end - end - - coeffs = zeros(n) - for j=1:nfld - for k=1:n - coeffs[k] = fld[k][j] - end - if length(fixed_coeffs) != 0 - #Logging.info("constrained problem, some coefficients are fixed") - N = Int[] # rest of coeffs - S = Int[] # fixed coeffs - for i = 1:n - if i in fixed_coeffs - push!(S, i) - else - push!(N, i) - end - end - lhs = M[N,N] - rhs = b[N,j] - M[N,S]*coeffs[S] - coeffs[N] = lhs \ rhs - else - coeffs[:] = M \ b[:,j] - end - for k=1:n - fld[k][j] = coeffs[k] - end - end - set_field(el, field, fld) - return -end - - diff --git a/src/equations.jl b/src/equations.jl index 16ca2da..25555d3 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -1,6 +1,8 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +# Functions to handle element level things -- integration, assembly, ... + abstract Equation abstract Assembly @@ -55,17 +57,61 @@ function initialize_local_assembly!(assembly::LocalAssembly, equation::Equation) end has_mass_matrix(equation::Equation) = false -get_mass_matrix(equation::Equation, ip, time) = nothing +function get_mass_matrix(equation::Equation, ip, time=Inf, problem=nothing) + get_mass_matrix(equation, ip, time) +end +function get_mass_matrix(equation::Equation, ip, time=Inf) + get_mass_matrix(equation, ip) +end +function get_mass_matrix(equation::Equation, ip) + nothing +end + has_stiffness_matrix(equation::Equation) = false -get_stiffness_matrix(equation::Equation, ip, time) = nothing +function get_stiffness_matrix(equation::Equation, ip, time=Inf, problem=nothing) + get_stiffness_matrix(equation, ip, time) +end +function get_stiffness_matrix(equation::Equation, ip, time=Inf) + get_stiffness_matrix(equation, ip) +end +function get_stiffness_matrix(equation::Equation, ip) + nothing +end + has_force_vector(equation::Equation) = false -get_force_vector(equation::Equation, ip, time) = nothing +function get_force_vector(equation::Equation, ip, time=Inf, problem=nothing) + get_force_vector(equation, ip, time) +end +function get_force_vector(equation::Equation, ip, time=Inf) + get_force_vector(equation, ip) +end +function get_force_vector(equation::Equation, ip) + nothing +end + has_residual_vector(equation::Equation) = false -get_residual_vector(equation::Equation, ip, time) = nothing +function get_residual_vector(equation::Equation, ip, time=Inf, problem=nothing) + get_residual_vector(equation, ip, time) +end +function get_residual_vector(equation::Equation, ip, time=Inf) + get_residual_vector(equation, ip) +end +function get_residual_vector(equation::Equation, ip) + nothing +end + has_potential_energy(equation::Equation) = false -get_potential_energy(equation::Equation, ip, time) = nothing +function get_potential_energy(equation::Equation, ip, time=Inf, problem=nothing) + get_potential_energy(equation, ip, time) +end +function get_potential_energy(equation::Equation, ip, time=Inf) + get_potential_energy(equation, ip) +end +function get_potential_energy(equation::Equation, ip) + nothing +end + get_element(equation::Equation) = equation.element -get_number_of_dofs(equation::Equation) = nothing get_integration_points(equation::Equation) = equation.integration_points @@ -85,19 +131,19 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, for ip in get_integration_points(equation) s = ip.weight*detJ(ip) if has_mass_matrix(equation) - assembly.mass_matrix += s*get_mass_matrix(equation, ip, time) + assembly.mass_matrix += s*get_mass_matrix(equation, ip, time, problem) end if has_stiffness_matrix(equation) - assembly.stiffness_matrix += s*get_stiffness_matrix(equation, ip, time) + assembly.stiffness_matrix += s*get_stiffness_matrix(equation, ip, time, problem) end if has_force_vector(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, "$unknown_field_name nodal load") - assembly.force_vector += element["$unknown_field_name nodal load"](time)[:] + assembly.force_vector += s*get_force_vector(equation, ip, time, problem) end end + # external loads -- if any nodal loads is defined add to force vector + if haskey(element, "$unknown_field_name nodal load") + assembly.force_vector += vec(element["$unknown_field_name nodal load"](time)) + end end # 2. variational / energy form - user has defined some potential energy / variational form @@ -116,14 +162,14 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, # external energy -- if any nodal loads is defined, decrease from potential energy if haskey(element, "$unknown_field_name nodal load") P = element["$unknown_field_name nodal load"](time) - assembly.potential_energy -= dot(P[:], df[:]) + assembly.potential_energy -= dot(vec(P), vec(df)) end if isa(assembly.potential_energy, Array) return assembly.potential_energy[1] end return assembly.potential_energy end - hessian, allresults = ForwardDiff.hessian(potential_energy, field[:], + hessian, allresults = ForwardDiff.hessian(potential_energy, vec(field), AllResults, cache=autodiffcache) assembly.stiffness_matrix += hessian assembly.force_vector -= ForwardDiff.gradient(allresults) # <--- minus explained in tutorial @@ -144,15 +190,14 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, end # external loads -- if any nodal loads is defined, remove from residual if haskey(element, "$unknown_field_name nodal load") - assembly.residual_vector -= element["$unknown_field_name nodal load"](time)[:] + assembly.residual_vector -= vec(element["$unknown_field_name nodal load"](time)) end return assembly.residual_vector end - jacobian, allresults = ForwardDiff.jacobian(residual_vector, field[:], + jacobian, allresults = ForwardDiff.jacobian(residual_vector, vec(field), AllResults, cache=autodiffcache) assembly.stiffness_matrix += jacobian assembly.force_vector -= ForwardDiff.value(allresults) # <-- minus explained in tutorial end end - diff --git a/src/field_utils.jl b/src/field_utils.jl new file mode 100644 index 0000000..301dcdf --- /dev/null +++ b/src/field_utils.jl @@ -0,0 +1,140 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# There are here for now until I figure a better place for them. + + +""" +Fit field s.t. || ∫ (Nᵢ(ξ)αᵢ - f(el, ξ)) dS || -> min! + +Parameters +---------- +f::Function + Needs to take (el::Element, xi::Vector) as argument +fixed_coeffs::Int[] + These coefficients are not changed during fitting -> constrained optimizatio +""" +function fit_field!(el::Element, field, f, fixed_coeffs=Int[]) + w = [ + 128/225, + (332+13*sqrt(70))/900, + (332+13*sqrt(70))/900, + (332-13*sqrt(70))/900, + (332-13*sqrt(70))/900] + xi = Vector[ + [0.0], + [ 1/3*sqrt(5 - 2*sqrt(10/7))], + [-1/3*sqrt(5 - 2*sqrt(10/7))], + [ 1/3*sqrt(5 + 2*sqrt(10/7))], + [-1/3*sqrt(5 + 2*sqrt(10/7))]] + n = get_number_of_basis_functions(el) + fld = get_field(el, field) + nfld = length(fld[1]) + #Logging.debug("dim of field $field: $nfld") + + M = zeros(n, n) + b = zeros(n, nfld) + for i=1:length(w) + detJ = get_detJ(el, xi[i]) + N = get_basis(el, xi[i]) + M += w[i]*N*N'*detJ + fi = f(el, xi[i]) + for j=1:nfld + b[:, j] += w[i]*N*fi[j]*detJ + end + end + + coeffs = zeros(n) + for j=1:nfld + for k=1:n + coeffs[k] = fld[k][j] + end + if length(fixed_coeffs) != 0 + # constrained problem, some coefficients are fixed + N = Int[] # rest of coeffs + S = Int[] # fixed coeffs + for i = 1:n + if i in fixed_coeffs + push!(S, i) + else + push!(N, i) + end + end + lhs = M[N,N] + rhs = b[N,j] - M[N,S]*coeffs[S] + coeffs[N] = lhs \ rhs + else + coeffs[:] = M \ b[:,j] + end + for k=1:n + fld[k][j] = coeffs[k] + end + end + set_field(el, field, fld) + return +end + + +""" +Fit field s.t. || ∫ ∂/∂ξ(∑Nᵢ(ξ)αᵢ)f(el, ξ) dS || -> min! +""" +function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[]) + w = [ + 128/225, + (332+13*sqrt(70))/900, + (332+13*sqrt(70))/900, + (332-13*sqrt(70))/900, + (332-13*sqrt(70))/900] + xi = Vector[ + [0.0], + [ 1/3*sqrt(5 - 2*sqrt(10/7))], + [-1/3*sqrt(5 - 2*sqrt(10/7))], + [ 1/3*sqrt(5 + 2*sqrt(10/7))], + [-1/3*sqrt(5 + 2*sqrt(10/7))]] + n = get_number_of_basis_functions(el) + fld = get_field(el, field) + nfld = length(fld[1]) + #Logging.debug("dim of field $field: $nfld") + + M = zeros(n, n) + b = zeros(n, nfld) + for i=1:length(w) + detJ = get_detJ(el, xi[i]) + dNdxi = get_dbasisdxi(el, xi[i]) + dNdX = dNdxi / detJ + M += w[i]*dNdX*dNdX'*detJ + fi = f(el, xi[i]) + for j=1:nfld + b[:, j] += w[i]*dNdX*fi[j]*detJ + end + end + + coeffs = zeros(n) + for j=1:nfld + for k=1:n + coeffs[k] = fld[k][j] + end + if length(fixed_coeffs) != 0 + #Logging.info("constrained problem, some coefficients are fixed") + N = Int[] # rest of coeffs + S = Int[] # fixed coeffs + for i = 1:n + if i in fixed_coeffs + push!(S, i) + else + push!(N, i) + end + end + lhs = M[N,N] + rhs = b[N,j] - M[N,S]*coeffs[S] + coeffs[N] = lhs \ rhs + else + coeffs[:] = M \ b[:,j] + end + for k=1:n + fld[k][j] = coeffs[k] + end + end + set_field(el, field, fld) + return +end diff --git a/src/heat.jl b/src/heat.jl index 8808d96..faf73b1 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -6,7 +6,66 @@ abstract HeatProblem <: Problem abstract HeatEquation <: Equation -### Plane heat problem + equations ### +### Formulation ### + +""" Heat equations. + +Formulation +----------- + +Field equation is: + + ρc∂u/∂t = ∇⋅(k∇u) + f + +Weak form is: find u∈U such that ∀v in V + + ∫k∇u∇v dx = ∫fv dx + ∫gv ds, + +where + + k = temperature thermal conductivity defined on volume + f = temperature load defined on volume + g = temperature flux defined on boundary + +References +---------- +https://en.wikipedia.org/wiki/Heat_equation + +""" +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) + basis = get_basis(element) + dbasis = grad(basis) + detJ = det(basis) + for ip in get_integration_points(equation) + w = ip.weight * detJ(ip) + N = basis(ip, time) + if haskey(element, "density") + ρ = basis("density", ip, time) + assembly.mass_matrix += w * ρ*N'*N + end + if haskey(element, "temperature thermal conductivity") + dN = dbasis(ip, time) + k = basis("temperature thermal conductivity", ip, time) + assembly.stiffness_matrix += w * k*dN'*dN + end + if haskey(element, "temperature load") + f = basis("temperature load", ip, time) + assembly.force_vector += w * N'*f + end + if haskey(element, "temperature flux") + g = basis("temperature flux", ip, time) + assembly.force_vector += w * N'*g + end + end +end + +### Problems ### type PlaneHeatProblem <: HeatProblem unknown_field_name :: ASCIIString @@ -23,7 +82,7 @@ function PlaneHeatProblem() return PlaneHeatProblem("temperature", 1, [], element_mapping) end - +### Equations ### """ Diffusive heat transfer for 4-node bilinear element. """ type DC2D4 <: HeatEquation @@ -32,7 +91,9 @@ type DC2D4 <: HeatEquation end function DC2D4(element::Quad4) integration_points = get_default_integration_points(element) - push!(element, FieldSet("temperature")) + if !haskey(element, "temperature") + element["temperature"] = FieldSet() + end DC2D4(element, integration_points) end Base.size(equation::DC2D4) = (1, 4) @@ -44,34 +105,10 @@ type DC2D2 <: HeatEquation end function DC2D2(element::Seg2) integration_points = get_default_integration_points(element) - push!(element, FieldSet("temperature")) + if !haskey(element, "temperature") + element["temperature"] = FieldSet() + end DC2D2(element, integration_points) end 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) - 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 - diff --git a/src/lagrange.jl b/src/lagrange.jl index 1707c42..421ed0b 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -15,7 +15,6 @@ function calculate_lagrange_basis(P, X) for i=1:nbasis A[i,:] = P(X[:, i]) end -# Logging.debug("Calculating inverse of A") invA = inv(A)' basis(xi) = (invA*P(xi))' dbasisdxi(xi) = (ForwardDiff.jacobian((xi) -> invA*P(xi), xi, cache=autodiffcache))' @@ -34,9 +33,6 @@ 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, dbasisdxi = calculate_lagrange_basis($P, $X) type $eltype <: CG connectivity :: Array{Int, 1} @@ -47,8 +43,6 @@ macro create_lagrange_element(element_name, element_description, X, P) $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 Base.size(el::Type{$eltype}) = Base.size($X) end end diff --git a/src/problems.jl b/src/problems.jl index 69e4abe..17a829e 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -5,22 +5,72 @@ abstract Problem abstract BoundaryProblem <: Problem abstract FieldProblem <: Problem +""" Return all equations beloging to this problem. """ function get_equations(problem::Problem) problem.equations end +""" Return the dimension of the unknown field of this problem. """ function get_unknown_field_dimension(problem::Problem) problem.unknown_field_dimension end +""" Return the name of the unknown field of this problem. """ function get_unknown_field_name(problem::Problem) problem.unknown_field_name end -""" Add new element to problem. """ +""" +Add new equation to problem. + +Parameters +---------- +problem +element + +Notes +----- +Equation is automatically created during process based on problem +element -> equation mapping and element type. +""" function Base.push!(problem::Problem, element::Element) element_type = typeof(element) equation_type = problem.element_mapping[element_type] push!(problem.equations, equation_type(element)) end +""" +Return the size of the problem, i.e. +maximum number of connectivity × unknown field dimension. +""" +function Base.size(problem::Problem) + mc = 0 + for equation in get_equations(problem) + element = get_element(equation) + mc = max(mc, get_connectivity(element)...) + end + dim = get_unknown_field_dimension(problem) + return (dim, dim*mc) +end + +""" Assign new equation mapping to problem for some element. + +Examples +-------- +>>> p = PlaneHeatProblem() +>>> p[Seg2] = DC2D2 + +""" +function Base.setindex!(problem::Problem, equation, element) + problem.element_mapping[element] = equation +end + +""" Return global degrees of freedom of element in matrix level. """ +function get_gdofs(problem::Problem, equation::Equation) + dim = get_unknown_field_dimension(problem) + element = get_element(equation) + conn = get_connectivity(element) + gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...)) + return gdofs +end + diff --git a/src/shape_functions.jl b/src/shape_functions.jl deleted file mode 100644 index e49c6c3..0000000 --- a/src/shape_functions.jl +++ /dev/null @@ -1,3 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - diff --git a/src/solvers.jl b/src/solvers.jl index 58532b2..6161b5f 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -5,17 +5,89 @@ abstract Solver +""" +Solve field equations for single element with some dofs fixed. This can be used +to test nonlinear element formulations. +""" +function solve!(equation::Equation, unknown_field_name::ASCIIString, + free_dofs::Array{Int, 1}, time::Number=Inf; + max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false) + element = get_element(equation) + x0 = element[unknown_field_name](-Inf) + x = zeros(prod(size(equation))) + dx = fill!(similar(x), 0.0) + la = initialize_local_assembly() + for i=1:max_iterations + calculate_local_assembly!(la, equation, unknown_field_name) + A = la.stiffness_matrix[free_dofs, free_dofs] + b = la.force_vector[free_dofs] + if dump_matrices + dump(full(A)) + dump(full(b)') + end + dx[free_dofs] = A \ b + x += dx + new_field = similar(x0, x) + new_field.time = time + new_field.increment = i + push!(element[unknown_field_name], new_field) + if norm(dx) < tolerance + return + end + end + Logging.err("Did not converge in $max_iterations iterations") +end + +""" +Solve field equations for a single problem with some dofs fixed. This can be used +to test nonlinear element formulations. Dirichlet boundary is assumed to be homogeneous +and degrees of freedom are eliminated. So if boundary condition is known in nodal +points and everything is zero this should be quite good. +""" +function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=Inf; + max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false) + ga = initialize_global_assembly(problem) + x = zeros(ga.ndofs) + dx = fill!(similar(x), 0.0) + field_name = get_unknown_field_name(problem) + dim = get_unknown_field_dimension(problem) + for i=1:max_iterations + calculate_global_assembly!(ga, problem) + A = ga.stiffness_matrix[free_dofs, free_dofs] + b = ga.force_vector[free_dofs] + if dump_matrices + dump(full(A)) + dump(full(b)') + end + dx[free_dofs] = lufact(A) \ full(b) + x += dx + for equation in get_equations(problem) + element = get_element(equation) + conn = get_connectivity(element) + gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...)) + old_field = element[field_name](Inf) + new_field = similar(old_field, full(x[gdofs])) + new_field.time = time + new_field.increment = i + push!(element[field_name], new_field) + end + if norm(dx) < tolerance + return + end + end + Logging.err("Did not converge in $max_iterations iterations") +end + """ Add new problem to solver. """ function add_problem!(solver::Solver, problem::Problem) push!(solver.problems, problem) end + function Base.push!(solver::Solver, problem::Problem) push!(solver.problems, problem) end -""" -Get all problems assigned to solver -""" +""" Get all problems assigned to solver. """ function get_problems(s::Solver) return s.problems end @@ -25,6 +97,7 @@ end type SimpleSolver <: Solver problems end + """ Default initializer. """ function SimpleSolver() SimpleSolver(Problem[]) @@ -33,34 +106,32 @@ end """ Call solver to solve a set of problems. -This is simple serial solver for demonstration purposes. It handles the most +This is a simple direct solver for demonstration purposes. It handles the common situation, i.e., some main field problem and it's Dirichlet boundary. + + Au + C'λ = f + Cu = g + """ -function call(solver::SimpleSolver, t) - problems = get_problems(solver) - problem1 = problems[1] - problem2 = problems[2] +function call(solver::SimpleSolver, time::Number=Inf) + p1, p2 = get_problems(solver) - # calculate order of degrees of freedom in global matrix - # and set the ordering to problems - dofmap = calculate_global_dofs(problem1) - assign_global_dofs!(problem1, dofmap) - assign_global_dofs!(problem2, dofmap) + ga1 = initialize_global_assembly(p1) + calculate_global_assembly!(ga1, p1) + ga2 = initialize_global_assembly(p2) + calculate_global_assembly!(ga2, p2) - # assemble problem 1 - A1 = sparse(get_lhs(problem1, t)...) - b1 = sparsevec(get_rhs(problem1, t)..., size(A1, 1)) - - # assemble problem 2 - A2 = sparse(get_lhs(problem2, t)...) - b2 = sparsevec(get_rhs(problem2, t)..., size(A2, 1)) - - # make one monolithic assembly + A1 = ga1.stiffness_matrix + b1 = ga1.force_vector + A2 = ga2.stiffness_matrix + b2 = ga2.force_vector + + # create a saddle point problem A = [A1 A2; A2' zeros(A2)] b = [b1; b2] # solve problem - nz = unique(rowvals(A)) + nz = unique(rowvals(A)) # here we remove any zero rows x = zeros(b) x[nz] = lufact(A[nz,nz]) \ full(b[nz]) @@ -68,26 +139,23 @@ function call(solver::SimpleSolver, t) x1 = x[1:length(b1)] x2 = x[length(b1)+1:end] - # check residual - R1 = A1*x1 - b1 - R2 = A2*x2 - b2 - println("Residual norm: $(norm(R1+R2))") - # update field for elements in problem 1 - for equation in get_equations(problem1) - gdofs = get_global_dofs(equation) + for equation in get_equations(p1) element = get_element(equation) - field_name = get_unknown_field_name(equation) # field we are solving - field = Field(t, full(x1[gdofs])[:]) + field_name = get_unknown_field_name(p1) + gdofs = get_gdofs(p1, equation) + element_solution = full(x1[gdofs]) + field = Field(time, element_solution) push!(element[field_name], field) end - # update field for elements in problem 2 - for equation in get_equations(problem2) - gdofs = get_global_dofs(equation) + # update field for elements in problem 2 (Dirichlet boundary) + for equation in get_equations(p2) element = get_element(equation) - field_name = get_unknown_field_name(equation) - field = Field(t, full(x2[gdofs])) + field_name = get_unknown_field_name(p2) + gdofs = get_gdofs(p2, equation) + element_solution = full(x2[gdofs]) + field = Field(time, element_solution) push!(element[field_name], field) end end diff --git a/src/types.jl b/src/types.jl index 5fcfbea..5b34528 100644 --- a/src/types.jl +++ b/src/types.jl @@ -7,7 +7,7 @@ using ForwardDiff """ Field is a fundamental data type which holds some values in some time t """ type Field{T} - time :: Float64 + time :: Number increment :: Int64 values :: T end @@ -106,12 +106,15 @@ type FieldSet fields :: Array{Field, 1} end """ Initializer for FieldSet. """ -function FieldSet(field_name) +function FieldSet(field_name::ASCIIString) FieldSet(field_name, []) end function FieldSet() FieldSet("unknown field", []) end +function FieldSet(fields::Array{Field, 1}) + FieldSet("unknown field", fields) +end """ Add new field to fieldset. """ function Base.push!(fs::FieldSet, field::Field) push!(fs.fields, field) @@ -130,6 +133,9 @@ end function Base.endof(fieldset::FieldSet) length(fieldset) end +function Base.convert(fieldset::Type{FieldSet}, field::Field) + FieldSet(Field[field]) +end """ Basis function. """ diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl new file mode 100644 index 0000000..267cfe1 --- /dev/null +++ b/test/test_elasticity.jl @@ -0,0 +1,26 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# unit tests for heat equations + +using FactCheck +using JuliaFEM: Quad4, Field, FieldSet, CPS4, get_basis, solve!, PlaneStressElasticityProblem + + +facts("test plane elasticity on single element, volume load") do + element = Quad4([1, 2, 3, 4]) + element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]])) + element["youngs modulus"] = FieldSet(Field(500.0)) + element["poissons ratio"] = FieldSet(Field(0.3)) + element["displacement load"] = FieldSet(Field(0.0, Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]])) + equation = CPS4(element) + free_dofs = [3, 4, 5, 6] + problem = PlaneStressElasticityProblem([equation]) + solve!(problem, free_dofs; max_iterations=10) + #solve!(equation, "displacement", free_dofs; max_iterations=10) + disp = get_basis(element)("displacement", [1.0, 1.0])[2] + Logging.info("displacement at tip: $disp") + # verified using Code Aster. + @fact disp --> roughly(-8.77303119819776E+00) +end + diff --git a/test/test_elements.jl b/test/test_elements.jl index 314e400..2f7cdb6 100644 --- a/test/test_elements.jl +++ b/test/test_elements.jl @@ -34,12 +34,14 @@ end facts("test adding fieldsets and fields to element") do el = MockElement([1, 2, 3, 4]) + fieldset = JuliaFEM.FieldSet("geometry") field1 = JuliaFEM.Field(0.0, [0.0, 0.0, 0.0, 0.0]) push!(fieldset, field1) field2 = JuliaFEM.Field(1.0, [1.0, 1.0, 1.0, 1.0]) push!(fieldset, field2) - push!(el, fieldset) + + el["geometry"] = fieldset fields = el["geometry"] @fact length(fields) --> 2 @fact fields[1] --> field1 @@ -57,13 +59,13 @@ facts("interpolation of fields in some function space") do fieldset6 = FieldSet("vector field 3", [Field(0.0, Vector[[1.0, 5.0, 9.0], [2.0, 6.0, 10.0], [3.0, 7.0, 11.0], [4.0, 8.0, 12.0]])]) fieldset7 = FieldSet("tensor field 1", [Field(0.0, Matrix[[1.0 5.0; 9.0 13.0], [2.0 6.0; 10.0 14.0], [3.0 7.0; 11.0 15.0], [4.0 8.0; 12.0 16.0]])]) - push!(element, fieldset1) - push!(element, fieldset2) - push!(element, fieldset3) - push!(element, fieldset4) - push!(element, fieldset5) - push!(element, fieldset6) - push!(element, fieldset7) + element["geometry"] = fieldset1 + element["constant scalar field"] = fieldset2 + element["scalar field"] = fieldset3 + element["vector field 1"] = fieldset4 + element["vector field 2"] = fieldset5 + element["vector field 3"] = fieldset6 + element["tensor field 1"] = fieldset7 xi = [0.0, 0.0] t = 0.0 diff --git a/test/test_global_assembly.jl b/test/test_global_assembly.jl new file mode 100644 index 0000000..6862590 --- /dev/null +++ b/test/test_global_assembly.jl @@ -0,0 +1,32 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM: Quad4, Seg2, FieldSet, Field, PlaneHeatProblem +using JuliaFEM: initialize_global_assembly, calculate_global_assembly! +using FactCheck + +facts("assemble a simple two element problem and solve") do + el1 = Quad4([1, 2, 3, 4]) + el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] + el1["temperature thermal conductivity"] = 6.0 + el1["temperature load"] = [12.0, 12.0, 12.0, 12.0] + el1["density"] = 10 + + el2 = Seg2([1, 2]) + el2["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]] + + # Boundary load, linear ramp 0 -> 600 at time 0 -> 1 + el2["temperature flux"] = FieldSet(Field[Field(0.0, 0.0), Field(1.0, 600.0)]) + + problem = PlaneHeatProblem() + push!(problem, el1) + push!(problem, el2) + + global_assembly = initialize_global_assembly(problem) + calculate_global_assembly!(global_assembly, problem) + free_dofs = [1, 2] + A = lufact(global_assembly.stiffness_matrix[free_dofs, free_dofs]) + b = full(global_assembly.force_vector)[free_dofs] + u = A \ b + @fact u --> roughly([101.0, 101.0]) +end diff --git a/test/test_heat.jl b/test/test_heat.jl new file mode 100644 index 0000000..5f48c04 --- /dev/null +++ b/test/test_heat.jl @@ -0,0 +1,43 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# unit tests for heat equations + +using FactCheck +using JuliaFEM: Seg2, Quad4, Field, FieldSet, DC2D4, initialize_local_assembly, calculate_local_assembly!, DC2D2 + +facts("tests on [0x1]x[0x1] domain") do + + # volume element + element = Quad4([1, 2, 3, 4]) + element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])) + element["temperature thermal conductivity"] = FieldSet(Field(0.0, 6.0)) + element["temperature load"] = FieldSet(Field(0.0, [12.0, 12.0, 12.0, 12.0])) + element["density"] = FieldSet(Field(0.0, 36.0)) + + # boundary element + boundary_element = Seg2([1, 2]) + boundary_element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [1.0, 0.0]])) + # linear ramp from 1 to 6 in time 0 to 1 + boundary_element["temperature flux"] = FieldSet(Field[Field(0.0, 0.0), Field(1.0, 6.0)]) + + # Set constant source f=12 with k=6. Accurate solution is + # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x) + equation = DC2D4(element) + la = initialize_local_assembly() + calculate_local_assembly!(la, equation, "temperature") + fdofs = [1, 2] + A = la.stiffness_matrix + b = la.force_vector + @fact A[fdofs, fdofs] \ b[fdofs] --> roughly([1.0, 1.0]) + + # Set constant flux g=6 on boundary. Accurate solution is + # u(x,y) = x which equals T=1 on boundary. + boundary_equation = DC2D2(boundary_element); + + calculate_local_assembly!(la, boundary_equation, "temperature") + b = la.force_vector + @fact A[fdofs, fdofs] \ b[fdofs] --> roughly([1.0, 1.0]) + +end + diff --git a/test/test_interpolation.jl b/test/test_interpolation.jl index 1b4dade..ef39013 100644 --- a/test/test_interpolation.jl +++ b/test/test_interpolation.jl @@ -4,28 +4,17 @@ using JuliaFEM: get_basis, grad, FieldSet, Field, Quad4 using FactCheck -element = Quad4([1, 2, 3, 4]) - -geometry_field = Field(0.0, Vector[]) # Create empty field at time t=0.0 -push!(geometry_field, [ 0.0, 0.0]) # push some values for field -push!(geometry_field, [ 1.0, 0.0]) -push!(geometry_field, [ 1.0, 1.0]) -push!(geometry_field, [ 0.0, 1.0]) -geometry_fieldset = FieldSet("geometry") # create fieldset "geometry" -push!(geometry_fieldset, geometry_field) # add field to fieldset -push!(element, geometry_fieldset) # add fieldset to element - -temperature_fieldset = FieldSet("temperature") -push!(temperature_fieldset, Field(0.0, [0.0, 0.0, 0.0, 0.0])) -push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0])) -push!(element, temperature_fieldset) - -displacement_fieldset = FieldSet("displacement") -push!(displacement_fieldset, Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])) -push!(displacement_fieldset, Field(1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])) -push!(element, displacement_fieldset) facts("basic continuum interpolations") do + + element = Quad4([1, 2, 3, 4]) + + element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] + element["temperature"] = ([0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 3.0, 4.0]) + element["displacement"] = ( + Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]], + Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]) + # from my old home works basis = get_basis(element) dbasis = grad(basis) diff --git a/test/test_solver.jl b/test/test_solver.jl new file mode 100644 index 0000000..e8f4156 --- /dev/null +++ b/test/test_solver.jl @@ -0,0 +1,71 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# test SimpleSolver + +using FactCheck +using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis + +""" Define Problem 1: + +- Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]} +- Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁ +""" +function get_heatproblem() + el1 = Quad4([1, 2, 3, 4]) + # these might look like normal values but believe me, they + # are fields with temporal and spatial dimension + el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] + el1["temperature thermal conductivity"] = 6.0 + el1["density"] = 36.0 + + el2 = Seg2([1, 2]) + el2["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]] + # Boundary load, linear ramp 0 -> 600 at time 0 -> 1 + # yet another simplification, if field is given as a tuple, + # multiple fields are created. there is 1 second time step between + # each field. So the following is basically same as + # fieldset = FieldSet("temperature flux") + # field1 = Field(0.0, 0.0) + # field2 = Field(1.0, 600.0) + # push!(fieldset, field1) + # push!(fieldset, field2) + # element["temperature flux"] = fieldset + el2["temperature flux"] = (0.0, 600.0) + + problem1 = PlaneHeatProblem() + push!(problem1, el1) + push!(problem1, el2) + return problem1 +end + +""" Define Problem 2: + - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂ +""" +function get_boundaryproblem() + el3 = Seg2([3, 4]) + el3["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]] + problem2 = DirichletProblem(1) + push!(problem2, el3) + return problem2 +end + +facts("test simplesolver") do + problem1 = get_heatproblem() + problem2 = get_boundaryproblem() + # Create a solver for a set of problems + solver = SimpleSolver() + push!(solver, problem1) + push!(solver, problem2) + # Solve problem at time t=1.0 and update fields + call(solver, 1.0) + # Postprocess. + # Interpolate temperature field along boundary of Γ₁ at time t=1.0 + xi = [0.0, -1.0] + el2 = get_element(problem1.equations[2]) + basis = get_basis(el2) + X = basis("geometry", xi, 1.0) + T = basis("temperature", xi, 1.0) + Logging.info("Temperature at point X = $X is T = $T") + @fact T --> roughly(100.0) +end