Unstructured grid generation

Outlines the procedure to create the necessary unstructured grid files for FVCOM, including open boundaries.

SMS

To avoid repetition, the manual has an extensive description of how to create an unstructured grid, interpolate bathymetry and export the relevant files in section 20.2.

N.B. In SMS version 10, to export the 'pts' file, select the interpolated bathymetry in the Mesh Data layer, right click on the actual interpolated bathymetry, select Export, change the file type to Generic ASCII file, optionally choose a file name (I suggest casename.dat as the output file name) and click Save. The format of this exported ASCII file can be read in with the MATLAB fvcom-toolbox function read_sms_mesh, and in turn exported to FVCOM format with write_FVCOM_bath. The 2dm file is automatically created once an unstructured grid has been created and the project saved.

Some tips we've discovered follow:

1. To have some control on the resolution within the domain, create arcs with specific node spacings, forcing the mesh to refine through those nodes.
2. The elements which have two nodes on the open boundaries should have their vertices be perpendicular to the boundary. This is quite difficult to achieve automatically, though a MATLAB script to take the positions of the open boundaries, calculate the orientation of a given pair of element edges surrounding the current node and calculate a new point within the model domain perpendicular to that orientation and spaced equal to the two surrounding element lengths should be relatively straightforward to create.
3. Poor quality elements (i.e. those which fail the quality criteria) can often be fixed automatically by selecting Elements > Relax. If no elements have been selected, the entire grid will be relaxed (which may cause previously good elements to fail the quality criteria in some way), otherwise, only the selected elements will be relaxed. Multiple executions of Elements > Relax may be necessary to improve the element quality.

Importing Shapefiles

1. Select the GIS module in SMS.
2. From the Data menu, select Add Shapefile Data. Select a shapefile and click Open.
3. To convert the shapefile to features from which a mesh can be generated, click on the select tool in the GIS module, and select the polygons you wish to use.
4. In the Mapping menu, select Shapes -> Feature Objects to convert the polygons/lines to features which can be edited as per usual in the Map module.

MATLAB fvcom-toolbox

To create the ASCII files necessary for FVCOM using the fvcom-toolbox, the procedure is as follows:

1. Read in an SMS mesh and bathymetry with read_sms_mesh.
2. Extract open boundary nodes (add_obc_nodes_list), assigning a value of 2 to the ObcType in add_obc_nodes_list's arguments (this is "clamped" after Beardsley and Haidvogel (1981)). See Table 6.1 in the FVCOM manual for the ObcType options.
3. Export the FVCOM grid ASCII file (_grd.dat) with write_FVCOM_grid.
4. Export the FVCOM bathymetry ASCII file (_dep.dat) with write_FVCOM_bath.
5. Export the FVCOM open boundary ASCII file (_obc.dat) with write_FVCOM_obc.

Fine tuning an SMS grid with the fvcom-toolbox

Page 77, Figure 6.2 in the FVCOM manual suggests that elements along an open boundary should be as close to right-angled triangles as possible. As such, fvcom-toolbox includes a routine called fix_inside_boundary which will read a given SMS unstructured grid and adjust the positions of the nodes adjacent to the open boundaries specified to ensure they are as close to right-angled as possible. The resulting grid has no additional quality checks performed, so it is often necessary to load the grid into SMS to ensure the elements all fulfil the quality criteria defined in the FVCOM manual. In order to facilitate this, a new fuction write_SMS_2dm has been created to export an unstructured grid to SMS format. The figure below shows the results for the optimisation of an SMS unstructured grid (in red) and the resulting adjusted grid (in blue).

/images/adjusted_triangulation.png

The fix_inside_boundary.m script is a bit brute force and has difficulties with the first and last nodes on an open boundary. This will possibly create overlapping elements at corners and between land and the open boundary. To mitigate this, make your open boundary more curved or edit fix_inside_boundary.m to check the angle between the two adjacent element edges to the current node for some minimum (i.e. if the angle is less than 90 degrees, skip this node).

Finally, the readjustment needed to correct overlapping elements can cause problems with the element table in the 2dm file. Thus, selecting a nodestring and running Nodestrings > Renumber in SMS is sometimes necessary. The bathymetry can then be interpolated and exported as necessary.

A sample MATLAB script showing how to read in SMS files and perform the optimisation is given below.

    % Read in an existing SMS mesh and adjust the nodes adjacent to the 
    % open boundary to be approximately normal to it.
    
    matlabrc
    close all
    clc
    
    global ftbverbose
    ftbverbose = 1; % be noisy
    
    %%%------------------------------------------------------------------------
    %%%                          INPUT CONFIGURATION
    %%%
    %%% Define the input parameters here.
    %%%
    %%%------------------------------------------------------------------------
    
    conf.base = '/data/medusa/pica/models/FVCOM/wavehub/run/';
    
    % Case name for the model inputs and outputs
    conf.casename = 'wavehub_v1';
    
    conf.coordType = 'cartesian'; % 'cartesian' or 'spherical'
    
    % Give some names to the boundaries. This must match the number of node
    % strings defined in SMS. Ideally, the order of the names should match the
    % order in which the boundaries were made in SMS.
    conf.boundaryNames = {'West'}; % Number is important!
    
    % Type of open boundary treatment (see Table 6.1 in the manual and mod_obcs.F).
    % 1 or 2 - Active (ASL): sea level is specified at the open boundary.
    % 3 or 4 - Clamped: zeta = 0 at the boundary (Beardsley and Haidvogel, 1981).
    % 5 or 6 - Implicit Gravity Wave Radiation.
    % 7 or 8 - Partial Clamped Gravity Wave Radiation (Blumberg and Kantha, 1985).
    % 9 or 10 - Explicit Orlanski Radiation (Orlanski, 1976; Chapman, 1985)
    conf.obc_type = 1;
    
    % Output file prefix.
    conf.outbase = fullfile(conf.base, 'input', conf.casename);
    
    
    %%%------------------------------------------------------------------------
    %%%                      END OF INPUT CONFIGURATION
    %%%------------------------------------------------------------------------
    
    %% Prepare the data
    
    % Read the input mesh and bathymetry. Also creates the data necessary for
    % the Coriolis correction in FVCOM.
    Mobj = read_sms_mesh(...
        '2dm', fullfile(conf.base, 'raw_data/', [conf.casename, '.2dm']),...
        'bath', fullfile(conf.base, 'raw_data/',[conf.casename, '.dat']),...
        'coordinate', conf.coordType, 'addCoriolis', true);
    
    % Parse the open boundary nodes and adjust accordingly.
    x3 = Mobj.x;
    y3 = Mobj.y;
    if Mobj.have_strings
        for i = 1:size(Mobj.read_obc_nodes,2)
            nodeList = double(cell2mat(Mobj.read_obc_nodes(i)));
            Mobj = add_obc_nodes_list(Mobj, nodeList, conf.boundaryNames{i}, conf.obc_type);
            % For each of the open boundaries, find the points inside the
            % boundary and move the closest existing point to that location.
            % This should make the boundaries more orthogonal (as suggested by
            % the FVCOM manual in Figure 6.2 on page 77 [as of 2013-03-11]).
            [x3, y3] = fix_inside_boundary(x3, y3, Mobj.read_obc_nodes{i});
        end
    end
    
    %% Write out the required files.
    
    % Make the output directory if it doesn't exist.
    if exist(conf.outbase, 'dir') ~= 7
        mkdir(conf.outbase)
    end
    
    % Save a new 2dm file from the adjusted points for final touching up in
    % SMS.
    write_SMS_2dm(fullfile(conf.outbase, [conf.casename, '_adjusted.2dm']), ...
        Mobj.tri, x3, y3, Mobj.read_obc_nodes)

Other resources for mesh generation

There are other tools for generating unstructured grid meshes that work with FVCOM. At PML we have used Matlab and Python software such as Distmesh, ADMESH, MESH2D, and OceanMesh2D. Other approaches worth noting are Shingle and Triangle