Files

2226 lines
488 KiB
Plaintext
Raw Permalink Normal View History

2015-11-24 03:06:56 +02:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 2d contacts\n",
2015-11-24 03:06:56 +02:00
"\n",
"Author: Jukka Aho\n",
"\n",
"Abstract: 2d tie contact.\n",
"\n",
"## Model 1: three body tie contact\n",
2015-11-24 03:06:56 +02:00
"\n",
"\n",
"![model](http://s4.postimg.org/u9yqeryul/Screenshot_from_2015_11_24_02_00_00.png)\n",
"\n",
"Each element is modelled as own \"body\" and they are connected using tie contacts. Segments 5-6 and 9-10 and 6-7 are slave surfaces, so node 6 or 9 is on at least two tie contacts as slave node. Moreover this model has dirichlet boundary $y=0$ at bottom of body 1 and $x=0$ on left. To get the accurate solution one needs to minimize \n",
"\\begin{equation}\n",
2015-11-24 11:45:51 +02:00
"\\frac{15}{2}u_1^4 + 60 u_1^3 + \\frac{15}{4}u_1^2 u_2^2 + 15 u_1^2 u_2 + 120 u_1^2 + 15 u_1 u_2^2 + 60 u_1 u_2 + \\frac{15}{2}u_2^4 + 60 u_2^3 + 120 u_2^2 + 50 u_2,\n",
2015-11-24 03:06:56 +02:00
"\\end{equation}\n",
"which gives approximate $u_1 = 0.0634862$ and $u_2 = -0.277183$ for the displacement of upper right corner.\n",
"[Wolfram](http://www.wolframalpha.com/input/?i=local+minimum+15*x^4%2F2+%2B+60*x^3+%2B+15*x^2*y^2%2F4+%2B+15*x^2*y+%2B+120*x^2+%2B+15*x*y^2+%2B+60*x*y+%2B+15*y^4%2F2+%2B+60*y^3+%2B+120*y^2+%2B+50*y)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM.Core: Element, Seg2, Quad4, PlaneStressElasticityProblem,\n",
" DirichletProblem, MortarProblem, DirectSolver"
2015-11-24 03:06:56 +02:00
]
},
{
"cell_type": "code",
"execution_count": 10,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nodes = Dict{Int64, Vector{Float64}}(\n",
" 1 => [0.0, 0.0], 2 => [2.0, 0.0],\n",
" 3 => [2.0, 1.0], 4 => [0.0, 1.0],\n",
" 5 => [0.0, 1.0], 6 => [1.0, 1.0],\n",
" 7 => [1.0, 2.0], 8 => [0.0, 2.0],\n",
" 9 => [1.0, 1.0], 10 => [2.0, 1.0],\n",
" 11 => [2.0, 2.0], 12 => [1.0, 2.0]);"
2015-11-24 03:06:56 +02:00
]
},
{
"cell_type": "code",
"execution_count": 11,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"connectivity = Dict{Int64, Vector{Int64}}(\n",
" 1 => [1, 2, 3, 4],\n",
" 2 => [5, 6, 7, 8],\n",
" 3 => [9, 10, 11, 12]);"
]
},
{
"cell_type": "code",
"execution_count": 13,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: number of elements: 3\n"
]
2015-11-24 03:06:56 +02:00
}
],
"source": [
"elements = Element[]\n",
"for c in values(connectivity)\n",
" element = Quad4(c)\n",
" element[\"geometry\"] = Vector{Float64}[nodes[i] for i in c]\n",
" element[\"youngs modulus\"] = 900.0\n",
" element[\"poissons ratio\"] = 0.25\n",
" push!(elements, element)\n",
"end\n",
"nelements = length(elements)\n",
"info(\"number of elements: $nelements\")"
2015-11-24 03:06:56 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create three bodies, each containing one element."
]
},
{
"cell_type": "code",
"execution_count": 14,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"body1 = PlaneStressElasticityProblem(\"body 1\")\n",
"body2 = PlaneStressElasticityProblem(\"body 2\")\n",
"body3 = PlaneStressElasticityProblem(\"body 3\")\n",
2015-11-24 03:06:56 +02:00
"push!(body1, elements[1])\n",
"push!(body2, elements[2])\n",
"push!(body3, elements[3]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Surface traction to the top of bodies 2 and 3:"
]
},
{
"cell_type": "code",
"execution_count": 15,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t2 = Seg2([8, 7])\n",
"t2[\"geometry\"] = Vector{Float64}[nodes[8], nodes[7]]\n",
"t2[\"displacement traction force\"] = Vector{Float64}[[0.0, -100.0], [0.0, -100.0]]\n",
"t3 = Seg2([12, 11])\n",
"t3[\"geometry\"] = Vector{Float64}[nodes[12], nodes[11]]\n",
"t3[\"displacement traction force\"] = Vector{Float64}[[0.0, -100.0], [0.0, -100.0]]\n",
"push!(body2, t2)\n",
"push!(body3, t3);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Boundary conditions: $x=0$ for left boundary."
]
},
{
"cell_type": "code",
"execution_count": 16,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dx1 = Seg2([1, 4])\n",
"dx1[\"geometry\"] = Vector[nodes[1], nodes[4]]\n",
"dx1[\"displacement 1\"] = 0.0\n",
"dx2 = Seg2([5, 8])\n",
"dx2[\"geometry\"] = Vector[nodes[5], nodes[8]]\n",
"dx2[\"displacement 1\"] = 0.0\n",
"bc1 = DirichletProblem(\"displacement\", 2)\n",
"push!(bc1, dx1)\n",
"push!(bc1, dx2);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$y=0$ for bottom of model"
]
},
{
"cell_type": "code",
"execution_count": 17,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dy1 = Seg2([1, 2])\n",
"dy1[\"geometry\"] = Vector[nodes[1], nodes[2]]\n",
"dy1[\"displacement 2\"] = 0.0\n",
"bc2 = DirichletProblem(\"displacement\", 2)\n",
"push!(bc2, dy1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Mortar boundary conditions: tie contact between body 1 and body 2"
]
},
{
"cell_type": "code",
"execution_count": 19,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]\n",
"using JuliaFEM.Core: calculate_normal_tangential_coordinates!\n",
2015-11-24 03:06:56 +02:00
"\n",
"master1 = Seg2([4, 3])\n",
"master1[\"geometry\"] = Vector[nodes[4], nodes[3]]\n",
"slave1 = Seg2([5, 6])\n",
"slave1[\"geometry\"] = Vector[nodes[5], nodes[6]]\n",
"slave1[\"master elements\"] = Element[master1]\n",
"calculate_normal_tangential_coordinates!(slave1, 0.0)\n",
"\n",
"#slave1[\"nodal ntsys\"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]\n",
2015-11-24 03:06:56 +02:00
"contact1 = MortarProblem(\"displacement\", 2)\n",
"push!(contact1, slave1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Tie contact between body 1 and body 3"
]
},
{
"cell_type": "code",
"execution_count": 20,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"slave2 = Seg2([9, 10])\n",
"slave2[\"geometry\"] = Vector[nodes[9], nodes[10]]\n",
"#slave2[\"nodal ntsys\"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]\n",
2015-11-24 03:06:56 +02:00
"slave2[\"master elements\"] = Element[master1]\n",
"calculate_normal_tangential_coordinates!(slave2, 0.0)\n",
2015-11-24 03:06:56 +02:00
"contact2 = MortarProblem(\"displacement\", 2)\n",
"push!(contact2, slave2);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Tie contact between body 2 and body 3"
]
},
{
"cell_type": "code",
"execution_count": 21,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"master2 = Seg2([6, 7])\n",
"master2[\"geometry\"] = Vector[nodes[6], nodes[7]]\n",
"slave3 = Seg2([9, 12])\n",
"slave3[\"geometry\"] = Vector[nodes[9], nodes[12]]\n",
"#slave3[\"nodal ntsys\"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)]\n",
"calculate_normal_tangential_coordinates!(slave3, 0.0)\n",
2015-11-24 03:06:56 +02:00
"slave3[\"master elements\"] = Element[master2]\n",
"contact3 = MortarProblem(\"displacement\", 2)\n",
"push!(contact3, slave3);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All defined. Solve it."
]
},
{
"cell_type": "code",
"execution_count": 22,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"solver = DirectSolver()\n",
"push!(solver, body1)\n",
"push!(solver, body2)\n",
"push!(solver, body3)\n",
"push!(solver, bc1)\n",
"push!(solver, bc2)\n",
"push!(solver, contact1)\n",
"push!(solver, contact2)\n",
"push!(solver, contact3);"
]
},
{
"cell_type": "code",
"execution_count": 23,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Starting solver DirectSolver\n",
2015-11-24 03:06:56 +02:00
"INFO: # of field problems: 3\n",
"INFO: # of boundary problems: 5\n",
"INFO: Starting iteration 1\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1: body 1\n",
"INFO: Assembling body 2: body 2\n",
"INFO: Assembling body 3: body 3\n",
"INFO: dim = 24\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: dirichlet boundary\n",
"INFO: Assembling boundary 3: mortar problem\n",
"INFO: Assembling boundary 4: mortar problem\n",
"INFO: Assembling boundary 5: mortar problem\n",
"INFO: Solving system\n",
"INFO: UMFPACK: solved in 0.3197059631347656 seconds. norm = 0.5357583756107197\n",
"INFO: timing info for iteration:\n",
"INFO: boundary assembly : 0.3747282028198242\n",
"INFO: field assembly : 2.646785020828247\n",
"INFO: dump matrices to disk : 9.5367431640625e-7\n",
"INFO: solve problem : 0.4615659713745117\n",
"INFO: update element data : 0.02040410041809082\n",
"INFO: non-linear iteration : 3.5035040378570557\n",
"INFO: Starting iteration 2\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1: body 1\n",
"INFO: Assembling body 2: body 2\n",
"INFO: Assembling body 3: body 3\n",
"INFO: dim = 24\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: dirichlet boundary\n",
"INFO: Assembling boundary 3: mortar problem\n",
"INFO: Assembling boundary 4: mortar problem\n",
"INFO: Assembling boundary 5: mortar problem\n",
"INFO: Solving system\n",
"INFO: UMFPACK: solved in 0.0003139972686767578 seconds. norm = 0.12311066326855769\n",
"INFO: timing info for iteration:\n",
"INFO: boundary assembly : 0.036119937896728516\n",
"INFO: field assembly : 0.0374150276184082\n",
"INFO: dump matrices to disk : 9.5367431640625e-7\n",
"INFO: solve problem : 0.06756091117858887\n",
"INFO: update element data : 0.00013899803161621094\n",
"INFO: non-linear iteration : 0.1412510871887207\n",
"INFO: Starting iteration 3\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1: body 1\n",
"INFO: Assembling body 2: body 2\n",
"INFO: Assembling body 3: body 3\n",
"INFO: dim = 24\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: dirichlet boundary\n",
"INFO: Assembling boundary 3: mortar problem\n",
"INFO: Assembling boundary 4: mortar problem\n",
"INFO: Assembling boundary 5: mortar problem\n",
"INFO: Solving system\n",
"INFO: UMFPACK: solved in 0.0003230571746826172 seconds. norm = 0.006976385449837204\n",
"INFO: timing info for iteration:\n",
"INFO: boundary assembly : 0.03788185119628906\n",
"INFO: field assembly : 0.036910057067871094\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 0.07133316993713379\n",
"INFO: update element data : 0.0001468658447265625\n",
"INFO: non-linear iteration : 0.1462879180908203\n",
"INFO: Starting iteration 4\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1: body 1\n",
"INFO: Assembling body 2: body 2\n",
"INFO: Assembling body 3: body 3\n",
"INFO: dim = 24\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: dirichlet boundary\n",
"INFO: Assembling boundary 3: mortar problem\n",
"INFO: Assembling boundary 4: mortar problem\n",
"INFO: Assembling boundary 5: mortar problem\n",
"INFO: Solving system\n",
"INFO: UMFPACK: solved in 0.00030493736267089844 seconds. norm = 2.1945519744339283e-5\n",
"INFO: timing info for iteration:\n",
"INFO: boundary assembly : 0.037918806076049805\n",
"INFO: field assembly : 0.036936044692993164\n",
"INFO: dump matrices to disk : 9.5367431640625e-7\n",
"INFO: solve problem : 0.07027888298034668\n",
"INFO: update element data : 0.00013899803161621094\n",
"INFO: non-linear iteration : 0.14528894424438477\n",
"INFO: Starting iteration 5\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1: body 1\n",
"INFO: Assembling body 2: body 2\n",
"INFO: Assembling body 3: body 3\n",
"INFO: dim = 24\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: dirichlet boundary\n",
"INFO: Assembling boundary 3: mortar problem\n",
"INFO: Assembling boundary 4: mortar problem\n",
"INFO: Assembling boundary 5: mortar problem\n",
"INFO: Solving system\n",
"INFO: UMFPACK: solved in 0.0003120899200439453 seconds. norm = 2.172141514107444e-10\n",
"INFO: timing info for iteration:\n"
2015-11-24 03:06:56 +02:00
]
},
{
"data": {
"text/plain": [
"(5,true)"
]
},
"execution_count": 23,
2015-11-24 03:06:56 +02:00
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: boundary assembly : 0.035440921783447266\n",
"INFO: field assembly : 0.0363919734954834\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 0.06740903854370117\n",
"INFO: update element data : 0.0001380443572998047\n",
"INFO: non-linear iteration : 0.13939404487609863\n",
"INFO: solver finished in 4.193101167678833 seconds.\n"
2015-11-24 03:06:56 +02:00
]
}
],
"source": [
"iterations, converged = call(solver, 0.0)"
]
},
{
"cell_type": "code",
"execution_count": 24,
2015-11-24 03:06:56 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: displacement at [2.0,2.0] = [0.06348623177789363,-0.2771830378556528]\n"
2015-11-24 03:06:56 +02:00
]
},
{
"data": {
"text/plain": [
"Test Passed\n",
" Expression: isapprox(u,[0.0634862,-0.277183],atol=1.0e-5)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
2015-11-24 03:06:56 +02:00
}
],
"source": [
"using JuliaFEM.Test\n",
"\n",
"@test converged\n",
"\n",
"X = elements[2](\"geometry\", [1.0, 1.0], 0.0)\n",
"u = elements[2](\"displacement\", [1.0, 1.0], 0.0)\n",
"info(\"displacement at $X = $u\")\n",
"@test isapprox(u, [0.0634862, -0.277183], atol=1.0e-5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Model 2: splitted beam, tie contact\n",
"\n",
"<img src=\"http://results.juliafem.org/splitted-2d-beam/2015-12-22-splitted-beam-mesh.png\">\n",
"\n",
2015-12-30 08:24:18 +02:00
"Put some load on the top, dx=dy=0 on left boundary."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
2015-12-31 12:35:58 +02:00
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM.Preprocess: parse_aster_med_file\n",
"using JuliaFEM.Core: PlaneStressLinearElasticityProblem, DirichletProblem,\n",
" get_connectivity, Quad4, Tri3, Seg2, LinearSolver,\n",
" update!, get_elements"
]
},
{
"cell_type": "code",
2015-12-31 08:40:29 +02:00
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
2015-12-30 08:24:18 +02:00
"create_problems (generic function with 1 method)"
]
},
2015-12-31 08:40:29 +02:00
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
2015-12-30 08:24:18 +02:00
"using JuliaFEM.Core: Element, MortarProblem, calculate_normal_tangential_coordinates!\n",
"\n",
2015-12-30 17:25:48 +02:00
"phi = 0\n",
2015-12-30 08:24:18 +02:00
"rmat = [cos(phi) -sin(phi); sin(phi) cos(phi)]\n",
"\n",
2015-12-30 08:24:18 +02:00
"function create_problems()\n",
"\n",
2015-12-30 08:24:18 +02:00
" mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/2d_beam/BEAM.med\")\n",
" \n",
" for (k, v) in mesh[\"nodes\"]\n",
" mesh[\"nodes\"][k] = rmat*v\n",
" end\n",
" \n",
" field_problem = PlaneStressLinearElasticityProblem()\n",
"\n",
2015-12-30 08:24:18 +02:00
" # field problems\n",
" joo = Dict(:QU4 => Quad4, :TR3 => Tri3)\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype in keys(joo) || continue\n",
" element = joo[eltype](elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
" element[\"youngs modulus\"] = 900.0\n",
" element[\"poissons ratio\"] = 0.25\n",
" push!(field_problem, element)\n",
" end\n",
"\n",
2015-12-30 08:24:18 +02:00
" # neumann boundary condition -1 on y direction\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :SE2 || continue\n",
" elset == :LOAD || continue\n",
" element = Seg2(elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
" # FIXME\n",
" # element[\"displacement traction force\"] = rmat*[0.0, -0.01]\n",
" f = rmat*[0.0, -0.01]\n",
" element[\"displacement traction force\"] = Vector{Float64}[f, f]\n",
" push!(field_problem, element)\n",
" end\n",
"\n",
2015-12-30 08:24:18 +02:00
" # boundary conditions\n",
" boundary_problem = DirichletProblem(\"displacement\", 2)\n",
"\n",
2015-12-30 08:24:18 +02:00
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :SE2 || continue\n",
" elset == :LEFT || continue\n",
" element = Seg2(elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
2015-12-30 17:25:48 +02:00
" # FIXME\n",
" element[\"displacement\"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])\n",
2015-12-30 08:24:18 +02:00
" push!(boundary_problem, element)\n",
" end\n",
"\n",
2015-12-30 08:24:18 +02:00
" info(\"created $(length(get_elements(field_problem))) field elements.\")\n",
" info(\"created $(length(get_elements(boundary_problem))) boundary elements.\")\n",
"\n",
2015-12-30 08:24:18 +02:00
" # Contact definition: contact pair is `LOWER_TO_UPPER <--> UPPER_TO_LOWER`:\n",
"\n",
2015-12-31 08:40:29 +02:00
" slave_surface = :LOWER_TO_UPPER\n",
" master_surface = :UPPER_TO_LOWER\n",
2015-12-30 08:24:18 +02:00
"\n",
2015-12-31 12:35:58 +02:00
" contact_problem = MortarProblem(\"displacement\", 2)\n",
2015-12-30 08:24:18 +02:00
" master_elements = JuliaFEM.Core.Element[]\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :SE2 || continue\n",
2015-12-31 08:40:29 +02:00
" elset == master_surface || continue\n",
2015-12-30 08:24:18 +02:00
" element = Seg2(elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
2015-12-30 17:25:48 +02:00
" element[\"displacement\"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])\n",
2015-12-30 08:24:18 +02:00
" push!(master_elements, element)\n",
2015-12-31 12:35:58 +02:00
" push!(contact_problem, element)\n",
2015-12-30 08:24:18 +02:00
" end\n",
"\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :SE2 || continue\n",
" elset == slave_surface || continue\n",
" element = Seg2(elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
" element[\"master elements\"] = master_elements\n",
" calculate_normal_tangential_coordinates!(element, 0.0)\n",
2015-12-30 17:25:48 +02:00
" element[\"displacement\"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])\n",
2015-12-30 08:24:18 +02:00
" push!(contact_problem, element)\n",
" end\n",
"\n",
" info(\"# of master elements: $(length(master_elements))\")\n",
" info(\"# of slave elements: $(length(contact_problem.elements))\")\n",
"\n",
" return field_problem, boundary_problem, contact_problem\n",
"\n",
"end"
]
},
{
"cell_type": "code",
2015-12-31 08:40:29 +02:00
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
2015-12-30 08:24:18 +02:00
"INFO: Found 5 element sets: LOWER_TO_UPPER, RIGHT, UPPER_TO_LOWER, LEFT, LOAD\n",
"INFO: created 104 field elements.\n",
"INFO: created 4 boundary elements.\n",
2015-12-31 08:40:29 +02:00
"INFO: # of master elements: 14\n",
2015-12-31 12:35:58 +02:00
"INFO: # of slave elements: 34\n",
"INFO: Starting solver divided_beam\n",
"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: plane stress linear elasticity\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 = 210\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: mortar problem\n",
"INFO: Solving system\n",
2015-12-31 12:35:58 +02:00
"INFO: UMFPACK: solved in 0.3899998664855957 seconds. norm = 1.049388621344361\n",
2015-12-30 08:24:18 +02:00
"INFO: timing info for iteration:\n",
2015-12-31 12:35:58 +02:00
"INFO: boundary assembly : 2.3869998455047607\n",
"INFO: field assembly : 1.4670000076293945\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 0.49900007247924805\n",
"INFO: update element data : 0.014999866485595703\n",
"INFO: non-linear iteration : 4.367999792098999\n",
2015-12-30 08:24:18 +02:00
"INFO: Starting iteration 2\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1: plane stress linear elasticity\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 = 210\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: mortar problem\n",
"INFO: Solving system\n",
2015-12-31 12:35:58 +02:00
"INFO: UMFPACK: solved in 0.0 seconds. norm = 1.8501897782230326e-13\n"
]
},
{
"data": {
"text/plain": [
2015-12-30 08:24:18 +02:00
"(2,true)"
]
},
2015-12-31 08:40:29 +02:00
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM.Core: DirectSolver\n",
2015-12-30 08:24:18 +02:00
"\n",
"field_problem, boundary_problem, contact_problem = create_problems()\n",
"\n",
"solver = DirectSolver()\n",
"solver.name = \"divided_beam\"\n",
"solver.method = :UMFPACK\n",
2015-12-30 17:25:48 +02:00
"solver.max_iterations = 2\n",
"solver.dump_matrices = false\n",
"push!(solver, field_problem)\n",
"push!(solver, boundary_problem)\n",
"push!(solver, contact_problem)\n",
"call(solver, 0.0)"
]
},
{
"cell_type": "code",
2015-12-31 08:40:29 +02:00
"execution_count": 4,
"metadata": {
"collapsed": false
},
2015-12-31 12:35:58 +02:00
"outputs": [],
"source": [
2015-12-30 08:24:18 +02:00
"import PyPlot"
]
},
{
"cell_type": "code",
2015-12-31 08:40:29 +02:00
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
2015-12-31 12:35:58 +02:00
"INFO: timing info for iteration:\n",
"INFO: boundary assembly : 0.10899996757507324\n",
"INFO: field assembly : 0.04700016975402832\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 0.06200003623962402\n",
"INFO: update element data : 0.015999794006347656\n",
"INFO: non-linear iteration : 0.23399996757507324\n",
"INFO: solver finished in 4.741999864578247 seconds.\n",
2015-12-31 08:40:29 +02:00
"INFO: displacement at tip: [-0.025032650050967196,-0.1947753325658252]\n"
]
2015-12-30 08:24:18 +02:00
},
{
"data": {
2015-12-31 08:40:29 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMQAAAF0CAYAAADMwuUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XdYVNfWwOEfQxWxIRJbNMYWJSoXY7AQYywx1qjXBGsIXo2KGFQQRQEVUKQoYiXGhjUaiTWWWBNbLJ+xRo0aS+wgUUSkznx/HAdFUREZBob1Ps88wJRz9sCaYc46a69tpNFoNAghhBBCCCGEEEIIUUSo9D0AIYQQQgghhBBCCCHykyTEhBBCCCGEEEIIIUSRIgkxIYQQQgghhBBCCFGkSEJMCCGEEEIIIYQQQhQpkhATQgghhBBCCCGEEEWKJMSEEEIIIYQQQgghRJEiCTEhhBBCCCGEEEIIUaRIQkwIIYQQQgghhBBCFCmSEBNCCCGEEEIIIYQQRYokxIQQQgghhBBCCCFEkaLzhNjhw4dxd3fHzs4OKysrqlatirOzM+fPn89yv6+//hqVSvXcpU6dOroeohBCCCGEEEIIIYQoQkx0vYOQkBAOHDjAF198Qf369bl58yYzZ87EwcGB33//HTs7u8z7mpubM3/+/CyPL1WqlK6HKIQQQgghhBBCCCGKECONRqPR5Q4OHDhAo0aNMDF5knu7cOEC9erVo3v37ixZsgRQKsR++uknEhISdDkcIYQQQgghhBBCCFHE6XzKZJMmTbIkwwBq1KhB3bp1OXv2bJbrNRoNarVakmJCCCGEEEIIIYQQQmf00lRfo9Fw+/ZtbGxsslyflJREyZIlKV26NGXLlsXd3Z2HDx/qY4hCCCGEEEIIIYQQwkDpvIdYdpYtW8aNGzcICgrKvK5ixYqMGjUKBwcH1Go1mzdvZvbs2Rw/fpzdu3djbGysj6EKIYQQQgghhBBCCAOj8x5izzp79iyOjo7Uq1ePPXv2YGRk9ML7BgcHM3bsWFasWIGzs/Nzt8fFxbF161beeecdihUrpsthCyGEEEIIIYQQQogC7NGjR1y+fJm2bds+NyvxWfmaELt16xbNmjUjIyOD33//nfLly7/0/snJyVhZWdGvXz/mzp373O3Lli2jT58+uhquEEIIIYQQQgghhChkli5dSu/evV96n3ybMnn//n3atWtHQkICe/bseWUyDMDCwgJra2vi4+Ozvf2dd94BlCdap06dvByuKKCGDx9ORESEvochRJ6T2BaGSmJbGCKJa2GoJLaFoZLYLjrOnDlDnz59MvNFL5MvCbHk5GQ6derEhQsX2L59O++9916OHvfgwQPi4uIoV65ctrdrp0nWqVMHBweHPBuvKLju378vf2thkCS2haGS2BaGSOJaGCqJbWGoJLaLnpy01dL5KpMZGRk4Oztz8OBBfvzxRxwdHZ+7T0pKCg8ePHju+sDAQAA+++wzXQ9TFBL379/X9xCE0AmJbWGoJLaFIZK4FoZKYlsYKoltkR2dV4h5enqyYcMGOnXqRFxcHEuXLs1ye58+fbh58yb/+c9/6NWrF7Vr1wZg69atbN68mXbt2vH555/repiikKhXr56+hyCETkhsC0MlsS0MkcS1MFQS28JQSWyL7Og8IXb8+HGMjIzYsGEDGzZsyHKbkZERffr0oUyZMnTq1Ilt27YRHR1NRkYGNWvWJDg4GC8vL10PUQghhBBCCCGEEEIUITpPiO3ateuV9ylVqhSLFy/W9VCEAejZs6e+hyCETkhsC0MlsS0MkcS1MFQS28JQSWyL7BhpNBqNvgeRW0ePHqVhw4b83//9nzTIE0IIIYQQQgghhCjCXidPpPOm+kLkpc6dO+t7CELohMS2MFQS28IQSVwLQyWxLQyVxLbIjiTERKHi7u6u7yEIoRMS28JQSWwLQyRxLQyVxLYwVBLbIjsyZVIIIYQQQgghhBBCFHoyZVIIIYQQQgghhBBCiBeQhJgQQgghhBBCCCGEKFIkISYKlbVr1+p7CELohMS2MFQS28IQSVwLQyWxLQyVxLbIjiTERKGyYsUKfQ9BCJ2Q2BaGSmJbGCKJa2GoJLaFoZLYFtmRpvpCCCGEEEIIIYQQotCTpvpCCCGEEEIIIYQQQryAJMSEEEIIIYQQQgghRJEiCTEhhBBCCCGEEEIIUaRIQkwUKq6urvoeghA6IbEtDJXEtjBEEtfCUElsC0MlsS2yIwkxUah8+umn+h6CEDohsS0MlcS2MEQS18JQSWwLQyWxLbIjq0wKIYQQQgghhBBCiEJPVpkUQgghhBBCCCGEEOIFJCEmhBBCCCGEEEIIIYoUSYiJQmXv3r36HoIQOiGxLQyVxLYwRBLXwlBJbAtDJbEtsiMJMVGohIaG6nsIQuiExLYwVBLbwhBJXAtDJbEtDJXEtsiONNUXhUpSUhKWlpb6HoYQeU5iWxgqiW1hiCSuhaGS2BaGSmK76JCm+sJgyZuYMFQS28JQSWwLQyRxLQyVxLYwVBLbIjuSEBNCCCGEEEIIIYQQRYokxIQQQgghhBBCCCFEkSIJMVGojBw5Ut9DEEInJLaFoZLYFoZI4loYKoltYagktkV2JCEmCpUqVaroewhC6ITEtjBUEtvCEElcC0MlsS0MlcS2yI6sMimEEEIIIYQQQgghCj1ZZVIIIYQQQgghhBBCiBeQhJgQQgghhBBCCCGEKFIkISYKlbNnz+p7CELohMS2MFQS28IQSVwLQyWxLQyVxLbIjiTERKHi7e2t7yEIoRMS28JQSWwLQyRxLQyVxLYwVBLbIjuSEBOFysyZM/U9BCF0QmJbGCqJbWGIJK6FoZLYFoZKYltkRxJiolCR5XKFoZLYFoZKYlsYIolrYagktoWhktgW2dFpQuzw4cO4u7tjZ2eHlZUVVatWxdnZmfPnzz933zNnzvDZZ59RokQJypYty1dffUVcXJwuhyeEEEIIIYQQQgghiiATXW48JCSEAwcO8MUXX1C/fn1u3rzJzJkzcXBw4Pfff8fOzg6Aa9eu0bx5c8qUKUNwcDAPHjwgPDyckydPcujQIUxNTXU5TCGEEEIIIYQQQghRhOi0QszT05MrV64wbdo0+vXrx9ixY9mzZw/p6elMnjw5836TJk3i0aNH7Ny5E3d3d3x8fFi1ahXHjx9n0aJFuhyiKGRCQkL0PQQhdEJiWxgqiW1hiCSuhaGS2BaGSmJbZEenCbEmTZpgYpK1CK1GjRrUrVs3y7KnMTExdOzYkcqVK2de16pVK2rVqsWqVat0OURRyCQlJel7CELohMS2MFQS28IQSVwLQyWxLQyVxLbITr431ddoNNy+fRsbGxsArl+/TmxsLB988MFz923UqBF//PFHfg9RFGATJkzQ9xCE0AmJbWGoJLaFIZK4FoZKYlsYKoltkZ18T4gtW7aMGzdu4OzsDMDNmzcBqFChwnP3rVChAvHx8aSlpeXrGIUQQgghhBBCCCGE4dJpU/1nnT17liFDhtC0aVNcXFwAePToEQDm5ubP3d/CwiLzPtJYX8T98w8no6IobmWFysQEVCo0oHw1MgIjI4wvX0ZjbPzkOpVK+dnICIyNSbOxwahYsecel3nftDRUiYnKdcbGYGyMRqVCo1JlbktdvDhof35mO6hUaDIy4PH48tqqVavYtWsXAQEBlCtXLs+3XxicPn2ayMhIRo4cSc2aNfU9HL3IyMhgzJgxlC9fnr59+2ZW3BY1CxYs4K+//sLX1xcrKyt9D0cv/vjjD+bMmUNAQADly5fX93D0IjExkYCAAOrVq0fXrl2LbCzMnDmTK1euMGrUqCL7nrBjxw5iYmLw9/cvsq+HhIQE/Pz8aNCgAZ07dy6SsaBWq1m3bh3Xrl2jR48eRfbz0vr169m4cSNBQUHY2trqezh6kZyczKJFi7C2tqZly5ZF8vUAMG7cOEqXLs2AAQOK7P/I9PR07t+/T+vWrYtsHIjsGWk0Gk1+7OjWrVs0a9aMjIwMfv/998wPKkeOHOHDDz9kyZIl9O7dO8tjvL29CQ8PJyUlJduE2NGjR2nYsCFvvfUWH374YZbbYmNjGTVqFF26dMm87pdffmHmzJmsX78+y32HDBmCg4MD//vf/7Jse/z48SxYsCDLi2bcuHF
"text/plain": [
2015-12-31 12:35:58 +02:00
"PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x000000002D7D0A58>)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
2015-12-30 08:24:18 +02:00
"function plot(field_problem, scaling_factor=30, time=0.0)\n",
"\n",
" fig = PyPlot.figure(figsize=(15, 4))\n",
"\n",
" g_tip = rmat*[100.0, 0.0]\n",
" u_tip = nothing\n",
" \n",
" for element in field_problem.elements\n",
" conn = get_connectivity(element)\n",
" X = element(\"geometry\", time)\n",
" u = element(\"displacement\", time)\n",
" x = X + scaling_factor*u\n",
" if isapprox(element(\"geometry\", [1.0, -1.0], time), g_tip)\n",
" u_tip = element(\"displacement\", [1.0, -1.0], time)\n",
" info(\"displacement at tip: $u_tip\")\n",
" end\n",
"\n",
" # undeformed\n",
" for i=1:length(X)\n",
" px1 = X[i][1]\n",
" py1 = X[i][2]\n",
" px2 = X[mod(i,length(X))+1][1]\n",
" py2 = X[mod(i,length(X))+1][2]\n",
" PyPlot.plot([px1, px2], [py1, py2], \"k-\", alpha=0.5)\n",
" end\n",
"\n",
" # deformed\n",
" for i=1:length(x)\n",
" px1 = x[i][1]\n",
" py1 = x[i][2]\n",
" px2 = x[mod(i,length(x))+1][1]\n",
" py2 = x[mod(i,length(x))+1][2]\n",
" PyPlot.plot([px1, px2], [py1, py2], \"r--\", alpha=0.5)\n",
" end\n",
" end\n",
"\n",
2015-12-30 08:24:18 +02:00
" PyPlot.axis(\"equal\")\n",
2015-12-30 17:25:48 +02:00
" #PyPlot.axis(\"off\")\n",
2015-12-23 17:39:53 +02:00
"\n",
2015-12-30 08:24:18 +02:00
" return u_tip\n",
"end\n",
"\n",
2015-12-30 17:25:48 +02:00
"u_tip = plot(field_problem, 30.0)\n",
"#PyPlot.xlim(95, 105)\n",
"#PyPlot.ylim(-10, 0)\n",
"PyPlot.grid()"
2015-12-23 17:39:53 +02:00
]
},
{
"cell_type": "code",
2015-12-31 08:40:29 +02:00
"execution_count": 6,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
2015-12-30 08:24:18 +02:00
"data": {
"text/plain": [
2015-12-31 08:40:29 +02:00
"true"
2015-12-30 08:24:18 +02:00
]
},
2015-12-31 08:40:29 +02:00
"execution_count": 6,
2015-12-30 08:24:18 +02:00
"metadata": {},
"output_type": "execute_result"
2015-12-23 17:39:53 +02:00
}
],
"source": [
2015-12-31 08:40:29 +02:00
"isapprox(u_tip, [-0.025032650050963334,-0.19477533256579582]) || warn(\"different result\")"
2015-12-23 17:39:53 +02:00
]
},
{
"cell_type": "code",
2015-12-31 08:40:29 +02:00
"execution_count": 7,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
2015-12-30 08:24:18 +02:00
"data": {
"text/plain": [
2015-12-31 08:40:29 +02:00
"true"
2015-12-30 08:24:18 +02:00
]
},
2015-12-31 08:40:29 +02:00
"execution_count": 7,
2015-12-30 08:24:18 +02:00
"metadata": {},
"output_type": "execute_result"
2015-12-23 17:39:53 +02:00
}
],
"source": [
2015-12-31 08:40:29 +02:00
"isapprox(norm(u_tip), 0.19637735038616433)"
2015-12-30 08:24:18 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Model 3: splitted beam, frictionless contact, small sliding\n",
2015-12-23 17:39:53 +02:00
"\n",
2015-12-30 08:24:18 +02:00
"Modification to the above: allow beams to slide in tangential direction without a friction. Global assembly operator can be modified by writing own functions `preprocess_assembly!` and `postprocess_assembly!`. Here we simply remove all kinematic constraints in tangential direction to achieve frictionless sliding between bodies."
2015-12-23 17:39:53 +02:00
]
},
{
"cell_type": "code",
2015-12-31 08:40:29 +02:00
"execution_count": 8,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
2015-12-30 08:24:18 +02:00
"INFO: Found 5 element sets: LOWER_TO_UPPER, RIGHT, UPPER_TO_LOWER, LEFT, LOAD\n",
"INFO: created 104 field elements.\n",
"INFO: created 4 boundary elements.\n",
2015-12-31 08:40:29 +02:00
"INFO: # of master elements: 14\n",
2015-12-31 12:35:58 +02:00
"INFO: # of slave elements: 34\n",
2015-12-23 17:39:53 +02:00
"INFO: Starting solver divided_beam_small_sliding_contact\n",
"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: plane stress linear elasticity\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 = 210\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
2015-12-30 08:24:18 +02:00
"INFO: Assembling boundary 2: mortar problem\n",
"INFO: assemble: doing postprocess for problem JuliaFEM.Core.BoundaryProblem{JuliaFEM.Core.MortarProblem} assembly\n",
"INFO: postprocess mortar assembly: remove contraints in tangent direction on boundary.\n",
"INFO: postprocess mortar assembly: done.\n",
2015-12-23 17:39:53 +02:00
"INFO: Solving system\n",
2015-12-31 12:35:58 +02:00
"INFO: UMFPACK: solved in 0.0 seconds. norm = 3.146559185661418\n",
2015-12-30 08:24:18 +02:00
"INFO: timing info for iteration:\n",
2015-12-31 12:35:58 +02:00
"INFO: boundary assembly : 0.7179999351501465\n",
"INFO: field assembly : 0.07800006866455078\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 0.09299993515014648\n",
"INFO: update element data : 0.0\n",
"INFO: non-linear iteration : 0.8889999389648438\n",
2015-12-30 08:24:18 +02:00
"INFO: Starting iteration 2\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1: plane stress linear elasticity\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 = 210\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: mortar problem\n",
"INFO: assemble: doing postprocess for problem JuliaFEM.Core.BoundaryProblem{JuliaFEM.Core.MortarProblem} assembly\n",
"INFO: postprocess mortar assembly: remove contraints in tangent direction on boundary.\n",
"INFO: postprocess mortar assembly: done.\n",
"INFO: Solving system\n",
2015-12-31 12:35:58 +02:00
"INFO: UMFPACK: solved in 0.0 seconds. norm = 3.5785351812097138e-12\n",
2015-12-30 08:24:18 +02:00
"INFO: timing info for iteration:\n",
2015-12-31 12:35:58 +02:00
"INFO: boundary assembly : 0.10900020599365234\n",
"INFO: field assembly : 0.06299996376037598\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 0.09299993515014648\n",
"INFO: update element data : 0.0\n",
"INFO: non-linear iteration : 0.2650001049041748\n",
"INFO: solver finished in 1.1540000438690186 seconds.\n",
2015-12-31 08:40:29 +02:00
"INFO: displacement at tip: [-0.03914822922685343,-0.592340359574268]\n"
2015-12-23 17:39:53 +02:00
]
},
{
"data": {
2015-12-31 08:40:29 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMQAAAF0CAYAAADMwuUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcVHX3wPEPw4gouCGaIonmguIa7hthZi5pWpJ74ZqiKAgKKosrogEiKW5puKA+uZYYRm6luaXhkqam/lxCcSVARECY+f1xc0tRlhkG8rxfL15Pz8y9d84wd3DumfM9x0ir1WoRQgghhBBCCCGEEOI1oTJ0AEIIIYQQQgghhBBCFCRJiAkhhBBCCCGEEEKI14okxIQQQgghhBBCCCHEa0USYkIIIYQQQgghhBDitSIJMSGEEEIIIYQQQgjxWpGEmBBCCCGEEEIIIYR4rUhCTAghhBBCCCGEEEK8ViQhJoQQQgghhBBCCCFeK5IQE0IIIYQQQgghhBCvFUmICSGEEEIIIYQQQojXit4TYkeOHMHV1ZV69ephbm6OjY0Nffr04fz5889sN2jQIFQq1XM/devW1XeIQgghhBBCCCGEEOI1otb3A8yZM4eDBw/yySef0LBhQ+Lj41mwYAH29vYcOnSIevXqPd62ePHiLF++/Jn9y5Qpo+8QhRBCCCGEEEIIIcRrxEir1Wr1+QAHDx6kWbNmqNVPcm8XLlygQYMGODk5sXr1akCpENu8eTPJycn6DEcIIYQQQgghhBBCvOb0vmSyVatWzyTDAGrWrImdnR1nz5595natVotGo5GkmBBCCCGEEEIIIYTQG4M01ddqtdy8eRNLS8tnbk9NTaV06dKULVuW8uXL4+rqyv379w0RohBCCCGEEEIIIYT4j9J7D7EXWbNmDdevX2fmzJmPb7OyssLb2xt7e3s0Gg3bt29n4cKFnDhxgp9++gljY2NDhCqEEEIIIYQQQggh/mP03kPs386ePUuLFi1o0KAB+/btw8jIKNttAwMD8fHxYd26dfTp0+e5++/cuUNMTAzVqlWjRIkS+gxbCCGEEEIIIYQQQhRiDx484PLly3Tq1Om5VYn/VqAJsRs3btCmTRuysrI4dOgQlSpVeun2aWlpmJubM2TIEJYuXfrc/WvWrGHgwIH6ClcIIYQQQgghhBBCFDGRkZEMGDDgpdsU2JLJpKQkunTpQnJyMvv27XtlMgzA1NQUCwsLEhISXnh/tWrVAOWJ1q1bV5fhvpbGjRtHaGioocMQRYycNyKv5NwReSHnjcgLOW9EXsh5I/JKzh2RF3Le6MaZM2cYOHDg43zRyxRIQiwtLY3u3btz4cIFdu7cSZ06dXK0371797hz5w4VKlR44f2PlknWrVsXe3t7ncX7uipTpoz8HkWuyXkj8krOHZEXct6IvJDzRuSFnDcir+TcEXkh541u5aStlt6nTGZlZdGnTx8OHz7Mhg0baNGixXPbpKenc+/evedunzFjBgCdO3fWd5hCCCGEEEIIIYQQ4jWh9woxT09PoqKi6N69O3fu3CEyMvKZ+wcOHEh8fDxvv/02/fv3x9bWFoCYmBi2b99Oly5d6NGjh77DFEIIIYQQQgghhBCvCb0nxE6cOIGRkRFRUVFERUU9c5+RkREDBw6kXLlydO/enR07drBy5UqysrKoVasWgYGBjB8/Xt8hCiGEEEIIIYQQQojXiN4TYnv27HnlNmXKlGHVqlX6DkW8Qr9+/QwdgiiC5LwReSXnjsgLOW9EXsh5I/JCzhuRV3LuiLyQ86bgGWm1Wq2hg8ir2NhYmjRpwm+//SbN54QQQgghhBBCCCFeY7nJE+m9qb4QQgghhBBCCCGEEIWJJMSEEEIIIYQQQgghxGtFEmJCCCGEEEIIIYQQ4rUiCTEhhBBCCCGEEEII8VqRhJgQQgghhBBCCCGEeK1IQkwIIYQQQgghhBBCvFYkISaEEEIIIYQQQgghXiuSEBNCCCGEEEIIIYQQrxVJiAkhhBBCCCGEEEKI14okxIQQQgghhBBCCCHEa0USYkIIIYQQQgghhBDitSIJMSGEEEIIIYQQQgjxWpGEmBBCCCGEEEIIIYR4rUhCTAghhBBCCCGEEEK8ViQhJoQQQgghhBBCCCFeK5IQE0IIIYQQQgghhBCvFUmICSGEEEIIIYQQQojXiiTEhBBCCCGEEEIIIcRrRRJiQgghhBBCCCGEEOK1IgkxIYQQQgghhBBCCPFakYSYEEIIIYQQQgghhHit6DUhduTIEVxdXalXrx7m5ubY2NjQp08fzp8//9y2Z86coXPnzpQqVYry5cvz2WefcefOHX2GJ4QQQgghhBBCCCFeQ2p9HnzOnDkcPHiQTz75hIYNGxIfH8+CBQuwt7fn0KFD1KtXD4C4uDgcHBwoV64cgYGB3Lt3j+DgYH7//Xd+/fVXihUrps8whRBCCCGEEEIIIcRrRK8JMU9PT5o1a4Za/eRh+vTpQ4MGDZg9ezarV68GYNasWTx48IBjx45hbW0NQPPmzenYsSMrVqxg+PDh+gxTCCGEEEIIIYQQQrxG9LpkslWrVs8kwwBq1qyJnZ0dZ8+efXzbpk2b6Nat2+NkGECHDh2oXbs269ev12eIQgghhBBCCCGEEOI1U+BN9bVaLTdv3sTS0hKAa9eucfv2bZo2bfrcts2aNePYsWMFHaIQQgghhBBCCCGE+A8r8ITYmjVruH79On369AEgPj4egMqVKz+3beXKlUlISODhw4cFGqMQQgghhBBCCCGE+O/Saw+xfzt79iyjR4+mdevWODs7A/DgwQMAihcv/tz2pqamj7eRxvpCFD537tzh/8aOpdL582RYWZHy9ttoVSowNkarUqFVqXhobs692rXh4UMsdu1Ca2ys/KjVyo+xMSb37qFVqUirWpWsMmWU+//ZH2NjskxM0BYrBpmZoNWCsTGodJvPX7JkCXv27GHu3LlYWVnp9NgFbdq0ady5c4fAwEDMzc0NHU6eZWRk4O7uTqlSpZgwYcLjyuKi6Pbt23h5edG2bVv69OlTpF+Xn376icWLFzNt2jRsbW0NHU6+fPXVVxw4cICQkBAsLCwMHU6eaTQafH19uXv3LgEBAUX6vZKZmcnYsWOpUqUKbm5uRfq9cvnyZfz9/WnSpAkDBgwo0q/Lr7/+ytq1axk/fvwzLU6KogULFnDw4EGCg4Nf+IV8UbJy5Ur++usvRowYQYUKFQwdTp5pNBpcXFywsbFh7NixRfp9D8rfsaSkJN57770i/b4X4nVgpNVqtQXxQDdu3KBNmzZkZWVx6NAhKlWqBMDRo0dp3rw5q1evZsCAAc/s4+XlRXBwMOnp6S9MiMXGxtKkSRMcHBwoU6bMM/f169ePfv366e8JCSFYHxHBgc8/p6VKRYJazaUSJTDWalEDKq0WY+COsTFRZmaYaDSMuXULY54vTW2YmYkJcLp4cf5+wXt9b4kSxBYvTuOUFDqnpACQBWQCGiMjTLVaams0ZAD7S5VC88/9GiMjsv757x9LluSBSkWr5GTKZWUp2xgZkanVcj89nYfp6RwBqhkZ0bRUKe6rVPxZogSpKhXpKhXpRkakAxgZ6eNXqTMajYZr166h0WgoX758kf1QqdVq+fvvv0lJSUGr1VKyZMki/WH//v373LlzB4AKFSpQsmRJA0eUN1qtlhs3bpCRkVHkX5P09HRu375NVlYWFhYWlCpVytAh5YlWqyUpKYnk5GS0Wi0lSpSgYsWKhg4rz9LT07lx4wYAFStWpESJEgaOKG80Gg137tx5/MVv2bJln/usWpT8/fffJCcnY2pqyhtvvGHocPIsMzOT69ev/yf+XdFqtdy+fZsHDx5gampKxYoVMSrkn1Gy8/DhQ65fvw5AiRIlsLS0RKXjLz4LUmZmJikpKSxevJi+ffsaOhwh/tPWrVvHunXrnrktKSmJvXv38ttvv2Fvb//S/QukQiwpKYkuXbqQnJzMvn37HifD4MlSyUdLJ58WHx9P+fLlX1kdFhoa+sonKoTQvcY1a3K7bl3ULi7UqlmTWtls9+G/b8jMxEijwSgjA6PMTDKTk9E8fEi1UqWoVqwYRhoNZGUp22g0NCpVisxSpTC5do2Sf/6p3JeVhVFmJkaZmagfPKDkzZug1fJ
2015-12-23 17:39:53 +02:00
"text/plain": [
2015-12-31 12:35:58 +02:00
"PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x00000000303F3F98>)"
2015-12-23 17:39:53 +02:00
]
},
"metadata": {},
2015-12-30 08:24:18 +02:00
"output_type": "display_data"
2015-12-23 17:39:53 +02:00
},
{
2015-12-31 08:40:29 +02:00
"data": {
"text/plain": [
"true"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
2015-12-23 17:39:53 +02:00
}
],
"source": [
2015-12-30 08:24:18 +02:00
"using JuliaFEM.Core: SparseMatrixCOO, Element, get_integration_points, get_jacobian, get_connectivity, add!\n",
"using JuliaFEM.Core: BoundaryAssembly, BoundaryProblem, get_elements\n",
"\n",
"function calculate_normal_tangential_coordinates(elements::Vector{Element}, time::Real)\n",
" P = SparseMatrixCOO()\n",
" field_dim = 2\n",
" for element in elements\n",
2015-12-31 12:35:58 +02:00
" haskey(element, \"normal-tangential coordinates\") || continue\n",
2015-12-30 08:24:18 +02:00
" for ip in get_integration_points(element, Val{2})\n",
" J = get_jacobian(element, ip, time)\n",
" w = ip.weight*norm(J)\n",
" nt = transpose(element(\"normal-tangential coordinates\", ip, time))\n",
" normal = nt[1,:]\n",
" tangent = nt[2,:]\n",
" for nid in get_connectivity(element)\n",
" ndofs = [2*(nid-1)+1, 2*(nid-1)+2]\n",
" add!(P, [2*(nid-1)+1], ndofs, normal)\n",
" add!(P, [2*(nid-1)+2], ndofs, tangent)\n",
" end\n",
" end\n",
" end\n",
" P = sparse(P)\n",
" for i=1:size(P,1)\n",
2015-12-30 17:25:48 +02:00
" n = norm(P[i,:])\n",
" if n > 0.0\n",
" P[i,:] = P[i,:] / n\n",
2015-12-30 08:24:18 +02:00
" end\n",
" end\n",
" return P\n",
"end\n",
"\n",
"import JuliaFEM.Core: postprocess_assembly!\n",
"\n",
"function JuliaFEM.Core.postprocess_assembly!(\n",
" assembly::BoundaryAssembly, problem::BoundaryProblem{MortarProblem},\n",
" time::Real)\n",
" info(\"postprocess mortar assembly: remove contraints in tangent direction on boundary.\")\n",
" C1 = sparse(assembly.C1)\n",
" C2 = sparse(assembly.C2)\n",
" dim = size(C1, 1)\n",
" P = calculate_normal_tangential_coordinates(get_elements(problem), time)\n",
" C1 = P*C1\n",
" C2 = P*C2\n",
" for i=2:2:dim\n",
" C1[i,:] = 0\n",
" C2[i,:] = 0\n",
" end\n",
" assembly.C1 = C1\n",
" assembly.C2 = C2\n",
" info(\"postprocess mortar assembly: done.\")\n",
"end\n",
"\n",
"\n",
"field_problem, boundary_problem, contact_problem = create_problems()\n",
2015-12-23 17:39:53 +02:00
"using JuliaFEM.Core: DirectSolver\n",
"solver = DirectSolver()\n",
"solver.name = \"divided_beam_small_sliding_contact\"\n",
"solver.method = :UMFPACK\n",
2015-12-30 17:25:48 +02:00
"solver.max_iterations = 10\n",
2015-12-23 17:39:53 +02:00
"solver.dump_matrices = false\n",
"push!(solver, field_problem)\n",
"push!(solver, boundary_problem)\n",
"push!(solver, contact_problem)\n",
2015-12-30 08:24:18 +02:00
"call(solver, 0.0)\n",
2015-12-23 17:39:53 +02:00
"\n",
2015-12-30 08:24:18 +02:00
"u_tip = plot(field_problem, 30)\n",
2015-12-31 08:40:29 +02:00
"isapprox(u_tip, [-0.03914822922685343,-0.592340359574268]) || warn(\"result changed\")"
2015-12-30 17:25:48 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Frictionless 2d contact with active and inactive nodes"
2015-12-23 17:39:53 +02:00
]
},
{
"cell_type": "code",
2015-12-31 12:35:58 +02:00
"execution_count": 76,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
2015-12-30 17:25:48 +02:00
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Found 5 element sets: LOWER_TO_UPPER, RIGHT, UPPER_TO_LOWER, LEFT, LOAD\n",
"INFO: created 104 field elements.\n",
"INFO: created 4 boundary elements.\n",
2015-12-31 08:40:29 +02:00
"INFO: # of master elements: 14\n",
2015-12-31 12:35:58 +02:00
"INFO: # of slave elements: 34\n",
2015-12-30 17:25:48 +02:00
"INFO: Starting solver divided_beam_inequality_constraints\n",
"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: plane stress linear elasticity\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 = 210\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: mortar problem\n",
2015-12-31 08:40:29 +02:00
"INFO: assemble: doing postprocess for problem JuliaFEM.Core.BoundaryProblem{JuliaFEM.Core.MortarProblem} assembly\n",
2015-12-30 17:25:48 +02:00
"INFO: postprocess mortar assembly: peforming PDASS.\n",
2015-12-31 12:35:58 +02:00
"INFO: weighted gap = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.0,2.5,0.0,0.0,-0.0,2.5,0.0,0.0,0.0,5.0,0.0,0.0,-0.0,5.0,0.0,0.0,-0.0,5.0,0.0,0.0,-0.0,5.0,0.0,0.0,-0.0,5.0,0.0,0.0,-0.0,5.0,0.0,0.0,0.0,5.0,0.0,0.0,0.0,5.0,0.0,5.0,0.0,5.0,-0.0,5.0,-0.0,5.0,-0.0,5.0,0.0,5.0,0.0,5.0,0.0,5.0,0.0,5.0,0.0,5.0,-0.0,5.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]\n",
"INFO: lambda = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]\n",
"INFO: dof 45: normal gap = 2.5, X=[100.0,10.0], C=0.0, la=0.0\n",
"INFO: contact dof 45 active (C = 0.0)\n",
"INFO: dof 49: normal gap = 2.5, X=[0.0,10.0], C=-1.0, la=0.0\n",
"INFO: dof 53: normal gap = 5.0, X=[5.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 57: normal gap = 5.0, X=[10.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 61: normal gap = 5.0, X=[15.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 65: normal gap = 5.0, X=[20.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 69: normal gap = 5.0, X=[25.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 73: normal gap = 5.0, X=[30.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 77: normal gap = 5.0, X=[35.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 81: normal gap = 5.0, X=[40.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 83: normal gap = 5.0, X=[95.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 85: normal gap = 5.0, X=[90.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 87: normal gap = 5.0, X=[85.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 89: normal gap = 5.0, X=[80.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 91: normal gap = 5.0, X=[75.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 93: normal gap = 5.0, X=[70.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 95: normal gap = 5.0, X=[65.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 97: normal gap = 5.0, X=[60.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 99: normal gap = 5.0, X=[55.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 101: normal gap = 5.0, X=[50.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 103: normal gap = 5.0, X=[45.0,10.0], C=0.0, la=0.0\n",
2015-12-30 17:25:48 +02:00
"INFO: postprocess mortar assembly: done.\n",
"INFO: Solving system\n",
2015-12-31 12:35:58 +02:00
"INFO: UMFPACK: solved in 0.0 seconds. norm = 5.426685089269588\n",
2015-12-31 08:40:29 +02:00
"INFO: timing info for iteration:\n",
2015-12-31 12:35:58 +02:00
"INFO: boundary assembly : 0.15599989891052246\n",
"INFO: field assembly : 0.06300020217895508\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 0.07800006866455078\n",
"INFO: update element data : 0.0\n",
"INFO: non-linear iteration : 0.2970001697540283\n",
2015-12-30 17:25:48 +02:00
"INFO: Starting iteration 2\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1: plane stress linear elasticity\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 = 210\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assembling boundary 1: dirichlet boundary\n",
"INFO: Assembling boundary 2: mortar problem\n",
2015-12-31 08:40:29 +02:00
"INFO: assemble: doing postprocess for problem JuliaFEM.Core.BoundaryProblem{JuliaFEM.Core.MortarProblem} assembly\n",
2015-12-30 17:25:48 +02:00
"INFO: postprocess mortar assembly: peforming PDASS.\n",
2015-12-31 12:35:58 +02:00
"INFO: weighted gap = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.33,0.0,0.0,0.0,-0.01,2.49,0.0,0.0,-0.07,4.96,0.0,0.0,-0.15,4.89,0.0,0.0,-0.23,4.77,0.0,0.0,-0.3,4.61,0.0,0.0,-0.37,4.41,0.0,0.0,-0.42,4.18,0.0,0.0,-0.47,3.92,0.0,0.0,-0.51,3.65,-0.67,0.2,-0.67,0.51,-0.67,0.82,-0.66,1.14,-0.66,1.45,-0.65,1.77,-0.64,2.09,-0.62,2.42,-0.6,2.73,-0.58,3.05,-0.55,3.35,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]\n",
"INFO: lambda = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,-0.0,0.0,0.0,0.74,-0.19,0.0,0.0,-0.0,0.0,0.0,0.0,0.0,-0.0,0.0,0.0,-0.0,0.0,0.0,0.0,0.0,-0.0,0.0,0.0,-0.0,0.0,0.0,0.0,0.0,-0.0,0.0,0.0,-0.0,0.0,0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-0.0,-0.0,0.0,0.0,-4.87,-4.26]\n",
"INFO: dof 45: normal gap = 0.0, X=[100.0,10.0], C=0.0, la=-0.0\n",
"INFO: contact dof 45 active (C = 0.0)\n",
"INFO: dof 49: normal gap = 2.493, X=[0.0,10.0], C=-1.0, la=-0.18516\n",
"INFO: dof 53: normal gap = 4.959, X=[5.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 57: normal gap = 4.885, X=[10.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 61: normal gap = 4.768, X=[15.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 65: normal gap = 4.606, X=[20.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 69: normal gap = 4.407, X=[25.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 73: normal gap = 4.179, X=[30.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 77: normal gap = 3.924, X=[35.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 81: normal gap = 3.646, X=[40.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 83: normal gap = 0.205, X=[95.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 85: normal gap = 0.514, X=[90.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 87: normal gap = 0.825, X=[85.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 89: normal gap = 1.138, X=[80.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 91: normal gap = 1.455, X=[75.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 93: normal gap = 1.774, X=[70.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 95: normal gap = 2.094, X=[65.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 97: normal gap = 2.415, X=[60.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 99: normal gap = 2.734, X=[55.0,10.0], C=0.0, la=0.0\n",
"INFO: dof 101: normal gap = 3.048, X=[50.0,10.0], C=0.0, la=-0.0\n",
"INFO: dof 103: normal gap = 3.353, X=[45.0,10.0], C=0.0, la=0.0\n",
2015-12-30 17:25:48 +02:00
"INFO: postprocess mortar assembly: done.\n",
"INFO: Solving system\n",
2015-12-31 12:35:58 +02:00
"INFO: UMFPACK: solved in 0.0 seconds. norm = 1.303211404108607e-12\n",
2015-12-31 08:40:29 +02:00
"INFO: timing info for iteration:\n",
2015-12-31 12:35:58 +02:00
"INFO: boundary assembly : 0.1399998664855957\n",
"INFO: field assembly : 0.06200003623962402\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 0.07800006866455078\n",
"INFO: update element data : 0.0\n",
"INFO: non-linear iteration : 0.2799999713897705\n",
"INFO: solver finished in 0.5770001411437988 seconds.\n",
"INFO: displacement at tip: [-0.036167518632111914,-0.4890644234713899]\n"
2015-12-30 17:25:48 +02:00
]
},
2015-12-23 17:39:53 +02:00
{
"data": {
2015-12-31 12:35:58 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABLkAAAF0CAYAAADPdesAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XdcleX7wPHPORymCCiCOHFmuDX3+jlSM7OpqaWZWZmoOcCNG0GU4d7m/poDzdDK1DRHrpxZWmqJIeKARJB5OOf3xw0CykEFFNDr/XqdF8p5xn0OZzzP9VzXdWuMRqMRIYQQQgghhBBCCCEKMW1+D0AIIYQQQgghhBBCiNySIJcQQgghhBBCCCGEKPQkyCWEEEIIIYQQQgghCj0JcgkhhBBCCCGEEEKIQk+CXEIIIYQQQgghhBCi0JMglxBCCCGEEEIIIYQo9CTIJYQQQgghhBBCCCEKPQlyCSGEEEIIIYQQQohCT4JcQgghhBBCCCGEEKLQkyCXEEIIIYQQQgghhCj0chTkOn78OIMGDaJGjRrY2tri6upK9+7duXjxYqblPv74Y7Ra7UM3Nze3PBm8EEIIIYQQQgghhBAAupys5Ofnx+HDh+nWrRu1a9fm+vXrzJs3j/r163PkyBFq1Khxf1lLS0uWL1+eaX17e/vcjVoIIYQQQgghhBBCiAw0RqPR+KQrHT58mIYNG6LTpcfILl26RK1atejatStr1qwBVCbXli1buHv3bt6NWAghhBBCCCGEEEKIB+SoXLFp06aZAlwAVapUoXr16ly4cCHT741GIwaDQQJdQgghhBBCCCGEEOKpybPG80ajkRs3blCiRIlMv4+Li8POzg4HBwccHR0ZNGgQ9+7dy6vdCiGEEEIIIYQQQgiRs55cWVm3bh3h4eF4e3vf/13p0qUZNWoU9evXx2Aw8P3337NgwQLOnDnDvn37MDMzy6vdCyGEEEIIIYQQQogXWI56cj3owoULNG7cmFq1anHgwAE0Go3JZX19fRk3bhzr16+ne/fuD91/+/Ztdu7cSYUKFbC2ts7t0IQQQgghhBBCCCFEIRUfH8+VK1fo2LHjQ9WDD8p1kCsiIoLmzZuTkpLCkSNHcHFxyXb5hIQEbG1t+eSTT1iyZMlD969bt45evXrlZkhCCCGEEEIIIYQQ4jmydu1aPvzww2yXyVW5YnR0NJ06deLu3bscOHDgkQEuACsrK4oXL05UVFSW91eoUAFQg3dzc8vN8ITI1rBhwwgKCsrvYQjxXJH3lRB5S95TQuQteU8JkffkfSWetvPnz9OrV6/78aLs5DjIlZCQQJcuXbh06RK7d+/m5Zdffqz1YmJiuH37Nk5OTlnen1ai6ObmRv369XM6PCEeyd7eXl5jQuQxeV8JkbfkPSVE3pL3lBB5T95X4ll5nJZWOZpdMSUlhe7du3P06FE2bdpE48aNH1omMTGRmJiYh34/depUAF577bWc7FoIIYQQQgghhBBCiIfkKJPLw8ODkJAQunTpwu3bt1m7dm2m+3v16sX169epV68eH3zwAdWqVQNg586dfP/993Tq1Im33nor96MXQgghhBBCCCGEEIIcBrnOnDmDRqMhJCSEkJCQTPdpNBp69epFsWLF6NKlC7t27WLVqlWkpKRQtWpVfH198fT0zJPBCyGEEEIIIYQQQggBOQxy7d2795HL2Nvbs3r16pxsXohnomfPnvk9BCGeO/K+EiJvyXtKiLwl7ykh8p68r0RBojEajcb8HkRGJ0+e5JVXXuHEiRPSvE4IIYQQQgghhBDiBfYkcaIcNZ4XQgghhBBCCCGEEKIgkSCXEEIIIYQQQgghhCj0JMglhBBCCCGEEEIIIQo9CXIJIYQQQgghhBBCiEJPglxCCCGEEEIIIYQQotCTIJcQQgghhBBCCCGEKPQkyCWEEEIIIYQQQgghCj0JcgkhhBBCCCGEEEKIQk+CXEIIIYQQQgghhBCi0JMglxBCCCGEEEIIIYQo9CTIJYQQQgghhBBCCCEKPQlyCSGEEEIIIYQQQohCT4JcQgghhBBCCCGEEKLQkyCXEEIIIYQQQgghhCj0JMglhBBCCCGEEEIIIQo9CXIJIYQQQgghhBBCiEJPglxCCCGEEEIIIYQQotCTIJcQQgghhBBCCCGEKPQkyCWEEEIIIYQQQgghCj0JcgkhhBBCCCGEEEKIQk+CXEIIIYQQQgghhBCi0JMglxBCCCGEEEIIIYQo9CTIJYQQQgghhBBCCCEKPQlyCSGEEEIIIYQQQohCT4JcQgghhBBCCCGEEKLQkyCXEEIIIYQQQgghhCj0JMglhBBCCCGEEEIIIQo9CXIJIYQQQgghhBBCiEJPglxCCCGEEEIIIYQQotCTIJcQQgghhBBCCCGEKPQkyCWEEEIIIYQQQgghCj0JcgkhhBBCCCGEEEKIQk+CXEIIIYQQQgghhBCi0JMglxBCCCGEEEIIIYQo9HT5PQBTzq1bh37/fpLt7IgrV+6h+4v/8AOkpIBWS5Hr1wEwAmg0nD5zhu/Cw+nz6aeUTr0PrRajRqP+rdGQYmVFYrFiRDdrhtHKKssxFLl6FcurV9H995/avkYDGg1GrYoNJtvZkeToSPzLL2e5vjYxEYs7d7D6+281tvt3aO9vI8nenmRnZ1IcHLLchll8PNq4OHSRkWBmBhoNhrTHodViNDfHYGmJ3t4edFn/OTXJyWiSksCYOgozM7WN1DEcPHyYr1auZOKkSbi6uma5DZEuMjKS8ePHU79+fd5++21KlCiR30Mq0PR6PaNHj6ZSpUp89NFH2Nra5veQCrxp06aRnJyMp6enPF+PYdWqVZw4cYJp06ZRtGjR/B5Ogbdnzx5WrVrFtGnTKJfF96vI7OrVq0yfPp2OHTvSrl07eU8+QmxsLGPHjqVixYr07t1bviMfwWAwsHHjRiIiIvjwww9xcnLK7yEVeAsXLuTs2bNMmTJFnq/HcPr0aTZu3MjgwYMpVapUfg+nwLt48SK+vr4MHDiQV155Jb+HUyjo9Xqio6N59dVX5TNfFAgao9FofPRiz87Jkyd55ZVXGGxrSxmdjsvm5vxoY/PQcl/euIEudegN9HpSwz4YDAb+Mhr5GWik0fCKVnv/PkhPXYvWaLhoZsa8EiWINREcGnTnDtXi4ymVkpLl/f9otfxmYcFqE2/misnJvHXvHo1iY7Ew9XjNzNhta8sxEydmHePiaBwXR+2EhCzvv6PRcCmXj+NCSgq/ABY6He8ClkYjRki/aTQYgDNmZuwrUoRfTYy1UUICDeLjqR0fjxEwZFjfCMQAF8zN2ergYHKsHeLiqJKQgEty8v39GzJsK1yn46K5OQft7LJc3y4lhTpJSbjFxWGW8XGkPgYDcEmn47KVFdcsLbPcRim9ntJJSZRLSgKNBn2G/RuMRm7GxXEpOZkrQC1ra6pYW6ttazSkaDSkpC4fo9Xyn5kZty2y/utrjUZsjEas9HqMgF6jwaDV3l9fD6Skjj0tIFkYJSUlcT012FyiRAmKFCmSzyMq2BITE4mIiADA0dFRTqgfISkpiRs3bmAwGOT5egyJiYncvHkTg8GAhYUFLi4uaDSaR6/4Art79y7/pV7scnZ2xtraOp9HVHAZjUb+++8/YmJiAHBwcMDe3j6fR1WwGY1Gbty4QWJiIvb29jiYuOgpFL1ez7Vr1wAoWrQoxYsXz+cRFXyRkZHExsZiYWFByZIl0RbiY8pn4datW8TFxaHT6ShWrBg2WZyHisz0ej2xsbEsWrSIHj165PdwxHMqLU504sQJ6tevn+2yBTaT62U/P6pWrUp9oNsjlr0HXLl4kW8XLiT51i0qFimC7d9/c9PMjIMtW1J7zJjMK6QGeyobjQRptSYDCNqkJDTx8cSkZkFpDAa0BoO602DAycyM1tbWtDRxQKJJTkYXFwfh4STr9ep3BgOk3QC3okWpWqIEn5n4kjaPjsY8MhLN1av319EYjWhSH4ODuTl1bG2ZVacORhPBlCKhoVheuYLF7dv3t5H2eO7FxBC6ezdF/vuPKlZWFO/aFW1yMhqjEQwGNd7UZZs7OVG3enWTmWs2V69i++efFD11So0xdRtp66dYWFCnRAlad+2KwcRBr9PBgxQ9dQrrf/9VmWdp20kNaCYWLUqsqyv
2015-12-23 17:39:53 +02:00
"text/plain": [
2015-12-31 12:35:58 +02:00
"PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x0000000033496E10>)"
2015-12-23 17:39:53 +02:00
]
},
2015-12-30 17:25:48 +02:00
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
2015-12-31 12:35:58 +02:00
" -0.0361675\n",
" -0.489064 "
2015-12-30 17:25:48 +02:00
]
},
2015-12-31 12:35:58 +02:00
"execution_count": 76,
2015-12-23 17:39:53 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
2015-12-30 17:25:48 +02:00
"function create_problems_2()\n",
"\n",
" mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/2d_beam/BEAM.med\")\n",
" \n",
" upper_body_nodes = Set()\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :TR3 || continue\n",
" push!(upper_body_nodes, elcon...)\n",
" end\n",
" for (nid, coords) in mesh[\"nodes\"]\n",
" nid in upper_body_nodes || continue\n",
2015-12-31 12:35:58 +02:00
" mesh[\"nodes\"][nid] = mesh[\"nodes\"][nid] + [0.0, 1.0]\n",
2015-12-30 17:25:48 +02:00
" end\n",
" \n",
" field_problem = PlaneStressLinearElasticityProblem()\n",
"\n",
" # field problems\n",
" joo = Dict(:QU4 => Quad4, :TR3 => Tri3)\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype in keys(joo) || continue\n",
" element = joo[eltype](elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
" element[\"youngs modulus\"] = 900.0\n",
" element[\"poissons ratio\"] = 0.25\n",
" push!(field_problem, element)\n",
" end\n",
"\n",
" # neumann boundary condition -1 on y direction\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :SE2 || continue\n",
" elset == :LOAD || continue\n",
" element = Seg2(elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
" # FIXME\n",
" # element[\"displacement traction force\"] = rmat*[0.0, -0.01]\n",
2015-12-31 12:35:58 +02:00
" f = rmat*[0.0, -0.0125]*1.5\n",
2015-12-30 17:25:48 +02:00
" element[\"displacement traction force\"] = Vector{Float64}[f, f]\n",
" push!(field_problem, element)\n",
" end\n",
"\n",
" # boundary conditions\n",
" boundary_problem = DirichletProblem(\"displacement\", 2)\n",
"\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :SE2 || continue\n",
" elset == :LEFT || continue\n",
" element = Seg2(elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
" element[\"displacement\"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])\n",
" push!(boundary_problem, element)\n",
" end\n",
"\n",
" info(\"created $(length(get_elements(field_problem))) field elements.\")\n",
" info(\"created $(length(get_elements(boundary_problem))) boundary elements.\")\n",
"\n",
" # Contact definition: contact pair is `LOWER_TO_UPPER <--> UPPER_TO_LOWER`:\n",
"\n",
2015-12-31 08:40:29 +02:00
" master_surface = :UPPER_TO_LOWER\n",
" slave_surface = :LOWER_TO_UPPER\n",
2015-12-30 17:25:48 +02:00
"\n",
2015-12-31 12:35:58 +02:00
" contact_problem = MortarProblem(\"displacement\", 2)\n",
2015-12-30 17:25:48 +02:00
" master_elements = JuliaFEM.Core.Element[]\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :SE2 || continue\n",
2015-12-31 08:40:29 +02:00
" elset == master_surface || continue\n",
2015-12-30 17:25:48 +02:00
" element = Seg2(elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
" element[\"displacement\"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])\n",
" push!(master_elements, element)\n",
2015-12-31 12:35:58 +02:00
" push!(contact_problem, element)\n",
2015-12-30 17:25:48 +02:00
" end\n",
"\n",
" for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
" eltype == :SE2 || continue\n",
" elset == slave_surface || continue\n",
" element = Seg2(elcon)\n",
" update!(element, \"geometry\", mesh[\"nodes\"])\n",
" element[\"master elements\"] = master_elements\n",
" element[\"displacement\"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])\n",
" calculate_normal_tangential_coordinates!(element, 0.0)\n",
" push!(contact_problem, element)\n",
" end\n",
"\n",
" info(\"# of master elements: $(length(master_elements))\")\n",
" info(\"# of slave elements: $(length(contact_problem.elements))\")\n",
"\n",
" return field_problem, boundary_problem, contact_problem\n",
"\n",
"end\n",
"\n",
2015-12-31 08:40:29 +02:00
"using JuliaFEM.Core: calculate_nodal_vector, Element\n",
2015-12-30 17:25:48 +02:00
"\n",
"function JuliaFEM.Core.postprocess_assembly!(\n",
" assembly::BoundaryAssembly, problem::BoundaryProblem{MortarProblem},\n",
" time::Real)\n",
" info(\"postprocess mortar assembly: peforming PDASS.\")\n",
2015-12-31 08:40:29 +02:00
" \n",
2015-12-30 17:25:48 +02:00
" C1 = sparse(assembly.C1)\n",
" C2 = sparse(assembly.C2)\n",
" dim = size(C1, 1)\n",
2015-12-31 12:35:58 +02:00
" elements = get_elements(problem)\n",
" P = -calculate_normal_tangential_coordinates(elements, time)\n",
" resize!(C1, 156, 156)\n",
" resize!(C2, 156, 156)\n",
" resize!(P, 156, 156)\n",
" la = calculate_nodal_vector(\"reaction force\", 2, elements, time)\n",
" X = calculate_nodal_vector(\"geometry\", 2, elements, time)\n",
" u = calculate_nodal_vector(\"displacement\", 2, elements, time)\n",
2015-12-31 08:40:29 +02:00
" x = X+u\n",
" #info(\"x = \", x)\n",
" #info(\"lambda = \", la)\n",
" gh = -C1*x # weighted gap vector\n",
" info(\"weighted gap = \", round(gh, 2))\n",
2015-12-31 12:35:58 +02:00
" info(\"lambda = \", round(la, 2))\n",
"\n",
" #=\n",
" info(\"size P = \", size(P))\n",
" info(\"size C1 = \", size(C1))\n",
" info(\"size C2 = \", size(C2))\n",
" info(\"size la = \", size(la))\n",
" info(\"size gh = \", size(gh))\n",
" info(\"displacement = \", round(u, 2))\n",
" =#\n",
"\n",
" #nz = sort(unique(rowvals(P)))\n",
" #dump(round(full(P[nz,nz]), 3))\n",
" \n",
2015-12-30 17:25:48 +02:00
" C1 = P*C1\n",
" C2 = P*C2\n",
2015-12-31 08:40:29 +02:00
" gh = P*gh\n",
" la = P*la\n",
2015-12-30 17:25:48 +02:00
" c = 1.0\n",
" # complementarity function\n",
2015-12-31 08:40:29 +02:00
" C = la - clamp(la - c*gh, 0, Inf)\n",
2015-12-31 12:35:58 +02:00
" C[abs(C) .< 1.0e-12] = 0.0\n",
" C[49] = -1 # node in support\n",
" #info(\"complementarity function = \", round(C, 2))\n",
2015-12-31 08:40:29 +02:00
" X2 = reshape(X, 2, round(Int, length(X)/2))\n",
" j = 0\n",
2015-12-30 17:25:48 +02:00
" for i=1:2:dim\n",
2015-12-31 08:40:29 +02:00
" j += 1\n",
2015-12-31 12:35:58 +02:00
" X2[1,j] == 0 && continue\n",
" info(\"dof $i: normal gap = $(round(gh[i], 3)), X=$(round(X2[:,j], 2)), C=$(round(C[i], 5)), la=$(round(la[i], 5))\")\n",
" #if C[i] > 0\n",
" if i == 45\n",
" info(\"contact dof $i active (C = $(C[i]))\")\n",
" #gh[i] = -gh[i]\n",
" else\n",
" C1[i, :] = 0\n",
" C2[i, :] = 0\n",
" gh[i] = 0\n",
" end\n",
2015-12-30 17:25:48 +02:00
" end\n",
" for i=2:2:dim\n",
2015-12-31 08:40:29 +02:00
" C1[i, :] = 0\n",
" C2[i, :] = 0\n",
2015-12-31 12:35:58 +02:00
" gh[i] = 0\n",
2015-12-30 17:25:48 +02:00
" end\n",
2015-12-31 08:40:29 +02:00
" #assembly.g = g\n",
2015-12-30 17:25:48 +02:00
" assembly.C1 = C1\n",
" assembly.C2 = C2\n",
2015-12-31 12:35:58 +02:00
" assembly.g = sparse(gh)\n",
2015-12-30 17:25:48 +02:00
" info(\"postprocess mortar assembly: done.\")\n",
"end\n",
"\n",
"field_problem, boundary_problem, contact_problem = create_problems_2()\n",
"using JuliaFEM.Core: DirectSolver\n",
"solver = DirectSolver()\n",
"solver.name = \"divided_beam_inequality_constraints\"\n",
"solver.method = :UMFPACK\n",
2015-12-31 12:35:58 +02:00
"solver.max_iterations = 5\n",
2015-12-30 17:25:48 +02:00
"solver.dump_matrices = false\n",
"push!(solver, field_problem)\n",
"push!(solver, boundary_problem)\n",
"push!(solver, contact_problem)\n",
"call(solver, 0.0)\n",
"\n",
2015-12-31 12:35:58 +02:00
"u_tip = plot(field_problem, 1)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4x4 Array{Float64,2}:\n",
" 0.340028 0.654971 0.587056 0.0\n",
" 0.28579 0.626868 0.0744826 0.0\n",
" 0.215966 0.159712 0.301785 0.0\n",
" 0.0 0.0 0.0 0.0"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = sparse(rand(3, 3))\n",
"resize!(A, 4, 4)\n",
"full(A)"
2015-12-30 17:25:48 +02:00
]
},
{
"cell_type": "code",
2015-12-31 08:40:29 +02:00
"execution_count": 34,
2015-12-30 17:25:48 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
2015-12-31 08:40:29 +02:00
" 0.0\n",
" 0.0"
2015-12-30 17:25:48 +02:00
]
},
2015-12-31 08:40:29 +02:00
"execution_count": 34,
2015-12-30 17:25:48 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
2015-12-31 08:40:29 +02:00
"contact_problem.elements[5](\"displacement\", [0.0], 0.0)"
2015-12-30 17:25:48 +02:00
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Int64,1}:\n",
" 137\n",
" 138\n",
" 139\n",
" 140"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JuliaFEM.Core.get_gdofs(contact_problem.elements[5], 2)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: assemble: doing postprocess for problem JuliaFEM.Core.BoundaryProblem{JuliaFEM.Core.MortarProblem} assembly\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Array(Float64,(0,2)) "
]
},
{
"data": {
"text/plain": [
"4x1 Array{Float64,2}:\n",
" 0.0\n",
" -2.0\n",
" 0.0\n",
" -2.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: postprocess mortar assembly: remove contraints in tangent direction on boundary.\n",
"INFO: postprocess mortar assembly: done.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"0x2 Array{Float64,2}\n"
]
}
],
"source": [
"function JuliaFEM.Core.postprocess_assembly!(\n",
" assembly::BoundaryAssembly, problem::BoundaryProblem{MortarProblem},\n",
" time::Real)\n",
" info(\"postprocess mortar assembly: remove contraints in tangent direction on boundary.\")\n",
" C1 = sparse(assembly.C1)\n",
" C2 = sparse(assembly.C2)\n",
" g = sparse(assembly.g)\n",
" n = round(Int, length(g)/2)\n",
" dump(reshape(full(g), 2, n)[:,60:end]')\n",
" dim = size(C1, 1)\n",
" P = calculate_normal_tangential_coordinates(get_elements(problem), time)\n",
" C1 = P*C1\n",
" C2 = P*C2\n",
" for i=2:2:dim\n",
" C1[i,:] = 0\n",
" C2[i,:] = 0\n",
" end\n",
" assembly.g = g\n",
" assembly.C1 = C1\n",
" assembly.C2 = C2\n",
" info(\"postprocess mortar assembly: done.\")\n",
"end\n",
"\n",
"el1 = Seg2([1, 2])\n",
"el1[\"geometry\"] = Vector{Float64}[[0.0, 0.0], [4.0, 0.0]]\n",
"el2 = Seg2([3, 4])\n",
"el2[\"geometry\"] = Vector{Float64}[[0.0, 1.0], [4.0, 1.0]]\n",
"el1[\"master elements\"] = Element[el2]\n",
"JuliaFEM.Core.calculate_normal_tangential_coordinates!(el1, 0.0)\n",
"p = MortarProblem(\"displacement\", 2)\n",
"push!(p, el1)\n",
"a = JuliaFEM.Core.assemble(p, 0.0)\n",
"full(a.g)"
2015-12-23 17:39:53 +02:00
]
},
2015-12-29 17:50:41 +02:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simple test case, sliding contact"
]
},
2015-12-23 17:39:53 +02:00
{
"cell_type": "code",
2015-12-29 17:50:41 +02:00
"execution_count": 1,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
2015-12-29 17:50:41 +02:00
"Dict{Int64,Array{Float64,1}} with 8 entries:\n",
" 7 => [2.0,0.0]\n",
" 4 => [2.0,4.0]\n",
" 2 => [2.0,6.0]\n",
" 3 => [2.0,2.0]\n",
2015-12-23 17:39:53 +02:00
" 5 => [0.0,0.0]\n",
2015-12-29 17:50:41 +02:00
" 8 => [4.0,0.0]\n",
2015-12-23 17:39:53 +02:00
" 6 => [2.0,0.0]\n",
" 1 => [0.0,2.0]"
]
},
2015-12-29 17:50:41 +02:00
"execution_count": 1,
2015-12-23 17:39:53 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM.Core: Seg2, Tri3, Quad4, Node, update!,\n",
" calculate_normal_tangential_coordinates!, PlaneStressLinearElasticityProblem, \n",
" ContactProblem, SmallSlidingContact, Element\n",
"using JuliaFEM.Core: PlaneStressLinearElasticityProblem, DirichletProblem,\n",
" get_connectivity, Quad4, Tri3, Seg2, LinearSolver,\n",
2015-12-29 17:50:41 +02:00
" update!, get_elements, BiorthogonalBasis, StandardBasis\n",
2015-12-23 17:39:53 +02:00
"nodes = Dict{Int64, Node}(\n",
" 1 => [0.0, 2.0],\n",
2015-12-29 17:50:41 +02:00
" 2 => [2.0, 6.0],\n",
" 3 => [2.0, 2.0],\n",
" 4 => [2.0, 4.0],\n",
2015-12-23 17:39:53 +02:00
" 5 => [0.0, 0.0],\n",
" 6 => [2.0, 0.0],\n",
2015-12-29 17:50:41 +02:00
" 7 => [2.0, 0.0],\n",
" 8 => [4.0, 0.0])"
2015-12-23 17:39:53 +02:00
]
},
{
"cell_type": "code",
2015-12-29 17:50:41 +02:00
"execution_count": 2,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
2015-12-29 17:50:41 +02:00
"basis = StandardBasis\n",
"\n",
"el1 = Quad4([5, 6, 3, 1])\n",
"el2 = Tri3([1, 3, 4])\n",
"el3 = Tri3([7, 8, 2])\n",
2015-12-23 17:39:53 +02:00
"del1 = Seg2([5, 6])\n",
2015-12-29 17:50:41 +02:00
"del2 = Seg2([7, 8])\n",
"sel1 = Seg2([4, 3])\n",
"sel2 = Seg2([3, 6])\n",
"mel1 = Seg2([2, 7])\n",
2015-12-23 17:39:53 +02:00
"update!([el1, el2, el3, sel1, sel2, mel1, del1, del2], \"geometry\", nodes)\n",
"for el in [el1, el2, el3]\n",
" el[\"youngs modulus\"] = 90.0\n",
" el[\"poissons ratio\"] = 0.25\n",
"end\n",
"for el in [del1, del2]\n",
" el[\"displacement\"] = 0.0\n",
"end\n",
"for el in [sel1, sel2]\n",
" el[\"master elements\"] = Element[mel1]\n",
" calculate_normal_tangential_coordinates!(el, 0.0)\n",
"end\n",
"prob = PlaneStressLinearElasticityProblem()\n",
"push!(prob, el1, el2, el3)\n",
2015-12-29 17:50:41 +02:00
"bc = DirichletProblem(\"displacement\", 2; basis=basis)\n",
2015-12-23 17:39:53 +02:00
"push!(bc, del1, del2)\n",
"con = ContactProblem(\"cont\", \"displacement\", 2;\n",
" contact_type=SmallSlidingContact)\n",
"push!(con, sel1, sel2);"
]
},
{
"cell_type": "code",
2015-12-29 17:50:41 +02:00
"execution_count": 3,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
2015-12-29 17:50:41 +02:00
"INFO: slave dofs of element: [7,8,5,6]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: slave dofs of element: [5,6,11,12]\n"
2015-12-23 17:39:53 +02:00
]
},
{
"data": {
"text/plain": [
2015-12-29 17:50:41 +02:00
"JuliaFEM.Core.BoundaryAssembly(JuliaFEM.Core.SparseMatrixCOO([7,5,7,5,7,5,7,5,8,6 … 5,11,6,12,6,12,6,12,6,12],[7,7,5,5,3,3,13,13,8,8 … 13,13,6,6,12,12,4,4,14,14],[0.21522,0.0105929,0.0105929,0.000521371,-0.147011,-0.00723572,-0.0788018,-0.00387854,0.21522,0.0105929 … -0.0109405,-0.222282,0.000521371,0.0105929,0.0105929,0.21522,-0.00017379,-0.00353096,-0.0109405,-0.222282]),JuliaFEM.Core.SparseMatrixCOO([7,5,7,5,7,5,7,5,8,6 … 5,11,6,12,6,12,6,12,6,12],[7,7,5,5,3,3,13,13,8,8 … 13,13,6,6,12,12,4,4,14,14],[0.21522,0.0105929,0.0105929,0.000521371,-0.147011,-0.00723572,-0.0788018,-0.00387854,0.21522,0.0105929 … -0.0109405,-0.222282,0.000521371,0.0105929,0.0105929,0.21522,-0.00017379,-0.00353096,-0.0109405,-0.222282]),JuliaFEM.Core.SparseMatrixCOO(Int64[],Int64[],Float64[]),JuliaFEM.Core.SparseMatrixCOO(Int64[],Int64[],Float64[]))"
2015-12-23 17:39:53 +02:00
]
},
2015-12-29 17:50:41 +02:00
"execution_count": 3,
2015-12-23 17:39:53 +02:00
"metadata": {},
"output_type": "execute_result"
},
{
2015-12-29 17:50:41 +02:00
"name": "stdout",
2015-12-23 17:39:53 +02:00
"output_type": "stream",
"text": [
2015-12-29 17:50:41 +02:00
"\n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n",
"normal tangential = \n",
"[1.0 0.0\n",
" 0.0 -1.0]\n"
2015-12-23 17:39:53 +02:00
]
}
],
"source": [
"ass1 = assemble(prob, 0.0)\n",
"dbc = assemble(bc, 0.0)\n",
"cbc = assemble(con, 0.0)"
]
},
{
"cell_type": "code",
2015-12-29 17:50:41 +02:00
"execution_count": 4,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
2015-12-29 17:50:41 +02:00
"8x8 Array{Float64,2}:\n",
" 92.0 -15.0 0.0 0.0 -74.0 15.0 0.0 -12.0\n",
" -15.0 62.0 0.0 0.0 15.0 -14.0 -18.0 0.0\n",
" 0.0 0.0 6.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 16.0 0.0 0.0 0.0 0.0\n",
" -74.0 15.0 0.0 0.0 110.0 -15.0 -18.0 12.0\n",
" 15.0 -14.0 0.0 0.0 -15.0 110.0 18.0 -48.0\n",
" 0.0 -18.0 0.0 0.0 -18.0 18.0 18.0 0.0\n",
" -12.0 0.0 0.0 0.0 12.0 -48.0 0.0 48.0"
2015-12-23 17:39:53 +02:00
]
},
2015-12-29 17:50:41 +02:00
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"K = full(ass1.stiffness_matrix)\n",
"dims = size(K)\n",
"K[1:8,1:8]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"16x16 Array{Float64,2}:\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 -12.0 0.0 24.0 0.0 6.0 0.0 0.0 0.0 6.0 0.0 -24.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 -12.0 0.0 24.0 0.0 6.0 0.0 0.0 0.0 6.0 0.0 -24.0 0.0 0.0\n",
" 0.0 0.0 -10.0 0.0 6.0 0.0 12.0 0.0 0.0 0.0 0.0 0.0 -8.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 -10.0 0.0 6.0 0.0 12.0 0.0 0.0 0.0 0.0 0.0 -8.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 -2.0 0.0 6.0 0.0 0.0 0.0 0.0 0.0 12.0 0.0 -16.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 -2.0 0.0 6.0 0.0 0.0 0.0 0.0 0.0 12.0 0.0 -16.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"
]
},
"execution_count": 6,
2015-12-23 17:39:53 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
2015-12-29 17:50:41 +02:00
"fixed = [9, 10, 11, 12, 13, 14, 15, 16]\n",
2015-12-23 17:39:53 +02:00
"ENV[\"COLUMNS\"] = 300\n",
2015-12-29 17:50:41 +02:00
"C1 = full(cbc.C1, dims...)\n",
2015-12-23 17:39:53 +02:00
"C1[abs(C1) .< 1.0e-9] = 0\n",
2015-12-29 17:50:41 +02:00
"#C1[fixed, fixed] = 0\n",
"C1*18"
2015-12-23 17:39:53 +02:00
]
},
{
"cell_type": "code",
2015-12-29 17:50:41 +02:00
"execution_count": 7,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
2015-12-29 17:50:41 +02:00
"16x16 Array{Float64,2}:\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"
2015-12-23 17:39:53 +02:00
]
},
2015-12-29 17:50:41 +02:00
"execution_count": 7,
2015-12-23 17:39:53 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
2015-12-29 17:50:41 +02:00
"D = full(cbc.D, dims...)\n",
"D[abs(D) .< 1.0e-9] = 0\n",
"D[fixed, fixed] = 0\n",
"D"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"16x16 Array{Float64,2}:\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 -12.0 0.0 24.0 0.0 6.0 0.0 0.0 0.0 6.0 0.0 -24.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 -10.0 0.0 6.0 0.0 12.0 0.0 0.0 0.0 0.0 0.0 -8.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 -2.0 0.0 6.0 0.0 0.0 0.0 0.0 0.0 12.0 0.0 -16.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C2 = full(cbc.C2, dims...)\n",
2015-12-23 17:39:53 +02:00
"C2[abs(C2) .< 1.0e-9] = 0\n",
2015-12-29 17:50:41 +02:00
"#C2[fixed, :] = 0\n",
"for i=2:2:16\n",
" C1[i,:] = 0\n",
" C2[i,:] = 0\n",
"end\n",
"C2*18"
2015-12-23 17:39:53 +02:00
]
},
{
"cell_type": "code",
2015-12-29 17:50:41 +02:00
"execution_count": 9,
2015-12-23 17:39:53 +02:00
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
2015-12-29 17:50:41 +02:00
"16x16 Array{Float64,2}:\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"
2015-12-23 17:39:53 +02:00
]
},
2015-12-29 17:50:41 +02:00
"execution_count": 9,
2015-12-23 17:39:53 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"D"
]
2015-12-29 17:50:41 +02:00
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"16x16 Array{Float64,2}:\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 12.0 0.0 6.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 12.0 0.0 6.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.0 0.0 12.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.0 0.0 12.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 12.0 0.0 6.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 12.0 0.0 6.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.0 0.0 12.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.0 0.0 12.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"DC1 = full(dbc.C1, dims...)\n",
"DC2 = full(dbc.C2, dims...)\n",
"DD = full(dbc.D, dims...)\n",
"DC1*18"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"32-element Array{Float64,1}:\n",
" 18.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" ⋮ \n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = zeros(16)\n",
"f[1] = 18.0\n",
"g = zeros(16)\n",
"A = [K (C1+DC1)'; (C2+DC2) (D+DD)]\n",
"b = [f; g]"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: zero line 17\n",
"INFO: zero line 18\n",
"INFO: zero line 19\n",
"INFO: zero line 20\n",
"INFO: zero line 22\n",
"INFO: zero line 24\n"
]
}
],
"source": [
"for i=1:size(A,1)\n",
" if sum(abs(A[i,:])) == 0\n",
" info(\"zero line $i\")\n",
" A[i,i] = 1.0\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x8 Array{Float64,2}:\n",
" 0.406155 0.659049 0.219683 0.439366 0.0 0.0 0.0 0.0\n",
" 0.149452 0.0 -0.102834 -0.0562157 0.0 0.0 0.0 0.0"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = A \\ b\n",
"u = reshape(u, 2, 16)\n",
"u[abs(u) .< 1.0e-9] = 0\n",
"disp = u[:,1:8]"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# tie contact\n",
"@assert isapprox(norm(vec(disp)), 0.8223229928825583)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# frictionless sliding\n",
"@assert isapprox(norm(vec(disp)), 0.9363129098080032)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x8 Array{Float64,2}:\n",
" 0.0 2.0 2.0 2.0 0.0 2.0 2.0 4.0\n",
" 2.0 6.0 2.0 4.0 0.0 0.0 0.0 0.0"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = [0.0 2; 2 6; 2 2; 2 4; 0 0; 2 0; 2 0; 4 0]'"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x8 Array{Float64,2}:\n",
" 0.406155 2.65905 2.21968 2.43937 0.0 2.0 2.0 4.0\n",
" 2.14945 6.0 1.89717 3.94378 0.0 0.0 0.0 0.0"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = X + disp"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Int64,1}:\n",
" 7\n",
" 8\n",
" 2"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eldofs = Dict()\n",
"eldofs[1] = [5, 6, 3, 1]\n",
"eldofs[2] = [1, 3, 4]\n",
"eldofs[3] = [7, 8, 2]"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVgAAAG+CAYAAADWcRzfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XdYFWfaP/DvAQTFgr1HYyHG1RhbDIqo58SKhahBsSDMpO5m3U3dbMpm07Mpm01b32R/G8Heu8YaCyJC1FhjjBo8aDCg2FBApJzfH8+SaEI5wMw8Z+Z8P9eVi6zAzO37xq/DM89z3zaXy+UCERFpzkd2AUREVsWAJSLSCQOWiEgnDFgiIp0wYImIdMKAJSLSCQOWiEgnDFgiIp0wYImIdMKAJSLSCQOWiEgnDFgiIp0wYImIdMKAJSLSCQOWiEgnDFgiIp0wYImIdMKAJUs4efIkoqKicP36ddmlEP2MAUuWcODAASxatAirVq2SXQrRzxiwZClxcXGySyD6GQOWLGXTpk04c+aM7DIq78YNYMUKICtLdiWkIQYsWYqvry9mzZolu4zKS0sDxo0DDh2SXQlpiAFLljJq1CjExcWhuLhYdimV43SKj7ffLrMK0hgDlixFURSkpqZi586dskupHKcT8PEBrl8HevQAvvtOdkWkAQYsWUpYWBg6duyImTNnyi6lcpxOoHVroG1b4NtvgW3bZFdEGmDAkqXYbDbExsZi6dKlyM7Oll2O+06dEssDtWsD994LbN0quyLSAAOWLCcmJgZ5eXlYvHix7FLc53T+sv5qtwPbtwNmW0em32DAkuW0bt0aw4YNM9ee2F8H7IULwOHDMisiDTBgyZIURUFSUhKOHTsmu5SK3bgBZGf/ErB9+wIBAVwmsAAGLFlSREQEGjZsaI6nWH9/4OpVIDpa/O+aNYF+/fiiywIYsGRJAQEBmDJlCmbPno3CwkLZ5VTMZgP8/H753w4HsGMHYIbaqUwMWLIsRVGQkZGBDRs2yC6l8iZPBhYulF0FVRMDliyrR48e6N69u/n2xAJA+/bAiBG3PtWS6TBgydIURcGaNWtw/vx52aWQF2LAkqVNmTIFPj4+mDt3ruxSyAsxYMnSGjVqhIiICMTFxcHlcskuh7wMA5YsT1EUHD58GPv27ZNdCnkZBixZ3tChQ9GqVSvPfNl18iQQGsruWRbFgCXL8/X1RUxMDBYsWIC8vDzZ5dzq5EkgKQkIDJRdCemAAUteITY2FpcvX8bKlStll3IrpxPw9QVatSr989nZQGwskJxsZFWkEQYseYXg4GCEhYV53jKB0wncdlvZ+13r1AHWrgW+/NLQskgbDFjyGoqi4KuvvkJaWprsUn7hdALt2pX9eR8fYNAg9iUwKQYseY3IyEgEBgZ61lDEkkbb5bHbgZQUICfHkJJIOwxY8hp16tTBxIkTER8f7zlDEW/uA1sWhwMoKAB27TKiItIQA5a8iqIoOHXqFHbs2CG7FCA3Fzh3ruKAvfNOoHlz9oc1IQYseZXQ0FAEBwd7xsuuwkLg5ZeBe+4p/+tsNq7DmhQDlryKzWaDqqpYtmwZrly5IreYevWAV18FOneu+GsdDmDvXkB2zVQpDFjyOtOmTUN+fj4WLVokuxT3DRsGvPUWByGaDAOWvE7Lli0xfPhwz1gmcFebNsBzzwENGsiuhCqBAUteSVEUpKSk4OjRo7JLIQtjwJJXGj16NBo1amSOoYhkWgxY8koBAQGYOnUq5syZg4KCAtnlkEUxYMlrKYqCzMxMrF+/XnYpZFEMWPJad999N3r27CnnZdfly6JN4fXrxt+bDMOAJa+mqirWrVuHzMxMY2+8a5dotJ2VZex9yVAMWPJqkyZNkjMU0ekEatQAWrSo3Pe5XMA//wkkJOhSFmmLAUterWHDhhg7dixmzpxp7FBEp1PsbfX1rdz32WzAF18A8+bpUhZpiwFLXk9RFBw9ehR79uwx7qbudNEqi93Oxi8mwYAlrzd48GC0bt3a2JddFTXaLo/DIWZ5nTmjaUmkPQYseT1fX1/ExsZiwYIFyM3NNeam1XmCHThQfGR3LY/HgCWCGIqYnZ2NFStW6H+za9fE7oGqBmzjxsDddzNgTYABSwSgQ4cOGDhwoDHLBCUzwaoasMAv67BGvpijSmPAEv2PqqrYunUrnE6nvjfq0kX0da2o0XZ57Hbg9Gkx04s8FgOW6H/Gjx+PunXrIj4+Xv+b1asH+PtX/fsHDAAmTxazushjMWCJ/qd27dqeNxSxLPXri72wnTrJroTKwYAluomqqkhLS8M2vkAiDTBgiW4SEhKCTp06mWvaAXksBizRTUqGIi5fvhyXL1+WXQ6ZHAOW6Feio6NRUFCAhQsXyi6FTI4BS/QrLVq0wIgRI7hMQNXGgCUqhaqq2LNnD44cOaLthRMSgHHjxGkusjwGLFEpRo4cicaNG2s/FPHAAeDLL4HAQO2uuX27mI5AHocBS1QKf39/REdHaz8U0ekE2rYFfDT8o/fqq8C772p3PdIMA5aoDKqq4vz581i3bp12F61OF62yOBziKbaoSNvrUrUxYInK0LVrV/Tu3Vvbl116BKzdLnobHDig7XWp2hiwROVQVRVffvklMjIytLmgHgHbp49Y0+WUA4/DgCUqR1RUFPz8/DBnzpzqX+zKFeDSpapPMiiLvz/Qvz/7w3ogBixRORo0aIBx48ZpMxSxpA2i1k+wgFgm2LmT3bU8DAOWqAKqquLYsWNISUmp3oUCA4FHHgE6dtSmsJs5HGJv7d692l+bqowBS1QBh8OBNm3aVP9lV3Aw8PnnYuSL1nr2BEJCgKtXtb82VRkDlqgCPj4+iI2NxcKFC5GTkyO7nNL5+QG7dwNDh8quhG7CgCVyQ2xsLK5evYrly5fLLoVMhAFL5IZ27drBbrezAQxVCgOWyE2qqmL79u1ITU2VXQqZBAOWyE3jxo1DvXr1jBmKSJbAgCVyU2BgIKKiohAfH48invsnNzBgiSpBVVWcOXMGWyt7LDUrC0hNBap7WIFMhQFLVAl9+vRB586dK/+ya+5coGtXfYr6tcxM4NAhY+5F5WLAElVCyVDEFStW4NKlS+5/Y0mTF5tNr9J+8ec/ixNjJB0DlqiSoqOjUVhYiAULFrj/TadO6dODoDR2uzgym51tzP2oTAxYokpq1qwZRo4cWbllAj3aFJbF4RDNt3fuNOZ+VCYGLFEVqKqKffv24ZA7a50ul7EB27Ej0KoV2xd6AAYsURWEh4ejadOm7g1FvHxZ/LhuVMDabOIplgErHQOWqApq1KiB6OhozJ07Fzdu3Cj/i0v6wGrdaLs8djuwfz9w8aJx96TfYMASVZGiKMjKysLatWvL/0I9G22XxeEQSxMJCcbdk36DAUtURV26dEGfPn0qftk1ciTw/ff69IEtS9u2QOfOQFqacfek3/CTXQCRmamqij/84Q84e/YsWrZsWfoX+fsDd9xhbGEAcPgw4Otr/H3pZ3yCJaqGqKgo+Pv7azMUUWsMV+kYsETVEBQUhPHjx2szFJEshwFLVE2qquL48eNISkqSXQp5GAYsUTUNGjQIt99+u3t7YsmrMGCJqqlkKOKiRYs8dygiScGAJdJAbGwscnJysHTpUtmlkAdhwBJpoG3btnA4HL/dE7tiBfDMM3KKKuFyARkZcmvwUgxYIo2oqoqEhAScPHnyl1/cuhXYsEFeUQDwwgtA375ya/BSDFgijYwdOxZBQUG3DkV0Oo3tQVCakBBRx6lTcuvwQgxYIo3UqlULkyZNunUoopFtCssyYIDosMXuWoZjwBJpSFVVpKenY/PmzWLt08hJBmVp0ADo2ZMBKwEDlkhDvXv3RpcuXcSe2AsXgJwc+QELiPaFW7dyqq3BGLBEGioZirhy5UpcOXhQ/KInBKzDAZw9C5w4IbsSr8KAJdLY1KlTUVxcjKT588UveELA9u8vmr9s3Sq7Eq/CgCXSWNOmTTFq1CisTEgAIiOBhg1llwTUrQv06QMkJ8uuxKuwHyyRDlRVxZgxY/D7JUvQ3WaTXY6wcqWxTb+JT7BEehgxYgSaN2/uWQ1gmjYFfPhH3kj8vzaRDvz8/H4eipifny+7HJKEAUukE0VRcPHiRaxevVp2KSQJA5ZIJ507d0b
"text/plain": [
"PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x000000002615EAC8>)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"(-0.1,6.1)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using PyPlot\n",
"\n",
"function get_coords(X)\n",
" return [X X[:, 1]]\n",
"end\n",
"\n",
"function plot(X, args...; kwargs...)\n",
" X1 = X[:, eldofs[1]]\n",
" X2 = X[:, eldofs[2]]\n",
" X3 = X[:, eldofs[3]]\n",
" PyPlot.plot(get_coords(X1)[1,:]', get_coords(X1)[2,:]', args...; kwargs...)\n",
" PyPlot.plot(get_coords(X2)[1,:]', get_coords(X2)[2,:]', args...; kwargs...)\n",
" PyPlot.plot(get_coords(X3)[1,:]', get_coords(X3)[2,:]', args...; kwargs...)\n",
" axis(\"off\")\n",
"end\n",
"figure(figsize=(4, 5))\n",
"plot(X, \"-k\")\n",
"plot(x, \"--r\")\n",
"xlim(-0.1, 4.1)\n",
"ylim(-0.1, 6.1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
2015-11-24 03:06:56 +02:00
}
],
"metadata": {
"kernelspec": {
2015-12-31 12:35:58 +02:00
"display_name": "Julia 0.4.1",
2015-11-24 03:06:56 +02:00
"language": "julia",
"name": "julia-0.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
2015-12-31 12:35:58 +02:00
"version": "0.4.1"
2015-11-24 03:06:56 +02:00
}
},
"nbformat": 4,
"nbformat_minor": 0
}