mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-10 13:17:42 +00:00
Merge branch 'master' of git://github.com/JuliaFEM/JuliaFEM.jl
This commit is contained in:
Olli Väinölä
@@ -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 @@
|
||||
"</Grid>\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": {
|
||||
|
||||
+17
-8
@@ -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)
|
||||
|
||||
+4
-4
@@ -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')
|
||||
|
||||
|
||||
Reference in New Issue
Block a user