JuliaFEM.jl/docs/tutorials/2015-12-10-3d-mortar-assembly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3d mortar assembly\n",
    "\n",
    "Author: Jukka Aho\n",
    "\n",
    "**Abstract**: Mortar assembly for 3d problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using JuliaFEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{JuliaFEM.Core.Element{E<:JuliaFEM.Core.AbstractElement},1}:\n",
       " JuliaFEM.Core.Element{JuliaFEM.Core.Tri3}([4,5,6],Dict{ASCIIString,JuliaFEM.Core.Field{A<:Union{JuliaFEM.Core.Continuous,JuliaFEM.Core.Discrete},B<:Union{JuliaFEM.Core.Constant,JuliaFEM.Core.Variable},C<:Union{JuliaFEM.Core.TimeInvariant,JuliaFEM.Core.TimeVariant}}}(\"geometry\"=>JuliaFEM.Core.Field{JuliaFEM.Core.Discrete,JuliaFEM.Core.Variable,JuliaFEM.Core.TimeInvariant}([[-1.0,1.0,0.1],[2.0,-0.5,0.1],[2.0,2.0,0.1]])))"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes = Vector{Float64}[\n",
    "    [0.0, 0.0, 0.0],\n",
    "    [3.0, 0.0, 0.0],\n",
    "    [0.0, 3.0, 0.0],\n",
    "    [-1.0,  1.0, 0.1],\n",
    "    [ 2.0, -0.5, 0.1],\n",
    "    [ 2.0,  2.0, 0.1]]\n",
    "\n",
    "mel = JuliaFEM.Core.Tri3([4, 5, 6])\n",
    "mel[\"geometry\"] = Vector{Float64}[nodes[4], nodes[5], nodes[6]]\n",
    "\n",
    "sel = JuliaFEM.Core.Tri3([1, 2, 3])\n",
    "sel[\"geometry\"] = Vector{Float64}[nodes[1], nodes[2], nodes[3]]\n",
    "R = [\n",
    "    0.0 1.0 0.0\n",
    "    0.0 0.0 1.0\n",
    "    1.0 0.0 0.0]\n",
    "sel[\"nodal ntsys\"] = Matrix{Float64}[R, R, R]\n",
    "sel[\"master elements\"] = JuliaFEM.Core.Element[mel]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(\n",
       "2x3 Array{Float64,2}:\n",
       " -1.0   2.0  -1.0\n",
       " -1.0  -1.0   2.0,\n",
       "\n",
       "2x3 Array{Float64,2}:\n",
       " -2.0   1.0  1.0\n",
       "  0.0  -1.5  1.0,\n",
       "\n",
       "2x6 Array{Float64,2}:\n",
       " -1.0   0.0   1.0  1.0  0.25  -1.0     \n",
       " -0.5  -1.0  -1.0  0.0  0.75   0.333333,\n",
       "\n",
       "3x3 Array{Int64,2}:\n",
       " 1  0  1\n",
       " 1  1  0\n",
       " 0  1  1,\n",
       "\n",
       "[0.039743589743589755,-0.1952991452991453])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "time = 0.0\n",
    "x0, Q = JuliaFEM.Core.create_auxiliary_plane(sel, time)\n",
    "S = Vector{Float64}[]\n",
    "for p in sel(\"geometry\", time)\n",
    "    push!(S, JuliaFEM.Core.project_point_to_auxiliary_plane(p, x0, Q))\n",
    "end\n",
    "\n",
    "M = Vector{Float64}[]\n",
    "for p in mel(\"geometry\", time)\n",
    "    push!(M, JuliaFEM.Core.project_point_to_auxiliary_plane(p, x0, Q))\n",
    "end\n",
    "\n",
    "#S = Matrix[x' for x in S]\n",
    "#S, M\n",
    "S = reshape([S...;], 2, 3)\n",
    "M = reshape([M...;], 2, 3)\n",
    "P, neighbours = JuliaFEM.Core.clip_polygon(S, M)\n",
    "C = JuliaFEM.Core.calculate_polygon_centerpoint(P)\n",
    "S, M, P, neighbours, C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: npts = 6\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAG7CAYAAABQNZVMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcVPX+x/HXYdhxFxAVcWOHUbFFzbzVLcvUXHJNzQXN0jLtV9piZXVLu5ZLZYuJmmaWWbm32nZzaRWVHcQVF0AFVPZhzu+P464g4DBnhvk8H4/7oJkDM2+6xJvvOef7/SqqqqoIIYQQDspJ7wBCCCGEnqQIhRBCODQpQiGEEA5NilAIIYRDkyIUQgjh0KQIhRBCODQpQiGEEA5NilAIIYRDkyIUQgjh0KQIhRBCODQpQiGEEA5NilAIIYRDkyIUQgjh0KQIhRBCODQpQiGEEA7NWe8AovbKzs7m1VenER//J4piQlWdiYy8meefn42Pj4/e8YQQAgBFNuYVNSErK4t+/W5hxIh0wsJAUcBshuRkWLGiLevWbZcyFELYBClCUSMmTx5DSMhHhIdfeSwhAVJTR/PWW0utH0wIIS4j1whFjYiP/5OwsKsfCwvTjgshhC2QIhQ1QlFMKMrVjzk5aceFEMIWSBGKGqGqzpR30t1s1o4LIYQtkCIUNSIy8maSkq5+LCkJgoP9rRtICCHKIUUoasTzz89mxYq2JCRoI0DQPiYkwOLF0LPndvLzy2lKIYSwIrlrVNSYc/MI//xzBU5OJsxmZ4KC/Lj//gwaNAB399Z07PgHrq4yjUIIoR8pQlHjtm3zp6TkMG5u/tx0UxI7d3bjzJmdANSrdwvt2/+IweCuc0ohhKOSU6PCqpyd6xAZuQFX12YAnDq1jZSUaOTvMSGEXqQIhdW5u/tjNG7AyckTgKysT9m//yV9QwkhHJYUodBF3bodCQ9fCWiTDQ8ceIVjxz7WN5QQwiFJEQrdeHv3pW3bOecfp6SMJTf3fzomEkI4IilCoSt//yk0a/YIAKpaSnx8fwoK0nROJYRwJFKEQleKohAY+A4NG94NgMl0kri4XpSWntQ5mRDCUUgRCt05OTkTEfE5np4RABQWphEffz9mc4nOyYQQjkCKUNgEZ+f6tGu3CRcXXwDy8n4lJWW8TKsQQtQ4KUJhM9zdW2I0rsfJSZtcn5m5jIMHZ+mcSghR20kRCptSr14nQkMvTKPYt286WVmrdEwkhKjtpAiFzfH1HUjr1hdGgklJo8jL265jIiFEbSZFKGxSQMDT+PlFA6CqxcTH96WwcJ/OqYQQtZEUobBJiqIQHPw+DRrcAUBpafbZaRW5OicTQtQ2UoTCZjk5uRIR8SUeHiEAFBQkkZg4CLO5VOdkQojaRIpQ2DQXl4a0a7cJZ+fGAOTkbCYt7VGZViGEsBgpQmHzPDzaEhm5FkVxBeDo0UUcOjTnGl8lhBCVI0Uo7EKDBrcSGrr0/OO9e6eRnb1Gx0RCiNpCilDYjSZNhtGq1UtnH6kkJQ3n1Km/9YwkhKgFpAiFXWnZ8kWaNBkBgNlcSHz8fRQVHdI5lRDCnkkRCruiKAohITHUr38rACUlx4iL643JdFrnZEIIeyVFKOyOk5MbERFrcHdvC0B+/m4SE4dgNpt0TiaEsEdShMIuubp6n51W0RCAkye/IT39CZ1TCSHskRShsFueniFERHyFojgDcPjwAjIy3tE5lRDC3kgRCrvWsOHtBAcvOv94z54pnDixScdEQgh7I0Uo7F7TpqMJCHju7CMzCQlDOHNml66ZhBD2Q4pQ1AqtW/8HH59BAJjN+cTF9aa4+IjOqYQQ9kCKUNQKiuJEaOgy6tbtBEBxcQZxcfdRVpavczIhhK2TIhS1hsHggdG4Dnf3VgCcObODxMThqGqZvsGEEDZNilDUKq6uTTAaN2Iw1APgxIl1pKc/rXMqIYQtkyIUtY6XVwQREV8ABgAyMuZw5MhCfUMJIWyWFKGolRo16k5w8HvnH6emPsrJk9/rmEgIYaukCEWt1azZePz9nzz7qIyEhEHk5yfomkkIYXukCEWt1rbtf/H27gdAWdkpdu/uRUlJps6phBC2RIpQ1GqKYiAsbAV16nQEoLj4AHFxfSkrK9Q5mRDCVkgRilrPYPDCaNyAm5s/AKdP/0Fy8ihU1axzMiGELZAiFA7Bza3Z2WkVdQDIzl7Nvn0v6JxKCGELpAiFw6hTpz3h4Z9x7sf+4MGZHD26VN9QQgjdSREKh9K4cS8CA+eff5yaOp6cnJ91TCSE0JsUoXA4/v6TaN78MQBU1URCwgAKClJ0TiWE0IsUoXBIbdvOo1GjngCYTDlnp1Uc1zmVEEIPUoTCITk5ORMe/hleXu0AKCpKJyGhP2Zzsc7JhBDWJkUoHJazc12Mxo24uvoBkJe3heTksaiqqnMyIYQ1SREKh+bu3oLIyA04OXkAkJX1CQcOvKJzKiGENUkRCodXr96NhIWtBBQA9u9/iczMT/QNJYSwGilCIQAfn360aTP7/OPk5Ghyc7fomEgIYS1ShEKc1aLFkzRtOh4AVS0hPr4fhYXpOqcSQtQ0KUIhzlIUhaCgBTRs2B0Ak+kEu3f3orQ0R+dkQoiaJEUoxEWcnFwID/8cT89wAAoLU0hIGIDZXKJzMiFETZEiFOIyLi4NMBo34uLiA0Bu7s+kpk6QaRVC1FJShEJchYdHayIj16MobgAcO7aEgwf/q3MqIURNkCIUohz163cmLGz5+cf79j1LVtYXOiYSQtQEKUIhKuDrO5jWrV87/zg5+UFOnfpDx0RCCEuTIhTiGgICnqVJk1EAmM1FxMX1obBwv76hhBAWI0UoxDUoikJIyIfUr38bAKWlWcTF9cZkytM5mRDCEqQIhagEJydXIiO/wsMjCICCggQSEgZjNpt0TiaEuF5ShEJUkotLI4zGTTg7NwIgJ+d79uyZJNMqhLBzUoRCVIGnZxCRkWtRFBcAjhz5gIyM+TqnEkJcDylCIaqoQYNuhIQsOf84Pf1Jjh9fp2MiIcT1kCIUohr8/EbQsuWLZx+pJCYO4/TpHbpmEkJUjxShENXUqtVL+Po+AIDZXEBc3H0UFWXonEoIUVVShEJUkzatYgn16t0CQEnJkbPTKs7onEwIURVShEJcB4PBncjItbi7twEgP38XSUkPoKplOicTQlSWFKEQ18nV1QejcRMGQ30ATpzYyJ49T+qcSghRWVKEQliAl1cokZFfoSjOABw+/BaHD7+rcyohRGVIEQphIQ0b/pvg4A/OP05Le5wTJ77WMZEQojKkCIWwoKZNx9KixdNnH5lJTBzCmTO7dc0khKiYFKEQFtamzUy8vQcAUFZ2hri43hQXH9U5lRCiPFKEQliYojgRFracunVvBqC4+BDx8X0oK8vXOZkQ4mqkCIWoAQaDJ5GR63BzCwDg9Om/SUp6EFU165xMCHE5KUIhaoibm9/ZaRV1ATh+fA179z6rcyohxOWkCIWoQXXqRBIRsRowAHDo0GyOHFmkbyghxCWkCIWoYY0a3UNQ0DvnH6elTeTkyc06JhJCXEyKUAgraN58Av7+UwBQVRMJCQPJz0/UOZUQAqQIhbCatm3fpHHj+wAoK8sjLq43JSVZOqcSQkgRCmElimIgLGwldepEAVBUtI/4+H6UlRXpnEwIxyZFKIQVOTvXwWjcgKtrMwBOndpOcvJomVYhhI6kCIWwMje35hiNG3Fy8gIgO3sV+/fP0DmVEI5LilAIHdStG0V4+KeAAsCBA69y7NhyfUMJ4aCkCIXQibf3fbRtO/f845SUceTm/qpjIiEckxShEDry959Ms2YTAVDVUuLj+1NQkKpzKiEcixShEDpSFIXAwLdo1KgHACZTDnFxvSgtPaFzMiEchxShEDpzcnImPHwVXl6RABQW7iE+vj9mc7HOyYRwDFKEQtgAZ+d6GI0bcXFpAkBe3m+kpDyEqqo6JxOi9pMiFMJGuLu3xGjcgJOTBwCZmR9z4MBrOqcSovaTIhTChtSrdxNhYR+ff7x//wtkZn6qYyIhaj8pQiFsjI/PANq0+e/5x8nJY8jL26ZjIiFqNylCIWxQixZT8fMbC4CqFhMf35fCwr06pxKidpIiFMIGKYpCcPD7NGjwbwBKS4+
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f170b2a8490>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(-2.1,2.1,-1.6,2.1)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyPlot\n",
    "fig = figure(figsize=(5.0, 5.0))\n",
    "plot(S[1, [1, 2, 3, 1]]', S[2, [1, 2, 3, 1]]', \"-yo\", lw=2.0)\n",
    "plot(M[1, [1, 2, 3, 1]]', M[2, [1, 2, 3, 1]]', \"-ro\", lw=2.0)\n",
    "#plot(P[1, [1, 2, 3, 4, 5, 6, 1]]', P[2, [1, 2, 3, 4, 5, 6, 1]]', \"-k\", lw=2.5)\n",
    "npts = size(P, 2)\n",
    "info(\"npts = $npts\")\n",
    "# loop over cells\n",
    "ips = JuliaFEM.Core.get_integration_points(JuliaFEM.Core.Tri3, Val{5})\n",
    "#info(\"ips = $ips\")\n",
    "SM = zeros(3, 3)\n",
    "SS = zeros(3, 3)\n",
    "for i=1:npts\n",
    "    xvec = [C[1], P[1, i], P[1, mod(i, npts)+1]]\n",
    "    yvec = [C[2], P[2, i], P[2, mod(i, npts)+1]]\n",
    "    plot(xvec, yvec, \"--ko\", lw=1.0)\n",
    "    X = hcat(xvec, yvec)'\n",
    "    geom = JuliaFEM.Core.Field(Vector{Float64}[X[:,i] for i=1:size(X,2)])\n",
    "    #info(\"geom = $geom\")\n",
    "    for ip in ips\n",
    "        dN = JuliaFEM.Core.get_dbasis(JuliaFEM.Core.Tri3, ip.xi)\n",
    "        J = sum([kron(dN[:,j], geom[j]') for j=1:length(geom)])\n",
    "        w = ip.weight*det(J)\n",
    "        x = vec(JuliaFEM.Core.get_basis(JuliaFEM.Core.Tri3, ip.xi)*geom)\n",
    "        plot(x[1], x[2], \".g\")\n",
    "        theta1 = JuliaFEM.Core.project_point_from_plane_to_surface(x, x0, Q, sel, time)\n",
    "        N1 = sel(theta1[2:3], time)\n",
    "        theta2 = JuliaFEM.Core.project_point_from_plane_to_surface(x, x0, Q, mel, time)\n",
    "        N2 = mel(theta2[2:3], time)\n",
    "        SS += w*N1'*N1\n",
    "        SM += w*N1'*N2\n",
    "    end\n",
    "end\n",
    "plot(C[1], C[2], \"ko\")\n",
    "xlim(-2.1, 2.1)\n",
    "ylim(-1.6, 2.1)\n",
    "axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: SS = \n",
      "[0.5123885459533705 0.29124871399177565 0.23957261659808454\n",
      " 0.29124871399177565 0.40967399691358797 0.23773469650206241\n",
      " 0.23957261659808454 0.23773469650206241 0.2491587362825836]\n",
      "INFO: SM = \n",
      "[0.4042245370370448 0.38539094650206535 0.2535943930041208\n",
      " 0.21609760802469608 0.37920524691358753 0.3433545524691425\n",
      " 0.24657600308642466 0.1835519547325143 0.2963380915637917]\n"
     ]
    }
   ],
   "source": [
    "info(\"SS = \\n$SS\")\n",
    "info(\"SM = \\n$SM\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.1",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}