{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem 1\n", "\n", "Consider the following two-dimensional diffision problem. Use the discretization below to solve for the" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "unknown diffusing concentrations at the nodes with $a = 1$ and $u_o = 100$ assuming a consistent unit system. The diffusion coefficient matrix is the identity matrix. Perform the following three tasks:\n", "\n", " 1. Solve this problem using a direct integration of the trianglur elements. (**20 points**)\n", " \n", " 2. Solve this problem by using a \"parent\" element mapping and Gauss integration on the rectangular element (i.e. ignore the triangular elements). Use a $2 \\times 2$ Gauss integration rule. (**20 points**)\n", "\n", " 3. Create effective plots to visualize your results from 1 and 2. (**5 points**)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note:** Submit a working version of your code to [Canvas](https://utexas.instructure.com/courses/1119539). Any supplemental material explaining your answer in part (b) can be turned in to me via hard copy or scanned and submitted to Canvas with your code." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**\n", "\n", "Below we will write the class and methods that will help us solve the problem." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import scipy.linalg\n", "import scipy.sparse\n", "import scipy.sparse.linalg\n", "\n", "class TwoDimFEM():\n", " \n", " def __init__(self, nodes, connect, C=np.array([[1., 0, 0], [0, 1., 0], [0, 0, 0]])):\n", " \"\"\"\n", " Initialize the model.\n", " \n", " input: nodes -> the nodal locations as a Numpy array of (x,y) pairs.\n", " input: connect -> the connectivitiy array\n", " output: the model problem object\n", " \n", " \"\"\"\n", " \n", " self.X = nodes[:,0]\n", " self.Y = nodes[:,1]\n", " #Subtract 1 to convert node numbers to be consistent with Numpy 0 indexing\n", " self.connect = (connect - 1)\n", " \n", " self.num_elem = len(self.connect)\n", " self.num_dof = len(self.X)\n", " \n", " self.Cmat = C\n", " \n", " #Allocate global tangent and r.h.s vector\n", " self.K = np.zeros((self.num_dof, self.num_dof), dtype=np.double)\n", " self.F = np.zeros(self.num_dof, dtype=np.double)\n", " \n", " \n", " def N(self, xi, eta):\n", " \"\"\"Compute linear shape functions in natural coordinates.\"\"\"\n", " \n", " return [1 / 4. - eta / 4. - xi / 4. + (eta * xi) / 4.,\n", " 1 / 4. - eta / 4. + xi / 4. - (eta * xi) / 4., \n", " 1 / 4. + eta / 4. + xi / 4. + (eta * xi) / 4.,\n", " 1 / 4. + eta / 4. - xi / 4. - (eta * xi) / 4.]\n", " \n", " \n", " def dNdxi(self, eta):\n", " \"\"\"Compute shape function derivatives with respect to xi\"\"\"\n", " \n", " return [-1 / 4. + eta / 4., \n", " 1 / 4. - eta / 4., \n", " 1 / 4. + eta / 4., \n", " -1 / 4. - eta / 4.]\n", " \n", " \n", " def dNdeta(self, xi):\n", " \"\"\"Compute shape function derivatives with respect to eta\"\"\"\n", " \n", " return [-1 / 4. + xi / 4.,\n", " -1 / 4. - xi / 4.,\n", " 1 / 4. + xi / 4.,\n", " 1 / 4. - xi / 4.]\n", " \n", " \n", " def compute_jacobian_matrix_and_inverse(self, xi, eta):\n", " \"\"\"\n", " Compute the Jacobian matrix, Det(J) and B for every element\n", " \"\"\"\n", " \n", " x = self.X\n", " y = self.Y\n", " con = self.connect\n", " \n", " #Understand we are broadcasting the dot product to every element\n", " J11 = np.dot(x[con], self.dNdxi(xi))\n", " J12 = np.dot(y[con], self.dNdxi(xi))\n", " J21 = np.dot(x[con], self.dNdeta(eta))\n", " J22 = np.dot(y[con], self.dNdeta(eta))\n", " \n", " #detJ is a vector containing the Jacobian determinate for every element\n", " self.detJ = J11 * J22 - J12 * J21\n", " \n", " self.Jinv11 = J22 / self.detJ\n", " self.Jinv12 = -J12 / self.detJ\n", " self.Jinv21 = -J21 / self.detJ\n", " self.Jinv22 = J11 / self.detJ\n", " \n", " \n", " def compute_B_matrix(self, xi, eta):\n", " \"\"\"Computes the B matrix for a given xi and eta\"\"\"\n", " \n", " #Returns detJ and Jinv components for this xi and eta\n", " self.compute_jacobian_matrix_and_inverse(xi, eta)\n", " \n", " \n", " Nmat = np.zeros((3, 4), dtype=np.double)\n", " Nmat[0,:] = self.dNdxi(xi)\n", " Nmat[1,:] = self.dNdeta(eta)\n", " Nmat[2,:] = self.N(xi, eta)\n", " \n", " zero = np.zeros(len(self.detJ))\n", " one = np.ones(len(self.detJ))\n", " \n", " Jmat = np.array([[self.Jinv11, self.Jinv12, zero],\n", " [self.Jinv21, self.Jinv22, zero],\n", " [ zero, zero, one]])\n", " \n", " #B = J * N\n", " return np.einsum('jk...,kl', Jmat, Nmat)\n", " \n", " \n", " def compute_stiffness_integrand(self, xi, eta):\n", " \"\"\"Computes the integrand of the stiffness matrix for this xi and eta\"\"\"\n", " \n", " Cmat = self.Cmat\n", " \n", " #Returns B\n", " self.Bmat = self.compute_B_matrix(xi, eta)\n", " \n", " #ke_{il} = B_{ji} C_{jk} B_{kl} \\det(J)\n", " return np.einsum('...ji,jk,...kl,...',self.Bmat,Cmat,self.Bmat,self.detJ)\n", " \n", " \n", " def integrate_element_matrices(self):\n", " \"\"\"Integrate stiffness matrix with Gauss integration\"\"\"\n", " \n", " #Use 2 x 2 Gauss integration\n", " wts = [1., 1.]\n", " pts = [-np.sqrt(1 / 3.), np.sqrt(1 / 3.)]\n", " \n", " ke = np.zeros((self.num_elem, 4, 4))\n", " \n", " for i in range(2):\n", " for j in range(2):\n", " ke += wts[i] * wts[j] * self.compute_stiffness_integrand(pts[i], pts[j])\n", " \n", " return ke\n", " \n", " \n", " def assemble(self):\n", " \"\"\"Assemble element stiffness into global\"\"\"\n", " \n", " #Check if using natural coordinate quads, if so we will integrate using Gauss integration\n", " if len(connect[0]) == 4:\n", " \n", " ke = self.integrate_element_matrices()\n", " \n", " for i in range(self.num_elem):\n", " \n", " idx_grid = np.ix_(self.connect[i], self.connect[i])\n", " self.K[idx_grid] += ke[i]\n", " \n", " #Check if using triangles, if so we will use the analytic stiffness matrix\n", " elif len(connect[0]) == 3:\n", " \n", " c00 = self.Cmat[2,2]\n", " c11 = self.Cmat[0,0]\n", " c12 = self.Cmat[0,1]\n", " c21 = self.Cmat[1,0]\n", " c22 = self.Cmat[1,1]\n", " \n", " ke = np.array([[1 / 12. * (c00 + 6 * (c11 + c12 + c21 + c22)),\n", " 1 / 24. * (c00 - 12 * (c11 + c21)),\n", " 1 / 24. * (c00 - 12 * (c12 + c22))],\n", " [1 / 24. * (c00 - 12 * (c11 + c12)), \n", " 1 / 12. * (c00 + 6 * c11), \n", " 1 / 24. * (c00 + 12 * c12)],\n", " [1 / 24. * (c00 - 12 * (c21 + c22)), \n", " 1 / 24. * (c00 + 12 * c21),\n", " 1 / 12. * (c00 + 6 * c22)]])\n", " \n", " for i in range(self.num_elem):\n", " \n", " idx_grid = np.ix_(self.connect[i], self.connect[i])\n", " self.K[idx_grid] += ke\n", " \n", " \n", " def apply_essential_bc(self, nodes, values):\n", " \"\"\"\n", " Modifies stiffness matrix and r.h.s. vector for essential b.c.'s\n", " at location of nodes with corresponding values.\n", " \"\"\"\n", " \n", " #Subtract 1 to bring the node numbers into alignment with Python indexing\n", " node_idx = nodes - 1\n", " \n", " row_replace = np.zeros(self.num_dof)\n", " \n", " for value_idx, node in enumerate(node_idx): \n", " \n", " \n", " self.K[node] = row_replace\n", " self.K[node,node] = 1\n", " \n", " self.F[node] = values[value_idx]\n", " \n", " \n", " def solve(self):\n", " \"\"\"Solve the global system of equations\"\"\"\n", " \n", " self.K = scipy.sparse.csr_matrix(self.K)\n", " \n", " return scipy.sparse.linalg.spsolve(self.K,self.F)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can solve the problem. I will answer part (2) first using the square elements with linear isoparametric interpolation." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 61.28417886, 53.07365574, 30.64208943, 0. ,\n", " 70.29995898, 60.88155036, 35.14997949, 0. ,\n", " 100. , 86.60254038, 50. , 0. ])" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#The node locations\n", "nodes = np.array([[0,0],[1,0],[2,0],[3,0],\n", " [0,1],[1,1],[2,1],[3,1],\n", " [0,2],[1,2],[2,2],[3,2]], dtype=np.double)\n", "\n", "#The connectivity array for the quadratic elements using standard node numbering scheme\n", "connect = np.array([[1, 2, 6, 5],\n", " [2, 3, 7, 6],\n", " [3, 4, 8, 7],\n", " [5, 6, 10, 9],\n", " [6, 7, 11, 10],\n", " [7, 8, 12, 11]], dtype=np.int64)\n", "\n", "#Create a boundary condition node set and vaules for the top side\n", "ns1 = np.array([9, 10, 11, 12], dtype=np.int64)\n", "val1 = np.cos(np.pi * np.array([0., 1., 2., 3.]) / 6.) * 100\n", "\n", "#Create a boundary condition node set for the right side\n", "ns2 = np.array([4, 8, 12], dtype=np.int64)\n", "val2 = np.zeros(len(ns2))\n", "\n", "#Instatiate the problem\n", "problem = TwoDimFEM(nodes, connect)\n", "#Assemble\n", "problem.assemble()\n", "#Apply boundary conditions\n", "problem.apply_essential_bc(ns1,val1)\n", "problem.apply_essential_bc(ns2,val2)\n", "#Solve\n", "u = problem.solve()\n", "u" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/john/miniconda3/lib/python3.5/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n", " if self._edgecolors == str('face'):\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAADpCAYAAADoFbhNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE5NJREFUeJzt3X+s3XV9x/Hni1KCirMhsKLci8VRGEyEutmLxQKbmpVu\n1LmRQIPR4B8alyqZmWOaJrDEP+bMNiMCIw4ITqTLcJKSQMQ5ftQBlwJtKbRVOtG1dVwgQqdUI8h7\nf5zvLYfTe8733PP9fb6vR3LD+fHt+Xy+33vvs5/zvd9eFBGYmWVxWNUTMLPmc0jMLDOHxMwyc0jM\nLDOHxMwyc0jMLLOBIZE0KeluSU9IelzSp/ps92VJT0raJmlZMVM1s7o6POX5l4C/iIitko4CHpH0\nnYjYObuBpNXASRGxVNIUcC1wVnFTNrO6GbgiiYinI2JrcvvnwE7gLT2brQFuSraZBhZJWlzAXM2s\npoY+RyJpCbAMmO556nhgT9f9vcBE1omZWXMMFZLkbc2twGXJyuSQTXru+7p7sxZJO0eCpIXAN4Gv\nR8Rtc2yyD5jsuj+RPNb7Oo6LWQNERO/CINXAkEgScD2wIyK+1GezjcA6YIOks4AXImJmrg2/+9bT\n5zu/Utz0wgwfWZTttM7b156Z02xe64ubtvGZlWcU8toAv7HyPSP/2c/ffDvrL7kg0/jPnbQi05/v\n5x+uuoZPf/LP+z6/4xcnFzLurM07F/R97ru3/g3vvfCKkV/7kfufGvnPpvnWVUtH+nNpK5KzgQ8B\nj0nakjz2OeAEgIi4LiLukLRa0m7gReDSkWbSUEUFxJprUETG1cCQRMT3GOI8SkSsy21GDTEOAcmy\nGslDUauRNEWvRtoo9RxJG5xx5BuG3rbsgKw4oZifpOcRkXNOH/0bsuiIvHv5uwp9/SxOPO3cqqeQ\nO4cEOPPIowY+X+Xq4+y3HlfZ2GnOeccpVU+hr3dP1TckbzvtvKqnkDuHZIBxePsyl7a+pSlD0edH\nijzRmoVDModxDQg4Ij4/UgyHpMs4B8SsSA4J7QlI21cjRWvjj31ntTYkbYnHLEfEb2uK1LqQtC0g\nUH1EbPy1JiRtDEhd1GE1YsUa+5C0PSBVr0bqEpEq/21NG4xtSNoeEHBExk1dryGBMQuJ42FWjbEI\niQNyKK9GrEyNDokDMjdH5LV8fqR4jQyJA2JWL40JieMxHK9GrAq1D4kDMjxH5FC+mrUctQ2JAzI/\nVUekrXx+pKN2IXFAmqmOq5FxUudrSKBGIXFARlf1asQRsUpD4nhkV3VE6sznR8pTSUgckPHR5tWI\nz4+8qtSQOCD5qno10uaI2GsN/T8Rt3pxRAbz25pyOSQNVHVEzHo5JDZvXo34/Egvh6Rhql6NOCLl\nq/s1JOCQNErVETHrxyGxoXk1Yv04JA1R9Wqk7hEpk8+PHCo1JJJukDQjaXuf58+TtF/SluRjff7T\nbLeqI9IEXo1Ua5gL0m4ErgK+NmCbeyNiTT5Tsrqp+2rEEaleakgiYpOkJSmbKZfZ2EF1WYXUPSJW\nD3lcIh/ACknbgH3AX0bEjhxet5XqEpCmKHs14vMjc8sjJI8CkxFxQNL5wG2A15rzUNd4eDVSvSZc\nQwI5hCQiftZ1+05J10g6OiJ+2rvtFzdtO3h7xQmLOfutx2UdvtHqGhBoRkR8biS7Z/dO89y+6cyv\nkzkkkhYDz0RESFoOaK6IAHxm5RlZh2u8OsdjVhMiYvk4dmKKYyemDt7f9dBVI71Oakgk3QKcCxwj\naQ9wBbAQICKuAy4EPiHpZeAAcPFIMxlzTQhIk1SxGvH5kf6G+anN2pTnrwauzm1GY6SJ8WjCasRv\naeqnNr+zdZw0MSDQjIhYPTkkOWlqPJrGq5F6ckgyGKd4eDUymM+PDOaQjGCcAgLNiUjbViNNuYYE\nHJKhjVs8zPLkkKQY94B4NWJ5cEjmMO7xmOWIDMfnR9I5JF3aEhCzvLU+JG2Nh1cjlqdWhqSt8ZjV\nlIhYc7QqJG0PCDQrInVYjfj8yHDGPiSOhzVRk64hgTEOiQNyKK9GrChjFRLHoz9HxIo0FiFxQKwI\nPj8yvMaGxPEYnlcjVrTGhcQBmZ8mRcSaqxEhcTzawauR5qptSByP7Jq0GqlbRKo8P9K0H/1CDUPi\ngOSjSRGx5qtFSByPdqvbasTmr9KQOCDF8GrEylZ6SByPYjUtInVcjfj6kfk7rMzBHJFiNS0iNj5K\nDYlZtzquRmw0DsmYaNpqxBEZLw7JGGhaROqs6vMjTbyGBBwSq4BXI+PHIWk4r0asDhySBmtiRLwa\nGU+pIZF0g6QZSdsHbPNlSU9K2iZpWb5TNCtH1edHmmyYFcmNwKp+T0paDZwUEUuBjwHX5jQ3G8Cr\nEauT1JBExCbg+QGbrAFuSradBhZJWpzP9GwujojVTR7nSI4H9nTd3wtM5PC6ZtYQeZ1sVc/9yOl1\nrYdXI8Wow/mRpl5DAvn8o719wGTX/YnksUN8/ubbD94+5/STOecdp+QwfHs0MSJWb8/unea5fdOZ\nXyePkGwE1gEbJJ0FvBARM3NtuP6SC3IYzpqkCauRNjt2YopjJ6YO3t/10FUjvU5qSCTdApwLHCNp\nD3AFsBAgIq6LiDskrZa0G3gRuHSkmZhZY6WGJCLWDrHNunymY/34bY3Vma9stdarw4nWpnNIrDA+\nP9IeDkkD+G2N1Z1DYmaZOSRmlplDUnNNfVvj8yPt4pDUWFMjYu3jkNSUI2JN4pDUkCNSHl9Dkg+H\nxHLn8yPt45DUjFcj7dTkXyEADkmtOCLWVA5JTTgi1mQOSQ2MU0R8fqSdHJKKjVNErL0ckgo5IjYu\nHJKKOCLV8zUk+XFIKjCuEfH5kfZySEo2rhGxdnNISuSI2LhySEoy7hHx25p2c0hKMO4RMXNICuaI\nWBs4JAVyRKwtHJKCtCkiTTw/4mtI8uWQFKBNETEDhyR3jojNV9N/Fwk4JLlyRKytHJKctDUiTTw/\nYvlzSMwss9SQSFolaZekJyVdPsfz50naL2lL8rG+mKnWV1tXI2azDh/0pKQFwFeA9wH7gM2SNkbE\nzp5N742INQXNsdYcEbP0FclyYHdE/CgiXgI2AB+YYzvlPrMGaHtEmnp+xNeQ5C8tJMcDe7ru700e\n6xbACknbJN0h6bQ8J1hXbY+IWbeBb23oRCLNo8BkRByQdD5wG9DMv6qG5IiYvVZaSPYBk133J+ms\nSg6KiJ913b5T0jWSjo6In/a+2Odvvv3g7XNOP5lz3nHKSJOukiNi4+TZvdM8t2868+ukheRhYKmk\nJcBPgIuAtd0bSFoMPBMRIWk5oLkiArD+kgsyT7hKjsirmnp+xF7r2Ikpjp2YOnh/10NXjfQ6A0MS\nES9LWgd8G1gAXB8ROyV9PHn+OuBC4BOSXgYOABePNJOac0TM+ktbkRARdwJ39jx2Xdftq4Gr859a\nfTgiZoP5ytYUjohZOodkAEdkbk0+P+JrSIrhkPThiJgNzyGZgyNiZRmH30UCDskhHBGz+XNIujgi\n6Zp8fsSK45AkHBGz0TkkOCJmWbU+JI6IWXatDokjMj9NPz/ia0iK09qQOCJm+WltSMwsP60MiVcj\nZvlqXUgckdE0/fyIFatVIXFEzIrRmpA4ImbFaUVIHJFs/LbG0ox9SBwRA19DUrSxDokjYlaOsQ2J\nI2J1Ny6/iwTGNCSOSH58fsSGMXYhcUTMyjdWIXFEzKoxNiFxRMyqMxYhcUSK4fMjNqzGh8QRsTS+\nhqR4jQ6JI2JWD40NiSNiVh+NDIkjUjyfH7H5aFxIHBGz+kkNiaRVknZJelLS5X22+XLy/DZJy/Kf\nZocjYlZPA0MiaQHwFWAVcBqwVtKpPdusBk6KiKXAx4BrC5prYR6Y3lz1FPqq89y2bb636in0Vee5\n/XDHPVVPIXdpK5LlwO6I+FFEvARsAD7Qs80a4CaAiJgGFklanPdEi1yNPPBQfb9Zq5jbsOdHHnv4\nvoJnMro6z+2pHfWN3KjSQnI8sKfr/t7ksbRtJrJP7VV+S2Oj8jUk5UgLSQz5Ohrxz6VyRMzqTxH9\nv+clnQVcGRGrkvufBV6JiC90bfNPwD0RsSG5vws4NyJmel4rt7iYWXEiondhkOrwlOcfBpZKWgL8\nBLgIWNuzzUZgHbAhCc8LvREZdXJm1gwDQxIRL0taB3wbWABcHxE7JX08ef66iLhD0mpJu4EXgUsL\nn7WZ1crAtzZmZsPI/crWOl3ANt+5STpP0n5JW5KP9SXN6wZJM5K2D9imqmM2cG5VHbNk7ElJd0t6\nQtLjkj7VZ7vSj90wc6vw6+1ISdOStiZzu7LPdsMft4jI7YPO25/dwBJgIbAVOLVnm9XAHcntKeDB\nPOeQcW7nARvLmE/PuCuBZcD2Ps9XcsyGnFslxywZ+zjgzOT2UcD3a/T1Nszcqjx2r0/+ezjwIDCV\n5bjlvSKpzQVsI84NDv1RduEiYhPw/IBNqjpmw8wNKjhmABHxdERsTW7/HNgJvKVns0qO3ZBzg+qO\n3YHk5hF0/mJ9pWeTeR23vENSiwvY+hhmbgGsSJZyd0g6rYR5DaOqYzaMWhyz5CeLy4DpnqcqP3YD\n5lbZsZN0mKStwAxwV0T0XkI9r+OW9uPf+ar8ArYBhhnjUWAyIg5IOh+4DajLv6ev4pgNo/JjJuko\n4FbgsuRv/0M26blf2rFLmVtlxy4iXgHOlPQm4FuSficinujZbOjjlveKZB8w2XV/kk7JBm0zkTxW\ntNS5RcTPZpd8EXEnsFDS0SXMLU1VxyxV1cdM0kLgm8DXI+K2OTap7Nilza3qY5eMux+4m84/zO02\nr+OWd0gOXsAm6Qg6F7Bt7NlmI/BhOHjl7JwXsBUgdW6SFktScns5nR+P/7SEuaWp6pilqvKYJeNe\nD+yIiC/12aySYzfM3Ko6dpKOkbQouf064P10zuF0m9dxy/WtTdT4ArZh5gZcCHxC0svAAeDiMuYm\n6RbgXOAYSXuAK+icAKv0mA0zNyo6ZomzgQ8Bj0nakjz2OeCE2flVeOxS50Z1x+7NwE3q/JqQw4B/\nTY7TyN+nviDNzDJr3K9aNLP6cUjMLDOHxMwyc0jMLDOHxMwyc0jMLDOHxMwyc0jMLDOHxMwyc0jM\nLDOHxMwyc0jMLDOHxMwyc0jMLDOHxMwyc0jMLDOHxMwyc0jMLDOHxMwyc0jMLDOHxMwyy/v/tFco\nSf6V92YDREQl/y/hRoUE4IOffLK0sX53xYmljQXwrlN/Xep4AKe97gelj3nM7vtLH/P/Nn2v1PEe\nv2VrqeMBvPfH20sfc5bf2phZZg6JmWXmkJhZZg6JmWXmkAzwwx33lD7mts33lj7mA9ObSx3vvse+\nX+p4AP/146dLH3PrL39e+phVcUgGeGpH+d/Ujz18X+ljPvBQySHZXv5Piu7/n5nSx9z2yxdLH7Mq\nDomZZeaQmFlmimjOxaK+stVssKqubG1USMysnvzWxswyc0jMLLNGhETSKkm7JD0p6fKCxrhB0oyk\n7V2PXSlpr6QtyceqHMc7UtK0pK2SHpd0ZfL40ZK+I+kHku6StCivMbvGXpDsz+3J/cL2M3n9RZJu\nlbRT0g5JU0Xup6RTuvZli6T9ki4r+PN5maTtyefysuSxXPexz9fonGNIWiLpF137ek22PUwREbX+\nABYAu4ElwEJgK3BqAeOsBJYB27seuwL4dIH79vrkv4cDDwJTwN8Bf5U8fjnwtwWM+2ngZmBjSft5\nE/DRrn19Uxn7mbz2YcD/ApNF7SfwdmA7cGTy9fod4Lfy3sc+X6NzjpF8v2zPMt58PpqwIlkO7I6I\nH0XES8AG4AN5DxIRm4Dn53iqsLPgEXEguXkEnUgGsIbONx7Jf/8kzzElTQCrgX/m1X0TBe2npDcB\nKyPiBoCIeDki9lPwfnZ5H52vnz0Ut5+/DUxHxC8j4tfAvcCfkfM+9vkaLes4DtSEkBwP7Om6vzd5\nrCyflLRN0vV5v82QdJikrcAMcFdEPAQsjojZyzBngMV5jgn8I/AZ4JWux4Li9vNE4FlJN0p6VNJX\nJb2B4vdz1sXALcntovbzcWBl8jbj9XRCPUE5+zhojBOTtzX3SHpPAWMf1ISQVPnz6WvpfCOcSWd5\n/Pd5vnhEvBIRZ9L5opuS9Pae54Mc91/SHwPPRMQWXvs3c5H7eTjwTuCaiHgn8CLw190b5L2fsyQd\nAVwA/FvyUCH7GRG7gC8AdwF30nn7/euebQrZxwFj/ASYjIhldN7KfkPSG4sauwkh2Ufn/e2sSTqr\nksJFxDORoPNWYHlB4+wH7gb+EJiRdByApDcDz+Q41ApgjaSn6Pwt/QeSvlbwfu4F9kbE7D/ouZVO\nWJ4ucD9nnQ88EhHPQrGfz4i4ISJ+LyLOpfP24wcU+7mcNecYEfGriHg+uf0o8N/A0gLGB5oRkoeB\npclZ6COAi4CNZQycfGJmfZDOCbW8XvuYrjPsrwPeD+yks28fSTb7CHBbXmNGxOciYjIiTqSz5P/P\niPhwkfsZEU8DeySdnDz0PuAJ4HYK2s8ua3n1bU3Rn8/fTP57AvCnwDco8HPZZc4xkq+vBcntt9GJ\nyA8LGL+jrLO6WT7o/M3yfTo/vflsQWPcQmc5+Cs652Q+CnwNeAzYlnyCFuc43unAo8lrbwfWJ48f\nDfwHnb/R7gIWFbS/5/HqT23+paj9TF7/DGBz8vr/TuenNoXuJ/AG4DngjV2PFfn5vI9OILcCv1/E\n53KOr9FL+41BJ2aPA1uAR4A/KuLraPbDl8ibWWZNeGtjZjXnkJhZZg6JmWXmkJhZZg6JmWXmkJhZ\nZg6JmWXmkJhZZv8P2xRnWkrmmk8AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "#Plot the results\n", "plot = plt.contourf(nodes[:,0].reshape(3,4), nodes[:,1].reshape(3,4), u.reshape(3,4),cmap=\"coolwarm\")\n", "plt.colorbar(plot, orientation='horizontal', shrink=0.6);\n", "plt.clim(0,100)\n", "plt.axes().set_aspect('equal')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will solve part (1) using the triangles and an analytic stiffness matrix" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#The node locations are the same, so we will updated the connectivities\n", "connect = np.array([[5, 1, 6],\n", " [2, 6, 1],\n", " [6, 2, 7],\n", " [3, 7, 2],\n", " [7, 3, 8],\n", " [4, 8, 3],\n", " [9, 5, 10],\n", " [6, 10, 5],\n", " [10, 6, 11],\n", " [7, 11, 6],\n", " [11, 7, 12],\n", " [8, 12, 7]], dtype=np.int64)\n", "\n", "#Instantiate the problem\n", "problem = TwoDimFEM(nodes, connect)\n", "#Assemble\n", "problem.assemble()\n", "#Apply boundary conditions\n", "problem.apply_essential_bc(ns1,val1)\n", "problem.apply_essential_bc(ns2,val2)\n", "#Solve\n", "u = problem.solve()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/john/miniconda3/lib/python3.5/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n", " if self._edgecolors == str('face'):\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAADpCAYAAADoFbhNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE11JREFUeJzt3X+s3XV9x/Hni1ICirMhsKLci+AojE6EutmLxQKbmpVu\n1LmRSKPR4B8alyqZiWOaJrDEP+bMNiMCIxsYnEiX4SQlKRHm+FEnXgq0pdBW6UTX1nGBCJ1SjSDv\n/XG+tz2c3nO+557v7/N9PZKTe873fO738/l+z7mv+/5+zqe3igjMzLI4quoBmFnzOUjMLDMHiZll\n5iAxs8wcJGaWmYPEzDIbGCSSJiXdK+kJSY9L+mSfdl+S9KSk7ZKWFTNUM6uro1Oefwn4i4jYJul4\n4BFJ90TErtkGklYDZ0TEEklTwA3A+cUN2czqZmBFEhFPR8S25P7PgV3AG3uarQFuSdpMA4skLS5g\nrGZWU0PPkUg6DVgGTPc8dQqwt+vxPmAi68DMrDmGCpLksuZ24MqkMjmiSc9jr7s3a5G0ORIkLQS+\nAXwtIu6Yo8l+YLLr8USyrXc/DhezBoiI3sIg1cAgkSTgJmBnRHyxT7ONwDpgg6TzgRciYmauht9+\n0znzHV8pbnlhhg8vyjat85a15+U0mlf7wubtfHrluYXsG+A3Vr5z5O/93K13sv4Dl2Yew3NnrMi8\nj15/f+31fOoTf973+Z2/ODP3Prtt2bWg73Pfvv2veddlV4+870e++9TI35vmm9cuGen70iqSC4AP\nAo9J2pps+yxwKkBE3BgRmyStlrQHeBG4YqSRNFRRAdIWRYRI1QaFyLgaGCQR8R2GmEeJiHW5jagh\nxiFAslQjeRjHEGmr1DmSNjj32NcO3bbsAFlxan0/Sb/wnGIvD7J4x/K3932u6MuaNKcvvajS/ovg\nIAHOO/b41DZVVSAXvOnkQvabRzVy4VvPGvl7i65G3jHVP0iq9ualF1c9hNw5SFKMwyVML1/SFKeN\n8yPgIOlrHAPEqr+syarIT2yycJD0GPcAcTViRXCQJMY9QMAhYsVpfZC0IUCsHG2dH4EWB0nbAsTV\nSPPnR+qsdUHStgCB6kPExl9rgqSNAVIXdahGrFhjHSQOj+qrkbqEiC9rijWWQeIA6ag6RNqkjInW\nuq4hgTELEgdIvdSlGrHijUWQOECOVHU14hBpl0YHiQNkblWHSN14fqR4jQwSB0i9ta0aafNCtFmN\nChIHSDpXI1aFRgSJA2Q4dQiRulUjvqwpR62DxAHSLHULEStPLYPEATJ/dahG2sjzIx21ChIHyGjq\nECKuRopV58VoUJMgcYA0W11DxPMj5ak0SBwg2dWhGjErPUgcHvmpQ4jUtRqxcg31n4jnxSEyXtoe\nIp5oPazUILH81KEaqTPPj5TLQdJAdQiROlcjDpHyOUjMLDMHScO4GhmsrGrE8yOv5iBpEIdIO9V9\nMRoMESSSbpY0I2lHn+cvlnRA0tbktj7/YZql89xIdYZZR/IV4FrgqwPa3B8Ra/IZks3F1YjVWWpF\nEhGbgedTmimf4dhcHCLpyqxGPD9ypDzmSAJYIWm7pE2SluawTzNrkDyWyD8KTEbEQUmXAHcAvljN\nQR0qEXA1YukyB0lE/Kzr/l2Srpd0QkT8tLftFzZvP3R/xamLueBNJ2ftfizVJUCg/iFi2Ty7b5rn\n9k9n3k/mIJG0GHgmIkLSckBzhQjAp1eem7W7sVanAGkKVyPZnDQxxUkTU4ce737o2pH2kxokkm4D\nLgJOlLQXuBpYCBARNwKXAR+X9DJwELh8pJG0VJ3Dw9XIkcqeaG3CGhIYIkgiYm3K89cB1+U2opao\nc4A0hauR+qjFX0hriyaFh6sRmw8HSQmaFCDQjBBxNVIvDpICNS1AbDAvROvPQZKzpoeHqxEbhYMk\nJ00PEGhGiFg9OUgyGIfwaBpXI/XkIBnBOAZIE6qRKkPE8yODOUjmYRwDxOqrKYvRwEGSqg3h4WrE\nsnKQ9NGGAIFmhIjVn4OkS1vCo2mqrkY8P5LOQUJ7A8TViOWl1UHS1gCB5oRI1dWIDad1QdLm8DAr\nSmuCxAFymKsRy9tYB4nD40hNCZG68ETrcMYySBwgzdf2aqRJi9FgzILEATKYqxErSuODxOExftpe\njTRRY4PEATI/rkbmz/Mjw2tUkDg8RtOkEHE10kyNCBIHiFm91TpIHCDZuRqxMtQuSBwe+WlSiNSN\n50fmpzZB4gBpN1cjhzVtDQlUHCQOj+I0qRpxiDRfJUHiADEbL6UGiQOkHK5GrGxHVT0Ay1eTQqSu\nPNE6fw4Sq4yrkfGRGiSSbpY0I2nHgDZfkvSkpO2SluU7RBuWqxGryjAVyVeAVf2elLQaOCMilgAf\nBW7IaWw2D00LEVcj4yU1SCJiM/D8gCZrgFuSttPAIkmL8xmeWbk8PzKaPOZITgH2dj3eB0zksF8b\nkquR8dHExWiQ32Sreh5HTvu1FE0LERtPeawj2Q9Mdj2eSLYd4XO33nno/oXnnMmFbz0rh+6tSVyN\n1Muz+6Z5bv905v3kESQbgXXABknnAy9ExMxcDdd/4NIcurNZrkYsq5MmpjhpYurQ490PXTvSflKD\nRNJtwEXAiZL2AlcDCwEi4saI2CRptaQ9wIvAFSONxKxinmgdXWqQRMTaIdqsy2c4Ns58WTO+vLLV\nzDJzkDSU50esThwkZpaZg8TMMnOQWCk80TreHCQN5PkRqxsHiZll5iAxw4vRsnKQNEwTL2s8PzL+\nHCRmNdHUPyEADhIzy4GDxMwyc5A0iOdHrK4cJGaWmYPEzDJzkDREEy9rrD0cJNZ6XoyWnYPECuOJ\n1vZwkJhZZg6SBvD8iNWdg8TMMnOQWCE8P9IuDhIzy8xBUnOeH7EmcJDUmEPEmsJBUlMOkXJ4MVo+\nHCQ11PQQ8URr+zhIaqbpIWKjafJfRwMHiZnlwEFSI65GrKlSg0TSKkm7JT0p6ao5nr9Y0gFJW5Pb\n+mKGOt7GJUQ8P9JORw96UtIC4MvAu4H9wBZJGyNiV0/T+yNiTUFjHHvjEiLWXmkVyXJgT0T8KCJe\nAjYA752jnXIfWUs4RGwcpAXJKcDersf7km3dAlghabukTZKW5jnAceYQsXEx8NKGTkikeRSYjIiD\nki4B7gB8oZxiHEOkafMjXoyWn7Qg2Q9Mdj2epFOVHBIRP+u6f5ek6yWdEBE/7d3Z526989D9C885\nkwvfetZIg266cQwRa6Zn903z3P7pzPtJC5KHgSWSTgN+ArwfWNvdQNJi4JmICEnLAc0VIgDrP3Bp\n5gGbWX5OmpjipImpQ493P3TtSPsZGCQR8bKkdcC3gAXATRGxS9LHkudvBC4DPi7pZeAgcPlII2kJ\nVyM2jtIqEiLiLuCunm03dt2/Drgu/6GNH4eIjSuvbC3JuIdI0yZaLV8OkhKMe4iYOUgK5hCxNnCQ\nFMghYm3hILHMmjg/4sVo+XKQFMTViLWJg6QADhGbj6b/dTRwkOTOIWJt5CDJkUPE2spBkpO2hkgT\nJ1otfw6SHLQ1RMxmOUgycoiYOUjMLAcOkgzaXo00dX7Ei9Hy5yAZUdtDxKybg2QEDhGzV3OQzJND\nxOxIDpJ5cIgc1tT5ESuGg2RIDhGz/hwkZpaZg2QIrkbMBnOQpHCImKVzkAzgEJlbkydavRitGA6S\nPhwiZsNzkMzBIWJlGYe/jgYOkiM4RMzmz0HSxSGSrsnzI1YcB4mZZeYgSbgaMRudgwSHiFlWqUEi\naZWk3ZKelHRVnzZfSp7fLmlZ/sMsjkNkeJ4fsX4GBomkBcCXgVXAUmCtpLN72qwGzoiIJcBHgRsK\nGmvuZkPkwektFY+kvzqPbfuW+6seQl9zja0ui9F+uPO+qoeQu7SKZDmwJyJ+FBEvARuA9/a0WQPc\nAhAR08AiSYtzH2nOuiuRBx+q7w9rncf22MMPVD2Evuo8tqd21jeAR5UWJKcAe7se70u2pbWZyD60\n4vhyxixfaUESQ+5HI36fmY0BRfT/mZd0PnBNRKxKHn8GeCUiPt/V5h+B+yJiQ/J4N3BRRMz07Mvh\nYtYAEdFbGKQ6OuX5h4Elkk4DfgK8H1jb02YjsA7YkATPC70hMurgzKwZBgZJRLwsaR3wLWABcFNE\n7JL0seT5GyNik6TVkvYALwJXFD5qM6uVgZc2ZmbDyH1la50XsKWNTdLFkg5I2prc1pc0rpslzUja\nMaBNVeds4NiqOmdJ35OS7pX0hKTHJX2yT7vSz90wY6vw/XaspGlJ25KxXdOn3fDnLSJyu9G5/NkD\nnAYsBLYBZ/e0WQ1sSu5PAd/LcwwZx3YxsLGM8fT0uxJYBuzo83wl52zIsVVyzpK+TwbOS+4fD3y/\nRu+3YcZW5bl7TfL1aOB7wFSW85Z3RVLnBWzDjA2O/Ci7cBGxGXh+QJPKFv0NMTao4JwBRMTTEbEt\nuf9zYBfwxp5mlZy7IccG1Z27g8ndY+j8Yn2lp8m8zlveQVLnBWzDjC2AFUkpt0nS0hLGNYw6L/qr\nxTlLPllcBkz3PFX5uRswtsrOnaSjJG0DZoC7I6J3CfW8zlvax7/zVecFbMP08SgwGREHJV0C3AHU\n5V+q1XXRX+XnTNLxwO3Alclv/yOa9Dwu7dyljK2ycxcRrwDnSXo98E1JvxMRT/Q0G/q85V2R7Acm\nux5P0kmyQW0mkm1FSx1bRPxstuSLiLuAhZJOKGFsaao6Z6mqPmeSFgLfAL4WEXfM0aSyc5c2tqrP\nXdLvAeBeOv8wt9u8zlveQXJoAZukY+gsYNvY02Yj8CE4tHJ2zgVsBUgdm6TFkpTcX07n4/GfljC2\nNFWds1RVnrOk35uAnRHxxT7NKjl3w4ytqnMn6URJi5L7xwHvoTOH021e5y3XS5uo8QK2YcYGXAZ8\nXNLLwEHg8jLGJuk24CLgREl7gavpTIBVes6GGRsVnbPEBcAHgcckbU22fRY4dXZ8FZ671LFR3bl7\nA3CLOn8m5CjgX5PzNPLPqRekmVlm/lOLZpaZg8TMMnOQmFlmDhIzy8xBYmaZOUjMLDMHiZll5iAx\ns8wcJGaWmYPEzDJzkJhZZg4SM8vMQWJmmTlIzCwzB4mZZeYgMbPMHCRmlpmDxMwyc5CYWWYOEjPL\nzEFiZpnl/T/tFUqS/+S92QARUcn/JdyoIAF43yeeLK2v311xeml9Abz97F+X2h/A0uN+UHqfJ+75\nbul9/t/m75Ta3+O3bSu1P4B3/XhH6X3O8qWNmWXmIDGzzBwkZpaZg8TMMnOQDPDDnfeV3uf2LfeX\n3ueD01tK7e+Bx75fan8A//Xjp0vvc9svf156n1VxkAzw1M7yf6gfe/iB0vt88KGSg2RH+Z8Uffd/\nZkrvc/svXyy9z6o4SMwsMweJmWWmiOYsFvXKVrPBqlrZ2qggMbN68qWNmWXmIDGzzBoRJJJWSdot\n6UlJVxXUx82SZiTt6Np2jaR9krYmt1U59nespGlJ2yQ9LumaZPsJku6R9ANJd0talFefXX0vSI7n\nzuRxYceZ7H+RpNsl7ZK0U9JUkccp6ayuY9kq6YCkKwt+Pa+UtCN5La9MtuV6jH3eo3P2Iek0Sb/o\nOtbrsx1hioio9Q1YAOwBTgMWAtuAswvoZyWwDNjRte1q4FMFHttrkq9HA98DpoC/Bf4y2X4V8DcF\n9Psp4FZgY0nHeQvwka5jfX0Zx5ns+yjgf4HJoo4TeAuwAzg2eb/eA/xW3sfY5z06Zx/Jz8uOLP3N\n59aEimQ5sCcifhQRLwEbgPfm3UlEbAaen+OpwmbBI+JgcvcYOiEZwBo6P3gkX/8kzz4lTQCrgX/m\n8LGJgo5T0uuBlRFxM0BEvBwRByj4OLu8m877Zy/FHedvA9MR8cuI+DVwP/Bn5HyMfd6jZZ3HgZoQ\nJKcAe7se70u2leUTkrZLuinvywxJR0naBswAd0fEQ8DiiJhdhjkDLM6zT+AfgE8Dr3RtC4o7ztOB\nZyV9RdKjkv5J0msp/jhnXQ7cltwv6jgfB1YmlxmvoRPUE5RzjIP6OD25rLlP0jsL6PuQJgRJlZ9P\n30DnB+E8OuXx3+W584h4JSLOo/Omm5L0lp7ngxyPX9IfA89ExFZe/Zu5yOM8GngbcH1EvA14Efir\n7gZ5H+csSccAlwL/lmwq5DgjYjfweeBu4C46l9+/7mlTyDEO6OMnwGRELKNzKft1Sa8rqu8mBMl+\nOte3sybpVCWFi4hnIkHnUmB5Qf0cAO4F/hCYkXQygKQ3AM/k2NUKYI2kp+j8lv4DSV8t+Dj3Afsi\nYvYf9NxOJ1ieLvA4Z10CPBIRz0Kxr2dE3BwRvxcRF9G5/PgBxb6Ws+bsIyJ+FRHPJ/cfBf4bWFJA\n/0AzguRhYEkyC30M8H5gYxkdJy/MrPfRmVDLa98nds2wHwe8B9hF59g+nDT7MHBHXn1GxGcjYjIi\nTqdT8v9nRHyoyOOMiKeBvZLOTDa9G3gCuJOCjrPLWg5f1hT9ev5m8vVU4E+Br1Pga9llzj6S99eC\n5P6b6YTIDwvov6OsWd0sNzq/Wb5P59ObzxTUx210ysFf0ZmT+QjwVeAxYHvyAi3Osb9zgEeTfe8A\n1ifbTwD+g85vtLuBRQUd78Uc/tTmX4o6zmT/5wJbkv3/O51PbQo9TuC1wHPA67q2Ffl6PkAnILcB\nv1/EaznHe/SKfn3QCbPHga3AI8AfFfE+mr15ibyZZdaESxszqzkHiZll5iAxs8wcJGaWmYPEzDJz\nkJhZZg4SM8vMQWJmmf0/3aU5hYxpWUkAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Plot the results\n", "plot = plt.contourf(nodes[:,0].reshape(3,4), nodes[:,1].reshape(3,4), u.reshape(3,4),cmap=\"coolwarm\")\n", "plt.colorbar(plot, orientation='horizontal', shrink=0.6);\n", "plt.clim(0,100)\n", "plt.axes().set_aspect('equal')" ] } ], "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.0" } }, "nbformat": 4, "nbformat_minor": 0 }