function p = jigglemesh(p,e,t,NVArgs)
%JIGGLEMESH Jiggle internal points of a triangular mesh.
%
%       P1=JIGGLEMESH(P,E,T) jiggles the triangular mesh by adjusting
%       the node point positions. The quality of the mesh normally
%       increases.
%
%       See PDETRIQ for more info on triangle quality.
%
%       The following Property/Value pairs are allowed:
%
%       Property    Value/{Default}           Description
%       ------------------------------------------------------------
%       Opt         {off}|mean|minimum        Optimization method
%       Iter        1 or 20 (see below)       Maximum number of iterations
%
%       Each mesh point that is not located on a mesh boundary is
%       moved toward the center of mass of the polygon formed by the
%       adjacent triangles. This process is repeated according to the
%       setting of the Opt and Iter variables:
%       - When Opt is set to off this process is repeated Iter times
%       (default: once).
%       - When Opt is set to mean the process is repeated until the
%       mean triangle quality does not significantly increase, or until
%       the bound Iter is reached (default: 20).
%       - When Opt is set to min the process is repeated until the
%       minimum triangle quality does not significantly increase, or until
%       the bound Iter is reached (default: 20).
%
%       See also INITMESH, PDETRIQ

%       Copyright 1994-2023 The MathWorks, Inc.

arguments
    p
    e
    t
    NVArgs.Opt {mustBeNonempty, mustBeMember(NVArgs.Opt,["off","mean","minimum"])} = "off"
    NVArgs.Iter (1,1) {mustBeNonnegative}
end

opt = NVArgs.Opt;
if ~isfield(NVArgs,"Iter") && strcmp(opt,"off")
    iter = 1;
elseif ~isfield(NVArgs,"Iter")
    iter = 20;
elseif isinf(NVArgs.Iter) && strcmp(opt,"off")
    error(message('pdelib:jigglemesh:OptOffIterInf'))
else
    iter = NVArgs.Iter;
end

% Determine interior non boundary points
j    = unique(e(1:2,:))';   % all points on subdomain boundaries
np   = size(p, 2);
i    = true(1, np);
i(j) = false;               % logical vector indicating if point is interior

if ~strcmp(opt,"off")
  q = pdetriq(p,t);
  if strcmp(opt,"minimum")
    q = min(q);
  else
    q = mean(q);
  end
end

for k = 1:iter
  X = sparse(t([1 2 3],:),t([2 3 1],:),p(1,t(1:3,:)),np,np);
  Y = sparse(t([1 2 3],:),t([2 3 1],:),p(2,t(1:3,:)),np,np);
  N = sparse(t([1 2 3],:),t([2 3 1],:),1,np,np);
  m = sum(N);
  X = sum(X)./m;
  Y = sum(Y)./m;
  p1 = p;
  p(1,i) = X(i);
  p(2,i) = Y(i);
  if ~strcmp(opt,"off")
    q1 = q;
    q = pdetriq(p,t);
    if strcmp(opt,"minimum")
      q = min(q);
    elseif strcmp(opt,"mean")
      q = mean(q);
    end
    if q < q1 % if triangle quality has decreased, revert and exit
      p = p1;
      break
    elseif q < q1+1e-4 % if triangle quality only slightly increased, exit
      break
    end
  end
end
end