nonlinear-kirchhoff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The non-linear Kirchhoff model\n",
    "\n",
    "## Model and energy\n",
    "\n",
    "## The isometry constraint\n",
    "\n",
    "## Discretisation\n",
    "\n",
    "## Linearisation of the nodal isometry constraint\n",
    "\n",
    "> For this, it is important ot realize that for the empleyed finite element space $W_h$, the nodal values of the discrete deformation ($u_h(z) : z \\in \\mathcal{N}_h$) are mutually independent variables in the minimization problem.\n",
    "\n",
    "## Discrete $H^2$ gradient flow\n",
    "\n",
    "# Implementation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Besides the work already done for [the linear case](linear-kirchhoff.ipynb), we need to compute the matrix $B_{k-1}$ which enforces the nodal isometry constraint for the solution $d_t y^k$ at step $k$.\n",
    "\n",
    "The system to solve is (Bartels, 2013):\n",
    "\n",
    "$$\\left(\\begin{array}{cc}\n",
    "  (1 + \\alpha \\tau) M^{\\top} A^{(2)} M & B^{\\top}_{k - 1}\\\\\n",
    "  B_{k - 1} & I\n",
    "\\end{array}\\right)  \\left(\\begin{array}{c}\n",
    "  d_t Y^k\\\\\n",
    "  \\Lambda\n",
    "\\end{array}\\right) = \\left(\\begin{array}{c}\n",
    "  - \\alpha M^{\\top} A^{(2)} MY^{k - 1} + \\tau F\\\\\n",
    "  0\n",
    "\\end{array}\\right)$$\n",
    "\n",
    "The local tensor matrix requires then the discrete gradient matrix $M$, the local tensor for the form $(\\nabla u, \\nabla v)$ in $P_2$ and the constraints matrix $B_{k-1}$. Note also that we now have vector valued functions and this will probably require tweaking some of the previous code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete isometry constraint\n",
    "\n",
    "\\\\[ \\left( \\begin{array}{ccccccccc}\n",
    "     \\quad & 2 y^{k - 1}_{i, 1} &  & \\quad & 2 y^{k - 1}_{i, 1} &  & \\quad & 2\n",
    "     y^{k - 1}_{i, 1} & \\\\\n",
    "     & y^{k - 1}_{i, 2} & y^{k - 1}_{i, 1} &  & y^{k - 1}_{i, 2} & y^{k -\n",
    "     1}_{i, 1} &  & y^{k - 1}_{i, 2} & y^{k - 1}_{i, 1}\\\\\n",
    "     & y^{k - 1}_{i, 2} & y^{k - 1}_{i, 1} &  & y^{k - 1}_{i, 2} & y^{k -\n",
    "     1}_{i, 1} &  & y^{k - 1}_{i, 2} & y^{k - 1}_{i, 1}\\\\\n",
    "     &  & 2 y^{k - 1}_{i, 2} &  &  & 2 y^{k - 1}_{i, 2} &  &  & 2 y^{k -\n",
    "     1}_{i, 2}\n",
    "   \\end{array}\\right)  \\left(\\begin{array}{c}\n",
    "     d_t y^k_i\\\\\n",
    "     d_t y^k_{i, 1}\\\\\n",
    "     d_t y^k_{i, 2}\\\\\n",
    "     d_t y^k_i\\\\\n",
    "     d_t y^k_{i, 1}\\\\\n",
    "     d_t y^k_{i, 2}\\\\\n",
    "     d_t y^k_i\\\\\n",
    "     d_t y^k_{i, 1}\\\\\n",
    "     d_t y^k_{i, 2}\n",
    "   \\end{array}\\right) \\\\]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "np.set_printoptions(precision=4, linewidth=130, threshold=5000, suppress=True)\n",
    "\n",
    "import nbimporter\n",
    "from interpolation import interpolate\n",
    "\n",
    "import matplotlib.pyplot as pl\n",
    "%matplotlib inline\n",
    "\n",
    "#info(parameters, True)\n",
    "\n",
    "def bitmap(A, rtol=1e-5, atol=1e-8, figsize=(10,10), cmap='binary', **kwargs):\n",
    "    \"\"\"Draw the number of non zeros of A.\"\"\"\n",
    "    if np.sum(A.shape) > 4000:\n",
    "        print(\"Matrix is too big (%d x %d)\" % A.shape)\n",
    "    pl.figure(figsize=figsize)\n",
    "    bmap = 1 - np.isclose(A, 0, rtol, atol).astype(np.int)\n",
    "    pl.imshow(bmap, cmap=cmap, **kwargs)\n",
    "    print(\"%.2f%% non zeros\" % (100*bmap.sum()/np.product(bmap.shape)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#path = \"/home/fenics/local/src/nonlinear-kirchhoff/build\"\n",
    "path = \"/home/fenics/.netbeans/remote/127.0.0.1/hpel620-Windows-x86_64/C/Users/Miguel/Devel/nonlinear-kirchhoff\"\n",
    "A = np.loadtxt(path + \"/A.data\")\n",
    "Bk = np.loadtxt(path + \"/Bk.data\")\n",
    "#D = np.loadtxt(path + \"/D.data\")\n",
    "D = np.identity(Bk.shape[0])\n",
    "Fk = np.loadtxt(path + \"/Fk.data\")\n",
    "Y0 = np.loadtxt(path + \"/y0.data\")\n",
    "Yk = np.loadtxt(path + \"/yk.data\")\n",
    "dtY = np.loadtxt(path + \"/dtY_L.data\")\n",
    "#P26 = np.loadtxt(path + \"/P26.data\")\n",
    "# Cut out the bottom extra chunk in the solution and RHS\n",
    "dtY = dtY[:-Bk.shape[0]]\n",
    "Fk = Fk[:-Bk.shape[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.shape, Bk.shape, Fk.shape, Y0.shape, Yk.shape, dtY.shape #, P26.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = np.loadtxt(path+\"/energy.data\")\n",
    "pl.figure(figsize=(16,4))\n",
    "_ = pl.plot(E)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring the stiffness matrix\n",
    "\n",
    "First a nice picture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bitmap(A, figsize=(5,5), interpolation='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The condition nunmber is **awful**, something is definitely wrong:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.cond(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bitmap(Bk[176:180,190:200], figsize=(10,4), interpolation='none')\n",
    "print(Bk[176:180,190:200])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = VectorFunctionSpace(UnitSquareMesh(10,10, \"right\"), \"DKT\", degree=3, dim=3)\n",
    "v2d = vertex_to_dof_map(W)\n",
    "d2v = (dof_to_vertex_map(W)/W.dofmap().num_entity_dofs(0)).astype(np.int)\n",
    "print((d2v.min(),d2v.max()))\n",
    "\n",
    "#idx = 44\n",
    "#for sub in range(3):\n",
    "#    print(v2d[9*idx + 3*sub + 0], v2d[9*idx + 3*sub + 1], v2d[9*idx + 3*sub + 2])\n",
    "\n",
    "#A[189:197,189:197], Fk[189:197], dtY[189:197]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reread full RHS and solve here\n",
    "Fk = np.loadtxt(path + \"/Fk.data\")\n",
    "dtY = np.loadtxt(path + \"/dtY_L.data\")\n",
    "\n",
    "assert np.all((D - np.identity(D.shape[0])) == 0), \"Padding matrix must be the identity\"\n",
    "M = np.block([[A, Bk.T],[Bk, D]])\n",
    "\n",
    "# Compare np's solution with ours\n",
    "dZ = np.linalg.solve(M,Fk)\n",
    "np.allclose(dZ, dtY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the components of the solution\n",
    "\n",
    "A quick visualization skipping Paraview...\n",
    "\n",
    "Recall that the $W^3$ space has 27 dofs in 3 groups of 9, one per subspace. Inside each group the dofs are: [evaluation, eval of dx, eval of dy] at each of the 3 vertices in turn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = VectorFunctionSpace(RectangleMesh(Point(-2,0), Point(2,1), 10, 10, \"right\"),\n",
    "                        \"DKT\", degree=3, dim=3)\n",
    "y0 = Function(W)\n",
    "yk = Function(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y0 = np.loadtxt(path + \"/y0.data\")\n",
    "Yk = np.loadtxt(path + \"/yk.data\")\n",
    "y0.vector().set_local(Y0)\n",
    "_ = yk.vector().set_local(Yk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_xx = np.array(sorted(W.mesh().coordinates()[:,0]))\n",
    "#_xx = np.arange(-2,-1.8,0.001)\n",
    "pl.figure(figsize=(12,6))\n",
    "for _y in np.linspace(0,1,6):\n",
    "    pl.plot(_xx, [y0(x,_y)[2] - 0 for x in _xx], label='$u_3^0(%.1f)$' % _y)\n",
    "for _y in np.linspace(0,1,6):\n",
    "    pl.plot(_xx, [yk(x,_y)[2] - 0. for x in _xx], label='$u_3^k(%.1f)$' % _y)\n",
    "_ = pl.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For later use, let us extract the dofs for evaluations and both partial derivatives and arrange them into arrays with three components, one per dimension of the range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dofs import plot_dofs, plot_field_at_dofs, extract_dofs_with_mask\n",
    "\n",
    "_f = extract_dofs_with_mask(W, np.arange(0,27,3))\n",
    "f = [_f[::3], _f[1::3], _f[2::3]]\n",
    "_dfdx = extract_dofs_with_mask(W, np.arange(1,27,3))\n",
    "dfdx = [_dfdx[::3], _dfdx[1::3], _dfdx[2::3]]\n",
    "_dfdy = extract_dofs_with_mask(W, np.arange(2,27,3))\n",
    "dfdy = [_dfdy[::3], _dfdy[1::3], _dfdy[2::3]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "_sf, _sdfdx, _sdfdy = set(_f), set(_dfdx), set(_dfdy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial condition has to have $\\nabla^T y^0 \\nabla y^0 = I_2$ at all nodes. If the following is not `(True, True, True)`, then we have a(nother) problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nvert = W.mesh().num_vertices()\n",
    "dxdx = np.dot(Y0[_dfdx], Y0[_dfdx]) / nvert\n",
    "dxdy = np.dot(Y0[_dfdx], Y0[_dfdy]) / nvert\n",
    "dydy = np.dot(Y0[_dfdy], Y0[_dfdy]) / nvert\n",
    "np.isclose(dxdx, 1), np.isclose(dxdy, 0), np.isclose(dydy,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dxdx, dxdy, dydy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How is it possible that dxdy is not 0? $\\partial_x y_0$ and $\\partial_y y_0$ should be orthogonal!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All iterates $y^k$ fulfill the boundary condition since updates are always made with $d_tY$ with zero value and derivative at the boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.allclose(Yk[dfdx[0]], 1), np.allclose(Yk[dfdy[1]], 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was bad! Let's see where it failed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.where(np.logical_not(np.isclose(Y0[dfdx[0]], 1.)))[0].size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl.figure(figsize=(16,6))\n",
    "pl.subplot(2,2,1)\n",
    "plot_dofs(W, np.where(np.logical_not(np.isclose(Y0[dfdy[0]], 0.)))[0])\n",
    "pl.subplot(2,2,2)\n",
    "plot_dofs(W, np.where(np.logical_not(np.isclose(Y0[dfdy[1]], 1.)))[0])\n",
    "pl.subplot(2,2,3)\n",
    "plot_dofs(W, np.where(np.logical_not(np.isclose(Y0[dfdy[2]], 0.)))[0])\n",
    "\n",
    "#plot_dofs(W, np.where(np.logical_not(np.isclose(Yk[dfdx[0]], 1.)))[0])\n",
    "#pl.subplot(2,2,4)\n",
    "#plot_dofs(W, np.where(np.logical_not(np.isclose(Yk[dfdy[1]], 1.)))[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Something is clearly off, let's plot the components of the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_component(Y, W, i, ii, **kwargs):\n",
    "    y = Function(FunctionSpace(W.mesh(), \"DKT\", 3))\n",
    "    vec = np.zeros_like(y.vector().array())\n",
    "    vec[i::3] = Yk[ii[i]]/1e8\n",
    "    _ = y.vector().set_local(vec)\n",
    "\n",
    "    P = FunctionSpace(W.mesh(), \"Lagrange\", 1)\n",
    "    p = Function(P)\n",
    "    vals = y.compute_vertex_values(W.mesh())\n",
    "\n",
    "    if np.any(np.isnan(vals)):\n",
    "        print(\"%.2f%% NaNs in vertex values for component %d, unable to plot\"\n",
    "              % (100*np.sum(np.isnan(vals))/len(vals), i))\n",
    "        return\n",
    "    d2v = dof_to_vertex_map(P)\n",
    "    _ = p.vector().set_local(vals[d2v])\n",
    "\n",
    "    plot(p, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(dtY, W, 0, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl.figure(figsize=(14,4))\n",
    "pl.subplot(1,2,1)\n",
    "plot_component(Y0, W, 1, dfdx)\n",
    "pl.subplot(1,2,2)\n",
    "plot_component(Y0, W, 2, dfdy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl.figure(figsize=(18,10))\n",
    "pl.subplot(3,3,1)\n",
    "plot_component(Yk, W, 0, f, title='$y_k^0$')\n",
    "pl.subplot(3,3,2)\n",
    "plot_component(Yk, W, 1, f, title='$y_k^1$')\n",
    "pl.subplot(3,3,3)\n",
    "plot_component(Yk, W, 2, f, title='$y_k^2$')\n",
    "pl.subplot(3,3,4)\n",
    "plot_component(Yk, W, 0, dfdx, title='$\\partial_x{y_k^0}$')\n",
    "pl.subplot(3,3,5)\n",
    "plot_component(Yk, W, 1, dfdx, title='$\\partial_x{y_k^1}$')\n",
    "pl.subplot(3,3,6)\n",
    "plot_component(Yk, W, 2, dfdx, title='$\\partial_x{y_k^2}$')\n",
    "pl.subplot(3,3,7)\n",
    "plot_component(Yk, W, 0, dfdy, title='$\\partial_y{y_k^0}$')\n",
    "pl.subplot(3,3,8)\n",
    "plot_component(Yk, W, 1, dfdy, title='$\\partial_y{y_k^1}$')\n",
    "pl.subplot(3,3,9)\n",
    "plot_component(Yk, W, 2, dfdy, title='$\\partial_y{y_k^2}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting the $P_2^{3 \\times 2}$ tensor\n",
    "\n",
    "Plots of the two local tensors for a 1x1 grid. What are those spurious entries on the top right corner?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bitmap(P26, figsize=(5,5), cmap='Purples', interpolation='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```c++\n",
    "  // position in array is destination row, value is source row:\n",
    "  int permutations[] = {0,6,1,7,2,8,3,9,4,10,5,11};\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T = TensorFunctionSpace(UnitSquareMesh(1,1), \"Lagrange\", 2, shape=(3,2))\n",
    "for sub in range(6):\n",
    "    print(T.sub(sub).dofmap().cell_dofs(0))\n",
    "\n",
    "dm = T.dofmap()\n",
    "dm.cell_dofs(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This replicates the contents of the UFL file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "domain = W.mesh()\n",
    "T = TensorFunctionSpace(domain, \"Lagrange\", 2, shape=(3,2))\n",
    "#W = VectorFunctionSpace(domain, \"DKT\", 3, dim=3)\n",
    "P = VectorFunctionSpace(domain, \"Lagrange\", 3, dim=3)\n",
    "\n",
    "p = TrialFunction(T)\n",
    "q = TestFunction(T)\n",
    "p22 = inner(nabla_grad(p), nabla_grad(q))*dx\n",
    "\n",
    "u = TrialFunction(W)\n",
    "v = TestFunction(W)\n",
    "dkt = inner(u,v)*dx\n",
    "\n",
    "f = Coefficient(W)\n",
    "force = inner(f,v)*dx\n",
    "\n",
    "# Define variational problem for projection\n",
    "g = Coefficient(P)\n",
    "Pg = TrialFunction(W)\n",
    "w = TestFunction(W)\n",
    "project_lhs = inner(w, Pg)*dx\n",
    "project_rhs = inner(w, g)*dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = assemble(p22)\n",
    "A.array().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bitmap(A.array())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The local tensor on a 1x1 grid\n",
    "\n",
    "Plots of the two local tensors for a 1x1 grid. What are those spurious entries on the top right corner?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Alocbmap = (np.round(Aloc, 5) != 0).astype(np.int)\n",
    "pl.imshow(Alocbmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aloc2bmap = (np.round(Aloc2, 5) != 0).astype(np.int)\n",
    "pl.imshow(Aloc2bmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "v = 2  # m/s\n",
    "r = 0.152 / 2  # m\n",
    "A = np.pi*r**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "LN = 10 + 50 * np.log (v) + 10 * np.log (A) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "LN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "W = VectorFunctionSpace(UnitSquareMesh(1,1, \"right\"), \"DKT\", degree=3, dim=3)\n",
    "T = TensorFunctionSpace(W.mesh(), \"Lagrange\", degree=2, shape=(3,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "msh = W.mesh()\n",
    "geo = msh.geometry()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(3):\n",
    "    dm = W.sub(i).dofmap()\n",
    "    print(dm.cell_dofs(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(3):\n",
    "    dm = W.sub(i).dofmap()\n",
    "    print(dm.cell_dofs(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(3):\n",
    "    dm = T.sub(i).dofmap()\n",
    "    print(dm.cell_dofs(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(3):\n",
    "    dm = T.sub(i).dofmap()\n",
    "    print(dm.cell_dofs(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "msh = RectangleMesh(Point(-2,0), Point(2,1), 400, 20, \"crossed\")\n",
    "V = FunctionSpace(msh, \"Lagrange\", 1)\n",
    "v = TestFunction(V)\n",
    "b = assemble(v*dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b.array().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "W = VectorFunctionSpace(msh, \"DKT\", degree=3, dim=3)\n",
    "w = TestFunction(W)\n",
    "\n",
    "p = Constant((1,0,0))\n",
    "\n",
    "a = assemble(inner(v, p)*dx)\n",
    "\n",
    "#c = assemble(inner(w,p)*dx)\n",
    "\n",
    "a,b #,c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.all(a.array()==b.array())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "set(a.array()[:333:3].round(6))"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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"
  },
  "toc": {
   "nav_menu": {
    "height": "167px",
    "width": "254px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "717px",
    "left": "0px",
    "right": "1019px",
    "top": "83px",
    "width": "261px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}