Specify Internal Heat Generation at multiple points

I'm currently working through the PDEs toolbox, trying to solve the heat conduction equation for a 2D surface with a diamond-shaped hole at the middle. Heat-generation is supposed to happen only along the edges of the diamond. I have a script that finds the coordinates for the nodes along the edges of the diamond. However, I am struggling to find a way to implement heat-generation at specific nodes in a matrix.

I have an array of the nodes where the heat-generation is meant to occur as well as an array specifying the heat-generation over position and time (heat generation is higher at the wide tips of the diamond).

Essentially, I am looking to specify the internal heating source at each point along the edge of the diamond. I have the time-dependent heat generation vector for each point along the edge of the diamond and the xy-coordinates for each node along the diamond's edge. How can I put this into the PDEs Toolbox internalHeatSource function?

I am using MATLAB R2017a and the PDEs Toolbox

If possible, would someone be able to explain it in these steps so I can understand the process better?

1. Constant internal heat generation at a point (x,y)
2. Time-dependent heat generation at a point (x,y) where we have the numerical data for the heat-gen based on the input time
3. Time-dependent heat generation at two points (the wide tips) given by the coordinates (x1,y1) and (x2,y2) where we have the numerical data for the heat-gen based on the input time.
4. Time-dependent heat generation at all points along the diamond edge where we have the time-dependent heat generation numerical data at each node along the edge and the (x,y) coordinates of each node.

Here's my code: EdgeNodes holds the (x,y) coordinates of all nodes along the edge of the diamond, while qdot is the heat generation over time correlating to each EdgeNode.

%This script executes the whole crack-healing analysis.
close all
clear all
clc

global a %[m] Crack half-length
global theta %[deg] Crack tip half-degree
global L %[m] Plate length
global W %[m] Plate width
global B %[kg/m^3] Mass density
global Cp %[J/kgC] Specific heat capacity
global k %[W/mC] Thermal conductivity
global r %[m] Radial distance from crack tip
global rho %[Ohm-m] Resistivity
global t %[s] Pulse time
global V %[V] Source pulse
global Tamb %[C] Ambient temperature

%For this example, we shall use steel.
a=2e-3;
theta=0.5; %Crack tip half-angle
L=10*a;
W=L; %Square plate
B=7861.2;
Cp=444;
k=64.6;
rho=0.14224e-6;
t=0:0.01:1;
V=100*1./(cosh(pi*t).^2);
Tamb=21;

b=a*tan(theta*pi/180); %[m] Crack maximum half-width.
c=a/cos(theta*pi/180); %[m] Crack hypoteneuse.
r=c/100:c/100:c; %Radius along crack edge

%Calculate KJ term for infinite plate (L>>a).
E=V/L; %[V/m] Electric field created by the pulse.
Jo=(1/rho).*E; %[A/m^2] Calculate applied current density.
KJ=Jo*sqrt(pi*a); %[A/(m^3/2)] Term used in calculation of current density around crack tip.

%%Calculate the current density along the crack edge.
%Now calculate the current density.
J=zeros(length(t),length(r)); %Initialize magnitude of current density.

for i=1:length(t)
    for j=1:length(r)
    J(i,j)=KJ(i)/sqrt(2*pi*r(j)); %[A/m^2] Calculate current density from tip.
    end
end
%Current density is J(time, distance r from tip).

qdotsample=rho.*J.^2; %[W] Internal heat generation due to Joule Heating.
%qdotsample(time,radial distance from crack tip)

thermalmodelT = createpde('thermal','transient');

%Create crack geometry model.
r1 = [3 4 -L/2 L/2 L/2 -L/2 W/2 W/2 -W/2 -W/2];
r2 = [2 4 0 a 0 -a b 0 -b 0];
gdm = [r1; r2]';
g = decsg(gdm,'R1-R2',['R1'; 'R2']');

geometryFromEdges(thermalmodelT,g);

figure
pdegplot(thermalmodelT,'EdgeLabels','on');
axis([-L/2 L/2 -W/2 W/2]);
title 'Crack Geometry with Edge Labels';
axis equal

%Create the mesh.
mesh=generateMesh(thermalmodelT,'Hmax',a/5);
figure
pdeplot(thermalmodelT);
axis equal
title 'Block With Finite Element Mesh Displayed'
[Nodes,Medges,Melements]=meshToPet(thermalmodelT.Mesh);

%Search for nodes within a given box to find nodes along crack edges.
eps=[a/10 b/10];
j=1;
for i=1:length(Nodes)
   if abs(Nodes(1,i))<(a+eps(1)) && abs(Nodes(2,i))<(b+eps(2))
       EdgeNodes(1,j)=Nodes(1,i);
       EdgeNodes(2,j)=Nodes(2,i);
       j=j+1;
   end
end

%Check to make sure gathered nodes are the correct ones.
dil=1.5; %Change how vertically wide the crack appears, for visual purposes only.
figure
scatter(EdgeNodes(1,:),EdgeNodes(2,:),'filled')
axis([-a-eps(1) a+eps(1) -b*dil b*dil]) 
title 'Crack Edge Nodes';
%The output should be a sideways diamond. If there are points beyond the
%diamond edges, then the node search code will need to be altered.

%Create a new heat-gen array based on the extracted nodes. For crack tip
%node, use the smallest value of r instead.
da=zeros(1,length(EdgeNodes));
db=zeros(1,length(EdgeNodes));
for i=1:length(EdgeNodes)
    da(i)=a-abs(EdgeNodes(1,i));
    db(i)=abs(EdgeNodes(2,i));
    dr(i)=sqrt((db(i)^2)+(da(i)^2));
end
%Check dr
%dr should hold two values of 0, two values greater than the crack length,
%and four values of each other node.
drSort=sort(dr,'descend');

%Assign heat generation values to radial points.
err=min(r); %Address the crack tip singularity.
qdot=zeros(length(t),length(dr));
for i=1:length(t)
    for j=1:length(dr)
       if dr(1,j)==0 %For the crack tip, select the value closest to the tip
           dr(1,j)=min(r);
           J(i,j)=KJ(i)/sqrt(2*pi*dr(1,j));
           qdot(i,j)=rho.*J(i,j).^2;
       else
           J(i,j)=KJ(i)/sqrt(2*pi*dr(1,j));
           qdot(i,j)=rho.*J(i,j).^2;
    end
    end
end
checknqd=sort(qdot(1,:),'descend');

thermalBC(thermalmodelT,'Edge',[1 2 7 8],'Temperature',Tamb);
thermalProperties(thermalmodelT,'ThermalConductivity',k,'MassDensity',B,'SpecificHeat',Cp);
thermalIC(thermalmodelT,Tamb);
internalHeatSource(thermalmodelT,1e9);
R=solve(thermalmodelT,t);
T=R.Temperature;

figure
pdeplot(thermalmodelT,'XYData',T(:,end),'Contour','on','ColorMap','hot');

[Christian McKinney]

Did you ever find the answer?