FEniCS-FE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What is Finite Element Theory?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Galerkin Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a bit of maths that shows you can legitimately model real world continuous problems as discrete, finite ones - and, even more, that, in certain ways, you can make it more accurate by making it bigger."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This turns it into a problem of lots (but finite quantities) of real numbers. Computers are perfect for these kind of problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image1](https://i2.wp.com/fenicsproject.org/wp-content/uploads/2016/10/stokesexample.png?resize=768%2C359&ssl=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://fenicsproject.org/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image2](https://i0.wp.com/fenicsproject.org/wp-content/uploads/2017/10/tc_vm.png?resize=372%2C280)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://fenicsproject.org/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image3](http://cbc.simula.no/pub/cm/figs/twist.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://cbc.simula.no/pub/cm/fenics.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image4](http://www.fenics-hpc.org/wp-content/uploads/2017/06/cropped-cropped-cyl_pylon_aoa22.56_trip_vortq01b.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://www.fenics-hpc.org/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image5](https://ecmiindmath.files.wordpress.com/2015/11/fenics-3.png?w=610)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://fenicsproject.org/pub/graphics/fenics-poster-2015.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![klein-sensitivity](data:image/jpeg;base64,/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEAAkGBxATEhAQDxIQEBUVFRUPDxAQFRUPDxAQFRUWFhUVFRUYHSggGBolHRUVITMhJSkrLi4uFx8zODMuNygtLisBCgoKDg0OGhAQGi0fHx8tLS0tKy0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLf/AABEIAOEA4QMBEQACEQEDEQH/xAAcAAABBQEBAQAAAAAAAAAAAAAAAgMEBQYHAQj/xABHEAABAwIDBAYGBwQHCQAAAAABAAIDBBEFITEGEkFRBxMiYXGBFEKRobHBIzJSYnKCkjNTotEWJENEY5PCFTVUc3SjsrPw/8QAGgEAAgMBAQAAAAAAAAAAAAAAAAQBAgMFBv/EADURAAICAQMBBQUHBQEBAQAAAAABAgMRBBIhMQUTIkFRMmFxkdEUI0KBobHBM0NS4fDxcjT/2gAMAwEAAhEDEQA/AO4oAEACABAAgAQAIAEAMVVZFGN6WSOMc5HBg9pKhtLqWjCU3iKyVUm2GGjI1lN+WRrvgqd7D1GVoNS/7cvkxpm3GGH+9wj8RLfiFCurfmS+ztUv7bHm7X4ac/TaXzlYD7CVPew9Sr0OpX9uXyY3PtphjBd1ZTn8DxIfY25R3sPUlaHUt42P5GexLpZoWXFMyapPA26mM+Ln9ofpWM9XCPvH6OwtRZ7eIr5/t9TI4r0pV779X1NMOG43rZLci59wfJoS0tZZL2Vg6tfYmlr/AKjcn8l+n1MtW7XVkhvJVVDu7rXtb+lpAHsWblbLrIajVpa+IVx+Wf3K1+KPJBLnEjQkkkeBOio62+Wy3exXCS+RNpdoqlhuyonZ+CV7fgUZmujZlJVS9qEX+S+hqcI6TK+MgSPZUN4tlaN63IPbY37zdaR1VsevInb2bpbOi2v3fR5/g6Js50g0dSWsefRpDkGSEdW48myaHlY2J4BO1aqE+OjOPqezLaluj4l7uv5r/wBNemTmggAQAIAEACABAAgAQAIAEACABAAgAQBlNqNvKWkJjF6iYZGKMgBh/wAR+jfAXPcl7dTCvjqzpaPsu7U+L2Y+r/hef7e85jje3+IT3HW+jsPqU94zbvk+sT4EDuSE9XZLpwel0/Y2nq5a3P3/AE6fuZSeUuJc4lzjq5xLnHxJzKw5fU6kalFYSwhl0ilRJ2JDLnlXSKvA2SrGfLHYob5lUlIuopcs8qKsNyaphW5csWu1SjwitlqyU1GtI5k9TkZ60q20y70W2RQ4l1PJIjes2i6kS4nrGSLKRMikWMkXUjoGxO38lOWw1RdLBoHZukgHdxc0fZ1HDSxZ0+rcHtn0/Y5+s0EbvFDiX6P6P3/P1Ow01QyRjZI3Ne1wDmOabtc06EFdVNNZR56UXFtNYaHVJUEACABAAgAQAIAEACABAAgAQByzb3b5xLqWgfZou2apYe048WxOGg+8NeHM87U6vHhh8z0nZnZKaVt6+C+v0+ZzMrnHp0hp6sjRDDyrokacrlJCSpM2ORR3VZSDoM19WGjdCvVXnlnP1WqxwinkkJTqjg487HJiLqSmT0FBKYthVWawZJjKykbolxLGRcmRLKRJKjKxYZN70bbWGnkbSzu+hkdZhOkErjr3McdeRN+ac0eo2vZLoc/X6XvY95H2l+q+q/7yOyLrHABAAgAQAIAEACABAAgAQAIAwPSjtOYWCjgdaSVt5nDWOE3Fhyc7MdwB5gpLWX93HaurO12PolbPvZrwx6e9/wCjkZC4+T1yG3BWRdMZeFdE5GXhXQbhpwV0Q5HjWobM3IVVS7je9RCO6QrqLtkSglkJNyuhGOEcGc3J5Y3dWM8nqCT1QWQ4wKrNYIlRNWUmMIsKOme82Y1z/wAIJt48kvZOMeZPBbKXUu6bAJz9YNZ+J2fuukZ6ypdOSrmkT49njxkb5NJ+aXetXlEo7oocOz54SNPi0j5qv21ehC1ETqWym0zWwRxVj/pGDq+sAc5r2jJribXBtYG/EX4rtabtaiUUpvD95x9RpW5uVfRmspqmOQb0b2vHNpDh7l1IWRmsxeV7hKUZReGsDquVBAAgAQAIAEACABADVXUNjY+WQ7rWNdI93JrQST7Aobwsloxcmorqz55xSvfUTS1En1pHF5H2Ro1vk0AeS87dY5zcj3empVNcYLyIbgqIYyNOCsidwy8K6DcMPCuiHIaKuUcxcTVDZXcVGKVF3W5JuiGEcXVXbplfdMCWQQGT0ILIWAoNEiywvC5Zj9G3Li85MHn/ACSt+ohV7T/I1jJLg1VDgMEdjKetdyOTB5cfNcq3WWT4jwv1LSsx14LdlQAN1oDRwAyCTcG+WKz1OOg/G8lUcfQVle2S4mo7mTMnaySyNT9mZXvGPCJUemkiysY5CXscHxucxw4tNj/93KsJW0y3RbXwNFblYlyjT4NtdmI6vLgJgLD84GniPZxXe0XbCliF3D9fqZT06lzX8voa5rgQCCCDmCMwRzC7yeRM9QAIAEACABAAgDE9KGOxR0slK2VnXS7rerBvIIt4F5cOALQRnzSursUa2s8nT7L08p3Rm14V5+/yOR2XAPXZEuClEZGnBWQbhl4WiDcMPCuirkMkK5RyPJX7rHHuRFZkkZW2bYNkWg2flmAkLmxtdmCQXOI57vLzW1uthU9qWWjlR087Fu6ZLul2PpxbffK891mNPlYn3pGfadv4Ul+pf7JFdW2WDNlqL90T4vf8il32hqP8v0RV0wQO2Yo/3ZH55Pm5C1+o/wAv0X0MZYj0I7dnqRh3t1zreq512+zitHrb5LGcBBTn7kPzV7WjdZYAZADIAeCzjS3yy07Y1rESOyclaOCQjO1sn06yccmDkWlOxaRrRXJZQRLeNZBNjhWqqIyPthVu5DIowrCenTJTI80C5t+lxyjSMsFjs3j5p3CGY3hJyJ/sSeI+7zHDXnd3sztB1PurOn7f6L2QVqyuv7nQQV6cSBAAgAQB45wAJJAAzJOQAQBybbzpPN3U+GutbsyVQsbniIuFvv8As4FJ26nyj8z0mg7F4Vl6+Efr9Pn6HMcPe50jnPJc43c5ziXOcScySdSufe/CdicVFxS4RcWSJLYkhSVyIcFZBkYeFdEZI0i0RRyGXLRFXIaqWb24zg5wB/CMz7gVaD25l6GF/jxD1Zfwy6DTkufKJs2TI5li4mE2PCdU2ClkhmaqV41mEYbmUtbiBJ3W+ZTtdKXLC25RW2JFatmISlknUgS9jM2XdJGskyC5pYlrFkFrTxJiBBPiiTMSCQ2JbKJAdUh1gMyxJW2ksmVdbBkVxNTTh5RrCWDUbBYuXsdSyG7oheO+rodLflNh4Fq73ZWq72vZLrH9iNRD8a8/3NauqLAgAQByPpT2xc9z6CldZjTu1UjdZHDWIH7I48zloDdDVajnZH8z1fYnZiwtRav/AJX8/T5nLnMSSZ6VxyP0As/yVLeYiGqjtlF/H+C3Shg2eEIKZEOUoMjEi0RVsiyrWJRyIzitUZuQne7bTyB9uQ/mpx4WC5mvcibHMsHEtKRJZULJwFZyFuqVHdisnllfW1Z+qPPuTFdXmwnPYsIisWrEpMfjWbM2i0oWJW1kYNBRxpTfhkYLuljWsbCMFpTsTUJkYJ8TE5CRUkNYm4MqBYtkiBqSNUnEkrayNcvU1cF0ymp630aohn0DXWk743ZP9xJ8QFz9FZ3N69BiK3RcTr4K9cIggDNbf48aSlcYzaWQ9TBza4jtP/KLnx3eaw1FvdQbH+zdL9ovUX0XL+n5nCjAuDvPdKaSwJNOp3h3o26PdLXcjY+By+NlZPcmhXVTzFP0f+ixCXMGBUGbEOVkRkYkWiKtkOZbRMpMhvctkjJyGJ5LFvgrxjlMtGXI5HOquBEpkhs6zcBacjySosPgpUMmaeOSIH8StcC8nkcY5VaM8EuAXKylwiNpfUEWi59sidpfUjUhNlXEuKZEbCuCzp01XaVwWEIT1dpVolMan67CrQvcTkJZKiJI1d8kFbWRJO+PBZGUxqPIrgXR22JjFb5OnbI1fW0dLITc9W1rjzczsO97SvWUS3VxfuF7VibRbrUzOP8ASXXGatMY+rAwRgcOseA95HkWD8i4naVvjUfQ9T2NX3dDn5yf6Lj6mZEC5e86UrANOjeU7wj1NLdrhzGXceHvWkLMNMpOW6LRHpJd5oPtHEFaWRwyIS3RTJCzBiHqyM2RpVpEo2QZyt4mUmQZXLeKF5SIlW7Jp77LatcslSGo5CrOISbJDJVm4mEmXWzmzVZXvcykYHbgvJI87kUd9A53M8gCVpXU5GNtiguSDjWEVNJMaerjMTwN4A2c17Do5jhk5uRzHIjUWROG0rDE1lDUZWDLNFnQhK2kqJfUhSFiLbS3p5gk5xI2FhBUhLygyrqLKnrG81TdOJR1MsYKsc1pDVOPUo4MmxVAXQp1qfmZuJJZIF1adUmUcRZK6ELUymCHVsyUWcoEZPG2ZFcLVxw8msHya/oxfegYPsyTAecjnf6ivQ6J/cRK6j2/kaxNGBwqrd1k9RIc9+aV4PcXuIHkLDyXlNbPN0viev0/gogvchbIEi5lZWHroFCmZ94R5YVpGQd4UE7eqmI0bJ228t71h8/NdCD7yvPmv+Rame2e31JIKyGmhLlKMmR5gtYmciuqExAXmV0pTMReTI9R9U+1aQ6kJkZjlo0ap5RYYbSSTSRwwtMkkjgyNg1c4/AcSeABKootvCMZtRWWfUOx2zsdBSx0zLFw7c0gFutmdbfd4cByAA4J2MVFYRyLJucssibf7JR4jTGPJszLvpZT6klvqn7jrAHyOoCJxUlgmqxwlk+baunkhkfBOwxSRndkjdk5p+Y4gjIggjJITg0dFNS5RZ4UJJDuwRyzkathY6Zw8QwFLOicuiL5jHq8DwxFwJaGm4JBvlYjIhYS06XtMYhXKXKQ82tnOhYPaVm6q16jEdLY/RChU1H7xv6T/NV2Vehf7DY/xL5CxX1TdHMPtCjuaX5MrLQahdHF/ND8O0s7Prtv3gqktBXP2WJ2wtq5shx69V80XNBtk3IPu3x09qSs7MkuYi+6qZpqDHmPtZwSv31LKyp9C7p60Hiuhpu0OcMXlDA7I8ELtV3qaMmjOY5HkUlrFlZJiaHosjtROPOeU+FrN+S7Og/oRIvfiNgnDE4ThrDbtfWv2vxcfevF6ttWSz6nq4y+7WPQt44Ug5Ck5nr4UKZl3hDniW0ZFlMo8aousYQPrDtMP3hw89E9prdksvo+pdvKKjDqveFjkRkQm7q9r46D9Fqsj7ya5Yl5IjyrSJkyvqQmIC80VUmqbQo+o7htWIpoJnDeEckcrmnMOaxwc5pB1BAIt3rSDxJFJLKaPpifYzCpRvOoqQ73a3mRMYXXzvdoBT2EctWTXmyTg+y9DSuL6WmhhcRYva3t7ptcbxzAyGXchJIJWSl1eS3UlAQBErcLp5i108EMxabsMsbJCw82lwNlGEyVJroLmLIonua1rWsa59mjdADQToPBSHVny9RtNgTck5knMknUlcK2WWz2NUMInsCXY1FDzQqM1ieOQjVESYraJbBFa8XsU3DEuJHmO2tHGpxtrWM8P4/9klQuew3Y4t8NPYqW6bK9UcaF0omgwnaVzSGyZd/Arj36BPmI3G2M+GbXD8Wa8DNKV3WUPEuhSyrHQ9xLtNK6veq2HAvjBo+jiO1E3vkmJ/zHD4ALuaD/APPEzteZGnThmcUZHuT1MRFtyeVoH3RI7dPmLHzXj+04bbpHoaJ7qV8C5gYuLJi1ssDj4kLIvuIFTEtoMtGZTVbU7WxiMjE4q7q6hxZlcB5HC5vf4X8126FvqW74FlNwllFjRVgeOR4hL2VODOrVarF7x6QKiCUSvqWreDFpoqJtU5HoJT6jIY5xDWNL3OIa1jQXOc4mwaAMySTay1isso3jk+tcEp5I6emjlO9IyGOORw0MjWAOOXeCnl0ORJpt4JqkgEACABADNZAJI5IybB7XMJGoDgRf3oJTw8nzdW4XNTSGnqGFkjeGe69oyD2H1mnn7bG4XCvrlCXJ6/TXwtjuiwalmNpjgVDRSPHBSjWMiJOFtE1TKx7s0ylwcvtiO7TP3YZZU7rhb1zysM8a0Oujuq21J8oE8D+H4g+IgXJb8FzL9Op/Ecqu8mbWhxISN1XLjuol7i9tfmjfbCEei7o9WSQHxLt74OC9Z2bNS06x7/3ELOpoU+UORbaQ9Rich9WdjJ28r26tw9rAfzLz/bFOXuOroLeNrJtE64XlbFhl71gllqdqhlCTZCqmKZ1YBMz+Itsr1M3hI5zj0n9Zf3Bo/hB+a9HpY/cr8xhLKG4pCCCDYq8op8MmMnF5Rb0eJB3ZfkefApOyhrmJ0KtTGfEuGSZ2XCyi8Gs4FSMOmmmjggYZJJHbkbG6uOupyAABJJyABK6NHi4RzNR4OWd42A6OIKDdnm3Z6q37S30cNxYiIHuNt45nPQGy6EIKJyLbnPjyN0rmIIAEACABAAgCuxzA6arj6qqjbI3Vp0fG77THjNp7wVEoqSwy8LJQeYvDOD7cYNJhlQ2JzuvhkaZIJHWEu6DZzX2yLhcZ2sQR3gczUaOK5jwd3S9pSkvGskClqmPHZOfI5Fc2dcodTr12wsXhZILVnk2UsESpYtYM1jMpZ/rAd9vano+yKdoPNMvgT6crJSwzx0kTWFNRllFDyRqXtj5kpj2F1pjduk5cEjfUprI9TZ5M630bYiHddDfW0zR/C4/+Cc7GnhSrfx/79BfVV7Xk3K7gqc96YqEmCnqm6wybjv8AlyjP+JrPaktdXvr+BrTPZLJm8CxAOaM143VUuLOu8WRyjRxuuFOks/CxCccMZqGLqOvKMsmexVmRSezbI0izl+Pt+neee6f4QPkvQ6R/dI6VSzDJGjctWiGhZKqZtEinrnsy1HIrOdMZG1epnDjqjUbE7TxUtXHUvbcAOjkHrBj7Xcw6XFh4i44qaN1MueUW1GzU17U8P3n0JhWKQVMYmp5GysPrNOh5OGrT3HNdOMlJZRwbK5VvbJYZMVigIAEACABAAgCo2k2jpqKPrKh2Zv1cTc5ZXDg1vzOQvmVSdkYLLNqaJ3SxFHzptvjs9dUmpl7IA3IYgbtijvcC/Fx1J4nkAAEnfvOm9L3S4KONxBuCQVVpPqEZOLyi7w/Fb9mTyckbdPjmJ1NPrFLwz+ZZSsBGSWi8MdfHKKKqhPWxjm4fFPQl4H8BbVzzUywdCl1I8zOJ6wraEzFocuryllEEebmsEaQeDYdHuLdXV05Jyc7qXd4k7I/iLT5KdL93qE/XgbtW+rPodxXoDmGe6QIA/Dq0HhEZfOMiQe9qzuWa2iV1OG4dVGMg8F57U0bkN0XODNvheJBwGa4NlUq5ZQ7OCsWUWpeCF09JqlLh9RCcGiqxOG4Kauqysoqmct2mhtNfmPgSm9HL7vB2NF4q37ipAsnC044PQVBg0eqCuDxSRgs8Bx+qo5BNSymN2W83WORo9V7dHD3jhYqYycXlBOKnHbI+ntm8WFVS09UG7nWxtkLNd1xHabfiAbi/FPJ5RyZx2ya9CyUlQQAIAEACAPnHG6+WqnlqJi4uc47rXf2UYJ3YwOAaDa3O51JK4N9znNtnsNPRGutRj/77yumprhZRswWsryitnpE1Gw51lOCMYyFplMXw0WmG1xHZdmErdVnlHR02px4ZE+eAFzH8jf3EfNYRm0nE21MMweB5zFRM4U4jL2K6kLSiNFabjPAzIVKLRQ5h9SWODmmxaQ5vc4Zj3hTLhpryHqeVhn0N/SCPl8f5Lu70czYxrbs2w6v/AOnlHmWEBTP2WRHqjgMGYXNjFTiT0JtJUuYciufqNKb1XOBpcOxkGwJXGt0zi8oeUoWrnqWhma4LejWyh4ZrKF7NNJcoxe2mDuIErAXbv1gMyGnjZdXTWwbex8MZ7Os2TcJcZ/cxhjTu46VkBBiVtwrKAkhSZNHlkEYPHKUQfUPR5BuYZh4508T/ANbQ/wD1J+PQ5NrzN/E0SkzBAAgBuedjBvPc1g03nENF/EozglJt4QmCrif+zkY/8Dg74FRlMlxlHqsHB9paHqqyri+zK5wHJsn0jB+l7V5/WR22s9hoJ95p4S92PlwVxYlcjMokaeFaxkK2QIE0CYjMSnAjmOy03ZMsYLGhn9V2iWth5o6FFuVtZOaeCwwI6mnZLAh6lCEokeRaIycSNItUCiNsdmrNDNXU0H+2p/3P8L013k/8Se7j6nZekV9sNrTzi3f1ODfmujb7DOXX7SOF0zMgubppeRMlySDGnHBNFDxoI0SF2lTLxm0Taavc3iuVbpceQ3XqpLqWEeKc0o6GnlDKtqn7SINZh9LKS4t3HHVzOz5kaH2LeGovh55+IzXPCxGXHo+f9lZNs6z1JR4Ob8wmY66XnE13Z9CI/Zx37yP3rVa5f4sq4Z9Bs7Pc5G+QJ/kr/bfSJHdx82j0YFF6znv8LMHzUfbJ+SSI21rzbNHU7SV7mNi9IlYxrQxjYz1QDGiwB3LXyA1UvV2y6y/gpGumPswX58/uVUNZUxP66KaVklwd9rjvG2m99odxuEQvcXlMY2KaxJZR3DYXaltdBvO3WzR2bURjIB3B7Qc911suRBHBdum1WRycLV6Z0Tx5PoX1ZVMiY+WVwaxjS97joGgXK0bwssXjFyaiurPnva/HJK2d0slw0XEMRzEUfIDTeNgSeJ7gLcW292Sz5eR6/R0x08VBdfN+r/7oZmVimLOqluiWOETucXh7nOORu4lxta2p8FhqucMWjBQk4pY8y1CTLNCJGqUzCaIkjFsmKTiRJI1qpC8ojYFld8hHgnRyXF+Wvgl3HkZsXeV+9A56EjmSgMvetEjJwI0j1okRtI735H4LRRyzWtcndf6Eu7/1BdXuzm98Xe37L4dXCwNoXuz+72r+5Ws9hmVftI4jQtuAuLRLEi81yT2wrrweUZMHUys45IG/RktZUmWTPRTlIWadF0z3qCk50tF1Nh1BWDTRorWApiqbzRWsUKRR3hopsV6Go701iz30NHejMENTUqtGwcgiNQV89LM2emduPbkb5se06se31mn5AixAK6GnvlB5RW+ELIbZot9o9t6qtYI3tZDE2znsjJd1sg+rvE+qDmG89b2FmdTq3ZHYuMiuj0cKpOx846GScbpZcDkLMsYkYrpnVomKw7KTxB+Si7mBNyxOL+P8Fy1yTaM2ekqDOSGJAtELSRGkatExeSI72rRMywLgcokjeqQ3LJYkK0VlC9kMNoYfKtFExcRh8iuolNpc7C4aanEKSG1wJBNJxAji7br9x3Q38wTFEMyM7pbK2z6ZXROORMWoxNBPAdJYpIj4PaW/NQ1lEp4eT59wm9hcWOhByIPEFedT22NDFi5L2GNdemXBgx8QJpFQ9GUSRIoUqWnEkUKNJ2RLIUKNIWRLo9FGufaXQr0VKOZrFnvoyjvBmDEvgUqY1BkCqjTEGb7ynqY05CRCzJ4K+teBZg4a+KYrWfEydRNRSgiMFqZ1SPCEHTpmJhykZ5/BTL2GMWyy4fF/sWQelsAxQeowZsS4qUYyQy9XRhJEd60Ri0IBViY9SLiDrEHmPgtKllF5rJCdKt1ExcRt0qsomTR2zoO2cMcMmISts6cdXBcZinabl353D2MaeKdphtWTk6y3dLavI6kthMEAcN2ow70bEKmLRr3ekxd7JSXG3cHb7fyrg6+Gy3d6jKe6KHqVbUWGTRYRsXRhIqPtiWmSB5kKykSOiBKWIshYp0hai6DqFzbUXR51C5dvBdHjoVipDEGQ6gLaDN4yKiqTdZtHkpa+cMBPHgnqoOTwM7lVDc+pnXTXN10VHCOXKzdLI6xyq0MVMWqnRqkNtP0jPBx+FvmrP2GMbs2RXx/glb6ywbMU16hooxW8owUY25yskYyQy9XRi0NkqxVEDF3dlp77e5MadcstKWCqMiawZOaNf0bbGSYjPd4c2mjINRJpvHURMP2jx5DPUi+sK8sR1WoUFhdT6WhiaxrWMAa1oDWtaLNa0CwAHAAJk4wtAAgDFdKGBGaBtTE28lPdxAzL4DbrBbiRYOHg4DVKayjva+OqNK5YeH5nP8NqAQCuFVY4vDNJIu4HLqV2mbRNjTSmVwSYwquRI+1qXmyUK3UlYXQbq59pdHhC5dxdEadyUiuTWJVVcqariMwjkosRq2tBJKfprcnhDsUoLdIyFbVGRxPDguxXWoLBzNRqHZIilamMR2ORVcRquQ+1yzaHq5jVO+8j3cANweWvvPuVprEEvXkb073TcvTgfL1TA0z0SIwVYvfVcFWeF6nBlIbc5WSMmNlytgoVmMv7LR3/ACTOmXLMb5YSNJ0f9G9TiBbLLvU9Le5mIs+UcoQdfxHId5Fk9GGTnW6pQ4XLPorCMLgpomU9NG2KNgs1rfeSTmSeJOZW6WDmSk5PLJiCAQAIAEAcj202YdRSGop2/wBWebua3+7PJ0twjJ0OgvbLK/G1+j/uQN4S3cPqQqKqBAzXOrva4ZLRbQyp2N5XBLjkWnfEYJDZFlK0nArrEpZYWSPTIkLbC6QxLOubOW58GkVkrqqqCtCsZrrM7iuLNYDc+A4ldGjTuQ3mNazIyNbVukN3acGrr11xgsIxlVqtSt0IPb8v36/kRHP7lskLPR6iPWD+Qgv8FbBK093+D+Q2ZgOKttbNo6e30PH4hwYM+Z4eSlU+bGIwlEmUbd1g71hY90jp0R2QAyIwaZPRIjaVbFiRV2lWw30YM2JL1ODNiC5TgozpHRLsfS1TZKyri67cl6qBkmcNw1pc4s0fm62dxkcrro6WtKOTj9oXNSUU/I7I0AZDIDIAaAJs5h6gAQAIAEACAEyxtcHNcA5rgWua4Atc0ixBB1CAOa7R7AyRF02HdtmrqVxs5nPqnHUfdPfY6BcvV9nRs8UOGbxt8pGZp8Ts4xyB0b25PjkBY9p5FpzC4dldtTw0abc9C0hrhzWffyXUjaSmVg5qktSydoo1g5rCVsmWUGNSYgOax2Sl1NY1srK3GWtBJcAt6tK30QxGCj1M3X4852UY/M7IeQXSq0aj7Q3TTbd/Sjx6vhf7/IpZHkm7iXHmU6kksI7Om7Lqre+zxy9/T8l9Rh7lokdCchh5V0hWchl5V0KzYy5XQtIcpIblVslhEQryywnktkl4x8xmTxwRusWu0z3HokUbSMixIo2kZPesUYKNgXowZtiS4mwaC4kgNaBdznHIADiSrRjl4MpSSPprY7BvQ6OnpjbeYy8pGYMryXyZ8RvONu6y6sI7Ukecus7ybl6l0rGYIAEACABAAgAQAIArsXwKlqhu1MLJbCzXEWkaPuvbZzfIhVlCMliSySpNdDE7QdHAZFLLQSz77Wl8dPIWyMeRnuNNg4E5gXJzskbezqpdFg3hc8pSOYQYzMQCC3zuFyJaWtPDOmtNJ9Gh04pOeLR7Sq/Z6jWOjn6ojy1UztZCPwi3vKuq610QxDQ59qXyX/pGc0am5PN2ZWqY/To6YPO3L9/I29yskdOJHe5aJGrY04q6F5MZcVdC8mNuVkYSEhl1OSihkmx2YLrB+Jl+IrJDlmW8Yi7mN9YrbSu89EijaRuFiRVwG4UJFGCjkBkRgzlI6R0NbJmeYV87foYXfQBwylqB6w5hmt/tW+yU1RX+JnM1mo42I7omjmAgAQAIAEACABAAgAQAIAEAcW6U9lDTSuroG/QSuvMBpBO46/heT5OJ5gJDV0Z8aOtodT+CX5GJa5cxo7UWekqBiI28qyGYEd5WiGIjLldFpMacroxkNkKyMWhO6pyU2ZHmNAzKo3ngJNRRDqqm+i2hXg51t+WRS9a4F3Nnm8jAbmKDlGCd4sPUYI3ig5Rgq5mz6Pdg58RkEjw6KlafpJtDJbWOK+p4F2gz45LWuvPLFL9SocLqfRVBRRwxxwwsEcbGhkbG6NaPj4po5TbbyyQggEACABAAgAQAIAEACABAAgBqqp2SMfHI0PY9pY9jhdrmuFiCOVkEp45Rw/bnYSWhLp6cOmpdSfrS03dJzZyf+rmedqNL+KJ2dJrlLwz6mQEgKQawdmE8njihDUZDT1dDEZDLgro0zkbIVijR5uqclGkjxzg3MoSbF7bowRAqaonTRMQrwcu29zI11oL5BSAoBQGT0BQGS1wTZ6sqzu0lPLNnYuaLRg/ekdZrfMqVFszndCPVnWtkOhtjC2XE3iUjMU0RIiB/xH5F/gABlqQto1JdRG3Vt8R4OsQQsY1rI2tY1oDWMYA1jWgWAaBkAOS1Em88scQAIAEACABAAgAQAIAEACABAAgAQAEIA59tR0W005dLRu9DkOZY0b1M898fqH8Nh3FYWaeEx2jXWVcPlHNcY2LxOmv1lM+Vo/tab+sMPfZo3wPFoSc9JNdOTr09p1S6vBm31DQSHHdIyLXdlwPIg6LHu5LyOjDVQfRiHTN5hCg/Q0+0x9RsTtJDW3e45BrRvOJ7gNVdVSfkZz10I+ZeYbsditR+xo5mC4G/OPR22PH6SxI8AVvHTPzOdd2pHyZsMJ6EZXWdXVbW846Zpef8x9rfpKZjSkcyzXuT4XzNVSdDuEMFnsnmP2pJXA/9vdHuV9iFnqrH5j56JMF/4d47xNNcfxKdiI+0WeoxJ0O4QdGTt72yuJ991HdxJ+02eouLohwcWvFM/wDFNIL9/ZIR3cSPtNnqXOHbBYTDbq6KC4NwZWmdwPO8hcVKikUds31ZomMAADQABkABYDyVjMUgAQAIAEACABAAgAQAIAEACABAAgAQAIAEACABAGO6Q/2f5f8AUqyNK+pw+m/3jJ5/+sLL8Q8/6SO67Bfsx+AfBi2QhPqaxSUBAAgAQAIAEACABAAgAQAIAEACABAAgAQB/9k=)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://www.dolfin-adjoint.org/en/latest/documentation/klein/klein.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image6](https://ecmiindmath.files.wordpress.com/2015/11/fenics-2.png?w=610)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://fenicsproject.org/pub/graphics/fenics-poster-2015.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image7](https://i2.wp.com/legato-team.eu/wp-content/uploads/2017/02/image_RF.png?resize=300%2C174)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image8](https://i.ytimg.com/vi/TZtX2fDIsHo/maxresdefault.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://www.youtube.com/watch?v=TZtX2fDIsHo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite Element, Volume, Difference... Infinite Finites?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finite element is the mathematician's approach. But that doesn't mean it's entirely bad."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image9](https://lh3.googleusercontent.com/Yjwlav26_7a3pmhi6elB8PDLgwJgfNBwaGmttX8CA513HoySlY7hIqzbLm9qB3GGOTCkce3d7xPCciLLyCog=w1366-h635)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://femtable.org/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No, please stay."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At a high level, if we have a continuous formula that applies over some domain, for example:\n",
    "* Maxwell's equations around a coax cable\n",
    "* Elasticity equations in a sheet of ice\n",
    "* Navier-Stokes equations in a water pipe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, multiple equations linked together,\n",
    "* potential flow and elasticity: hydroelasticity\n",
    "* fluid flow and heat equation: re-entry\n",
    "* magnetohydronamics: most existing equations rolled together for star mechanics and general masochism"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerically, there are different ways of approaching each of these. For example, conventional wisdom says \"Finite Volume\" for fluids, \"Finite Element\" for solids. Historically, that was true. However, it turns out, with a modern approach, Finite Volume and Finite Difference are special cases of Finite Element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That then means you need to pick the right finite element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...but that involves a lot of maths, so just ask someone."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most common finite elements are _Lagrange_ finite elements. They're also the most intuitive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![lagrange-o1-1d](./triangle.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a 2D, order 1 Lagrange finite element. Or, as it's also known, a triangle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Essentially, you define a 2D function by making a patch-work of these over your domain(s) and assigning values at each corner of a triangle. So that tells you the value at a series of points - the value at any point inside a triangle (or on an edge) is based on the values at the corners like so: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![linear-triangle](triangle-function.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Put these together and require that overlapping corners on the neighbouring elements match, and you get a continuous, piecewise linear function:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![cont-lin-piecewise](triangle-tesselation.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Different types of finite elements work in different ways. They may be different shapes, involve higher order polynomials instead of being linear (or more exotic things), they may specify derivatives, second derivatives, or vectors of function values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, each element (triangle) is defined now by the location of its corners, and the function values at each corner. Correspondingly our whole surface is defined by a finite set of numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![approximation](triangle-approximation.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That mesh of elements can also be made finer and finer, to better approximate any 2D function you want:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![fin-ele.png](approximation.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting all of those corner values, over all the mesh, into a single large vector, gives us a way of representing the whole function as a string of numbers, which a computer can process. Finite Element theory also gives us ways of turning differential operators into matrices, so we can approximate a partial differential equation (PDE) of continuous functions with a matrix problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\begin{eqnarray}\n",
    "\\frac{d}{dt} \\, f(x) &= 2 f(x) \\\\\n",
    "\\mathrm{approximated...} & \\\\\n",
    "\\frac{d}{dt} \\tilde{f} &= 2 \\tilde{f}(x) \\,  \\\\\n",
    "\\implies M \\, \\bf{F} &= 2 \\bf{F} \\\\\n",
    "\\textrm{Solve with a computer...} & \\\\\n",
    "\\bf{F} &= \\begin{bmatrix} 1 \\\\ 3 \\\\ 5 \\\\ \\vdots \\\\ \\end{bmatrix} \\\\\n",
    "\\textrm{Gives} & \\\\\n",
    "\\tilde{f}(x) &= \\textrm{FEalgorithm}(\\bf{F}, x)\n",
    "\\end{eqnarray} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make $\\tilde{f}$ a better and better approximation of $f$, but this normally involves making the matrix problem the computer solves bigger and bigger, so it's often a trade off of accuracy vs time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational Form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finite element actually works by integration. To solve a differential problem, we solve a mathematically equivalent \"variational problem\". This is important because this is the actually maths that the computer will discretize, solve and present back to you."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why? Because\n",
    "\n",
    "1. a differential equation involves finding the slope of a function, and that doesn't really exist where the triangles in $\\tilde{f}$ meet\n",
    "2. we want to get an approximate solution, so we average error over an area (like a cell), rather than trying to get all points exact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, I'm not going to dwell on that now, but it may help explain why the equations we enter are a bit different to the equations we ultimately want to solve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, if we want to solve Poisson's equation (which comes up in electrostatics, and a few other places). It looks like (in 2D):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{d^2}{d x^2} f + \\frac{d^2}{d y^2} f = g\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for some _known_ function $g$ and some _unknown_ function $f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make it more intuitive, this is saying that the \"divergence of $f$\" is $g$, or that $g$ is the flux entering/leaving each point of $f$, or that where-ever $g$ is big, $f$ may have a pointy basin (to which things flood), and if it's very negative, $f$ may have a pointy peak (from which things leave). If you think about this for a while, you might realise that, if we are trying to define a function only by an equation of its derivatives, then we also need some more information..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are boundary conditions. It turns out that, if we say $f$ is defined on a square, know it has (some vaguely believable) flux $g$ and we can say it's flat at the square's boundary (zero normal derivative), then we have a properly mathematically defined problem. We can find $f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is exciting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from fenics import *\n",
    "\n",
    "%matplotlib notebook\n",
    "\n",
    "# Create mesh and define function space\n",
    "mesh = UnitSquareMesh(3, 3)\n",
    "V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "\n",
    "# Define variational problem\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "f = Expression(\"sin(2 * pi * x[0]) + cos(2 * pi * x[1])\", degree=1)\n",
    "a = dot(grad(u), grad(v))*dx\n",
    "L = f*v*dx\n",
    "\n",
    "# Compute solution\n",
    "u = Function(V)\n",
    "solve(a == L, u)\n",
    "plot(u)\n",
    "\n",
    "u.rename(\"u\", \"Solution Variable\")\n",
    "outfile = File('u3.pvd')\n",
    "outfile << u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although, 2D is less fun than 3D..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "from fenics import *\n",
    "from mshr import Sphere, generate_mesh\n",
    "\n",
    "centre = Point(0.0, 0.0, 0.0)\n",
    "sphere = Sphere(centre, 1)\n",
    "mesh = generate_mesh(sphere, 10)\n",
    "plot(mesh)\n",
    "\n",
    "f = File(\"f.pvd\")\n",
    "f << mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "\n",
    "# Define variational problem\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "f = Constant(9.81)\n",
    "a = dot(grad(u), grad(v))*dx\n",
    "L = f*v*dx\n",
    "\n",
    "class SphereBoundary(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return on_boundary\n",
    "    \n",
    "bc = DirichletBC(V, Constant(0.0), SphereBoundary())\n",
    "\n",
    "# Compute solution\n",
    "u = Function(V)\n",
    "\n",
    "solve(a == L, u, bc)\n",
    "u.rename(\"u\", \"Solution Variable\")\n",
    "outfile = File('u.pvd')\n",
    "outfile << u\n",
    "plot(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Credit: FEniCS Project official demos\n",
    "\n",
    "from fenics import *\n",
    "\n",
    "# Load mesh from file\n",
    "mesh = Mesh(\"../demo/documented/navier-stokes/lshape.xml.gz\")\n",
    "\n",
    "# Define function spaces (P2-P1)\n",
    "V = VectorFunctionSpace(mesh, \"Lagrange\", 2)\n",
    "Q = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "\n",
    "# Define trial and test functions\n",
    "u = TrialFunction(V)\n",
    "p = TrialFunction(Q)\n",
    "v = TestFunction(V)\n",
    "q = TestFunction(Q)\n",
    "\n",
    "# Set parameter values\n",
    "dt = 0.01\n",
    "T = 3\n",
    "nu = 0.01\n",
    "\n",
    "# Define time-dependent pressure boundary condition\n",
    "p_in = Expression(\"sin(3.0*t)\", t=0.0, degree=1)\n",
    "\n",
    "# Define boundary conditions\n",
    "noslip  = DirichletBC(V, (0, 0),\n",
    "                      \"on_boundary && \\\n",
    "                       (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \\\n",
    "                       (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))\")\n",
    "inflow  = DirichletBC(Q, p_in, \"x[1] > 1.0 - DOLFIN_EPS\")\n",
    "outflow = DirichletBC(Q, 0, \"x[0] > 1.0 - DOLFIN_EPS\")\n",
    "bcu = [noslip]\n",
    "bcp = [inflow, outflow]\n",
    "\n",
    "# Create functions\n",
    "u0 = Function(V)\n",
    "u1 = Function(V)\n",
    "p1 = Function(Q)\n",
    "\n",
    "# Define coefficients\n",
    "k = Constant(dt)\n",
    "f = Constant((0, 0))\n",
    "\n",
    "# Tentative velocity step\n",
    "F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \\\n",
    "     nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx\n",
    "a1 = lhs(F1)\n",
    "L1 = rhs(F1)\n",
    "\n",
    "# Pressure update\n",
    "a2 = inner(grad(p), grad(q))*dx\n",
    "L2 = -(1/k)*div(u1)*q*dx\n",
    "\n",
    "# Velocity update\n",
    "a3 = inner(u, v)*dx\n",
    "L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx\n",
    "\n",
    "# Assemble matrices\n",
    "A1 = assemble(a1)\n",
    "A2 = assemble(a2)\n",
    "A3 = assemble(a3)\n",
    "\n",
    "# Use amg preconditioner if available\n",
    "prec = \"amg\" if has_krylov_solver_preconditioner(\"amg\") else \"default\"\n",
    "\n",
    "# Use nonzero guesses - essential for CG with non-symmetric BC\n",
    "parameters['krylov_solver']['nonzero_initial_guess'] = True\n",
    "\n",
    "# Create files for storing solution\n",
    "ufile = File(\"results/velocity.pvd\")\n",
    "pfile = File(\"results/pressure.pvd\")\n",
    "\n",
    "# Time-stepping\n",
    "t = dt\n",
    "while t < T + DOLFIN_EPS:\n",
    "\n",
    "    # Update pressure boundary condition\n",
    "    p_in.t = t\n",
    "\n",
    "    # Compute tentative velocity step\n",
    "    begin(\"Computing tentative velocity\")\n",
    "    b1 = assemble(L1)\n",
    "    [bc.apply(A1, b1) for bc in bcu]\n",
    "    solve(A1, u1.vector(), b1, \"bicgstab\", \"default\")\n",
    "    end()\n",
    "\n",
    "    # Pressure correction\n",
    "    begin(\"Computing pressure correction\")\n",
    "    b2 = assemble(L2)\n",
    "    [bc.apply(A2, b2) for bc in bcp]\n",
    "    [bc.apply(p1.vector()) for bc in bcp]\n",
    "    solve(A2, p1.vector(), b2, \"bicgstab\", prec)\n",
    "    end()\n",
    "\n",
    "    # Velocity correction\n",
    "    begin(\"Computing velocity correction\")\n",
    "    b3 = assemble(L3)\n",
    "    [bc.apply(A3, b3) for bc in bcu]\n",
    "    solve(A3, u1.vector(), b3, \"bicgstab\", \"default\")\n",
    "    end()\n",
    "\n",
    "    # Plot solution\n",
    "    # plot(p1, title=\"Pressure\", rescale=True)\n",
    "    # plot(u1, title=\"Velocity\", rescale=True)\n",
    "\n",
    "    # Save to file\n",
    "    ufile << u1\n",
    "    pfile << p1\n",
    "\n",
    "    # Move to next time step\n",
    "    u0.assign(u1)\n",
    "    t += dt\n",
    "    print(\"t =\", t)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While this may seem long, it's because this contains _everything_ specific to the physics. The best comparison with other tools, is looking at the source files for the FE tool's Navier-Stokes solver: here's Elmer's closest equivalent at 1652 lines... https://raw.githubusercontent.com/tehnick/elmerfem/master/fem/src/NavierStokes.src\n",
    "\n",
    "FEniCS is able to bundle all the physics into a couple dozen lines (and Navier-Stokes is definitely at the complex end), because it has a deep understanding of the mathematical theory - this allows it provide just as much abstraction needed to be generic, but no more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Isn't this way more complicated?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not necessarily - firstly this is doing something different to the traditional approach (e.g. Elmer, COMSOL, Abaqus...), where the software recognises the general type of equation and picks a standard set of configurations, including elements, to solve it (potentially with some parameters). NASTRAN, PATRAN (engineering industry standard) let you pick an element, but don't really give you enough information about it... (apparently the numerics you're using are their trade secret)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting the settings wrong can lead to unexpected problems..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![tricell](tricell.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![sleipner](http://www-users.math.umn.edu/~arnold/disasters/sleipner.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://www-users.math.umn.edu/~arnold/disasters/sleipner.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is why having access to the actual code, even if you need to get help reading it, can avoid months of battling configuration, undiagnosable weird numerical problems, or expensive fees to licensed consultants, or worse (^^^), when the answer is in the maths on your laptop."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ALSO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://www.comsol.com/model/microwave-heating-of-a-cancer-tumor-30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having researched this (and built a comparable model in FEniCS), the research papers which used this often got the maths completely wrong, or misread the manual and copied the text in a nonsensical way... not ideal for reproducability... which meant cross-correlating half-a-dozen papers to understand the physics, and separate out the mistakes. However, our FEniCS and Elmer models were able to reproduce COMSOL results, and the process highlighted where the limitations were, and the things COMSOL had left out of the manual..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How do I get a black box?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power of Python..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn an official demo into a solver\n",
    "def solve_navier_stokes(mesh, bcu, bcp, dt, T, nu, ufile, pfile):\n",
    "    # Define function spaces (P2-P1)\n",
    "    V = VectorFunctionSpace(mesh, \"Lagrange\", 2)\n",
    "    Q = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "\n",
    "    # Define trial and test functions\n",
    "    u = TrialFunction(V)\n",
    "    p = TrialFunction(Q)\n",
    "    v = TestFunction(V)\n",
    "    q = TestFunction(Q)\n",
    "\n",
    "    # Define time-dependent pressure boundary condition\n",
    "    p_in = Expression(\"sin(3.0*t)\", t=0.0, degree=1)\n",
    "\n",
    "\n",
    "    # Create functions\n",
    "    u0 = Function(V)\n",
    "    u1 = Function(V)\n",
    "    p1 = Function(Q)\n",
    "\n",
    "    # Define coefficients\n",
    "    k = Constant(dt)\n",
    "    f = Constant((0, 0))\n",
    "\n",
    "    # Tentative velocity step\n",
    "    F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \\\n",
    "         nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx\n",
    "    a1 = lhs(F1)\n",
    "    L1 = rhs(F1)\n",
    "\n",
    "    # Pressure update\n",
    "    a2 = inner(grad(p), grad(q))*dx\n",
    "    L2 = -(1/k)*div(u1)*q*dx\n",
    "\n",
    "    # Velocity update\n",
    "    a3 = inner(u, v)*dx\n",
    "    L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx\n",
    "\n",
    "    # Assemble matrices\n",
    "    A1 = assemble(a1)\n",
    "    A2 = assemble(a2)\n",
    "    A3 = assemble(a3)\n",
    "\n",
    "    # Use amg preconditioner if available\n",
    "    prec = \"amg\" if has_krylov_solver_preconditioner(\"amg\") else \"default\"\n",
    "\n",
    "    # Use nonzero guesses - essential for CG with non-symmetric BC\n",
    "    parameters['krylov_solver']['nonzero_initial_guess'] = True\n",
    "\n",
    "    # Time-stepping\n",
    "    t = dt\n",
    "    while t < T + DOLFIN_EPS:\n",
    "\n",
    "        # Update pressure boundary condition\n",
    "        p_in.t = t\n",
    "\n",
    "        # Compute tentative velocity step\n",
    "        begin(\"Computing tentative velocity\")\n",
    "        b1 = assemble(L1)\n",
    "        [bc.apply(A1, b1) for bc in bcu]\n",
    "        solve(A1, u1.vector(), b1, \"bicgstab\", \"default\")\n",
    "        end()\n",
    "\n",
    "        # Pressure correction\n",
    "        begin(\"Computing pressure correction\")\n",
    "        b2 = assemble(L2)\n",
    "        [bc.apply(A2, b2) for bc in bcp]\n",
    "        [bc.apply(p1.vector()) for bc in bcp]\n",
    "        solve(A2, p1.vector(), b2, \"bicgstab\", prec)\n",
    "        end()\n",
    "\n",
    "        # Velocity correction\n",
    "        begin(\"Computing velocity correction\")\n",
    "        b3 = assemble(L3)\n",
    "        [bc.apply(A3, b3) for bc in bcu]\n",
    "        solve(A3, u1.vector(), b3, \"bicgstab\", \"default\")\n",
    "        end()\n",
    "\n",
    "        # Save to file\n",
    "        ufile << u1\n",
    "        pfile << p1\n",
    "\n",
    "        # Move to next time step\n",
    "        u0.assign(u1)\n",
    "        t += dt\n",
    "        print(\"t =\", t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In real-life we would normally have the above solver in one Python file, and then be able to run the code below in a much shorter, tidier, less confusing script:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is now my Python script\n",
    "\n",
    "# Load mesh from file\n",
    "mesh = Mesh(\"../demo/documented/navier-stokes/lshape.xml.gz\")\n",
    "\n",
    "# Set parameter values\n",
    "dt = 0.01\n",
    "T = 3\n",
    "nu = 0.01\n",
    "\n",
    "# Create files for storing solution\n",
    "ufile = File(\"results/velocity.pvd\")\n",
    "pfile = File(\"results/pressure.pvd\")\n",
    "\n",
    "\n",
    "# Define boundary conditions\n",
    "noslip  = DirichletBC(V, (0, 0),\n",
    "                      \"on_boundary && \\\n",
    "                       (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \\\n",
    "                       (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))\")\n",
    "inflow  = DirichletBC(Q, p_in, \"x[1] > 1.0 - DOLFIN_EPS\")\n",
    "outflow = DirichletBC(Q, 0, \"x[0] > 1.0 - DOLFIN_EPS\")\n",
    "bcu = [noslip]\n",
    "bcp = [inflow, outflow]\n",
    "\n",
    "solve_navier_stokes(mesh, bcu, bcp, dt, T, nu, ufile, pfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FEniCS demos...\n",
    "\n",
    "import os\n",
    "\n",
    "demos = os.listdir('../demo/documented') + os.listdir('../demo/undocumented')\n",
    "from IPython.display import display, Markdown\n",
    "display(Markdown('* ' + '\\n* '.join(demos)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find Python code [here](/tree/demo) and explanatory documentation [here](https://fenics.readthedocs.io/projects/dolfin/en/stable/demos.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The documented demos provide a walk-through of the finite element theory showing that the solver is mathematically sensible, and why it is the way it is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next Time\n",
    "\n",
    "* Which solvers/physics?\n",
    "* When?\n",
    "* Any questions for devs?\n",
    "* Using a development environment\n",
    "* Take a look at Garth Well's materials on Microsoft Azure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Challenge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In remaining time, create a new notebook and try to get one of the demos to run."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stretch Goal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get that working and are feeling adventurous, make a function, like the black box Navier-Stokes approach above, so you can pass any desired parameters as arguments, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "from fenics import *\n",
    "\n",
    "param1 = 2.3\n",
    "param2 = 84.1\n",
    "bc = DirichetBC(...)\n",
    "mesh = Mesh('....') # or maybe a sphere, box, etc.\n",
    "outfile = File('outvar.pvd')\n",
    "\n",
    "solve_blackboxphysics(param1, param2, bc, outfile)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}