heatsink class and model-heatsink.py

The 3-dimensional plot that "model-heatsink.py" draws didn't show the fins; this was due to a slight problem in the HeatSink class, not in "model-heatsink.py" itself. The cube was created using "np.empty()", and never got initialized to the "surround" temperature. Most or all of the empty points happen to wind up with values higher than the "surround" value (not implausible, just bad luck). Thus, when the "coords" list is built up, every point in the cube is added to the list.

As this code shows, the fix is simply to create the cube using "np.full()" instead, filling it with "surround". Then the "coords" list only includes the points added by "form_base()" and "form_fins()".

heatsink.py

# -*- coding: utf-8 -*-
"""
Created on Thu Apr 15 09:35:33 2021

@author: rmontant
"""
import numpy as np

def pdf(dx):
    return ( np.exp(-(dx**2)/2) / np.sqrt(2*np.pi) )
#--------

class HeatSink:
    surround = -1e100
    initial_temperature = 20

    def __init__(self, nfins=4, finwidth=4, gapwidth=6, baseheight=5):
        self.nfins = nfins
        self.finwidth = finwidth
        self.gapwidth = gapwidth
        self.baseheight = baseheight
        self.ngaps = self.nfins - 1

        self.width = nfins * finwidth + self.ngaps * gapwidth
        self.length = self.width
        self.height = int( 1.2 * self.width ) # 1.2 seems like a good number

        ###
        ### 2021-04-22 afternoon:
        ### The cube MUST be initialized to the "surround" value here!
        ### (because it never happens anywhere else...)
        ###
        self.cube = np.full(
            (self.width + 2, self.length + 2, self.height + 2), self.surround )

        self.form_base()
        self.form_fins()
        self.build_pointarray()

        self.tmax = self.initial_temperature
    #--------

    def form_base(self):
        self.cube[1:self.width-1, 1:self.length-1, 1:self.baseheight+1] \
            = self.initial_temperature

    def form_fins(self):
        for f in range(self.nfins):
            fstart = 1 + f * (self.gapwidth + self.finwidth)
            fend = fstart + self.finwidth
            self.cube[ \
                fstart:fend, 1:self.length-1 , self.baseheight+1:self.height-1 \
            ] = self.initial_temperature
    #--------

    def build_pointarray(self):
        coords = []
        for x in range(self.width):
            for y in range(self.length):
                for z in range(self.height):
                    if self.cube[x,y,z] > self.surround:
                        coords.append( (x,y,z) )
        self.coords = np.array( coords )

    def set_hotspot(self, w, l, t, scale):
        rmax = ((self.width**2 + self.length**2)**0.5) * scale
        for x in range(1, self.width-1):
            for y in range(1, self.length-1):
                r = ( (w-x)**2 + (l-y)**2 )**0.5
                s = pdf(r/rmax) / pdf(0)

                self.cube[x,y,0] = s * t \
                    + (1-s) * self.cube[x,y,0]
        if t > self.tmax:
            self.tmax = t
#--------

model-heatsink.py

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 20 10:06:10 2021

@author: rmontant
"""
import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import Normalize
from mpl_toolkits.mplot3d import Axes3D

import heatsink

def draw_contour(axis, surface, tmin=None, tmax=None):
    if not tmin:
        tmin = np.min(surface)
    if not tmax:
        tmax = np.max(surface)
    normal = Normalize(tmin, tmax)
    pcolormap = axis.pcolormesh(surface, cmap=cm.afmhot, norm=normal)
    axis.contour(surface, cmap=cm.hot)
    plt.colorbar(pcolormap, ax=axis)
#--------

def rgb(v):
    return (np.sqrt(v), v, 0.9*v, 0.05)
#--------

def draw_heatsink(axis, heatsink):
    xcoords = [ p[0] for p in heatsink.coords]
    ycoords = [ p[1] for p in heatsink.coords]
    zcoords = [ p[2] for p in heatsink.coords]
    colors = [ rgb(0.6 * z / heatsink.height) for z in zcoords]
    axis.scatter(xcoords, ycoords, zcoords, color=colors, marker=',')
    axis.set_xlabel('width')
    axis.set_ylabel('length')
    axis.set_zlabel('height')
#--------

def main(argv=[__name__]):
    htsnk = heatsink.HeatSink(nfins=5, finwidth=6, gapwidth=9, baseheight=10)

    nhotspots = 3
    for i in range(nhotspots):
        temperature = 100 + 10 * np.random.random()
        w = np.random.randint(1, htsnk.width-1)
        l = np.random.randint(1, htsnk.length-1)
        scale = 0.07 + 0.06 * np.random.random()
        print('hotspot {:d} is {:.3} at ({:d}, {:d})' \
              .format(i, scale, w, l))
        htsnk.set_hotspot(w, l, temperature, scale)

    fig = plt.figure(figsize=(8,8))
    a1 = fig.add_subplot(2, 2, 1)
    draw_contour(a1, htsnk.cube[:,:,0], \
                 tmin=htsnk.initial_temperature,
                 tmax=htsnk.tmax
    )

    a2 = fig.add_subplot(2,2, 2)
    surface = htsnk.cube[:, htsnk.length//2, :]
    draw_contour(a2, surface,
                 tmin=htsnk.initial_temperature,
                 tmax=htsnk.tmax,
    )

    a4 = fig.add_subplot(2,2, 4, projection='3d')
    draw_heatsink(a4, htsnk)

    prows, pcols = 4, 4
    fig2, axes = plt.subplots(nrows=prows, ncols=pcols, figsize=(8,8))
    for i in range(prows):
        for j in range(pcols):
            zval = int( htsnk.height * (i*pcols + j) / (prows*pcols) )
            surface = htsnk.cube[ :, :, zval]
            draw_contour(axes[i][j], surface,
                          tmin=htsnk.initial_temperature,
                          tmax=htsnk.tmax,
            )

    plt.show()
#--------

if __name__ == '__main__':
    sys.exit(main(sys.argv))
#--------