Skip to content

Commit 3f7d2ff

Browse files
com-pycom-py
authored
Jupyter notebook version added
1 parent 2cc9026 commit 3f7d2ff

File tree

3 files changed

+268
-0
lines changed

3 files changed

+268
-0
lines changed

‎ch07/Program_7.3_laplacefem.ipynb

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"#\n",
10+
"# Program 7.3: Laplace equation by FEM (laplacefem.ipynb)\n",
11+
"# J Wang, Computational modeling and visualization with Python\n",
12+
"#\n",
13+
"\n",
14+
"from scipy.linalg import solve\n",
15+
"from mpl_toolkits.mplot3d import Axes3D\n",
16+
"import matplotlib.pyplot as plt, numpy as np\n",
17+
"%matplotlib notebook\n",
18+
"\n",
19+
"def mesh(L, N): # generate mesh\n",
20+
" elm, bv = [], [] # elements: $[[n_1,n_2,n_3],..]$, bndry value\n",
21+
" x, y = np.linspace(0,L,N+1), np.linspace(0,L,N+1) # same x,y grids\n",
22+
" ndn = np.arange((N+1)*(N+1)).reshape(N+1,N+1) # node number \n",
23+
" node = [[xi, yj] for xi in x for yj in y] # nodes \n",
24+
" for j in range(N):\n",
25+
" for i in range(N):\n",
26+
" elm.append([ndn[i,j], ndn[i+1,j+1], ndn[i,j+1]]) # upper \n",
27+
" elm.append([ndn[i,j], ndn[i+1,j], ndn[i+1,j+1]]) # lower \n",
28+
"\n",
29+
" ip = ndn[1:-1,1:-1].flatten() # internal nodes \n",
30+
" bp = np.delete(ndn, ip) # boundary nodes=all-internal \n",
31+
" for p in bp:\n",
32+
" bv.append((node[p][0]*node[p][1])**2) # boundary values\n",
33+
" return node, elm, bp, ip, bv, x, y\n",
34+
"\n",
35+
"def fem_mat(node, elm): # fills matrix $A_{ij} = \\langle \\nabla \\phi_i \\cdot \\nabla \\phi_j \\rangle $\n",
36+
" A = np.zeros((len(node),len(node)))\n",
37+
" for e in elm:\n",
38+
" (x1,y1), (x2,y2), (x3,y3) = node[e[0]], node[e[1]], node[e[2]]\n",
39+
" beta, gama = [y2-y3, y3-y1, y1-y2], [x3-x2, x1-x3, x2-x1]\n",
40+
" ar = 2*(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))\n",
41+
" for i in range(3):\n",
42+
" for j in range(i,3):\n",
43+
" A[e[i],e[j]] += (beta[i]*beta[j] + gama[i]*gama[j])/ar\n",
44+
" if (i != j): A[e[j],e[i]] = A[e[i],e[j]] # symmetry\n",
45+
" return A\n",
46+
" \n",
47+
"L, N = 1.0, 40 # length of square, number of intervals\n",
48+
"node, elm, bp, ip, bv, x, y = mesh(L, N) # generate mesh \n",
49+
"ip.sort() # sort ip, just in case \n",
50+
"A, b = fem_mat(node,elm), np.zeros(len(ip)) # build matrices\n",
51+
"for j in range(len(ip)): \n",
52+
" b[j] = np.dot(A[ip[j], bp], bv) # boundary condition \n",
53+
"\n",
54+
"A = np.delete(A, bp, axis=0) # delete rows specified by bp \n",
55+
"A = np.delete(A, bp, axis=1) # delete cols specified by bp \n",
56+
"u = solve(A, -b) # solve \n",
57+
"\n",
58+
"u = np.concatenate((u, bv)) # combine internal+boundary values \n",
59+
"all = np.concatenate((ip, bp)) # internal+boundary nodes\n",
60+
"idx = np.argsort(all) # index sort nodes \n",
61+
"u = np.take(u, idx) # now u[n] is the value at node n \n",
62+
"u = np.reshape(u, (N+1, N+1)) # reshape grid for graphing \n",
63+
"x, y = np.meshgrid(x, y) \n",
64+
"\n",
65+
"plt.figure()\n",
66+
"ax = plt.subplot(111, projection='3d')\n",
67+
"ax.plot_surface(x, y, u, rstride=1, cstride=1,\n",
68+
" linewidth=0, cmap=plt.cm.jet)\n",
69+
"ax.set_xlabel('x'), ax.set_ylabel('y'), ax.set_zlabel('V')\n",
70+
"plt.figure()\n",
71+
"plt.subplot(111, aspect='equal')\n",
72+
"plt.contour(x, y, u, 26)\n",
73+
"plt.show()"
74+
]
75+
},
76+
{
77+
"cell_type": "code",
78+
"execution_count": null,
79+
"metadata": {
80+
"collapsed": true
81+
},
82+
"outputs": [],
83+
"source": []
84+
}
85+
],
86+
"metadata": {
87+
"kernelspec": {
88+
"display_name": "Python 3",
89+
"language": "python",
90+
"name": "python3"
91+
},
92+
"language_info": {
93+
"codemirror_mode": {
94+
"name": "ipython",
95+
"version": 3
96+
},
97+
"file_extension": ".py",
98+
"mimetype": "text/x-python",
99+
"name": "python",
100+
"nbconvert_exporter": "python",
101+
"pygments_lexer": "ipython3",
102+
"version": "3.6.1"
103+
}
104+
},
105+
"nbformat": 4,
106+
"nbformat_minor": 2
107+
}

‎ch07/Program_7.4_laplacerbf.ipynb

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"#\n",
10+
"# Program 7.4: Laplace equation by RBF (laplacerbf.ipynb)\n",
11+
"# J Wang, Computational modeling and visualization with Python\n",
12+
"#\n",
13+
"\n",
14+
"from scipy.linalg import solve\n",
15+
"import matplotlib.pyplot as plt, numpy as np\n",
16+
"%matplotlib notebook\n",
17+
"\n",
18+
"def initialize(L, R, N=22): # calc nodes and boundary values\n",
19+
" inode, bnode, b = [], [], [] # intnl, bndry nodes, bndry values\n",
20+
" x, y = np.linspace(0,L,N+1), np.linspace(0,L,N+1) # same x,y grids\n",
21+
" for j in range(N+1):\n",
22+
" for i in range(N+1):\n",
23+
" r2 = (x[i]-L/2.)**2 + (y[j]-L/2.)**2\n",
24+
" if (i==0 or i==N or j==0 or j==N or r2<=R*R):\n",
25+
" bnode.append([x[i], y[j]]) # bndry node \n",
26+
" b.append(1. if r2<=R*R else 0.)\n",
27+
" else: \n",
28+
" inode.append([x[i], y[j]])\n",
29+
" ni, node = len(inode), inode+bnode # combine nodes \n",
30+
" return np.array(node), [0.]*ni+b, ni, len(node)\n",
31+
" \n",
32+
"def phi(i, r): # Gaussian $\\phi_i$ and $\\nabla^2 \\phi_i$\n",
33+
" r2, e2 = (node[i,0]-r[0])**2 + (node[i,1]-r[1])**2, eps*eps\n",
34+
" f = np.exp(-e2*r2)\n",
35+
" return f, 4*e2*f*(e2*r2-1) \n",
36+
"\n",
37+
"def rbf_mat(ni, nt): # fills matrix $A_{ij} = \\nabla^2_j \\phi_i$ or $\\phi_{ij}$\n",
38+
" A = np.zeros((nt, nt))\n",
39+
" for j in range(nt):\n",
40+
" f, df = phi(np.arange(nt), node[j]) # vector operation \n",
41+
" A[j,:] = df if j < ni else f\n",
42+
" return A \n",
43+
"\n",
44+
"def usum(a, x, y): # calc u at (x,y)\n",
45+
" u = 0.\n",
46+
" for i in range(len(a)):\n",
47+
" u += a[i]*phi(i, [x, y])[0] # [0]= first return value \n",
48+
" return u\n",
49+
" \n",
50+
"L, R, eps = 1.0, 0.25, 20. # sizes of square, disk; shape parameter\n",
51+
"\n",
52+
"node, b, ni, nt = initialize(L, R) # ni, nt = num. intrnl/tot nodes\n",
53+
"A = rbf_mat(ni, nt)\n",
54+
"ai = solve(A, b) # solve \n",
55+
"\n",
56+
"x = np.linspace(0, L, 41) # same x,y plotting grids\n",
57+
"x, y = np.meshgrid(x, x)\n",
58+
"u = usum(ai, x, y)\n",
59+
"\n",
60+
"plt.figure()\n",
61+
"img = plt.imshow(u, cmap=plt.cm.jet) # plot as image \n",
62+
"plt.colorbar(img), plt.axis('off')\n",
63+
"plt.show()"
64+
]
65+
},
66+
{
67+
"cell_type": "code",
68+
"execution_count": null,
69+
"metadata": {
70+
"collapsed": true
71+
},
72+
"outputs": [],
73+
"source": []
74+
}
75+
],
76+
"metadata": {
77+
"kernelspec": {
78+
"display_name": "Python 3",
79+
"language": "python",
80+
"name": "python3"
81+
},
82+
"language_info": {
83+
"codemirror_mode": {
84+
"name": "ipython",
85+
"version": 3
86+
},
87+
"file_extension": ".py",
88+
"mimetype": "text/x-python",
89+
"name": "python",
90+
"nbconvert_exporter": "python",
91+
"pygments_lexer": "ipython3",
92+
"version": "3.6.1"
93+
}
94+
},
95+
"nbformat": 4,
96+
"nbformat_minor": 2
97+
}

‎ch07/Program_7.8_dipole.ipynb

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"#\n",
10+
"# Program 7.8: Angular distribution of dipole radiation (dipole.ipynb)\n",
11+
"# J Wang, Computational modeling and visualization with Python\n",
12+
"#\n",
13+
"\n",
14+
"import numpy as np\n",
15+
"import matplotlib.pyplot as plt\n",
16+
"from mpl_toolkits.mplot3d import Axes3D\n",
17+
"%matplotlib notebook\n",
18+
"\n",
19+
"theta = np.linspace(0, np.pi, 21)\n",
20+
"phi = np.linspace(0, 2*np.pi, 41)\n",
21+
"theta, phi = np.meshgrid(phi, theta)\n",
22+
"x = np.sin(theta)*np.cos(phi)\n",
23+
"y = np.sin(theta)*np.sin(phi)\n",
24+
"z = np.sin(theta)*np.cos(theta) # z=radial times cos(theta) \n",
25+
"\n",
26+
"plt.figure()\n",
27+
"ax = plt.subplot(111, projection='3d', aspect=.5)\n",
28+
"ax.plot_surface(x, y, z, rstride=1, cstride=1, color='w')\n",
29+
"plt.axis('off')\n",
30+
"plt.show()"
31+
]
32+
},
33+
{
34+
"cell_type": "code",
35+
"execution_count": null,
36+
"metadata": {
37+
"collapsed": true
38+
},
39+
"outputs": [],
40+
"source": []
41+
}
42+
],
43+
"metadata": {
44+
"kernelspec": {
45+
"display_name": "Python 3",
46+
"language": "python",
47+
"name": "python3"
48+
},
49+
"language_info": {
50+
"codemirror_mode": {
51+
"name": "ipython",
52+
"version": 3
53+
},
54+
"file_extension": ".py",
55+
"mimetype": "text/x-python",
56+
"name": "python",
57+
"nbconvert_exporter": "python",
58+
"pygments_lexer": "ipython3",
59+
"version": "3.6.1"
60+
}
61+
},
62+
"nbformat": 4,
63+
"nbformat_minor": 2
64+
}

0 commit comments

Comments
 (0)