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');