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');endif length(mesh) < 1  error('First argument must be simplicial complex');endglobal GM_SIMPCOMP_TYPE_CODEtypecode = double(mesh{0});if ~strcmpi(typecode, GM_SIMPCOMP_TYPE_CODE)  error('First argument must be simplicial complex');endgdim = double(mesh{1});if gdim ~= 2  error('Only 2D simplicial complexes supported')enddi = 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));endso1 = 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'];  endendfor 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'];  endendif 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');endif any(unewtri_int > cmaplength)  error('Some field values go above the maximum field range');endcolors = [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% ------------------------------------------------------------------