From d9b9edd6dd26976653b29fc59605d2e56d7c0dd3 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 2 Dec 2015 08:36:17 +0200 Subject: [PATCH 1/2] little tweaking of assembly --- src/assembly.jl | 25 +++++++++++++++++-------- src/directsolver.jl | 8 ++++---- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/src/assembly.jl b/src/assembly.jl index 39061b6..b1dfb92 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -47,9 +47,10 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int}) # empty assembly to release memory for factorization empty!(assembly.stiffness_matrix) empty!(assembly.force_vector) + gc() if dim < 100000 - # no need to do any reduction of matrix size at all + # no need to do any reduction of matrix size at all, just \ it. return CAssembly([], all_dofs, Matrix{Float64}(), K, f, spzeros(0, 0), spzeros(0,1)) end @@ -64,9 +65,10 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int}) fb = f[boundary_dofs] F = cholfact(K[interior_dofs, interior_dofs]) - K = spzeros(0, 0) + K = 0 + gc() + -#= if dim < 100000 # for small problems we don't need to care about memory usage Kd = Kib' * (F \ Kib) @@ -78,31 +80,38 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int}) for bi in 1:nb mod(bi, p) == 0 && info("Reduction: ", round(Int, bi/nb*100), " % done") C = full(F \ Kib[:, bi]) - for bj in 1:nb + for bj in bi:nb d = Kib[:, bj] - Kd[bj,bi] = dot(C[rowvals(d)], nonzeros(d)) + @inbounds Kd[bj,bi] = dot(C[rowvals(d)], nonzeros(d)) end end + Kd += tril(Kd, -1)' end Kc = spzeros(dim, dim) Kc[boundary_dofs, boundary_dofs] = Kbb - Kd -=# +#= # this is slightly faster but uses more memory chunks = round(Int, dim/3000) info("Reduction is done in $chunks chunks.") nb = length(boundary_dofs) kk = round(Int, collect(linspace(0, nb, chunks+1))) sl = [kk[j]+1:kk[j+1] for j=1:length(kk)-1] - Kc = spzeros(dim, dim) + Kd = zeros(Float64, nb, nb) + #Kd = SharedArray(Float64, nb, nb) for (k,sli) in enumerate(sl) b1 = boundary_dofs[sli] Sc = F \ Kib[:,sli] for slj in sl b2 = boundary_dofs[slj] - Kc[b2,b1] = Kbb[slj,sli] - Kib[:,slj]'*Sc + #Kc[b2,b1] = Kbb[slj,sli] - Kib[:,slj]'*Sc + Kd[slj, sli] = Kib[:,slj]'*Sc end info("Reduction: ", round(k/chunks*100, 0), " % done") end + + Kc = spzeros(dim, dim) + Kc[boundary_dofs, boundary_dofs] = Kbb - Kd +=# fc = spzeros(dim, 1) fc[boundary_dofs] = fb - Kib' * (F \ fi) diff --git a/src/directsolver.jl b/src/directsolver.jl index e32838f..4148dfb 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -163,16 +163,16 @@ function call(solver::DirectSolver, ::Type{Val{:noreduce}}, time::Number=0.0) tic(timing, "solution of system") # solve increment for linearized problem nz = unique(rowvals(A)) # take only non-zero rows - sol = zeros(b) - sol[nz] = A[nz,nz] \ full(b[nz]) + sol = zeros(length(b)) + sol[nz] = full(A[nz,nz]) \ full(b[nz]) #info("solution vector before reconstruction") #info(full(sol)') la = sol[dim+1:end] #for assembly in assemblies # reconstruct!(assembly, sol) #end - la = vec(full(la)) - sol = vec(full(sol)) +# la = vec(full(la)) +# sol = vec(full(sol)) #info("la = ", la') #info("sol = ", sol') From 071307a8f96e96819f2a6d5fe18aea09bc9bea15 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 2 Dec 2015 08:36:53 +0200 Subject: [PATCH 2/2] elasticity notebook --- ...2015-06-25-elasticity-solver-example.ipynb | 711 +++++++++++++----- 1 file changed, 514 insertions(+), 197 deletions(-) diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index 49a08bc..076f47c 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -59,11 +59,14 @@ "cell_type": "code", "execution_count": 1, "metadata": { - "collapsed": false + "collapsed": true }, "outputs": [], "source": [ - "using JuliaFEM" + "using JuliaFEM\n", + "using JuliaFEM.Core: Element, Tri3, Tet4, Tri6, Tet10\n", + "using JuliaFEM.Core: ElasticityProblem, DirichletProblem\n", + "using JuliaFEM.Core: DirectSolver" ] }, { @@ -79,38 +82,55 @@ "text": [ "INFO: Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\",\"ELSET\"]\n", "INFO: Parsing elements\n", - "INFO: 31437 elements found\n", - "INFO: Creating ELSET PISTON\n" + "INFO: 216652 elements found\n", + "INFO: Creating ELSET PISTON\n", + "INFO: Parsing elements\n", + "INFO: 24783 elements found\n" ] }, { - "data": { - "text/plain": [ - "Dict{Any,Any} with 3 entries:\n", - " \"nodes\" => Dict{Any,Any}(2843=>[-27.41338,-2.18093,-16.93988],4495=>[-14.0…\n", - " \"elements\" => Dict{Any,Any}(43367=>[4492,4602,4494,6450],35510=>[2796,2798,80…\n", - " \"elsets\" => Dict{Any,Any}(\"BC3\"=>[3975,3976,3977,3978,3979,3980,3981,3982,3…" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + " " + ] }, { "name": "stderr", "output_type": "stream", "text": [ - "INFO: Parsing elements\n", - "INFO: 5894 elements found\n", "INFO: Creating element set BC1\n", "INFO: Creating element set BC2\n", - "INFO: Creating element set BC3\n" + "INFO: Creating element set BC3\n", + "INFO: model loaded.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.760825 seconds (44.46 M allocations: 1.301 GB, 30.79% gc time)\n" ] } ], "source": [ - "model = open(JuliaFEM.parse_abaqus, \"../geometry/piston/piston_8789_P1.inp\")\n", - "model" + "@time begin\n", + " # Linear models\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_8789_P1.inp\")\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_16436_P1.inp\")\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_27343_P1.inp\")\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_45510_P1.inp\")\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_75470_P1.inp\")\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/wrench/wrench_128903_P1.inp\")\n", + "\n", + " # Quadratic models\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_19611_P2.inp\")\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_55950_P2.inp\")\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_107168_P2.inp\")\n", + " #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_345757_P2.inp\")\n", + "\n", + " info(\"model loaded.\")\n", + "end" ] }, { @@ -119,12 +139,207 @@ "metadata": { "collapsed": false }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: 241435 elements created.\n", + "INFO: all ready for solver\n" + ] + } + ], "source": [ "perm = Dict{Int64, Int64}()\n", "for (i, nid) in enumerate(keys(model[\"nodes\"]))\n", " perm[nid] = i\n", - "end" + "end\n", + "\n", + "elements = Dict{Int64, Element}()\n", + "for (elid, node_ids) in model[\"elements\"]\n", + " connectivity = Int64[perm[nid] for nid in node_ids]\n", + " coords = Vector{Float64}[model[\"nodes\"][nid] for nid in node_ids]\n", + " elmap = Dict(3 => Tri3, 4 => Tet4, 6 => Tri6, 10 => Tet10)\n", + " for (dim, eltype) in elmap\n", + " if length(coords) == dim\n", + " element = eltype(connectivity)\n", + " element[\"geometry\"] = coords\n", + " element[\"id\"] = elid\n", + " elements[elid] = element\n", + " end\n", + " end\n", + "end\n", + "\n", + "info(\"$(length(elements)) elements created.\")\n", + "\n", + "problem = ElasticityProblem()\n", + "#for elid in model[\"elsets\"][\"WRENCH\"]\n", + "for elid in model[\"elsets\"][\"PISTON\"]\n", + " elements[elid][\"youngs modulus\"] = 210.0e3\n", + " elements[elid][\"poissons ratio\"] = 0.3\n", + " push!(problem, elements[elid])\n", + "end\n", + "\n", + "traction = Vector{Float64}[[0.0, 0.0, -10.0] for i in 1:6]\n", + "#traction = Vector{Float64}[[100.0, 100.0, 100.0] for i in 1:3]\n", + "#traction = Vector{Float64}[[0.0, -100.0, 20.0] for i in 1:3]\n", + "for elid in model[\"elsets\"][\"BC1\"]\n", + " elements[elid][\"displacement traction force\"] = traction\n", + " push!(problem, elements[elid])\n", + "end\n", + "\n", + "bc = DirichletProblem(\"displacement\", 3)\n", + "for elid in model[\"elsets\"][\"BC2\"]\n", + " elements[elid][\"displacement\"] = 0.0\n", + " push!(bc, elements[elid])\n", + "end\n", + "\n", + "solver = DirectSolver()\n", + "push!(solver, problem)\n", + "push!(solver, bc)\n", + "info(\"all ready for solver\")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "empty!(model)\n", + "gc()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: # of field problems: 1\n", + "INFO: # of boundary problems: 1\n", + "INFO: Starting iteration 1\n", + "INFO: Assembling boundary problems...\n", + "INFO: Assemble: 10 % done\n", + "INFO: Assemble: 20 % done\n", + "INFO: Assemble: 30 % done\n", + "INFO: Assemble: 40 % done\n", + "INFO: Assemble: 50 % done\n", + "INFO: Assemble: 60 % done\n", + "INFO: Assemble: 70 % done\n", + "INFO: Assemble: 80 % done\n", + "INFO: Assemble: 90 % done\n", + "INFO: Assemble: 100 % done\n", + "INFO: # of interface dofs: 19848\n", + "INFO: Assembling field problems...\n", + "INFO: Assembling body 1...\n", + "INFO: Assemble: 10 % done\n" + ] + }, + { + "ename": "LoadError", + "evalue": "LoadError: InterruptException:\nwhile loading In[5], in expression starting on line 1", + "output_type": "error", + "traceback": [ + "LoadError: InterruptException:\nwhile loading In[5], in expression starting on line 1", + "" + ] + } + ], + "source": [ + "call(solver, 0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: Assemble: 10 % done\n", + "INFO: Assemble: 20 % done\n", + "INFO: Assemble: 30 % done\n", + "INFO: Assemble: 40 % done\n", + "INFO: Assemble: 50 % done\n", + "INFO: Assemble: 60 % done\n", + "INFO: Assemble: 70 % done\n", + "INFO: Assemble: 80 % done\n", + "INFO: Assemble: 90 % done\n", + "INFO: Assemble: 100 % done\n", + "INFO: 1350 boundary dofs\n", + "INFO: Assemble: 10 % done\n", + "INFO: Assemble: 20 % done\n", + "INFO: Assemble: 30 % done\n", + "INFO: Assemble: 40 % done\n", + "INFO: Assemble: 50 % done\n", + "INFO: Assemble: 60 % done\n", + "INFO: Assemble: 70 % done\n", + "INFO: Assemble: 80 % done\n", + "INFO: Assemble: 90 % done\n", + "INFO: Assemble: 100 % done\n", + "INFO: 26367 field dofs\n", + "INFO: 25017 interior dofs\n", + "INFO: Assembly: 17.84435486793518\n", + "INFO: det = Inf\n", + "INFO: Factorization: 1.87113618850708\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 19.775338 seconds (152.31 M allocations: 8.037 GB, 10.16% gc time)\n" + ] + } + ], + "source": [ + "using JuliaFEM.Core: assemble\n", + "\n", + "for elid in model[\"elsets\"][\"PISTON\"]\n", + " elements[elid][\"displacement\"] = (0.0 => Vector{Float64}[zeros(Float64, 3) for i = 1:4])\n", + "end\n", + "\n", + "function do_assembly()\n", + " t0 = time()\n", + " boundary_assembly = sum(map((p)->assemble(p, 0.0), solver.boundary_problems))\n", + " boundary_dofs = unique(boundary_assembly.stiffness_matrix.I)\n", + " info(\"$(length(boundary_dofs)) boundary dofs\")\n", + " field_assembly = sum(map((p)->assemble(p, 0.0), solver.field_problems))\n", + " field_dofs = unique(field_assembly.stiffness_matrix.I)\n", + " info(\"$(length(field_dofs)) field dofs\")\n", + " interior_dofs = setdiff(field_dofs, boundary_dofs)\n", + " info(\"$(length(interior_dofs)) interior dofs\")\n", + " info(\"Assembly: \", time()-t0)\n", + "\n", + " t0 = time()\n", + " K = sparse(field_assembly.stiffness_matrix)\n", + " f = sparse(field_assembly.force_vector)\n", + " @assert maximum(abs(1/2*(K + K') - K)) < 1.0e-6\n", + " K = 1/2*(K + K')\n", + " Kii = K[interior_dofs, interior_dofs]\n", + " info(\"det = \", det(Kii))\n", + " F = cholfact(Kii)\n", + " Kbb = K[boundary_dofs, boundary_dofs]\n", + " Kib = K[interior_dofs, boundary_dofs]\n", + " fi = f[interior_dofs]\n", + " info(\"Factorization: \", time()-t0)\n", + " dim = size(K, 1)\n", + " return dim, F, Kbb, Kib, fi, boundary_dofs, interior_dofs\n", + "end\n", + "\n", + "@time dim, F, Kbb, Kib, fi, boundary_dofs, interior_dofs = do_assembly();" ] }, { @@ -133,34 +348,107 @@ "metadata": { "collapsed": false }, + "outputs": [ + { + "data": { + "text/plain": [ + "58833x1 sparse matrix with 0 Float64 entries:" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Kc = spzeros(dim, dim)\n", + "fc = spzeros(dim, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "metadata": { + "collapsed": false + }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "INFO: 37331 elements created.\n" + "INFO: 10 % done.\n", + "INFO: 20 % done.\n", + "INFO: 30 % done.\n", + "INFO: 40 % done.\n", + "INFO: 50 % done.\n", + "INFO: 60 % done.\n", + "INFO: 70 % done.\n", + "INFO: 80 % done.\n", + "INFO: 90 % done.\n", + "INFO: 100 % done.\n", + "INFO: norm(2) = 535881.6985084609\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 20.500665 seconds (13.36 M allocations: 3.063 GB, 3.34% gc time)\n" ] } ], "source": [ - "using JuliaFEM: Element, Tri3, Tet4\n", - "elements = Dict{Int64, Element}()\n", - "for (elid, node_ids) in model[\"elements\"]\n", - " connectivity = Int64[perm[nid] for nid in node_ids]\n", - " coords = Vector{Float64}[model[\"nodes\"][nid] for nid in node_ids]\n", - " if length(coords) == 3\n", - " element = Tri3(connectivity)\n", - " element[\"geometry\"] = coords\n", - " elements[elid] = element\n", - " elseif length(coords) == 4\n", - " element = Tet4(connectivity)\n", - " element[\"geometry\"] = coords\n", - " elements[elid] = element\n", - " else\n", - " warn(\"unknown element dim = $(length(connectivity))\")\n", + "@time begin\n", + " nb = length(boundary_dofs)\n", + " p = round(Int, nb/10)\n", + " Kd = zeros(nb, nb)\n", + " Vd = zeros(nb)\n", + " for bi in 1:nb\n", + " mod(bi, p) == 0 && info(round(Int, bi/nb*100), \" % done.\")\n", + " C = full(F \\ Kib[:, bi])\n", + " fill!(Vd, 0.0)\n", + " for bj in bi:nb\n", + " d = Kib[:, bj]\n", + " Kd[bj,bi] = dot(C[rowvals(d)], nonzeros(d))\n", + " Vd[bj] = dot(C[rowvals(d)], nonzeros(d))\n", + " end\n", + " Kd[:, bi] = Vd\n", " end\n", - "end\n", - "info(\"$(length(elements)) elements created.\")" + " Kd += tril(Kd, -1)'\n", + " Kc = spzeros(dim, dim)\n", + " Kc[boundary_dofs, boundary_dofs] = Kbb - Kd\n", + " info(\"norm(2) = \", norm(Kd))\n", + " @assert isapprox(norm(Kd), 535881.698508461)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "79868.97400906119" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "maximum(abs(1/2*(Kd + Kd') - Kd))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "30% memory usage before start. 100% CPU usage. 85 % maximum memory usage." ] }, { @@ -173,7 +461,7 @@ { "data": { "text/plain": [ - "37331" + "get_slices (generic function with 1 method)" ] }, "execution_count": 5, @@ -182,74 +470,17 @@ } ], "source": [ - "31437+5894" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "using JuliaFEM: ElasticityProblem, DirichletProblem\n", - "\n", - "piston = ElasticityProblem()\n", - "for elid in model[\"elsets\"][\"PISTON\"]\n", - " elements[elid][\"youngs modulus\"] = 210.0e3\n", - " elements[elid][\"poissons ratio\"] = 0.3\n", - " push!(piston, elements[elid])\n", + "function get_slices(nb, n)\n", + "# n = 3\n", + "# nb = 25\n", + " kk = round(Int, collect(linspace(0, nb, n+1)))\n", + " return [kk[j]+1:kk[j+1] for j=1:length(kk)-1]\n", "end" ] }, { "cell_type": "code", - "execution_count": 7, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "load = Vector{Float64}[[0.0, 0.0, -10.0], [0.0, 0.0, -10.0], [0.0, 0.0, -10.0]]\n", - "for elid in model[\"elsets\"][\"BC1\"]\n", - " elements[elid][\"displacement traction force\"] = load\n", - " push!(piston, elements[elid])\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "encastre = DirichletProblem(\"displacement\", 3)\n", - "for elid in model[\"elsets\"][\"BC2\"]\n", - " elements[elid][\"displacement\"] = 0.0\n", - " push!(encastre, elements[elid])\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "using JuliaFEM: DirectSolver\n", - "solver = DirectSolver()\n", - "push!(solver, piston)\n", - "push!(solver, encastre);" - ] - }, - { - "cell_type": "code", - "execution_count": 10, + "execution_count": 16, "metadata": { "collapsed": false }, @@ -258,64 +489,59 @@ "name": "stderr", "output_type": "stream", "text": [ - "INFO: # of field problems: 1\n", - "INFO: # of boundary problems: 1\n", - "INFO: Starting iteration 1\n", - "INFO: # of dofs: 26367, # of interface dofs: 1350\n", - "INFO: solved. length of solution vector = 52734\n", - "INFO: timing info for non-linear iteration:\n", - "INFO: boundary assembly : 1.4367611408233643\n", - "INFO: field assembly : 17.380300998687744\n", - "INFO: create sparse matrices : 0.38140082359313965\n", - "INFO: solution of system : 3.159921169281006\n", - "INFO: update element data : 3.3958511352539062\n", - "INFO: non-linear iteration : 25.75425100326538\n", - "INFO: Starting iteration 2\n", - "INFO: # of dofs: 26367, # of interface dofs: 1350\n", - "INFO: solved. length of solution vector = 52734\n", - "INFO: timing info for non-linear iteration:\n", - "INFO: boundary assembly : 0.1028587818145752\n", - "INFO: field assembly : 15.321197032928467\n", - "INFO: create sparse matrices : 0.11679887771606445\n", - "INFO: solution of system : 2.1216089725494385\n", - "INFO: update element data : 3.484022855758667\n", - "INFO: non-linear iteration : 21.146512031555176\n", - "INFO: Starting iteration 3\n", - "INFO: # of dofs: 26367, # of interface dofs: 1350\n", - "INFO: solved. length of solution vector = 52734\n", - "INFO: timing info for non-linear iteration:\n", - "INFO: boundary assembly : 0.09190106391906738\n", - "INFO: field assembly : 14.772294044494629\n", - "INFO: create sparse matrices : 0.12086606025695801\n", - "INFO: solution of system : 2.260561943054199\n", - "INFO: update element data : 3.4714930057525635\n", - "INFO: non-linear iteration : 20.717132806777954\n", - "INFO: Starting iteration 4\n", - "INFO: # of dofs: 26367, # of interface dofs: 1350\n", - "INFO: solved. length of solution vector = 52734\n", - "INFO: timing info for non-linear iteration:\n", - "INFO: boundary assembly : 0.09454894065856934\n", - "INFO: field assembly : 15.800831079483032\n", - "INFO: create sparse matrices : 0.10474300384521484\n", - "INFO: solution of system : 2.3751449584960938\n", - "INFO: update element data : 3.471247911453247\n", - "INFO: non-linear iteration : 21.846541166305542\n", - "INFO: solver finished in 90.86189484596252 seconds.\n" + "INFO: 10.0 % done\n", + "INFO: 20.0 % done\n", + "INFO: 30.0 % done\n", + "INFO: 40.0 % done\n", + "INFO: 50.0 % done\n", + "INFO: 60.0 % done\n", + "INFO: 70.0 % done\n", + "INFO: 80.0 % done\n", + "INFO: 90.0 % done\n", + "INFO: 100.0 % done\n" ] }, { - "data": { - "text/plain": [ - "(4,true)" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + " 47" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: 3.525162568997601e6\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + ".470587 seconds (251.48 k allocations: 10.678 GB, 5.06% gc time)\n" + ] } ], "source": [ - "solver(0.0)" + "@time begin\n", + " chunks = 10\n", + " nb = length(boundary_dofs)\n", + " kk = round(Int, collect(linspace(0, nb, chunks+1)))\n", + " sl = [kk[j]+1:kk[j+1] for j=1:length(kk)-1]\n", + " Kc = spzeros(dim, dim)\n", + " for (k, sli) in enumerate(sl)\n", + " b1 = boundary_dofs[sli]\n", + " Sc = F \\ Kib[:,sli]\n", + " for slj in sl\n", + " b2 = boundary_dofs[slj]\n", + " Kc[b2,b1] = Kbb[slj,sli] - Kib[:,slj]'*Sc\n", + " end\n", + " info(round(k/chunks*100, 0), \" % done\")\n", + " end\n", + "end\n", + "info(norm(Kc, 1))\n", + "@assert isapprox(norm(Kc, 1), 3.525162568997601e6)" ] }, { @@ -327,7 +553,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 20, "metadata": { "collapsed": false }, @@ -340,15 +566,15 @@ "\n" ] }, - "execution_count": 11, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "xdoc, xmodel = JuliaFEM.xdmf_new_model()\n", - "temporal_collection = JuliaFEM.xdmf_new_temporal_collection(xmodel)\n", - "grid = JuliaFEM.xdmf_new_grid(temporal_collection; time=0)" + "xdoc, xmodel = JuliaFEM.Postprocess.xdmf_new_model()\n", + "temporal_collection = JuliaFEM.Postprocess.xdmf_new_temporal_collection(xmodel)\n", + "grid = JuliaFEM.Postprocess.xdmf_new_grid(temporal_collection; time=0)" ] }, { @@ -360,7 +586,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 21, "metadata": { "collapsed": false }, @@ -395,7 +621,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 22, "metadata": { "collapsed": false }, @@ -406,35 +632,13 @@ "true" ] }, - "execution_count": 13, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "JuliaFEM.xdmf_new_mesh(grid, X, elmap)" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "988949" - ] - }, - "execution_count": 14, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "JuliaFEM.xdmf_save_model(xdoc, \"/tmp/piston.xmf\")" + "JuliaFEM.Postprocess.xdmf_new_mesh(grid, X, elmap)" ] }, { @@ -446,7 +650,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 23, "metadata": { "collapsed": false }, @@ -457,13 +661,13 @@ "(3,8789)" ] }, - "execution_count": 15, + "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "using JuliaFEM: get_connectivity\n", + "using JuliaFEM.Core: get_connectivity\n", "\n", "u = zeros(3, nnodes)\n", "\n", @@ -480,7 +684,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 24, "metadata": { "collapsed": false }, @@ -491,18 +695,18 @@ "true" ] }, - "execution_count": 16, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "JuliaFEM.xdmf_new_field(grid, \"Displacement\", \"nodes\", u)" + "JuliaFEM.Postprocess.xdmf_new_field(grid, \"Displacement\", \"nodes\", u)" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 25, "metadata": { "collapsed": false }, @@ -510,17 +714,130 @@ { "data": { "text/plain": [ - "1573751" + "1525505" ] }, - "execution_count": 17, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "JuliaFEM.xdmf_save_model(xdoc, \"/tmp/piston.xmf\")" + "JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/piston_8789_P1.xmf\")" ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "3x10 Array{Float64,2}:\n", + " 0.0582911 0.339694 0.431768 … 0.124351 0.346492 0.000617283\n", + " 0.0900323 0.22174 0.433814 0.140503 0.275036 0.0044649 \n", + " 0.0838928 -0.148224 0.00262483 -0.11171 0.0843743 -0.0175412 " + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "u[:, 1:10]" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1x3 Array{Float64,2}:\n", + " -0.0471355 -0.468786 -0.222571" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "minimum(u, 2)'" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1x3 Array{Float64,2}:\n", + " 0.511168 0.563004 0.66328" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "maximum(u, 2)'" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "38.52232721814439" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "piston_8789_P1_solution_norm = 38.52232721814439" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "@assert isapprox(norm(vec(u)), piston_8789_P1_solution_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] } ], "metadata": {