function [u,res]=pdenonlin(varargin) % b,p,e,t,c,a,f,varargin
%PDENONLIN Solve nonlinear PDE problem.
%
%       [U,RES]=PDENONLIN(B,P,E,T,C,A,F) solves the nonlinear
%       PDE problem -div(c*grad(u))+a*u=f, on a mesh described by P, E, and T,
%       with boundary conditions given by the function name B, using damped
%       Newton iteration with the Armijo-Goldstein line search strategy.
%
%       [U,RES]=PDENONLIN(PDEM,C,A,F...) Allows the boundary conditions B 
%       and mesh (P,E,T) to be defined by PDEM - a PDEModel object that 
%       contains this data.
%
%       The solution u is represented as the MATLAB column vector U.
%       See ASSEMPDE for details. RES is the 2-norm of the Newton step
%       residuals.
%
%       B describes the boundary conditions of the PDE problem. B
%       can either be a Boundary Condition Matrix or the name of Boundary
%       MATLAB-file. See PDEBOUND for details.
%
%       The optional arguments TOL and U0 signify a relative tolerance
%       parameter and an initial solution respectively.
%
%       The coefficients C, A, F of the PDE problem can be given in a
%       wide variety of ways. In the context of PDENONLIN the
%       coefficients can depend on U. The coefficients cannot depend
%       on T, the time.  See ASSEMPDE for details.
%
%       Fine tuning parameters may be sent to the solver in the form of
%       Parameter-Value pairs. It is possible to choose between computing
%       the full Jacobian, a "lumped" approximation or a fixed point
%       iteration and supply an initial guess of the solution.
%       The maximum number of iterations, the minimal damping factor of
%       the search vector, the size and norm of the residual at termination
%       can be controlled. These adjustments are made by setting one or more
%       of the following property/value pairs:
%
%       Property     Value/{Default}                Description
%       ----------------------------------------------------------------------------
%       Jacobian    'lumped'|{'full'}(3D)|      Jacobian approximation
%                            {'fixed'}(2D)
%       U0          string|numeric {0}          Set initial guess
%       Tol         numeric scalar {1e-4}       Acceptable residual
%       MaxIter     numeric scalar {25}         Maximum iterations
%       MinStep     numeric scalar {1/2^16}     Minimal damping factor
%       Report      'on'|{'off'}                Write convergence info.
%       Norm        numeric|{Inf}|'energy'      Residual norm
%
%       See also: ASSEMPDE, PDEBOUND, CREATEPDE, pde.PDEModel
%
%       Diagnostics: If the Newton-iteration does not converge, the
%       error messages 'Too many iterations' or 'Stepsize too small' are
%       displayed.
%
%       Note: This function does not support quadratic triangular elements
%

%       Copyright 1994-2020 The MathWorks, Inc.

% Error checks
if nargin > 0
    [varargin{:}] = convertStringsToChars(varargin{:});
end
nargs = nargin;
narginchk(4,21)

meshedpde=false; 
if isa(varargin{1}, 'pde.PDEModel') 
    mypde = varargin{1}; 
    if isa(mypde.Mesh, 'pde.FEMesh')
      meshedpde=true;
    end
end

offset=0;
if meshedpde
    mypde = varargin{1};    
    msh = mypde.Mesh;
    [p,e,t] = msh.meshToPet();
    b = pde.PDEModel(mypde.PDESystemSize);  
    b.BoundaryConditions = mypde.BoundaryConditions;
    offset=3;
    nargs = nargs+3;  
else
    b = varargin{1};  
    p = varargin{2};  
    e = varargin{3};   
    t = varargin{4};   
end
quadraticTriNotSupported(t);
 
if rem(nargs+1,2) || (nargs<7)
   error(message('pde:pdenonlin:nargin'))
end

c = varargin{5-offset};   
a = varargin{6-offset};   
f = varargin{7-offset};   


[dim, np] = size(p);
jacobianDefault='fixed';
if (dim==3)
    jacobianDefault='full';
end


% Recover Param/Value pairs from argument list
argsToPass = {};
if nargs > 7
argsToPass = varargin((8-offset):end);
end
parser = inputParser;
parser.addParameter('Jacobian', jacobianDefault,@(x) validateattributes(x,{'char','string'},{'nonempty','scalartext'}) );
parser.addParameter('Tol',1E-4,@(x) validateattributes(x,{'double'},{'nonsparse','real','finite','positive','scalar'}));
parser.addParameter('U0',0);
parser.addParameter('MaxIter',25,@(x) validateattributes(x,{'double'},{'nonsparse','real','finite','positive','scalar'}));
parser.addParameter('MinStep',2^(-16),@(x) validateattributes(x,{'double'},{'nonsparse','real','finite','positive','scalar'}));
parser.addParameter('Report', 'off',@(x) validateattributes(x,{'char','string'},{'nonempty','scalartext'}) );
parser.addParameter('Norm', 'inf',@(x) validateattributes(x,{'char','string','double'},{'nonempty'}) );
parser.parse(argsToPass{:});

jacobian = validatestring(parser.Results.Jacobian,{'lumped','full','fixed'});
tol = parser.Results.Tol;
u0 = parser.Results.U0;
maxiter = parser.Results.MaxIter;
MinStep = parser.Results.MinStep;
report=0;
reportOnOff = validatestring(parser.Results.Report,{'on','off'});
if ~isempty(reportOnOff) && startsWith(reportOnOff,'on')
    report = 1;
end
usenorm = parser.Results.Norm;
if ischar(usenorm)
    usenorm = validatestring(usenorm,{'inf','energy'});
end

if isnumeric(usenorm) && (~isscalar(usenorm) )
    error(message('pde:pdenonlin:InvalidNorm'));
end

if(dim==3 && ~strcmp(jacobian, 'full'))
  error(message('pde:pdenonlin:JacobianInvalidString3D'));
end

% Find an initial guess that satisfies the BC
thePde = getPde( p, e, t, b, f, u0, []);
N = thePde.PDESystemSize;
u = pdeuxpd(p,u0,N);
[K,M,F,Q,G,H,R]=assempde(b,p,e,t,c,a,f,u);
if ~any(any(H))
  H=sparse(0,size(H,2));
  R=sparse(0,1);
end

if(all(u==0))
  % if the initial guess is all zeros, use the stiffness
  % matrix as the initial Jacobian to improve the guess
  u = [K+M+Q H'; H zeros(size(H,1))] \ [F+G; R];
else
  u=[u;H'\((F+G)-(K+M+Q)*u)];
end

if ~all(isfinite(u),'all')
  error(message('pde:pdenonlin:InvalidInitGuess'))
end

global pdenonb pdenonp pdenone pdenont pdenonc pdenona pdenonf pdenonn
pdenonb=b;pdenonp=p;pdenone=e;pdenont=t;pdenonc=c;pdenona=a;pdenonf=f;
pdenonn=np*N;
global  pdenonK

r=pderesid(0,u);
K=pdenonK; % pdenonK may change after call to pderesid or numjac
if ischar(usenorm)
  if strcmp(usenorm,'energy')
    nr2=r'*K*r;
  else
    nr2=norm(r,usenorm)^2;
  end
else
  nr2=(norm(r,usenorm)/length(r)^(1/usenorm))^2;
end
res=sqrt(nr2);
if report
  disp(['Iteration     Residual     Step size  Jacobian: ',jacobian])
  fprintf('%4i%20.4e\n',0,res);
end


if strcmp(jacobian,'full') && (sqrt(nr2)> tol)
  %                                       c, a,            f,           u
  [a1,a2,a3,a4,a5,a6,a7]=assempde(b,p,e,t,0,ones(N*N,1), zeros(N,1), ones(np*N,1));
  if ~any(any(a6))
    a6=sparse(0,size(a6,2));
  end
  S=([a2 a6';a6 zeros(size(a6,1))]~=0);
  thresh = calcThresh(u, np*N);
  pdeResidual = str2func('pderesid');
  [dfdu0,fac,G]=numjac(pdeResidual,0,u,r,thresh,[],0,S,[]);
  first_time=1;
end

tmp=u;cor=zeros(size(u));iter=0;     % Gauss-Newton iteration
while sqrt(nr2)> tol
  iter=iter+1;
  if iter>maxiter,error(message('pde:pdenonlin:TooManyIter')),end
  J=K;
  if strcmp(jacobian,'lumped')
    % Find dc/du da/du df/du
    v=full(u(1:np*N));
    du=(max(abs(v))-min(abs(v)))/max(100,sqrt(np));
    if du==0, du=1e-5; end
    [dcc,daa,dff]=pdetxpd(p,t,du+v,c,a,f);
    [dccm,daam,dffm]=pdetxpd(p,t,v-du,c,a,f);
    dcc=(dcc-dccm)/2/du;
    daa=(daa-daam)/2/du;
    dff=(dff-dffm)/2/du;
    [J1,J3,a1]=thePde.assema(p,t,dcc,-dff,zeros(N,1));
    [a2,J2,a3]=thePde.assema(p,t,0,daa,zeros(N,1));
    J(1:np*N,1:np*N)=K(1:np*N,1:np*N)+spdiags((J1+J2)*v,0,np*N,np*N)+J3;
  elseif strcmp(jacobian,'full')
    if ~first_time
      thresh = calcThresh(u, np*N);
      [J,fac,G]=numjac(pdeResidual,0,u,r,thresh,fac,0,S,G);
    else
      J=dfdu0;
    end
    first_time=0;
  end
  alfa=1;
  lastwarn('','');
  warnState = warning('off', 'MATLAB:singularMatrix');
  cor=-J\r;
  [~,wid] = lastwarn;
  if(strcmp(wid,'MATLAB:singularMatrix'))
    warning(warnState);
    error(message('pde:pdenonlin:SingularJacobian'))
  end
  warning(warnState);
  
  while(true)
    if alfa<MinStep, error(message('pde:pdenonlin:SmallStepSize')), end
    tmp=u+alfa*cor;
    rr=pderesid(0,tmp);
    if ischar(usenorm)
      if strcmp(usenorm,'energy')
        nrr=rr'*K*rr;
      else
        nrr=norm(rr,usenorm)^2;
      end
    else
      nrr=(norm(rr,usenorm)/length(rr)^(1/usenorm))^2;
    end
    if nr2-nrr < alfa/2*nr2
      alfa=alfa/2;
    else
      u=tmp;
      r=rr;
      nr2=nrr;
      K=pdenonK;
      if report
        fprintf('%4i%20.4e%12.7f\n',iter,sqrt(nr2),alfa);
      end
      break
    end
  end
  res=[res;sqrt(nr2)];
end

u=full(u(1:pdenonn));
clear global pdenonb pdenonp pdenone pdenont pdenonc pdenona pdenonf
clear global  pdenonK pdenonF pdenonn

end

function thresh = calcThresh(u, numDofs)
nu = size(u,1);
thresh=zeros(nu,1);
threshFact = 1e-6;
if all(u(1:numDofs)==0)
  thresh(1:numDofs)=sqrt(eps)*ones(numDofs,1);
else
  thresh(1:numDofs)=threshFact*max(abs(u(1:numDofs)))*ones(numDofs,1);
end
if numDofs<nu
  if all(u(numDofs+1:nu)==0)
    thresh(numDofs+1:nu)=sqrt(eps)*ones(nu-numDofs,1);
  else
    thresh(numDofs+1:nu)=threshFact*max(abs(u(numDofs+1:nu)))*ones(nu-numDofs,1);
  end
end

end

