function output=linearize (self,inputOptions)
%LINEARIZE  Linearize thermal model.
%
%   LINEARIZE extracts a sparse linear model of a thermal model. 
%   For transient thermal models, LINEARIZE uses time t=0. Use 
%   LINEARIZEINPUT to specify the inputs of the linear model that 
%   correspond to external forcing, such as heat sources. Use 
%   LINEARIZEOUTPUT to specify the outputs of the linear model that
%   correspond to the selected regions of the geometry.
%
%   SYS = LINEARIZE(thermalmodel) returns the linearized model SYS 
%   in SPARSS form. This syntax requires Control System Toolbox. 
%   Use SYS.InputName and SYS.OutputGroup to locate the various 
%   inputs and output groups.
%
%   MX = LINEARIZE(thermalmodel,'OutputType','matrices') returns the
%   corresponding A,B,C,D,E or M,K,B,F matrices in a structure MX.
%
%   Example:
%      tmodel = createpde('thermal','transient');
%      geometryFromEdges(tmodel,@crackg);
%      generateMesh(tmodel);
%      thermalProperties(tmodel,'ThermalConductivity',1,...
%         'MassDensity',1,...
%         'SpecificHeat',1);
%      thermalBC(tmodel,'Edge',6,'Temperature',20,'Label','TempBC');
%      thermalBC(tmodel,'Edge',1,'HeatFlux',-10,'Label','FluxBC');
%      internalHeatSource(tmodel,25,'Label','HeatSource');
%      thermalIC(tmodel,0);

%      % Create inputs for HeatSource, TempBC and FluxBC on crackg
%      linearizeInput(tmodel,'HeatSource');
%      linearizeInput(tmodel,'TempBC');
%      linearizeInput(tmodel,'FluxBC');
%      % Measure temperature of Edge 2
%      linearizeOutput(tmodel,'Edge',2);
%      sys = linearize(tmodel);
%
%   See also pde.ThermalModel/linearizeInput, pde.ThermalModel/linearizeOutput

%    Copyright 2021 The MathWorks, Inc.

arguments
    self pde.ThermalModel
    inputOptions.OutputType {mustBeTextScalar}
end

state.time=0;
data = assembleFEMatrices(self,state);
K = data.K;
H = data.H;
R = data.R;
M = data.M;

% Get A and E
if isequal(nnz(H),0)
    A=-K;
    E=M;
else
    A = -[K H';H zeros(size(R,1))];
    E = blkdiag(M,zeros(size(R,1)));
end

%Inputs B: loads and BCs
[B,InputName,InputGroup]=getB(self,A,K,H);

% Outputs C: regions of interest
[C,OutputName,OutputGroup]=getC(self,A);

if isfield(inputOptions,'OutputType')
    if ~startsWith('matrices',inputOptions.OutputType,'IgnoreCase',true)
      error(message('pde:linearize:invalidVal',inputOptions.OutputType));
   end
   output.A=A;
   output.B=B;
   output.C=C;
   output.E=E;
else
   D = sparse(size(C,1),size(B,2));
   output = sparss(A,B,C,D,E,...
      'InputName',InputName, 'InputGroup', InputGroup,...
      'OutputName',OutputName, 'OutputGroup', OutputGroup);
end
end

function [bCstb,InputName,InputGroup]=getB(userModel,A,K,H)
nBCRows=size(K,1);
inputs=userModel.LinearizeInputs;
numInputs=numel(inputs);
InputName= cell(numInputs,1);
InputGroup = struct;

state.time=0;
bCstb = sparse(size(A,1),numInputs);

for i=1:numInputs
   try
      actBC=findThermalBC(userModel.BoundaryConditions,'Label',inputs(i).Label);
      if ~isempty(actBC.Temperature)
         [~,nodesH]=find(H);
         nodes=findNodes(userModel.Mesh,'region',inputs(i).RegionType,inputs(i).RegionID);
         idx=find(ismember(nodesH',nodes));
         bCstb(nBCRows+idx,:)=0;
         bCstb(nBCRows+idx,i)=1;

      elseif ~isempty(actBC.HeatFlux)
         model=createModelWithSingleLoad(userModel,actBC);
         DATA = assembleFEMatrices(model,'G',state);
         Bcol = DATA.G;
         bCstb(1:nBCRows,i)=Bcol;
      end
   catch
      actBC=findHeatSource(userModel.HeatSources,'Label',inputs(i).Label);
      if isprop(actBC,'HeatSource')
         model=createModelWithSingleLoad(userModel,actBC);
         DATA = assembleFEMatrices(model,'F',state);
         Bcol = DATA.F;
         bCstb(1:nBCRows,i)=Bcol;
      end
   end
   InputName{i} = inputs(i).Label;
end

end

function [C,OutputName,OutputGroup]=getC(userModel,A)
numOutputs=numel(userModel.LinearizeOutputs);
f = strings(numOutputs,1);
c = cell(numOutputs,1);
ny=0;
ng = 0;
OutputName = strings(0,1);
iy=[];
for i=1:numOutputs
    regionType=userModel.LinearizeOutputs(i).RegionType;
    regionID=userModel.LinearizeOutputs(i).RegionID;
    
    try
        nodes = findNodes(userModel.Mesh,'region',regionType,regionID);
    catch
        error(message('pde:linearize:outdatedMesh',regionID));
    end
    
    nn=numel(nodes);
    iy=[iy,nodes];
    yName = string(regionType) + regionID;
    
    if nn==1
        OutputName=[OutputName; yName];
    else
        ng = ng+1;
        f(ng) = yName;
        c{ng} = ny+1:ny+nn;
        OutputName=[OutputName; strings(nn,1)];
    end
    ny = ny+nn;
end
C = sparse(1:ny,iy,ones(1,ny),ny,size(A,1));
OutputGroup = cell2struct(c(1:ng,:),f(1:ng,:),1);
end

function pdem=createModelWithSingleLoad(userModel,bc)
modelType='thermal';

pdem = createpde(modelType,userModel.AnalysisType);
pdem.Geometry=userModel.Geometry;
pdem.Mesh=userModel.Mesh;
pdem.MaterialProperties=userModel.MaterialProperties;

if any(ismember(fieldnames(bc),'HeatSource'))
    internalHeatSource(pdem,1,bc.RegionType,bc.RegionID);
else
    thermalBC(pdem,bc.RegionType,bc.RegionID,'HeatFlux',1);
end
end

