diff --git a/notebooks/2015-12-18-normal-tangential-coordinates.ipynb b/notebooks/2015-12-18-normal-tangential-coordinates.ipynb
new file mode 100644
index 0000000..63c8c63
--- /dev/null
+++ b/notebooks/2015-12-18-normal-tangential-coordinates.ipynb
@@ -0,0 +1,2973 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Normal-tangential coordinates\n",
+ "\n",
+ "Author: Jukka Aho\n",
+ "\n",
+ "**Abstract**: Use of normal tangential coordinate system"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Array{Float64,1}"
+ ]
+ },
+ "execution_count": 1,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using JuliaFEM.Core: Quad4, Seg2, Seg3, PlaneStressLinearElasticityProblem, DirichletProblem, update!, DirectSolver\n",
+ "typealias Node Vector{Float64}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "using Gadfly\n",
+ "set_default_plot_size(10cm, 10cm)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Element normal definition in 2d"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAF6CAYAAAAXoJOQAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeXhU1fkH8O9770wWICKiKAhkEhBwgZnJgLZ1Cy6tCFRtjVWsFq0/la2ta1tbK2oXsS6tgvtW94JWWwQ3ZIJW65KZuTeIRZDMJCAIioiBJCQz5/39MQlFBR0gOSfhvp/n4SnJ3Lnnmy+Z0+uZO/cCQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCdEnMTFVVqb6mc3hJPF7Xz3QGL5G+9dLZt6VroK4uFov5fD48YjqHlxDxE6YzeIn0rZfOvmWiz1F9fYSJEDedw0uIOGY6g5dI33pJ30IIIdoNmQ6wu6b+ZOrtID4FwIDmvJb97rnnnk+2t93kiZNDFughBnoSaBkrnjDr4Vnrcx2Hma1EInV0WVnJonYLL75WIpEsD4dLKk3n8ArpWy+dfXf5pRsimu1LtxwBYPPXbWeB7meFa2c9NKsEpByy6LqdGScWi9mWRVfvVlixk6xrTSfwFulbL319d/mJ/vaHbn/tL4/es+brtpl0/qRSAP1nPTzrWQBA2roL4IqdGaemJqKI8OSuJxU7i0g9bjqDl0jfeunsu8tP9LnwMQYAWAWAAWBt09o6AL2nTZuWn+s+zjiDMsFg4L4Oiii2IxQqudt0Bi+RvvXS2bdP10BdQSKx4iDLsu9g5jdCoZJrXLfu+4CaBvD9778fmDNkSPKfRFZ+JqMuLysrdROJ5N8si/q1tHxyMtCrm99vP8VMH4RCxZNcN/ktgK4HaG4wWHyb46TOJ8JZStGN4XDxy4lE6s+WhRDgOy8Y7L/KcVLPAdgSCgV+WF1dW8rMdwN4KxgM/DaRSI21LPxCKTwUDgceSyRSV1gWvgvgl8FgIO44qQeIMGCvvXh8Q0N3Xzrd8AyAZDAYuDCRqBtlWeqPzDw/FCq5NZFI/sSy6McA3xwMlrzgOMkbiCiSyagLyspKaxOJ5D8ty+JgsPjUeLym2Lat+5ipKhQq/rXrJk8C6DKAHgkGix92nNSlRBijlHVVODzwnUQida9lIeDzdTtt48Ytme7dM/9iRl0oFPip66bKAMwA8GIwGLjJcWp/TMQ/UQp/CYcD8xKJ1B8sC4cT0UUjRhTXOE7qH0QYFgwGDnHdVf2B9IMAEsFg4MpEovZEy+IrleLHwuGSh1y39ucAjwP46mCw5E3Hqb2LiAe1tGRO79lTNTU0+J8D6MNgsHhiPF4TtG3rJmZaEAoVz3Dd5FkAnQ/QbcFg8VzHSV1HhG8rlZkcDg9ankik5gDoHg4HTo7H6/rZtvobM6pDocBlrps6DsCvATwZDAbud5zUVCKcohSuCYcDb7huahaAIT6f/aN16/p/vs8+tc8z46NQKHBOdfXK4cyZWwBEg8HAHx2n9gwi/j+leFY4XPKs69ZeA/BRgJoaDJa+7zi1TwLoFQoVf2/x4pr9lbIeVYqXhMMlv0gkkuWWRb9hpjmhUPE9rpucBNAPAOu6YHDga46Tuo0IB7e05J0VifRdX11d+5JS+DgcDkyIx2sPsW3+KzNeJUIjgCSAiwHcFQwGnnac1G+JcGwmQz8vKyt+z3VTjwHoM2JE8XdjsTW9/f7mJ5ixNBQKTHPduqMB9Tul8Ew4HLjDcWovJOIKpfgP4XBJpevW3grwYZalfjx8eOlax6l9EVCfhUIlP3LdmqGANVMpfj0cLpmeSCRPtSyawkz3hkLFs103+WuAjiOyLx0xYsBix0k9TIS+n35aPKZ377oiZp6tFJaHw4HJiUTqO5aFa5npX6FQ8e2um/opgDMBuiEYLH7FcVI3ESGYyVg/KSsbuDqRSM0nQkMoFDg9Fls52OfL3MmMN0OhwNWuWzse4J8B/EAwWPJEIlF7pWXxiUrxFeFwieM4tQ8Scf9u3VrGrV+P/MJC/9PMXBMKlVzkuskjAPq9UjwvHC75SyKRnGhZdDaz+nMoVPqS4ySjRJS2rPT5w4cPXuk4qbkA0qFQ4DTXTZU0Nfk+P+KI/jm/j+gJUydO2XThhRfuu73HJp0/qXTqxCnr0Prm89Rzppa0fp2zqqoqv+OkFrRDVJGjRKJW3vjWSPrWS2ffe/TSzeSJk0+acu6U3nc+cGcNiFdNmzj5NABgmycR46md2VckEskoxdd3TFKxfeoa0wm8RfrWS1/fXf/0yvOm3A3GWAD9APoIwOKZD838HgBMnTjlA2XxxDseuOPf086dFmZL/Q2MXiD8l8g6+/YHb//YbHohhBCdRjTKPsepvdN0Di9JJFL3ms7gJdK3Xjr73qOXbtpTUVGMAD7IdA5voSGmE3iL9K2X9N0pLV26tMh0Bi+RvvWSvvWSvoUQQgjd5PRK/eR0P72kb73k9MpOqL4+wgC+9lILon0R8WrTGbxE+tZL+hZCCCF0Y2Zq/Ri/0CQerxtpOoOXSN966exblm5yFIvFfMy40XQOLyHim01n8BLpWy+dfctEn6OamogC6FXTObyESFWazuAl0rde0rcQQgihGzNbrps803QOL0kkUhNMZ/AS6VsvnX3L0k2OYrGYzUwXmM7hLXSR6QTeIn3rpa9vmehzFIlEMsy41XQOL7Es3GQ6g5dI33pJ30IIIYRus2ez7Ti1vzOdw0scJ3mt6QxeIn3rpbNvWbrJUWlpzAL4GNM5vITZKjedwUukb7109i0TfY4ikUgayEw2ncNLLCtzoekMXiJ96yV9CyGEELpVVVX5Xbf2WdM5vMRxUs+ZzuAl0rdeOvuWpZudwMw9TGfwEmaSO/BoJH3rJX0LIYQQujEzVVWl+prO4SXxeF0/0xm8RPrWS2ffsnSTo1gs5vP58IjpHF5CxE+YzuAl0rdeOvv26Rqoq6uvj3CvXnXLTefwFl5mOoG37GTfb3AhGtAXFvoCOAAW+oGxP4C9AewNQq+tf2f0AlDQ+kw/gFze79oMoBmEDWA0t369CUALCJ+C8SkYG7b+3cIGKHwKG+uwGR/iZPp8p34e7fT9fpOugYQQXUwV+/EZSmCjFIxSEEqBrX+KkZ3EGcA6ENYCWAXGOjA2APgMwGewtvm7QgMAQKEF2Qn76/nRHRnkgdALjDxY6A6gBxj5rd/bB4R9Wv9PZJ9t/vRB9v9MGgGsBvARgDVgrIGFJBRqQFiBTajBeGpov8I6L5noc8TMViKROrqsrETbndu9LpFIlofDJZWmc+zxmAmVGHRy3qYz52/poUAYDuBQEIaCYQFYCaAGjBpQ6/8yUgA+RE+sw0hqMZp/exbw/vCjDxj9wegDxoEg9AW1/p8WUAIgH8BHW38m4H0wlsHGcjRhGb5Hmzsyos7fb1m6yVEsFrN9vn2vBiATvTbWtQCONZ1ijxPl/iCMgsIoEA5HJUYCsN/P5DMszIHCO2A8CGAJjqcPTcfdJSfQWgBrASze7uPT2cJoHIhM63+pEAaDEQThdCgMQR4KEeWVAJaBsQyEahBcMBZjNH3zf43kRN/vt0z0OaqpiaihQ2ufNJ3DS4jU46Yz7BEW8cFQOBbY+ucAMN4H4U0Ac5DBpfgUS54e8uEFoVDJ3WbDajKdFKZjJbL/tfLFgzdmwkIMhI0hYAwB4WAAE8C4AYwiRHkFABeAA0I1mhHDd2l1rkMzswXg8E8+qV/OzEcBeIuoY/+rSJZuhNjTvMp9oTAOjOORndj3A5BAdkKrhA+v42jaYDRjV8RMeAUl8CEIRhBAEEAIQADAhwDeAeFtMN5BGu/gRNr41V3woZvTuD9NKFnWACopAAptrOlu4QIieqejostEn6PZs9keOrT2vGAwcJ/pLF7hOMmLPHOEubsWchCE8QC+D2AkgHcBvAhGJZrwWi5noEjfu+hV3g9pjAJhFIBRAA4HsC+yyz5vA/g3gNeTo7HiAAX3vjUo+UMdfJvTQJ4FXDIAmUv7Y22hhcOI6LOOiCgTfY6qqqr8Pt++z4dCgRNMZ/GKRKJ2UThcLGv0O7KIw8jgLBAqQOgDxkIA85DGfJxIdTu7O+m7HUU5AMbhIHwbwJEAwn5C09E9UfB6PewtCgBn598CG+yMRMPgPFzg85EsD5vEzFY8npQXgUaJRLLcdIZOZyEPRZSvQZSXIsqbUMmPoZLHI8oF3/zkryd9d6AXufuMOr7pvKXcRFFmfOnPsx9zg1L8XkcNL0f0QnR2/+YiNOMsEC4AIQjgBTCeQDPmdvQpgKL9MPOpdVvw8KA30SPNX5x7nVHYPLwQT9s2/cRUPgEgGmWf49TeaTqHlyQSqXtNZzAqyt9ClO9DlOsR5TgqeRpe4H06ajjP993BmHmvRsXrfl3DGbuSVdvR/MXvM9dn+PNG5oEdNbacXpmjoqIYAfseZDqHt9AQ0wm0m8/5KMDZsPDz1g/5PAbG0RhNTscP7sG+NSKiz5n57MsG4KGz+qDXe5uRP7Q7mgbko7GQcb4NrDKd0fPk6pX6eepqii9yn9a1949QyS4q+TzM53ydETzVtyHM7Ktn7rOF+azPmlpuaGY+h5n7MXOHHnTLGr0QJkV5MIBfAZgAYCEYt+I4esVwKtEOTvvud/ukm/0PcPZ01zSIFlIPunju3Lnar68jlynOUVVVld9xUgtM5/CSRKJ2z73cRJSHIcqPIPsRfR8YYYymcSYn+T26bwO2NPssgO/ahM39M4VqEBT3xCb1q7bHdfYta/Q5yl6muHaN6RxeQsQ5f6y8y1jIh4LwWwCnAXgUGRyKE6jGdCxgD+3boPmV8z8C0HZf2PS40eMWKebhbY/r7FuWboTQIcoHgHEdCOeA8DAYf8JoSpmOJfQoLy8vKEL3t5npN88tem6u7vFl6SZHzEyOs0LOStDIdWuGms6w2+bzXojyDQCWwwLDwhCU00WdcZLfI/ruhMrLy3090P1xBby47SSvs2+Z6HMUi8V8gH2H6RxeopR9j+kMu4yZUMnnoBDvgnA0FL6HcroIx9JK09F2pEv33UlVVFTYPdDtUQbWzaucd8W2j+nsW9boc1RTE1FDhtS9ajqHlxCpStMZdskCHoFKzAJQDMYlGE1Pm46Uiy7bd+dFjes2P0BETc9Vzpv0lQc19i1r9EK0l+xVDG8G4XQANwK4AaOpyXQsYca448YdDsVvEbCGAQUABETnVs47x3Q2sQPMbLlu8kzTObwkkUhNMJ0hZ1H+EaK8GlF+Awt4hOk4u6JL9b0H0Nm3rNHnKBaL2cx0gekc3kIXmU7Qhpl7MfPRzDyWmYcxsw0AiPK+iPJsAPcCuA6LcBROoGqjYXdZ5+nbG/T1LWv0OYpEIpnq6pTcdEQrNn4TDGbu3gx8rz6NyzYpHLAhDaufHw15Fl495C1e+V4DrgAQBTAEo+kj03l3j/m+vUVf37JGL8TXaGQuTzXi3gdXo/jp9bCbFVCSD04DmfgmkI/wi01H4U4QsemsQuyILN3kaPZsth2n9nemc3iJ4ySv1TneuGPHvTKufOyqcceO/XBc+diXTy0/KdCiMMr9HP3vWgPfigZYK5tgvbYR9vsN8L84Am79UYjvKZO87r69TmffMtHnqLQ0ZgF8jOkcXsJslescj/x03nOV8/o/t2hefxC/k4Z1w2ct6FvbDHye/t9//TKAQAFUAaEYwJU6M3Yk3X17nc6+ZaLPUSQSSQOZyaZzeIllZS7UOd7cBXPrAGD69OkEWD4iWtfdxqZ8+uoSZ74F7uVHPYAtOjN2JN19e53OvmWNXohtjDt23CsgPgJAitM4uvzR5y646UPMYAZ93JI9mu+fD/ysP9JT+uPlboQpRJQ0nVsI0Q6ytxJMPWI6h5c4TuoJE+NGIhH/+PKxNw255rF3KMqNZyzhl5Zu5vXPr+fGxz7iprc+5/rVTfxqC/OxzLzHHCyZ6turdPYtp1fmqPVWgnKHKY2Yycgdj2K/rKIl3TYMyNh5ZRZh3N8PwX8agWF9/ThCEYryCcuoBW/7gNW0h7wRC5jr26ukbyE0O+3403qffPzJQ/Aa90JULSx4bsPHR024RG40I/YI8mZsjuSesfrpvIdpmtLdG/fq93T3dR+u2ze5+NvH3Dnpte4r3/PUNUnknrF6Sd+dkNxKUD+tt7ar5G8jyp9gId+G2a2XN/AYuZWgXjr7liP6HNXXRxig5aZzeAsv2+VnAhUM7JfTxgt5LBgvArgex9HPcAZldnXcrm3X+xa7QvoWYpcxcDADDQwM/MaNF/IZiPLniPKPNUQTQnRmzGzF48ljTefwkkQiWb4rz2PgSQa++W5gC/l8RLkeC/nkXRlnT7OrfYtdI313QrJGr9+urGEyEGRgMwNf/0bXQv4dKnk9KnnULgfcw8gavV6yRt8J1dREFBGeNJ3DS4jU47vwtN8BuI+A1TvcIso3gDAVhBNQTu/scsA9zC72LXaR9C3ELmBgVOvR/AHbe3z69OnWAXe8/goqeS0WclB3PiFEJzd7Ntuum5I7TGnkOMmdugMPA/9i4OYdbrCg5dr85zYoLOAhux1uD7SzfYvdo7NvWbrJUWlpzGKG3DNWI2Yr53tqMnAEgHIAf9ruBlH+pa0yVx7+yNU87vfj7hh/7Nin2yflnmNn+ha7T/ruhOSsG/125qwEBp5n4IbtPriQJ6OSN/qeb/j2uPKxm9or355GzgLRS/oWYicwcDQDGxno/ZUHo3wmKnkzonxCeXm5TyZ64UWydJOjaJR9rpu60XQOL3Gc1C05bnoNgNsIWP+F71byaAD3ATgdo0lOjf0GO9G3aAc6+5aJPkdFRTFiRpnpHF7CTJFv3Ca7Ll+GL78J+wqXAXgWwE9RTs8DQGVlZRqMTEVFRV77p+36culbtB+dfctEn6NIJJJOp+Gpqxmaxkxn5bDZ75E9mv9s63eiPBgWXoDCHzCa/v6FrS2e1fhxQ3xc+Vg5wv+SHPsW7UT6FiIHDJzAwHoGem79ZpT3RpQXYyHfYzCaEJ2KHNHnSC6BoF8OHxG/DsAtBGwEAMxmG4S/A/gQBLmR+06SSyDopbNvuZXgTiAiOWNDIyKu39FjDIwBMAjAX7d+sw9uQ/ZTsUdiNKU7POAe5uv6Fu1P+hbiazBADFQxcOXWb0b5YlTyWrzCxQajCSG6MmYmx1khH53XyHVrhm7v+wyMZ+AjBroDyJ5Gmb3c8JFaA+5hdtS36Bg6+5Y1+hzFYjEfYH/zNc5Fu1HK/sobqgwQgGsB3EjAZrzCB4LxBIArcBy9rj3kHmR7fYuOo7NvmehzVFMTUQC9ajqHlxCpyu18+zQAfQDcgSr2w8LfwZiP0XSX3nR7nh30LTqI9C3EdjBgMfAuA9MAAJX8V0S5CvM533A0IcSegJmtRCI11nQOL3Hd2vHbfs3AjxioYyAflVyBKH8ib762ny/3LTqWzr5l6SZHsVjMJsIlpnN4iVK4vO3vnD0V+DoAf6Qol4BxPwjn4XiqNZdwz7Jt36Lj6exbJvocRSKRDBHfZzqHt/Dd23xxJgD/Fef8+REAjwJ4GOU010yuPdUX+hYdTl/fpGugjjJ54uSQBXqIgZ4EWsaKJ8x6eNb6L283ZeKUpZT9qHwm+x368cyHZlbqTbv7Wj/uPwHAAwRsMZ1Hh9aj+f8CmEGVfDAYx6MRR+Bk8sTPL8Tu6vJH9Bbofla4dtZDs0pAyiGLrtvRtuSzRs18aFb/7J+dm+Rnz2bbcWp/t9uBd98gZD8otIKBnzFQaDpQR3Gc5LWtfz0HAPV+dv1KMC5CBhNkkm9/2/QtNNDZd5ee6CedP6kUQP9ZD896FgCQtu4CuKIjxiotjVkAH9MR+94ZBMSRnewnAZgIYCUD0xnY22iwDsBslTOQB+CalfsNuPXTnvvcD+C3OIHeM51tT8RslZvO4CU6++7SE72PMQDAKgAMAGub1tYB6D1t2rTtnm7HaVU5deKUFVMnTr7zyvPPL9qZsSKRSJpom4/cG0SAImAugAiA8wCMQ/ZslBsY2MdsuvbDTJch+39mjaWP15wAgovyba5tI9pVa99CE519e+aiZpShMTMfmZmadva0IuSpextUwY3IHhVvVVW1YqDPZ/2GmdxwOHBHIpEsJ8JZRHiGiF5wnFQoe+d2dXMoNGiZ46SuB7jPhg2BKX37vl/Y1JR/E0ArQ6HA7xcvrhmRydAUZloUDgceTySSpxJhDBEeDAZL3nSc1KUAD7Vt/t3w4aVrHSc1E0BzKBS41HVX9WduuRqgxaFQYGZ1dc0xStHZzPTPcDgwP5FITiTCt+MZ+kskEphb7SRH9vr7I5/0v+Ga48A85bPR3019NOWK54edfuKVrlt3GHNmGmD9OxQqfsR1677PnBnLbD0cDhe/7rq1P2dWh6TTNH3kyMAax0n+lRkcDpf8Ih6v62dZmWuIrPeCweK/JhK1RxKpc4ms54LB4rmOU3sOoI4CrNtCoeIlrlt7DbPql07nX9K79xa1cSP+yow14XDJ9FgsdbBt8y+Y8UY4XPK3RCI1loi/T2Q/GgwOfM1xaqcB6jDAd10oNOBDx0ndYjd+ng9g/MOnXfyyz2/9+IG9PrrhXOrL8Xjq25bFE5kxPxwu+WcikTqbiI+xLN/MESMGLHac1G8BHuD3d7+ssbFui8/XeyZA60KhwNWuWzOUmS4F8GYoVPJgdXVqjFJ8qlJ4vKysZFEikZpMxMF0Wv1h5MhBdY6TuokIhcFgYEo8vnw/y/L9HqBloVDgZsepORygnwL0QigUeMZ1k2cyYzQz7gyHSxzXTV3FzMVNTekr6uoGbx4yJHUHM30SDgd+E4utHGzb6SuI6J1gMHCf4yS/B+AHRPT3YDCw0HGSFwMIM+NP4XBJynGSMwBrr1CoeNJbb63qnZ/f8kdma0U4XHxjPF430rIy/0dkvRwMFj/lOLVnAOp4Iro7GAzEHSf5KwAllmX9avjwgZ+5buouABtCoZJfVVfXliqlfglYsVBo4D3xeO0JlqUqlLLmlJUVL3Cc2gsBFbFt+8bhwweucJzaPwFqn2AwcLHjpHoSYQaAZChUcoPrpsqY+SIAC0Ohkr+7bu3pzOpEZvu+cHjgO4lE6goiHrxli/+qI47ov95xau8k4vpgMHBlIpEMEOHXRHCCwZI7HSc1GuAzAfwjFCp50XVTP2Xmw5nVTeHwoOWJROoPRLzvsmWByQMHftC9oMD3Z2aqC4cDf4jHa4KWRZOJUBkMljzhOKnTAD6JCA8EgyVvOU7qMoCHEGWuDgYHr3Pd1CxmNIZCgcsXL/5gQCZj/5aZqsPhwKx4PHmsZWGCZdGzI0YEnnec5HkAvkXEtwSDpe87Tuo6gPdPp9dP9fv75zE33gLQqlAocH119crhSqWnMtOr4XDgsUQieQoRTlaKHiorC/yHKHOU4yT/z+/HNYceWvKR46RuI6JMMFh8ieOsPJA5vTkcLvkM7aBLT/Rpwko7e1RPAHj/gv0HArz+9ttv/8r67cxHZiYB4PbHbv/8Z+dOvost+sptvHr2VGsbG30zbFs1AMCWLel4t255dZmM+jQaZR9z7Rjbtn65aZNvDQDYdvoeZr+/vByZOXOGNgwbVjdDKWtL9rkFH+TltcxobvZ/DgB5eb5XMxlVbduF6wAgncYTeXlW4ZYtH38KAES4OZNRCgB8vo3rMpkeW3M0NGScbt3yVmUy6lMAYE7Ps+28VxsarDUAQHbmvg1n/eRv/f94dRJ+/wl7vfHajJ6VL09hIH/jBRP+Wnf/YzPS6eb67Dj8GpH1bl5e48cA0NLCs/PyrMLCwk3rs/vGrW195OfXf5LJ9JjR3MyNAFBY2Fjd3Fw4A8AGAFCq+QWfL+/15mb/R9ln2PdbFvIikb5NALB4cd0MopaWbB92qkcPtTWHbdMbAP23oWHLJwDg96s5mYw1z7Y//zibg2/b78F7/9Xi8286b9qsU4bY6qpD1JbHAaBbt6Z3m5sLZ9i29RkAZDJ5L+bltfynudn+KPsz+h4iyuQtWbJfY0XFfmrbHJs3++t69FAzmFs2ZZ/L/7Ft6/2mpuZPsv9O/I9MxnqhZ8/02uzPqG63bcvK/n3jZz7ffjPSaWrKZu7xXibTuDVHS0v+gry8lreZ89dmfwb7b5al8uvqBm+uqIBavNia0dycTgPAPvukV9XXW1tzMOMt27aWNzen12f74WeY7ZcKC5vXZMfmO3w+2NnfhwM3FhbWbc2Rn1+wNJNpnNHYaG9s7fqVwkJUteUAfI9Ylspfv35gPRFxdXXt1hxFRWp1fb01gyiz2XFST6TT6Yt9Pl9NOt2WI/NPZntBjx6Z1dnf2/SdeXk+HxFxNMqbevf+3+88UPC+ZW3ZmqOpyY4WFiJu2wXrsj+j9Zhtc0FDw4EbAcCy8Od0WmUAoEeP9JrGxrwZRJnN2ee2xLp1y6tt+50H0nMtyx8tKsLq7GskfRez319RATVnzuDN2772lCpY7vP977Xn99uLMhnltr32lLKe8Pm4oLl5/YbWX/WblMq+9goKeF1jo7X1tdfcnE5065a3si0HUWYekX/Rpk1W2xxwL7PfH4lE0nPmQA0bVjeD2W7OPtdekZentp0DXstk1OKCguxrD7DKLQuXNjZ+vD67b9yqlOJs5s8/XrLkkAxE1tTzJsenTZz8AwCYMnHKjdN+MmXr9WgmT5x80pRzp/S+8vzzi6ZdMK0/AEybNi1/ysQpD06ZOOXBnRmnq12PnoGjGFjAQBMDdzNwoOlMO4OBwpaevbb8fNpfXkcl/9N0Hi+Q69HrpbPvLn965bRzp4XZUn8DoxcI/yWyzr79wds/BoCpE6d8oCye6E/76zJWej6AXmAoJizKs/N+duv9t35qOP4uKy8vL+jB3dcy8BcAw4i4hNieMnfR3He23Y6BowD8EsCJAP4G4HrKvq/RqTHwiw1FvX7e+5/r9+Y0HYrv0mrTmYQQHrB06dKdegO3I5WXlxeMKx/LY48dezQAjBs97uRx5ePm72h7Br7DwFwGtjDwMAOD9aXdOQx0U0Rrf/S7v29AlMJEHBUAACAASURBVC8ynccrOtPvtxfo7LtLn3WjU1VVlb+pqeAZ0zm+gNEwb9G81wCAbHoXQOmONiXgDQLGAzgcQAGAJa0Tfme8xv6U1AEBrj7+VMIi3Gs6jFc0NhY+ZzqDl+jsWyb6HNXXRxig5aZzbIsILW1/T6fTisDf+OY6AS4BZwAY2fqtdxmYzUCnuOkEA0Vp2/ebS6bc2uvybp+9jOmkTGfyDl5mOoG3SN/iG5SXlxeMLx+79dSrMeVj+o8vH/vBzu6HgUNbj+ybWyf8g9s36c5psXy/dUtHbEKUrzeZQ4g9iRzR54iZrXg8eazpHO2NgCUEnAtgGLKnTTqta/llurMwsHfGtq/63XnX1SMff0wkkuW6M3iZ9K2X9N0JdbXTK3cVAwEG/spAY+uEP0rX2Gv2OeCmfx92ZBpRPgmQ0/10k7710tm3HNHnKHsrQZ5nOkdHIyBFwM+RPcKvAfAqAy8zcERHjsvAPj03b5x22w9/FsNoegEALIvlMsQaSd96Sd+i02Bg/9Zr6DQw8G8GRnfEOO/3H3LvohHHZLCID+qI/QshxDeaPZtt101dYDqHKQz0aZ3wN7dO+Me31743de++f0N+Yfqcqx5+bNvvZ68rJHSRvvXS2bcs3eSotDRmMeNM0zlMIWAdAb8CUAxgAYCnWyf83b7v5bsDD7vrnaGjMo+cfM7F236f2Zqwu/sWuZO+9ZK+O6E99aybXcVA79br4H/KwOsMjOdduKTG24eNGrA5v5ua8vOZN375MTkrQS/pWy/pW3QZDBQx8EsG1jOQYKBiZyb8l8tOqFw04pjPMZvtjswphBDfKBpln+umvnLUKbIY6MHAzxlYw4Cby4T/6JizD2nMK1CXT/rzL7b3uOOkvnIpadFxpG+9dPYta/Q5KiqKEbP+DxF1FQRsIuCvyF4s7QFkr6pZzcC5DGz3aL33+vV/e/vgwzfcdOcVf9ne48wU6bjE4sukb72k706ImamqKtXXdI6ugoF8Bi7k7D1tl7RO+FuvxXPTjy4b1ZhXwNf/+Lfn7Wgf8XhdPz1pBSB96yZ9iz0GA3mtk/wHDKxonfx980eNWRwNln9kOp8QQmzllUsgdBQGChmYxsDKFX1LaxrzCviWiku+9vQy+Ui+XtK3Xjr77tL3jNWNiDaZztBVEdAI4HYG7jnrmicW9P14zep/Xn3q41/7HOJ6TfEEpG/dpG+x54pyAFFuwkIOmo4ihFfIWTc5YmZynBWd8W5MXQvjKgDP4Thyv2lT163pFDdD8QrpWy+dfctEn6NYLOYD7DtM5+jSXuaBsPATAL/PZXOl7Hs6OJHYhvStl86+ZaLPUX19hIkQN52jS/PhUjBewWhyctmciGMdHUn8j/Stl/Qt9jwvch9U8mZE+VumowjhNXJEnyNmthKJ1FjTObosP34BxjsYTW/m+hTXrd3tK2OK3EnfeunsWyb6HMViMZsIl5jO0SXN571AmATCjJ15mlK4vKMiia+SvvXS2bdM9DmKRCIZIr7PdI4uqRD/B+B9lNPzO/dEvrtD8ogdkL71kr7FniLKPkQ5hSh79qYtQpgmR/Q5mj2bbcdJytLNziKcBgIBeGpnn+q6KVlK0Ej61ktn3zLR56i0NGYBJG/G7izGJWDcjtGU3tmnKkXy5qBG0rdeOvuWiT5HkUgkTYQrTefoUrKnUh6CRuzSB0OY6bJ2TiS+hvStl/Qt9gxRfhRRvt10DCGEyEk0yj7HST1iOkeXsYD3R5SbsIB3+fpAjpN6oj0jia8nfeuls29ZuslRUVGMAMgdpnJl4SIAUZxAy3Z1F8wkd+DRSPrWS/rupJYuXVpkOkOXUMV+VPKHWMjf353dSN96Sd966exbjuhF+9uEcWA04xPMMx1FCCETfc6qqqr8TU0Fz5jO0SUwJgG4D2dQZnd209hY+Fw7JRI5kL710tm3TPQ5qq+PMEDLTefo9LJvvh6DDNrhchG8y+v7YldI33pJ36KrquQZqOSnTccQQvyPHNHniJnJdVNlpnN0alXsBzARjHvbY3fxeN3I9tiPyI30rZfOvjt8ov/RmNPPPuPkioZTTz117+18v+m0007r3dEZ2kMsFvMx40bTOTq1z3EyGFuwCC+1x+6I+Ob22I/IjfStl86+O3yi/xybnwKwKa8579xtv6+ILgLRU88888z6js7QHmpqIgpgOYvk6xAmgvEoppNqj91ZFs9tj/2I3Ejfeunsm3QMcsaY028A0fjZ8+ccCgA/HP/DYXbG+i8rHD3nhTn/1pFBdLAoHwCgDhYOxbHyprUQnYmWNXqLfXcBGFZxUsVRAGCn7YsAvNuVJvnZs9l23dQFpnN0WoSzALzdnpO84yQvaq99iW8mfeuls28tE/2TLzyZAvA8WXzRmDFj8kF8LjO61N1VSktjFjPk5hk7wjgXQLteC4jZmtCe+xNfT/rWS2ff+s66UXQHQKcXUdFFAPLR0r6TQkeLRCIZZtxqOkentIBHgDAEaTzZnru1LNzUnvsTX0/61ktn31rW6AFg+vTp1n/fXrKcgf4EPPL3+XNkGWRPsZBvBNAfx5EcEQrRCWk7op8+fbpSwL0A8gjqLl3jtpdolH2um5LTK79sOlsgnAkL7X7JVcdJ3dLe+xQ7Jn3rpbNvrR+YsggHAog9Of/pKp3jtoeiohgxQz4w9WXH4GgA3dADL7T3rpkp0t77FDsmfeuls2+fjkFOPfXUvf1N/sOYcT4R/1THmO0tEomkY7Hac0zn6HQsTADwNEZSS3vvmpnOau99ih2TvvXS2bee8+hPrlgKYADAD8+e/9RkAKxjXNHB5nM+CrEGwA8wmipNxxFCiN1SVVXld93aZ03n6FQW8lhEeTVms90Ru3eclFw2VyPpWy+dfctFzcSuI1QAeHp3rzsvhBCiM5rNeYjyZ6jko01HEUJ8PTmizxEzk+OsGGI6R6exH44D0IBKvN5RQ7huzdCO2rf4KulbL519y0Sfo1gs5gPsO0zn6ER+CMaz7XWlyu1Ryr6no/Ytvkr61ktn3zLR56i+PsJEiJvO0SlE2QfgVAAdeicpIo515P7FF0nfeknfonOr5NGI8rqOOttGCNG+5Ig+R8xsJRKpsaZzdAqMUwDM6+izbVy3dnxH7l98kfStl86+ZaLPUSwWs4lwiekcncSpYDzT0YMohcs7egzxP9K3Xjr7lok+R5FIJEPE95nOYVyUQyDshwK83PGDcZe6Z0HXJ33rJX2LzirKv0VU7i0qRFciR/Q5mj2bbcdJytINMBaAlo9uu25KlhI0kr710tm3TPQ5Ki2NWQB5+83YBbw/gFFogZYjeqVI3hzUSPrWS2ffMtHnKBKJpIlwpekcRtkYA6Aa36XVOoZjpst0jCOypG+9pG/ROUV5DhbydaZjCCF2jhzR5ygaZZ/jpLrUDc3bVfbTsCeC2/9OUjviOKl2vz2h2DHpWy+dfctEn6OiohgB6Gs6hzEK3wEhg/V4S9eQzNRP11hC+tZNZ98y0edo5MiRLQUFTaeZzmGMhRPBWKjz2vOFhY3jdI0lpG/ddPYtE73I1YlgHR+SEkK0N5noc1RVVeVvairo8I/9d0ov8D4ARoLxos5hGxsL5dZ2GknfeunsWyb6HNXXRxjAGtM5jCjA8SCswPFUq3NYItZyGqfIkr71kr5F57KQ70GUZ5qOIYTYNXJEnyNmJtdNlZnOYYSF40F4Rfew8XjdSN1jepn0rZfOvmWiz1EsFvMx40bTObRbxAPAKIaNSt1DE/HNusf0MulbL519y0Sfo5qaiAJ4nukc2ikcByCBo2mD7qEtS66SqZP0rZf0LTqPhfwAovxn0zGEELtOjuhzxMyW6ybPNJ1DO0I5gEUmhk4kUhNMjOtV0rdeOvuWiT5HsVjMZqYLTOfQ6hUuBjAQabxmJgBdZGZcr5K+9dLXt0z0OYpEIhlm3Go6h1aEYwE4OJE2mhjesnCTiXG9SvrWS/oWnUMl348oyy+jEF2cHNHnKBpln+umvHV6JeMYsKllG8BxUreYGtuLpG+9dPYtE32OiopixAzvfGDqFT4QjFL48aqpCMwUMTW2F0nfeunsWyb6HEUikTSQmWw6hzYWvg3C+ybOn98awcpcaGpsL5K+9dLZN+kaqKNMnjg5ZIEeYqAngZax4gmzHp61fle3E60W8i0g7IXRHjvTSIg9UJc/ordA97PCtbMemlUCUg5ZtN17mua63Y5UVVX5Xbf22fZJ3QVYOBKMN0xGcJyUXDZXI+lbL519d+mJftL5k0oB9J/18KzsBJy27gK4Yle3E63mcjcwwgBeNx1FCLH7fKYD7A4fYwADqwAwAKxtWlu3f/c+vadNm5Z/++23b9nZ7aqqVnfLy2s5mIg+Gz584Ip4fPl+Pl/eQNtWHx56aMlHrpu6pLq6NpKX17hs2LBh9a5bdxgR5w8fPjBeWQm7d++6YEsLN0Qigf9WVa3omZfnG9zcnP545MhBdY6z8kDLUgc0NtqpI47ov95xVgyxLF9RUZFaUlJS0uQ4qbBSSpWVlbpvvLGysEcPdUhLi7UxEhnwQVuO5mZePXJkYE0ikQzYttW7oaF5+be+ddDnjlN7qGWhYPjwgYk5c0DDhtWF0mk0lpUVv9eWI53OfFJWVlobj9f18/m4b3Ozv3bkyH6fJBIrDrJt315tORKJZOiaxoayf23p9tncvdasLK6ujWQy6c/D4UHLq6pW75uX11KcTtOasrKBq+PxmmKfz963LUc8XnuIz4fCpUsHOgAwbFhdSCk0hULFS958c/le3brlHZTJqPXhcEmqqirVNy+P+rXliMVWDvb7Vc9Nm6z3vvOdAY3xeE0QwNUAkEwmC+rrrUOVSteHQoOWvfXWqt6FhZlAW46qqhUD8/J8+zU3pz8YOXLQxlgsdbDfT93Wrx/olpdDLV5cF2amLcHgwHeXLl1a1NxcOKQtx5IlyQMyGevAdLq5rqzsoI8XL64bxMx7Nzf7/ztyZL+GeLwmaNuWHQwG4suXL89vbMw7rC3HkiUr98lkVIlS1keh0IAPFy/+YACzv08mo1aEwyWfteVobv64OhKJpBcvritry7FkyboemUzjUGb+NBgMJNtyMLesDAYHr2vLwZy/NBg8YPPixTUjlLJ8wWAgvmTJkrxMpsdw5symYLD0/erq2l4ASpnttcFg/1VtOQDUjBhRvKG6OjkMsLrb9qbFhx56aHN1dW0EsJpHjBiw2HU/6k60ZRiADSNGFI9bvLhmf2a7P1Fm1fDhpWurq2tLAfSy7cL3Dz20z6bq6pXDAZU3YkRxrKqqyp+Xt98IQG0eMaJkaVsOopZ1w4cPXum6q/oTZfa3bSt56KEDPnXdmqFEdo+2HK6bKrMslR4+vLS6LUfba891P+hD5B+wzWuvhIj2acux7WsvFov58vL2G9H22kskknvbtjXoy6+9thxtr73CwuZ3DzrooC2umyrLZFSmrKzUbZsDvvzaa8vR9tprmwO2fe1VVsLq3bsu+OXXXluOttde2xzAnLmkuro2su0cwMwcDpc4b7yxsjAvb2165MiRLe0xV3bpI/r25vNleimlKpTiI7Pf8Q9WSlU0N2MYMxMzfTf7deF+AKCUGqOUqqishL3PPmvzlVIVPh+dCABEdj+lVIVtW61n6qQPU0pVFBSkiwGA2XesUqpi/fr8Hq3D/8Cy7FMAIC8vs7dSqsKy0kcBgG37SpVSFX4/Ds5uao3Kfl3QJ/s1n6SUqnjvvfd8/fuvylNKVRCp7wKA308HZPdlRbK5+NDWfQcAwLLsY5RSFZ9/ntkr+zhOy4AnAHgzkN9UlN2XfXQ2R7pEKVUB8CHZr62R2Z/JOiC7L/U9pVRF//6r8sLhD3yt234PAPLzrf1bc4xszXVItq+m0tZ9H62UqujZM92zNdcpzHwuAGzYkG7LcQwAFBSki1t/hsOy/252RClVkZdn9c1+nf132meftfmxWMzO5lAnAUBLS7c+2X3hiOzXdHD238k3CADS6cxR2VyZXtkc1ngAPwCATz8t7JHdlzUaAJqb0wOz+0oPzz7XLst+bR2Y/ZlwvFKqokePHgVz5sBq/d06Ofvczfu29vMtANiyBUOzj/sGA0Amw9/O/oyb9sn+rlnjlKLTs8/t1i37XDou+2+WGZB9bnpEdlt/qDVn/+y+6DilVEVjY0EhM5NSqoI5PS7782/urZSqyGT42/F4Xb9Mxj6oNceQ1hzfUkpVbNmyqXd23+mx2eczpdM9C7PPpeOzv9N0YPa5/lDrtiOyz20Z0Pr46Oz2PbpnH6fT02lrfOvro1drju9kH/MNzj4XQ1t/x49QSlWk0w2trz0+WSlVMWcOrMLCgfnZf0M6ITuO1fbaC2f7SQ9XSlW0tGSKs/vylSulKj79tHCb1571/ey/d7rttXdk62tvUPa5dHB2X9bh2a+79Wl9PY1RSlXEYjG7qGhNXnZcnJh9HVt9s79LdqR134cppSoKCzNtr71TlFIVGzaki1pznGZZ1ikA0LNnumc63bMQIrskM3XilHVofVN56jlTS1q/3qXtvk5VVZXfcVIL2iN3pxflZ7CQrzIdI5GoNXKNHa+SvvXS2XeXPqK/84E7a0C8atrEyacBANs8iRhPtT0+eeLkk6acO6X3N22Xi/r6CBMh3r4/QSdFOBwW/mM8BnHMdAYvkb710tl3l57oAYAy9k8ZNH3qT6asJCAEy7qm7TELNJN9fPA3bZeL0aMpHQwGrmzv/J3OS9wPCgegAcZf9KFQ4FLTGbxE+tZL+u6EmNlKJFJjTefocAv5FET5v6ZjAIDr1o43ncFLpG+9dPbd5Y/odYnFYjYRLjGdo8MRRoHwtukYAKAULjedwUukb7109i0TfY5qaiKKCE+aztHhCKPAeMd0DAAgUo+bzuAl0rde0rcwg5lQyeuxkI8wHUUI0X7kiD5Hs2ez7TjJPXvpphKDAPQAwTUdBQBcNyVLCRpJ33rp7Fsm+hyVlsYsgPb0N2MjYLyH0dRkOggAKEXy5qBG0rdeOvuWiT5HkUgkTYQ9+/RKRgQwf1plG2a6zHQGL5G+9ZK+hRlRfgmVPM10DCFE+5Ij+hxFo+xznNo7TefoYCEAjukQbRKJ1L2mM3iJ9K2Xzr679NUrdSoqihGw70Gmc+yKceVjbwf4FIAGEKz95lbO/eQrG73E/QDsi4bO8UZsFg0xncBbpG+99PUtR/Q5GjlyZEtBQdNppnPsCmbMbslkjgCw+Yvf5x4NLfydhgxPu34obt3bj083nYjDmXkfQ1G/oLCwcZzpDF4ifeuls285oveAeYvmvQYA48r/d9IQM1vNGZyabMY1HzSi/5sb4RtcAGt5Ex4b2g03MPPdRNRgLLQQot3IEX2Oqqqq/E1NBc+YztGO+jQp/P6uD1E64T3kz1sPX6we1p/qsN/GDK4FYHyZqrGxUG5tp5H0rZfOvmWiz1F9fYQBrDGdox0dvq4F3R/6CLQ5k71OPwOYvQ60bDN8zQo/MpwPRLzadAYvkb710tm3LN3kaPRoSgM4x3SOdrTXhkzrvRW/5NM0AMIYAEZvPhIKBc4yOb7XSN966exbjuhzxMzkuqky0zna0YpAIdKDCsEW/W++H5APDO4GZTH+YTIcAMTjdSNNZ/AS6VsvnX3LRJ+jWCzmY8aNpnPsivHHjr17XPnYVQC6Aap63LEnvwgg1tPG25cPQMNxPaFCPaCO3AuZSwageWA+HJ+Fe0znJuKbTWfwEulbL519y9JNjmpqImrIkNQ80zl2xdxF8y768veICMx85Wn74qoTeuGgZY3oEchH475+rO5m448APjYQ9Qssi+eazuAl0rdeOvsmXQOJzomZfQDKABwI4BMAMTmtUgjhScxsuW7yTNM5vCSRSE0wncFLpG+9dPYta/Q5isViNjNdYDqHt9BXlpxER5K+9dLXt0z0OYpEIhlm3Go6h5dYFm4yncFLpG+9pG8hhBBCt2iUfa6b6pKnV3ZVjpO6xXQGL5G+9dLZtyzd5KioKEbM2JM+MNXpMVPEdAYvkb710tm3TPQ5ikQiaSAz2XQOL7GszIWmM3iJ9K2X9C2EEELoVlVV5Xfd2mdN5/ASx0nJZXM1kr710tm3LN3sBGbuYTqDlzBTkekMXiJ96yV9CyGEELoxM1VVpfqazuEl8XhdP9MZvET61ktn37J0k6NYLObz+fCI6RxeQsRPmM7gJdK3Xjr7lok+R/X1ESZC3HQOLyHimOkMXiJ96yV9CyGEELoxsxWPJ481ncNLEolkuekMXiJ966Wzb1m6yVEsFrMti642ncNbrGtNJ/AW6VsvfX3LRJ+jmpqIIsKTpnN4CZF63HQGL5G+9ZK+hRBCCN1mz2bbcZKXmM7hJa6butx0Bi+RvvXS2bcs3eSotDRmATTWdA4vUYrGm87gJdK3Xjr7lok+R5FIJKMUX286h7eoa0wn8BbpWy/pWwghhNArGmWf49TeaTqHlyQSqXtNZ/AS6VsvnX3L0k2OiopiBPBBpnN4Cw0xncBbpG+9pO9OaenSpXL9aI2kb72kb72kbyGEEEK3qqoqv+OkFpjO4SWJRO0i0xm8RPrWS2ffskafo/r6CANYYzqHlxDxatMZvET61kv6FkIIIXRjZnLdVJnpHF4Sj9eNNJ3BS6RvvXT2LUs3OYrFYj5m3Gg6h5cQ8c2mM3iJ9K2Xzr5los9RTU1EAfSq6RxeQqQqTWfwEulbL+lbCCGE0I2ZLddNnmk6h5ckEqkJpjN4ifStl86+ZekmR7FYzGamC0zn8Ba6yHQCb5G+9dLXt0z0OYpEIhlm3Go6h5dYFm4yncFLpG+9pG8hhBBCt+ytBGt/ZzqHlzhO8lrTGbxE+tZLZ9+ydJOj7K0E+RjTObyE2So3ncFLpG+9dPYtE32OIpFIGshMNp3DSywrc6HpDF4ifeslfQshhBC6VVVV+V239lnTObzEcVLPmc7gJdK3Xjr7lqWbncDMPUxn8BJmkjvwaCR96yV9CyGEELoxM1VVpfqazuEl8XhdP9MZvET61ktn37J0k6NYLObz+fCI6RxeQsRPmM7gJdK3Xjr7lok+R/X1ESZC3HQOLyHimOkMXiJ96yV9CyGEELoxsxWPJ481ncNLEolkuekMXiJ966Wzb1m6yVEsFrMti642ncNbLLn2ilbSt176+paJPkc1NRFFhCdN5/ASIvW46QxeIn3rJX0LIYQQumUvU5y8xHQOL3Hd1OWmM3iJ9K2Xzr5l6SZH2csU01jTObxEKRpvOoOXSN966exbJvocRSKRjFJ8vekc3qKuMZ3AW6RvvaRvIYQQQq9olH2OU3un6Rxekkik7jWdwUukb7109i1LNzkqKooRwAeZzuEtNMR0Am+RvvWSvjulpUuXyvWjNZK+9ZK+9ZK+hRBCCN2qqqr8jpNaYDqHlyQStYtMZ/AS6VsvnX3LGn2O6usjDGCN6RxeQsSrTWfwEulbL519k66BOsrkiZNDFughBnoSaBkrnjDr4Vnrv7zdlIlTlhLQE0Am+x368cyHZlbqTSuEEPp1+SN6C3Q/K1w766FZJSDlkEXX7Whb8lmjZj40q3/2z85N8sxMjrNC3iXXyHVrhprO4CXSt146++7SE/2k8yeVAug/6+FZzwIA0tZdAFd0xFixWMwH2Hd0xL7F9ill32M6g5dI33rp7LtLT/Q+xgAAqwAwAKxtWlsHoPe0adPyt7c9p1Xl1IlTVkydOPnOK88/f6dObaqpiSiAXt3t0CJnRKrSdAYvkb710tm3T9dAplGGxsx8ZGZq2tnTipCn7m1QBTcCmLTtNrFY6mDbxiMAFoVCgcsSidoKIv4lM88Mh+khx0n5HSdVZVl00YgRxTHHST0DYIDfv/k7mUyP7krxy0R4PxgMnO26dUczq1sBfioUKrkhkUhNIcJ5zLgmHA7My37Klkel05kfjBw5qM5xUq8DaAqFAscnEisOIrKfYOZ/h8Mlv3Cc1A8AXEWEO4PBwP2JRHI6EY0D1ORQqPRtx0k9BSDQ0GAfbVlN/oIC/0KAPgiFis9MJFLfIcJtAJ4NhQK/d5zkxQBdQGRdFwwO/JfjpGYC+BYRKoLBQDKRqF1kWczBYKC8urq2VCmeDeA/oVBgWiKRPIWIrgbonlCo+B7HSf0WwKlKYVpZWeA/jpP8O0CDCgqaRm/atFfG52t+lRnJcDhQ4brJI5hpFhHNDQaLr3Xd1AXMuBjAH0KhwDOOk/wrQEcyZ84Khwctd5zUQgB+ANfE4zXFlmU9DdA7oVDxJMepHQfwdCK+PxgsudN1k79mph8SWZcEgwNfSyRSjxNhCDOf0KNHunHzZv/rAFaGQoHT4vG6kZal7gIwPxQK/M5xkucBNIWIbggGi59ynNQtAI6xLP7xiBElSx0n+RIR9QgGA99x3VX9mdPPEiEeDAYurK5OjVEK1wN4KBQKzEwkaq8k4jOY+fJwuKTScVKPADh4yxbf9xoaDtzYq1ftmwCvDoVKvu84qTCAe5nxYjgc+I3r1p7LzD8j4puCwZInE4nUn4kwGqCfhELFSxwn9TzAvUOhksOrqlJ9fT7MBeCGQoGfJhK1JxLxn5j50XC45C+Ok7oMwFlE9MtgsPgVx6l9EODhROmTR4wY9LHr1r5DhLXBYGBsdfXK4UplHgR4QShU8qtEIjWBCJcy49ZwOPCY49T+CeATbVudP3x4abXjpJ4DcEAwWDyqunrFfsy++UT0bjBYPNF1U8cx40YiPBkMBm5KJJK/IKIfA+qqUKj0JcdJ3QcgpJT1/bKygasdJ/k2YH0aChWfFI/XHmJZ/DAzouFw4ArHSf4IoCuI6LZgsPhh1039nhknAfi/UCiQSCSS/ySiAzdsKP5Wnz6r9mppybwE4L+hUOCceDx5rGXRzQDNCYWKZzhOaiqAiZaFq0eMCDzvOMm7AYrYdvq04cMHr3Td1BsANgeDgRNdt2Yos/UYEb0WDBZf8xer2QAABO1JREFU4rq1pzPzrwCeFQqVPOg4yWsBGstsTQqHB77jOKl/ABjYvXvLkRs3WgU+n/0KgOWhUOAsx6k9CuC/EOEfwWDgj66bnMRMPyWia4PB4rnMtJ/jpKqY+fRwuCTlOKnXALSEQoHjYrGVg4maNpaVHfT/7d3fa1tlHMfxz3PStV3H8Fe1NyoOquImzJILUUH0ZlpnsYOmFcU03Yay9GQXirsc/gFjCEk6V4YEf4GJSG9sqxcKeqOMsVGEwZQ5WYeICrKpo22Sx4usgqXNTlj6nGbn/bpMnuZ8+zlfnoSnJd/fmrL/NeNFXPFT6bRkXpekqqxvPXsuVvW+zRXyPZKs/4q/TTH7Xa6Qv6ve6xxMpp+ueuZorpDvC3pta603N3dheOfObQwfceT06Qsv9fXdx3AGR8jbLZd5t9Qn+lxhYkLS/87J/bH0fCaV3pMtTHxqY/aAZ/XJ8nPpVPpZUzUnt7RdXbzqbbkleyI7n8lkOipXqklJZxq59qlTp2Jtbd37JaZMuWNek8TG4wx5u+Uu75Y+o5ckU4ntszJv+aPjF430iDzvv6/+9GRyts0+tFDeepstV2f91Pgle7n6o6RN7bH2Nxq5TjwerxhjTzT9F0Ad9njYFUQLebtF3gAAuFUbJfjz4bDriJLaH77gCnm75TLvlj+6caU2StA+GXYdUWKt91TYNUQJebvlMm82+oDi8XhZqqTDriNKPK/yatg1RAl5u0XeAAAAAAAAAAAAAAAAwEbWUl9q5hKTq9wImnPQdaiPvnbLH/WzMvYFSfcsti/dOTk5+ftq69a7v/k/+jW4mlwVdUFzbuR+YG30tVvGmGJbeelRSX/XW7fe/c1GvwqXk6uiLGjO3I/mIEf3soXsN29/MPlLvTUu7gsb/SpcTq6KsqA5N3o/sDr6emNy0d9s9DfIVEx/rpC/3yx5fcaY269NrgJaGn19c2mpwSPrZeXkqrKx52K1d1kjyfZ09twr2T+y2ezCyp/NvZ/7SZKyH2YvH0ym37GeOeq2+tZVNroYJOeg61BfIznS1+646G8+0as2uSpXyPfmCvneicLE7LF3j52XsfOZVHqPJNmYPWBWTK4aT47fcWjv3q2Z/Zm7JSmTyXRUPJO0DU6uirKgOV9vHYKhrzcWl/3NRr8GV5Oroi5Iztdbh+Doa7f8sfHjfmp8XlJX+2L7nJ/yP19+jv4GAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAANiATdgFAqxscHLy1Y3HTGUlffTxdGpOkgYGBrs2VzpOSzhanS0PhVoio88IuAGh1U1NTfxpVR6z08kj/UFKSuiqdeUmb7aL2hVweoFjYBQA3g+9/OHtpR+/2f2TMkR0PbO+UTNpUvf7iF8XzYdcGcHQDNI9JPJeYNdIuGfNm8bPikbALAiSOboCmSfQnuo3Vw5Iqsnow7HqAZWz0QHMYY/SeNfrVWu2SbGpk99CLYRcFSGz0QFMM7x4+JOmJqlcdKc2UvpTsYWvNZOKZRG/YtQGc0QM3KPF84jFT1ddWGi1Nlz669nDtvN6q+4r+enxmZmYh1CIBAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAQIf8Cdp3e7YQUGScAAAAASUVORK5CYII=",
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ],
+ "text/html": [
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ "Plot(...)"
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "el = Seg3([1, 2, 3])\n",
+ "el[\"geometry\"] = Node[[-0.2, -0.1], [1.0, 0.8], [0.2, 0.7]]\n",
+ "x = Float64[]\n",
+ "y = Float64[]\n",
+ "for xi in linspace(-1, 1)\n",
+ " X = el(\"geometry\", [xi], 0.0)\n",
+ " push!(x, X[1])\n",
+ " push!(y, X[2])\n",
+ "end\n",
+ "l1 = layer(x=x, y=y, Geom.line)\n",
+ "\n",
+ "X1 = el(\"geometry\", [-1.0], 0.0)\n",
+ "X2 = el(\"geometry\", [1.0], 0.0)\n",
+ "X3 = el(\"geometry\", [0.0], 0.0)\n",
+ "l2 = layer(x=[X1[1],X2[1],X3[1]], y=[X1[2],X2[2],X3[2]], label=[\"1\",\"2\",\"3\"],\n",
+ " Geom.point, Geom.label(;hide_overlaps=false))\n",
+ "\n",
+ "xi = [-0.5]\n",
+ "X = el(\"geometry\", xi, 0.0)\n",
+ "J = JuliaFEM.Core.get_jacobian(el, xi, 0.0)\n",
+ "t = vec(J)/norm(J)\n",
+ "n = [cos(pi/2) -sin(pi/2); sin(pi/2) cos(pi/2)]*t\n",
+ "a = 0.3\n",
+ "l3 = layer(x=[X[1], X[1]+a*t[1]], y=[X[2], X[2]+a*t[2]], label=[\"\",\"t\"],\n",
+ " Geom.label(;hide_overlaps=false), Geom.line, Theme(default_color=colorant\"red\"))\n",
+ "l4 = layer(x=[X[1], X[1]+a*n[1]], y=[X[2], X[2]+a*n[2]], label=[\"\",\"n\"],\n",
+ " Geom.label(;hide_overlaps=false), Geom.line, Theme(default_color=colorant\"red\"))\n",
+ "\n",
+ "p = plot(l1, l2, l3, l4, Coord.Cartesian(xmin=-0.5, ymin=-0.5, xmax=1.0, ymax=1.0))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Model defined in cartesian coordinates"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Model defined using cartesian coordinates. In general set individual components of vector values using index number, 1=x, 2=y, 3=z."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "INFO: Starting solver simple 2d block\n",
+ "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: block\n",
+ "INFO: dim = 8\n",
+ "INFO: Assembling boundary problems...\n",
+ "INFO: Assembling boundary 1: dirichlet boundary conditions\n",
+ "INFO: Solving system\n",
+ "INFO: UMFPACK: solved in 0.2082831859588623 seconds. norm = 0.16197088596792505\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "(1,true)"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "INFO: timing info for iteration:\n",
+ "INFO: boundary assembly : 0.1141359806060791\n",
+ "INFO: field assembly : 1.1088120937347412\n",
+ "INFO: dump matrices to disk : 1.1920928955078125e-6\n",
+ "INFO: solve problem : 0.38747596740722656\n",
+ "INFO: update element data : 0.030972003936767578\n",
+ "INFO: non-linear iteration : 1.6414198875427246\n",
+ "INFO: solver finished in 1.7901818752288818 seconds.\n"
+ ]
+ }
+ ],
+ "source": [
+ "geometry = Dict{Int64, Node}(\n",
+ " 1 => [0.0, 0.0],\n",
+ " 2 => [1.0, 0.0],\n",
+ " 3 => [1.0, 1.0],\n",
+ " 4 => [0.0, 1.0])\n",
+ "el1 = Quad4([1, 2, 3, 4])\n",
+ "el2 = Seg2([3, 4]) # load\n",
+ "el3 = Seg2([1, 2]) # symmetry 2\n",
+ "el4 = Seg2([4, 1]) # symmetry 1\n",
+ "update!([el1, el2, el3, el4], \"geometry\", geometry)\n",
+ "el1[\"youngs modulus\"] = 900.0\n",
+ "el1[\"poissons ratio\"] = 0.25\n",
+ "el2[\"displacement traction force 2\"] = -100.0\n",
+ "el3[\"displacement 2\"] = 0.0\n",
+ "el4[\"displacement 1\"] = 0.0\n",
+ "problem = PlaneStressLinearElasticityProblem(\"block\")\n",
+ "boundary = DirichletProblem(\"dirichlet boundary conditions\", \"displacement\", 2)\n",
+ "push!(problem, el1, el2)\n",
+ "push!(boundary, el3, el4)\n",
+ "solver = DirectSolver(\"simple 2d block\")\n",
+ "push!(solver, problem)\n",
+ "push!(solver, boundary)\n",
+ "solver.method = :UMFPACK\n",
+ "solver.nonlinear_problem = false\n",
+ "call(solver, 0.0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Test Passed\n",
+ " Expression: isapprox(el1(\"displacement\",[1.0,1.0],0.0),[1 / 36,-1 / 9])"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using JuliaFEM.Test\n",
+ "@test isapprox(el1(\"displacement\", [1.0, 1.0], 0.0), [1/36, -1/9])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "0.11453071182271282"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "norm(el1(\"displacement\", [1.0, 1.0], 0.0))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAF6CAYAAAAXoJOQAAAABmJLR0QA/wD/AP+gvaeTAAAFlElEQVR4nO3YsU0DQRRF0WdMEe6/NiNRBNYSABIF2Frp6pxkfviiG8xl2zEAst7/3Y/TVgDwCte/49j2deIQAF7jvu14O3sFAK8l9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxB32Xb83h9nDgHg6W7brv9DD0CQrxuAuPff97Ht88whADzdbdt1+/m6uZ+7BYAXuG87fN0AxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxF22Hdse2z5P3gLAc922Xf9CD0DUN3uKFY9lvC/DAAAAAElFTkSuQmCC",
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ],
+ "text/html": [
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ "Compose.Context(Measures.BoundingBox{Tuple{Measures.Length{:w,Float64},Measures.Length{:h,Float64}},Tuple{Measures.Length{:w,Float64},Measures.Length{:h,Float64}}}((0.0w,0.0h),(1.0w,1.0h)),Nullable{Compose.UnitBox{S,T,U,V}}(),Nullable{Compose.Rotation{P<:NTuple{N,Measures.Measure}}}(),Nullable{Compose.Mirror}(),Compose.ListNode{Compose.Container}(Compose.Context(Measures.BoundingBox{Tuple{Measures.Length{:w,Float64},Measures.Length{:h,Float64}},Tuple{Measures.Length{:w,Float64},Measures.Length{:h,Float64}}}((0.0w,0.0h),(1.0w,1.0h)),Nullable{Compose.UnitBox{S,T,U,V}}(),Nullable{Compose.Rotation{P<:NTuple{N,Measures.Measure}}}(),Nullable{Compose.Mirror}(),Compose.ListNull{Compose.Container}(),Compose.ListNode{Compose.Form{P<:Compose.FormPrimitive}}(Compose.Form{Compose.SimplePolygonPrimitive{Tuple{Measures.Length{:cx,Float64},Measures.Length{:cy,Float64}}}}([Compose.SimplePolygonPrimitive{Tuple{Measures.Length{:cx,Float64},Measures.Length{:cy,Float64}}}([(0.0cx,0.0cy),(1.0277777777777777cx,0.0cy),(1.0277777777777777cx,0.8888888888888888cy),(0.0cx,0.8888888888888888cy)])],symbol(\"\")),Compose.ListNull{Compose.Form{P<:Compose.FormPrimitive}}()),Compose.ListNode{Compose.Property{P<:Compose.PropertyPrimitive}}(Compose.Property{Compose.StrokePrimitive}([Compose.StrokePrimitive(RGBA{Float64}(0.0,0.0,0.0,1.0))]),Compose.ListNode{Compose.Property{P<:Compose.PropertyPrimitive}}(Compose.Property{Compose.FillPrimitive}([Compose.FillPrimitive(RGBA{Float64}(0.0,0.0,0.0,0.0))]),Compose.ListNode{Compose.Property{P<:Compose.PropertyPrimitive}}(Compose.Property{Compose.LineWidthPrimitive}([Compose.LineWidthPrimitive(1.0mm)]),Compose.ListNull{Compose.Property{P<:Compose.PropertyPrimitive}}()))),0,false,false,false,false,nothing,nothing,0.0,symbol(\"\")),Compose.ListNull{Compose.Container}()),Compose.ListNode{Compose.Form{P<:Compose.FormPrimitive}}(Compose.Form{Compose.SimplePolygonPrimitive{Tuple{Measures.Length{:cx,Float64},Measures.Length{:cy,Float64}}}}([Compose.SimplePolygonPrimitive{Tuple{Measures.Length{:cx,Float64},Measures.Length{:cy,Float64}}}([(0.0cx,0.0cy),(1.0cx,0.0cy),(1.0cx,1.0cy),(0.0cx,1.0cy)])],symbol(\"\")),Compose.ListNull{Compose.Form{P<:Compose.FormPrimitive}}()),Compose.ListNode{Compose.Property{P<:Compose.PropertyPrimitive}}(Compose.Property{Compose.StrokePrimitive}([Compose.StrokePrimitive(RGBA{Float64}(0.0,0.0,0.0,1.0))]),Compose.ListNode{Compose.Property{P<:Compose.PropertyPrimitive}}(Compose.Property{Compose.FillPrimitive}([Compose.FillPrimitive(RGBA{Float64}(0.0,0.0,0.0,0.0))]),Compose.ListNode{Compose.Property{P<:Compose.PropertyPrimitive}}(Compose.Property{Compose.LineWidthPrimitive}([Compose.LineWidthPrimitive(1.0mm)]),Compose.ListNull{Compose.Property{P<:Compose.PropertyPrimitive}}()))),0,false,false,false,false,nothing,nothing,0.0,symbol(\"\"))"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using Compose\n",
+ "\n",
+ "X = el1(\"geometry\")\n",
+ "u = el1(\"displacement\", 0.0)\n",
+ "x = X+u\n",
+ "\n",
+ "root = context()\n",
+ "p1 = polygon([tuple(X[i][1], X[i][2]) for i=1:4])\n",
+ "p2 = polygon([tuple(x[i][1], x[i][2]) for i=1:4])\n",
+ "undeformed = compose(root, p1, linewidth(1mm), fill(nothing), stroke(\"black\"))\n",
+ "deformed = compose(root, p2, linewidth(1mm), fill(nothing), stroke(\"black\"))\n",
+ "compose(undeformed, deformed)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Using normal tangential coordinates\n",
+ "\n",
+ "Sometimes it's more convenient to set boundary conditions using normal-tangential coordinates."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "rot (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "rot(ϕ) = [cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "ename": "LoadError",
+ "evalue": "LoadError: MethodError: `*` has no method matching *(::Function, ::Array{Float64,1})\nClosest candidates are:\n *(::Any, ::Any, !Matched::Any, !Matched::Any...)\n *{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},S}(!Matched::Union{DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}}, ::Union{DenseArray{S,1},SubArray{S,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}})\n *{TA,TB}(!Matched::Base.LinAlg.AbstractTriangular{TA,S<:AbstractArray{T,2}}, ::Union{DenseArray{TB,1},DenseArray{TB,2},SubArray{TB,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},SubArray{TB,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}})\n ...\nwhile loading In[9], in expression starting on line 10",
+ "output_type": "error",
+ "traceback": [
+ "LoadError: MethodError: `*` has no method matching *(::Function, ::Array{Float64,1})\nClosest candidates are:\n *(::Any, ::Any, !Matched::Any, !Matched::Any...)\n *{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},S}(!Matched::Union{DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}}, ::Union{DenseArray{S,1},SubArray{S,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}})\n *{TA,TB}(!Matched::Base.LinAlg.AbstractTriangular{TA,S<:AbstractArray{T,2}}, ::Union{DenseArray{TB,1},DenseArray{TB,2},SubArray{TB,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},SubArray{TB,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}})\n ...\nwhile loading In[9], in expression starting on line 10",
+ ""
+ ]
+ }
+ ],
+ "source": [
+ "using JuliaFEM.Core: Quad4, Seg2, Seg3, PlaneStressLinearElasticityProblem, DirichletProblem, update!, DirectSolver\n",
+ "typealias Node Vector{Float64}\n",
+ "using JuliaFEM.Core: calculate_normal_tangential_coordinates!\n",
+ "\n",
+ "# 2d rotation matrix\n",
+ "#ϕ = -pi/4\n",
+ "ϕ = 0\n",
+ "rmat(ϕ) = [cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)]\n",
+ "\n",
+ "geometry = Dict{Int64, Node}(\n",
+ " 1 => rmat*[0.0, 0.0],\n",
+ " 2 => rmat*[1.0, 0.0],\n",
+ " 3 => rmat*[1.0, 1.0],\n",
+ " 4 => rmat*[0.0, 1.0])\n",
+ "el1 = Quad4([1, 2, 3, 4])\n",
+ "el2 = Seg2([3, 4]) # load\n",
+ "el3 = Seg2([1, 2]) # symmetry 2\n",
+ "el4 = Seg2([4, 1]) # symmetry 1\n",
+ "update!([el1, el2, el3, el4], \"geometry\", geometry)\n",
+ "calculate_normal_tangential_coordinates!([el2, el3, el4], 0.0)\n",
+ "\n",
+ "el1[\"youngs modulus\"] = 900.0\n",
+ "el1[\"poissons ratio\"] = 0.25\n",
+ "# traction force in normal direction, (i.e. pressure load)\n",
+ "el2[\"displacement traction force N\"] = 100.0\n",
+ "# support sides in normal direction\n",
+ "el3[\"displacement 2\"] = 0.0\n",
+ "el4[\"displacement 1\"] = 0.0\n",
+ "problem = PlaneStressLinearElasticityProblem(\"block\")\n",
+ "boundary = DirichletProblem(\"dirichlet boundary conditions\", \"displacement\", 2)\n",
+ "push!(problem, el1, el2)\n",
+ "push!(boundary, el3, el4)\n",
+ "solver = DirectSolver(\"rotated 2d block\")\n",
+ "push!(solver, problem)\n",
+ "push!(solver, boundary)\n",
+ "solver.method = :UMFPACK\n",
+ "solver.nonlinear_problem = false\n",
+ "call(solver, 0.0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "7x7 Array{Float64,2}:\n",
+ " 0.333333 0.0 0.0 0.0 0.0 0.0 0.166667\n",
+ " 0.0 0.333333 0.0 0.166667 0.0 0.0 0.0 \n",
+ " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n",
+ " 0.0 0.166667 0.0 0.333333 0.0 0.0 0.0 \n",
+ " 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 \n",
+ " 0.166667 0.0 0.0 0.0 0.0 0.0 0.333333"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "full(JuliaFEM.Core.assemble(boundary, 0.0).stiffness_matrix)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "8x1 sparse matrix with 2 Float64 entries:\n",
+ "\t[6, 1] = -50.0\n",
+ "\t[8, 1] = -50.0"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "K = sparse(JuliaFEM.Core.assemble(problem, 0.0).stiffness_matrix)\n",
+ "f = sparse(JuliaFEM.Core.assemble(problem, 0.0).force_vector)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 264,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Array(Float64,(8,8)) 8x8 Array{Float64,2}:\n",
+ " 4.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0\n",
+ " 0.0 4.0 0.0 1.0 0.0 0.0 0.0 0.0\n",
+ " 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 1.0 0.0 2.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\n",
+ " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 1.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0\n",
+ " 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ "Array(Float64,(8,8)) 8x8 Array{Float64,2}:\n",
+ " 4.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0\n",
+ " 0.0 4.0 0.0 1.0 0.0 0.0 0.0 0.0\n",
+ " 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 1.0 0.0 2.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\n",
+ " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 1.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0\n",
+ " 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ "Array(Float64,(8,8)) 8x8 Array{Float64,2}:\n",
+ " 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\n",
+ " 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\n",
+ " 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\n",
+ " 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\n"
+ ]
+ }
+ ],
+ "source": [
+ "# normal version\n",
+ "m = 1/6*[2 1; 1 2]*6\n",
+ "C1 = spzeros(8, 8)\n",
+ "C1[[1,3],[1,3]] += m\n",
+ "C1[[2,4],[2,4]] += m\n",
+ "C1[[1,7],[1,7]] += m\n",
+ "C1[[2,8],[2,8]] += m\n",
+ "\n",
+ "D = spzeros(8, 8)\n",
+ "C2 = spzeros(8, 8)\n",
+ "C2[[1,3],[1,3]] += m\n",
+ "C2[[2,4],[2,4]] += m\n",
+ "C2[[1,7],[1,7]] += m\n",
+ "C2[[2,8],[2,8]] += m\n",
+ "\n",
+ "#C1[:, [3,4]] = (rot(pi/2)'*C1[:, [3,4]]')'\n",
+ "#C2[:, [3,4]] = (rot(pi/2)'*C2[:, [3,4]]')'\n",
+ "#C1[[3,4], :] = rot(pi/2)'*C1[[3,4], :]\n",
+ "#C2[[3,4], :] = rot(pi/2)'*C2[[3,4], :]\n",
+ "\n",
+ "# eliminate tangential constraint\n",
+ "C1[:, 3] = 0\n",
+ "C2[:, 3] = 0\n",
+ "#C1[3, :] = 0\n",
+ "#C2[3, :] = 0\n",
+ "\n",
+ "C1[:, 8] = 0\n",
+ "C2[:, 8] = 0\n",
+ "#C1[8, :] = 0\n",
+ "#C2[8, :] = 0\n",
+ "\n",
+ "dump(round(full(C1), 3))\n",
+ "dump(round(full(C2), 3))\n",
+ "dump(round(full(D), 3))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 371,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Array(Float64,(8,8)) 8x8 Array{Float64,2}:\n",
+ " 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 2.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0\n",
+ " 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0\n",
+ " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0\n",
+ "Array(Float64,(8,8)) 8x8 Array{Float64,2}:\n",
+ " 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 0.0 1.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\n",
+ " 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\n",
+ " 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0\n",
+ " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ "Array(Float64,(8,8)) 8x8 Array{Float64,2}:\n",
+ " 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\n",
+ " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0\n",
+ " 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 1.0\n"
+ ]
+ }
+ ],
+ "source": [
+ "# biorthogonal version\n",
+ "m = [1 0; 0 1]\n",
+ "#m = 1/6*[2 1; 1 2]\n",
+ "C1 = spzeros(8, 8)\n",
+ "C1[[1,3],[1,3]] += m\n",
+ "C1[[2,4],[2,4]] += m\n",
+ "C1[[1,7],[1,7]] += m\n",
+ "C1[[2,8],[2,8]] += m\n",
+ "\n",
+ "C1[[3,5],[3,5]] += m\n",
+ "C1[[4,6],[4,6]] += m\n",
+ "C1[[7,5],[7,5]] += m\n",
+ "C1[[8,6],[8,6]] += m\n",
+ "\n",
+ "D = spzeros(8, 8)\n",
+ "C2 = spzeros(8, 8)\n",
+ "#C2[[1,3],[1,3]] += m\n",
+ "#C2[[2,4],[2,4]] += m\n",
+ "#C2[[1,7],[1,7]] += m\n",
+ "#C2[[2,8],[2,8]] += m\n",
+ "#C2[[3,4],:] = rot(pi/2)*C2[[3,4],:]\n",
+ "\n",
+ "# working\n",
+ "#=\n",
+ "C2[1, 1] = 1.0\n",
+ "C2[2, 2] = 1.0\n",
+ "D[3, 3] = 1.0\n",
+ "C2[4, 4] = 1.0\n",
+ "C2[7, 7] = 1.0\n",
+ "D[8, 8] = 1.0\n",
+ "=#\n",
+ "\n",
+ "tangents = Vector{Float64}[\n",
+ " [1.0, 0.0],\n",
+ " [1.0, 0.0],\n",
+ " [0.0, 1.0],\n",
+ " [0.0, 1.0]]\n",
+ "\n",
+ "normals = Vector{Float64}[rot(pi/2)*t for t in tangents]\n",
+ "\n",
+ "#t1 = [1, 0]\n",
+ "#n1 = rot(pi/2)*t1\n",
+ "# both N and T fixed\n",
+ "C2[1, [1,2]] = normals[1]\n",
+ "C2[2, [1,2]] = tangents[1]\n",
+ "\n",
+ "#t2 = [1, 0]\n",
+ "#n2 = rot(pi/2)*t2\n",
+ "# normal direction fixed, tangential allowed to move\n",
+ "C2[3, [3,4]] = normals[2]\n",
+ "D[4, [3,4]] = tangents[2]\n",
+ "#C1[4, :] = 0\n",
+ "\n",
+ "#t3 = [0, 1]\n",
+ "#n3 = rot(pi/2)*t3\n",
+ "# allowed to move in n and t directions\n",
+ "#D[5, [5,6]] = normals[3]\n",
+ "#D[6, [5,6]] = tangents[3]\n",
+ "D[5,5] = 1.0\n",
+ "D[6,6] = 1.0\n",
+ "\n",
+ "#t4 = [0, -1]\n",
+ "#n4 = rot(pi/2)*t4\n",
+ "# normal fixed, tangential free\n",
+ "C2[7, [7,8]] = normals[4]\n",
+ "D[8, [7,8]] = tangents[4]\n",
+ "#C2[6, :] = 0\n",
+ "\n",
+ "#D[4,:] = 0\n",
+ "#D[8,:] = 0\n",
+ "#C1[4,:] = 0\n",
+ "#C1[:,4] = 0\n",
+ "#C1[8,:] = 0\n",
+ "#C1[:,8] = 0\n",
+ "\n",
+ "\n",
+ "dump(round(full(C1), 3))\n",
+ "dump(round(full(C2), 3))\n",
+ "dump(round(full(D), 3))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 372,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "INFO: size of A = (16,16)\n",
+ "INFO: total dim = 16\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "2x8 Array{Float64,2}:\n",
+ " 0.0 0.0277778 0.0277778 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.0 0.0 -0.111111 -0.111111 -25.0 -25.0 0.0 0.0"
+ ]
+ },
+ "execution_count": 372,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "INFO: solving\n",
+ "INFO: solved\n"
+ ]
+ }
+ ],
+ "source": [
+ "g = spzeros(8, 1)\n",
+ "A = [K C1'; C2 D]\n",
+ "b = [f; g; h]\n",
+ "\n",
+ "info(\"size of A = $(size(A))\")\n",
+ "E = spzeros(2, 16)\n",
+ "E[2,16] = 1.0\n",
+ "h = spzeros(2, 1)\n",
+ "#A = [A E'; E spzeros(2,2)]\n",
+ "#b = [b; h]\n",
+ "#full(A)\n",
+ "\n",
+ "ENV[\"COLUMNS\"] = 300\n",
+ "\n",
+ "dim = size(A, 1)\n",
+ "info(\"total dim = $dim\")\n",
+ "\n",
+ "nz1 = sort(unique(rowvals(A)))\n",
+ "nz2 = sort(unique(rowvals(A')))\n",
+ "A = A[nz1, nz2]\n",
+ "b = b[nz1]\n",
+ "u = zeros(dim)\n",
+ "T = full(A)\n",
+ "T[abs(T) .< 1.0e-9] = 0\n",
+ "info(\"solving\")\n",
+ "sol = T \\ full(b)\n",
+ "info(\"solved\")\n",
+ "\n",
+ "u[nz1] = sol\n",
+ "u[abs(u).<1.0e-9] = 0.0\n",
+ "u = reshape(full(u), 2, round(Int, length(u)/2))\n",
+ "full(u)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 373,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2-element Array{Float64,1}:\n",
+ " 0.0277778\n",
+ " -0.111111 "
+ ]
+ },
+ "execution_count": 373,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "u[:,3]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 374,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "INFO: displacement on tip: [0.027777777777777783,-0.11111111111111115], magnitude = 0.11453071182271282\n"
+ ]
+ }
+ ],
+ "source": [
+ "u_tip = el1(\"displacement\", [1.0, 1.0], 0.0)\n",
+ "info(\"displacement on tip: $u_tip, magnitude = $(norm(u_tip))\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 375,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Test Passed\n",
+ " Expression: isapprox(norm(u_tip),norm([1 / 36,-1 / 9]))"
+ ]
+ },
+ "execution_count": 375,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using JuliaFEM.Test\n",
+ "@test isapprox(norm(u_tip), norm([1/36, -1/9]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 376,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "16x16 Array{Int64,2}:\n",
+ " 440 150 -260 -30 -220 -150 40 30 2 0 0 0 0 0 0 0\n",
+ " 150 440 30 40 -150 -220 -30 -260 0 2 0 0 0 0 0 0\n",
+ " -260 30 440 -150 40 -30 -220 150 0 0 2 0 0 0 0 0\n",
+ " -30 40 -150 440 30 -260 150 -220 0 0 0 2 0 0 0 0\n",
+ " -220 -150 40 30 440 150 -260 -30 0 0 0 0 2 0 0 0\n",
+ " -150 -220 -30 -260 150 440 30 40 0 0 0 0 0 2 0 0\n",
+ " 40 -30 -220 150 -260 30 440 -150 0 0 0 0 0 0 2 0\n",
+ " 30 -260 150 -220 -30 40 -150 440 0 0 0 0 0 0 0 2\n",
+ " 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
+ " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
+ " 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0\n",
+ " 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0\n",
+ " 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0\n",
+ " 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0\n",
+ " 0 0 0 0 0 0 -1 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 1"
+ ]
+ },
+ "execution_count": 376,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "round(Int, T)"
+ ]
+ },
+ {
+ "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
+}
diff --git a/src/directsolver.jl b/src/directsolver.jl
index 040f7c8..45fc641 100644
--- a/src/directsolver.jl
+++ b/src/directsolver.jl
@@ -141,9 +141,10 @@ end
""" Call solver to solve a set of problems. """
function call(solver::DirectSolver, time::Number=0.0)
+ info("Starting solver $(solver.name)")
info("# of field problems: $(length(solver.field_problems))")
info("# of boundary problems: $(length(solver.boundary_problems))")
- (length(solver.field_problems) != 0) || error("no field problems defined.")
+ (length(solver.field_problems) != 0) || error("no field problems defined for solver, use push!(solver, problem, ...) to define field problems.")
timing = Dict{ASCIIString, Float64}()
tic(timing, "solver")
@@ -202,7 +203,7 @@ function call(solver::DirectSolver, time::Number=0.0)
info("Assembling field problems...")
field_assembly = Assembly()
for (i, problem) in enumerate(solver.field_problems)
- info("Assembling body $i...")
+ info("Assembling body $i: $(problem.name)")
append!(field_assembly, assemble(problem, time))
end
K = sparse(field_assembly.stiffness_matrix)
@@ -217,7 +218,7 @@ function call(solver::DirectSolver, time::Number=0.0)
info("Assembling boundary problems...")
boundary_assembly = Assembly()
for (i, problem) in enumerate(solver.boundary_problems)
- info("Assembling boundary $i...")
+ info("Assembling boundary $i: $(problem.name)")
append!(boundary_assembly, assemble(problem, time))
end
C = sparse(boundary_assembly.stiffness_matrix, dim, dim)
diff --git a/src/elements.jl b/src/elements.jl
index 13bfe20..3614e1a 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -243,34 +243,48 @@ end
""" Calculate local normal-tangential coordinates for element. """
function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::Real)
- proj(u, v) = dot(v, u) / dot(u, u) * u
ntcoords = Matrix[]
refcoords = get_reference_element_coordinates(E)
x = element("geometry", time)
for xi in refcoords
dN = get_dbasis(E, xi)*x
- normal = cross(dN[:,1], dN[:,2])
- normal /= norm(normal)
- u1 = normal
- j = indmax(abs(u1))
- v2 = zeros(3)
- v2[mod(j,3)+1] = 1.0
- u2 = v2 - proj(u1, v2)
- u3 = cross(u1, u2)
- tangent1 = u2/norm(u2)
- tangent2 = u3/norm(u3)
- push!(ntcoords, [normal tangent1 tangent2])
+ n, m = size(dN)
+ @assert n != m # if n == m -> this is not manifold
+ if m == 1 # plane case
+ tangent = dN / norm(dN)
+ normal = [-tangent[2] tangent[1]]'
+ push!(ntcoords, [normal tangent])
+ elseif m == 2
+ normal = cross(dN[:,1], dN[:,2])
+ normal /= norm(normal)
+ u1 = normal
+ j = indmax(abs(u1))
+ v2 = zeros(3)
+ v2[mod(j,3)+1] = 1.0
+ u2 = v2 - dot(u1, v2) / dot(v2, v2) * v2
+ u3 = cross(u1, u2)
+ tangent1 = u2/norm(u2)
+ tangent2 = u3/norm(u3)
+ push!(ntcoords, [normal tangent1 tangent2])
+ else
+ error("calculate_normal_tangential_coordinates!(): n=$n, m=$m")
+ end
end
element["normal-tangential coordinates"] = ntcoords
end
-
-""" Pick values from nodes and set to element according to connectivity. """
-function update(element::Element, field_name::ASCIIString, data::Union{Vector, Dict})
- element[field_name] = [data[i] for i in get_connectivity(element)]
-end
-function update(elements::Vector{Element}, field_name::ASCIIString, data::Union{Vector, Dict})
-# info("update $field_name for $(length(elements)) elements.")
+function calculate_normal_tangential_coordinates!{E}(elements::Vector{Element{E}}, time::Real)
for element in elements
- update(element, field_name, data)
+ calculate_normal_tangential_coordinates!(element, time)
+ end
+end
+
+""" Pick values from nodes and set to element according to connectivity. """
+function update!(element::Element, field_name::ASCIIString, data::Union{Vector, Dict})
+ element[field_name] = [data[i] for i in get_connectivity(element)]
+end
+function update!(elements::Vector{Element}, field_name::ASCIIString, data::Union{Vector, Dict})
+# info("update $field_name for $(length(elements)) elements.")
+ for element in elements
+ update!(element, field_name, data)
end
end
diff --git a/src/fields.jl b/src/fields.jl
index c0449dd..9814849 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -180,8 +180,12 @@ for op = (:+, :*, :/, :-)
@eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data)
@eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data)
@eval ($op)(field1::DCTI, field2::DCTI) = ($op)(field1.data, field2.data)
- @eval ($op)(field::DCTI, k) = ($op)(field.data, k)
- @eval ($op)(k, field::DCTI) = ($op)(field.data, k)
+ @eval ($op)(field::DCTI, k::Number) = ($op)(field.data, k)
+ @eval ($op)(k::Number, field::DCTI) = ($op)(field.data, k)
+end
+
+function Base.(:+)(f1::DVTI, f2::DVTI)
+ return DVTI(f1.data + f2.data)
end
function Base.vec(field::DVTI)
diff --git a/src/linear_elasticity.jl b/src/linear_elasticity.jl
index a57245b..e32c71d 100644
--- a/src/linear_elasticity.jl
+++ b/src/linear_elasticity.jl
@@ -118,6 +118,23 @@ function assemble!{E<:CG, P<:PlaneStressLinearElasticityProblem}(assembly::Assem
L = w*T*N*norm(J)
add!(assembly.force_vector, gdofs, vec(L))
end
+ for dim in 1:problem.dim
+ if haskey(element, "displacement traction force $dim")
+ T = element("displacement traction force $dim", ip, time)
+ ldofs = gdofs[dim:problem.dim:end]
+ L = w*T*N*norm(J)
+ add!(assembly.force_vector, ldofs, vec(L))
+ end
+ end
+ if haskey(element, "displacement traction force N")
+ # surface pressure
+ p = zeros(2)
+ p[1] = element("displacement traction force N", ip, time)
+ R = element("normal-tangential coordinates", ip, time)
+ T = R'*p
+ L = w*T*N*norm(J)
+ add!(assembly.force_vector, gdofs, vec(L))
+ end
end
end
diff --git a/test/test_mortar.jl b/test/test_mortar.jl
index 980f0fe..3add883 100644
--- a/test/test_mortar.jl
+++ b/test/test_mortar.jl
@@ -6,7 +6,7 @@ module MortarTests
using JuliaFEM.Test
using JuliaFEM.Core: Element, Seg2, Quad4, Tri3, Hex8, MortarProblem, Assembly, assemble!,
- get_connectivity, update
+ get_connectivity, update!
using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem, DirectSolver
# 2d stuff
@@ -727,7 +727,7 @@ function test_3d_problem()
l2u = Quad4([5, 6, 7, 8])
u2l = Quad4([9, 10, 11, 12])
elements = Element[el1, el2, sym121, sym131, sym132, sym231, sym232, force, l2u, u2l]
- update(elements, "geometry", nodes)
+ update!(elements, "geometry", nodes)
el1["youngs modulus"] = el2["youngs modulus"] = 900.0
el1["poissons ratio"] = el2["poissons ratio"] = 0.25
sym121["displacement 3"] = 0.0
@@ -859,8 +859,8 @@ end
mel7 = Quad4([18, 19, 23, 22])
mel8 = Quad4([19, 20, 24, 23])
mel9 = Quad4([20, 21, 25, 24])
- update(Element[sel1, sel2, sel3, sel4, mel1, mel2, mel3,
- mel4, mel5, mel6, mel7, mel8, mel9], "geometry", nodes)
+ update!(Element[sel1, sel2, sel3, sel4, mel1, mel2, mel3,
+ mel4, mel5, mel6, mel7, mel8, mel9], "geometry", nodes)
calculate_normal_tangential_coordinates!(sel1, 0.0)
calculate_normal_tangential_coordinates!(sel2, 0.0)
calculate_normal_tangential_coordinates!(sel3, 0.0)