Description of savesurfpoly

0001 function savesurfpoly(v,f,holelist,regionlist,p0,p1,fname,forcebox)
0002 %
0003 % savesurfpoly(v,f,holelist,regionlist,p0,p1,fname)
0004 %
0005 % save a set of surfaces into poly format (for tetgen)
0006 %
0007 % author: Qianqian Fang, <q.fang at neu.edu>
0008 % date: 2007/11/21
0009 %
0010 % input:
0011 %      v: input, surface node list, dimension (nn,3)
0012 %         if v has 4 columns, the last column specifies mesh density near each node
0013 %      f: input, surface face element list, dimension (be,3)
0014 %      holelist: list of holes, each hole is represented by an internal point
0015 %      regionlist: list of regions, similar to holelist
0016 %      p0: coordinate of one of the end of the bounding box
0017 %      p1: coordinate for the other end of the bounding box
0018 %      fname: output file name
0019 %      forcebox: non-empty: add bounding box, []: automatic
0020 %                if forcebox is a 8x1 vector, it will be used to
0021 %                specify max-edge size near the bounding box corners
0022 %
0023 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0024 %
0025 dobbx=0;
0026 if(nargin>=8)
0027     dobbx=~isempty(forcebox) & all(forcebox);
0028 end
0029 
0030 if(~iscell(f) && size(f,2)==4)
0031     faceid=f(:,4);
0032     f=f(:,1:3);
0033 end
0034 
0035 if(~iscell(f))
0036     edges=surfedge(f);
0037 else
0038     edges=[];
0039 end
0040 bbxnum=0;
0041 
0042 nodesize=[];
0043 if(size(v,2)==4)
0044    nodesize=v(:,4);
0045    v=v(:,1:3);
0046 end
0047 node=v;
0048 loopid=[];
0049 loopvert={};
0050 loopnum=1;
0051 if(~isempty(edges))
0052     loops=extractloops(edges);
0053     if(length(loops)<3)
0054         error('degenerated loops detected');
0055     end
0056     seg=[0,find(isnan(loops))];
0057     segnum=length(seg)-1;
0058     newloops=[];
0059     for i=1:segnum
0060        if(seg(i+1)-(seg(i)+1)==0)
0061            continue;
0062        end
0063        oneloop=loops(seg(i)+1:seg(i+1)-1);
0064        if(oneloop(1)==oneloop(end)) oneloop(end)=[]; end
0065        newloops=[newloops nan bbxflatsegment(node,oneloop)];
0066     end
0067     loops=[newloops nan];
0068 
0069     seg=[0,find(isnan(loops))];
0070     segnum=length(seg)-1;
0071     bbxnum=6;
0072     loopcount=zeros(bbxnum,1);
0073     loopid=zeros(segnum,1);
0074     for i=1:segnum     % walk through the edge loops
0075         subloop=loops(seg(i)+1:seg(i+1)-1);
0076         if(isempty(subloop))
0077             continue;
0078         end
0079         loopvert{loopnum}=subloop;
0080         loopnum=loopnum+1;
0081         boxfacet=find(sum(abs(diff(v(subloop,:))))<1e-8); % find a flat loop
0082         if(length(boxfacet)==1)   % if the loop is flat along x/y/z dir
0083             bf=boxfacet(1);    % no degeneracy allowed
0084             if(sum(abs(v(subloop(1),bf)-p0(bf)))<1e-2)
0085                 loopcount(bf)=loopcount(bf)+1;
0086                 v(subloop,bf)=p0(bf);
0087                 loopid(i)=bf;
0088             elseif(sum(abs(v(subloop(1),bf)-p1(bf)))<1e-2)
0089                 loopcount(bf+3)=loopcount(bf+3)+1;
0090                 v(subloop,bf)=p1(bf);
0091                 loopid(i)=bf+3;
0092             end
0093         end
0094     end
0095 end
0096 
0097 if(dobbx && isempty(edges))
0098     bbxnum=6;
0099     loopcount=zeros(bbxnum,1);    
0100 end
0101 
0102 if(dobbx || ~isempty(edges))
0103     nn=size(v,1);
0104     boxnode=[p0;p1(1),p0(2:3);p1(1:2),p0(3);p0(1),p1(2),p0(3);
0105               p0(1:2),p1(3);p1(1),p0(2),p1(3);p1;p0(1),p1(2:3)];
0106     boxelem=[
0107         4 nn nn+3 nn+7 nn+4;   % x=xmin
0108         4 nn nn+1 nn+5 nn+4;   % y=ymin
0109         4 nn nn+1 nn+2 nn+3;   % z=zmin
0110         4 nn+1 nn+2 nn+6 nn+5; % x=xmax
0111         4 nn+2 nn+3 nn+7 nn+6; % y=ymax
0112         4 nn+4 nn+5 nn+6 nn+7];% z=zmax
0113 
0114     node=[v;boxnode];
0115 end
0116 
0117 node=[(0:size(node,1)-1)',node];
0118 
0119 fp=fopen(fname,'wt');
0120 fprintf(fp,'#node list\n%d 3 0 0\n',length(node));
0121 fprintf(fp,'%d %.16f %.16f %.16f\n',node');
0122 
0123 if(~iscell(f))
0124     fprintf(fp,'#facet list\n%d 1\n',length(f)+bbxnum+length(loopvert));
0125     elem=[3*ones(length(f),1),f-1];
0126     if(~isempty(elem))
0127         if(exist('faceid','var') && length(faceid)==size(elem,1))
0128             fprintf(fp,'1 0 %d\n%d %d %d %d\n',[faceid(:) elem]');
0129         else
0130             fprintf(fp,'1 0\n%d %d %d %d\n',elem');
0131         end
0132     end
0133     if(~isempty(loopvert))
0134         for i=1:length(loopvert)     % walk through the edge loops
0135             subloop=loopvert{i}-1;
0136             fprintf(fp,'1 0 %d\n%d',i,length(subloop));
0137             fprintf(fp,'\t%d',subloop);
0138             fprintf(fp,'\n');
0139         end
0140     end
0141 else % if the surface is recorded as a cell array
0142     totalplc=0;
0143     for i=1:length(f)
0144         if(~iscell(f{i}))
0145             totalplc=totalplc+size(f{i},1);
0146         else
0147             totalplc=totalplc+size(f{i}{1},1);
0148         end
0149     end
0150     fprintf(fp,'#facet list\n%d 1\n',totalplc+bbxnum);
0151     for i=1:length(f)
0152         plcs=f{i};
0153         faceid=-1;
0154         if(iscell(plcs)) % if each face is a cell, use plc{2} for face id
0155             if(length(plcs)>1)
0156                 faceid=plcs{2};
0157             end
0158             plcs=plcs{1};
0159         end
0160         for row=1:size(plcs,1);
0161          plc=plcs(row,:);
0162          if(any(isnan(plc))) % we use nan to separate outter contours and holes
0163             holeid=find(isnan(plc));
0164             if(faceid>0) 
0165                 fprintf(fp,'%d %d %d\n%d',length(holeid)+1,length(holeid),faceid,holeid(1)-1);
0166             else
0167                 fprintf(fp,'%d %d\n%d',length(holeid)+1,length(holeid),holeid(1)-1);
0168             end
0169             fprintf(fp,'\t%d',plc(1:holeid(1)-1)-1);
0170             fprintf(fp,'\t1\n');
0171             for j=1:length(holeid)
0172                 if(j==length(holeid))
0173                     fprintf(fp,'%d',length(plc(holeid(j)+1:end)));
0174                 fprintf(fp,'\t%d',plc(holeid(j)+1:end)-1);
0175                 else
0176                     fprintf(fp,'%d',length(plc(holeid(j)+1:holeid(j+1)-1)));
0177                     fprintf(fp,'\t%d',plc(holeid(j)+1:holeid(j+1)-1)-1);
0178                 end
0179                 fprintf(fp,'\t1\n');
0180             end
0181             for j=1:length(holeid)
0182                 if(j==length(holeid))
0183                     fprintf(fp,'%d %.16f %.16f %.16f\n',j,mean(node(plc(holeid(j)+1:end),2:4)));
0184                 else
0185                     fprintf(fp,'%d %.16f %.16f %.16f\n',j,mean(node(plc(holeid(j)+1:holeid(j+1)-1),2:4)));
0186                 end
0187             end
0188          else
0189              if(faceid>0)
0190                  fprintf(fp,'1 0 %d\n%d',faceid,length(plc));
0191              else
0192                  fprintf(fp,'1 0\n%d',length(plc));
0193              end
0194              fprintf(fp,'\t%d',plc-1);
0195              fprintf(fp,'\t1\n');
0196          end
0197         end
0198     end
0199 end
0200 
0201 if(dobbx || ~isempty(edges))
0202     for i=1:bbxnum
0203         fprintf(fp,'%d %d 1\n',1+loopcount(i),loopcount(i));
0204         fprintf(fp,'%d %d %d %d %d\n',boxelem(i,:));
0205         if(~isempty(edges) && loopcount(i) &&~isempty(find(loopid==i, 1)))
0206             endid=find(loopid==i);
0207             for k=1:length(endid)
0208                 j=endid(k);
0209                 subloop=loops(seg(j)+1:seg(j+1)-1);
0210                 fprintf(fp,'%d ',length(subloop));
0211                 fprintf(fp,'%d ',subloop-1);
0212                 fprintf(fp,'\n');
0213             end
0214             for k=1:length(endid)
0215                 j=endid(k);
0216                 subloop=loops(seg(j)+1:seg(j+1)-1);
0217                 fprintf(fp,'%d %.16f %.16f %.16f\n',k,internalpoint(v,subloop)); %mean(v(subloop,:)));
0218             end
0219         end
0220     end
0221 end
0222 
0223 if(size(holelist,1))
0224         fprintf(fp,'#hole list\n%d\n',size(holelist,1));
0225         for i=1:size(holelist,1)
0226                 fprintf(fp,'%d %.16f %.16f %.16f\n', i, holelist(i,:));
0227         end
0228 else
0229     fprintf(fp,'#hole list\n0\n');
0230 end
0231 
0232 if(size(regionlist,1))
0233     fprintf(fp,'#region list\n%d\n',size(regionlist,1));
0234     if(size(regionlist,2)==3)
0235         for i=1:size(regionlist,1)
0236             fprintf(fp,'%d %.16f %.16f %.16f %d\n', i, regionlist(i,:),i);
0237         end
0238     elseif(size(regionlist,2)==4)
0239         for i=1:size(regionlist,1)
0240             fprintf(fp,'%d %.16f %.16f %.16f %d %.16f\n', i, regionlist(i,1:3),i,regionlist(i,4));
0241         end
0242     end
0243 end
0244 fclose(fp);
0245 
0246 if(~isempty(nodesize))
0247     if(size(nodesize,1)+size(forcebox(:),1)==size(node,1))
0248         nodesize=[nodesize;forcebox(:)];
0249     end
0250     fid=fopen(regexprep(fname,'\.poly$','.mtr'),'wt');
0251     fprintf(fid,'%d 1\n',size(nodesize,1));
0252     fprintf(fid,'%.16f\n',nodesize);
0253     fclose(fid);
0254 end
savesurfpoly

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
