gmq11read.m

matlab有限元分析工具,比经较全面的一个手册,请大家下载呀

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 file


fi = fopen(filename, 'r');
if fi < 0
  error(sprintf('File %s not found', filename))
end
sbrep = fread(fi);
fclose(fi);

% strip comments

sbrep = 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)];
end
ca1 = split_bracket_list(sbrep, 2);

% process header info

if ~strcmpi(ca1{1},'brep')
  error('Keyword ''brep'' missing from file header');
end
ca2 = 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')
end
faces = 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 dim
for 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};
end
if ~strcmpi(faces, 'nil')
  error('Brep does not terminate with nil');
end
b = 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
  return
end
leftdelim = (s1 == '(') | (s1 == '<');
rightdelim = (s1 == ')') | (s1 == '>');
c = cumsum(leftdelim - rightdelim);
if any(c < 0) | c(end) ~= 0
  error('Unbalanced delimiters')
end

issp1 = 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);
  end
end


function 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');
end
split = split_helper(s1(2:length(s1) - 1));
if length(split) ~= expected_size
  error('Wrong number of entries in bracket list');
end

function split = split_paren_list(string1, expected_size)
s1 = gm_strim(string1);
if s1(1) ~= '(' | s1(end) ~= ')'
  error('Parenthesized list expected but delimiters ( ) not found');
end
split = 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
% ------------------------------------------------------------------