From 095ff09a4827fabd564612f02e70c5b4a8dc2bc4 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 7 Oct 2015 09:32:38 +0300 Subject: [PATCH] mixed problem solution --- .../2015-08-29-developing-juliafem.ipynb | 501 +++++++++++++----- 1 file changed, 380 insertions(+), 121 deletions(-) diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb index bb8ca11..3c9077f 100644 --- a/notebooks/2015-08-29-developing-juliafem.ipynb +++ b/notebooks/2015-08-29-developing-juliafem.ipynb @@ -188,14 +188,14 @@ "name": "stderr", "output_type": "stream", "text": [ - "30-Sep 09:08:39:INFO:root:Testing element MyQuad4\n", - "30-Sep 09:08:39:INFO:root:number of basis functions in this element: 4\n", - "30-Sep 09:08:39:INFO:root:Initializing element\n", - "30-Sep 09:08:40:INFO:root:Element dimension: 2\n", - "30-Sep 09:08:40:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n", - "30-Sep 09:08:40:INFO:root:Pushing field to element.\n", - "30-Sep 09:08:40:INFO:root:Interpolating scalar field at [0.0,0.0]\n", - "30-Sep 09:08:40:INFO:root:Value: 2.5\n" + "07-Oct 09:28:37:INFO:root:Testing element MyQuad4\n", + "07-Oct 09:28:37:INFO:root:number of basis functions in this element: 4\n", + "07-Oct 09:28:37:INFO:root:Initializing element\n", + "07-Oct 09:28:37:INFO:root:Element dimension: 2\n", + "07-Oct 09:28:37:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n", + "07-Oct 09:28:37:INFO:root:Pushing field to element.\n", + "07-Oct 09:28:38:INFO:root:Interpolating scalar field at [0.0,0.0]\n", + "07-Oct 09:28:38:INFO:root:Value: 2.5\n" ] }, { @@ -244,7 +244,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "30-Sep 09:08:40:INFO:root:Element MyQuad4 passed tests.\n" + "07-Oct 09:28:38:INFO:root:Element MyQuad4 passed tests.\n" ] } ], @@ -762,20 +762,22 @@ "name": "stderr", "output_type": "stream", "text": [ - "30-Sep 09:08:43:DEBUG:root:Problem (matrix) dimension: 4\n" + "07-Oct 09:28:41:DEBUG:root:total dofs: 4\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dofmap: Dict(4=>[3],2=>[1],3=>[2],5=>[4])\n" ] }, { "data": { "text/plain": [ - "(\n", - "4x4 Array{Float64,2}:\n", - " 4.0 -1.0 -2.0 -1.0\n", - " -1.0 4.0 -1.0 -2.0\n", - " -2.0 -1.0 4.0 -1.0\n", - " -1.0 -2.0 -1.0 4.0,\n", - "\n", - "[300.0,300.0,0.0,0.0])" + "4x1 sparse matrix with 2 Float64 entries:\n", + "\t[1, 1] = 300.0\n", + "\t[2, 1] = 300.0" ] }, "execution_count": 23, @@ -784,13 +786,14 @@ } ], "source": [ - "using JuliaFEM: set_global_dofs!, get_global_dofs, add_element!, get_equations, get_matrix_dimension\n", + "using JuliaFEM: get_connectivity, set_global_dofs!, get_global_dofs\n", + "using JuliaFEM: add_element!, get_equations, get_matrix_dimension\n", "\n", "# create elements and add necessary properties like connectivity and geometry\n", - "el1 = Quad4([1, 2, 3, 4])\n", + "el1 = Quad4([2, 3, 4, 5])\n", "new_field!(el1, :Geometry, Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n", "new_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n", - "el2 = Seg2([1, 2])\n", + "el2 = Seg2([2, 3])\n", "new_field!(el2, :Geometry, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n", "new_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n", "\n", @@ -798,25 +801,92 @@ "add_element!(problem, el1)\n", "add_element!(problem, el2)\n", "\n", - "# set global dofs for equations\n", - "set_global_dofs!(problem)\n", - "\n", - "n = get_matrix_dimension(problem)\n", - "\n", - "# integrate and assembly\n", - "t = 1.0\n", - "A = zeros(n, n)\n", - "b = zeros(n)\n", - "for eq in get_equations(problem)\n", - " dofs = get_global_dofs(eq)\n", - " if has_lhs(eq)\n", - " A[dofs, dofs] += integrate_lhs(eq, t)\n", + "function JuliaFEM.get_connectivity(pr::Problem)\n", + " conn = Int[]\n", + " for eq in get_equations(pr)\n", + " el = get_element(eq)\n", + " append!(conn, get_connectivity(el))\n", " end\n", - " if has_rhs(eq)\n", - " b[dofs] += integrate_rhs(eq, t)\n", + " conn = unique(conn)\n", + " return conn\n", + "end\n", + "\n", + "\"\"\"\n", + "Calculate global dofs for equations, maybe using some bandwidth\n", + "minimizing or fill reducing algorithm\n", + "\"\"\"\n", + "function calculate_global_dofs(pr::Problem)\n", + " conn = get_connectivity(pr)\n", + " dim = get_dimension(typeof(pr))\n", + " ndofs = dim*length(conn)\n", + " Logging.debug(\"total dofs: $ndofs\")\n", + "\n", + " mconn = maximum(conn)\n", + " gdofs = reshape(collect(1:mconn), dim, mconn)\n", + " dofmap = Dict{Int64, Array{Int64, 1}}()\n", + " for (i, c) in enumerate(conn)\n", + " dofmap[c] = gdofs[:, i]\n", + " end\n", + " return dofmap\n", + "end\n", + "\n", + "\"\"\"\n", + "Assign global dofs for equations.\n", + "\"\"\"\n", + "function assign_global_dofs!(pr::Problem, dofmap)\n", + " for eq in get_equations(pr)\n", + " el = get_element(eq)\n", + " c = get_connectivity(el)\n", + " #gdofs = [dofmap[ci] for ci in c]\n", + " gdofs = Int64[]\n", + " for ci in c\n", + " append!(gdofs, dofmap[ci])\n", + " end\n", + " set_global_dofs!(eq, gdofs)\n", " end\n", "end\n", - "A, b" + "\n", + "dofmap = calculate_global_dofs(problem)\n", + "println(\"dofmap: $dofmap\")\n", + "assign_global_dofs!(problem, dofmap)\n", + "\n", + "function get_lhs(pr::Problem, t::Float64)\n", + " I = Int64[]\n", + " J = Int64[]\n", + " V = Float64[]\n", + " dim = get_dimension(typeof(pr))\n", + " for eq in filter(has_lhs, get_equations(pr))\n", + " dofs = get_global_dofs(eq)\n", + " lhs = integrate_lhs(eq, t)\n", + " for (li, i) in enumerate(dofs)\n", + " for (lj, j) in enumerate(dofs)\n", + " push!(I, i)\n", + " push!(J, j)\n", + " push!(V, lhs[li, lj])\n", + " end\n", + " end\n", + " end\n", + " return I, J, V\n", + "end\n", + "\n", + "function get_rhs(pr::Problem, t::Float64)\n", + " I = Int64[]\n", + " V = Float64[]\n", + " dim = get_dimension(typeof(pr))\n", + " for eq in filter(has_rhs, get_equations(pr))\n", + " dofs = get_global_dofs(eq)\n", + " rhs = integrate_rhs(eq, t)\n", + " for (li, i) in enumerate(dofs)\n", + " push!(I, i)\n", + " push!(V, rhs[li])\n", + " end\n", + " end\n", + " return I, V\n", + "end\n", + "\n", + "t = 1.0\n", + "A = sparse(get_lhs(problem, t)...)\n", + "b = sparsevec(get_rhs(problem, t)..., size(A, 1))" ] }, { @@ -829,9 +899,9 @@ { "data": { "text/plain": [ - "2-element Array{Float64,1}:\n", - " 100.0\n", - " 100.0" + "2x1 sparse matrix with 2 Float64 entries:\n", + "\t[1, 1] = 100.0\n", + "\t[2, 1] = 100.0" ] }, "execution_count": 24, @@ -864,7 +934,7 @@ { "data": { "text/plain": [ - "InterfaceProblem" + "DirichletProblem" ] }, "execution_count": 25, @@ -873,15 +943,17 @@ } ], "source": [ - "type InterfaceProblem <: Problem\n", + "abstract BoundaryProblem <: Problem\n", + "\n", + "type DirichletProblem <: BoundaryProblem\n", " equations :: Array{Equation, 1}\n", "end\n", - "InterfaceProblem() = InterfaceProblem([])" + "DirichletProblem() = DirichletProblem([])" ] }, { "cell_type": "code", - "execution_count": 48, + "execution_count": 26, "metadata": { "collapsed": false }, @@ -892,7 +964,7 @@ "get_equation (generic function with 4 methods)" ] }, - "execution_count": 48, + "execution_count": 26, "metadata": {}, "output_type": "execute_result" } @@ -910,110 +982,62 @@ " global_dofs :: Array{Int64, 1}\n", " fieldval :: Function\n", "end\n", - "#function DBC2D2(el::Seg2, args...)\n", "function DBC2D2(el::Seg2)\n", " integration_points = [\n", " IntegrationPoint([-sqrt(1/3)], 1.0),\n", " IntegrationPoint([+sqrt(1/3)], 1.0)]\n", " new_field!(el, \"reaction force\")\n", " fieldval(X, t) = 0.0\n", - "# if length(args) != 0\n", - "# for (k, v) in args\n", - "# if k == :fieldval\n", - "# fieldval = v\n", - "# end\n", - "# end\n", - "# end\n", " DBC2D2(el, integration_points, [], fieldval)\n", "end\n", "\n", "function JuliaFEM.get_lhs(eq::DBC2D2, ip, t)\n", " el = get_element(eq)\n", " h = get_basis(el)(ip.xi)\n", - " println(\"basis h: $h\")\n", - " out = h*h'\n", - " return out\n", + " return h*h'\n", "end\n", "\n", "function JuliaFEM.get_rhs(eq::DBC2D2, ip, t)\n", " el = get_element(eq)\n", " h = get_basis(el)\n", - " #f = eq.fieldval\n", - " #X = interpolate(el, :Geometry, ip.xi, t)\n", - " #return h(ip.xi, t)*f(X, t)\n", - " return [0.0, 0.0]\n", + " f = eq.fieldval\n", + " X = interpolate(el, :Geometry, ip.xi, t)\n", + " return h(ip.xi)*f(X, t)\n", "end\n", "JuliaFEM.has_lhs(eq::DBC2D2) = true\n", "JuliaFEM.has_rhs(eq::DBC2D2) = true\n", - "JuliaFEM.get_dimension(pr::Type{InterfaceProblem}) = 1\n", - "JuliaFEM.get_equation(pr::Type{InterfaceProblem}, el::Type{Seg2}) = DBC2D2" + "JuliaFEM.get_dimension(pr::Type{DirichletProblem}) = 1\n", + "JuliaFEM.get_equation(pr::Type{DirichletProblem}, el::Type{Seg2}) = DBC2D2" ] }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { - "name": "stderr", - "output_type": "stream", - "text": [ - "30-Sep 09:15:17:DEBUG:root:Problem (matrix) dimension: 2\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "lhs\n" - ] - }, - { - "ename": "LoadError", - "evalue": "LoadError: MethodError: `call` has no method matching call(::JuliaFEM.Basis, ::Array{Float64,1}, ::Float64)\nClosest candidates are:\n BoundsError(!Matched::Any...)\n TypeVar(!Matched::Any...)\n TypeConstructor(!Matched::Any...)\n ...\nwhile loading In[49], in expression starting on line 17", - "output_type": "error", - "traceback": [ - "LoadError: MethodError: `call` has no method matching call(::JuliaFEM.Basis, ::Array{Float64,1}, ::Float64)\nClosest candidates are:\n BoundsError(!Matched::Any...)\n TypeVar(!Matched::Any...)\n TypeConstructor(!Matched::Any...)\n ...\nwhile loading In[49], in expression starting on line 17", - "", - " in get_lhs at In[48]:32", - " in integrate at /home/jukka/.julia/v0.4/JuliaFEM/src/equations.jl:69", - " in integrate_lhs at /home/jukka/.julia/v0.4/JuliaFEM/src/equations.jl:36", - " [inlined code] from In[49]:21", - " in anonymous at no file:0" - ] + "data": { + "text/plain": [ + "1-element Array{JuliaFEM.Equation,1}:\n", + " DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(:Geometry=>JuliaFEM.Field[JuliaFEM.Field{Array{Array{T,1},1}}(0.0,1,Array{T,1}[[0.0,0.0],[0.0,1.0]])],symbol(\"reaction force\")=>JuliaFEM.Field[])),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Any,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Any,Any}())],Int64[],fieldval)" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ "# create elements and add necessary properties like connectivity and geometry\n", - "el3 = Seg2([1, 2])\n", - "new_field!(el2, :Geometry, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n", + "el3 = Seg2([4, 5])\n", + "new_field!(el3, :Geometry, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n", "\n", - "bc1 = InterfaceProblem()\n", - "add_element!(bc1, el3)\n", - "\n", - "# set global dofs for equations\n", - "set_global_dofs!(bc1)\n", - "\n", - "n = get_matrix_dimension(bc1)\n", - "\n", - "# integrate and assembly\n", - "t = 1.0\n", - "A = zeros(n, n)\n", - "b = zeros(n)\n", - "for eq in get_equations(bc1)\n", - " dofs = get_global_dofs(eq)\n", - " println(\"lhs\")\n", - " if has_lhs(eq)\n", - " A[dofs, dofs] += integrate_lhs(eq, t)\n", - " end\n", - " println(\"rhs\")\n", - " if has_rhs(eq)\n", - " b[dofs] += integrate_rhs(eq, t)\n", - " end\n", - "end\n", - "A2, b2" + "bc1 = DirichletProblem()\n", + "#get_equation(typeof(bc1), typeof(el3))\n", + "#DBC2D2(el3)\n", + "add_element!(bc1, el3)" ] }, { @@ -1026,7 +1050,9 @@ { "data": { "text/plain": [ - "DBC2D2" + "4x1 sparse matrix with 2 Float64 entries:\n", + "\t[3, 1] = 0.0\n", + "\t[4, 1] = 0.0" ] }, "execution_count": 28, @@ -1035,7 +1061,21 @@ } ], "source": [ - "JuliaFEM.get_equation(typeof(bc1), typeof(el3))" + "assign_global_dofs!(bc1, dofmap)\n", + "\n", + "# integrate and assembly\n", + "t = 1.0\n", + "A2 = sparse(get_lhs(bc1, t)...)\n", + "b2 = sparsevec(get_rhs(bc1, t)...)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false + }, + "source": [ + "Now we have two problems defined, " ] }, { @@ -1048,7 +1088,11 @@ { "data": { "text/plain": [ - "DBC2D2(JuliaFEM.Seg2([1,2],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.Field[])),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Any,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Any,Any}())],Int64[],fieldval)" + "4x4 Array{Float64,2}:\n", + " 4.0 -1.0 -2.0 -1.0\n", + " -1.0 4.0 -1.0 -2.0\n", + " -2.0 -1.0 4.0 -1.0\n", + " -1.0 -2.0 -1.0 4.0" ] }, "execution_count": 29, @@ -1057,7 +1101,7 @@ } ], "source": [ - "DBC2D2(el3)" + "full(A)" ] }, { @@ -1070,7 +1114,11 @@ { "data": { "text/plain": [ - "JuliaFEM.Seg2([1,2],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.Field[]))" + "4x1 Array{Float64,2}:\n", + " 300.0\n", + " 300.0\n", + " 0.0\n", + " 0.0" ] }, "execution_count": 30, @@ -1079,8 +1127,219 @@ } ], "source": [ - "el3" + "full(b)" ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "4x4 Array{Float64,2}:\n", + " 0.0 0.0 0.0 0.0 \n", + " 0.0 0.0 0.0 0.0 \n", + " 0.0 0.0 0.333333 0.166667\n", + " 0.0 0.0 0.166667 0.333333" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "full(A2)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "4x1 Array{Float64,2}:\n", + " 0.0\n", + " 0.0\n", + " 0.0\n", + " 0.0" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "full(b2)" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8x8 sparse matrix with 24 Float64 entries:\n", + "\t[1, 1] = 4.0\n", + "\t[2, 1] = -1.0\n", + "\t[3, 1] = -2.0\n", + "\t[4, 1] = -1.0\n", + "\t[1, 2] = -1.0\n", + "\t[2, 2] = 4.0\n", + "\t[3, 2] = -1.0\n", + "\t[4, 2] = -2.0\n", + "\t[1, 3] = -2.0\n", + "\t[2, 3] = -1.0\n", + "\tâ‹®\n", + "\t[8, 3] = 0.166667\n", + "\t[1, 4] = -1.0\n", + "\t[2, 4] = -2.0\n", + "\t[3, 4] = -1.0\n", + "\t[4, 4] = 4.0\n", + "\t[7, 4] = 0.166667\n", + "\t[8, 4] = 0.333333\n", + "\t[3, 7] = 0.333333\n", + "\t[4, 7] = 0.166667\n", + "\t[3, 8] = 0.166667\n", + "\t[4, 8] = 0.333333" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Atot = [A A2; A2 zeros(A2)]" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8x1 sparse matrix with 4 Float64 entries:\n", + "\t[1, 1] = 300.0\n", + "\t[2, 1] = 300.0\n", + "\t[7, 1] = 0.0\n", + "\t[8, 1] = 0.0" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "btot = [b; b2]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Problem is that now are total matrix Atot has zero rows which needs to be removed." + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8x8 Array{Float64,2}:\n", + " 4.0 -1.0 -2.0 -1.0 0.0 0.0 0.0 0.0 \n", + " -1.0 4.0 -1.0 -2.0 0.0 0.0 0.0 0.0 \n", + " -2.0 -1.0 4.0 -1.0 0.0 0.0 0.333333 0.166667\n", + " -1.0 -2.0 -1.0 4.0 0.0 0.0 0.166667 0.333333\n", + " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", + " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", + " 0.0 0.0 0.333333 0.166667 0.0 0.0 0.0 0.0 \n", + " 0.0 0.0 0.166667 0.333333 0.0 0.0 0.0 0.0 " + ] + }, + "execution_count": 36, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "full(Atot)" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Unique rows: [1,2,3,4,7,8]\n" + ] + }, + { + "data": { + "text/plain": [ + "8x1 Array{Float64,2}:\n", + " 100.0\n", + " 100.0\n", + " 0.0\n", + " 0.0\n", + " 0.0\n", + " 0.0\n", + " 600.0\n", + " 600.0" + ] + }, + "execution_count": 37, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "r = unique(rowvals(Atot))\n", + "println(\"Unique rows: $r\")\n", + "xtot = zeros(btot)\n", + "F = lufact(Atot[r,r])\n", + "s = full(btot[r])\n", + "xtot[r] = F \\ s\n", + "full(xtot)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] } ], "metadata": {