(Disclaimer: The physics here may not be totally sound. Fluxes are just differences between triangle temperatures times edge lengths.)
In [ ]:
%pylab inline
In [ ]:
import numpy as np
#import time
import matplotlib.pyplot as plt
from IPython.display import clear_output
#from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
#import math
#import pylab
#import copy
#import JMesh
In [ ]:
# Setting up a 2D mesh and boundary conditions
vert = [ [0, 0], [0, 1], [0.25, 0.5], [0.5, 0], [1, 1], [1, 0] ] # (x, y) pairs
tri = [ [0, 2, 1], [0, 3, 2], [1, 2, 4], [2, 3, 4], [3, 5, 4] ] # Vertex indices
nbr = [ [2, 1, 2, -1], [2, -1, 3, 0], [2, 0, 3, -1], [3, 1, 4, 2], [1, -1, -1, 3] ]; # Count, neighbour triangle (-1 if none.) Note: CCW indexing from corner 0.
In [ ]:
# To refine the handwritten mesh defined above a little bit
"""
random.seed(0)
import JMesh
#JMesh.showTemp(vert, tri, np.random.rand(len(tri)))
for i in range(3):
JMesh.refineAll(vert, tri, nbr)
JMesh.improveBySwapping(vert, tri, nbr, 10000)
#JMesh.showTemp(vert, tri, np.random.rand(len(tri)))
print "vert=", vert
print "tri=", tri
print "nbr=", nbr
""";
In [ ]:
vert= [[0, 0], [0, 1], [0.25, 0.5], [0.5, 0], [1, 1], [1, 0], [0.8333333333333333, 0.3333333333333333], [0.5833333333333333, 0.5], [0.41666666666666663, 0.8333333333333333], [0.25, 0.16666666666666666], [0.08333333333333333, 0.5], [0.6111111111111109, 0.5555555555555555], [0.75, 0.7222222222222221], [0.41666666666666663, 0.611111111111111], [0.1111111111111111, 0.2222222222222222], [0.19444444444444445, 0.38888888888888884], [0.25, 0.611111111111111], [0.16666666666666663, 0.7777777777777777], [0.6388888888888888, 0.27777777777777773], [0.36111111111111105, 0.3888888888888889], [0.4444444444444444, 0.2222222222222222], [0.027777777777777776, 0.5], [0.25, 0.05555555555555555], [0.47222222222222215, 0.9444444444444443], [0.9444444444444443, 0.4444444444444444], [0.7777777777777777, 0.1111111111111111]]
tri= [[2, 10, 15], [10, 2, 16], [1, 0, 21], [0, 3, 22], [4, 1, 23], [5, 4, 24], [3, 5, 25], [16, 8, 17], [21, 0, 14], [18, 6, 11], [11, 7, 18], [14, 9, 15], [15, 10, 21], [21, 14, 15], [13, 2, 19], [19, 7, 13], [13, 7, 11], [11, 8, 13], [17, 1, 21], [16, 17, 21], [21, 10, 16], [19, 2, 15], [15, 9, 19], [17, 8, 23], [23, 1, 17], [22, 3, 20], [20, 9, 22], [16, 2, 13], [13, 8, 16], [24, 6, 25], [25, 5, 24], [24, 4, 12], [18, 3, 25], [25, 6, 18], [12, 4, 23], [23, 8, 12], [19, 9, 20], [20, 3, 18], [18, 7, 19], [19, 20, 18], [12, 8, 11], [11, 6, 24], [24, 12, 11], [14, 0, 22], [22, 9, 14]]
nbr= [[3, 1, 12, 21], [3, 0, 27, 20], [2, -1, 8, 18], [2, -1, 25, 43], [2, -1, 24, 34], [2, -1, 31, 30], [2, -1, 30, 32], [3, 28, 23, 19], [3, 2, 43, 13], [3, 33, 41, 10], [3, 16, 38, 9], [3, 44, 22, 13], [3, 0, 20, 13], [3, 8, 11, 12], [3, 27, 21, 15], [3, 38, 16, 14], [3, 15, 10, 17], [3, 40, 28, 16], [3, 24, 2, 19], [3, 7, 18, 20], [3, 12, 1, 19], [3, 14, 0, 22], [3, 11, 36, 21], [3, 7, 35, 24], [3, 4, 18, 23], [3, 3, 37, 26], [3, 36, 44, 25], [3, 1, 14, 28], [3, 17, 7, 27], [3, 41, 33, 30], [3, 6, 5, 29], [3, 5, 34, 42], [3, 37, 6, 33], [3, 29, 9, 32], [3, 31, 4, 35], [3, 23, 40, 34], [3, 22, 26, 39], [3, 25, 32, 39], [3, 10, 15, 39], [3, 36, 37, 38], [3, 35, 17, 42], [3, 9, 29, 42], [3, 31, 40, 41], [3, 8, 3, 44], [3, 26, 11, 43]]
len(tri)
In [ ]:
#vert= [[0, 0], [0, 1], [0.25, 0.5], [0.5, 0], [1, 1], [1, 0], [0.8333333333333333, 0.3333333333333333], [0.5833333333333333, 0.5], [0.41666666666666663, 0.8333333333333333], [0.25, 0.16666666666666666], [0.08333333333333333, 0.5], [0.6111111111111109, 0.5555555555555555], [0.75, 0.7222222222222221], [0.41666666666666663, 0.611111111111111], [0.1111111111111111, 0.2222222222222222], [0.19444444444444445, 0.38888888888888884], [0.25, 0.611111111111111], [0.16666666666666663, 0.7777777777777777], [0.6388888888888888, 0.27777777777777773], [0.36111111111111105, 0.3888888888888889], [0.4444444444444444, 0.2222222222222222], [0.027777777777777776, 0.5], [0.25, 0.05555555555555555], [0.47222222222222215, 0.9444444444444443], [0.9444444444444443, 0.4444444444444444], [0.7777777777777777, 0.1111111111111111], [0.20370370370370372, 0.14814814814814814], [0.12037037037037036, 0.09259259259259259], [0.7685185185185184, 0.574074074074074], [0.7962962962962961, 0.4444444444444444], [0.5925925925925924, 0.7037037037037035], [0.4814814814814814, 0.2962962962962963], [0.5277777777777777, 0.38888888888888884], [0.5277777777777777, 0.16666666666666663], [0.35185185185185175, 0.25925925925925924], [0.5462962962962963, 0.8333333333333333], [0.7407407407407407, 0.8888888888888888], [0.75, 0.2407407407407407], [0.6388888888888888, 0.1296296296296296], [0.898148148148148, 0.7222222222222221], [0.9074074074074072, 0.18518518518518517], [0.8518518518518517, 0.2962962962962963], [0.36111111111111105, 0.6851851851851851], [0.3055555555555555, 0.5740740740740741], [0.31481481481481477, 0.14814814814814814], [0.39814814814814814, 0.09259259259259259], [0.21296296296296294, 0.9074074074074072], [0.35185185185185175, 0.8518518518518517], [0.2685185185185185, 0.31481481481481477], [0.2685185185185185, 0.4259259259259259], [0.12037037037037036, 0.537037037037037], [0.14814814814814814, 0.6296296296296295], [0.0648148148148148, 0.7592592592592592], [0.4814814814814814, 0.6666666666666666], [0.5370370370370369, 0.5555555555555555], [0.4537037037037036, 0.5], [0.34259259259259256, 0.5], [0.11111111111111112, 0.37037037037037035], [0.10185185185185186, 0.4629629629629629], [0.18518518518518517, 0.2592592592592592], [0.6111111111111109, 0.4444444444444443], [0.6944444444444443, 0.38888888888888884], [0.046296296296296294, 0.24074074074074073], [0.27777777777777773, 0.7407407407407406], [0.7592592592592592, 0.037037037037037035], [0.9814814814814814, 0.48148148148148145], [0.4907407407407407, 0.9814814814814814], [0.25, 0.018518518518518517], [0.009259259259259259, 0.5], [0.19444444444444442, 0.537037037037037], [0.17592592592592593, 0.4629629629629629]]
#tri= [[12, 11, 28], [6, 24, 29], [11, 12, 30], [23, 8, 35], [24, 6, 41], [16, 2, 43], [8, 23, 47], [21, 10, 50], [7, 11, 54], [10, 21, 58], [11, 7, 60], [3, 5, 64], [5, 4, 65], [4, 1, 66], [0, 3, 67], [1, 0, 68], [2, 16, 69], [63, 16, 42], [35, 12, 36], [43, 13, 42], [42, 16, 43], [70, 2, 69], [53, 13, 54], [38, 18, 33], [33, 3, 38], [36, 4, 66], [66, 23, 35], [35, 36, 66], [43, 2, 56], [56, 13, 43], [65, 4, 39], [48, 15, 59], [59, 9, 48], [52, 17, 46], [46, 1, 52], [55, 13, 56], [56, 19, 55], [51, 17, 52], [33, 18, 31], [31, 20, 33], [30, 53, 54], [54, 11, 30], [61, 18, 37], [37, 6, 61], [45, 20, 44], [40, 5, 65], [54, 13, 55], [55, 7, 54], [38, 25, 37], [37, 18, 38], [27, 22, 26], [26, 14, 27], [67, 3, 45], [32, 7, 55], [55, 19, 32], [64, 25, 38], [38, 3, 64], [67, 22, 27], [27, 0, 67], [53, 8, 42], [42, 13, 53], [59, 15, 57], [33, 20, 45], [45, 3, 33], [31, 19, 34], [34, 20, 31], [68, 0, 62], [26, 22, 44], [44, 9, 26], [59, 14, 26], [26, 9, 59], [68, 62, 57], [36, 12, 39], [39, 4, 36], [60, 7, 32], [67, 45, 44], [44, 22, 67], [63, 17, 51], [51, 16, 63], [70, 69, 50], [62, 0, 27], [27, 14, 62], [61, 60, 32], [32, 18, 61], [61, 6, 29], [49, 19, 56], [56, 2, 49], [30, 12, 35], [59, 57, 62], [62, 14, 59], [57, 15, 70], [70, 58, 57], [65, 24, 41], [41, 40, 65], [51, 50, 69], [69, 16, 51], [66, 1, 46], [60, 61, 29], [64, 5, 40], [40, 25, 64], [48, 19, 49], [49, 15, 48], [58, 70, 50], [50, 10, 58], [46, 17, 63], [63, 47, 46], [47, 23, 66], [66, 46, 47], [70, 15, 49], [49, 2, 70], [34, 9, 44], [44, 20, 34], [68, 57, 58], [58, 21, 68], [31, 18, 32], [32, 19, 31], [34, 19, 48], [48, 9, 34], [41, 6, 37], [28, 11, 60], [60, 29, 28], [37, 25, 40], [40, 41, 37], [35, 8, 53], [53, 30, 35], [42, 8, 47], [47, 63, 42], [52, 1, 68], [39, 12, 28], [65, 39, 28], [51, 52, 68], [50, 51, 68], [68, 21, 50], [65, 28, 29], [29, 24, 65]]
#nbr= [[3, 2, 119, 128], [3, 4, 134, 84], [3, 0, 87, 41], [3, 6, 123, 26], [3, 1, 118, 92], [3, 16, 28, 20], [3, 3, 106, 125], [3, 9, 103, 132], [3, 10, 41, 47], [3, 7, 113, 103], [3, 8, 74, 119], [2, -1, 98, 56], [2, -1, 30, 45], [2, -1, 96, 25], [2, -1, 52, 58], [2, -1, 66, 127], [3, 5, 95, 21], [3, 78, 20, 126], [3, 87, 72, 27], [3, 29, 60, 20], [3, 17, 5, 19], [3, 109, 16, 79], [3, 60, 46, 40], [3, 49, 38, 24], [3, 63, 56, 23], [3, 73, 13, 27], [3, 106, 3, 27], [3, 18, 25, 26], [3, 5, 86, 29], [3, 35, 19, 28], [3, 12, 73, 129], [3, 101, 61, 32], [3, 70, 117, 31], [3, 37, 104, 34], [3, 96, 127, 33], [3, 46, 29, 36], [3, 85, 54, 35], [3, 77, 33, 130], [3, 23, 114, 39], [3, 65, 62, 38], [3, 124, 22, 41], [3, 8, 2, 40], [3, 83, 49, 43], [3, 118, 84, 42], [3, 62, 111, 75], [3, 98, 12, 93], [3, 22, 35, 47], [3, 53, 8, 46], [3, 55, 121, 49], [3, 42, 23, 48], [3, 57, 67, 51], [3, 69, 81, 50], [3, 14, 63, 75], [3, 74, 47, 54], [3, 36, 115, 53], [3, 99, 48, 56], [3, 24, 11, 55], [3, 76, 50, 58], [3, 80, 14, 57], [3, 123, 125, 60], [3, 19, 22, 59], [3, 31, 90, 88], [3, 39, 44, 63], [3, 52, 24, 62], [3, 115, 116, 65], [3, 111, 39, 64], [3, 15, 80, 71], [3, 50, 76, 68], [3, 110, 70, 67], [3, 89, 51, 70], [3, 68, 32, 69], [3, 66, 88, 112], [3, 18, 128, 73], [3, 30, 25, 72], [3, 10, 53, 82], [3, 52, 44, 76], [3, 67, 57, 75], [3, 104, 37, 78], [3, 95, 17, 77], [3, 21, 94, 102], [3, 66, 58, 81], [3, 51, 89, 80], [3, 97, 74, 83], [3, 114, 42, 82], [3, 43, 1, 97], [3, 100, 36, 86], [3, 28, 109, 85], [3, 2, 18, 124], [3, 61, 71, 89], [3, 81, 69, 88], [3, 61, 108, 91], [3, 102, 112, 90], [3, 134, 4, 93], [3, 122, 45, 92], [3, 131, 79, 95], [3, 16, 78, 94], [3, 13, 34, 107], [3, 82, 84, 120], [3, 11, 45, 99], [3, 121, 55, 98], [3, 116, 85, 101], [3, 108, 31, 100], [3, 91, 79, 103], [3, 7, 9, 102], [3, 33, 77, 105], [3, 126, 107, 104], [3, 6, 26, 107], [3, 96, 105, 106], [3, 90, 101, 109], [3, 86, 21, 108], [3, 117, 68, 111], [3, 44, 65, 110], [3, 71, 91, 113], [3, 9, 132, 112], [3, 38, 83, 115], [3, 54, 64, 114], [3, 64, 100, 117], [3, 32, 110, 116], [3, 4, 43, 122], [3, 0, 10, 120], [3, 97, 133, 119], [3, 48, 99, 122], [3, 93, 118, 121], [3, 3, 59, 124], [3, 40, 87, 123], [3, 59, 6, 126], [3, 105, 17, 125], [3, 34, 15, 130], [3, 72, 0, 129], [3, 30, 128, 133], [3, 37, 127, 131], [3, 94, 130, 132], [3, 113, 7, 131], [3, 129, 120, 134], [3, 1, 92, 133]]
#len(tri)
In [ ]:
# Setting up boundary conditions to be constant all around the edges
def setBoundary(t, n, u, boundaryValue):
for i in range(len(t)):
if n[i][0] != 3:
# This is a boundary triangle
u[i] = boundaryValue
In [ ]:
def initializeSolution(t, n, u, initValue):
for i in range(len(t)):
if n[i][0] == 3:
u[i] = initValue
In [ ]:
def showTemp(vert, tri, u, counter):
fig=plt.figure()
ax=fig.add_subplot(111)
title("frame #" + str(counter))
patches = []
for i in range(len(tri)):
corners = vert[tri[i][0]], vert[tri[i][1]], vert[tri[i][2]]
polygon = Polygon(corners, True)
patches.append(polygon)
colors = u
p = PatchCollection(patches, cmap=matplotlib.cm.jet, alpha=0.7)
p.set_array(np.array(colors))
ax.add_collection(p)
plt.colorbar(p)
plt.show()
In [ ]:
u = np.zeros( len(tri) )
initializeSolution(tri, nbr, u, 0.0)
setBoundary(tri, nbr, u, 0.5)
showTemp(vert, tri, u, 0)
In [ ]:
# Heat diffusion constant
kappa = 0.3
# Maximum time step size (constrained by CFL)
#dt = (dx*dx)/(4*kappa) # ok when dx==dy
#dt = (dx*dx*dy*dy)/(dx*dx+dy*dy) / (2*kappa) # if dx!=dy
dt = 10.0
In [ ]:
# Should really be called edgeLength...
def faceArea(vert, tr, face): # List of vertices, a single triangle, face index in triangle
v0 = np.array( vert[ tr[face] ] ) # So that we can apply subtraction to the two lists
v1 = np.array( vert[ tr[(face+1)%3] ] )
#print "v0 =", v0, ", v1 =", v1
#print "diff =", v1-np.array(v0)
return math.sqrt(pylab.dot(v1-v0, v1-v0))
In [ ]:
def assembleA(vert, tri, nbr):
n = len(tri)
A = np.zeros( (n, n) )
# Assemble the matrix one cell at a time.
for i in range(n):
A[i, i] = 1;
# Looping over neighbours.
#print "i =", i, ", #neighbours =", nbr[i][0]
for k in range(3):
j = nbr[i][1+k]
if (j!=-1):
# Flux across face between the two cells:
f = faceArea(vert, tri[i], k)*dt*kappa
A[i, i] += f
A[i, j] -= f
# Apply Dirichlet boundary conditions, # u_i^n = u_i^n-1, for i=0 and i=n-1.
for i in range(n):
if nbr[i][0] != 3:
#print "Resetting row ", i
A[ i, : ] = np.zeros(n)
A[ i, i ] = 1
return A;
In [ ]:
pl2D = False # 2D image plot
In [ ]:
A = assembleA(vert, tri, nbr)
# Visualizing A and inv(A):
if False:
fig2 = plt.figure(figsize=(16, 4))
plt.subplot(121); # Denne går inn i current figure, som nå er fig2?
plt.imshow(A, interpolation='nearest');
plt.colorbar();
plt.subplot(122)
plt.imshow(pylab.inv(A), interpolation='nearest')
qqq=plt.colorbar() # Trenger å putte den i en variabel for å unngå utskrift av "handle"
#plt.sca(ax) # Setting current axis to 'ax' and current figure to the parent figure of 'ax'.
#plt.figure(fig.number); # Setting the current figure
u = np.zeros( len(tri) )
initializeSolution(tri, nbr, u, 0.0)
setBoundary(tri, nbr, u, 1.0)
iterations = 50
if not(pl2D):
fig, ax = plt.subplots()
for k in range(iterations):
v = np.linalg.solve(A, u);
if not(pl2D):
n = len(tri)
x = np.linspace(0.0, n-1, n) + 0.5
plot(x, u, 'r.-', label='$u^{k+1}$')
pylab.ylim([-0.5, 1.5])
#time.sleep(0.2)
clear_output()
display(fig)
ax.cla()
title("Time = " + str(dt*k))
legend()
else:
clear_output()
vv = v[0]
v[0] = 0.0 # Just to fool the image colour scaling
showTemp(vert, tri, v, k)
v[0] = vv
u, v = v, u;
if not(pl2D):
plt.close()
1D plot above: Note the boundary triangles to the left, they are kept constant. The rest of the points are more or less randomly ordered triangles in the mesh
2D plot above: Triangle #0 is just used to force the image to be scaled betwen 0 and 1. The rest of the triangles are coloured according to the solution.