Skip to content

bast/smeshing

cpp-version
Switch branches/tags
Code

Latest commit

 

Git stats

Files

Permalink
Failed to load latest commit information.
Type
Name
Latest commit message
Commit time
img
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
https://travis-ci.org/bast/smeshing.svg?branch=master https://coveralls.io/repos/github/bast/smeshing/badge.svg?branch=master

Simple mesh generator in Python (fork of PyDistMesh)

img/example.jpg

This version of the code is based on PyDistMesh developed by Bradley M. Froehle. This code generates unstructured triangular meshes using signed distance functions. Like DistMesh and PyDistMesh upon which this work is based, this code is distributed under the GNU GPL.

I am in the process of rewriting this code to Rust: https://github.com/bast/smeshing-rs

Nice things about this code

  • Individual components live in separate libraries.
  • A lot of effort was invested in avoiding quadratic scaling.
  • Optimization is fully relaxed.
  • Delaunay is performed at every step.
  • Good memory profile (hopefully, please report if not).
  • Gives the user a lot of flexibility to define a distance-dependent resolution.

Known issues

  • Code uses shared-memory parallelization but the load leveling is not optimal and the scaling has not been studied in detail.

Recommended citation

If you use this code in a program or publication, please cite the following references:

The algorithm: - `P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB, SIAM

Review, Volume 46 (2), pp. 329-345, June 2004 <http://persson.berkeley.edu/distmesh/persson04mesh.pdf>`__

The present code:

@misc{smeshing,
  author    = {Bast, Radovan},
  title     = {Smeshing: Mesh generator for unstructured triangular grids},
  month     = {4},
  year      = {2020},
  publisher = {Zenodo},
  version   = {v0.2.0},
  doi       = {10.5281/zenodo.3774421},
  url       = {https://doi.org/10.5281/zenodo.3774421}
}

References

The DistMesh algorithm is described in the following two references:

Installation

Installation using virtualenv

virtualenv venv
source venv/bin/activate

pip install --process-dependency-links git+https://github.com/bast/smeshing.git

smesh --help

Installing dependencies for development

virtualenv venv
source venv/bin/activate
pip install -r requirements.txt

Installation on the Stallo supercomputer

#!/bin/bash

#SBATCH --account=your-account
#SBATCH --job-name=install
#SBATCH --ntasks=1
#SBATCH --time=0-01:00:00
#SBATCH --partition short
#SBATCH --mem-per-cpu=1000MB
#SBATCH --mail-type=ALL

module purge
module load foss/2016b
module load Python/3.5.2-foss-2016b
module load CMake/3.7.1-foss-2016b
module load libffi/3.2.1-foss-2016b

cd ${SLURM_SUBMIT_DIR}

python3 -m venv venv
source venv/bin/activate

python --version

export CC=gcc
export CXX=g++
export FC=gfortran

pip install --process-dependency-links git+https://github.com/bast/smeshing.git

Running tests

py.test -vv smeshing/*.py

How to run the code

Launching the code

The code is launched using the smesh script. Example:

$ smesh --boundary=/home/user/smeshing/data/happy-bear/boundary.txt \
        --islands=/home/user/smeshing/data/happy-bear/islands.txt \
        --config=/home/user/smeshing/data/happy-bear/config.yml \
        --output=data.txt

For an explanation of the options try:

$ smesh --help

Usage: smesh [OPTIONS]

Options:
  --boundary TEXT           File containing boundary data.
  --islands TEXT            Island file names (it is possible to use
                            wildcards).
  --resolution-fields TEXT  File name(s) containing resolution fields (it is
                            possible to use wildcards).
  --config TEXT             Read configuration from this file.
  --output TEXT             Write output to this file.
  --restart TEXT            Restart from this file.
  --help                    Show this message and exit.

You can take the files here as a starting point: https://github.com/bast/smeshing/tree/master/data/happy-bear

Example run script for the Stallo supercomputer

#!/bin/bash

#SBATCH --account=your-account
#SBATCH --job-name=smesh
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=20
#SBATCH --exclusive
#SBATCH --time=0-00:10:00
#SBATCH --partition short
#SBATCH --mem-per-cpu=500MB
#SBATCH --mail-type=ALL

# load a couple of modules
module purge
module load foss/2016b
module load Python/3.5.2-foss-2016b
module load libffi/3.2.1-foss-2016b

# use all available threads for shared-memory parallelization
export OMP_NUM_THREADS=${SLURM_TASKS_PER_NODE}

# compile the custom functions
cd ${SLURM_SUBMIT_DIR}
g++ -O3 -shared -fpic custom_functions.cpp -o libcustom_functions.so

# this will define the custom functions for the meshing code
export LD_PRELOAD=${SLURM_SUBMIT_DIR}/libcustom_functions.so

# load the virtual environment that contains the installation
source /home/user/smeshing/venv/bin/activate

# start the actual code
smesh --boundary=${SLURM_SUBMIT_DIR}/boundary.txt \
      --islands=${SLURM_SUBMIT_DIR}/islands.txt \
      --config=${SLURM_SUBMIT_DIR}/config.yml \
      --output=${SLURM_SUBMIT_DIR}/data.txt

exit 0

How to provide polygon data for the boundary and islands

Boundary polygon data has to be in a separate file from island data but both are given in the same format. Island data polygons can be all in one file, or in multiple files. Each polygon starts with one line specifying the number of points, followed by the polygon points, each point in one line. First and last point of the polygon have the same coordinates.

As an example, this file contains two polygons, one with 5 points, one with 4 points:

5
0.0 0.0
1.0 0.0
1.0 1.0
0.0 1.0
0.0 0.0
4
5.0 0.0
6.0 0.0
6.0 1.0
5.0 0.0

It would be equally fine to split this file into two files if you prefer.

Configuration

Configuration is given in YAML format. You can name the configuration file as you like, for instance config.yml. The order of keywords does not matter and you can add comments as in this example:

# number of grid points
num_grid_points: 5000

# number of all boundary and coastline interpolation points
# these will not be part of the grid points
# instead of num_interpolation_points you can also provide
# interpolation_step_length using the same units as the coordinates of your data
num_interpolation_points: 1000

# number of iterations
num_iterations: 100

Restart

It is possible to restart a calculation if you provide --restart=/path/to/restart/file.

How is the resolution computed?

The resolution is computed as minimum of the resolution function and the optional resolution field (below). But this needs some explanation and background so let's start simple:

One approach would be to define the resolution as the distance to islands and the boundary but this would lead to a couple of problems:

  • Resolution would decrease to zero close to polygons and lead to numerical problems.
  • We would see many grid points at the boundary.
  • We would treat the entire coastline and all islands on the same footing but perhaps some portions are scientifically more interesting and require a finer mesh than others?

So we decided to make two improvements:

  • We wanted to be able to make some polygon points more "attractive" for a finer mesh than others. For this we made it possible to define islands not only as points with x and y coordinates but also to give a coefficient which we can use in the resolution function. For this we introduced function g. This function also allows to cap the minimum to avoid zero resolution.
  • We wanted to make it possible for the resolution function to depend non-linearly on the distance to the nearest coastal point. For this we introduced function h.

User has the possibility to express g and h and the resolution is defined as the sum of both (see below).

One problem remains: this approach does not allow to have finer grid depending on local features which are not related to the coastline, such as water depth or other local data. For this we introduced the optional resolution field (see below).

If the resolution field is provided, the code will take the resolution field value in the closest point to the reference point, the code will also compute g + h at the reference point, and use the minimum of both values as resolution.

How to express the resolution function

Grid points move depending on forces and forces depend on the resolution. You have to define the resolution yourself by writing a C++ file, compiling it, and feeding it to the meshing algorithm using LD_PRELOAD. To get you started, here is an example custom_functions.cpp:

// provides std::max
#include <algorithm>

// The resolution is expressed as distance using the same (arbitrary)
// units as used by the boundary and polygons - this means that
// larger resolution number means that points are farther apart.

// Resolution in point r is defined as min(f(r, p)), where the miminum
// is taken over all boundary points p for a particular boundary point p,
// f is given as f(r, p) = g(d(r, p)) + h(c_p).
// d(r, p) is the distance r to p and function h(c_p) depends on
// coefficients c_p of a boundary point p. The number of coefficients
// per point and their meaning can be freely specified and interpreted.

// Below you are asked to specify functions g and h.
// You have two restrictions:
// 1) You have to respect is that g should not decrease for an increasing d.
//    In other words, for an increasing distance the resolution should not
//    decrease.
// 2) The sum g + h should never become zero since the code will divide by
//    the distance.

// This function only depends on the distance to a boundary point but not
// on coefficients at the boundary point.
double g_function(const double distance)
{
    // this is to make sure we do not end up with zero distance
    // and then try to divide by zero later
    double result = std::max(0.5, distance);

    return result;
}

// The code will give you all coefficients for a point in h_function
// and then you can use and combine them freely.
double h_function(const double coefficients[])
{
    // in this example we simply return the first coefficient
    return coefficients[0];
}

To see how this file can be compiled and provided to the meshing script, please have a look at the run script example.

Adding resolution fields

Sometimes the resolution should not only be dependent on the distance and the boundary coefficients but also on local features. For this you can provide resolution fields with the --resolution-fields flag. Point it to a file or files that contain the following format:

N
x1 y1 r1
x2 y2 r2
...
xN yN rN

The code will then use this field to interpolate a local resolution for each of the resolution fields. The resolution for a grid point is then given as the minimum taken over all resolution fields and the distance-dependent resolution provided by g_function + h_function.

Design choices

Why do we need to provide islands and the boundary separately?

  • We compute view vectors for nearest neighbor polygon points in view. For the boundary they point to the "inside". For islands they point to the "outside".
  • During the computation we need to figure out whether points are inside or outside of polygons. We want grid points to be inside the boundary but outside islands.

Postprocessing

The repository contains a tiny script which can be used to plot the generated grid:

python smeshing/plot.py data/happy-bear/result.txt example.png

About

Mesh generator for unstructured triangular grids.

Topics

Resources

License

Stars

Watchers

Forks

Packages

No packages published