femwind/add_terrain_to_mesh.m

   1 function XX=add_terrain_to_mesh(X, kind, how, val)
   2 % XX=add_terrain_to_mesh(X, kind, how, val)
   3 % Modify 3D mesh vertically to follow terrain
   4 % in:
   5 %      kind  'hill'  add hill
   6 %            numeric uniform shift up
   7 %      how   'shift' move up each column the same
   8 %            'squash' keep top flat
   9 %      val   relative height of the hill
  10 % out:
  11 %      XX    modified mesh
  12
  13 check_mesh(X);
  14
  15 x = X{1}(:,:,1); y=X{2}(:,:,1);z=X{3};
  16 kmax=size(X{1},3);
  17 if ischar(kind),
  18     switch kind
  19         case 'hill'
  20             cx = mean(x(:));
  21             cy = mean(y(:));
  22             hz = max(z(:))-min(z(:));
  23             rx = mean(abs((x(:)-cx)));
  24             ry = mean(abs((y(:)-cy)));
  25             a = ((x-cx)./rx).^2 + ((y-cy)./ry).^2 ;
  26             t = hz*exp(-a*2)*val;
  27         case 'xhill'
  28             cx = mean(x(:));
  29             cy = mean(y(:));
  30             hz = max(z(:))-min(z(:));
  31             rx = mean(abs((x(:)-cx)));
  32             a = ((x-cx)./rx).^2;
  33             t = hz*exp(-a*2)*val;
  34         otherwise
  35             error('add_terrain_to_mesh: unknown kind')
  36     end
  37 elseif isnumeric(kind),
  38     t=kind;
  39 else
  40     error('kind must be string or numeric')
  41 end
  42 switch how
  43     case {'shift','s'}
  44         disp('shifting mesh by terrain vertically')
  45         XX=X;
  46         for k=1:kmax
  47             XX{3}(:,:,k)=X{3}(:,:,k)+t;
  48         end
  49     case {'compress','c','squash'}
  50         if any(any(t > X{3}(:,:,end))),
  51             error('shift values are too large to be compressed')
  52         end
  53         disp('compressing mesh keeping top unchanged')
  54         XX=X;
  55         for k=1:kmax
  56             XX{3}(:,:,k)=X{3}(:,:,k)+t*(XX{3}(1,1,kmax) - XX{3}(1,1,k)) / XX{3}(1,1,kmax);
  57         end
  58     otherwise
  59         error('add_terrain_to_mesh: unknown how')
  60 end
  61 check_mesh(XX)
  62 end