clear all
clc

diam_bottom = 50.0;
diam_top = 83.0;
wall_thickness = 3.0;
bottom_thickness = 6.0;
top_level = 105.0;
fluid_level = 66.0;
floor_level = 0.0;

inner_diam = @(height) diam_bottom - 2*wall_thickness + (diam_top - diam_bottom)/(top_level - bottom_thickness)*(height - bottom_thickness);

diam_fluid = inner_diam(fluid_level);


g=[2 -diam_bottom/2 diam_bottom/2  floor_level floor_level  1  0;        % #vertices,v_1x, v_2x, v_1y, v_2y, subdomain_left, subdomain_right
   2 diam_bottom/2 diam_top/2  floor_level top_level  1  0;
   2 diam_top/2 diam_top/2-wall_thickness  top_level top_level  1  0;
   2 diam_top/2-wall_thickness diam_fluid/2  top_level fluid_level  1  3;
   2 diam_fluid/2 diam_bottom/2-wall_thickness  fluid_level bottom_thickness 1 2;
   2 diam_bottom/2-wall_thickness -(diam_bottom/2-wall_thickness) bottom_thickness bottom_thickness 1 2;
   2 -(diam_bottom/2-wall_thickness) -diam_fluid/2 bottom_thickness fluid_level 1 2;
   2 -diam_fluid/2 -(diam_top/2-wall_thickness) fluid_level top_level 1 3;
   2 -(diam_top/2-wall_thickness) -diam_top/2 top_level top_level 1 0;
   2 -diam_top/2 -diam_bottom/2 top_level floor_level 1 0;

   % surface between fluid and air
   2 -diam_fluid/2 diam_fluid/2 fluid_level fluid_level 3 2;

   % top (air) surface
   2 -(diam_top/2-wall_thickness) diam_top/2-wall_thickness top_level top_level 0 3
   ]';


[p,e,t] = initmesh(g,'hmax',10); 
pdemesh(p,e,t)


ascii_write_mesh( p, t, e, mfilename);

ascii_write_subdomains( p, t, e, mfilename);