gmplot.m

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

function gmplot(mesh,soln,solnrange,cmap);
% gmplot(mesh, soln {, solnrange {, cmap}})
%  Takes a simplicial complex and a solution vector computed
%   by gmfem for example and makes a color plot of the solution.
%  Simpcomplex have intrinsic dimension of 2, and can be either 2 or 3D.
%  It does "manual" interpolated shading, i.e., it breaks each triangle
%  into 9 little triangles and colors each little triangle flat.
%  This is because interpolated colors in matlab don't seem to
%  well on my printer.
%  Colormap is jet unless a cmap argument is given.
%  Third optional argument is an upper and lower bound for the
%  the solution range; this range is used to map the solution into
%  the color map.  This is useful if you want to make several plots
%  and fix the absolute meaning of the colors.  If omitted, the
%  mapping automatically goes from the lowest color to the highest
%  color in the colormap.


if ~isa(mesh,'zba')
  error('First argument must be simplicial complex');
end

if length(mesh) < 1
  error('First argument must be simplicial complex');
end

global GM_SIMPCOMP_TYPE_CODE
typecode = double(mesh{0});
if ~strcmpi(typecode, GM_SIMPCOMP_TYPE_CODE)
  error('First argument must be simplicial complex');
end

gdim = double(mesh{1});
if gdim ~= 2
  error('Only 2D simplicial complexes supported')
end

di = double(mesh{2});

xy = mesh{4}(1:di, :)';
trilist0 = zba([]);
[scrap,numtoplev] = size(mesh{7});
for i = 0 : numtoplev - 1
  trilist0 = [trilist0; mesh{7}{1,i}'];
end
[numvtx,scrap] = size(xy);
[numtri,scrap] = size(trilist0);


vtx_global_to_ordinal = ...
   sparse((double(mesh{4}(0,:)))'+1,ones(numvtx,1), (0:numvtx-1)');
trilist = zba(full(vtx_global_to_ordinal(double(trilist0)+1)));


[numvtx1, n] = size(soln);
if numvtx1 ~= numvtx | n ~= 1
  error('Size of soln argument is incompatible with simpcomplex')
end

[engine, filename, filestat, notify, format, ...
 thickness, default_color, beziersub] = gmvizgui_get_state(di);

if nargin < 4
  cmap = jet(64);
end


%% Refine the mesh by a factor of 3.

refine = 3;
h = 1 / refine;

th = 1/3;
tth = 2/3;

re1 = refine - 1;
re2 = refine - 2;

x1 = double(xy(trilist(:,0),0));
x2 = double(xy(trilist(:,1),0));
x3 = double(xy(trilist(:,2),0));

y1 = double(xy(trilist(:,0),1));
y2 = double(xy(trilist(:,1),1));
y3 = double(xy(trilist(:,2),1));

if di == 3
  z1 = double(xy(trilist(:,0),2));
  z2 = double(xy(trilist(:,1),2));
  z3 = double(xy(trilist(:,2),2));
end

so1 = double(soln(trilist(:,0)));
so2 = double(soln(trilist(:,1)));
so3 = double(soln(trilist(:,2)));

pts = [];
simps = [];
fval = [];
for i = 0 : re1
  for j = 0 : re1 - i;
    k = re1 - i - j;
    trix1 = ((i+1) * x1 + j * x2 + k * x3) * h;
    trix2 = (i * x1 + (j+1) * x2 + k * x3) * h;
    trix3 = (i * x1 + j * x2 + (k+1) * x3) * h;
    triy1 = ((i+1) * y1 + j * y2 + k * y3) * h;
    triy2 = (i * y1 + (j+1) * y2 + k * y3) * h;
    triy3 = (i * y1 + j * y2 + (k+1) * y3) * h;
    if di == 3
      triz1 = ((i+1) * z1 + j * z2 + k * z3) * h;
      triz2 = (i * z1 + (j+1) * z2 + k * z3) * h;
      triz3 = (i * z1 + j * z2 + (k+1) * z3) * h;
    end
    func = ((i+th) * so1 + (j+th) * so2 + (k+th) * so3) * h;
    [base,scrap] = size(pts);
    if di == 3
      pts = [pts;[trix1,triy1,triz1];[trix2,triy2,triz2];[trix3,triy3,triz3]];
    else
      pts = [pts;[trix1,triy1];[trix2,triy2];[trix3,triy3]];
    end
    simps = [simps;...
         [(base:base+numtri-1)', (base+numtri:base+2*numtri-1)', ...
            (base+2*numtri:base+3*numtri-1)']];
    fval = [fval,func'];
  end
end

for i = 0 : re2
  for j = 0 : re2 - i;
    k = re2 - i - j;
    trix1 = ((i+1) * x1 + (j+1) * x2 + k * x3) * h;
    trix2 = (i * x1 + (j+1) * x2 + (k+1) * x3) * h;
    trix3 = ((i+1) * x1 + j * x2 + (k+1) * x3) * h;
    triy1 = ((i+1) * y1 + (j+1) * y2 + k * y3) * h;
    triy2 = (i * y1 + (j+1) * y2 + (k+1) * y3) * h;
    triy3 = ((i+1) * y1 + j * y2 + (k+1) * y3) * h;
    if di == 3
      triz1 = ((i+1) * z1 + (j+1) * z2 + k * z3) * h;
      triz2 = (i * z1 + (j+1) * z2 + (k+1) * z3) * h;
      triz3 = ((i+1) * z1 + j * z2 + (k+1) * z3) * h;
    end
    [base,scrap] = size(pts);
    if di == 3
      pts = [pts;[trix1,triy1,triz1];[trix2,triy2,triz2];[trix3,triy3,triz3]];
    else
      pts = [pts;[trix1,triy1];[trix2,triy2];[trix3,triy3]];
    end
    simps = [simps;...
         [(base:base+numtri-1)', (base+numtri:base+2*numtri-1)', ...
            (base+2*numtri:base+3*numtri-1)']];
    fval = [fval,func'];
  end
end



if nargin < 3
  umin = min(fval);
  umax = max(fval);
else
  umin = solnrange(1);
  umax = solnrange(2);
end

% Get colors corresponding to u.

spread = umax - umin;
if spread == 0
  spread = 1;
end

[newnumtri,scrap] = size(simps);
[cmaplength, scrap] = size(cmap);
multiplier = (cmaplength - 1) / spread;
unewtri_int = floor((fval - umin) * multiplier) + 1;
if any(unewtri_int < 1)
  error('Some field values go below the minimum field range');
end

if any(unewtri_int > cmaplength)
  error('Some field values go above the maximum field range');
end

colors = [cmap(unewtri_int,1), cmap(unewtri_int,2), cmap(unewtri_int,3), ...
          ones(newnumtri,1)];

bbox = 1e307 * [ones(di,1),-ones(di,1)];
bbox = gmviz_update_bbox(bbox, pts);


if strcmpi(engine,'local')
  handle1 = gmviz_init_matlabplot(bbox);
  caxis([umin,umax]);
  gmviz_matlabplot(handle1, bbox, pts, simps, colors, 2, thickness);
  gmviz_terminate_matlabplot(handle1);
else
  handle1 = gmviz_init_vrmlplot(bbox, filename, filestat);
  gmviz_vrmlplot(handle1, bbox, pts, simps, colors, 2, thickness);
  gmviz_terminate_vrmlplot(handle1, filename, notify);
end


% ------------------------------------------------------------------
% 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
% ------------------------------------------------------------------