function [u,p,e,t] = adaptmesh(varargin)
%ADAPTMESH Adaptive mesh generation and PDE solution.
%
%       [U,P,E,T]=ADAPTMESH(G,B,C,A,F,P1,V1,...) performs adaptive mesh
%       generation and PDE solution.  The large number of possible input
%       option is handled using property-value pair arguments.  The first
%       five arguments G, B, C, A, and F are not optional.
%
%       The function produces a solution U to the elliptic scalar PDE problem
%       -div(c*grad(u))+a*u=f where the problem geometry and boundary
%       conditions given by G and B.
%
%       The solution u is represented as the MATLAB column vector U.
%       See ASSEMPDE for details.
%
%       G describes the geometry of the PDE problem. G can
%       either be a Decomposed Geometry Matrix or the name of Geometry
%       MATLAB-file. See either DECSG or PDEGEOM for details.
%
%       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 adapted triangular mesh of the PDE problem is given by the triangle
%       data P, E, and T.  See either INITMESH or PDEGEOM for details.
%
%       The coefficients C, A and F of the PDE problem can
%       be given in a wide variety of ways.  See ASSEMPDE for details.
%
%       Valid property/value pairs include
%
%       Prop.   Value/{Default}         Description
%       ----------------------------------------------------------------------
%       Maxt    Positive integer {Inf}  Maximum number of new triangles
%       Ngen    Positive integer {10}   Maximum number of triangle generations
%       Mesh    P1, E1, T1              Initial mesh
%       Tripick {pdeadworst}|pdeadgsc   Triangle selection method
%       Par     Numeric {0.5}           Function parameter
%       Rmethod {longest}|regular       Triangle refinement method
%       Nonlin  on | off                Use nonlinear solver
%       Toln    numeric {1e-3}          Nonlinear tolerance
%       Init    string|numeric          Nonlinear initial solution value
%       Jac     {fixed}|lumped|full     Nonlinear solver Jacobian calculation
%       Norm    Numeric {Inf}           Nonlinear solver residual norm
%       MesherVersion   {preR2013a}|R2013a      Meshing algorithm used
%
%       Par is passed to the tripick function. Normally it is used as
%       tolerance of how well the solution fits the equation. No more than
%       Ngen successive refinements are attempted. Refinement is also
%       stopped when the number of triangles in the mesh exceeds the
%       Maxt.
%
%       P1, E1, and T1 are the input mesh data. This triangle mesh is used
%       as a starting mesh for the adaptive algorithm. If no initial mesh
%       is provided, the result of a call to INITMESH with no options is
%       used as initial mesh.
%
%       The triangle pick method is a user-definable triangle selection
%       method.  Given the error estimate computed by the function PDEJMPS,
%       the triangle pick method selects the triangles to be refined
%       in the next triangle generation. The function is called using the
%       arguments P, T, CC, AA, FF, U, ERRF, and PAR.  P and T represent the
%       current generation of triangles, CC, AA, FF are the current
%       coefficients for the PDE problem, expanded to triangle midpoints,
%       U is the current solution, ERRF is the computed error estimate,
%       and PAR, the function parameter, given to ADAPTMESH as optional
%       argument. The matrices CC, AA, FF, and ERRF all have NT columns,
%       where NT is the current number of triangles.
%       The number of rows in CC, AA, and FF are exactly the same as the
%       input arguments C, A, and F. ERRF has one row for each equation
%       in the system. There are two standard triangle selection methods
%       in the PDE Toolbox - PDEADWORST and PDEADGSC.
%       PDEADWORST selects triangles where ERRF exceeds a fraction
%       (default: 0.5) of the worst value. PDEADGSC selects triangles
%       using a relative tolerance criterion.
%
%       The refinement method is either 'longest' or 'regular'.
%       See REFINEMESH for details.
%
%       Also nonlinear PDE problems can be solved by the adaptive algorithms.
%       For nonlinear PDE problems, the 'Nonlin' parameter must be set
%       to 'on'. The nonlinear tolerance Toln and nonlinear initial value
%       U0 are passed to the nonlinear solver. See PDENONLIN for details.
%
%       The MesherVersion parameter selects the specific method used to
%       produce a mesh. Option R2013a selects a Delaunay-based algorithm that is
%       faster and more robust at dealing with non-standard geometries. The
%       default option, preR2013a, selects the method used in previous versions
%       of the toolbox.
%
%       See also ASSEMPDE, PDEBOUND, PDEGEOM, INITMESH, REFINEMESH, PDENONLIN

%       Copyright 1994-2023 The MathWorks, Inc.

narginchk(5,Inf)
[g,b,c,a,f] = varargin{1:5}; % deal out first five arguments
varargin = varargin(6:end);
if rem(length(varargin),2)
    error(message('pde:adaptmesh:KeyWithoutValue'))
end
meshInd = find(cellfun(@(x) any(strcmpi(x,'Mesh')),varargin));
for i = meshInd(end:-1:1)
    if isa(varargin{i+1},'cell') && isequal(size(varargin{i+1}), [1 3]) ...
            && all(cellfun(@isnumeric,varargin{i+1})) % next arg is {p,e,t}
        continue
    elseif i+3>length(varargin) % fewer than three args after 'Mesh'
        error(message('pde:adaptmesh:InvalidMeshValue'))
    elseif all(cellfun(@isnumeric,varargin(i+1:i+3))) % next three args are numeric
        varargin{i+1} = varargin(i+1:i+3);
        varargin(i+2:i+3) = [];
    else
        error(message('pde:adaptmesh:InvalidMeshValue'))
    end
end
try
    [u,p,e,t] = iadaptmesh(g,b,c,a,f,varargin{:});
catch err
    throw(err)
end
end

function [u,p,e,t] = iadaptmesh(g,b,c,a,f,NVArgs)
arguments
    g {mustBeA(g, ["double", "char", "string", "function_handle", "pde.AnalyticGeometry"])}
    b {mustBeA(b, ["double", "char", "string", "function_handle", "pde.PDEModel"])}
    c {mustBeA(c, ["double", "char", "string", "function_handle"])}
    a {mustBeA(a, ["double", "char", "string", "function_handle"])}
    f {mustBeA(f, ["double", "char", "string", "function_handle"])}
    NVArgs.Maxt (1,1) {mustBePositive} = Inf
    NVArgs.Ngen (1,1) {mustBePositive} = 10
    NVArgs.Mesh (1,3) cell
    NVArgs.Tripick {mustBeNonempty, mustBeTextScalar} = "pdeadworst"
    NVArgs.Par (1,1) {mustBeInRange(NVArgs.Par,0,1)} = 0.5
    NVArgs.Rmethod {mustBeNonempty, mustBeMember(NVArgs.Rmethod, ["regular", "longest"])} = "longest"
    NVArgs.Nonlin {mustBeNonempty, mustBeMember(NVArgs.Nonlin, ["off", "on"])} = "off"
    NVArgs.Toln (1,1) {mustBePositive} = 1e-4
    NVArgs.Init {mustBeA(NVArgs.Init, ["double", "char", "string"])} = 0
    NVArgs.Jac {mustBeNonempty, mustBeMember(NVArgs.Jac, ["fixed", "lumped", "full"])} = "fixed"
    NVArgs.Norm (1,1) {mustBeNumeric, mustBeNonNan, mustBeReal, mustBePositiveOrNegInf} = Inf
    NVArgs.MesherVersion {mustBeNonempty, mustBeMember(NVArgs.MesherVersion, ["preR2013a", "R2013a"])} = "preR2013a"
end

if isa(b, 'pde.PDEModel') && isa(b.Mesh, 'pde.FEMesh')
    error(message('pde:adaptmesh:noMeshedPDEModelSupport'));
end

Maxt = NVArgs.Maxt;
Ngen = min(NVArgs.Ngen,intmax);
if isfield(NVArgs,'Mesh')
    [p,e,t] = NVArgs.Mesh{:};
else
    [p,e,t] = initmesh(g, 'MesherVersion', NVArgs.MesherVersion);
end
Tripick = NVArgs.Tripick;
Par     = NVArgs.Par;
Rmethod = NVArgs.Rmethod;
nonl    = strcmpi(NVArgs.Nonlin, "on");
Toln    = NVArgs.Toln;
u       = NVArgs.Init;
Jac     = NVArgs.Jac;
norm    = NVArgs.Norm;

alfa    = 0.15;
beta    = 0.15;
mexp    = 1;

for gen = 1:Ngen+1
    fprintf('Number of triangles: %g\n',size(t,2))
    if nonl
        u = pdenonlin(b,p,e,t,c,a,f,'jacobian',Jac,'U0',u,'tol',Toln,'norm',norm);
    else
        u = assempde(b,p,e,t,c,a,f);
    end

    if any(isnan(u))
        error(message('pde:adaptmesh:NaNinSolution'))
    end

    % Expand values
    [cc,aa,ff] = pdetxpd(p,t,u,c,a,f);

    errf = pdejmps(p,t,cc,aa,ff,u,alfa,beta,mexp);

    i = feval(Tripick,p,t,cc,aa,ff,u,errf,Par);

    if isempty(i)
        fprintf('\nAdaption completed.\n')
        break
    elseif size(t,2)>Maxt
        fprintf('\nMaximum number of triangles obtained.\n');
        break
    elseif gen==Ngen+1
        fprintf('\nMaximum number of refinement passes obtained.\n');
        break
    end

    tl=i';
    u = reshape(u,size(p,2),[]);
    [p,e,t,u] = refinemesh(g,p,e,t,u,tl,Rmethod);
    u = u(:);
end
end % END iadaptmesh

function mustBePositiveOrNegInf(Norm)
if Norm~=-Inf
    mustBePositive(Norm)
end
end % END mustBePositiveOrNegInf