gridgen.c

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

        while (isspace(s[0]))
            s++;

        if (isspace(s[len]))
            s[len] = 0;
    } while (strcmp(key, s) != 0);

    if (s == NULL)
        return 0;

    if (fseek(fp, fpos + (s - buf) + len, 0))
        quit("%s: no characters found after \"%s\"\n", fname, key);

    return 1;
}

/* This is modified readkeyprm_s() from `sjwlib'.
 */
static int prm_read(char* fname, FILE* fp, char* key, char* p)
{
    char line[BUFSIZE];
    char* s = line;
    char* r = p;

    if (!key_find(fname, fp, key))
        return 0;

    if (fgets(line, BUFSIZE, fp) != line)
        quit("%s: key \"%s\": read failed: %s\n", fname, key, strerror(errno));

    for (s = line; *s && isspace(*s); s++);

    for (r = p; *s && (*s != '\n'); *r++ = *s++);

    *r = 0;

    if (gg_verbose)
        fprintf(stderr, "-> %s = \"%s\"\n", key, p);

    return 1;
}

static gridgen* gridgen_init(void)
{
    gridgen* gg = malloc(sizeof(gridgen));

    gg->prmfname = NULL;
    gg->datafname = NULL;
    gg->sigmafname = NULL;
    gg->rectfname = NULL;
    gg->out = NULL;
#if !defined (GRIDGEN_STANDALONE) && defined(HAVE_GRIDNODES_H)
    gg->gn = NULL;
#endif
    gg->fsigma = NULL;
    gg->frect = NULL;
    gg->nrect = NULL;
    gg->xrect = NULL;
    gg->yrect = NULL;
    gg->xmin = DBL_MAX;
    gg->xmax = -DBL_MAX;
    gg->ymin = DBL_MAX;
    gg->ymax = -DBL_MAX;
    gg->thin = THIN_DEF;
    gg->checksimplepoly = CHECKSIMPLEPOLY_DEF;
    gg->reversed = 0;
    gg->nx = -1;
    gg->ny = -1;
    gg->ngridpoints = 0;
    gg->gridpoints = NULL;
    gg->eps = EPS_DEF;
    gg->nnodes = NNODES_DEF;
    gg->newton = NEWTON_DEF;
    gg->vertices = NULL;
    gg->d = NULL;
    gg->nquadrilaterals = 0;
    gg->lastqi = 0;
    gg->nz = 0;
    gg->zs = NULL;
    gg->betas = NULL;
    gg->sc = NULL;
    gg->nsigmas = NULL;
    gg->sigmas = NULL;
    gg->ws = NULL;
    gg->As = NULL;
    gg->Bs = NULL;
    gg->ncorners = -1;
    gg->rzs = NULL;
    gg->rvids = NULL;
    gg->newxmin = DBL_MAX;
    gg->newxmax = -DBL_MAX;
    gg->newymin = DBL_MAX;
    gg->newymax = -DBL_MAX;
    gg->newbetas = NULL;
    gg->newsc = NULL;
    gg->newAs = NULL;
    gg->newBs = NULL;
    gg->newzs = NULL;
    gg->newrzs = NULL;
    gg->nppe = NPPE_DEF;
    gg->nppq = gg->nppe * 4 + 5;
    gg->nqivertices = NULL;
    gg->qivertices = NULL;
    gg->lastqid = -1;
    gg->diagonal = -1;

    return gg;
}

static gridgen* gridgen_create(char* prmfname)
{
    gridgen* gg = gridgen_init();
    char buf[BUFSIZE];
    FILE* prm = NULL;
    FILE* data = NULL;
    int nold;

    prm = gg_fopen(prmfname, "r");
    gg->prmfname = prmfname;

    gg->vertices = vertlist_create();

    if (!prm_read(prmfname, prm, "input", buf))
        quit("%s: key \"input\" not found\n", prmfname);
    gg->datafname = strdup(buf);
    data = gg_fopen(buf, "r");

    if (gg_verbose)
        fprintf(stderr, "reading:\n");
    vertlist_read(gg->vertices, data);

    if (gg->vertices->n == 0)
        quit("%s: no XY points found\n", buf);
    if (gg->vertices->n < NMIN)
        quit("%s: %d points found -- too few (should be >= %d)\n", buf, gg->vertices->n, NMIN);
    if (gg->vertices->n > NMAX)
        quit("%s: %d points found -- too many (should be <= %d)\n", buf, gg->vertices->n, NMAX);
    if (gg_verbose)
        fprintf(stderr, "  %d vertices after read\n", gg->vertices->n);

    if (gg_verbose > 1)
        vertlist_print(gg->vertices, stderr);

    gg->reversed = 0;
    if (vertlist_area(gg->vertices) < 0.0) {
        if (gg_verbose)
            fprintf(stderr, "  clockwise, reversing\n");
        vertlist_change_order(&gg->vertices);
        gg->reversed = 1;
    } else if (gg_verbose)
        fprintf(stderr, "  counterclockwise ok\n");

    vertlist_find_minmax(gg->vertices, &gg->xmin, &gg->xmax, &gg->ymin, &gg->ymax);

    nold = gg->vertices->n;
    if (prm_read(prmfname, prm, "thin", buf))
        gg->thin = atoi(buf);
    if (gg->thin) {
        if (gg_verbose)
            fprintf(stderr, "thinning:\n");
        vertlist_thin(gg->vertices, (gg->xmax - gg->xmin) / BIGDOUBLE, (gg->ymax - gg->ymin) / BIGDOUBLE);

        if (gg->vertices->n < NMIN)
            quit("%s: %d vertices after thinning -- too little (should be >= %d)\n", buf, gg->vertices->n, NMIN);

        if (gg_verbose)
            fprintf(stderr, "  %d vertices after thinning\n", gg->vertices->n);
    }

    if (prm_read(prmfname, prm, "checksimplepoly", buf))
        gg->checksimplepoly = atoi(buf);
    if (gg->checksimplepoly) {
        double* x = NULL;
        double* y = NULL;

        if (gg_verbose)
            fprintf(stderr, "checking for self-intersections:\n");
        vertlist_toxy(gg->vertices, &x, &y);
        if (!issimplepolygon(gg->vertices->n, x, y)) {
            fprintf(stderr, "  Beware that code for testing a polyline on self-intersections is still in beta\n");
            fprintf(stderr, "  stage. You may try to skip the test by adding parameter \"checksimplepoly 0\"\n");
            fprintf(stderr, "  to your paremeter file. (If true, you should get a segfault during\n");
            fprintf(stderr, "  triangulation.)\n");
            quit("%s: not a simple polyline\n", gg->datafname);
        }
        free(x);
        free(y);
    }

    if (gg_verbose > 1 && gg->vertices->n != nold)
        vertlist_print(gg->vertices, stderr);

    if (prm_read(prmfname, prm, "nnodes", buf))
        gg->nnodes = atoi(buf);
    if (gg->nnodes < NNODES_MIN)
        gg->nnodes = NNODES_MIN;
    if (gg->nnodes > NNODES_MAX)
        gg->nnodes = NNODES_MAX;
    if (gg_verbose)
        fprintf(stderr, "nnodes = %d\n", gg->nnodes);

    if (prm_read(prmfname, prm, "newton", buf))
        gg->newton = atoi(buf);

    if (prm_read(prmfname, prm, "precision", buf))
        gg->eps = atof(buf);
    if (gg->eps < EPS_MIN)
        gg->eps = EPS_MIN;
    if (gg->eps > EPS_MAX)
        gg->eps = EPS_MAX;
    if (gg_verbose)
        fprintf(stderr, "precision = %3g\n", gg->eps);

    if (prm_read(prmfname, prm, "grid", buf)) {
        FILE* gridfile = gg_fopen(buf, "r");
        vertlist* l = vertlist_create();

        if (gridfile == NULL)
            quit("could not open \"%s\": %s\n", buf, strerror(errno));

        vertlist_read(l, gridfile);
        gg->gridpoints = vertlist_tozdouble(l);
        gg->ngridpoints = l->n;

        fclose(gridfile);
        vertlist_destroy(l);
    }

    if (gg->gridpoints == NULL) {
        if (prm_read(prmfname, prm, "nx", buf))
            gg->nx = atoi(buf);
        else
            gg->nx = N_DEF;
        if (prm_read(prmfname, prm, "ny", buf))
            gg->ny = atoi(buf);
        else
            gg->ny = N_DEF;
        if (gg->nx < N_MIN)
            gg->nx = N_MIN;
        if (gg->nx > N_MAX)
            gg->nx = N_MAX;
        if (gg->ny < N_MIN)
            gg->ny = N_MIN;
        if (gg->ny > N_MAX)
            gg->ny = N_MAX;

        if (gg_verbose)
            fprintf(stderr, "going to generate %dx%d grid\n", gg->nx, gg->ny);
    }

    if (prm_read(prmfname, prm, "output", buf))
        gg->out = gg_fopen(buf, "w");
    else
        gg->out = stdout;

    if (prm_read(prmfname, prm, "sigmas", buf)) {
        gg->sigmafname = strdup(buf);
        gg->fsigma = fopen(buf, "r");
    }

    if (prm_read(prmfname, prm, "rectangle", buf)) {
        gg->rectfname = strdup(buf);
        gg->frect = gg_fopen(buf, "w");
    }

    if (prm_read(prmfname, prm, "ppe", buf))
        gg->nppe = atoi(buf);
    gg->nppq = gg->nppe * 4 + 5;

    fclose(data);
    fclose(prm);

    return gg;
}

#if !defined (GRIDGEN_STANDALONE) && defined(HAVE_GRIDNODES_H)

static gridgen* gridgen_create2(int nbdry, double xbdry[], double ybdry[], double beta[], int ul, int nx, int ny, int ngrid, double xgrid[], double ygrid[], int nnodes, int newton, double precision, int checksimplepoly, int thin, int nppe, int verbose, int* nsigmas, double** sigmas, int* nrect, double** xrect, double** yrect)
{
    gridgen* gg = gridgen_init();
    int nold;

    gridgen_setverbose(verbose);

    if (gg_verbose)
        fprintf(stderr, "reading:\n");
    gg->vertices = vertlist_create2(nbdry, xbdry, ybdry, beta);
    vertlist_setfirstnode(gg->vertices, ul);

    if (gg->vertices->n < NMIN)
        quit("%d points found -- too few (should be >= %d)\n", gg->vertices->n, NMIN);
    if (gg->vertices->n > NMAX)
        quit("%d points found -- too many (should be <= %d)\n", gg->vertices->n, NMAX);
    if (gg_verbose)
        fprintf(stderr, "  %d vertices after read\n", gg->vertices->n);

    if (gg_verbose > 1)
        vertlist_print(gg->vertices, stderr);

    gg->reversed = 0;
    if (vertlist_area(gg->vertices) < 0.0) {
        if (gg_verbose)
            fprintf(stderr, "  clockwise, reversing\n");
        vertlist_change_order(&gg->vertices);
        gg->reversed = 1;
    } else if (gg_verbose)
        fprintf(stderr, "  counterclockwise ok\n");

    vertlist_find_minmax(gg->vertices, &gg->xmin, &gg->xmax, &gg->ymin, &gg->ymax);

    nold = gg->vertices->n;
    gg->thin = thin;
    if (gg->thin) {
        if (gg_verbose)
            fprintf(stderr, "thinning:\n");
        vertlist_thin(gg->vertices, (gg->xmax - gg->xmin) / BIGDOUBLE, (gg->ymax - gg->ymin) / BIGDOUBLE);

        if (gg->vertices->n < NMIN)
            quit("%d vertices after thinning -- too little (should be >= %d)\n", gg->vertices->n, NMIN);

        if (gg_verbose)
            fprintf(stderr, "  %d vertices after thinning\n", gg->vertices->n);
    }

    gg->checksimplepoly = checksimplepoly;
    if (gg->checksimplepoly) {
        double* x = NULL;
        double* y = NULL;

        if (gg_verbose)
            fprintf(stderr, "checking for self-intersections:\n");
        vertlist_toxy(gg->vertices, &x, &y);
        if (!issimplepolygon(gg->vertices->n, x, y)) {
            fprintf(stderr, "  Beware that code for testing a polyline on self-intersections is still in beta\n");
            fprintf(stderr, "  stage. You may try to skip the test by adding parameter \"checksimplepoly 0\"\n");
            fprintf(stderr, "  to your paremeter file. (If true, you should get a segfault during\n");
            fprintf(stderr, "  triangulation.)\n");
            quit("%s: not a simple polyline\n", gg->datafname);
        }
        free(x);
        free(y);
    }

    if (gg_verbose > 1 && gg->vertices->n != nold)
        vertlist_print(gg->vertices, stderr);

    gg->nnodes = nnodes;
    if (gg->nnodes < NNODES_MIN)
        gg->nnodes = NNODES_MIN;
    if (gg->nnodes > NNODES_MAX)
        gg->nnodes = NNODES_MAX;
    if (gg_verbose)
        fprintf(stderr, "nnodes = %d\n", gg->nnodes);

    gg->newton = newton;

    gg->eps = precision;
    if (gg_verbose)
        fprintf(stderr, "precision = %3g\n", gg->eps);

    if (ngrid > 0) {
        int i;

        gg->ngridpoints = ngrid;
        gg->gridpoints = malloc(ngrid * sizeof(zdouble));

        for (i = 0; i < ngrid; ++i)
            gg->gridpoints[i] = xgrid[i] + I * ygrid[i];
    }

    if (gg->gridpoints == NULL) {
        gg->nx = nx;
        gg->ny = ny;
        if (gg->nx < N_MIN)
            gg->nx = N_MIN;
        if (gg->nx > N_MAX)
            gg->nx = N_MAX;
        if (gg->ny < N_MIN)
            gg->ny = N_MIN;
        if (gg->ny > N_MAX)
            gg->ny = N_MAX;

        if (gg_verbose)
            fprintf(stderr, "going to generate %dx%d grid\n", gg->nx, gg->ny);
    }

    if (gg->gridpoints == NULL)
        gg->gn = gridnodes_create(gg->nx, gg->ny, NT_NONE);
    else
        gg->gn = gridnodes_create(gg->ngridpoints, 1, NT_NONE);

    /*
     * Storage coordinates of the working stuff
     */
    gg->nsigmas = nsigmas;
    gg->sigmas = sigmas;
    gg->nrect = nrect;
    gg->xrect = xrect;
    gg->yrect = yrect;

    return gg;
}

#endif

static void get_rpoints(gridgen* gg)
{
    vertlist* l = gg->vertices;
    vertnode* now = l->first;
    int n = l->n;
    double sum, mult;
    int i, count;

    if (gg_verbose) {
        fprintf(stderr, "getting marked vertices:\n");
        fflush(stderr);
    }

    /*
     * count corner markers 
     */
    for (i = 0, count = 0, sum = 0.0; i < n; ++i) {
        point* p = &now->p;

        if (p->z != 0) {
            if (fabs(p->z) > 2.0)
                quit("vertex %d (%.15g,%.15g): |beta| = %.15g > 2\n", i, p->x, p->y, fabs(p->z));
            count++;
            sum += p->z;
        }
        now = now->next;
    }