surfboolean

 [newnode,newelem,newelem0]=surfboolean(node1,elem1,op2,node2,elem2,op3,node3,elem3,...)

 merge two or more triangular meshes and resolve intersecting elements
 
 author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>

 input:
      node: node coordinates, dimension (nn,3)
      elem: triangle surfaces (ne,3)
      op:  a string of a boolean operator, possible op values include
           'union' or 'or': the outter surface of the union of the enclosed space
           'inter' or 'and': the surface of the domain contained by both meshes
           'diff' or '-': the surface of the domain in mesh 1 excluding that of
                   mesh 2
           'all' or 'xor' or '+': the output contains 4 subsurfaces, identified by the 4th
                  column of newelem:
                    1: mesh 1 outside of mesh 2
                    2: mesh 2 outside of mesh 1
                    3: mesh 1 inside of mesh 2
                    4: mesh 2 inside of mesh 1
                  you can use newelem(find(mod(newelem(:,4),2)==1),:) to
                  get mesh 1 cut by mesh 2, or newelem(find(mod(newelem(:,4),2)==0),:) 
                  to get mesh 2 cut by mesh 1;
           'first': combine 1 and 3 from the output of 'all'
           'second': combine 2 and 4 from the output of 'all'
           'self': test for self-intersections; only the first mesh is
                   tested; other inputs are ignored.
           'decouple': separate two shells and make sure there is no intersection;
                   the input surfaces must be closed and ordered from outer to inner

 output:
      newnode: the node coordinates after boolean operations, dimension (nn,3)
      newelem: tetrahedral element or surfaces after boolean operations (nn,4) or (nhn,5)
      newelem0: when the operator is 'self', return the intersecting
               element list in terms of the input node list (experimental)

 example:

   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
   [node2,face2,elem2]=meshabox([0 0 0]+5,[10 10 10]+5,1,1);
   [newnode,newface]=surfboolean(node1,face1,'union',node2,face2);
   plotmesh(newnode,newface);
   figure;
   [newnode,newface]=surfboolean(node1,face1,'diff',node2,face2);
   plotmesh(newnode,newface,'x>5');

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function [newnode,newelem,newelem0]=surfboolean(node,elem,varargin)
0002 %
0003 % [newnode,newelem,newelem0]=surfboolean(node1,elem1,op2,node2,elem2,op3,node3,elem3,...)
0004 %
0005 % merge two or more triangular meshes and resolve intersecting elements
0006 %
0007 % author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
0008 %
0009 % input:
0010 %      node: node coordinates, dimension (nn,3)
0011 %      elem: triangle surfaces (ne,3)
0012 %      op:  a string of a boolean operator, possible op values include
0013 %           'union' or 'or': the outter surface of the union of the enclosed space
0014 %           'inter' or 'and': the surface of the domain contained by both meshes
0015 %           'diff' or '-': the surface of the domain in mesh 1 excluding that of
0016 %                   mesh 2
0017 %           'all' or 'xor' or '+': the output contains 4 subsurfaces, identified by the 4th
0018 %                  column of newelem:
0019 %                    1: mesh 1 outside of mesh 2
0020 %                    2: mesh 2 outside of mesh 1
0021 %                    3: mesh 1 inside of mesh 2
0022 %                    4: mesh 2 inside of mesh 1
0023 %                  you can use newelem(find(mod(newelem(:,4),2)==1),:) to
0024 %                  get mesh 1 cut by mesh 2, or newelem(find(mod(newelem(:,4),2)==0),:)
0025 %                  to get mesh 2 cut by mesh 1;
0026 %           'first': combine 1 and 3 from the output of 'all'
0027 %           'second': combine 2 and 4 from the output of 'all'
0028 %           'self': test for self-intersections; only the first mesh is
0029 %                   tested; other inputs are ignored.
0030 %           'decouple': separate two shells and make sure there is no intersection;
0031 %                   the input surfaces must be closed and ordered from outer to inner
0032 %
0033 % output:
0034 %      newnode: the node coordinates after boolean operations, dimension (nn,3)
0035 %      newelem: tetrahedral element or surfaces after boolean operations (nn,4) or (nhn,5)
0036 %      newelem0: when the operator is 'self', return the intersecting
0037 %               element list in terms of the input node list (experimental)
0038 %
0039 % example:
0040 %
0041 %   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
0042 %   [node2,face2,elem2]=meshabox([0 0 0]+5,[10 10 10]+5,1,1);
0043 %   [newnode,newface]=surfboolean(node1,face1,'union',node2,face2);
0044 %   plotmesh(newnode,newface);
0045 %   figure;
0046 %   [newnode,newface]=surfboolean(node1,face1,'diff',node2,face2);
0047 %   plotmesh(newnode,newface,'x>5');
0048 %
0049 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0050 %
0051 
0052 allinputs=varargin;
0053 opt=struct;
0054 if(length(allinputs)>0 && isstruct(allinputs{end}))
0055     opt=allinputs{end};
0056     allinputs{end}=[];
0057 end
0058 len=length(varargin);
0059 newnode=node;
0060 newelem=elem;
0061 if(len>0 && mod(len,3)~=0)
0062    error('you must give operator, node and element in a triplet form');
0063 end
0064 
0065 try
0066     exename=evalin('caller','ISO2MESH_SURFBOOLEAN');
0067 catch
0068     try
0069         exename=evalin('base','ISO2MESH_SURFBOOLEAN');
0070     catch
0071         exename='cork';
0072     end
0073 end
0074 isgts=0;
0075 if(strcmp(exename,'gtsset'))
0076     isgts=1;
0077 end
0078 
0079 exesuff=fallbackexeext(getexeext,exename);
0080 randseed=hex2dec('623F9A9E'); % "U+623F U+9A9E"
0081 if(~isempty(getvarfrom({'caller','base'},'ISO2MESH_RANDSEED')))
0082     randseed=getvarfrom({'caller','base'},'ISO2MESH_RANDSEED');
0083 end
0084 
0085 for i=1:3:len
0086    op=varargin{i};
0087    no=varargin{i+1};
0088    el=varargin{i+2};
0089    opstr=op;
0090    if(strcmp(op,'or'))   opstr='union'; end
0091    if(strcmp(op,'xor'))  opstr='all';   end
0092    if(strcmp(op,'and'))
0093         if(isgts)
0094            opstr='inter'; 
0095         else
0096        opstr='isct'; 
0097     end
0098    end
0099    if(strcmp(op,'-'))    opstr='diff';  end
0100    if(strcmp(op,'self')) opstr='inter -s';  end
0101    if(isgts && (strcmp(op,'first') || strcmp(op,'second') || strcmp(op,'+')))
0102        opstr='all';
0103    end
0104    if(strcmp(opstr,'all') && isgts==0)
0105        opstr='resolve';
0106    end
0107    tempsuff='gts';
0108    if(isgts==0)
0109        tempsuff='off';
0110    end
0111    deletemeshfile(mwpath(['pre_surfbool*.' tempsuff]));
0112    deletemeshfile(mwpath('post_surfbool.off'));
0113    if(strcmp(opstr,'all'))
0114       deletemeshfile(mwpath('s1out2.off'));
0115       deletemeshfile(mwpath('s1in2.off'));
0116       deletemeshfile(mwpath('s2out1.off'));
0117       deletemeshfile(mwpath('s2in1.off'));
0118    end
0119    if(strcmp(op,'decouple'))
0120        if(exist('node1','var')==0)
0121           node1=node;
0122           elem1=elem;
0123           newnode(:,4)=1;
0124           newelem(:,4)=1;
0125        end
0126        opstr=[' --decouple-inin 1 --shells 2']; %-q
0127        saveoff(node1(:,1:3),elem1(:,1:3),mwpath('pre_decouple1.off'));
0128        if(size(no,2)~=3)
0129            opstr=['-q --shells ' num2str(no)];
0130            cmd=sprintf('cd "%s" && "%s%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
0131                mwpath('pre_decouple1.off'),opstr);
0132        else
0133            saveoff(no(:,1:3),el(:,1:3),mwpath('pre_decouple2.off'));
0134            cmd=sprintf('cd "%s" && "%s%s" "%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
0135                mwpath('pre_decouple1.off'),mwpath('pre_decouple2.off'),opstr);
0136        end
0137    elseif(strcmp(op,'decoupleout'))
0138        if(exist('node1','var')==0)
0139           node1=node;
0140           elem1=elem;
0141           newnode(:,4)=1;
0142           newelem(:,4)=1;
0143        end
0144        opstr=[' --decouple-outout 1 --shells 2']; %-q
0145        saveoff(node1(:,1:3),elem1(:,1:3),mwpath('pre_decouple1.off'));
0146        if(size(no,2)~=3)
0147            opstr=['-q --shells ' num2str(no)];
0148            cmd=sprintf('cd "%s" && "%s%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
0149                mwpath('pre_decouple1.off'),opstr);
0150        else
0151            saveoff(no(:,1:3),el(:,1:3),mwpath('pre_decouple2.off'));
0152            cmd=sprintf('cd "%s" && "%s%s" "%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
0153                mwpath('pre_decouple1.off'),mwpath('pre_decouple2.off'),opstr);
0154        end
0155    elseif(strcmp(op,'separate'))
0156        if(exist('node1','var')==0)
0157           node1=node;
0158           elem1=elem;
0159           newnode(:,4)=1;
0160           newelem(:,4)=1;
0161        end
0162        opstr=[' --shells 2']; %-q
0163        saveoff(node1(:,1:3),elem1(:,1:3),mwpath('pre_decouple1.off'));
0164        if(size(no,2)~=3)
0165            opstr=['-q --shells ' num2str(no)];
0166            cmd=sprintf('cd "%s" && "%s%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
0167                mwpath('pre_decouple1.off'),opstr);
0168        else
0169            saveoff(no(:,1:3),el(:,1:3),mwpath('pre_decouple2.off'));
0170            cmd=sprintf('cd "%s" && "%s%s" "%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
0171                mwpath('pre_decouple1.off'),mwpath('pre_decouple2.off'),opstr);
0172        end       
0173    else
0174        if(isgts)
0175            savegts(newnode(:,1:3),newelem(:,1:3),mwpath(['pre_surfbool1.' tempsuff]));
0176            savegts(no(:,1:3),el(:,1:3),mwpath(['pre_surfbool2.' tempsuff]));
0177            cmd=sprintf('cd "%s" && "%s%s" %s "%s" "%s" -v > "%s"',mwpath,mcpath('gtsset'),exesuff,...
0178                opstr,mwpath('pre_surfbool1.gts'),mwpath('pre_surfbool2.gts'),mwpath('post_surfbool.off'));
0179        else
0180            saveoff(newnode(:,1:3),newelem(:,1:3),mwpath(['pre_surfbool1.' tempsuff]));
0181            saveoff(no(:,1:3),el(:,1:3),mwpath(['pre_surfbool2.' tempsuff]));
0182            cmd=sprintf('cd "%s" && "%s%s" %s%s "%s" "%s" "%s" -%d',mwpath,mcpath(exename),exesuff,'-',...
0183                opstr,mwpath(['pre_surfbool1.' tempsuff]),mwpath(['pre_surfbool2.' tempsuff]),mwpath('post_surfbool.off'),randseed);
0184        end
0185    end
0186    [status outstr]=system(cmd);
0187    if(status~=0 && strcmp(op,'self')==0)
0188        error(sprintf('surface boolean command failed:\n%s\nERROR: %s\n',cmd,outstr));
0189    end
0190    if(status~=0 && strcmp(op,'self') && ~isempty(strfind(outstr,'(new_ear): assertion failed')))
0191        fprintf(1,'no self-intersection was found! (ignore the above error)\n');
0192        newnode=[];
0193        newelem=[];
0194        newelem0=[];
0195        return;
0196    end
0197    if(strcmp(opstr,'all'))
0198       % tag the 4 piceses of meshes, this tag do not propagate to the next boolean operation
0199       [nnode nelem]=readoff(mwpath('s1out2.off'));
0200       newelem=[nelem ones(size(nelem,1),1)];
0201       newnode=[nnode ones(size(nnode,1),1)];
0202 
0203       [nnode nelem]=readoff(mwpath('s1in2.off'));
0204       newelem=[newelem; nelem+size(newnode,1) 3*ones(size(nelem,1),1)];
0205       newnode=[newnode; nnode 3*ones(size(nnode,1),1)];
0206 
0207       [nnode nelem]=readoff(mwpath('s2out1.off'));
0208       newelem=[newelem; nelem+size(newnode,1) 2*ones(size(nelem,1),1)];
0209       newnode=[newnode; nnode 2*ones(size(nnode,1),1)];
0210 
0211       [nnode nelem]=readoff(mwpath('s2in1.off'));
0212       newelem=[newelem; nelem+size(newnode,1) 4*ones(size(nelem,1),1)];
0213       newnode=[newnode; nnode 4*ones(size(nnode,1),1)];
0214       if(isgts)
0215           if(strcmp(op,'first'))
0216               newelem=newelem(find(mod(newelem(:,4),2)==1),:);
0217               [newnode,nelem]=removeisolatednode(newnode,newelem(:,1:3));
0218               newelem=[nelem newelem(:,4)];
0219           elseif(strcmp(op,'second'))
0220               newelem=newelem(find(mod(newelem(:,4),2)==0),:);
0221               [newnode,nelem]=removeisolatednode(newnode,newelem(:,1:3));
0222               newelem=[nelem,newelem(:,4)];
0223           end
0224       end
0225    elseif(strcmp(op,'decouple'))
0226       [newnode,newelem]=readoff(mwpath('pre_decouple1_fixed.off')); %[node1,elem1]
0227       %newelem=[newelem;elem1+size(newnode,1) (i+1)*ones(size(elem1,1),1)];
0228       %newnode=[newnode;node1 (i+1)*ones(size(node1,1),1)];
0229    elseif(strcmp(op,'separate'))
0230        [newnode,newelem]=readoff(mwpath('pre_decouple1_fixed.off'));
0231    elseif(strcmp(op,'decoupleout'))
0232       [newnode,newelem]=readoff(mwpath('pre_decouple1_fixed.off')); %[node1,elem1]
0233    else
0234       [newnode,newelem]=readoff(mwpath('post_surfbool.off'));
0235       if(strcmp(op,'self'))
0236           fprintf(1,'a total of %d self-intersecting elements were found\n',size(newelem,1));
0237           if(nargout>=3)
0238               [found,newelem0]=ismember(newnode,node,'rows');
0239               if(~all(found))
0240                   error('self intersecting elements contain new nodes');
0241               end
0242               newelem0=newelem0(newelem);
0243           end
0244           return;
0245       end
0246    end
0247 end