{
"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
}