gridgen.c

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

    if (gg_verbose)
        fprintf(stderr, "  image region: %d corners\n", count);

    if (count < 3)
        quit("less than 3 corners in the image region\n", count);
    if (fabs(fabs(sum) - 4.0) > EPS)
        quit("sum of beta_i = %.15f != 4\n", sum);

    /*
     * normalize sum of image region corner angles if a small discrepancy
     * exists 
     */
    mult = sum / 4.0;
    if (mult != 1.0)
        for (i = 0; i < n; ++i) {
            now->p.z *= mult;
            now = now->next;
        }

    gg->ncorners = count;
    gg->rzs = malloc(sizeof(zdouble) * count);
    gg->rvids = malloc(sizeof(int) * count);
    gg->newrzs = malloc(sizeof(zdouble) * count);

    /*
     * some extra vertices may be inserted, so we will get `rvids' later 
     */
    for (i = 0; i < gg->ncorners; ++i)
        gg->rvids[i] = -1;

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

        if (p->z != 0.0) {
            gg->rzs[count] = p->x + I * p->y;
            count++;
        }
        now = now->next;
    }

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

static double* calculate_betas(vertlist* l)
{
    int n = l->n;
    double* betas = malloc(n * sizeof(double));
    vertnode* now = l->first;
    int i;

    point* pprev = &now->prev->p;
    point* pnow = &now->p;
    point* pnext = &now->next->p;

    double dx1 = pnow->x - pprev->x;
    double dy1 = pnow->y - pprev->y;
    double dx2 = pnext->x - pnow->x;
    double dy2 = pnext->y - pnow->y;

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

    for (i = 0; i < n; ++i) {

        double arg = carg(-(dx1 + I * dy1) / (dx2 + I * dy2));

        if (arg < 0)
            arg += M_2PI;
        betas[i] = arg / M_PI - 1.0;

        if (gg_verbose > 1)
            fprintf(stderr, "  %d: %f\n", i, betas[i]);

        now = now->next;

        pprev = pnow;
        pnow = pnext;
        pnext = &now->next->p;
        dx1 = dx2;
        dy1 = dy2;
        dx2 = pnext->x - pnow->x;
        dy2 = pnext->y - pnow->y;
    }

    return betas;
}

static quadrilateral* set_quadrilaterals(delaunay* d)
{
    quadrilateral* out = malloc((d->npoints - 3) * sizeof(quadrilateral));
    istack* diags = istack_create();
    istack* neighbours = istack_create();
    int* edges = d->edges;
    int maxdiff = d->npoints - 1;
    int n = 0;
    int i, ii;

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

    for (i = 0, ii = 0; i < d->nedges; ++i, ii += 2) {
        int diff = abs(edges[ii] - edges[ii + 1]);

        if (diff != 1 && diff != maxdiff) {
            istack_push(diags, edges[ii]);
            istack_push(diags, edges[ii + 1]);
            n++;
        }
    }

    if (gg_verbose)
        fprintf(stderr, "  %d diagonals\n", n);

    for (i = 0; i < n; ++i) {
        quadrilateral* q = &out[i];
        int vid1 = istack_pop(diags);
        int vid2 = istack_pop(diags);
        int i1, i2;
        triangle* triangles[2];
        int nt = 0;

        q->id = i;
        for (i1 = 0; i1 < d->n_point_triangles[vid1]; ++i1) {
            int tid1 = d->point_triangles[vid1][i1];

            for (i2 = 0; i2 < d->n_point_triangles[vid2]; ++i2) {
                int tid2 = d->point_triangles[vid2][i2];

                if (tid1 == tid2) {
                    triangles[nt] = &d->triangles[tid1];
                    q->tids[nt] = tid1;
                    nt++;
                }
            }
        }
        q->vids[0] = vid1;
        q->vids[2] = vid2;
        for (ii = 0; ii < 3; ++ii) {
            int vid = triangles[0]->vids[ii];

            if (vid != vid1 && vid != vid2)
                q->vids[1] = vid;
        }
        for (ii = 0; ii < 3; ++ii) {
            int vid = triangles[1]->vids[ii];

            if (vid != vid1 && vid != vid2)
                q->vids[3] = vid;
        }

        for (ii = 0; triangles[0]->vids[ii] != vid1; ++ii);
        if (triangles[0]->vids[(ii + 1) % 3] == vid2) {
            int tmp = q->vids[1];

            q->vids[1] = q->vids[3];
            q->vids[3] = tmp;
            tmp = q->tids[0];
            q->tids[0] = q->tids[1];
            q->tids[1] = tmp;
        }
    }

    /*
     * set neighbours 
     */
    n = d->npoints - 3;
    for (i = 0; i < n; ++i) {
        quadrilateral* q = &out[i];
        int t0 = q->tids[0];
        int t1 = q->tids[1];
        int ii;

        istack_reset(neighbours);

        for (ii = 0; ii < n; ++ii) {
            quadrilateral* qq = &out[ii];
            int tt0 = qq->tids[0];
            int tt1 = qq->tids[1];

            if (ii == i)
                continue;

            if (t0 == tt0 || t0 == tt1 || t1 == tt0 || t1 == tt1)
                istack_push(neighbours, ii);
        }

        q->nneighbours = neighbours->n;
        q->neighbours = malloc(q->nneighbours * sizeof(int));
        memcpy(q->neighbours, neighbours->v, q->nneighbours * sizeof(int));
    }

    if (gg_verbose > 1)
        for (i = 0; i < n; ++i) {
            quadrilateral* q = &out[i];
            int ii;

            fprintf(stderr, "  %d:  (%d,%d,%d,%d) (%d,%d)", i, q->vids[0], q->vids[1], q->vids[2], q->vids[3], q->tids[0], q->tids[1]);
            for (ii = 0; ii < q->nneighbours; ++ii)
                fprintf(stderr, (ii == 0) ? " (%d" : ",%d", q->neighbours[ii]);
            fprintf(stderr, ")\n");
        }

    istack_destroy(neighbours);
    istack_destroy(diags);

    return out;
}

static double* calculate_cis(int nquads, quadrilateral* quads, delaunay* d)
{
    double* cis = malloc(nquads * sizeof(double));
    point* points = d->points;
    int i;

    if (gg_verbose) {
        fprintf(stderr, "calculating log|ro|:\n");
        fflush(stderr);
    }

    for (i = 0; i < nquads; ++i) {
        quadrilateral* q = &quads[i];
        point* p1 = &points[q->vids[0]];
        point* p2 = &points[q->vids[1]];
        point* p3 = &points[q->vids[2]];
        point* p4 = &points[q->vids[3]];
        zdouble z1 = p1->x + I * p1->y;
        zdouble z2 = p2->x + I * p2->y;
        zdouble z3 = p3->x + I * p3->y;
        zdouble z4 = p4->x + I * p4->y;

        cis[i] = log(cabs((z4 - z1) * (z2 - z3) / (z3 - z4) / (z1 - z2)));
        if (gg_verbose > 1)
            fprintf(stderr, "  %d: %.10g\n", i, cis[i]);
    }

    return cis;
}

static void process_quadri(quadrilateral qs[], int qindex, int do_second_triangle, double x[], zdouble w[], int done[], int* count)
{
    quadrilateral* q = &qs[qindex];
    double ro = -exp(x[qindex]);
    int i;

    if (do_second_triangle) {
        zdouble* pa = &w[q->vids[0]];
        zdouble* pb = &w[q->vids[1]];
        zdouble* pc = &w[q->vids[2]];

        if (*pc != *pb) {
            zdouble h = ro * (*pb - *pa) / (*pc - *pb);

            w[q->vids[3]] = (h * *pc + *pa) / (h + 1.0);
        } else
            /*
             * OK, crowding occurs here. The principal point with CRDT is
             * that it occurs for quadrilaterals far from the starting one
             * for this (local) mapping and can be safely ignored. 
             */
            w[q->vids[3]] = *pc;
    } else {
        zdouble* pa = &w[q->vids[2]];
        zdouble* pb = &w[q->vids[3]];
        zdouble* pc = &w[q->vids[0]];

        if (*pc != *pb) {
            zdouble h = ro * (*pb - *pa) / (*pc - *pb);

            w[q->vids[1]] = (h * *pc + *pa) / (h + 1.0);
        } else
            /*
             * see the comment above 
             */
            w[q->vids[1]] = *pc;
    }
    done[qindex] = 1;
    (*count)++;

    for (i = 0; i < q->nneighbours; ++i) {
        int nid = q->neighbours[i];
        quadrilateral* q0 = &qs[nid];

        if (done[nid])
            continue;
        if (q0->tids[0] == q->tids[0] || q0->tids[0] == q->tids[1])
            process_quadri(qs, nid, 1, x, w, done, count);
        else
            process_quadri(qs, nid, 0, x, w, done, count);
    }
}

static void calculate_fi(int nq, quadrilateral qs[], int qindex, double x[], zdouble w[])
{
    int count = 0;
    int* done = calloc(nq, sizeof(int));

    /*
     * startup 
     */
    double ro = exp(x[qindex]);
    double alpha = atan(sqrt(ro));
    double sin_alpha = sin(alpha);
    double cos_alpha = cos(alpha);
    quadrilateral* q = &qs[qindex];

    w[q->vids[0]] = cos_alpha + I * sin_alpha;
    w[q->vids[1]] = -cos_alpha + I * sin_alpha;
    w[q->vids[2]] = -cos_alpha - I * sin_alpha;

    /*
     * go 
     */
    process_quadri(qs, qindex, 1, x, w, done, &count);  /* recursive */

    free(done);
}

static zdouble integrate(swcr* sc, int vid, zdouble w[])
{
    return sc_w2z(sc, w[vid], vid, ZZERO, ZZERO, -1, 1.0 + 0.0 * I, w);
}

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

    calculate_fi(nq, qs, qindex, x, w);

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

        dzetas[i] = integrate(sc, vid, w);
    }

    *B = (z2 - z0) / (dzetas[2] - dzetas[0]);
    *A = z0 - *B * dzetas[0];
}

static void F(double x[], double f[], void* p)
{
    gridgen* gg = (gridgen*) p;
    int nq = gg->nquadrilaterals;
    int nz = gg->nz;
    zdouble dzetas[4];
    int i;

    for (i = 0; i < nq; ++i) {
        zdouble dzeta0, dzeta1, dzeta2, dzeta3;

        calculate_dzetas(nq, gg->quadrilaterals, nz, gg->zs, i, gg->sc, x, &gg->ws[i * nz], &gg->As[i], &gg->Bs[i], dzetas);

        dzeta0 = dzetas[0];
        dzeta1 = dzetas[1];
        dzeta2 = dzetas[2];
        dzeta3 = dzetas[3];
        f[i] = log(cabs((dzeta3 - dzeta0) * (dzeta1 - dzeta2) / (dzeta2 - dzeta3) / (dzeta0 - dzeta1))) - gg->cis[i];
        if (isnan(f[i])) {
            if (gg_verbose && sc_issingular(gg->sc)) {
                fprintf(stderr, "  You have encountered a (very rare) case when a pre-vertice in an arising\n");
                fprintf(stderr, "  Schwarz-Christoffel integral coincides with the point of crowding. You may\n");
                fprintf(stderr, "  try to overcome this by either:\n");
                fprintf(stderr, "    -- deleting the file with the initial sigma values (the one pointed by\n");
                fprintf(stderr, "       \"sigmas\" entry in the parameter file);\n");
                fprintf(stderr, "    -- changing \"newton 1\" to \"newton 0\" or vice versa in the parameter file;\n");
                fprintf(stderr, "    -- slightly moving the boundary polygon vertices.\n");
            }
            quit("F(): NaN detected\n");
        }
    }
}

static void find_sigmas(gridgen* gg, func F)
{
    int n = gg->nquadrilaterals;
    double* x = malloc(n * sizeof(double));
    double* f = malloc(n * sizeof(double));
    double* w = NULL;
    double error = DBL_MAX;
    double error_prev;
    int count = 0;
    int nlastbad = 0;

    if (gg_verbose) {
        fprintf(stderr, "solving for sigmas:\n");
        fflush(stderr);
    }

    if (gg->sigmas != NULL && *gg->nsigmas == n && gg->sigmas != NULL && *gg->sigmas != NULL)
        memcpy(x, *gg->sigmas, n * sizeof(double));
    else if (gg->fsigma == NULL || (int) fread(x, sizeof(double), n, gg->fsigma) != n)
        memcpy(x, gg->cis, n * sizeof(double));

    gg->ws = calloc(n * gg->vertices->n, sizeof(zdouble));
    gg->As = calloc(gg->vertices->n, sizeof(zdouble));
    gg->Bs = calloc(gg->vertices->n, sizeof(zdouble));

    if (gg->newton == 1) {
        double* ww;
        int i;

        w = calloc(n * n, sizeof(double));

        for (i = 0, ww = w; i < n; ++i, ww += n + 1)
            ww[0] = -1.0;
    } else if (gg->newton > 1)
        w = calloc(n * M * 2, sizeof(double));

    if (gg_verbose == 1)
        fprintf(stderr, "  ");

    do {
        int i;

        error_prev = error;
        error = 0.0;

        if (gg->newton) {
            /*
             * newton method with broyden update