function mesh_oops_test ( )

%*****************************************************************************80
%
%% mesh_oops_test() demonstrates MESH2D on some example problems.
%
%  Modified:
%
%    27 October 2011
%
%  Author:
%
%    Darren Engwirda
%
  clc
%
%  I have wasted considerable time reading inadequate MATLAB diagnostics
%  and unhelpful QHULL documentation, just to try to suppress warnings
%  about TSEARCH being retired (and I can't figure out how to use MATLAB's
%  suggested DelaunayTri/pointLocation replacement properly!) and the
%  warnings about duplicate points being fed into DelaunayN, which are
%  NOT helped by trying the UNIQUE function.  So...TURN THE WARNINGS OFF,
%  THEY JUST MAKE ME ANGRY.  (JVB, 27 October 2011.)
%
  warning ( 'off' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'mesh_oops_test:\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Demonstrate mesh2d() for creating meshes.\n' );

%*****************************************************************************80
%  #1) Circle, with 24 points on the circumference.
%*****************************************************************************80
  fprintf ( 1, '\n' );
  fprintf ( 1, '  EXAMPLE 1: The Circle\n' );
  fprintf ( 1, '  We start with 24 points on the circumference.\n' );
  fprintf ( 1, '  The element size function is generated automatically to try\n' );
  fprintf ( 1, '  to adequately resolve the geometry.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This means that the mesh size is related to the length of \n' );
  fprintf ( 1, '  the line segments used to define the geometry.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );

  dtheta = pi / 12;
  theta  = (-pi:dtheta:(pi-dtheta))';
  node   = [cos(theta) sin(theta)];

  [p,t] = mesh2d ( node );

  [ nv, ~ ] = size ( node );
  [ np, ~ ] = size ( p );
  [ nt, ~ ] = size ( t );
  fprintf ( 1, ...
    '  %d boundary vertices input, %d nodes and %d triangles created\n', ...
    nv, np, nt );

  fprintf ( 1, 'The following example is the same as the last, \n' );
  fprintf ( 1, 'but with more lines used to represent the circle.\n' );
  pause ( 5 )

%*****************************************************************************80
%  #2) Circle, with 150 points on the circumference.
%*****************************************************************************80
  dtheta = pi/75;
  theta  = (-pi:dtheta:(pi-dtheta))';
  node   = [cos(theta) sin(theta)];

  plot ( node(:,1), node(:,2), 'r.', 'MarkerSize', 20 );
  axis ( 'equal' );
  grid ( 'on' );
  title ( '#1: 150 Boundary vertices' );

  [p,t] = mesh2d(node);

  [ nv, ~ ] = size ( node );
  [ np, ~ ] = size ( p );
  [ nt, ~ ] = size ( t );
  fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );

%*****************************************************************************80
%                           Square
%*****************************************************************************80

fprintf ( 1, 'It is often necessary to specify the element size in some locations.\n' );                                               ...
pause ( 5 )

close all

node = [
    0   0
    1   0
    1   1
    0   1
        ];

hdata.fun = @hfun1;

[p,t] = mesh2d(node,[],hdata);

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );
%*****************************************************************************80
%                       Sliver regions
%*****************************************************************************80

fprintf ( 1, 'Mesh2d can now deal with very fine "sliver" geometry features.\n' );

pause ( 5 )

close all



        node  = [0 0
                 3 0
                 3 3
                 0 3
                 0.1 1
                 0.11 1
                 0.11 2
                 0.1 2
                 ];
        cnect = [1 2
                 2 3
                 3 4
                 4 1
                 5 6
                 6 7
                 7 8
                 8 5
                 ];

        [p,t] = mesh2d(node,cnect);

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );
%*****************************************************************************80
%                    Cylinder in crossflow
%*****************************************************************************80

fprintf ( 1, 'The following is a mesh used to simulate the flow past a cylinder. \n' );

pause ( 5 )

close all

theta = 0:pi/100:(2*pi-pi/100);
x     = cos(theta)/2;
y     = sin(theta)/2;

node = [ x'  y'
        -5  -10
        25  -10
        25   10
        -5   10 
       ];

n = size(node,1)-4;

cnect = [(1:n-1)' (2:n)'
          n        1
          n+1      n+2
          n+2      n+3
          n+3      n+4
          n+4      n+1
        ]; 

hdata = [];
hdata.fun = @hfun2;

[p,t] = mesh2d(node,cnect,hdata);

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );
%*****************************************************************************80
%                         Airfoil + flap
%*****************************************************************************80

fprintf ( 1, 'The following is a mesh used to simulate the flow past an airfoil/flap configuration. \n' );

pause ( 5 )

close all

wing = [
 1.00003  0.00126  
 0.99730  0.00170  
 0.98914  0.00302  
 0.97563  0.00518  
 0.95693  0.00812  
 0.93324  0.01176  
 0.90482  0.01602  
 0.87197  0.02079  
 0.83506  0.02597  
 0.79449  0.03145  
 0.75070  0.03712  
 0.70417  0.04285  
 0.65541  0.04854  
 0.60496  0.05405  
 0.55335  0.05924  
 0.50117  0.06397  
 0.44897  0.06811  
 0.39733  0.07150  
 0.34681  0.07402  
 0.29796  0.07554  
 0.25131  0.07597  
 0.20738  0.07524  
 0.16604  0.07320  
 0.12732  0.06915  
 0.09230  0.06265  
 0.06203  0.05382  
 0.03730  0.04324  
 0.01865  0.03176  
 0.00628  0.02030  
 0.00015  0.00956  
 0.00000  0.00000  
 0.00533 -0.00792  
 0.01557 -0.01401  
 0.03029 -0.01870  
 0.04915 -0.02248  
 0.07195 -0.02586  
 0.09868 -0.02922  
 0.12954 -0.03282  
 0.16483 -0.03660  
 0.20483 -0.04016  
 0.24869 -0.04283  
 0.29531 -0.04446  
 0.34418 -0.04510  
 0.39476 -0.04482  
 0.44650 -0.04371  
 0.49883 -0.04188  
 0.55117 -0.03945  
 0.60296 -0.03655  
 0.65360 -0.03327  
 0.70257 -0.02975  
 0.74930 -0.02607  
 0.79330 -0.02235  
 0.83407 -0.01866  
 0.87118 -0.01512  
 0.90420 -0.01180  
 0.93279 -0.00880  
 0.95661 -0.00621  
 0.97543 -0.00410  
 0.98901 -0.00254  
 0.99722 -0.00158  
 0.99997 -0.00126   
 ];

flap = rotate(0.4*wing,10);
flap = move(flap,0.95,-0.1);
wing = rotate(wing,5);
wing = move(wing,0,0.05);

wall = [
      -1      -3
       4      -3
       4       3
      -1       3 
       ];

nwing = size(wing,1);
nflap = size(flap,1);
nwall = size(wall,1);

cwing = [(1:nwing-1)', (2:nwing)'; nwing, 1];
cflap = [(1:nflap-1)', (2:nflap)'; nflap, 1];
cwall = [(1:nwall-1)', (2:nwall)'; nwall, 1];

cnect = [
        cwing
        cflap+nwing
        cwall+nflap+nwing
        ];

node = [wing; flap; wall];

hdata = [];
hdata.edgeh = [(1:size(cnect,1)-4)', 0.005*ones(size(cnect,1)-4,1)];

options.dhmax = 0.25;
options.stats = false;

[p,t] = mesh2d(node,cnect,hdata,options);

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );
%*****************************************************************************80
%                        Lake Superior
%*****************************************************************************80

fprintf ( 1, 'The following is a mesh of Lake Superior and is a standard test of mesh algorithms. \n');

pause ( 5 )

close all

%       From:
%       Computer Solutions Europe AB
%       Bjornnasvagen 21
%       S-113 47 Stockholm
%       SWEDEN
%       Tel:    +46 8 15 30 22
%       Fax:    +46 8 15 76 35
%       E-mail: info@comsol.se
%       URL:    http://www.comsol.se/

p1=[
-8.9154147   1.6615920	
-8.8847064   1.5538653
-8.7270571   1.4396306
-8.6780329   1.3310253
-8.2763869   1.3642085
-7.7589437   1.5521208
-7.3399811   1.5866412
-7.2064877   1.7400111
-6.9871610   1.7838988
-6.9190874   1.6751018
-6.6594362   1.8237880
-6.4705078   2.0285109
-6.3982459   2.0257761
-6.1776033   2.1236043
-5.9683636   1.9041636
-6.0500631   1.6420917
-6.2751037   1.4382319
-6.2083449   1.2768143
-6.3593873   1.1233594
-6.4379782   0.9672913
-6.4389687   0.9408290
-6.3308569   0.9103191
-5.9216680   1.1606956
-5.8853530   1.1594465
-5.5306989   0.8827500
-5.4571118   0.9068911
-5.2776995   0.8218321
-4.6882658   1.0696109
-4.4325551   1.1157840
-3.9867787   1.5284986
-3.6232783   1.5733560
-3.5156848   1.5181578
-3.1151141   1.6162773
-2.6031886   1.9253444
-2.4201615   2.0814581
-2.0576784   2.1825510
-1.9469493   2.3929843
-1.6564457   2.6016249
-1.4385855   2.8113568
-1.1857180   2.9681802
-0.8965760   3.2839699
-0.4656316   3.4939341
-0.0715676   3.6520868
 0.3221563   3.5994536
 0.6089023   3.4943945
 0.6092879   3.3885375
 0.3942456   3.3878929
 0.2866788   3.4141384
 0.0717038   3.3345153
 0.0896864   3.2286668
 0.1076576   3.1757488
-0.1436793   3.0169900
-0.3235837   2.8584538
-0.4676923   2.7529346
-0.9391942   2.0668111
-1.0495193   1.7498742
-1.0882363   1.3531376
-1.0158452   1.3262275
-0.9424102   1.4845979
-0.6152380   1.7477533
-0.5061946   1.9061816
-0.2890742   2.0115316
-0.5066670   1.7473959
-0.6160022   1.5360392
-0.3979720   1.8000363
-0.1808403   1.8525961
 0.2894343   1.7998173
 0.5431103   1.6681146
 0.9058832   1.5373304
 0.9432856   1.3258124
 1.0534256   1.1147328
 1.1641909   0.8507940
 1.3109203   0.6930545
 1.3117251   0.5871977
 1.4214690   0.5351345
 1.9317215   0.5400872
 2.2226054   0.5965222
 2.3340522   0.4391798
 2.4073552   0.4137156
 2.4449497   0.3348350
 2.5909063   0.3369578
 2.6249835   0.4962917
 2.7343488   0.4979874
 2.8454551   0.3938967
 2.9366470   0.3954205
 2.9201944   0.2892562
 2.9931921   0.2905022
 3.3135431   0.7199108
 3.3863589   0.7213257
 3.9960329   1.1315466
 4.0705567   1.0538628
 4.6490632   1.1744760
 4.9039980   1.1550887
 5.4495532   1.1450567
 5.9850466   1.4278414
 6.1653805   1.4607402
 6.3476152   1.4409063
 6.6355959   1.5048353
 6.5309084   1.3948027
 6.4623697   1.2862027
 6.4803085   0.8098819
 6.4094979   0.7542354
 6.4100090   0.7012564
 6.4882671   0.5981826
 6.7068647   0.6065329
 6.8546905   0.5593290
 6.9296957   0.5093111
 7.3625916   0.6332439
 7.4195041   0.5826584
 7.4810256   0.4262426
 7.6622000   0.4606821
 7.7339418   0.4903887
 7.6927729   0.5946095
 7.6199031   0.5913803
 7.5788199   0.6956243
 7.4649955   0.7967085
 7.4604098   0.9025555
 7.6668095   1.1767640
 7.7709460   1.2875089
 7.7298813   1.3917132
 7.6707366   1.4951115
 7.6356734   1.4670289
 7.6246141   1.3074531
 7.5532118   1.2777792
 7.5280861   1.4357493
 7.5345391   1.7011603
 7.7131429   1.7621690
 7.8603274   1.7158544
 7.9228745   1.9308907
 7.7348899   2.0813402
 7.3713821   2.1182087
 7.2698318   1.9547751
 7.1885795   2.1633926
 7.1152657   2.1868122
 7.0441691   2.1573400
 6.9336956   2.2058041
 6.9676141   2.2602054
 7.0648074   2.5292990
 7.0897106   2.7954185
 7.2917362   3.1221593
 7.3230152   3.2295634
 7.2074219   3.4101256
 7.0987913   3.4319941
 7.0562590   3.5892427
 6.8349672   3.7390777
 6.5426013   3.8863685
 6.3538119   4.1440578
 6.6493349   4.7917061
 6.5385166   4.8932842
 6.7052993   5.1650313
 6.7009764   5.2708721
 5.7441174   5.2345817
 5.5686945   5.1755508
 4.7886070   5.1509481
 4.4636817   5.3536511
 4.2075761   5.6646354
 3.6075780   7.1858876
 3.3873644   7.6576363
 3.2804384   7.7613007
 3.1762847   7.7327191
 3.0381962   7.6506117
 2.8286881   7.6467634
 2.8598071   7.8590921
 2.6845531   7.9089847
 2.5468001   7.8008409
 2.4437671   7.6933492
 2.3390419   7.6917801
 2.1983122   7.7691861
 2.1616242   7.9010276
 2.0918985   7.9000964
 2.0235167   7.7933405
 1.9540682   7.7660064
 1.5705080   7.7353016
 1.5002116   7.7875578
 1.2547448   7.9442315
 0.7309650   8.1526824
 0.4522013   8.2575018
 0.1738657   8.3098811
-0.0694533   8.5215153
-0.1737495   8.4157382
-0.2084645   8.4422444
-0.3476150   8.3630957
-0.4518238   8.3898231
-0.5558115   8.4695476
-0.5215963   8.3106439
-0.5221188   8.1518584
-0.4879606   7.9400337
-0.2092654   7.8335660
-0.2097517   7.4630660
-0.5257651   7.0403597
-0.5604465   7.1463352
-0.7716285   6.9354925
-0.9473082   6.8835005
-0.9833619   6.7249256
-1.1593398   6.6731632
-1.2287946   6.7795396
-1.1574372   6.9378046
-0.6994015   7.4115277
-0.5933112   7.7286747
-0.6973173   7.8878837
-0.9759166   7.9422776
-1.1175557   7.6256239
-1.1194069   7.3609825
-1.2951661   7.2564438
-1.5483098   6.4117610
-1.9366363   6.3100700
-1.8637764   6.5209552
-1.7582839   6.5197737
-1.7202965   6.7840348
-1.5427314   6.9939696
-1.5417155   7.0998256
-1.5762336   7.1530937
-1.9985176   6.9989869
-2.3506579   6.8978706
-2.5645144   6.6893724
-2.5351899   6.3183276
-2.7562921   5.7395724
-3.7980607   5.0191439
-4.3373509   4.7676755
-4.9834948   4.5740477
-5.9256532   4.1282314
-7.0382285   3.1643631
-7.4517045   2.7580044
-8.2033903   2.1245258
-8.7235361   1.8639192
];

p2=[
-6.0222359   1.4026510
-5.9152018   1.5578835
-5.8726343   1.5299217
-5.7920161   1.5536639
-5.8345898   1.5816055
-5.8318773   1.6609928
-5.6635570   1.7083083
-5.5966615   1.7590827
-5.6808641   1.8413432
-5.7929717   1.7391469
-6.0818599   1.5637339
-6.0865686   1.4314222
];

p3=[
-5.7324455   1.9225485
-5.5683675   1.9594866
-5.5335991   1.9159672
-5.5318777   1.9688923
-5.4498160   2.0457126
-5.4914569   2.1000333
-5.7297704   2.0019357
-5.7677013   1.9502324
];

p4=[
-5.2706487   2.2254789
-5.2257103   2.2770640
-5.2196558   2.3563407
-5.1911650   2.4614151
-5.2924029   2.4539566
-5.3166868   2.3699543
-5.3193455   2.2852741
-5.3444930   2.1748120
];

p5=[
-2.6285090   4.8904087
-2.5231726   4.8093509
-2.4524965   4.7818143
-2.1664250   4.9101487
-1.8814086   4.9860466
-2.0588911   4.9881727
-2.0937130   5.0415487
-1.9150532   5.1452403
-1.7011793   5.2487523
-1.5239785   5.2470098
-1.1676824   5.5087407
-1.2022868   5.6148538
-0.9185046   5.7718016
-0.8828910   5.8245348
-0.9532116   5.8778612
-1.3075201   5.7215257
-1.5563954   5.5649116
-1.6276655   5.5126727
-1.7340945   5.4872998
-1.8405570   5.4619959
-2.0900008   5.3326514
-2.3756722   5.1777288
-2.5904890   5.0486276
];

p6=[
-0.1745905   7.6482740
-0.1048759   7.4629629
-0.0174793   7.4629295
-0.0349355   7.5687895
 0.0349007   7.7275752
 0.2267794   7.7806613
 0.3839078   7.7280333
 0.5583180   7.7550129
 0.6632222   7.7024853
 0.8201638   7.7296796
 0.7325480   7.8086476
 0.6273773   7.9405222
 0.4880418   7.9135694
 0.1393711   7.9922753
 0.0348428   7.9922181
-0.1743299   7.8864526
];

p7=[
 4.5964504   4.4039327
 4.8814214   4.4120893
 4.9866866   4.4681990
 5.1274771   4.5254771
 5.1614199   4.5795088
 5.0886198   4.6302255
 4.9818982   4.6269721
 4.8364965   4.7285905
 4.5535316   4.6675729
 4.3785074   4.5569131
 4.3813116   4.4510625
];


n1 = size(p1,1);
n2 = size(p2,1);
n3 = size(p3,1);
n4 = size(p4,1);
n5 = size(p5,1);
n6 = size(p6,1);
n7 = size(p7,1);

c1 = [(1:n1-1)', (2:n1)'; n1, 1];
c2 = [(1:n2-1)', (2:n2)'; n2, 1];
c3 = [(1:n3-1)', (2:n3)'; n3, 1];
c4 = [(1:n4-1)', (2:n4)'; n4, 1];
c5 = [(1:n5-1)', (2:n5)'; n5, 1];
c6 = [(1:n6-1)', (2:n6)'; n6, 1];
c7 = [(1:n7-1)', (2:n7)'; n7, 1];

node  = [p1; p2; p3; p4; p5; p6; p7];
cnect = [
        c1
        c2+n1
        c3+n2+n1
        c4+n3+n2+n1
        c5+n4+n3+n2+n1
        c6+n5+n4+n3+n2+n1
        c7+n6+n5+n4+n3+n2+n1
        ];

% Make mesh
[p,t] = mesh2d(node,cnect);

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );


fprintf ( 1, 'The following shows the influence of gradient limiting on the size function. \n' );

pause ( 5 )

close all

%
%  FAILURE: call DSEARCH in TINTERP.
%
% Solver options
%options.dhmax = 0.1;

% Make mesh
%[p,t] = mesh2d(node,cnect,[],options);

%[ nv, ~ ] = size ( node );
%[ np, ~ ] = size ( p );
%[ nt, ~ ] = size ( t );
%fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );

%*****************************************************************************80
% User defined size functions
%*****************************************************************************80

fprintf ( 1, 'The following example shows how the element size can be controlled.\n' );

pause ( 5 )

close all

hdata = [];
options.output = false;
options.stats = false;

node = [0 0; 1 0; 1 1; 0 1];              % Simple square example
% 
% [p,t] = mesh2d(node,[],[],options);       % Auto size fun only
% 
% figure
% subplot(2,2,1)
% patch('faces',t,'vertices',p,'facecolor','none','edgecolor','b');
% axis equal off;
% title('Automatic size fun only')
% 
% 
% hdata.hmax = 0.1;                         % Global size of 0.1
% 
% [p,t] = mesh2d(node,[],hdata,options);
% 
% subplot(2,2,3)
% patch('faces',t,'vertices',p,'facecolor','none','edgecolor','b');
% axis equal off;
% title('Global size fun with h0=0.1')


hdata.hmax = 0.1;                         % Global size of 0.1
hdata.edgeh = [1,0.05];                   % Boundary layer on bottom edge

[p,t] = mesh2d(node,[],hdata,options);

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );

figure
subplot(2,2,2)
patch('faces',t,'vertices',p,'facecolor','none','edgecolor','b');
axis equal off;
title('Additional boundary layer function on bottom edge')

% hdata = [];
% hdata.fun = @hfun1;
% 
% [p,t] = mesh2d(node,[],hdata,options);
% 
% subplot(2,2,4)
% patch('faces',t,'vertices',p,'facecolor','none','edgecolor','b');
% axis equal off; hold on;
% title('User defined funcion centred at [0.25,0.75]')

%*****************************************************************************80
% Refinement
%*****************************************************************************80

fprintf ( 1, 'The following example shows how an existing mesh can be refined.\n' );

pause ( 5 )

close all

hdata = [];
options.output = false;
options.stats = false;
hdata.hmax = 0.1;

node = [0 0; 1 0; 1 1; 0 1];              % Simple square example

[p,t] = mesh2d(node,[],hdata,options);    % Auto size fun only

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );

figure
subplot(2,2,1)
patch('faces',t,'vertices',p,'facecolor','none','edgecolor','b');
axis equal off; hold on;
title('Uniform mesh from MESH2D')

node = [0.3,0.3; 0.7,0.3; 0.7,0.7; 0.3,0.7];

in = inpoly(p,node);
ti = sum(in(t),2)>0;
[p,t] = refine(p,t,ti);

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );

subplot(2,2,2)
patch('faces',t,'vertices',p,'facecolor','none','edgecolor','b');
axis equal off; hold on;
title('Mesh refined in centre region using REFINE')

[p,t] = smoothmesh(p,t);

[ nv, ~ ] = size ( node );
[ np, ~ ] = size ( p );
[ nt, ~ ] = size ( t );
fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );

subplot(2,2,3)
patch('faces',t,'vertices',p,'facecolor','none','edgecolor','b');
axis equal off; hold on;
title('Smoothed mesh using SMOOTHMESH')
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'mesh_oops_test:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function p = move(p,xm,ym)

%*****************************************************************************80
%
% Move a node set p by [xm,ym]

n = size(p,1);
p = p + [xm*ones(n,1), ym*ones(n,1)];

end      % move()

function p = rotate(p,A)

%*****************************************************************************80
%
% Rotate a node set p by A degrees.

A = A*pi/180;
T = [ cos(A), sin(A)
     -sin(A), cos(A)];
p = (T*p')';

end      % rotate()
function h = hfun1(x,y)

%*****************************************************************************80
%
% User defined size function for square

h = 0.01 + 0.1*sqrt( (x-0.25).^2+(y-0.75).^2 );

end      % hfun1()

function h = hfun2(x,y)

%*****************************************************************************80
%
% User defined size function for cylinder

h1 = inf*ones(size(x));
in = (x>=0)&(x<=25)&(y>=-3)&(y<=3);
h1(in) = 0.2;

r = sqrt(x.^2+y.^2);
h2 = inf*h1;
h2(r<=3) = 0.02 + 0.05*r(r<=3);

h = min(h1,h2);

end      % hfun2()
