function g=generate_mesh(x,y,z,e);
%
% generate_mesh::generate a finite element mesh based on bathymetry data(step by step)
% function g=generate_mesh(x,y,z,(optional) e );
% This is a (crude) interface to use a bathymetry data to generate a
% traingular 2D finite element grid using triangle by J. Shewchuk.
% The string constants triangle_path and showme_path must be changed 
% to point to your triangle and showme executables.
%
% Inputs:  There are three ways you can input bathymetry data:
%
%         1) As gridded bathymetry data: 
%            >>g=generate_mesh(x,y,z) 
%               where
%               z (MxN) grid of bathymetric depth.
%               x (1xN) x coordinate of columns of z 
%               y (Mx1) y coordinate of rows of z.
%
%         2) As scattererd bathymetry data with a pre-defined triangulation
%           (the triangulation is used for interpolation and contouring).
%            >>g=generate_mesh(x,y,z,e); 
%               where
%               x (Nx1) x coordinate of depth measurements
%               y (Nx1) y coordinate of depth measurements
%               z (Nx1) depth measeurements
%               e (Kx3) verticies numbers for triangles in x and y.
%
%         3) As scattererd bathymetry data with no pre-defined triangulation  
%            >>g=generate_mesh(x,y,z); 
%               where
%               x (Nx1) x coordinate of depth measurements
%               y (Nx1) y coordinate of depth measurements
%               z (Nx1) depth measeurements,
%         In this case delaunay triangulation is used for interpolation. 
%
%         4) As vector of bathymetry data bases  
%            >>g=generate_mesh(batvect);
%            where 
%               batvect(1).z (M1xN1) grid of bathymetric depth.
%               batvect(1).x (1xN1) x coordinate of columns of z 
%               batvect(1).y (M1x1) y coordinate of rows of z.
%               batvect(2).z (M2xN2) grid of bathymetric depth.
%               batvect(2).x (1xN2) x coordinate of columns of z 
%               batvect(2).y (M2x1) y coordinate of rows of z.
% ...
%               batvect(j).z (MjxNj) grid of bathymetric depth.
%               batvect(j).x (1xNj) x coordinate of columns of z 
%               batvect(j).y (Mjx1) y coordinate of rows of z.
% ...
%               batvect(k).x (Nkx1) x coordinate of depth measurements
%               batvect(k).y (Nkx1) y coordinate of depth measurements.
%               batvect(k).z (Nkx1) list of bathymetric depth.
%               batvect(k).e (Kkx3) verticies numbers for triangles in x and y.
%          in case 4 the earlier data bases are take priority over the
%          later ones.  so if you have a high resolution local bathymetry
%          embedded in a course resolution broad coverage database list the
%          course/ broad data base last.
%       
% ...
%Outputs:  g finite element grid structure with fields g.x,g.y,g.z,g.e,g.bnd
%          output is also written in standard file formats to files 
%			Dartmouth FE model grid files for final triangle grid
%				[nml_fe_root,'.nod']
%				[nml_fe_root,'.ele']
%				[nml_fe_root,'.bat']
%			
%			Triangle grid files .node and .ele for each refinement of the grid
%			and the .poly file for the boundary pslg. 
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
% Note to use multiple bathymetry sets:  Send the bathymetry

%triangle [-prq__a__AcevngBPNEIOXzo_YS__iFlsCQVh] input_file
%    -p  Triangulates a Planar Straight Line Graph (.poly file).
%    -r  Refines a previously generated mesh.
%    -q  Quality mesh generation.  A minimum angle may be specified.
%    -a  Applies a maximum triangle area constraint.
%    -A  Applies attributes to identify elements in certain regions.
%    -c  Encloses the convex hull with segments.
%    -e  Generates an edge list.
%    -v  Generates a Voronoi diagram.
%    -n  Generates a list of triangle neighbors.
%    -g  Generates an .off file for Geomview.
%    -B  Suppresses output of boundary information.
%    -P  Suppresses output of .poly file.
%    -N  Suppresses output of .node file.
%    -E  Suppresses output of .ele file.
%    -I  Suppresses mesh iteration numbers.
%    -O  Ignores holes in .poly file.
%    -X  Suppresses use of exact arithmetic.
%    -z  Numbers all items starting from zero (rather than one).
%    -o2 Generates second-order subparametric elements.
%    -Y  Suppresses boundary segment splitting.
%    -S  Specifies maximum number of added Steiner points.
%    -i  Uses incremental method, rather than divide-and-conquer.
%    -F  Uses Fortune's sweepline algorithm, rather than d-and-c.
%    -l  Uses vertical cuts only, rather than alternating cuts.
%    -s  Force segments into mesh by splitting (instead of using CDT).
%    -C  Check consistency of final mesh.
%    -Q  Quiet:  No terminal output except errors.
%    -V  Verbose:  Detailed information on what I'm doing.
%    -h  Help:  Detailed instructions for Triangle.
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
% Keston Smith, Ata Bilgili, May 2001.
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%
% Update by Ata Bilgili, Dartmouth College, 10/29/2003.
%   - Added Collins bandwidth reduction algorithm. Made Cuthill-Mackee default.
%   - Added a quick wireframe plot of the mesh just before bathymetry interpolation for diagnostic
%     mesh visualization before bathymetry interpolation that might take a long time.
%   - Switched from Matlab Delaunay triangulation to Triangle one in the xyz2bat.m routine. The Matlab
%     one had a buggy behaviour with large data set triangulation and was causing crossing triangles.
%   - Added possibility to turn exact arithmetic off (EXARS, EXAR) for scattered bathymetry data set
%     triangulation on Intel platforms.
%   - Added control "if" sentence to the refinement section in the case Triangle does not finish. Now it 
%     goes back to the previous grid instead of just giving an error message and quitting.
%   - Changed string comparison (if bat.battype=="") in the OA input section to strcmpi.
%   - Removed final interpolation of the bathymetry into g.z. Instead a temporary variable gzfinal is 
%     used now to carry over the already interpolated g.z values in the final section.
%   - Similarly to the above, eliminated intermediate interpolations in the refinement section if
%     user wants to "redo" or "go back n steps" during refinement. The bathymetric map is now saved
%     for each iteration in a gtemp(num2str(k)) array. This speeds up the refinement process.
%   - Added default values to interparm structure elements, clevs in contouring and to many other input.
%   - Added keyboard input error handling to keep Matlab from quitting if it does not receive
%     the correct input. This is done for switches used the most. There are many that did not get
%     fixed. This is an ongoing project.
%
%XXXXXXXXXXXXXXXXXXXXXXXXXXXX<Verify inputs are consistent>XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

%triangle_path='/afs/thayer.dartmouth.edu/data/ocean/MSHGN/abilgili/meshwork/triangle/triangle'
%triangle_path='triangle'
%triangle_path='C:\Progra~1\CygWin\usr\local\triangle\triangle.exe'
triangle_path='/afs/isis/depts/marine/opnml/bin/triangle'

disp(['Before you begin, check that the Triangle path is correct and make sure you have the OPNML '...
        ,'finite element toolbox in your Matlab path.']);
disp(['    ']);
disp(['Turning exact arithmetic off may be necessary for Triangle to finish on Intel architecture.'])
disp(['This is not good practice and should only be performed if Triangle does not finish during'])
disp(['bathymetric structure triangulation with scattered bathymetry...']) 
disp(' ')

EXAR=input(['Would you like to turn off exact arithmetic for scattered bathymetric dataset interpolation 1(Yes)/0(No) (Defaults to 0)? : ']);
    if isempty(EXAR);
     EXAR=0;
    end
    
disp(['     ']);
% Reminder: ispc.m/isunix.m routines may be implemented later for the above after more tests. AB

    if nargin ==3,
    bat=xyz2bat(EXAR,triangle_path,x,y,z);
    clear x y z;
    
elseif nargin==1
    bat=x;
    clear x y z;
elseif nargin==4
    bat=xyz2bat(EXAR,triangle_path,x,y,z,e);
    clear x y z e;
end
[pslg,interparm,nml_fe_root,triangle_root,min_depth]=get_battri_parameters(bat);
%
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX<Begin manual editing of PSLG>XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%XXXXX<Edit a PSLG to define domain boundaries and constrain mesh creation process>XXXXXXXXXXXXXXXXX
%XXXXXXXXXXXXXXXXXXXXXXXXX<and compute the first cut triangulation>XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%
donepslgedit=0;
donerefinemesh=0;
k=0;
constraint=[];
while ~donerefinemesh,
    while ~donepslgedit,
        subplot(1,1,1);
        pslg=editpslg(pslg,bat);
        [g,donepslgedit]=firstcutmesh(pslg,bat,triangle_root,triangle_path,min_depth,interparm);
    end
%XXXXXXXXXXXXXXXXXXXXXXXX<Finished manual editing of PSLG>XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%XXXXXXXXXXXXXXXXXXXXXXXXXX<Begin mesh refinement loop>XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    k=1;
    donerefinemesh=0;
    [g,k,donerefinemesh,gzfinal,constraint]=refinemeshscript(g,k,donerefinemesh,triangle_path,triangle_root,bat,min_depth,interparm,constraint);
    donepslgedit=0;
end

%XXXXXXXXXXXXXXXXXXXXXXX<Finished with mesh refinement loop>XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%XXXXXXXXXXXXXXXXXX<Output NML standard finite element mesh files>XXXXXXXXXXXXXXXXXXXXXXXXXX

g=triangle2fe_struct2([triangle_root,'.',int2str(k)]);

%
% Reestablish final g.z without interpolating again. AB
%
g.z=gzfinal;

% Reduce the bandwidth of the mesh.
%
disp(' ')
bwmethod=input(['Choose bandwidth reduction algorithm: 0 for Cuthill-McKee (Fast, less reduction, default) or 1 for Collins (Much slower, more reduction): ']);
if isempty(bwmethod)==1
    bwmethod=0;
end
g=reduce_bwidth(g,bwmethod);
%
% Write nml standard mesh files
%
write_nod(g.x,g.y,[nml_fe_root,'.nod']);	
write_bat(g.z,[nml_fe_root,'.bat']);	
write_ele(g.e,[nml_fe_root,'.ele']);
%
% Yeehaa!...
%