{"id":365,"date":"2019-01-29T01:19:52","date_gmt":"2019-01-29T09:19:52","guid":{"rendered":"https:\/\/gantovnik.com\/bio-tips\/?p=365"},"modified":"2021-11-12T15:53:15","modified_gmt":"2021-11-12T23:53:15","slug":"solution-of-laplace-equation-using-fem","status":"publish","type":"post","link":"https:\/\/gantovnik.com\/bio-tips\/2019\/01\/solution-of-laplace-equation-using-fem\/","title":{"rendered":"#51 Solution of Laplace equation using FEM"},"content":{"rendered":"<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nimport os\r\nimport matplotlib.pyplot as plt\r\nimport numpy as np\r\nfrom numpy import linspace\r\nfrom scipy.spatial import Delaunay\r\nfrom numpy import cross\r\nfrom scipy.sparse import dok_matrix\r\nimport matplotlib.cm as cm\r\nos.chdir(r'D:\\projects\\wordpress\\ex51')\r\nos.getcwd()\r\nxmin = 0 ; xmax = 1 ; nXpoints = 10\r\nymin = 0 ; ymax = 1 ; nYpoints = 10\r\nhorizontal = linspace(xmin,xmax,nXpoints)\r\nvertical = linspace(ymin,ymax,nYpoints)\r\ny, x = np.meshgrid(horizontal, vertical)\r\nvertices = np.array(&#x5B;x.flatten(),y.flatten()])\r\ntriangulation = Delaunay(vertices.T)\r\nindex2point = lambda index: triangulation.points&#x5B;index]\r\nall_centers = index2point(triangulation.vertices).mean(axis=1)\r\ntrngl_set=triangulation.vertices\r\nplt.triplot(vertices&#x5B;0],vertices&#x5B;1],triangles=trngl_set)\r\nplt.savefig(&quot;example51_1.png&quot;, dpi=100)\r\nplt.show()\r\npoints=triangulation.points.shape&#x5B;0]\r\nstiff_matrix=dok_matrix((points,points))\r\nfor triangle in triangulation.vertices:\r\n    helper_matrix=dok_matrix((points,points))\r\n    pt1,pt2,pt3=index2point(triangle)\r\n    area=abs(0.5*cross(pt2-pt1,pt3-pt1))\r\n    coeffs=0.5*np.vstack((pt2-pt3,pt3-pt1,pt1-pt2))\/area\r\n    #helper_matrix&#x5B;triangle,triangle] = \\\r\n    np.array(np.mat(coeffs)*np.mat(coeffs).T)\r\n    u=None\r\n    u=np.array(np.mat(coeffs)*np.mat(coeffs).T)\r\n    for i in range(len(triangle)):\r\n        for j in range(len(triangle)):\r\n            helper_matrix&#x5B;triangle&#x5B;i],triangle&#x5B;j]] = u&#x5B;i,j]\r\n    stiff_matrix=stiff_matrix+helper_matrix\r\n\r\nallNodes = np.unique(trngl_set)\r\nboundaryNodes = np.unique(triangulation.convex_hull)\r\nNonBoundaryNodes = np.array(&#x5B;])\r\nfor x in allNodes:\r\n    if x not in boundaryNodes:\r\n        NonBoundaryNodes = np.append(NonBoundaryNodes,x)\r\n\r\nNonBoundaryNodes = NonBoundaryNodes.astype(int)\r\nnbnodes = len(boundaryNodes) # number of boundary nodes\r\nFbVals=np.zeros(&#x5B;nbnodes,1]) # Values on the boundary\r\nFbVals&#x5B;(nbnodes-nXpoints+1):-1]=np.ones(&#x5B;nXpoints-2, 1])\r\ntotalNodes = len(allNodes)\r\nstiff_matrixDense = stiff_matrix.todense()\r\nstiffNonb = stiff_matrixDense&#x5B;np.ix_(NonBoundaryNodes,NonBoundaryNodes)]\r\nstiffAtb = stiff_matrixDense&#x5B;np.ix_(NonBoundaryNodes,boundaryNodes)]\r\nU=np.zeros(&#x5B;totalNodes, 1])\r\nU&#x5B;NonBoundaryNodes] = np.linalg.solve( - stiffNonb,stiffAtb * FbVals )\r\nU&#x5B;boundaryNodes] = FbVals\r\nX = vertices&#x5B;0]\r\nY = vertices&#x5B;1]\r\nZ = U.T.flatten()\r\nfrom mpl_toolkits.mplot3d import axes3d\r\nfig = plt.figure()\r\nax = fig.add_subplot(111, projection='3d')\r\nsurf = ax.plot_trisurf(X, Y, Z, cmap=cm.jet, linewidth=0)\r\nfig.colorbar(surf)\r\nax.set_xlabel('X',fontsize=14)\r\nax.set_ylabel('Y',fontsize=14)\r\nax.set_zlabel(r&quot;$\\phi$&quot;,fontsize=16)\r\nplt.savefig(&quot;example51_2.png&quot;, dpi=100)\r\nplt.show()\r\nfrom numpy import pi, sinh, sin, cos\r\ndef f(x,y):\r\n    return sum( 2*(1.0\/(n*pi) - cos(n*pi)\/(n*pi))*(sinh(n*pi*x)\/ \\\r\n                   sinh(n*pi))*sin(n*pi*y) for n in range(1,200))\r\nZe = f(X,Y)\r\nZdiffZe = Ze - Z\r\nplt.plot(ZdiffZe)\r\nplt.savefig(&quot;example51_3.png&quot;, dpi=100)\r\nplt.show()\r\nplt.close()\r\n<\/pre>\n<p><img data-recalc-dims=\"1\" decoding=\"async\" src=\"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_1.png?resize=600%2C400&#038;ssl=1\" alt=\"example61_1\" width=\"600\" height=\"400\" class=\"alignnone size-full wp-image-366\" srcset=\"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_1.png?w=600&amp;ssl=1 600w, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_1.png?resize=300%2C200&amp;ssl=1 300w\" sizes=\"(max-width: 600px) 100vw, 600px\" \/><\/p>\n<p><img data-recalc-dims=\"1\" decoding=\"async\" src=\"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_2.png?resize=600%2C400&#038;ssl=1\" alt=\"example61_2\" width=\"600\" height=\"400\" class=\"alignnone size-full wp-image-367\" srcset=\"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_2.png?w=600&amp;ssl=1 600w, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_2.png?resize=300%2C200&amp;ssl=1 300w\" sizes=\"(max-width: 600px) 100vw, 600px\" \/><\/p>\n<p><img data-recalc-dims=\"1\" decoding=\"async\" src=\"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_3.png?resize=600%2C400&#038;ssl=1\" alt=\"example61_3\" width=\"600\" height=\"400\" class=\"alignnone size-full wp-image-368\" srcset=\"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_3.png?w=600&amp;ssl=1 600w, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example61_3.png?resize=300%2C200&amp;ssl=1 300w\" sizes=\"(max-width: 600px) 100vw, 600px\" \/><\/p>\n","protected":false},"excerpt":{"rendered":"<p>import os import matplotlib.pyplot as plt import numpy as np from numpy import linspace from scipy.spatial import Delaunay from numpy import cross from scipy.sparse import dok_matrix import matplotlib.cm as cm os.chdir(r&#8217;D:\\projects\\wordpress\\ex51&#8242;) os.getcwd() xmin = 0 ; xmax = 1 ; nXpoints = 10 ymin = 0 ; ymax = 1 ; nYpoints = 10 horizontal [&hellip;]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"nf_dc_page":"","_et_pb_use_builder":"","_et_pb_old_content":"","_et_gb_content_width":"","_lmt_disableupdate":"yes","_lmt_disable":"","jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[37,2],"tags":[],"class_list":["post-365","post","type-post","status-publish","format-standard","hentry","category-fem","category-python"],"modified_by":"gantovnik","jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"jetpack_shortlink":"https:\/\/wp.me\/p8bH0k-5T","jetpack_likes_enabled":true,"jetpack-related-posts":[{"id":1742,"url":"https:\/\/gantovnik.com\/bio-tips\/2023\/01\/204-mandelbrot-fractal-using-python-2-2-2-2-2-2-2-2-2\/","url_meta":{"origin":365,"position":0},"title":"#330 Gradient method using python","author":"gantovnik","date":"2023-01-03","format":false,"excerpt":"golden.py [code language=\"python\"] import math as mt def golden(f,a,b,tol=1.0e-10): # Golden section method for determining x # that minimizes the scalar function f(x). # The minimum must be bracketed in (a,b). c1 = (mt.sqrt(5.0)-1.0)\/2.0 c2 = 1.0 - c1 nIt = int(mt.ceil(mt.log(tol\/abs(a-b))\/mt.log(c1))) # first step x1 = c1*a + c2*b\u2026","rel":"","context":"In &quot;optimization&quot;","block_context":{"text":"optimization","link":"https:\/\/gantovnik.com\/bio-tips\/category\/optimization\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex330.png?resize=350%2C200&ssl=1","width":350,"height":200},"classes":[]},{"id":1738,"url":"https:\/\/gantovnik.com\/bio-tips\/2023\/01\/204-mandelbrot-fractal-using-python-2-2-2-2-2-2-2-2\/","url_meta":{"origin":365,"position":1},"title":"#329 Golden section method using python","author":"gantovnik","date":"2023-01-03","format":false,"excerpt":"golden.py [code language=\"python\"] import math as mt def golden(f,a,b,tol=1.0e-10): # Golden section method for determining x # that minimizes the scalar function f(x). # The minimum must be bracketed in (a,b). c1 = (mt.sqrt(5.0)-1.0)\/2.0 c2 = 1.0 - c1 nIt = int(mt.ceil(mt.log(tol\/abs(a-b))\/mt.log(c1))) # first step x1 = c1*a + c2*b\u2026","rel":"","context":"In &quot;optimization&quot;","block_context":{"text":"optimization","link":"https:\/\/gantovnik.com\/bio-tips\/category\/optimization\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex329.png?resize=350%2C200&ssl=1","width":350,"height":200},"classes":[]},{"id":260,"url":"https:\/\/gantovnik.com\/bio-tips\/2019\/01\/interpolation\/","url_meta":{"origin":365,"position":2},"title":"Interpolation","author":"gantovnik","date":"2019-01-15","format":false,"excerpt":"[code language=\"python\"] import os import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import interp1d os.chdir(r'D:\\data\\scripts\\wordpress\\ex50') os.getcwd() A,nu,k=10,4,2 def f(x,A,nu,k): return A*np.exp(-k*x)*np.cos(2*np.pi*nu*x) xmax,nx=0.5,8 x=np.linspace(0,xmax,nx) y=f(x,A,nu,k) f_nearest=interp1d(x,y,kind='nearest') f_linear =interp1d(x,y) f_cubic =interp1d(x,y,kind='cubic') x2=np.linspace(0,xmax,100) plt.plot(x,y,'o',label='datapoints') plt.plot(x2,f(x2,A,nu,k),label='exact') plt.plot(x2,f_nearest(x2),label='nearest') plt.plot(x2, f_linear(x2), label='linear') plt.plot(x2, f_cubic(x2), label='cubic') plt.legend() plt.show() plt.tight_layout() plt.savefig(\"example50.png\", dpi=100) plt.show() plt.close() [\/code]","rel":"","context":"In &quot;python&quot;","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"example50","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example50-1.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example50-1.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example50-1.png?resize=525%2C300 1.5x"},"classes":[]},{"id":1182,"url":"https:\/\/gantovnik.com\/bio-tips\/2021\/11\/204-mandelbrot-fractal-using-python-2-2\/","url_meta":{"origin":365,"position":3},"title":"#206 Triangulation using matplotlib","author":"gantovnik","date":"2021-11-26","format":false,"excerpt":"[code language=\"python\"] import matplotlib.pyplot as plt from matplotlib.tri import Triangulation from matplotlib.patches import Polygon import numpy as np def update_polygon(tri): if tri == -1: points = [0, 0, 0] else: points = triang.triangles[tri] xs = triang.x[points] ys = triang.y[points] polygon.set_xy(np.column_stack([xs, ys])) def on_mouse_move(event): if event.inaxes is None: tri = -1\u2026","rel":"","context":"In &quot;python&quot;","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2021\/11\/ex206.png?resize=350%2C200&ssl=1","width":350,"height":200},"classes":[]},{"id":2007,"url":"https:\/\/gantovnik.com\/bio-tips\/2024\/01\/401-add-labels-to-line-plot-in-python\/","url_meta":{"origin":365,"position":4},"title":"#401 Add labels to line plot in python","author":"gantovnik","date":"2024-01-05","format":false,"excerpt":"[code language=\"python\"] import matplotlib.pyplot as plt import numpy as np plt.clf() # using some dummy data for this example n = 25 xs = np.linspace(1, 100, num=n) ys = np.random.normal(loc=3, scale=0.4, size=n) # 'bo-' means blue color, round points, solid lines plt.plot(xs,ys,'bo-') # zip joins x and y coordinates in\u2026","rel":"","context":"In &quot;matplotlib&quot;","block_context":{"text":"matplotlib","link":"https:\/\/gantovnik.com\/bio-tips\/category\/matplotlib\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2024\/01\/ex401.png?fit=640%2C480&ssl=1&resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2024\/01\/ex401.png?fit=640%2C480&ssl=1&resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2024\/01\/ex401.png?fit=640%2C480&ssl=1&resize=525%2C300 1.5x"},"classes":[]},{"id":654,"url":"https:\/\/gantovnik.com\/bio-tips\/2020\/09\/107-tsa-data-2020-vs-2019\/","url_meta":{"origin":365,"position":5},"title":"#107: TSA Data 2020 vs 2019","author":"gantovnik","date":"2020-09-28","format":false,"excerpt":"#107: TSA Data 2020 vs 2019 The TSA has started to publish the daily volume of passengers going through checkpoints on its website. The data set also includes the numbers from 2019 in order to measure the impact of travel as a result of COVID-19. https:\/\/www.tsa.gov\/coronavirus\/passenger-throughput [code language=\"python\"] from bs4\u2026","rel":"","context":"In &quot;python&quot;","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2020\/09\/ex107.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2020\/09\/ex107.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2020\/09\/ex107.png?resize=525%2C300 1.5x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2020\/09\/ex107.png?resize=700%2C400 2x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2020\/09\/ex107.png?resize=1050%2C600 3x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2020\/09\/ex107.png?resize=1400%2C800 4x"},"classes":[]}],"_links":{"self":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts\/365","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/comments?post=365"}],"version-history":[{"count":0,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts\/365\/revisions"}],"wp:attachment":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/media?parent=365"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/categories?post=365"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/tags?post=365"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}