function [newnode, newelem, newface] = meshrefine(node, elem, varargin)
%
% [newnode,newelem,newface]=meshrefine(node,elem,face,opt)
%
% refine a tetrahedral mesh by adding new nodes or constraints
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input parameters:
%      node: existing tetrahedral mesh node list
%      elem: existing tetrahedral element list
%      face: (optional) existing tetrahedral mesh surface triangle list
%      opt:  options for mesh refinement:
%        if opt is a Nx3 array, opt is treated as a list of new nodes to
%          be inserted into the mesh. the new nodes must be located on the
%          surface or inside the original mesh. external nodes are
%          discarded, unless the opt.extcmdopt is specified.
%        if opt is a vector with a length that equals to that of node,
%          it will be used to specify the desired edge-length at each node;
%          setting a node value to 0 will by-pass the refinement at this node
%        if opt is a vector with a length that equals to that of elem,
%          it will be used as the desired maximum element volume of each
%          tetrahedron; setting to 0 will by-pass the refinement of that element.
%        if opt is a struct, it can have the following fields:
%          opt.newnode: same as setting opt to an Nx3 array
%          opt.reratio: radius-edge ratio, by default, iso2mesh uses 1.414
%          opt.maxvol: maximum element volume
%          opt.sizefield: a vector specifying either the desired edge-length
%              at each node, or the maximum volume constraint within each
%              tetrahedron, see above for details.
%          opt.extcmdopt: by default, meshrefine can only insert nodes
%              that are inside the original mesh. if one prefers to insert
%              nodes that are outside of the original mesh, one can define
%              this parameter to specify the meshing option (for tetgen)
%              for the extended domain, i.e. the convex hull including
%              both the original and the external nodes. If not defined,
%              '-Y' option is used by default (prevent tetgen from
%              inserting new nodes on the surface).
%          opt.extlabel: when external nodes are inserted, the new elements
%              will be assigned with an element label to group them
%              together, by default, this label is 0, unless opt.extlabel
%              is given
%          opt.extcorelabel: when external nodes are inserted, par of the
%              new elements share the polyhedra between the inserted nodes,
%              these special elements will be marked by opt.extcorelabel,
%              otherwise the label will be set to -1
%
% outputs:
%      newnode: node coordinates of the tetrahedral mesh
%      newelem: element list of the tetrahedral mesh
%      newface: mesh surface element list of the tetrahedral mesh
%             the last column denotes the boundary ID
%
% examples:
%
%     [node,face,elem]=meshasphere([0 0 0],24,1,2);
%     elem(:,5)=1;
%
%     % inserting nodes that are inside the original mesh
%     innernodes=double([1 1 1; 2 2 2; 3 3 3]);
%     [newno,newel]=meshrefine(node,elem,innernodes);
%     all(ismember(round(innernodes*1e10)*1e-10,round(newno*1e10)*1e-10,'rows'))
%     plotmesh(newno,[],newel,'x>-3')
%
%     % inserting nodes that are external to the original mesh
%     extnodes=double([-5 -5 25;-5 5 25;5 5 25;5 -5 25]);
%     [newno,newel]=meshrefine(node,elem,struct('newnode',extnodes,'extcmdopt','-Y'));
%     figure;
%     plotmesh(newno,[],newel,'x>-3')
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

exesuff = getexeext;
exesuff = fallbackexeext(exesuff, 'tetgen');

newpt = [];
sizefield = [];
if (size(node, 2) == 4)
    sizefield = node(:, 4);
    node = node(:, 1:3);
end
opt = struct;
if (length(varargin) == 1)
    face = [];
    if (isstruct(varargin{1}))
        opt = varargin{1};
    else
        if (length(varargin{1}) == size(node, 1) || length(varargin{1}) == size(elem, 1))
            sizefield = varargin{1};
        else
            newpt = varargin{1};
        end
    end
elseif (length(varargin) >= 2)
    face = varargin{1};
    if (isstruct(varargin{2}))
        opt = varargin{2};
    else
        if (length(varargin{2}) == size(node, 1) || length(varargin{1}) == size(elem, 1))
            sizefield = varargin{2};
        else
            newpt = varargin{2};
        end
    end
else
    error('meshrefine requires at least 3 inputs');
end
if (isstruct(opt) && isfield(opt, 'newnode'))
    newpt = opt.newnode;
end
if (isstruct(opt) && isfield(opt, 'sizefield'))
    sizefield = opt.sizefield;
end

% call tetgen to create volumetric mesh
deletemeshfile(mwpath('pre_refine.*'));
deletemeshfile(mwpath('post_refine.*'));

moreopt = '';
setquality = 0;
if (isstruct(opt) && isfield(opt, 'reratio'))
    moreopt = [moreopt sprintf(' -q %.10f ', opt.reratio)];
    setquality = 1;
end
if (isstruct(opt) && isfield(opt, 'maxvol'))
    moreopt = [moreopt sprintf(' -a %.10f ', opt.maxvol)];
end

externalpt = [];
if (isstruct(opt) && isfield(opt, 'extcmdopt'))
    isinside = tsearchn(node(:, 1:3), elem(:, 1:4), newpt(:, 1:3));
    externalpt = newpt(isnan(isinside), :);
    newpt = newpt(find(isnan(isinside) == 0), :);
end

if (~isempty(newpt))
    if (size(newpt, 1) < 4)
        newpt = [newpt; repmat(newpt(1, :), 4 - size(newpt, 1), 1)];
    end
    savetetgennode(newpt, mwpath('pre_refine.1.a.node'));
    moreopt = ' -i ';
end

if (~isempty(sizefield))
    if (length(sizefield) == size(node, 1))
        fid = fopen(mwpath('pre_refine.1.mtr'), 'wt');
        fprintf(fid, '%d 1\n', size(sizefield, 1));
        fprintf(fid, '%.16f\n', sizefield);
        fclose(fid);
        moreopt = [moreopt ' -qa '];
    else
        fid = fopen(mwpath('pre_refine.1.vol'), 'wt');
        fprintf(fid, '%d\n', size(sizefield, 1));
        fprintf(fid, '%d\t%.16f\n', [(1:size(sizefield, 1))' sizefield]');
        fclose(fid);
        moreopt = [moreopt ' -qa '];
    end
end

if (size(elem, 2) == 3 && setquality == 0)
    if (~isempty(newpt))
        error('inserting new point can not be used for surfaces');
    end
    nedge = savegts(node, elem, mwpath('pre_refine.gts'));
    exesuff = fallbackexeext(getexeext, 'gtsrefine');
elseif (size(elem, 2) == 3)
    savesurfpoly(node, elem, [], [], [], [], mwpath('pre_refine.poly'));
else
    savetetgennode(node, mwpath('pre_refine.1.node'));
    savetetgenele (elem, mwpath('pre_refine.1.ele'));
end

fprintf(1, 'refining the input mesh ...\n');

if (size(elem, 2) == 3 && setquality == 0)
    if (isstruct(opt) && isfield(opt, 'scale'))
        moreopt = sprintf('%s -n %d ', moreopt, round(nedge * opt.scale));
    else
        error('you must give opt.scale value for refining a surface');
    end
end
if (isstruct(opt) && isfield(opt, 'moreopt'))
    moreopt = [moreopt opt.moreopt];
end

if (size(elem, 2) == 3 && setquality == 0)
    status = system([' "' mcpath('gtsrefine') exesuff '" ' moreopt ' < "' ...
                     mwpath('pre_refine.gts') '" > "' mwpath('post_refine.gts') '"']);
    if (status)
        error('gtsrefine command failed');
    end
    [newnode, newelem] = readgts(mwpath('post_refine.gts'));
    newface = newelem;
elseif (size(elem, 2) == 3)
    status = system([' "' mcpath('tetgen') exesuff '" ' moreopt ' -p -A "' mwpath('pre_refine.poly') '"']);
    if (status)
        error('tetgen command failed');
    end
    [newnode, newelem, newface] = readtetgen(mwpath('pre_refine.1'));
elseif (~isempty(moreopt))
    status = system([' "' mcpath('tetgen') exesuff '" ' moreopt ' -r "' mwpath('pre_refine.1') '"']);
    if (status)
        error('tetgen command failed');
    end
    [newnode, newelem, newface] = readtetgen(mwpath('pre_refine.2'));
else
    newnode = node;
    newelem = elem;
    newface = face;
end

if (~isempty(externalpt)) % user request to insert nodes that are outside of the original mesh
    % create a mesh including the external points
    externalpt = unique(externalpt, 'rows');
    allnode = [newnode; externalpt];

    % define the convexhull as the external surface
    try
        outface = convhull(allnode, 'simplify', false);
    catch
        if (isoctavemesh)
            outface = volface(delaunayn(allnode));
        else
            outface = convhulln(allnode);
        end
    end
    outface = sort(outface, 2);
    face = volface(newelem(:, 1:4));
    inface = sort(face(:, 1:3), 2);

    % define the surface that bounds the newly extended convex hull space
    bothsides = removedupelem([outface; inface]);

    % define a seed point to avoid meshing the interior space
    holelist = surfseeds(newnode, face(:, 1:3));

    % mesh the extended space
    ISO2MESH_TETGENOPT = jsonopt('extcmdopt', '-Y', opt);
    try
        if (size(bothsides, 1) >= size(inface, 1))
            [no, el] = surf2mesh(allnode, bothsides, [], [], 1, 10, [], holelist);
        else
            [no, el] = surf2mesh(allnode, bothsides, [], [], 1, 10);
        end
    catch
        bothsides = maxsurf(finddisconnsurf(bothsides), allnode);
        if (size(bothsides, 1) >= size(inface, 1))
            [no, el] = surf2mesh(allnode, bothsides, [], [], 1, 10, [], holelist);
        else
            [no, el] = surf2mesh(allnode, bothsides, [], [], 1, 10);
        end
    end
    [isinside, map] = ismember(round(no * 1e10) * 1e-10, round(allnode * 1e10) * 1e-10, 'rows');
    snid = [length(newnode) + 1:length(allnode)];

    if (~isempty(map == 0))
        oldsize = size(allnode, 1);
        allnode = [allnode; no(map == 0, :)];
        map(map == 0) = oldsize + 1:size(allnode, 1);
    end
    % merge the external space with the original mesh
    el2 = map(el(:, 1:4));

    % label all new elements with -1
    if (size(newelem, 2) == 5)
        el2(:, 5) = jsonopt('extlabel', 0, opt);
        % search elements that contain source(s) and save their id(s)
        iselm = ismember(el2(:, 1:4), snid);
        el2(sum(iselm, 2) >= 3, 5) = jsonopt('extcorelabel', -1, opt);
    end

    % merge nodes/elements and replace the original ones
    newnode = allnode;
    newelem = [newelem; el2];
end

% read in the generated mesh

fprintf(1, 'mesh refinement is complete\n');