3d tie contact with multiple boundaries, shared nodes etc.

This commit is contained in:
Jukka Aho
2015-12-17 22:28:15 +02:00
parent 869c36ac70
commit d2d7b431bb
+301
View File
@@ -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": [
"<img src=\"http://results.juliafem.org/2015-12-17-superblock.png\"\n",
" style=\"width:80%; height:auto; vertical-align:middle;\">"
]
},
{
"cell_type": "code",
"execution_count": null,