# Last edited on 2022-09-05 11:51:24 by stolfi

def fusion_path_plan(MVS,nx,ny,nz,dx,dy,dz,dtmax):
  # Given: 
  #   {MVS} = set of all traces to fabricate, by layer.
  #   {nx,ny,nz} = total voxels in object's box, along each axis.
  #   {dx,dy,dz} = voxel size along each axis.
  #   {dtmax} = max time step for accurate enough simulation.
  #   {c} = thermal conductivity of solid metal, {J/s/m^3/K}.
  #   {f} = specific heat capacity of solid metal, {J/m^3/K}.
  
  tsim = 0 # Simulated time.
  nzs = 30 # Max number of layers retained in heat map.
  
  # Tomogram indices are {[iz][iy][ix]}:
  C = allocate_tomogram(nx,ny,nzs) # Voxel occupancy (0 = 100% air, 1 = 100% metal).
  T = allocate_tomogram(nx,ny,nzs) # Temperature of each voxel at time {tsim}.
  L = allocate_tomogram(nx,ny,nzs) # WORK: estimated net energy flow into pixel.
  
  for iz in range(nz):
    MV = MVS[iz] # Traces in layer {iz}.
    EL = chop_traces(MV,dtmax) # Set of those traces chopped into small enough pieces.
    shift_down(T) # Cyclic shift of the {nzs} layers of {T}.
    shift_down(C) # Cyclic shift of the {nzs} layers of {C}.
    izs = nzs-1 # Layer of tomograms corresponding to layer {iz} of object.
    fill_layer(C, izs, 0)   # Set current layer to all air.
    fill_layer(T, izs, NaN) # Temperature of empty voxels is irrelevant.
    while not size(EL) == 0:
      el = select_trace(EL, C, T, izs); # Select next trace, based on heat map.
      simulate_deposition(el, C,T, dx,dy,dz, izs); # Add molten metal to voxels of {C} and {T} along {el}.
      dt = duration(el) # Time taken to fabricate {el}.
      simulate_heat_transfer(C,T,L, dx,dy,dz, c,f, dt) # Integrate the discretized heat equation.
      tsim = tsim + dt # Advance the simulated time.
      
def simulate_heat_transfer(C,T,L, dx,dy,dz, c,f, dt):
  # Computes the heat flow and updates {T} of filled voxels over the time step {dt}.
  # Uses {L} as work area.
  
  nxs, nys, nzs = tomogram_size(C)
  
  # Compute the net energy accumulation per volume (J/s/m^3) in each voxel:
  for izs in range(nzs):
    for iy in range(nys):
      for ix in range(nxs):
        L[izs][iys][ixs] = net_energy_flow(C,T, ixs,iys,izs, dx,dy,dz, c) # Net energy accum in voxel.

  # Update the temperature of each voxel:
  for izs in range(nzs):
    for iys in range(ny):
      for ixs in range(nx):
        if C[izs][iys][ixs] > 0:
          T[izs][iys][ixs] = T[izs][iys][ixs] + L[izs][iys][ixs]/f*dt
        else:
          T[izs][iys][ixs] = NaN

def net_energy_flow(C,T, ixs,iys,izs, dx,dy,dz, c):
  # Computes the net energy accumulation per volume (J/s/m^3) at the center of
  # the voxel with indices {[izs][iys][ixs]}.
  # Result is relevant only if {C[izs][iys][ixs]} is nonzero.
  
  Ltot = 0
  if C[izs][iys][ixs] > 0:
    nxs, nys, nzs = tomogram_size(C)
    Ti = T[izs][iys][ixs] # Temperature at center of voxel.
    for axis in range(3):
      for dir in (-1, +1):
        # Estimate net energy flow per volume across face perp to {axis} in direction {dir}:
        jxyzs = [ixs, iys, izs]
        jxyzs[axis] = jxyzs[axis] - dir
        jxs, jys, jzs = jxyzs
        if jxs >= 0 and jxs < nxs and jys >= 0 and jys < nys and jzs >= 0 and jzs < nzs and C[jzs][jys][jxs] > 0:
          Tj = T[jzs][jys][jxs]
          dL = c*(Tj - Ti) # Energy flow across face, divided by voxel volume.
          Ltot = Ltot + dL
  return Ltot