function thermalLoad = generateThermalBoundaryLoad(PDESystemSize,coefstruct,subdomainIDforBoundary,Tref)
% generateThermalBoundaryLoad - Internal function to compute thermal
% effect on the boundary term.

% Copyright 2018-2023 The MathWorks, Inc.

if  any(isnan(subdomainIDforBoundary))  % a boundary of a single domain or a non-interface boundary of a multidomain model.
    if isscalar(coefstruct.Coefficients.CTE)
        nonNaNIdx =1; % Material packing found only one material even if the geomety is multidomain.
    else
        nonNaNIdx = subdomainIDforBoundary(~isnan(subdomainIDforBoundary));
    end
    CTE = coefstruct.Coefficients.CTE{nonNaNIdx};
    E = coefstruct.Coefficients.YoungsModulus{nonNaNIdx};
    nu =  coefstruct.Coefficients.PoissonsRatio{nonNaNIdx};
    thermalExpConstant = CTE*E/(1-2*nu);
    
    % Convention is such that if the first index is NaN then the edge is
    % oriented such that the normal is pointing inward, so flip the
    % sign to make it an outward normal. This issue is specific to
    % geometry created using decsg.
    if isnan(subdomainIDforBoundary(1))
        thermalExpConstant = -thermalExpConstant;
    end
else %multidomain, interface boundary.
    if isscalar(coefstruct.Coefficients.CTE)
        thermalExpConstant = 0; % interface is separating same material on either side, resultant thermal effect at this interface is zero.
    else
        primarySubdomain =  subdomainIDforBoundary(1);
        secondarySubdomain = subdomainIDforBoundary(2);
        
        CTE1 = coefstruct.Coefficients.CTE{primarySubdomain};
        E1 = coefstruct.Coefficients.YoungsModulus{primarySubdomain};
        nu1 =  coefstruct.Coefficients.PoissonsRatio{primarySubdomain};
        CTE2 = coefstruct.Coefficients.CTE{secondarySubdomain};
        E2 = coefstruct.Coefficients.YoungsModulus{secondarySubdomain};
        nu2 =  coefstruct.Coefficients.PoissonsRatio{secondarySubdomain};
        
        thermalExpConstant1 = CTE1*E1/(1-2*nu1); % back/left subdomain, normal is outward for this subdomain.
        thermalExpConstant2 = CTE2*E2/(1-2*nu2);
        
        thermalExpConstant = thermalExpConstant1 - thermalExpConstant2;
    end
end

T = coefstruct.Coefficients.Tintrp{1};
TimeStep = coefstruct.Coefficients.TimeStep{1};

    function g = gThermal(region,~)
        
        if (PDESystemSize == 3)
            qPoints = [region.x;region.y;region.z];
            normals = [region.nx;region.ny;region.nz];
        else
            qPoints = [region.x; region.y];
            normals = [region.nx;region.ny];            
        end
        
        if isnumeric(T)
            g = thermalExpConstant.*ones(PDESystemSize,1).*(T-Tref).*normals;
        else
            if isa(T,'pde.TransientThermalResults')
                if isempty(TimeStep)
                    TimeStep = numel(T.SolutionTimes);
                end
                Tgp = T.interpolateTemperature(qPoints,TimeStep);
            else
                Tgp = T.interpolateTemperature(qPoints);
            end
            
            if any(isnan(Tgp))
                nanIdx = isnan(Tgp);
                coordinatesOfNaNNodes = qPoints(:,nanIdx);
                nearestNodesInSourceMesh = T.Mesh.findNodes('nearest',coordinatesOfNaNNodes);
                Tgp(nanIdx) = T.Temperature(nearestNodesInSourceMesh);
            end
            
            g = thermalExpConstant.*ones(PDESystemSize,1).*(Tgp'-Tref).*normals;
        end
    end

thermalLoad = @gThermal;
end
