diff --git a/notebooks/2015-12-13-3d-tie-contact.ipynb b/notebooks/2015-12-13-3d-tie-contact.ipynb index e93c475..45b6d8b 100644 --- a/notebooks/2015-12-13-3d-tie-contact.ipynb +++ b/notebooks/2015-12-13-3d-tie-contact.ipynb @@ -1000,6 +1000,307 @@ " style=\"width:80%; height:auto; vertical-align:middle;\">" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Several contacts with shared nodes" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: Found 12 element sets: BLOCK4_TO_BLOCK3, SYM12, SYM23, BLOCK4_TO_BLOCK2, BLOCK3_TO_BLOCK4, BLOCK1_TO_BLOCK3, SYM13, LOAD, BLOCK1_TO_BLOCK2, BLOCK2_TO_BLOCK4, BLOCK2_TO_BLOCK1, BLOCK3_TO_BLOCK1\n" + ] + }, + { + "data": { + "text/plain": [ + "Dict{ASCIIString,Any} with 2 entries:\n", + " \"nodes\" => Dict(306=>[0.0,0.8469606634303125,0.4191161978569679],29=>[…\n", + " \"connectivity\" => Dict(2843=>(:TE4,:OTHER,[631,574,595,513]),1316=>(:TR3,:SYM…" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using JuliaFEM\n", + "using JuliaFEM.Preprocess: parse_aster_med_file\n", + "using JuliaFEM.Core: LinearElasticityProblem, DirichletProblem, get_connectivity, Tri3, Tet4, LinearSolver, update\n", + "mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/SUPERBLOCK.med\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: created 2878 elements.\n", + "INFO: created 403 boundary elements.\n" + ] + } + ], + "source": [ + "# interior elements are of type HE8\n", + "field_problem = LinearElasticityProblem()\n", + "\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " eltype == :TE4 || continue\n", + " element = Tet4(elcon)\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " element[\"youngs modulus\"] = 900.0\n", + " element[\"poissons ratio\"] = 0.25\n", + " #info(\"lower: add element with connectivity $elcon\")\n", + " push!(field_problem, element)\n", + "end\n", + "\n", + "# Neumann boundary condition, traction force -100 on Z direction for element set LOAD\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :TR3) && (elset == :LOAD) || continue\n", + " element = Tri3(elcon)\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " element[\"displacement traction force 3\"] = -100.0\n", + " push!(field_problem, element)\n", + "end\n", + "info(\"created $(length(field_problem.elements)) elements.\")\n", + "\n", + "# boundary conditions\n", + "boundary_problem = DirichletProblem(\"displacement\", 3)\n", + "\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :TR3) && (elset in [:SYM23, :SYM12, :SYM13]) || continue\n", + " element = Tri3(elcon)\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " if elset == :SYM23\n", + " element[\"displacement 1\"] = 0.0\n", + " elseif elset == :SYM12\n", + " element[\"displacement 3\"] = 0.0\n", + " elseif elset == :SYM13\n", + " element[\"displacement 2\"] = 0.0\n", + " end\n", + " push!(boundary_problem, element)\n", + "end\n", + "info(\"created $(length(boundary_problem.elements)) boundary elements.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Contact pairs: BLOCK1 <-> BLOCK3, BLOCK1 <-> BLOCK2, BLOCK2 <-> BLOCK4, BLOCK3 <-> BLOCK4" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "using JuliaFEM.Core: Element, MortarProblem, calculate_normal_tangential_coordinates!\n", + "\n", + "function create_contact(master_surface, slave_surface)\n", + "\n", + " master_elements = JuliaFEM.Core.Element[]\n", + " for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :TR3) && (elset == master_surface) || continue\n", + " element = Tri3(elcon)\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " push!(master_elements, element)\n", + " end\n", + "\n", + " contact_problem = MortarProblem(\"displacement\", 3)\n", + " for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :TR3) && (elset == slave_surface) || continue\n", + " element = Tri3(reverse(elcon))\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " element[\"master elements\"] = master_elements\n", + " calculate_normal_tangential_coordinates!(element, 0.0)\n", + " push!(contact_problem, element)\n", + " end\n", + "\n", + " return contact_problem\n", + "end\n", + "\n", + "tie1 = create_contact(:BLOCK1_TO_BLOCK3, :BLOCK3_TO_BLOCK1)\n", + "tie2 = create_contact(:BLOCK1_TO_BLOCK2, :BLOCK2_TO_BLOCK1)\n", + "tie3 = create_contact(:BLOCK2_TO_BLOCK4, :BLOCK4_TO_BLOCK2)\n", + "tie4 = create_contact(:BLOCK3_TO_BLOCK4, :BLOCK4_TO_BLOCK3);" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: # of field problems: 1\n", + "INFO: # of boundary problems: 5\n", + "INFO: Starting iteration 1\n", + "INFO: Assembling field problems...\n", + "INFO: Assembling body 1...\n", + "INFO: Assembly: 10.0 % done. \n", + "INFO: Assembly: 20.0 % done. \n", + "INFO: Assembly: 30.0 % done. \n", + "INFO: Assembly: 40.0 % done. \n", + "INFO: Assembly: 50.0 % done. \n", + "INFO: Assembly: 60.0 % done. \n", + "INFO: Assembly: 70.0 % done. \n", + "INFO: Assembly: 80.0 % done. \n", + "INFO: Assembly: 90.0 % done. \n", + "INFO: Assembly: 100.0 % done. \n", + "INFO: dim = 2562\n", + "INFO: Assembling boundary problems...\n", + "INFO: Assembling boundary 1...\n", + "INFO: Assembly: 10.0 % done. \n", + "INFO: Assembly: 20.0 % done. \n", + "INFO: Assembly: 30.0 % done. \n", + "INFO: Assembly: 40.0 % done. \n", + "INFO: Assembly: 50.0 % done. \n", + "INFO: Assembly: 60.0 % done. \n", + "INFO: Assembly: 70.0 % done. \n", + "INFO: Assembly: 80.0 % done. \n", + "INFO: Assembly: 90.0 % done. \n", + "INFO: Assembly: 100.0 % done. \n", + "INFO: Assembling boundary 2...\n", + "INFO: Assembling boundary 3...\n", + "INFO: Assembling boundary 4...\n", + "INFO: Assembling boundary 5...\n", + "INFO: dumping matrices to disk, file = matrices_tie_contact_3d_host_1_iteration_1.jld\n", + "INFO: Solving system\n", + "INFO: UMFPACK: solved in 0.2858870029449463 seconds. norm = 1.9710688787586395\n", + "INFO: timing info for iteration:\n", + "INFO: boundary assembly : 10.795902967453003\n", + "INFO: field assembly : 1.356456995010376\n", + "INFO: dump matrices to disk : 1.1161930561065674\n", + "INFO: solve problem : 0.4241180419921875\n", + "INFO: update element data : 0.05763602256774902\n", + "INFO: non-linear iteration : 13.75032901763916\n", + "INFO: solver finished in 13.92406702041626 seconds.\n", + "INFO: nid near corner = 502\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777977912483327,0.02776944079303928,-0.11111675469425412]\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777977912483327,0.02776944079303928,-0.11111675469425412]\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777977912483327,0.02776944079303928,-0.11111675469425412]\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777977912483327,0.02776944079303928,-0.11111675469425412]\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777977912483327,0.02776944079303928,-0.11111675469425412]\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777977912483327,0.02776944079303928,-0.11111675469425412]\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777977912483327,0.02776944079303928,-0.11111675469425412]\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777977912483327,0.02776944079303928,-0.11111675469425412]\n" + ] + } + ], + "source": [ + "using JuliaFEM.Core: DirectSolver\n", + "solver = DirectSolver()\n", + "solver.name = \"tie_contact_3d\"\n", + "solver.method = :UMFPACK\n", + "solver.nonlinear_problem = false\n", + "solver.max_iterations = 1\n", + "solver.dump_matrices = true\n", + "push!(solver, field_problem)\n", + "push!(solver, boundary_problem)\n", + "push!(solver, tie1, tie2, tie3, tie4)\n", + "call(solver, 0.0)\n", + "\n", + "using JuliaFEM.Test\n", + "\n", + "nid = 0\n", + "for (nid, coords) in mesh[\"nodes\"]\n", + " if isapprox(coords, [1.0, 1.0, 1.0])\n", + " info(\"nid near corner = $nid\")\n", + " break\n", + " end\n", + "end\n", + "nid\n", + "\n", + "known_value = [1/36, 1/36, -1/9]\n", + "\n", + "for element in field_problem.elements\n", + " i = indexin([nid], get_connectivity(element))[1]\n", + " i != 0 || continue\n", + " X = element(\"geometry\", 0.0)\n", + " u = element(\"displacement\", 0.0)\n", + " info(\"displacement X = $(X[i]), u = $(u[i])\")\n", + " @test isapprox(u[i], known_value, atol=1.0e-3)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: XDFM: ndim = 2562\n" + ] + } + ], + "source": [ + "xdoc, xmodel = JuliaFEM.Postprocess.xdmf_new_model()\n", + "coll = JuliaFEM.Postprocess.xdmf_new_temporal_collection(xmodel)\n", + "grid = JuliaFEM.Postprocess.xdmf_new_grid(coll; time=0.0)\n", + "\n", + "Xg = Dict{Int64, Vector{Float64}}()\n", + "ug = Dict{Int64, Vector{Float64}}()\n", + "for element in field_problem.elements\n", + " conn = get_connectivity(element)\n", + " X = element(\"geometry\", 0.0)\n", + " u = element(\"displacement\", 0.0)\n", + " for (i, c) in enumerate(conn)\n", + " Xg[c] = X[i]\n", + " ug[c] = u[i]\n", + " end\n", + "end\n", + "perm = sort(collect(keys(Xg)))\n", + "nodes = Vector{Float64}[Xg[i] for i in perm]\n", + "disp = Vector{Float64}[ug[i] for i in perm]\n", + "elements = []\n", + "for el in field_problem.elements\n", + " isa(el, JuliaFEM.Core.Element{JuliaFEM.Core.Tet4}) || continue\n", + " push!(elements, (:Tet4, get_connectivity(el)))\n", + "end\n", + "#elements\n", + "JuliaFEM.Postprocess.xdmf_new_mesh!(grid, nodes, elements)\n", + "JuliaFEM.Postprocess.xdmf_new_nodal_field!(grid, \"displacement\", disp)\n", + "JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/superblock.xmf\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "" + ] + }, { "cell_type": "code", "execution_count": null,