function u = solveStationary(self,femodel)
% solveStationary - Internal function for solving a time-independent PDE.
%

% Copyright 2015-2023 The MathWorks, Inc.

K = femodel.K;
A = femodel.A;
F = femodel.F;
Q = femodel.Q;
G = femodel.G;
H = femodel.H;
R = femodel.R;

[null,orth]=pdenullorth(H);
if size(orth,2)==0
    ud=zeros(size(K,2),1);
else
    ud=full(orth*((H*orth)\R));
end

KK=K+A+Q;
FF=null'*((F+G)-KK*ud);
KK=null'*KK*null;
KK = self.checkForSymmetry(KK);

% Check rigid motion in case of structural static analysis.
if isprop(self.SolverOptions,'RigidMotionCheck') && self.SolverOptions.RigidMotionCheck && ~isempty(null)
    Dim = size(self.Mesh.Nodes,1);
    %Construct unit translation matrix.
    unitDisp = ones(size(self.Mesh.Nodes,2),1);
    utest = kron(eye(Dim),unitDisp);

    %Use max diag element to scale tolerance. This is cheaper than
    %computing norm. The factor 100 is to account for the fact that max
    %diag element will be lower than norm.
    maxDiagK = max(diag(KK))*100;
    
    %Scaling factor for tol to account for norm of translation matrix.
    utestNorm = sqrt(size(self.Mesh.Nodes,2)*3/3);
    
    tol = eps(maxDiagK*utestNorm);
    
    reactionQR = qr(KK*null'*utest);
    reactionR = triu(reactionQR);
    
    %Singular values of R.
    d = svd(reactionR(1:Dim,:));
    if any(d < tol)
        error(message('pde:StructuralModel:possibleRigidbodyMode'))
    end
end

if size(null,2)==0
    u=ud;
else
    u=null*(KK\FF)+ud;
end

end
