function [u, conductorIDs] = solveDCConduction(self,coefstruct,constCoeffsBCs)
% Internal only.
% Identify conductors. If no conductors are found, results will be an
% all zero solution.

%   Copyright 2022 The MathWorks, Inc.

numSubdomains = numel(coefstruct.Coefficients.c);
conductorsFlag = cellfun(@(c) isa(c,'function_handle') || c ~= 0, coefstruct.Coefficients.c);
subdomains = 1:numSubdomains;
conductorIDs = subdomains(conductorsFlag);
domainRegion = 'Face';
if size(self.Mesh.Nodes,1) == 3 % ~isTwoD
    domainRegion = 'Cell';
end
conductorDOF = findNodes(self.Mesh,'region',domainRegion,conductorIDs);
u = zeros(size(self.Mesh.Nodes,2),1);

if isempty(conductorDOF)
    return;
end

[p,e,t] = self.Mesh.meshToPet();
u0 = pdeuxpd(p,0,self.PDESystemSize);
femodel = pde.DiscretizedPDEModel(self,p,e,t,coefstruct,u0,constCoeffsBCs);

% For linear conduction update femodel directly to retain  only
% conduction DoFs.
femodel.K = femodel.K(conductorDOF,conductorDOF);
femodel.A = femodel.A(conductorDOF,conductorDOF);
femodel.F = femodel.F(conductorDOF);
femodel.Q = femodel.Q(conductorDOF,conductorDOF);
femodel.G = femodel.G(conductorDOF);
femodel.H = femodel.H(:,conductorDOF);
femodel.R = femodel.R;
uConductor = self.solveStationary(femodel);
u(conductorDOF) = uConductor;
end