gridgen.c

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

             */
            if (count == 0)
                F(x, f, gg);
            else
                broyden_update(F, n, x, f, w, gg);
        } else {
            /*
             * simple iterations 
             */
            F(x, f, gg);
            for (i = 0; i < n; ++i)
                x[i] -= f[i];
        }

        for (i = 0; i < n; ++i)
            error += f[i] * f[i];
        error = sqrt(error);

        if (error < error_prev)
            nlastbad = 0;
        else
            nlastbad++;

        if (gg_verbose < 2 && gg_verbose + nlastbad > 2) {
            fprintf(stderr, "\n  error[%d] = %.3g, error[%d] = %.3g", count, error_prev, count + 1, error);
            fflush(stderr);
        } else if (gg_verbose > 1) {
            fprintf(stderr, "  %d: error = %.3g\n", count, error);
            fflush(stderr);
        } else if (gg_verbose) {
            fprintf(stderr, ".");
            fflush(stderr);
        }

        count++;
    } while (error > gg->eps && (!gg->newton || error_prev > gg->eps));

    if (gg_verbose)
        fprintf(stderr, "\n");
    if (gg_verbose > 1) {
        int i;

        fprintf(stderr, "  sigmas:\n");
        for (i = 0; i < n; ++i)
            fprintf(stderr, "  %d: %.15g\n", i, x[i]);
        fflush(stderr);
    }

    if (count > 1 && gg->sigmafname != NULL) {
        if (gg->fsigma != NULL)
            fclose(gg->fsigma);
        gg->fsigma = gg_fopen(gg->sigmafname, "w");
        if (gg_verbose)
            fprintf(stderr, "  saving sigmas to \"%s\":\n", gg->sigmafname);
        if ((int) fwrite(x, sizeof(double), n, gg->fsigma) != n)
            quit("could not save sigmas to \"%s\": %s\n", gg->sigmafname, strerror(errno));
        fclose(gg->fsigma);
        gg->fsigma = NULL;
    }

    if (w != NULL)
        free(w);
    free(f);

    if (gg->sigmas != NULL) {
        if (*gg->sigmas != NULL)
            free(*gg->sigmas);
        *gg->nsigmas = n;
        *gg->sigmas = x;
    } else
        free(x);
}

static void quadrilaterals_destroy(int n, quadrilateral* quads)
{
    int i;

    if (quads == NULL)
        return;
    for (i = 0; i < n; ++i)
        free(quads[i].neighbours);
    free(quads);
}

static void gridgen_destroy(gridgen* gg)
{
    if (gg == NULL)
        return;
    if (gg->zs != NULL)
        free(gg->zs);
    if (gg->ws != NULL)
        free(gg->ws);
    if (gg->As != NULL)
        free(gg->As);
    if (gg->Bs != NULL)
        free(gg->Bs);
    sc_destroy(gg->sc);
    if (gg->newbetas != NULL)
        free(gg->newbetas);
    if (gg->betas != NULL)
        free(gg->betas);
    quadrilaterals_destroy(gg->nquadrilaterals, gg->quadrilaterals);
    if (gg->cis != NULL)
        free(gg->cis);
    delaunay_destroy(gg->d);
    vertlist_destroy(gg->vertices);
    if (gg->datafname != NULL)
        free(gg->datafname);
    sc_destroy(gg->newsc);
    if (gg->newzs != NULL)
        free(gg->newzs);
    if (gg->newAs != NULL)
        free(gg->newAs);
    if (gg->Bs != NULL)
        free(gg->newBs);
    if (gg->gridpoints != NULL)
        free(gg->gridpoints);
    if (gg->out != stdout && gg->out != NULL)
        fclose(gg->out);
    if (gg->sigmafname != NULL)
        free(gg->sigmafname);
    if (gg->fsigma != NULL)
        fclose(gg->fsigma);
    if (gg->rectfname != NULL)
        free(gg->rectfname);
    if (gg->frect != NULL)
        fclose(gg->frect);
    if (gg->rzs != NULL)
        free(gg->rzs);
    if (gg->rvids != NULL)
        free(gg->rvids);
    if (gg->newrzs != NULL)
        free(gg->newrzs);
    if (gg->nqivertices != NULL) {
        free(gg->nqivertices);
        free(gg->qivertices);
    }
    free(gg);
}

static void gridgen_check(gridgen* gg)
{
    int n = gg->vertices->n;
    double* betas = gg->betas;
    double sum = 0.0;
    int i;

    if (gg_verbose) {
        fprintf(stderr, "checking input:\n");
        fflush(stderr);
    }

    for (i = 0; i < n; ++i)
        sum += betas[i];

    if (fabs(sum + 2.0) > EPS)
        quit("gridgen_check(): sum of betas = %.15f (should be 2)\n", sum);

    if (gg_verbose > 1)
        fprintf(stderr, "  sum of betas = 2.0 + %.3g\n", sum + 2.0);
}

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

    if (gg_verbose) {
        fprintf(stderr, "getting image region corner vertex indices:\n");
        fflush(stderr);
    }

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

        if (p->z != 0.0 && !isnan(p->z)) {
            gg->rvids[count] = i;
            count++;
        }
        now = now->next;
    }

    if (gg_verbose > 1) {
        zdouble* rzs = gg->rzs;

        for (i = 0; i < gg->ncorners; ++i)
            fprintf(stderr, "  %d: (%.10g,%.10g) %d\n", i, creal(rzs[i]), cimag(rzs[i]), gg->rvids[i]);
        fflush(stderr);
    }
}

static void set_newbetas(gridgen* gg)
{
    point* points = gg->d->points;
    int n = gg->ncorners;
    int i;

    if (gg_verbose) {
        fprintf(stderr, "setting new betas:\n");
        fflush(stderr);
    }

    gg->newbetas = calloc(gg->d->npoints, sizeof(double));

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

        if (points[index].z != 0)
            gg->newbetas[index] = -points[index].z * 0.5;
    }
}

static int get_first_quadrilateral(gridgen* gg)
{
    int qid = -1;
    istack* qids = istack_create();
    int n = gg->nquadrilaterals;
    int vid = gg->rvids[0];
    int vid1 = (vid + 1) % gg->d->npoints;
    int i;

    for (i = 0; i < n; ++i) {
        int* vids = gg->quadrilaterals[i].vids;

        if (vid == vids[0] || vid == vids[1] || vid == vids[2] || vid == vids[3])
            istack_push(qids, i);
    }

    assert(qids->n > 0);

    if (qids->n == 1)
        qid = qids->v[0];
    else
        for (i = 0; i < qids->n; ++i) {
            int* vids = gg->quadrilaterals[qids->v[i]].vids;

            if (vid1 == vids[0] || vid1 == vids[1] || vid1 == vids[2] || vid1 == vids[3]) {
                qid = qids->v[i];
                break;
            }
        }

    assert(qid >= 0);

    istack_destroy(qids);

    return qid;
}

/* Finds index (0..3) of a vertex within quadrilateral "newqid" not 
 * contained in the quadrilateral "qid".
 */
static int find_newvertex(int nq, quadrilateral qs[], int qid, int pqid)
{
    int* pvids = qs[pqid].vids;
    int* vids = qs[qid].vids;
    int i, j;

    for (i = 0; i < 4; ++i) {
        int vid = vids[i];

        for (j = 0; j < 4; ++j)
            if (vid == pvids[j])
                break;
        if (j == 4)
            return vid;
    }
    return -1;
}

static void calculate_zs(int nq, quadrilateral qs[], int nz, zdouble zs[], int qindex, int vid, swcr* sc, zdouble w[], int* count, int done[], zdouble A[], zdouble B[])
{
    quadrilateral* q = &qs[qindex];
    int* vids = q->vids;
    zdouble z0 = zs[vids[0]];
    zdouble z2 = zs[vids[2]];
    int i;

    if (*count) {
        zdouble dzeta0 = integrate(sc, vids[0], &w[qindex * nz]);
        zdouble dzeta2 = integrate(sc, vids[2], &w[qindex * nz]);
        zdouble dzeta = integrate(sc, vid, &w[qindex * nz]);

        B[qindex] = (z2 - z0) / (dzeta2 - dzeta0);
        A[qindex] = z0 - B[qindex] * dzeta0;

        zs[vid] = A[qindex] + B[qindex] * dzeta;
    }

    (*count)++;
    done[qindex] = 1;

    for (i = 0; i < q->nneighbours; ++i) {
        int nid = q->neighbours[i];
	int vid;

        if (done[nid])
            continue;

	vid = find_newvertex(nq, qs, nid, qindex);
	calculate_zs(nq, qs, nz, zs, nid, vid, sc, w, count, done, A, B);
    }
}

static void find_zminmax(int n, zdouble zs[], double* xmin, double* xmax, double* ymin, double* ymax)
{
    int i;

    *xmin = DBL_MAX;
    *xmax = -DBL_MAX;
    *ymin = DBL_MAX;
    *ymax = -DBL_MAX;

    for (i = 0; i < n; ++i) {
        zdouble z = zs[i];
        double x = creal(z);
        double y = cimag(z);

        if (x < *xmin)
            *xmin = x;
        if (x > *xmax)
            *xmax = x;
        if (y < *ymin)
            *ymin = y;
        if (y > *ymax)
            *ymax = y;
    }
}

static void align_newzs(int n, zdouble rzs[], int rvids[], int nz, zdouble zs[], double betas[], double eps)
{
    double eps5 = eps * 5.0;
    double beta;
    int i;

    for (i = 0; i < n; ++i) {
        if (fabs(creal(rzs[i])) < eps5)
            rzs[i] = I * cimag(rzs[i]);
        else if (fabs(1.0 - creal(rzs[i])) < eps5)
            rzs[i] = 1.0 + I * cimag(rzs[i]);
        if (fabs(cimag(rzs[i])) < eps5)
            rzs[i] = creal(rzs[i]);
        else if (fabs(1.0 - cimag(rzs[i])) < eps5)
            rzs[i] = creal(rzs[i]) + I * 1.0;
    }

    for (i = 0; i < n; ++i)
        zs[rvids[i]] = rzs[i];

    for (i = 0, beta = 0.0; i < n; ++i) {
        int vid0 = rvids[i];
        int vid1 = rvids[(i + 1) % n];
        zdouble z = zs[vid0];
        int j, jj;

        if (fabs(fmod(beta, 0.5)) < eps) {
            int k = rint(fabs(beta / 0.5));

            if (k % 2 == 0) {
                double x = creal(z);

                for (j = 0, jj = (vid0 + 1) % nz; j < (vid1 + nz - vid0) % nz; ++j, jj = (jj + 1) % nz)
                    zs[jj] = x + I * cimag(zs[jj]);
            } else {
                double y = cimag(z);

                for (j = 0, jj = (vid0 + 1) % nz; j < (vid1 + nz - vid0) % nz; ++j, jj = (jj + 1) % nz)
                    zs[jj] = creal(zs[jj]) + I * y;
            }
        }
        beta += betas[vid1];
    }

    for (i = 0; i < n; ++i)
        rzs[i] = zs[rvids[i]];
}

static void write_zarray(int n, zdouble* z, FILE* f)
{
    int i;

    for (i = 0; i < n; ++i)
        fprintf(f, "%f %f\n", creal(z[i]), cimag(z[i]));
    fflush(f);
}

static void calculate_newzs(gridgen* gg)
{
    int nz = gg->nz;
    int nq = gg->nquadrilaterals;
    int* done = calloc(nq, sizeof(int));
    int count = 0;
    int vid0 = gg->rvids[0];
    int qid0;
    quadrilateral* q0 = NULL;
    int i;

    if (gg_verbose) {
        fprintf(stderr, "calculating image region:\n");
        fflush(stderr);
    }

    gg->newzs = calloc(nz, sizeof(zdouble));
    gg->newAs = malloc(nq * sizeof(zdouble));
    gg->newBs = malloc(nq * sizeof(zdouble));

    qid0 = get_first_quadrilateral(gg);
    q0 = &gg->quadrilaterals[qid0];
    if (gg_verbose > 1) {
        int* vids = q0->vids;

        fprintf(stderr, "  start from quadrilateral %d (%d,%d,%d,%d)\n", qid0, vids[0], vids[1], vids[2], vids[3]);
    }

    /*
     * startup 
     */
    {
        zdouble* A = gg->newAs;
        zdouble* B = gg->newBs;
        int* vids = q0->vids;
        int vid1 = (vid0 + 1) % nz;
        int qvid0 = -1;
        int qvid1 = -1;
        zdouble dzetas[4];

        for (i = 0; i < 4; ++i) {
            int vid = vids[i];

            dzetas[i] = integrate(gg->newsc, vid, &gg->ws[qid0 * nz]);