function b = gmq11read(filename);% gmq11read: Read in a brep in QMG 1.1 ascii format%  b = gmq11read(filename);global GM_BREP_TYPE_CODE;% read in filefi = fopen(filename, 'r');if fi < 0  error(sprintf('File %s not found', filename))endsbrep = fread(fi);fclose(fi);% strip commentssbrep = char([10,sbrep']);while 1  p = findstr(sbrep, char([10,'#']));  if length(p) == 0 | length(sbrep) < 2    break  end  p = p(1);  p1 = findstr(sbrep(p + 1:end), char(10)) + p;  if length(p1) == 0 | length(sbrep) <= p    p1 = length(sbrep);  else    p1 = p1(1);  end  sbrep = [sbrep(1:p),sbrep(p1 + 1:end)];endca1 = split_bracket_list(sbrep, 2);% process header infoif ~strcmpi(ca1{1},'brep')  error('Keyword ''brep'' missing from file header');endca2 = split_bracket_list(ca1{2}, 3);gdim = sscanf(ca2{1},'%d');di = sscanf(ca2{2}, '%d');if length(gdim) ~= 1 | length(di) ~= 1 | ...     gdim < 0 | gdim > di | di < 2 | di > 3  error('intrinsic or embedded dimension out of range')endfaces = ca2{3};cplist = [];segtable = {};facecount = 0;allfacenames = {};allfspecs = [];newbrep0 = cell(6 + gdim, 1);newbrep0{1} = GM_BREP_TYPE_CODE;newbrep0{2} = gdim;newbrep0{3} = di;newbrep0{4} = {};% loop over faces of dimension dimfor dim = 0 : gdim  faces = split_bracket_list(faces, 2);  thislev = faces{1};  thislevfaces = split_paren_list(thislev);  newlev = {};  numface = length(thislevfaces);  for faceind = 1 : numface    thisface = split_bracket_list(thislevfaces{faceind}, 4);    facename = thisface{1};    facecount = facecount + 1;    allfacenames{facecount} = facename;    allfspecs(facecount,1:2) = [dim, faceind-1];    pvlist = split_paren_list(thisface{2});    chlist = split_paren_list(thisface{3});    iblist = split_paren_list(thisface{4});    newpv = {};    pointfound = 0;    newpvcount = 0;    % process property-value list.  Must save 'point' property    % for faces of dimension 0.  Else scrap point, affine_coef, and    % affine_rhs properties.  Also, reformat color property.    for j = 1 : length(pvlist)      thispv = split_bracket_list(pvlist{j},2);      if strcmpi(thispv{1}, 'point')        if dim == 0          if pointfound            error('Point property occurs twice in vertex');          end          pointfound = 1;          cval = split_paren_list(thispv{2});          if length(cval) ~= di            error('Number of coordinate entries does not match embedded dim');          end          cp = zeros(di,1);          for k = 1 : di            v = sscanf(cval{k}, '%f');            if length(v) ~= 1              error('Could not convert coordinate entry in point prop-val to double');            end            cp(k) = v;          end          cplist = [cplist,cp];        end      elseif strcmpi(thispv{1},'affine_coef')      elseif strcmpi(thispv{1},'affine_rhs')      elseif strcmpi(thispv{1},'color')        val1 = split_bracket_list(thispv{2}, 4);        newpvcount = newpvcount + 1;        newpv{1,newpvcount} = 'color';        newpv{2,newpvcount} = ['(',val1{1},' ', val1{2},' ', val1{3},' ', ...                           val1{4}, ')'];      else        newpvcount = newpvcount + 1;        newpv{1,newpvcount} = thispv{1};        newpv{2,newpvcount} = thispv{2};      end    end    if dim == 0 & ~pointfound      error('No point property for vertex');    end    chcount = length(chlist);    newchlist = cell(1,chcount);    newchlistn = zeros(1,chcount);    % process child and ib lists by seeking names in table.    for j = 1 : chcount;      i = strmatch(chlist{j}, allfacenames, 'exact');      if length(i) == 0        disp(sprintf('For face named %s, subface named %s', facename, chlist{j}));        error('Face named as subface not found')      end      if length(i) > 1        disp(sprintf('For face named %s, subface named %s', facename, chlist{j}));        error('Face named as subface occurs is the name of multiple faces')      end      subfdim = allfspecs(i,1);      subfind = allfspecs(i,2);      if subfdim ~= dim - 1        disp(sprintf('For face named %s, subface named %s', facename, chlist{j}));        error('Subface is not one dimension lower')      end      newchlist{j} = chlist{j};      newchlistn(j) = subfind;    end    ibcount = 0;    newiblist = {};    for j = 1 : length(iblist)      i = strmatch(iblist{j}, allfacenames, 'exact');      if length(i) == 0        disp(sprintf('For face named %s, subface named %s', facename, chlist{j}));        error('Face named as subface not found')      end      if length(i) > 1        disp(sprintf('For face named %s, subface named %s', facename, chlist{j}));        error('Face named as subface occurs is the name of multiple faces')      end      subfdim = allfspecs(i,1);      subfind = allfspecs(i,2);      if subfdim > dim - 1        disp(sprintf('For face named %s, subface named %s', facename, chlist{j}));        error('Subface is not one dimension or more dimensions lower')      elseif subfdim == dim - 1        chcount = chcount + 1;        newchlist{chcount} = iblist{j};        newchlistn(chcount) = subfind;        chcount = chcount + 1;        newchlist{chcount} = iblist{j};        newchlistn(chcount) = subfind;      else        ibcount = ibcount + 1;        newiblist{ibcount} = iblist{j};      end    end      % prepare bezier list    if dim == 0      bezlist = {'vertex'; []; [faceind - 1]};    elseif dim == 1      % In 1-D case, sort endpoints along edge and generate      % linear segments.      if chcount == 0 | rem(chcount,2) ~= 0        error('Edge must have an even number of boundaries');      end      maxdelta = -1;      bestd = -1;      newchlistn1 = newchlistn + ones(1,chcount);      for j = 1 : di        lb = min(cplist(j, newchlistn1));        ub = max(cplist(j, newchlistn1));        if ub - lb > maxdelta          maxdelta = ub - lb;          bestd = j;        end      end      [scrap,pos] = sort(cplist(bestd, newchlistn1));      bezlist = cell(2,chcount / 2);      seglist = zeros(chcount / 2, 2);      for j = 1 : 2 : chcount        k = (j + 1) / 2;        bezlist{1,k} = 'bezier_curve';        bezlist{2,k} = 1;        bezlist{3,k} = [newchlistn(pos(j)),newchlistn(pos(j+1))];         seglist(k,:) = [newchlistn(pos(j)),newchlistn(pos(j+1))];       end      segtable{faceind} = seglist;    elseif dim == 2 & di == 3           % In 2-D case, call Delaunay triangulation routine to      % make triangular patches      smallcoordcount = 0;      smallcoordlist = zba([]);      smallseglist = [];      [scrap,numcp] = size(cplist);      big2smallmap = zeros(numcp,1);      small2bigmap = [];      for j = 1 : chcount        subfind = newchlistn(j);        [numseg,scrap] = size(segtable{subfind + 1});        for k = 1 : numseg          thisseg = zeros(1,2);          for m = 1 : 2            ep = segtable{subfind + 1}(k,m) + 1;            if big2smallmap(ep) == 0              smallcoordcount = smallcoordcount + 1;              big2smallmap(ep) = smallcoordcount;              small2bigmap(smallcoordcount) = ep;              smallcoordlist(smallcoordcount - 1,:) = cplist(:,ep)';            end            thisseg(m) = big2smallmap(ep) - 1;          end          smallseglist = [smallseglist;thisseg];        end      end      for j = 1 : ibcount        ep = newibist(2,j) + 1;        if big2smallmap(ep) == 0          smallcoordcount = smallcoordcount + 1;          big2smallmap(ep) = smallcoordcount;          small2bigmap(smallcoordcount) = ep;          smallcoordlist(smallcoordcount - 1,:) = cplist(:,ep)';        end      end      triangles = gm_polytri(smallcoordlist, smallseglist);      [numtri, scrap] = size(triangles);      bezlist = cell(3, numtri);      for j = 1 : numtri        bezlist{1,j} = 'bezier_triangle';        bezlist{2,j} = 1;        bezlist{3,j} = small2bigmap(triangles(j,:)+[1,1,1]) - [1,1,1];      end    else      bezlist = {};    end    newlev{1,faceind} = facename;    newlev{2,faceind} = newpv;    newlev{3,faceind} = newchlist;    newlev{4,faceind} = newiblist;    newlev{5,faceind} = bezlist;  end  if dim == 0    newbrep0{5} = cplist;  end  newbrep0{6 + dim} = newlev;  faces = faces{2};endif ~strcmpi(faces, 'nil')  error('Brep does not terminate with nil');endb = zba(newbrep0);% Function to split a list, taking into account delimiters% () and <>function split = split_helper(string1)s1 = gm_strim(string1);split = {};if length(s1) == 0  returnendleftdelim = (s1 == '(') | (s1 == '<');rightdelim = (s1 == ')') | (s1 == '>');c = cumsum(leftdelim - rightdelim);if any(c < 0) | c(end) ~= 0  error('Unbalanced delimiters')endissp1 = isspace(s1);prev = 1;count = 1;while 1  m = find(~issp1(prev:end));  if length(m) == 0    return  end  st = prev + m(1) - 1;  m2 = find(issp1(st:end) & c(st:end) == 0);  if length(m2) == 0    split{count} = s1(st:end);    return;  else    split{count} = s1(st:st + m2(1) - 2);    count = count + 1;    prev = st + m2(1);  endendfunction split = split_bracket_list(string1, expected_size)s1 = gm_strim(string1);if s1(1) ~= '<' | s1(end) ~= '>'  error('Angle bracket list expected but delimiters < > not found');endsplit = split_helper(s1(2:length(s1) - 1));if length(split) ~= expected_size  error('Wrong number of entries in bracket list');endfunction split = split_paren_list(string1, expected_size)s1 = gm_strim(string1);if s1(1) ~= '(' | s1(end) ~= ')'  error('Parenthesized list expected but delimiters ( ) not found');endsplit = split_helper(s1(2:length(s1) - 1));% ------------------------------------------------------------------% Copyright (c) 1999 by Cornell University.  All rights reserved% See the accompanying file 'Copyright' for authorship information,% the terms of the license governing this software, and disclaimers% concerning this software.% ------------------------------------------------------------------% This file is part of the QMG software.  % Version 2.0 of QMG, release date September 3, 1999% ------------------------------------------------------------------