1. # -*- coding: utf-8 -*-
2. """
3. Created on Thu Aug  6 18:27:54 2020
4.  
5. Chapter 6 PDE: Excercise 2: Metallic plate and heat sources
6.  
7. 2D Elliptic PDE
8. """
9.  
10. """
11. Heat distribution in a metallic plate
12. size: 1m x 1m
13. distribution:
14.    d2T/d2x + d2T/d2y = f
15.    f = -Q/k (constant)
16. """
17. from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import
18. import numpy as np
19. import scipy.sparse as sp
20. from scipy.sparse.linalg import spsolve
21. import matplotlib.pyplot as plt
22. from matplotlib import cm
23. from matplotlib.ticker import LinearLocator, FormatStrFormatter
24.  
25.  
26. # define the space
27. x0, xf , nX = 0.0, 1.0, 100
28. y0, yf , nY = 0.0, 1.0, 100
29.  
30. # Vectors
31. X = np.linspace(x0, xf, nX)
32. Y = np.linspace(y0, yf, nY)
33. h = X[1] - X[0]
34. h2 = 1/(h**2)
35. dh = 1/(2*h) # 1/2dx and 1/2dy
36.  
37. # Size of the system
38. N = nX*nY
39. A = sp.lil_matrix((N,N))
40. b = np.zeros(N)
41.  
42. # Right hand side of the equation
43. Q = 100 # kW/m3 heat generation
44. k = 0.01 # W/mC thermal conductivity
45. f = -Q/k
46.  
47. # nested loop for the node mapping
48. for i in range(nX):
49.     for j in range(nY):
50.         # mapping relation center
51.         K = j*nY + i
52.         # neighbours p:plus m:minus
53.         KpN = (j+1)*nY + i              
54.         KpN2 = (j+2)*nY + i            
55.         KmN = (j-1)*nY + i              
56.         KmN2 = (j-2)*nY + i            
57.         Kp1 = j*nY + (i+1)              
58.         Kp2 = j*nY + (i+2)              
59.         Km1 = j*nY + (i-1)              
60.         Km2 = j*nY + (i-2)              
61.        
62.         # Check boundary conditions
63.         # Fill matrix
64.         # K,K is the diagonal value
65.         # Robin Boundary Condition
66.         if (i == 0):
67.             A[K,K] = -3*dh - 5
68.             A[K,Kp1] = 2/h
69.             A[K,Kp2] = -dh
70.             b[K] = -125.0  # degrees Celsius
71.            
72.         elif (i == nX-1):
73.             A[K,K] = -3*dh + 5
74.             A[K,Km1] = 2/h
75.             A[K,Km2] = -dh
76.             b[K] = 125.0  # degrees Celsius
77.            
78.         elif (j == 0):
79.             A[K,K] = -3*dh - 5
80.             A[K,KpN] = 2/h
81.             A[K,KpN2] = -dh
82.             b[K] = -125.0  # degrees Celsius
83.            
84.         elif (j == nY-1):
85.             A[K,K] = -3*dh + 5
86.             A[K,KmN] = 2/h
87.             A[K,KmN2] = -dh
88.             b[K] = 125.0  # degrees Celsius
89.            
90.         # Not Boundary Conditions
91.         else:
92.             A[K,K] = -4*h2
93.             A[K,Kp1] = h2
94.             A[K,Km1] = h2
95.             A[K,KpN] = h2
96.             A[K,KmN] = h2
97.             b[K] = f
98.  
99.  
100.  
101. # nested loop for the Corner mapping
102. # The loop is needed again to use the definition of K
103. for i in range(nX):
104.     for j in range(nY):
105.         # mapping relation
106.         K = j*nY + i
107.         # neighbours p:plus m:minus
108.         KpN = (j+1)*nY + i              
109.         KmN = (j-1)*nY + i              
110.         Kp1 = j*nY + (i+1)              
111.         Km1 = j*nY + (i-1)              
112.        
113.         # Corner 1: 0,0; 2: 1,1; 3: -1,0; 4: -1,-1
114.         if (i == 0 and j == 0):
115.             b[K] = (b[Kp1]+b[KpN]+b[11])/3
116.            
117.         elif (i == (nX-1) and j == 0):
118.             b[K] = (b[Km1]+b[KpN]+b[18])/3
119.            
120.         elif (i == 0 and j == (nY-1)):
121.             b[K] = (b[Kp1]+b[KmN]+b[81])/3
122.            
123.         elif (i == (nX-1) and j == (nY-1)):
124.             b[K] = (b[Km1]+b[KmN]+b[88])/3
125.  
126.  
127. # Convert the LIL matrix to a CSR matrix
128. A = A.tocsr()
129.  
130. # Solve using spsolve
131. z = spsolve(A, b)
132. Z = z.reshape(nX, nY)
133. plt.contourf(X, Y, Z, cmap='inferno')
134. plt.colorbar()
135. plt.xlabel('X')
136. plt.ylabel('Y')
137.  
138. # 3d plot
139. X, Y = np.meshgrid(X, Y)
140. fig = plt.figure()
141. ax = fig.gca(projection='3d')
142. surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1,
143.                        cmap=cm.plasma, linewidth=0, antialiased=False)
144. fig.colorbar(surf, shrink=0.5, aspect=5)
145.  
146.  
147. plt.show()