gridgen.c

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

        }
        for (i = 0; i < 4; ++i) {
            if (vids[i] == vid0)
                qvid0 = i;
            if (vids[i] == vid1)
                qvid1 = i;
        }
        B[qid0] = -I * 1.0 / (dzetas[qvid1] - dzetas[qvid0]);
        A[qid0] = 1.0 * I - B[qid0] * dzetas[qvid0];

        for (i = 0; i < 4; ++i)
            gg->newzs[vids[i]] = A[qid0] + B[qid0] * dzetas[i];
    }

    calculate_zs(nq, gg->quadrilaterals, nz, gg->newzs, qid0, -1, gg->newsc, gg->ws, &count, done, gg->newAs, gg->newBs);       /* recursive 
                                                                                                                                 */

    for (i = 0; i < gg->ncorners; ++i) {
        int vid = gg->rvids[i];

        gg->newrzs[i] = gg->newzs[vid];
    }

    if (gg_verbose > 1) {
        fprintf(stderr, "image vertices:\n");
        for (i = 0; i < nz; ++i)
            fprintf(stderr, "  %d: (%.15g, %.15g)\n", i, creal(gg->newzs[i]), cimag(gg->newzs[i]));
        fprintf(stderr, "image corners:\n");
        for (i = 0; i < gg->ncorners; ++i)
            fprintf(stderr, "  %d: (%.15g, %.15g)\n", i, creal(gg->newrzs[i]), cimag(gg->newrzs[i]));
        fprintf(stderr, "aligning:\n");
    }

    align_newzs(gg->ncorners, gg->newrzs, gg->rvids, gg->nz, gg->newzs, gg->newbetas, gg->eps);

    if (gg_verbose > 1) {
        fprintf(stderr, "image vertices:\n");
        for (i = 0; i < nz; ++i)
            fprintf(stderr, "  %d: (%.15g, %.15g)\n", i, creal(gg->newzs[i]), cimag(gg->newzs[i]));
        fprintf(stderr, "image corners:\n");
        for (i = 0; i < gg->ncorners; ++i)
            fprintf(stderr, "  %d: (%.15g, %.15g)\n", i, creal(gg->newrzs[i]), cimag(gg->newrzs[i]));
    }

    find_zminmax(gg->ncorners, gg->newrzs, &gg->newxmin, &gg->newxmax, &gg->newymin, &gg->newymax);
    if (gg_verbose)
        fprintf(stderr, "  conformal modulus = %.5g\n", (gg->newxmax - gg->newxmin) / (gg->newymax - gg->newymin));

    if (gg->frect != NULL) {
        if (gg_verbose)
            fprintf(stderr, "saving image region:\n");
        write_zarray(gg->ncorners, gg->newrzs, gg->frect);
    }
    if (gg->nrect != NULL) {
        *gg->nrect = gg->ncorners;
        if (*gg->xrect != NULL)
            free(*gg->xrect);
        *gg->xrect = malloc(sizeof(double) * gg->ncorners);
        if (*gg->yrect != NULL)
            free(*gg->yrect);
        *gg->yrect = malloc(sizeof(double) * gg->ncorners);

        for (i = 0; i < gg->ncorners; ++i) {
            (*gg->xrect)[i] = creal(gg->newrzs[i]);
            (*gg->yrect)[i] = cimag(gg->newrzs[i]);
        }
    }

    free(done);
}

static int z2q_simple(gridgen* gg, zdouble z, zdouble zs[], double eps)
{
    int nq = gg->nquadrilaterals;
    int i, qi;

    for (i = 0, qi = gg->lastqi; i < nq; ++i, qi = (qi + 1) % nq) {
        quadrilateral* q = &gg->quadrilaterals[qi];

        if (in_quadrilateral(&z, zs, q->vids, eps)) {
            gg->lastqi = qi;
            return qi;
        }
    }

    return -1;
}

static int z2q(gridgen* gg, zdouble z, zdouble zs[], double eps)
{
    if (gg->nppe > 0) {
        int nq = gg->nquadrilaterals;
        int i, qi;

        for (i = 0, qi = gg->lastqi; i < nq; ++i, qi = (qi + 1) % nq) {
            if (in_poly(z, gg->nqivertices[qi], &gg->qivertices[qi * gg->nppq], eps)) {
                gg->lastqi = qi;
                return qi;
            }
        }
    } else
        return z2q_simple(gg, z, zs, eps);

    return -1;
}

static void output_point(gridgen* gg, double x, double y)
{
#if !defined(GRIDGEN_STANDALONE) && defined(HAVE_GRIDNODES_H)
    if (gg->gn == NULL) {
#endif
        if (!isnan(x))
            fprintf(gg->out, "%.10g %.10g\n", x, y);
        else
            fprintf(gg->out, "NaN NaN\n");
        fflush(gg->out);
#if !defined(GRIDGEN_STANDALONE) && defined(HAVE_GRIDNODES_H)
    } else
        gridnodes_readnextpoint(gg->gn, x, y);
#endif
}

static void map_point(gridgen* gg, zdouble z)
{
    zdouble zz = NaN;
    double eps = gg->eps;
    quadrilateral* qs = gg->quadrilaterals;
    int nz = gg->nz;
    zdouble* zs = gg->newzs;
    zdouble* ws = gg->ws;

    int status;
    zdouble w;

    int qid = z2q(gg, z, zs, eps);
    int vid = -1;
    int vidmin = -1;
    double distmin = DBL_MAX;
    int* vids;
    int i;

    if (qid < 0)
        quit("map_point(): could not map z = (%.10g,%.10g) to a quadrilateral\n", creal(z), cimag(z));

    vids = qs[qid].vids;

    /*
     * find the closest vertex
     */
    for (i = 0; i < 4; ++i) {
        double dist = cabs(z - zs[vids[i]]);

        if (dist < eps) {
            w = gg->ws[qid * nz + vids[i]];
            vid = vids[i];
            status = 0;
            goto know_w;
        }
        if (dist < distmin) {
            vidmin = vids[i];
            distmin = dist;
        }
    }

    /*
     * (Note that because we eliminated the possibility of z being the nearest
     * vertex, w can not be a transform of a polygon vertex. The only danger
     * is that because of approximate nature of z2q() we could get a wrong
     * quadrilateral.)
     */
    w = sc_z2w(gg->newsc, z, zs[vidmin], ws[nz * qid + vidmin], gg->newAs[qid], ZZERO, gg->newAs[qid], ZZERO, gg->newBs[qid], &ws[qid * nz], eps * 10.0, &status);

  know_w:
    if (status < 2) {
        zz = sc_w2z(gg->sc, w, vid, ZZERO, gg->As[qid], -1, gg->Bs[qid], &ws[qid * nz]);
        if (status != 0)
            (void) z2q_simple(gg, zz, gg->zs, eps);
        output_point(gg, creal(zz), cimag(zz));
    } else {
        output_point(gg, NaN, NaN);
    }
    if (gg_verbose && gg->out != stdout) {
        if (gg_verbose > 1)
            fprintf(stderr, "(%.10g,%.10g)(%d)->(%.10f,%.10f)->(%.10g,%.10g)\n", creal(z), cimag(z), qid, creal(w), cimag(w), creal(zz), cimag(zz));
        else {
            if (status == 0)
                fprintf(stderr, ".");
            else if (status == 1)
                fprintf(stderr, "o");
            else
                fprintf(stderr, "-");
        }
        fflush(stderr);
    }
}

static void align_onboundaries(gridgen* gg, zdouble* z)
{
    double eps = gg->eps;
    double x = creal(z[0]);
    double y = cimag(z[0]);

    if (fabs(x - gg->newxmin) < eps)
        x = gg->newxmin;
    else if (fabs(x - gg->newxmax) < eps)
        x = gg->newxmax;
    if (fabs(y - gg->newymin) < eps)
        y = gg->newymin;
    else if (fabs(y - gg->newymax) < eps)
        y = gg->newymax;

    z[0] = x + I * y;
}

static void map_quadrilaterals(gridgen* gg)
{
    int nppe = gg->nppe;
    int nq = gg->nquadrilaterals;

    /*
     * (there is 1-to-1 correspondence between inner edges of the
     * triangulation and the diagonals of the quadrilaterals) 
     */
    zdouble* zm_all = malloc(nppe * nq * sizeof(zdouble));

    int maxdiff = gg->nz - 1;
    double eps = gg->eps;
    int nz = gg->nz;
    zdouble* zs = gg->zs;
    zdouble* ws = gg->ws;
    hashtable* e2q = NULL;      /* internal edge to quadrilateral */

    int i, j;

    if (gg_verbose)
        fprintf(stderr, "mapping quadrilaterals (nppe = %d):\n", gg->nppe);
    if (gg_verbose == 1)
        fprintf(stderr, "  ");

    for (i = 0; i < nq; ++i) {
        quadrilateral* q = &gg->quadrilaterals[i];
        zdouble* zm = &zm_all[nppe * i];
        zdouble z0, z1;

        if (gg_verbose > 1)
            fprintf(stderr, "  %d:\n", i);
        /*
         * end points of a diagonal 
         */
        z0 = gg->zs[q->vids[0]];
        z1 = gg->zs[q->vids[2]];

        for (j = 0; j < nppe; ++j) {
            zdouble z = (z0 * (nppe - j) + z1 * (j + 1)) / (nppe + 1);
            int vid = (j < nppe / 2) ? q->vids[0] : q->vids[2];
            int status = 0;

            zdouble w = sc_z2w(gg->sc, z, zs[vid], ws[nz * i + vid], gg->As[i], ZZERO, gg->As[i], ZZERO, gg->Bs[i], &ws[i * nz], eps * 10.0, &status);

            assert(status < 2); /* I can see no reason why the mapping could
                                 * fail here. Bail out. */

            zm[j] = sc_w2z(gg->newsc, w, -1, ZZERO, gg->newAs[i], -1, gg->newBs[i], &ws[i * nz]);

            if (gg_verbose == 1) {
                if (status == 0)
                    fprintf(stderr, ".");
                else if (status == 1)
                    fprintf(stderr, "o");
                else
                    fprintf(stderr, "-");
            }
        }

        if (gg_verbose > 1) {
            z0 = gg->newzs[q->vids[0]];
            z1 = gg->newzs[q->vids[1]];

            fprintf(stderr, "    (%.15g %.15g)\n", creal(z0), cimag(z0));
            for (j = 0; j < nppe; ++j)
                fprintf(stderr, "    (%.15g %.15g)\n", creal(zm[j]), cimag(zm[j]));
            fprintf(stderr, "    (%.15g %.15g)\n", creal(z1), cimag(z1));
        }
    }
    if (gg_verbose == 1)
        fprintf(stderr, "\n");

    /*
     * a bit of muscle building
     */
    e2q = ht_create_i2(gg->nquadrilaterals);
    for (i = 0; i < nq; ++i) {
        quadrilateral* q = &gg->quadrilaterals[i];
        int key[2];

        if (q->vids[0] < q->vids[2]) {
            key[0] = q->vids[0];
            key[1] = q->vids[2];
        } else {
            key[0] = q->vids[2];
            key[1] = q->vids[0];
        }

        /*
         * make sure that the edge generates a unique key
         */
        assert(ht_insert(e2q, &key, q) == NULL);
    }

    /*
     * storing the found images
     */
    gg->nqivertices = calloc(nq, sizeof(int));
    gg->qivertices = malloc(nq * gg->nppq * sizeof(zdouble));
    for (i = 0; i < nq; ++i) {
        quadrilateral* q = &gg->quadrilaterals[i];
        zdouble* qivertices = &gg->qivertices[i * gg->nppq];
        zdouble* zm;
        int edge;

        for (edge = 0; edge <4; ++edge) {
            int vid0 = q->vids[edge];
            int vid1 = q->vids[(edge +1) %4];
            int diff = abs(vid0 - vid1);

            /*
             * store the vertex
             */
            qivertices[gg->nqivertices[i]] = gg->newzs[vid0];
            gg->nqivertices[i]++;

            /*
             * if necessary, store the rest 
             */
            if (diff != 1 && diff != maxdiff) {
                int key[2];
                quadrilateral* q1;
                int qid;
                int j;

                if (vid0 < vid1) {
                    key[0] = vid0;
                    key[1] = vid1;
                } else {
                    key[0] = vid1;
                    key[1] = vid0;
                }
                q1 = ht_find(e2q, &key);
                assert(q1 != NULL);     /* found in the table */
                qid = q1->id;

                zm = &zm_all[nppe * qid];
                if (vid0 == gg->quadrilaterals[qid].vids[0]) {
                    for (j = 0; j < nppe; ++j) {
                        qivertices[gg->nqivertices[i]] = zm[j];
                        gg->nqivertices[i]++;
                    }
                } else {
                    for (j = nppe - 1; j >= 0; --j) {
                        qivertices[gg->nqivertices[i]] = zm[j];
                        gg->nqivertices[i]++;
                    }
                }
            }
        }
        /*
         * close up the polygon
         */
        qivertices[gg->nqivertices[i]] = gg->newzs[q->vids[0]];
        gg->nqivertices[i]++;
    }

    ht_destroy(e2q);
    free(zm_all);
}

static void generate_grid(gridgen* gg)
{
    int count = 0;
    int i, j;

    if (gg_verbose) {
        fprintf(stderr, "generating grid:\n");
        fflush(stderr);
    }

    ode_silent = 1;
    ode_stoponnan = 0;

    if (gg_verbose == 1 && gg->out != stdout)
        fprintf(stderr, "  ");

    if (gg->gridpoints == NULL) {
        int nx = gg->nx;
        int ny = gg->ny;
        double dzx = (gg->newxmax - gg->newxmin) / (nx - 1);
        double dzy = (gg->newymax - gg->newymin) / (ny - 1);
        zdouble z0 = gg->newxmin + I * gg->newymin;

        if (gg->out != NULL)
            fprintf(gg->out, "## %d x %d\n", nx, ny);

        for (j = 0; j < ny; ++j) {
            zdouble zy = z0 + I * dzy * j;

            for (i = 0; i < nx; ++i) {
                zdouble z = zy + dzx * i;

                align_onboundaries(gg, &z);
                if (in_poly(z, gg->ncorners, gg->newrzs, gg->eps)) {
                    map_point(gg, z);
                    count++;
                } else
                    output_point(gg, NaN, NaN);
            }
        }
    } else {
        double dzx = gg->newxmax - gg->newxmin;
        double dzy = gg->newymax - gg->newymin;
        zdouble z0 = gg->newxmin + I * gg->newymin;
        int n = gg->ngridpoints;
        zdouble* zs = (zdouble*) gg->gridpoints;

        for (i = 0; i < n; ++i) {
            zdouble z = z0 + creal(zs[i]) * dzx + I * cimag(zs[i]) * dzy;

            align_onboundaries(gg, &z);
            if (in_poly(z, gg->ncorners, gg->newrzs, gg->eps)) {