octave_packages/fpl-1.2.0/pdesurf.m

   1 ##  Copyright (C) 2010  Carlo de Falco
   2 ##
   3 ##  This file is part of:
   4 ##         FPL - Fem PLotting package for octave
   5 ##
   6 ##  FPL is free software; you can redistribute it and/or modify
   7 ##  it under the terms of the GNU General Public License as published by
   8 ##  the Free Software Foundation; either version 2 of the License, or
   9 ##  (at your option) any later version.
  10 ##
  11 ##  FPL is distributed in the hope that it will be useful,
  12 ##  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 ##  GNU General Public License for more details.
  15 ##
  16 ##  You should have received a copy of the GNU General Public License
  17 ##  along with FPL; If not, see <http://www.gnu.org/licenses/>.
  18 ##
  19 ##  author: Carlo de Falco     <cdf _AT_ users.sourceforge.net>
  20
  21 ## -*- texinfo -*-
  22 ## @deftypefn {Function File} {@var{h} =} pdesurf (@var{p}, @var{t}, @var{u})
  23 ##
  24 ## Plot a 3D surface given node or element data on a triangular mesh.
  25 ## @var{p}, @var{t} are the mesh vertices and connectivity, @var{u} node or element data.
  26 ##
  27 ## @seealso{fpl_dx_write_field, fpl_vtk_write_field}
  28 ##
  29 ## @end deftypefn
  30
  31 function h = pdesurf (p, t, u)
  32
  33   ## Check input
  34   if nargin!=3
  35     error("pdesurf: wrong number of input parameters");
  36   endif
  37
  38   nel = columns (t);
  39   npt = columns (p);
  40   if (numel (u) != npt && numel (u) != nel)
  41     error("pdesurf: wrong data size");
  42   endif
  43
  44   hs = ishold ();
  45
  46 ### node data
  47   if (numel (u) == npt)
  48     ## normalize data
  49     c  = colormap;
  50     uc = sum (u(t(1:3, :)), 1)/3;
  51     uc = floor ((rows (c)-1)*(uc - min (uc))/(max (uc) - min (uc))) + 1;
  52     H = patch ('Faces', t(1:3, :)',
  53                'Vertices', [p(1,:)', p(2,:)', u],
  54                'FaceVertexCData', c(uc(:), :));
  55 ### triangle data
  56   elseif (numel (u) == nel)
  57     tri = reshape (1:3*nel, 3, [])';
  58     pt(:, 1) = p(1, t(1:3, :));
  59     pt(:, 2) = p(2, t(1:3, :));
  60     pt(:, 3) = repmat (u(:)', 3, 1)(:);
  61     ## normalize data
  62     c  = colormap;
  63     uc = floor ((rows (c)-1)*(u - min (u))/(max (u) - min (u))) + 1;
  64     H = patch ('Faces', tri,
  65                'Vertices', pt,
  66                'FaceVertexCData', c(uc(:), :),
  67                'FaceColor', 'flat');
  68   endif
  69
  70   if (nargout == 1)
  71     h = H;
  72   endif
  73
  74   if (hs)
  75     hold on;
  76   else
  77     hold off;
  78   endif
  79
  80 endfunction
  81
  82 %!demo
  83 %! msh = msh2m_structured_mesh ([0:.05:1], [0:.1:1], 1, 1:4, 'random');
  84 %! x = msh.p(1,:)'; xm = sum(x(msh.t(1:3,:)),1)/3;
  85 %! y = msh.p(2,:)'; ym = sum(y(msh.t(1:3,:)),1)/3;
  86 %! pdesurf (msh.p, msh.t, x.*(1-x).*y.*(1-y))
  87 %! title ('node data')
  88 %! view(3)
  89 %! figure ()
  90 %! pdesurf (msh.p, msh.t, xm.*(1-xm).*ym.*(1-ym))
  91 %! title ('element data')
  92 %! view(3)