2015-12-14 02:09:33 +02:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tie contact 3d\n",
"\n",
"Author(s): Jukka Aho"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solid block\n",
"\n",
2015-12-14 17:49:36 +02:00
"Geometry and mesh:\n",
"\n",
"<div class=\"main_block\" style=\"max-width: 800px;\">\n",
" <div class=\"inner_block\" style=\"display: inline-block; float:left; width:50%;\">\n",
" <img src=\"http://results.juliafem.org/2015-12-13-box_geom.png\"\n",
" style=\"width:100%; height:auto; vertical-align:middle;\">\n",
" </div>\n",
" <div class=\"inner_block\" style=\"display: inline-block; float:left; width:50%;\">\n",
" <img src=\"http://results.juliafem.org/2015-12-13-box_bcs.png\"\n",
" style=\"width:100%; height:auto; vertical-align:middle;\">\n",
" </div>\n",
"</div>"
2015-12-14 02:09:33 +02:00
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM\n",
"using JuliaFEM.Preprocess: parse_aster_med_file\n",
2015-12-14 17:49:36 +02:00
"using JuliaFEM.Core: LinearElasticityProblem, DirichletProblem, get_connectivity, Quad4, Hex8, DirectSolver"
2015-12-14 02:09:33 +02:00
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Found 4 element sets: SYM23, SYM12, SYM13, LOAD\n"
]
},
{
"data": {
"text/plain": [
"Dict{ASCIIString,Any} with 2 entries:\n",
" \"nodes\" => Dict(2=>[0.0,0.0,0.0],11=>[0.0,0.3333333333333333,1.0],39=>…\n",
" \"connectivity\" => Dict(68=>(:QU4,:OTHER,[45,47,48,46]),2=>(:SE2,:OTHER,[9,10]…"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/mesh.med\")\n",
"#mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/twoelem_box.med\")"
]
},
{
"cell_type": "code",
2015-12-14 17:49:36 +02:00
"execution_count": 5,
2015-12-14 02:09:33 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: created 36 elements.\n"
]
}
],
"source": [
"# interior elements are of type HE8\n",
"field_problem = LinearElasticityProblem()\n",
"for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :HE8 || continue\n",
" element = Hex8(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mesh[\"nodes\"][i] for i in get_connectivity(element)]\n",
" element[\"youngs modulus\"] = 900.0\n",
" element[\"poissons ratio\"] = 0.25\n",
" push!(field_problem, element)\n",
"end\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 == :QU4) && (elset == :LOAD) || continue\n",
" element = Quad4(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mesh[\"nodes\"][i] for i in get_connectivity(element)]\n",
" element[\"displacement traction force\"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]\n",
" push!(field_problem, element)\n",
"end\n",
"info(\"created $(length(field_problem.elements)) elements.\")"
]
},
{
"cell_type": "code",
2015-12-14 17:49:36 +02:00
"execution_count": 6,
2015-12-14 02:09:33 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: created 27 boundary elements.\n"
]
}
],
"source": [
"# boundary conditions\n",
"boundary_problem = DirichletProblem(\"displacement\", 3)\n",
"for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" (eltype == :QU4) && (elset in [:SYM23, :SYM12, :SYM13]) || continue\n",
" element = Quad4(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mesh[\"nodes\"][i] for i in get_connectivity(element)]\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": "code",
2015-12-14 17:49:36 +02:00
"execution_count": 7,
2015-12-14 02:09:33 +02:00
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
2015-12-14 17:49:36 +02:00
"INFO: # of field problems: 1\n",
"INFO: # of boundary problems: 1\n",
"INFO: Starting iteration 1\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: dim = 192\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1...\n",
"INFO: dumping matrices to disk, file = matrices_block_host_1_iteration_1.jld\n",
"INFO: Solving system\n",
"INFO: CHOLMOD: all dofs = 192\n",
"INFO: CHOLMOD: interior dofs = 144\n",
"INFO: CHOLMOD: boundary dofs = 48\n",
"INFO: CHOLMOD: displacement on boundary solved.\n",
"INFO: CHOLMOD: homogeneous dirichlet boundary\n",
"INFO: CHOLMOD: LDLt factorization done in 0.0012049674987792969 seconds\n",
"INFO: CHOLMOD: solved in 0.05479288101196289 seconds. norm = 0.5879447357921324\n",
"INFO: timing info for iteration:\n",
"INFO: boundary assembly : 0.09608697891235352\n"
2015-12-14 02:09:33 +02:00
]
},
{
"data": {
"text/plain": [
2015-12-14 17:49:36 +02:00
"(1,true)"
2015-12-14 02:09:33 +02:00
]
},
2015-12-14 17:49:36 +02:00
"execution_count": 7,
2015-12-14 02:09:33 +02:00
"metadata": {},
"output_type": "execute_result"
2015-12-14 17:49:36 +02:00
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: field assembly : 1.0157549381256104\n",
"INFO: dump matrices to disk : 0.9569120407104492\n",
"INFO: solve problem : 0.48387885093688965\n",
"INFO: update element data : 0.01578807830810547\n",
"INFO: non-linear iteration : 2.568441867828369\n",
"INFO: solver finished in 2.7200798988342285 seconds.\n"
]
2015-12-14 02:09:33 +02:00
}
],
"source": [
2015-12-14 17:49:36 +02:00
"solver = DirectSolver()\n",
"solver.name = \"block\"\n",
"solver.nonlinear_problem = false\n",
"#solver.method = :UMFPACK\n",
"solver.dump_matrices = true\n",
"push!(solver, field_problem)\n",
"push!(solver, boundary_problem)\n",
"call(solver, 0.0)"
2015-12-14 02:09:33 +02:00
]
},
{
"cell_type": "code",
2015-12-14 17:49:36 +02:00
"execution_count": 8,
2015-12-14 02:09:33 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: nid near corner = 7\n"
]
},
{
"data": {
"text/plain": [
"7"
]
},
2015-12-14 17:49:36 +02:00
"execution_count": 8,
2015-12-14 02:09:33 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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"
]
},
{
"cell_type": "code",
2015-12-14 17:49:36 +02:00
"execution_count": 9,
2015-12-14 02:09:33 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
2015-12-14 17:49:36 +02:00
"INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777794,0.02777777777777784,-0.11111111111111129]\n",
"INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777794,0.02777777777777784,-0.11111111111111129]\n"
2015-12-14 02:09:33 +02:00
]
}
],
"source": [
"using JuliaFEM.Test\n",
"known_value = [1/36, 1/36, -1/9]\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",
2015-12-14 17:49:36 +02:00
" @test isapprox(u[i], known_value)\n",
2015-12-14 02:09:33 +02:00
"end"
]
},
{
"cell_type": "code",
2015-12-14 17:49:36 +02:00
"execution_count": 10,
2015-12-14 02:09:33 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: XDFM: ndim = 192\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.Hex8}) || continue\n",
" push!(elements, (:Hex8, 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/foobar2.xmf\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<img src=\"http://results.juliafem.org/2015-12-13-box_results.png\" width=300px style=\"float:left\">"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Block divived to two parts\n",
"\n",
2015-12-14 17:49:36 +02:00
"Block is now divided to two parts and meshes are tied using mortar method.\n",
"\n",
"<div class=\"main_block\" style=\"max-width: 800px;\">\n",
" <div class=\"inner_block\" style=\"display: inline-block; float:left; width:50%;\">\n",
" <img src=\"http://results.juliafem.org/2015-12-13-divided-block-geometry.png\"\n",
" style=\"width:100%; height:auto; vertical-align:middle;\">\n",
" </div>\n",
" <div class=\"inner_block\" style=\"display: inline-block; float:left; width:50%;\">\n",
" <img src=\"http://results.juliafem.org/2015-12-13-divided-block-both-parts.png\"\n",
" style=\"width:100%; height:auto; vertical-align:middle;\">\n",
" </div>\n",
"</div>"
2015-12-14 02:09:33 +02:00
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Found 4 element sets: LSYM23, LSYM12, LSYM13, LOWER_BLOCK_TO_UPPER_BLOCK\n",
"INFO: Found 4 element sets: USYM23, UPPER_BLOCK_TO_LOWER_BLOCK, USYM13, LOAD\n"
]
},
{
"data": {
"text/plain": [
"Dict{ASCIIString,Any} with 2 entries:\n",
" \"nodes\" => Dict(2=>[0.0,0.0,1.0],11=>[0.0,0.6666666666666666,1.0],39=>…\n",
2015-12-14 17:49:36 +02:00
" \"connectivity\" => Dict(68=>(:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[26,4,14,42]),2=…"
2015-12-14 02:09:33 +02:00
]
},
"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, Quad4, Hex8, LinearSolver\n",
"mesh_lower = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/LOWER_BLOCK.med\")\n",
"mesh_upper = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/UPPER_BLOCK.med\")"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1x75 Array{Int64,2}:\n",
" 1 2 3 4 5 6 7 8 9 10 11 12 … 67 68 69 70 71 72 73 74 75"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sort(collect(keys(mesh_lower[\"nodes\"])))'"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1x48 Array{Int64,2}:\n",
" 1 2 3 4 5 6 7 8 9 10 11 12 … 40 41 42 43 44 45 46 47 48"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sort(collect(keys(mesh_upper[\"nodes\"])))'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Need to renumber nodes."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Dict{Int64,Int64} with 48 entries:\n",
" 2 => 77\n",
" 11 => 86\n",
" 39 => 114\n",
" 46 => 121\n",
" 25 => 100\n",
" 42 => 117\n",
" 29 => 104\n",
" 8 => 83\n",
" 20 => 95\n",
" 14 => 89\n",
" 31 => 106\n",
" 33 => 108\n",
" 18 => 93\n",
" 26 => 101\n",
" 35 => 110\n",
" 17 => 92\n",
" 44 => 119\n",
" 4 => 79\n",
" 37 => 112\n",
" 45 => 120\n",
" 13 => 88\n",
" 30 => 105\n",
" 1 => 76\n",
" 47 => 122\n",
" 32 => 107\n",
" ⋮ => ⋮"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"new_node_numbering = Dict{Int64, Int64}()\n",
"upper_node_ids = sort(collect(keys(mesh_upper[\"nodes\"])))\n",
"for (k, old_node_id) in enumerate(upper_node_ids)\n",
" new_node_numbering[old_node_id] = 75+k\n",
"end\n",
"new_node_numbering"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1x48 Array{Int64,2}:\n",
" 76 77 78 79 80 81 82 83 84 85 … 117 118 119 120 121 122 123"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function renumber!(mesh, node_numbering)\n",
" old_nodes = mesh[\"nodes\"]\n",
" new_nodes = typeof(old_nodes)()\n",
" for (node_id, node_coords) in old_nodes\n",
" new_node_id = node_numbering[node_id]\n",
" new_nodes[new_node_id] = node_coords\n",
" end\n",
" mesh[\"nodes\"] = new_nodes\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" new_elcon = [node_numbering[node_id] for node_id in elcon]\n",
" mesh[\"connectivity\"][elid] = (eltype, elset, new_elcon)\n",
" end\n",
" return mesh\n",
"end\n",
"\n",
"renumber!(mesh_upper, new_node_numbering)\n",
"sort(collect(keys(mesh_upper[\"nodes\"])))'"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Dict{Int64,Tuple{Symbol,Symbol,Array{Int64,1}}} with 92 entries:\n",
2015-12-14 17:49:36 +02:00
" 68 => (:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[101,79,89,117])\n",
" 2 => (:SE2,:OTHER,[77,84])\n",
" 89 => (:HE8,:OTHER,[114,102,103,115,121,118,119,123])\n",
" 11 => (:SE2,:OTHER,[90,80])\n",
" 39 => (:QU4,:USYM13,[93,106,90,80])\n",
" 46 => (:QU4,:LOAD,[96,109,108,95])\n",
" 85 => (:HE8,:OTHER,[94,93,114,116,107,106,121,120])\n",
" 25 => (:SE2,:OTHER,[99,83])\n",
" 55 => (:QU4,:OTHER,[95,113,112,96])\n",
" 42 => (:QU4,:USYM13,[106,92,81,90])\n",
" 66 => (:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[116,88,76,94])\n",
" 58 => (:QU4,:OTHER,[113,101,100,112])\n",
" 29 => (:SE2,:OTHER,[79,101])\n",
" 59 => (:QU4,:OTHER,[87,79,101,113])\n",
" 8 => (:SE2,:OTHER,[88,76])\n",
" 74 => (:QU4,:OTHER,[99,83,103,119])\n",
" 90 => (:HE8,:OTHER,[89,117,101,79,105,122,113,87])\n",
" 57 => (:QU4,:OTHER,[112,100,83,99])\n",
" 20 => (:SE2,:OTHER,[96,95])\n",
" 78 => (:HE8,:OTHER,[104,120,122,105,85,110,108,86])\n",
" 14 => (:SE2,:OTHER,[92,91])\n",
" 31 => (:SE2,:OTHER,[103,102])\n",
" 70 => (:QU4,:OTHER,[98,119,118,97])\n",
" 33 => (:QU4,:USYM23,[88,104,84,76])\n",
" 52 => (:QU4,:LOAD,[111,92,91,110])\n",
2015-12-14 02:09:33 +02:00
" ⋮ => ⋮"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mesh_upper[\"connectivity\"]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: created 59 elements.\n"
]
}
],
"source": [
"# interior elements are of type HE8\n",
"field_problem = LinearElasticityProblem()\n",
"\n",
"for (elid, (eltype, elset, elcon)) in mesh_lower[\"connectivity\"]\n",
" eltype == :HE8 || continue\n",
" element = Hex8(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mesh_lower[\"nodes\"][i] for i in get_connectivity(element)]\n",
" element[\"youngs modulus\"] = 900.0\n",
" element[\"poissons ratio\"] = 0.25\n",
2015-12-14 17:49:36 +02:00
" #info(\"lower: add element with connectivity $elcon\")\n",
2015-12-14 02:09:33 +02:00
" push!(field_problem, element)\n",
"end\n",
"\n",
"for (elid, (eltype, elset, elcon)) in mesh_upper[\"connectivity\"]\n",
" eltype == :HE8 || continue\n",
" element = Hex8(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n",
" element[\"youngs modulus\"] = 900.0\n",
" element[\"poissons ratio\"] = 0.25\n",
2015-12-14 17:49:36 +02:00
" #info(\"upper: add element with connectivity $elcon\")\n",
2015-12-14 02:09:33 +02:00
" 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_upper[\"connectivity\"]\n",
" (eltype == :QU4) && (elset == :LOAD) || continue\n",
" element = Quad4(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n",
2015-12-14 17:49:36 +02:00
" element[\"displacement traction force\"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]\n",
2015-12-14 02:09:33 +02:00
" push!(field_problem, element)\n",
"end\n",
"info(\"created $(length(field_problem.elements)) elements.\")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: created 44 boundary elements.\n"
]
}
],
"source": [
"# boundary conditions\n",
"boundary_problem = DirichletProblem(\"displacement\", 3)\n",
"\n",
"for (elid, (eltype, elset, elcon)) in mesh_lower[\"connectivity\"]\n",
" (eltype == :QU4) && (elset in [:LSYM23, :LSYM12, :LSYM13]) || continue\n",
" element = Quad4(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mesh_lower[\"nodes\"][i] for i in get_connectivity(element)]\n",
" if elset == :LSYM23\n",
" element[\"displacement 1\"] = 0.0\n",
" elseif elset == :LSYM12\n",
" element[\"displacement 3\"] = 0.0\n",
" elseif elset == :LSYM13\n",
" element[\"displacement 2\"] = 0.0\n",
" end\n",
" push!(boundary_problem, element)\n",
"end\n",
"\n",
"for (elid, (eltype, elset, elcon)) in mesh_upper[\"connectivity\"]\n",
" (eltype == :QU4) && (elset in [:USYM23, :USYM12, :USYM13]) || continue\n",
" element = Quad4(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n",
" if elset == :USYM23\n",
" element[\"displacement 1\"] = 0.0\n",
" elseif elset == :USYM12\n",
" element[\"displacement 3\"] = 0.0\n",
" elseif elset == :USYM13\n",
" element[\"displacement 2\"] = 0.0\n",
" end\n",
" push!(boundary_problem, element)\n",
"end\n",
"\n",
"info(\"created $(length(boundary_problem.elements)) boundary elements.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Creating tie contact\n",
"- define slave element surface (the one where integration happend)\n",
"- define potential master elements for slave elements"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM.Core: Element, MortarProblem, calculate_normal_tangential_coordinates!\n",
"\n",
"mortar_side = mesh_upper\n",
"nonmortar_side = mesh_lower\n",
"\n",
"master_elements = JuliaFEM.Core.Element[]\n",
"for (elid, (eltype, elset, elcon)) in mortar_side[\"connectivity\"]\n",
" (eltype == :QU4) && (elset == :UPPER_BLOCK_TO_LOWER_BLOCK) || continue\n",
" element = Quad4(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[mortar_side[\"nodes\"][i] for i in get_connectivity(element)]\n",
" push!(master_elements, element)\n",
"end\n",
"\n",
"contact_problem = MortarProblem(\"displacement\", 3)\n",
"for (elid, (eltype, elset, elcon)) in nonmortar_side[\"connectivity\"]\n",
" (eltype == :QU4) && (elset == :LOWER_BLOCK_TO_UPPER_BLOCK) || continue\n",
" element = Quad4(elcon)\n",
" element[\"geometry\"] = Vector{Float64}[nonmortar_side[\"nodes\"][i] for i in get_connectivity(element)]\n",
" element[\"master elements\"] = master_elements\n",
" calculate_normal_tangential_coordinates!(element, 0.0)\n",
" push!(contact_problem, element)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: # of field problems: 1\n",
"INFO: # of boundary problems: 2\n",
"INFO: Starting iteration 1\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: dim = 369\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1...\n",
"INFO: Assembling boundary 2...\n",
"INFO: dumping matrices to disk, file = matrices_tie_contact_3d_host_1_iteration_1.jld\n",
"INFO: Solving system\n",
2015-12-14 17:49:36 +02:00
"INFO: UMFPACK: solved in 0.18637418746948242 seconds. norm = 0.7276121358181449\n",
"INFO: timing info for iteration:\n"
2015-12-14 02:09:33 +02:00
]
},
{
"data": {
"text/plain": [
2015-12-14 17:49:36 +02:00
"(1,true)"
2015-12-14 02:09:33 +02:00
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
2015-12-14 17:49:36 +02:00
"INFO: boundary assembly : 1.6767330169677734\n",
"INFO: field assembly : 1.1170330047607422\n",
"INFO: dump matrices to disk : 1.0959408283233643\n",
"INFO: solve problem : 0.3340458869934082\n",
"INFO: update element data : 0.022793054580688477\n",
"INFO: non-linear iteration : 4.246567964553833\n",
"INFO: solver finished in 4.391422986984253 seconds.\n"
2015-12-14 02:09:33 +02:00
]
}
],
"source": [
"using JuliaFEM.Core: DirectSolver\n",
"solver = DirectSolver()\n",
"solver.name = \"tie_contact_3d\"\n",
"solver.method = :UMFPACK\n",
2015-12-14 17:49:36 +02:00
"solver.nonlinear_problem = false\n",
2015-12-14 02:09:33 +02:00
"solver.max_iterations = 1\n",
"solver.dump_matrices = true\n",
"push!(solver, field_problem)\n",
"push!(solver, boundary_problem)\n",
"push!(solver, contact_problem)\n",
"call(solver, 0.0)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: nid near corner = 82\n",
2015-12-14 17:49:36 +02:00
"INFO: displacement X = [1.0,1.0,1.0], u = [0.02865803993668297,0.02865803993668297,-0.11229776250350555]\n"
2015-12-14 02:09:33 +02:00
]
},
{
"ename": "LoadError",
"evalue": "LoadError: test failed: isapprox(-(u[i]),known_value)\n in expression: isapprox(-(u[i]),known_value)\nwhile loading In[11], in expression starting on line 14",
"output_type": "error",
"traceback": [
"LoadError: test failed: isapprox(-(u[i]),known_value)\n in expression: isapprox(-(u[i]),known_value)\nwhile loading In[11], in expression starting on line 14",
"",
" in error at ./error.jl:21",
" in default_handler at test.jl:30",
" in do_test at test.jl:53",
" [inlined code] from In[11]:20",
" in anonymous at no file:0"
]
}
],
"source": [
"using JuliaFEM.Test\n",
"\n",
"nid = 0\n",
"for (nid, coords) in mesh_upper[\"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)\n",
"end"
]
},
{
"cell_type": "code",
2015-12-14 17:49:36 +02:00
"execution_count": 28,
2015-12-14 02:09:33 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: XDFM: ndim = 369\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.Hex8}) || continue\n",
" push!(elements, (:Hex8, 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/foobar3.xmf\");"
]
},
{
2015-12-14 17:49:36 +02:00
"cell_type": "markdown",
2015-12-14 02:09:33 +02:00
"metadata": {
"collapsed": true
},
2015-12-14 17:49:36 +02:00
"source": [
"<img src=\"http://results.juliafem.org/2015-12-14-divided-block-results.png\"\n",
" style=\"width:30%; height:auto; vertical-align:middle;\">"
]
2015-12-14 02:09:33 +02:00
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.4.2",
"language": "julia",
"name": "julia-0.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.4.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
