JuliaFEM.jl/notebooks/2015-12-13-3d-tie-contact.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tie contact 3d\n",
    "\n",
    "Author(s): Jukka Aho"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solid block\n",
    "\n",
    "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>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using JuliaFEM\n",
    "using JuliaFEM.Preprocess: parse_aster_med_file\n",
    "using JuliaFEM.Core: LinearElasticityProblem, DirichletProblem, get_connectivity, Quad4, Hex8, DirectSolver"
   ]
  },
  {
   "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",
   "execution_count": 5,
   "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",
   "execution_count": 6,
   "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",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: # of field problems: 1\n",
      "INFO: # of boundary problems: 1\n",
      "INFO: Starting iteration 1\n",
      "INFO: Assembling 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"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1,true)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "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"
     ]
    }
   ],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: nid near corner = 7\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "7"
      ]
     },
     "execution_count": 8,
     "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",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "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"
     ]
    }
   ],
   "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",
    "    @test isapprox(u[i], known_value)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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",
    "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-17-divided-block-both-parts.png\"\n",
    "         style=\"width:100%; height:auto; vertical-align:middle;\">\n",
    "  </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Found 6 element sets: SYM23, SYM12, UPPER_TO_LOWER, LOAD, SYM13, LOWER_TO_UPPER\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Dict{ASCIIString,Any} with 2 entries:\n",
       "  \"nodes\"        => Dict(11=>[1.0,0.0,0.5],134=>[0.5,0.25,1.0],158=>[1.0,0.5,0.…\n",
       "  \"connectivity\" => Dict(288=>(:HE8,:OTHER,[176,177,180,179,167,168,171,170]),3…"
      ]
     },
     "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,\n",
    "                     Quad4, Hex8, LinearSolver, update\n",
    "mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/BLOCKS.med\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: created 107 elements.\n",
      "INFO: created 59 boundary elements.\n"
     ]
    }
   ],
   "source": [
    "# interior elements are of type HE8\n",
    "field_problem = LinearElasticityProblem()\n",
    "\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    eltype == :HE8 || continue\n",
    "    element = Hex8(elcon)\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    element[\"youngs modulus\"] = 900.0\n",
    "    element[\"poissons ratio\"] = 0.25\n",
    "    #info(\"lower: add element with connectivity $elcon\")\n",
    "    push!(field_problem, element)\n",
    "end\n",
    "\n",
    "# Neumann boundary condition, traction force -100 on Z direction for element set LOAD\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    (eltype == :QU4) && (elset == :LOAD) || continue\n",
    "    element = Quad4(elcon)\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    element[\"displacement traction force 3\"] = -100.0\n",
    "    push!(field_problem, element)\n",
    "end\n",
    "info(\"created $(length(field_problem.elements)) elements.\")\n",
    "\n",
    "# boundary conditions\n",
    "boundary_problem = DirichletProblem(\"displacement\", 3)\n",
    "\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    (eltype == :QU4) && (elset in [:SYM23, :SYM12, :SYM13]) || continue\n",
    "    element = Quad4(elcon)\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    if elset == :SYM23\n",
    "        element[\"displacement 1\"] = 0.0\n",
    "    elseif elset == :SYM12\n",
    "        element[\"displacement 3\"] = 0.0\n",
    "    elseif elset == :SYM13\n",
    "        element[\"displacement 2\"] = 0.0\n",
    "    end\n",
    "    push!(boundary_problem, element)\n",
    "end\n",
    "info(\"created $(length(boundary_problem.elements)) boundary elements.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 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": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: using PyPlot.mesh in module Main conflicts with an existing identifier.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAG7CAYAAABEoAutAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt0VOW9N/Dv3plhkiEXEAhMQEuMjYCXhEBSlJyIeI66qByPNIhpZUHl4nuaHgERwbSAAo1IoUFsWAUvaVCbYlNKz/L1Lb7U0zeF5YJc0dJEEANIMmEUIRcnGWZmP+8fk4zmPkPm+vD9rJXF4pnZ2b/n90z2N3tndqIIIQSIiIgkoAa7ACIiIl9hqBERkTQYakREJA2GGhERSYOhRkRE0mCoERGRNBhqREQkDV2wCxjIl19+iUOHDmHixImIiooKdjlERBQE7e3tOHv2LB544AGMHj16wOeGdKgdOnQIjz/+eLDLICKiEPDWW2/hRz/60YDPCelQmzhxIgDXRCZPnhyUGlatWoWCgoKg7DtUXO89uN7nD7AH1/v8geD2oLa2Fo8//rg7EwYS0qHWdclx8uTJSEtLC0oNcXFxQdt3qLjee3C9zx9gD673+QOh0QNPfgzFN4oQEZE0GGpERCQNhhoREUmDoTaInJycYJcQdNd7D673+QPswfU+fyB8eqCE8t9Tq6qqwrRp01BZWRn0H1ASEVFweJMFPFMjIiJpMNSIiEgaDDUiIpIGQ42IiKTBUCMiImkw1IiISBoMNSIikgZDjYiIpMFQIyIiaTDUiIhIGgw1IiKSBkONiIikwVAjIiJpMNSIiEgaDDUiIpIGQ42IiKTBUCMiImkw1IiISBpehdrXX3+NjRs34sEHH8QNN9wAVVVRXFzs8fZXrlzB8uXLMWbMGERHR2P27Nmorq72umgiIqK+eBVqX3zxBTZv3oxPPvkEqampAABFUTzaVtM0fP/730dJSQmeeuopbNu2DRaLBbNmzcKnn37qfeVEREQ96Lx5ckJCApqamhAfH4/Kykqkp6d7vG1paSk+/PBDlJaWYt68eQCARx99FMnJydi4cSPefvtt7yonIiLqwasztWHDhiE+Ph4AIITwakelpaUYN26cO9AAYPTo0Xj00Ufx5z//GXa73avPR0RE1JNXZ2pDUV1djbS0tF7j6enp2Lt3L06dOoXbbrstUOUMyGKxYHf+y7hwog6qXYOmVzEhZRJ+krfCHeo0MIvFgvz8Ypw48QXsdj30ejtSUsYgL2+Rz3tYX1+PHbm5aKmthcHhgE2nQ+zkyVhdWIjExESf7iuQAtnDQAnkWsnYv0Crr69Hfu7TuFL7CYY5nLiqi8CIybcir/BXIfu1FbBQM5vNmDVrVq9xk8kEAGhsbAyJULNYLFg7bynus45FVvRtUHUqNKeGhrImrK1YipcOvMYviEFYLBbMm7cVVusyREffCp1OhdOpoaysDhUVW3HgwDqf9bC+vh7rs7KwxmbDHZGRUCMioDmd+KiiAuuzsrC5rCxkv/gGEsgeBkog10rG/gVafX09nsl6AP/LNgxJkaOhRriOhZ9WfIpnsh7A9rJDIfm1FbC39Hd0dMBgMPQaj4yMBAC0t7cHqpQB7c5/GbOt4zAhdjxU1dUeVVVxY2wCZlvHYnf+y0GuMPTl5xfDal2K2NjJ3XoYGzsFVusS5Od7/o7ZwezIzcUzNhtSjMZu+0o1GrHaZsOO3Fyf7SuQAtnDQAnkWsnYv0DLz30aT9oM+K5xdLceJhtHYbnNgPzcp4NcYd8CdqYWFRUFm83Wa7yjo8P9eCi4cKIOWdHfOmO0WgGnEwAwHjH4f0dqcLHGHKTqwsOZI/UwIQ5Kq6tPTicQEeF6bDhG4MyRz33Wwys1H+M2VQ/nVQcgADjsEBE6QFFwu6pHe83HYblegexhoLR0rZXN4Rpw2IEIHaC61uqKD9fqzJHPB+lffdj1L9Cu1JzELeoNgM0GQECn2eE0GCHUCHw3ciSu1H4S7BL7FLBQM5lMaGxs7DVuNrteWAkJCf1uu2rVKsTFxXUby8nJQU5Ojm+LBKDaNai6zhNYqxUjTh3v9vgo0YTalXt8vl+ZzD79MRTlVQCApgFXrwLDhgGd3+xBiFM+66F6qRlORYUTAIQGnfMqnKoOQnW9tMcLLSzXK5A9DBTjIGul+nCtZp8+NUj/Pg67/gVa3KVLMCpXewwCjqgYqKqKYQ6nX/ZbUlKCkpKSbmPNzc0ebx+wUEtNTcXf//53CCG63dt27NgxDB8+HMnJyf1uW1BQ0OebTPxB07uuG6uq6j5Daxk/GVqkEZoQuKR+gsk7nwxILeFqx1ILNG0ZFEVFeztw4QIwYQIQFQUIoUFVv8IPfdRDbc7riNA0qKoCXLVDbb0KJWY4hF4PIQQaVDUs1yuQPQwUa9daKQpg/2atMEwPTRPQfLhWO5Z+NUj/LGHXv0BrnvM7WLUboCoKFMdVGFovwaHqAbjuO76qi/DLfvs6YamqqsK0adM82t4vodbU1IQrV67glltugU7n2kV2djZKS0tx4MAB/OAHPwAAfPnll/jDH/6AuXPnQq/X+6MUr01ImYSGsibcGPvNmaMWaQRiYtDQ0oAJmakYm2oKYoWhLykzEWVlzYiNmYJWABdVIDoKcMYALS0nkZV5o896OCL1DpysqEBqpBEQTkABIvQqYNCh2mpFVNr0sFyvQPYwUGK71spoBPCttRqmw0dWK0b4cK2SMm8cpH+JYde/QBuRehs+rfgUycZR+OYGLtcJyemOyxgx/dZglTYgr98o8utf/xpbtmzBG2+8AQD47//+b2zZsgVbtmxBS0sLAGDdunWYMmVKt8uN2dnZmDFjBn784x9j8+bN2L17N2bNmgUhBF544QUfTWfofpK3Ah8Ym/B5SwO0znvxNCHweUsDPjBexE/yVgS5wtCXl7cIRuNraGk5CSE0AK7vjltaTsJofB15eYt8tq/VhYXYbjCg2mrttl7VVit2GAxYXVjos30FUiB7GCiBXCsZ+xdoeYW/wh6DDZ9YL3Vbr0+sl7DXYENe4a+CXGHfFOHlXdSJiYk4d+6ca+POy4hdlxTr6+tx00034cc//jH27dvn/n+XK1euYM2aNTh48CDa29uRkZGB7du393tpseuUs7KyMmCXH4Fv3ad2pAajTp/Bpe8mYUJmKu9T80LXPUJHjlzG6dPx+O53LcjMHOnX+9RsNTWYeOkSzo4aBUNqqjT3qQWih4HStVZtNTUYcekSrowahWg/rZWM/Qs0931qNScRd+kSmkeNwojU2wJ+n5pXWSBCWGVlpQAgKisrg7L/pupG8T/3bBRN1Y1B2b8MqquFuOce179y7SxwpJxWACclZf8CLNjHQm+ygH96hoiIpMFQIyIiaTDUiIhIGgw1IiKSBkONiIikwVAjIiJpMNSIiEgaDDUiIpIGQ42IiKTBUCMiImkw1IiISBoMNSIikgZDjYiIpMFQIyIiaTDUiIhIGgw1IiKSBkONiIikwVAjIiJpMNSIiEgaDDUiIpIGQ42IiKTBUCMiImkw1IiISBoMNSIikgZDjYiIpMFQIyIiaTDUiIhIGgw1IiKSBkONiIikwVAjIiJpMNSIiEgaDDUiIpIGQ42IiKTBUCMiImkw1IiISBoMNSIikgZDjYiIpMFQIyIiaTDUiIhIGgw1IiKSBkONiIikwVAjIiJpMNSIiEgaDDUiIpIGQ42IiKTBUCMiImkw1IiISBoMNSIikgZDjYiIpMFQIyIiaTDUiIhIGgw1IiKSBkONiIik4VWo2Ww2rF27FgkJCTAajZgxYwYOHz7s0baVlZV46KGHYDKZEBMTg5SUFLzyyivQNO2aCiciIurJq1BbvHgxCgoKsHDhQuzatQsRERGYM2cOjh49OuB2lZWVuPvuu3H+/HmsW7cOv/rVr3DzzTdjxYoVePrpp4c0ASIioi46T594/Phx7N+/H9u3b3cH0cKFC3H77bfj2WefHTDY9uzZA1VVUVZWhhEjRgAAli1bhlmzZuG3v/0tdu7cOcRpEBEReXGmVlpaCp1Oh+XLl7vHDAYDlixZgg8//BANDQ39btvS0gKDwYC4uLhu4+PGjYPRaLyGsomIiHrzONSqq6uRnJyM6OjobuPp6ekAgJq
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f5cf2a922d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(-0.1,1.1)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using JuliaFEM.Core: Element, MortarProblem, calculate_normal_tangential_coordinates!\n",
    "\n",
    "using PyPlot\n",
    "\n",
    "#mortar_surface = :LOWER_TO_UPPER\n",
    "#slave_surface = :UPPER_TO_LOWER\n",
    "mortar_surface = :UPPER_TO_LOWER\n",
    "slave_surface = :LOWER_TO_UPPER\n",
    "\n",
    "fig1 = figure(figsize=(5, 5))\n",
    "master_elements = JuliaFEM.Core.Element[]\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    (eltype == :QU4) && (elset == mortar_surface) || continue\n",
    "    element = Quad4(elcon)\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    X = element(\"geometry\", 0.0)\n",
    "    x = [X[mod(i, 4)+1][1] for i=1:5]\n",
    "    y = [X[mod(i, 4)+1][2] for i=1:5]\n",
    "    plot(x, y, \"-bo\", alpha=0.5)\n",
    "    push!(master_elements, element)\n",
    "end\n",
    "\n",
    "contact_problem = MortarProblem(\"displacement\", 3)\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    (eltype == :QU4) && (elset == slave_surface) || continue\n",
    "    element = Quad4(reverse(elcon))\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    element[\"master elements\"] = master_elements\n",
    "    X = element(\"geometry\", 0.0)\n",
    "    x = [X[mod(i, 4)+1][1] for i=1:5]\n",
    "    y = [X[mod(i, 4)+1][2] for i=1:5]\n",
    "    plot(x, y, \"-ro\", alpha=0.5)\n",
    "    calculate_normal_tangential_coordinates!(element, 0.0)\n",
    "    push!(contact_problem, element)\n",
    "end\n",
    "xlim(-0.1, 1.1)\n",
    "ylim(-0.1, 1.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: # of field problems: 1\n",
      "INFO: # of boundary problems: 2\n",
      "INFO: Starting iteration 1\n",
      "INFO: Assembling field problems...\n",
      "INFO: Assembling body 1...\n",
      "INFO: Assembly: 10.0 % done. \n",
      "INFO: Assembly: 20.0 % done. \n",
      "INFO: Assembly: 30.0 % done. \n",
      "INFO: Assembly: 40.0 % done. \n",
      "INFO: Assembly: 50.0 % done. \n",
      "INFO: Assembly: 60.0 % done. \n",
      "INFO: Assembly: 70.0 % done. \n",
      "INFO: Assembly: 80.0 % done. \n",
      "INFO: Assembly: 90.0 % done. \n",
      "INFO: Assembly: 100.0 % done. \n",
      "INFO: dim = 567\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",
      "INFO: UMFPACK: solved in 0.21409296989440918 seconds. norm = 1.0506403593090503\n",
      "INFO: timing info for iteration:\n",
      "INFO: boundary assembly       : 1.9287428855895996\n",
      "INFO: field assembly          : 1.2453219890594482\n",
      "INFO: dump matrices to disk   : 1.1484379768371582\n",
      "INFO: solve problem           : 0.34978604316711426\n",
      "INFO: update element data     : 0.023519039154052734\n",
      "INFO: non-linear iteration    : 4.695832014083862\n",
      "INFO: solver finished in 4.851454019546509 seconds.\n",
      "INFO: nid near corner = 95\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777818,0.02777777777777777,-0.11111111111111108]\n"
     ]
    }
   ],
   "source": [
    "using JuliaFEM.Core: DirectSolver\n",
    "solver = DirectSolver()\n",
    "solver.name = \"tie_contact_3d\"\n",
    "solver.method = :UMFPACK\n",
    "solver.nonlinear_problem = false\n",
    "solver.max_iterations = 1\n",
    "solver.dump_matrices = true\n",
    "push!(solver, field_problem)\n",
    "push!(solver, boundary_problem)\n",
    "push!(solver, contact_problem)\n",
    "call(solver, 0.0)\n",
    "\n",
    "using JuliaFEM.Test\n",
    "\n",
    "nid = 0\n",
    "for (nid, coords) in mesh[\"nodes\"]\n",
    "    if isapprox(coords, [1.0, 1.0, 1.0])\n",
    "        info(\"nid near corner = $nid\")\n",
    "        break\n",
    "    end\n",
    "end\n",
    "nid\n",
    "\n",
    "known_value = [1/36, 1/36, -1/9]\n",
    "\n",
    "for element in field_problem.elements\n",
    "    i = indexin([nid], get_connectivity(element))[1]\n",
    "    i != 0 || continue\n",
    "    X = element(\"geometry\", 0.0)\n",
    "    u = element(\"displacement\", 0.0)\n",
    "    info(\"displacement X = $(X[i]), u = $(u[i])\")\n",
    "    @test isapprox(u[i], known_value)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777818,0.02777777777777777,-0.11111111111111108]\n",
      "INFO: XDFM: ndim = 567\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/blocks.xmf\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<img src=\"http://results.juliafem.org/2015-12-17-divided-block-results.png\"\n",
    "         style=\"width:80%; height:auto; vertical-align:middle;\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using JuliaFEM\n",
    "using JuliaFEM.Preprocess: parse_aster_med_file\n",
    "using JuliaFEM.Core: LinearElasticityProblem, DirichletProblem, get_connectivity, Tri3, Tet4, LinearSolver, update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Found 6 element sets: SYM23, SYM12, UPPER_TO_LOWER, LOAD, SYM13, LOWER_TO_UPPER\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Dict{ASCIIString,Any} with 2 entries:\n",
       "  \"nodes\"        => Dict(11=>[0.0,0.3333333333333333,0.0],134=>[0.0,0.869378844…\n",
       "  \"connectivity\" => Dict(306=>(:TE4,:OTHER,[53,24,51,83]),1316=>(:TE4,:OTHER,[2…"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/BLOCKS_TET4.med\")\n",
    "#med = JuliaFEM.Preprocess.MEDFile(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/BLOCKS_TET4.med\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: created 955 elements.\n",
      "INFO: created 116 boundary elements.\n"
     ]
    }
   ],
   "source": [
    "# interior elements are of type HE8\n",
    "field_problem = LinearElasticityProblem()\n",
    "\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    eltype == :TE4 || continue\n",
    "    element = Tet4(elcon)\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    element[\"youngs modulus\"] = 900.0\n",
    "    element[\"poissons ratio\"] = 0.25\n",
    "    #info(\"lower: add element with connectivity $elcon\")\n",
    "    push!(field_problem, element)\n",
    "end\n",
    "\n",
    "# Neumann boundary condition, traction force -100 on Z direction for element set LOAD\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    (eltype == :TR3) && (elset == :LOAD) || continue\n",
    "    element = Tri3(elcon)\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    element[\"displacement traction force 3\"] = -100.0\n",
    "    push!(field_problem, element)\n",
    "end\n",
    "info(\"created $(length(field_problem.elements)) elements.\")\n",
    "\n",
    "# boundary conditions\n",
    "boundary_problem = DirichletProblem(\"displacement\", 3)\n",
    "\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    (eltype == :TR3) && (elset in [:SYM23, :SYM12, :SYM13]) || continue\n",
    "    element = Tri3(elcon)\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    if elset == :SYM23\n",
    "        element[\"displacement 1\"] = 0.0\n",
    "    elseif elset == :SYM12\n",
    "        element[\"displacement 3\"] = 0.0\n",
    "    elseif elset == :SYM13\n",
    "        element[\"displacement 2\"] = 0.0\n",
    "    end\n",
    "    push!(boundary_problem, element)\n",
    "end\n",
    "info(\"created $(length(boundary_problem.elements)) boundary elements.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: using PyPlot.mesh in module Main conflicts with an existing identifier.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAG7CAYAAABEoAutAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzsnXd0FFUbxp+ZbAohJBBDh4QiKL2G3hUpotIEEVSUKgLSm/QiCAYUBKUIHwgiiAiKIohIkx4CmNBbooQWWgjpO/f748mkbpIts9kkzu8cDrC7M3On3fe+XRJCCOjo6Ojo6OQDZEcPQEdHR0dHRyt0oaajo6Ojk2/QhZqOjo6OTr5BF2o6Ojo6OvkGXajp6Ojo6OQbdKGmo6Ojo5Nv0IWajo6Ojk6+weDoAWRFREQEdu3ahXLlyqFAgQKOHo6Ojo6OjgOIiYnBjRs30K5dO/j4+GT521wt1Hbt2oU+ffo4ehg6Ojo6OrmA9evXo3fv3ln+JlcLtXLlygHgiVSpUsUhYxg5ciQWLVrkkGPnFv7r1+C/fv6Afg3+6+cPOPYanD9/Hn369EmWCVmRq4WaanKsUqUK6tat65AxeHl5OezYuYX/+jX4r58/oF+D//r5A7njGpjjhtIDRXR0dHR08g26UNPR0dHRyTfoQk1HR0dHJ9+gC7Vs6NWrl6OH4HD+69fgv37+gH4N/uvnD+SdayDl5n5qp06dQr169RAYGOhwB6WOjo6OjmOwRBbompqOjo6OTr5BF2o6Ojo6OvkGXajp6Ojo6OQbdKGmo6Ojo5Nv0IWajo6Ojk6+QRdqOjo6Ojr5Bl2o6ejo6OjkG3ShpqOjo6OTb9CFmo6Ojo5OvkEXajo6Ojo6+QZdqOno6Ojo5Bt0oaajo6Ojk2/QhZqOjo6OTr5BF2o6Ojo6OvkGXajp6Ojo6OQbdKGmo6Ojo5Nv0IWajo6Ojk6+QRdqOjo6Ojr5BouE2tOnTzFt2jS0b98e3t7ekGUZa9euNXv7R48eYeDAgShatCg8PDzQpk0bBAUFWTxoHR0dHR0dU1gk1O7du4dZs2bh4sWLqF27NgBAkiSztlUUBS+//DI2btyI4cOHY/78+bh79y5atWqFK1euWD5yHR0dHR2ddBgs+XGpUqVw+/ZtFCtWDIGBgfD39zd72y1btuDIkSPYsmULunbtCgDo0aMHKleujGnTpmHDhg2WjVxHR0dHRycdFmlqLi4uKFasGABACGHRgbZs2YISJUokCzQA8PHxQY8ePbB9+3YkJCRYtD8dHR0dHZ305FigSFBQEOrWrZvhc39/f0RHR+PSpUs5NZRsOXDgAJr6lEI7gys6yy5oZ3BFU59SOHDggKOHlmc4cOAAfHxawWDoAFnuBoOhA3x8WtnlGoaEhKBn7droWqQI3vT0RNciRdCzdm2EhIRofqyc5Pr16+jYcSjKl38bZcsOQPnyb6Njx6G4fv26o4dmNTl5r06cOAE/vw5wc+sEF5fX4ebWCX5+HXDixAnNj5VfCQkJwcu1m6NtkbLo5FkWbYuUxcu1m+fud0tYyYkTJ4QkSWLt2rVm/b5gwYKif//+GT7/5ZdfhCRJYvfu3Rm+CwwMFABEYGCgtcO0mP3794u2krP4Cy7isVRIPHHyEo+lQuIQXEVbyVns378/x8aSV9m/f7+QpPYCOCEkySicnISQJKMAjgtJaq/pNQwODhadCxUSgW5uwliokBBeXsJYqJA46eYmOhcqJIKDgzU7Vk5y7do1UaZMb1G06GlRtqxR+PkJUbasURQtGiTKlOktrl275ughWkxO3qvjx48LZ+eXhSyfFM7ORuHqKoSzs1E4OZ0Qzs4vi+PHj2t2rPxKcHCweKFQSfGD2/MisFBzEeTVUgQWai62uFURLxQqmaPvliWyIMc0tdjYWLi6umb43M3NDQAQExOTU0PJkold38A0IaOmUwHIkgQIAVmWUcvJDVOEjIld33D0EHM9XbtOhRAz4eRUH5IkQwhAlmU4OflDiJno2nWqZsea2bs3JickoK6ra5r7Vc/VFR8lJGBm796aHSsn+eCDAMTFjYG7ey3IMl9TWZbh7l4bcXGj8cEHAQ4eoeWkuVeybNd71b37VCjKdDg714MspzyDBkN9KMo0dO+u3TOYXxnXezDeTyiMCq7F0zyDFV2LYXCCF8b1HuzgEZrGokARWyhQoADi4uIyfB4bG5v8fW7A49F91JAofN2UpwAAI5wgIKEeFJS7fw/7Os535BBzPf3u34XAXsD4Z5rPJQACAtL925pdw7LBF1FDKEhIiIWTMR5CcoKQnQAANRUFsSEX8uT98v/zJOqLssCT3UiIBwwGIGlegRAC0p9577xS3ytZGCEpChTZAEiS5veq1z/3IMQfkJS9UJLc/5KU9AwKAemf6Dx3/XKaCsF/o54oAjkhDK5KNOJlN9z3KMfvnIsiPjR3msFzTKiVLFkS4eHhGT6/desWAEZWZsbIkSPh5eWV5rNevXqhV69e2g4SQAFFJK9KjDDACYmQYUQiXADJAEXEoUjHxpofNz9xZOevAJogfSiRBE4sQvyKNzW6htF7JciSMyQhIBkBCQqEJEEYXOAEoIBQ8uT9OnX4MGS5CaKjgQQJcBJAQTdePwBQlCB0zWPnlXyvFAVyQgIgCciygDA4a36vju76A5LcBEIARlWoCUB2UgXbYfTKY9cvp3m81xmx8ELhhLtwE7GIF27J38myDFejfY67ceNGbNy4Me1YHj82e/scE2q1a9fGwYMHuUpKldt27NgxFCxYEJUrV85020WLFpkMMrEHMbIERVEgyzKMkhOEEJChwIB4xAkXPHCSUWto8xwZS17l8IgCMBqbIn0cEgWaAicnd82u4e0prpCNRshGIyBLgJMTJEUBJAWKwYB4J+c8eb+CA1biyZOmiIIMd08gOhooAMDLkzmfTk4r89x53Z7iCjkxEXJsAu+VJEOSJcDVAEVRNL1XR8c4IyGhKYSQAVqluSBQACcnBU6GT/Lc9ctpbk1xg2vUExRQohEtFUSCU4r7SFEUxDnZ57imFJZTp06hXr16Zm1vF5/a7du3ceHCBSQmJiZ/1r17d9y5cwdbt25N/iwiIgLff/89XnnlFTg7O9tjKBYTVfgZ/C3iU30iIw4FICDjEmJhcHd32NjyCm5uiQBOJWsVkqQKNAAIQuHC2vlPnf38EJSQACgKD1KgAODkBMTHIyguDs5+fpodKyd59lkvPHkSAmdnoFAhwN0diIkB4uOB2NgzqFLF09FDtBhnPz+cj43lg+DmRpuq0QgIgVMJCZreK09PDyhKMCSuc5L/BoDExBAUL24nNSMfUc7VHY8Tb+GB/AyinIuk+e5awj24+JVzzMCywWKh9sUXX2D27NlYvXo1AOCnn37C7NmzMXv2bERGRgIAJkyYgKpVq6YxN3bv3h2NGjXCu+++i1mzZmHZsmVo1aoVhBCYMWOGRqdjO3O3focZkoIgYywUoQAAFADHIGMegNkl/YCvv3boGHMzP/0EVKw4B8A0CHECvHpI+vskgM/Rq5d2QQ5TN2zAbGdnhCQk8EiSBMXFBUEAvkhMxIxhwzQ7Vk5SufJYAEvg6hoEIRQUKgRIkoIHD87BxSUAS5eOdvQQLWZunz7Yoig4K8tQZBlwcoKiKAiMjcUcZ2dM1agAw86dQNmyH0OS5kOSTkAkvcf8+xSAdejYMe8F2uQkj7/ZjnfivPCF5IxAZ5eUuVBRcCXuDr5yfoz5G75y8CgzwdLQynLlyglJkoQkSUKWZSHLcvK/Q0NDhRBC9O3bN83/VR4+fCj69+8vfHx8RMGCBUXr1q2zDNF0REi/EAxJb/JMSdFeMoiucBIvObmIxt4lRfuXDoo1z30sYhq2FOLjj4UwGnN0XLmdw4eFeOEFIUaPFmLfvv3imWdaCuAdIUm9hJNTe/HMMy3F8OEnRMuWQixdqt1xg8+eFYsKFRLD3dxEr0KFRJfChUXPGjXE5R49hHjxRSGOHtXuYDnA/v1CtGwpxKef3hAdOnwgypV7S5Qp01+UKDFAFCnyuZg7N8zRQ7ScoCAhXnhBhPbvL3rUrCm6FC4sehUqJN4zGET/4sU1Cw/ftUuINm2EGDFCiCNHjgtf3/bC2flVAfQRzs6vCF/f9mLKlL9Fy5Z
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f7024d3c2d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(-0.1,1.1)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using JuliaFEM.Core: Element, MortarProblem, calculate_normal_tangential_coordinates!\n",
    "\n",
    "using PyPlot\n",
    "\n",
    "#mortar_surface = :LOWER_TO_UPPER\n",
    "#slave_surface = :UPPER_TO_LOWER\n",
    "mortar_surface = :UPPER_TO_LOWER\n",
    "slave_surface = :LOWER_TO_UPPER\n",
    "\n",
    "fig1 = figure(figsize=(5, 5))\n",
    "master_elements = JuliaFEM.Core.Element[]\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    (eltype == :TR3) && (elset == mortar_surface) || continue\n",
    "    element = Tri3(elcon)\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    X = element(\"geometry\", 0.0)\n",
    "    x = [X[mod(i, 3)+1][1] for i=1:4]\n",
    "    y = [X[mod(i, 3)+1][2] for i=1:4]\n",
    "    plot(x, y, \"-bo\", alpha=0.5)\n",
    "    push!(master_elements, element)\n",
    "end\n",
    "\n",
    "contact_problem = MortarProblem(\"displacement\", 3)\n",
    "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n",
    "    (eltype == :TR3) && (elset == slave_surface) || continue\n",
    "    element = Tri3(reverse(elcon))\n",
    "    update(element, \"geometry\", mesh[\"nodes\"])\n",
    "    element[\"master elements\"] = master_elements\n",
    "    X = element(\"geometry\", 0.0)\n",
    "    x = [X[mod(i, 3)+1][1] for i=1:4]\n",
    "    y = [X[mod(i, 3)+1][2] for i=1:4]\n",
    "    plot(x, y, \"-ro\", alpha=0.5)\n",
    "    calculate_normal_tangential_coordinates!(element, 0.0)\n",
    "    push!(contact_problem, element)\n",
    "end\n",
    "xlim(-0.1, 1.1)\n",
    "ylim(-0.1, 1.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: # of field problems: 1\n",
      "INFO: # of boundary problems: 2\n",
      "INFO: Starting iteration 1\n",
      "INFO: Assembling field problems...\n",
      "INFO: Assembling body 1...\n",
      "INFO: Assembly: 10.0 % done. \n",
      "INFO: Assembly: 20.0 % done. \n",
      "INFO: Assembly: 30.0 % done. \n",
      "INFO: Assembly: 40.0 % done. \n",
      "INFO: Assembly: 50.0 % done. \n",
      "INFO: Assembly: 60.0 % done. \n",
      "INFO: Assembly: 70.0 % done. \n",
      "INFO: Assembly: 80.0 % done. \n",
      "INFO: Assembly: 90.0 % done. \n",
      "INFO: Assembly: 100.0 % done. \n",
      "INFO: dim = 807\n",
      "INFO: Assembling boundary problems...\n",
      "INFO: Assembling boundary 1...\n",
      "INFO: Assembly: 10.0 % done. \n",
      "INFO: Assembly: 20.0 % done. \n",
      "INFO: Assembly: 30.0 % done. \n",
      "INFO: Assembly: 40.0 % done. \n",
      "INFO: Assembly: 50.0 % done. \n",
      "INFO: Assembly: 60.0 % done. \n",
      "INFO: Assembly: 70.0 % done. \n",
      "INFO: Assembly: 80.0 % done. \n",
      "INFO: Assembly: 90.0 % done. \n",
      "INFO: Assembly: 100.0 % done. \n",
      "INFO: Assembling boundary 2...\n",
      "INFO: dumping matrices to disk, file = matrices_tie_contact_3d_host_1_iteration_1.jld\n",
      "INFO: Solving system\n",
      "INFO: UMFPACK: solved in 0.23568201065063477 seconds. norm = 1.2500614850006646\n",
      "INFO: timing info for iteration:\n",
      "INFO: boundary assembly       : 3.2558910846710205\n",
      "INFO: field assembly          : 1.247041940689087\n",
      "INFO: dump matrices to disk   : 1.1785790920257568\n",
      "INFO: solve problem           : 0.3804500102996826\n",
      "INFO: update element data     : 0.03557395935058594\n",
      "INFO: non-linear iteration    : 6.097559213638306\n",
      "INFO: solver finished in 6.260515928268433 seconds.\n",
      "INFO: nid near corner = 120\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n",
      "INFO: displacement X = [1.0,1.0,1.0], u = [0.027777777777777707,0.027777777777777592,-0.11111111111111081]\n"
     ]
    }
   ],
   "source": [
    "using JuliaFEM.Core: DirectSolver\n",
    "solver = DirectSolver()\n",
    "solver.name = \"tie_contact_3d\"\n",
    "solver.method = :UMFPACK\n",
    "solver.nonlinear_problem = false\n",
    "solver.max_iterations = 1\n",
    "solver.dump_matrices = true\n",
    "push!(solver, field_problem)\n",
    "push!(solver, boundary_problem)\n",
    "push!(solver, contact_problem)\n",
    "call(solver, 0.0)\n",
    "\n",
    "using JuliaFEM.Test\n",
    "\n",
    "nid = 0\n",
    "for (nid, coords) in mesh[\"nodes\"]\n",
    "    if isapprox(coords, [1.0, 1.0, 1.0])\n",
    "        info(\"nid near corner = $nid\")\n",
    "        break\n",
    "    end\n",
    "end\n",
    "nid\n",
    "\n",
    "known_value = [1/36, 1/36, -1/9]\n",
    "\n",
    "for element in field_problem.elements\n",
    "    i = indexin([nid], get_connectivity(element))[1]\n",
    "    i != 0 || continue\n",
    "    X = element(\"geometry\", 0.0)\n",
    "    u = element(\"displacement\", 0.0)\n",
    "    info(\"displacement X = $(X[i]), u = $(u[i])\")\n",
    "    @test isapprox(u[i], known_value, atol=1.0e-5)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: XDFM: ndim = 807\n"
     ]
    }
   ],
   "source": [
    "xdoc, xmodel = JuliaFEM.Postprocess.xdmf_new_model()\n",
    "coll = JuliaFEM.Postprocess.xdmf_new_temporal_collection(xmodel)\n",
    "grid = JuliaFEM.Postprocess.xdmf_new_grid(coll; time=0.0)\n",
    "\n",
    "Xg = Dict{Int64, Vector{Float64}}()\n",
    "ug = Dict{Int64, Vector{Float64}}()\n",
    "for element in field_problem.elements\n",
    "    conn = get_connectivity(element)\n",
    "    X = element(\"geometry\", 0.0)\n",
    "    u = element(\"displacement\", 0.0)\n",
    "    for (i, c) in enumerate(conn)\n",
    "        Xg[c] = X[i]\n",
    "        ug[c] = u[i]\n",
    "    end\n",
    "end\n",
    "perm = sort(collect(keys(Xg)))\n",
    "nodes = Vector{Float64}[Xg[i] for i in perm]\n",
    "disp = Vector{Float64}[ug[i] for i in perm]\n",
    "elements = []\n",
    "for el in field_problem.elements\n",
    "    isa(el, JuliaFEM.Core.Element{JuliaFEM.Core.Tet4}) || continue\n",
    "    push!(elements, (:Tet4, get_connectivity(el)))\n",
    "end\n",
    "#elements\n",
    "JuliaFEM.Postprocess.xdmf_new_mesh!(grid, nodes, elements)\n",
    "JuliaFEM.Postprocess.xdmf_new_nodal_field!(grid, \"displacement\", disp)\n",
    "JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/blocks_tet4.xmf\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://results.juliafem.org/divided-block-tet4/2015-12-17-tet4-results.png\"\n",
    "         style=\"width:80%; height:auto; vertical-align:middle;\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}