function Tintrp = interpolateThermalBodyLoad(self,bodyLoadStr)
% interpolateThermalBodyLoad - Internal function to interpolate thermal
% results on to the analysis mesh.

% Copyright 2018-2022 The MathWorks, Inc.

T = bodyLoadStr.Coefficients.Temperature{1};
TimeStep = bodyLoadStr.Coefficients.TimeStep{1};

if isa(T,'pde.TransientThermalResults')
    if isempty(TimeStep)
        TimeStep = numel(T.SolutionTimes);
    end
end

inputLoadMesh = T.Mesh;
analysisMesh = self.Mesh;

if T.IsTwoD
    inputDomainSize = area(inputLoadMesh);
    analysisDomainSize = area(analysisMesh);
    
else
    inputDomainSize = volume(inputLoadMesh);
    analysisDomainSize = volume(analysisMesh);
end

domainSizeDifference = abs(inputDomainSize - analysisDomainSize)/analysisDomainSize;
if domainSizeDifference > 0.05
    error(message('pde:bodyLoadSpecification:geometricSizeMismatch'));
end

if isa(T,'pde.TransientThermalResults')
    intrpTemperature = T.interpolateTemperature(analysisMesh.Nodes,TimeStep);
else
    intrpTemperature = T.interpolateTemperature(analysisMesh.Nodes);
end

nanIdx = isnan(intrpTemperature);

if any(nanIdx) % Replace NaN with nearest nodal values.
    coordinatesOfNaNNodes = analysisMesh.Nodes(:,nanIdx);
    nearestNodesInSourceMesh = inputLoadMesh.findNodes('nearest',coordinatesOfNaNNodes);
    intrpTemperature(nanIdx) = T.Temperature(nearestNodesInSourceMesh);
end

% Create a results object
tempModel = createpde('thermal');
tempModel.Mesh = analysisMesh;
geometryFromMesh(tempModel,analysisMesh.Nodes,analysisMesh.Elements);
tempModel.thermalProperties('ThermalConductivity',1);
Tintrp = pde.SteadyStateThermalResults(tempModel,intrpTemperature);

end

