gridgen.c

gridgen是一款强大的网格生成程序

/******************************************************************************
 *
 * File:           gridgen.c
 *
 * Created:        19/03/2001
 *
 * Author:         Pavel Sakov
 *                 CSIRO Marine Research
 *
 * Purpose:        Code for generation of 2D orthogonal grids.
 *
 * Description:    This program is based on the CRDT algorithm as described
 *                 in:
 *                 Tobin D. Driscoll and Stephen A. Vavasis, "Numerical
 *                 conformal mapping using cross-ratios and Delaunay
 *                 triangulation", SIAM J. Sci. Comp., 1998, 19(6), 1783-1803.
 *                 See http://www.math.udel.edu/~driscoll/pubs/crdt.pdf.
 *                 
 * Revisions:      02/11/2006 PS Minor structural changes associated with
 *                            the introduction of libgridgen.a. Changed the
 *                            structure name "gridspec" to "gridgen".
 *                 15/02/2007 PS Introduced a new function gridgen_create2().
 *                            Some associated structural changes.
 *
 *****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <ctype.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include <assert.h>
#include <errno.h>
#include "gridgen.h"
#include "ode.h"
#include "istack.h"
#include "delaunay.h"
#include "vertlist.h"
#include "swcr.h"
#include "broyden.h"
#include "geom.h"
#include "issimplepoly.h"
#include "hash.h"
#include "version.h"
#include "config.h"
#include "nan.h"

#define BUFSIZE 1024
#define NMIN 3                  /* minimal number of vertices */
#define NMAX 10000              /* maximal number of vertices */
#define M_2PI (M_PI * 2.0)
#define EPS 1.0e-6
#define EPS_DEF 1.0e-12         /* default precision */
#define EPS_MAX 1.0e-3
#define EPS_MIN 1.0e-14
#define NNODES_MIN 5
#define NNODES_MAX 20
#define NNODES_DEF 10           /* default number of nodes */
#define N_DEF 25
#define N_MIN 2                 /* min. number of nodes in grid edge */
#define N_MAX 10001             /* max. number of nodes in grid edge */
#define BIGDOUBLE 1.0e+6
#define THIN_DEF 1              /* thin input vertices by default */
#define CHECKSIMPLEPOLY_DEF 1   /* check input for self-intersections */
#define NEWTON_DEF 1            /* newton method by default */
#define M 4                     /* number of stored iterations in generalised
                                 * Broyden update */
#define NPPE_DEF 3              /* number of points per internal edge in
                                 * polygon images */
#define ZZERO 0.0 + 0.0 * I

int gg_verbose = 0;
int tr_verbose = 0;

/* The diagonal goes from vids[0] to vids[2];
 * tids[0] corresponds to vids[0], vids[1] and vids[2];
 * tids[1] corresponds to vids[4], vids[2], vids[3].
 */
typedef struct {
    int id;
    int vids[4];
    int tids[2];

    /*
     * quadrilaterals which share with this one one of the triangles 
     */
    int nneighbours;
    int* neighbours;
} quadrilateral;

typedef struct {
    char* prmfname;
    char* datafname;
    char* sigmafname;
    char* rectfname;
    FILE* out;
#if !defined (GRIDGEN_STANDALONE) && defined(HAVE_GRIDNODES_H)
    gridnodes* gn;
#endif
    FILE* fsigma;
    FILE* frect;
    int* nrect;
    double** xrect;
    double** yrect;

    double xmin;
    double ymin;
    double xmax;
    double ymax;

    int thin;
    int checksimplepoly;
    int reversed;

    int nx;
    int ny;

    int ngridpoints;
    zdouble* gridpoints;

    /*
     * specifications 
     */
    double eps;
    int nnodes;
    int newton;

    /*
     * z 
     */
    vertlist* vertices;
    delaunay* d;
    int nquadrilaterals;
    quadrilateral* quadrilaterals;
    int lastqi;

    int nz;
    zdouble* zs;

    /*
     * z<->w 
     */
    double* betas;
    swcr* sc;
    double* cis;
    int* nsigmas;
    double** sigmas;
    zdouble* ws;
    zdouble* As;
    zdouble* Bs;

    int ncorners;               /* number of corners in the mapped polygonal; 
                                 * 4 by default */

    /*
     * z1 
     */
    zdouble* rzs;
    int* rvids;

    double newxmin;
    double newymin;
    double newxmax;
    double newymax;

    /*
     * w<->z1 
     */
    double* newbetas;
    swcr* newsc;
    zdouble* newAs;
    zdouble* newBs;
    zdouble* newzs;
    zdouble* newrzs;

    /*
     * What we have is a number of mappings, each working best in a certain
     * quadrilateral of the complex plane of the original polygon. It looses
     * a bit of precision for the neighbour quadrilaterals (but still works 
     * OK for them), and can deteriorate substantially or completely for points
     * in distant quadrilaterals.
     *
     * Now, let one map a point from the complex plane of the image to the 
     * complex plane of the original. For best results, one needs to find the
     * best mapping to be used. To find it, one needs to know which
     * quadrilateral in the original plane the inverse transformation of a 
     * point in the image belongs to.
     *
     * The problem is that one does not know the images of the quadrilaterals, 
     * only the vertice images.
     *
     * To get coarse images of the quadrilaterals, we put gg->nppe points on 
     * the diagonal of each quadrilateral, and map them to the image region. 
     * The original quadrilaterals would then effectively transform into a
     * number of non-intersecting polygons in the image plane.
     *
     * (Each interior edge of the triangulation is a diagonal. By mapping
     * all diagonals, we map all interior edges of all quadrilaterals. In
     * this way, we ensure that each internal edge is mapped only once.
     * We then need to be able to get effectively the diagonal id for a given
     * edge of a quadrilateral. This is done via hashtable that allows to get 
     * the diagonal id from the edge vertex ids.)
     *
     * After that, almost any point in the image space can be easily mapped to
     * a certain quadrilateral. In rare cases when we will get a neighbour of
     * the right quadrilateral instead of the quadrilateral itself (because of
     * approximating the boundary with a polyline) this still will be good
     * enough to map the point.
     *
     * The mapping is controled by parameter `nppe'. To switch it off, set
     * nppe = 0 in the parameter file. In this case, coarser images of 
     * quadrilaterals will be used that are formed by straight lines
     * connecting the quadrilateral vertex images.
     */

    /*
     * quasi-images of quadrilaterals
     */
    int nppe;                   /* number of points per internal edge */
    int nppq;                   /* number of points allocated per
                                 * quadrilateral (= nppe * 4 + 5) */
    int* nqivertices;           /* number of vertices in this image
                                 * [nquadrilaterals] */
    zdouble* qivertices;        /* vertices [nquadrilaterals * nppq] */

    /*
     * id of quadrilateral containing the last valid grid point
     */
    int lastqid;

    /*
     * temporal stuff -- some output of interest to me. Ignore it. 
     */
    int diagonal;
} gridgen;

void gridgen_setverbose(int verbose)
{
    if (verbose <= 0) {
        gg_verbose = 0;
        tr_verbose = 0;
    } else if (verbose == 1) {
        gg_verbose = 1;
        tr_verbose = 0;
    } else if (verbose == 2) {
        gg_verbose = 2;
        tr_verbose = 1;
    } else if (verbose >= 3) {
        gg_verbose = 2;
        tr_verbose = 2;
    }
}

void gridgen_printversion(void)
{
    printf("gridgen version %s\n", gridgen_version);
}

void gridgen_printhelpalg(void)
{
    printf("  `gridgen' generates quasi-orthogonal grids for polygonal regions by CRDT\n");
    printf("algorithm. It is based on the Schwarz-Christoffel formula that allows to\n");
    printf("transform conformally a unit circumcircle into a simply connected polygon with\n");
    printf("specified vertex angles.\n");
    printf("  In the transformation defined by Schwarz-Christoffel formula, each vertex of\n");
    printf("the polygon corresponds to a point on the circumcircle (\"prevertex\"). By\n");
    printf("specifying different angles, it is possible to transform the same set of\n");
    printf("prevertices into different polygons. A conjunction of one forward and one\n");
    printf("backward transformation with the same prevertices but different vertex angles in\n");
    printf("the Schwarz-Christoffel formula yields a single conformal mapping between the\n");
    printf("original polygon and a new (\"image\") polygon.\n");
    printf("  The underlying idea for using conformal mappings is that these mappings do\n");
    printf("preserve local angles. Therefore, if the image polygon is set to be a rectangle\n");
    printf("or a set of intersecting rectangles, one can easily put an orthogonal grid on it\n");
    printf("and ...provided he can do the inverse transform... this orthogonal grid would\n");
    printf("yield a (quasi)-orthogonal grid in the original polygon.\n");
    printf("  The main obstacle to this idyllic picture is a potential loss of precision on\n");
    printf("the way, when a number of vertices degenerate into a single prevertex\n");
    printf("(\"crowding\"). The effect of crowding is particularly strong for polygons with\n");
    printf("long (thin) arms, and this is where the CRDT algorithm comes into play.\n");
    printf("  The essence of the CRDT algorithm is in splitting the original polygon into a\n");
    printf("number of intersecting quadrilaterals and generating not one but a number of\n");
    printf("conformally equivalent transformations (\"embeddings\") of the original polygon\n");
    printf("to the unit circumcircle, each of them attributed with a certain quadrilateral\n");
    printf("and working numerically best for this and adjacent quadrilaterals.\n");
    printf("  All these details are transparent to a user of `gridgen', who has to define\n");
    printf("only the original and the image polygons. The image polygon is defined in a\n");
    printf("semi-implicit manner, by specifying new angles (but not edge lengths) for each\n");
    printf("vertex. In the simplest case, a user must specify which four vertices are to\n");
    printf("become the vertices of the image rectangle by marking them with \"1\" in the\n");
    printf("input file specifying the original polygon. After that, `gridgen' will calculate\n");
    printf("images of the nodes of a rectangular grid with specified parameters (or images\n");
    printf("of a custom set of points) within the image polygon in the original polygon, and\n");
    printf("this will complete the task.\n");
    printf("  There are 4 main stages in generating a quasi-orthogonal grid:\n");
    printf("  -- preprocessing of the input polygon by removing duplicate vertices and\n");
    printf("     inserting new ones when necessary (to avoid thin triangles);\n");
    printf("  -- calculation of the conformal mappings that transfer the input polygon\n");
    printf("     into a unit circumcircle;\n");
    printf("  -- calculation of conformal mappings of a unit circumcircle into a polygon\n");
    printf("     that transfer the prevertices into vertices of a polygon with specified\n");
    printf("     vertex angles;\n");
    printf("  -- inverse mapping of specified points (grid nodes) from the image polygon\n");
    printf("     into the original one.\n");
}

void gridgen_printhelpprm(void)
{
    printf("  `gridgen': Generates grid for a polygonal region by CRDT algorithm.\n");
    printf("  Parameter file format:\n");
    printf("    input <input polygon data file>\n");
    printf("    [output <output grid data file>]\n");
    printf("    [grid <input grid data file>]\n");
    printf("    [nx <number of nodes in X direction>] (default = 25)\n");
    printf("    [ny <number of nodes in Y direction>] (25)\n");
    printf("    [nnodes <number of nodes in Gauss-Jacobi quadrature>] (12)\n");
    printf("    [precision <precision>] (1e-10)\n");
    printf("    [thin {0|1}] (1)\n");
    printf("    [checksimplepoly {0|1}] (1)\n");
    printf("    [newton {0|1}] (1)\n");
    printf("    [sigmas <intermediate backup file>]\n");
    printf("    [rectangle <output image polygon file>]\n");
    printf("    [nppe <number of points per internal edge>] (3)\n");
    printf("  Input polygon data file format:\n");
    printf("    <x_1> <y_1> [<beta_1>[*]]\n");
    printf("    <x_2> <y_2> [<beta_2>[*]]\n");
    printf("    ...\n");
    printf("    <x_n> <y_n> [<beta_n>[*]]\n");
    printf("  Input grid data file format:\n");
    printf("    <x_1> <y_1>\n");
    printf("    ...\n");
    printf("    <x_n> <y_n>\n");
    printf("    where 0 <= x_i, y_i <= 1\n");
    printf("  Remarks:\n");
    printf("    1. The input polygon must be simply connected.\n");
    printf("    2. Values `beta_i' in the input polygon data file (see above) define angles\n");
    printf("       of the image region corners:\n");
    printf("       beta_i = (1 - angle_i / pi) * 2; 0 by default,\n");
    printf("       where angle_i -- value of the interior angle at the i_th vertex. To get\n");
    printf("       a valid polygon, the sum of beta_i must be equal to 4. The image of the\n");
    printf("       very first vertex marked by \"*\" is placed in point (0,1), while the\n");
    printf("       edge connecting it with the next vertex is set directed towards\n");
    printf("       -i*infty. If not specified, beta_i is set to 0.\n");
    printf("    3. A uniform grid is specified by `nx' and `ny' parameters. This may be\n");
    printf("       overridden by specifying a custom grid from an input grid file. Such a\n");
    printf("       file should contain the desired node coordinates in index space scaled\n");
    printf("       to a 1x1 square.\n");
    printf("    4. The `nnodes' parameter affects precision and computation time; it is\n");
    printf("       advised to be set to the same value as -log10(<precision>) or slightly\n");
    printf("       higher.\n");
    printf("    5. If `newton 1' specified, the nonlinear system for sigmas is solved using\n");
    printf("       Gauss-Newton solver with Broyden update; otherwhile simple iterations\n");
    printf("       are used.\n");
    printf("    6. If no output file specified, the results are written to the standard\n");
    printf("       output.\n");
    printf("    7. If `sigmas' specified, the sigmas (an intermediate solution of nonlinear\n");
    printf("       system) will be read from/saved to a binary file. This may save a fair\n");
    printf("       bit of time in the case of multiple runs of the program.\n");
    printf("    8. If `rectangle' specified, the image polygon vertex coordinates will be\n");
    printf("       written to the corresponding text file.\n");
    printf("    9. Parameter `nppe' controls the coarseness of mapping of quadrilaterals.\n");
    printf("       Large values can take some time; small values can lead to failing\n");
    printf("       to map some points in difficult cases when quadrilateral images are\n");
    printf("       strongly distorted.\n");
    printf("  Acknowledgments. This program uses the following public code/algorithms:\n");
    printf("    1. CRDT algorithm by Tobin D. Driscoll and Stephen A. Vavasis -- for\n");
    printf("       conformal mapping.\n");
    printf("    2. `Triangle' by Jonathan Richard Shewchuk -- for Delaunay triangulation.\n");
    printf("    3. `SCPACK' by Lloyd N. Trefethen -- for Schwarz-Christoffel transform.\n");
    printf("    4. `DOPRI8' by E. Hairer, S.P. Norsett, G. Wanner -- for solving ODEs.\n");
    printf("    5. Shamos-Hoey algorithm implementation by softSurfer -- for testing a\n");
    printf("       polyline on self-intersections.\n");
    printf("\n");
    printf("  Good luck!\n");
    printf("\n");
    printf("  Pavel Sakov\n");
    printf("  April-May 2001 - January 2006\n");
}

static void quit(char* format, ...)
{
    va_list args;

    fflush(stdout);

    fprintf(stderr, "error: ");
    va_start(args, format);
    vfprintf(stderr, format, args);
    va_end(args);
    exit(1);
}

static FILE* gg_fopen(char* fname, char* mode)
{
    FILE* f = fopen(fname, mode);

    if (f == NULL)
        quit("could not open \"%s\" in \"%s\" mode: %s\n", fname, mode, strerror(errno));

    return f;
}

/* This is modified skipToKeyEnd() from `sjwlib'.
 */
static int key_find(char* fname, FILE* fp, char* key)
{
    char buf[BUFSIZE];
    int len = strlen(key);
    int fpos;
    char* s;
    int rewound = 0;

    do {
        fpos = ftell(fp);
        s = fgets(buf, BUFSIZE, fp);
        if (s == NULL) {
            if (!rewound) {
                fpos = 0;
                if (fseek(fp, fpos, 0))
                    quit("%s: could not rewind: %s\n", fname, strerror(errno));
                s = fgets(buf, BUFSIZE, fp);
                rewound = 1;
            } else
                return 0;
        }

        if (s == NULL)
            break;