# clone the code
$ git clone git@gitlab.com:oceanbox/meshrefine.git
$ cd meshrefine

# this is only needed if you don't have Rust installed
$ nix-env -iA nixos.cargo nixos.rustc

# finally, we compile the code
$ cargo build --release

This compiles a binary called meshrefine into $PWD/target/release/.

How to run it

The binary has the following options:

USAGE:
    meshrefine [OPTIONS]

OPTIONS:
    -a, --min-angle <MIN_ANGLE>
    -b, --boundary <BOUNDARY>
    -d, --delta-t <DELTA_T>                              [default: 0.1]
    -f, --fscale <FSCALE>                                [default: 1]
    -h, --help                                           Print help information
    -i, --islands <ISLANDS>
        --ld <LD>                                        [default: 10000]
        --lr <LR>                                        [default: 4000]
    -m, --num-repair-steps <NUM_REPAIR_STEPS>
    -n, --num-relaxation-steps <NUM_RELAXATION_STEPS>    [default: 0]
    -o, --output-file <OUTPUT_FILE>
    -r, --resolution-field <RESOLUTION_FIELD>
        --rmax <RMAX>                                    [default: 800]
        --rw <RW>                                        [default: 400]
    -s, --scaling <SCALING>                              [default: 10]
    -v, --vertices <VERTICES>
    -V, --version                                        Print version information

Algorithm

Basically a reimplementation of https://www.cs.cmu.edu/~quake/triangle.html with a customized triangle acceptance function based on the resolution.

Advantages of the algorithm:

Preserves the boundary perfectly
Guarantees a minimum angle
Takes seconds to minutes on a normal computer (code is not parallelized yet)

Disadvantages:

Currently does not look at area differences
May produce too many triangles (no cutoff for too small triangles, we can add that if we want that)

Steps:

The code starts by running a constrained Delaunay algorithm to preserve boundary edges.
The code then identifies all bad triangles. Bad triangles are triangles with too small minimum angle or where longest side is longer than the resolution at the triangle centroid.
It then starts repairing triangles starting from the worst (this is defined as the triangle with the largest circumradius, we can redefine that).
It repairs by creating smaller triangles which are again checked for their quality and if they are bad, they are placed on a priority queue for repair.
Code continues until no bad triangles are left for repair.

Lower bound on triangle size:

The code accepts triangles if their longest side is 0.2 times shortest edge of input polygons. This is to not get too many too tiny triangles when trying to make sure that all of them are above certain minimum angle.

Mesh improvement using mesh relaxation

After the mesh refinement steps (above), if --num-relaxation-steps is set and larger than 0, the code will run a series of iterations to relax vertex coordinates according to forces acting on vertices. For these steps the code uses --delta-t and --fscale and this algorithm corresponds to the Distmesh algorithm with the modification that all boundary vertices are fixed.

Visualize the result

This code produces a grid file, for instance grid.txt which I look at in my browser using https://gitlab.com/oceanbox/makai.

Resolution field

If a resolution field file is provided, the code will interpolate the data using bilinear interpolation using the library bilinear and resolution in a point is then defined as the minimum of the resolution coming from the distance function and the resolution approximated using the interpolation.

The file format for the resolution field file:

ox oy       # origin
sx sy       # step along x and y
xi yi z     # index along x (counting starts with 0), index along y (counting starts with 0), and value
xi yi z
...         # as many of these as we like

Only non-negative indices are allowed.
We don't have to give all values in places where the matrix has "holes". This can be relevant "on land" where we might not have meaningful data.
The rows with indices and values can come in any order.
It does not matter how many spaces are between numbers.

Example for a 3x2 grid:

When interpolating, the code first finds the right "tile", then checks that we have data available for all 4 corners. If yes, then the value is interpolated and used. If one or more of the 4 corner values are missing, then this tile is not considered and does not contribute to the resolution.

How to triangulate a set of vertices

Maybe you have polygons for boundary and islands and a set of vertices "in water" and want to Delaunay-triangulate them without any mesh refinement and without any mesh relaxation.

You can do it like this:

#!/usr/bin/env bash

set -euf -o pipefail

binary=$PWD/../../target/release/meshrefine

time $binary --boundary=boundary.txt \
             --islands=islands.txt \
             --vertices=vertices.txt \
             --num-repair-steps=0 \
             --output-file=grid.txt

Here we have set the number of repair steps to 0 and we provide a --vertices file which contains a list of vertices (x and y coordinates on each line).

You can find an example for this in folder examples/delaunay-only.