swcr.c

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

/******************************************************************************
 *
 * File:           swcr.c
 *
 * Created:        19/03/2001
 *
 * Author:         Pavel Sakov
 *                 CSIRO Marine Research
 *
 * Purpose:        Calculation of Schwarz-Christoffel transform
 *
 * Description:    Calculates Schwarz-Christoffel transform as described in:
 *                 Lloyd N. Trefethen, "Numerical Computation of the
 *                 Schwartz-Christoffel transformation", SIAM J. Sci. Stat.
 *                 Comput., 1980, 1(1), 82-102.
 *                 A fair amount of mission-critical code has been borrowed 
 *                 from SCPACK (ported from FORTRAN), 
 *                 see http://www.netlib.org/conformal.
 *
 * Revisions:      None
 *
 *****************************************************************************/

#include "config.h"
#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include "zode.h"
#include "swcr.h"
#include "nan.h"

#define NEWTON_NITER_MAX 20
#define INTEGRATION_NINT_MAX 1000
#define WABS_OK_MAX 1.1
#define WABS_NEWTON_MAX 2.0
#define ZZERO 0.0 + 0.0 * I

#define min(x,y) ((x) < (y) ? (x) : (y))

#if HAVE_TGAMMA
double tgamma(double);
#elif HAVE_LGAMMA
double lgamma(double);
#endif

struct swcr {
    int n;
    int nq;
    double* betas;
    double* nodes;
    double* weights;
    double* work;
    int singular;               /* flag */
};

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

    va_start(args, format);
    vfprintf(stderr, format, args);
    va_end(args);
    exit(1);
}

/* Borrowed from "sclibdbl.f" (SCPACK).
 * This procedure supplies the coefficients a(j), b(j) of the recurrence
 * relation
 *
 *   b[j] p[j](x) = (x - a[j]) p[j-1](x) - b[j-1] p[j-2](x)
 *
 * for the various classical (normalized) orthogonal polynomials, and the
 * zero-th moment
 *
 *   muzero = integral w(x) dx
 *
 * of the given polynomial's weight function w(x). Since the polynomials are
 * ortho-normalized, the tridiagonal matrix is guaranteed to be symmetric.
 */
static double class(int n, double alpha, double beta, double* b, double* a)
{
    int nm1 = n - 1;
    double ab = alpha + beta;
    double abi = ab + 2.0;
    double abi2;
    double bbaa = beta * beta - alpha * alpha;
    double a1 = alpha + 1;
    double b1 = beta + 1;

#if HAVE_TGAMMA
    double muzero = pow(2.0, ab + 1.0) * tgamma(a1) * tgamma(b1) / tgamma(abi);
#elif HAVE_LGAMMA
    double muzero = pow(2.0, ab + 1.0) * exp(lgamma(a1) + lgamma(b1) - lgamma(abi));
#endif
    int i;

    a[0] = (beta - alpha) / abi;
    b[0] = sqrt(4.0 * a1 * b1 / (ab + 3.0) / abi / abi);

    for (i = 1; i < nm1; ++i) {
        int i1 = i + 1;

        abi = 2.0 * i1 + ab;
        abi2 = abi * abi;
        a[i] = bbaa / (abi - 2.0) / abi;
        b[i] = sqrt(4.0 * i1 * (alpha + i1) * (beta + i1) * (ab + i1) / (abi2 - 1.0) / abi2);
    }

    abi += 2.0;
    a[nm1] = bbaa / (abi - 2.0) / abi;
    b[nm1] = 0.0;

    return muzero;
}

/* Borrowed from "sclibdbl.f" (SCPACK).
 * This is a modified version of the eispack routine imtql2(). It finds the
 * eigenvalues and first components of the eigenvectors of a symmetric
 * tridiagonal matrix by the implicit ql method.
 */
static void imtql2(int n, double* d, double* e, double* z)
{
    int nm1 = n - 1;
    double b, c, f, g, p, r, s;
    int i, l;

    if (n <= 1)
        return;

    e[nm1] = 0.0;

    for (l = 0; l < n; ++l) {
        int j = 0;
        int m, mml;

        while (1) {
            for (m = l; m < nm1; ++m)
                if (fabs(e[m]) < UROUND * (fabs(d[m]) + fabs(d[m + 1])))
                    break;
            p = d[l];

            if (m == l)
                goto continue_l;

            if (j == 30)
                quit("error: imtql2(): number of iterations > 30\n");

            j++;

            g = (d[l + 1] - p) / e[l] / 2.0;
            r = sqrt(g * g + 1.0);
            g = d[m] - p + e[l] / (g + copysign(r, g));
            s = 1.0;
            c = 1.0;
            p = 0.0;
            mml = m - l;

            for (i = m - 1; i >= l; --i) {
                f = s * e[i];
                b = c * e[i];
                if (fabs(f) >= fabs(g)) {
                    c = g / f;
                    r = sqrt(c * c + 1.0);
                    e[i + 1] = f * r;
                    s = 1.0 / r;
                    c *= s;
                } else {
                    s = f / g;
                    r = sqrt(s * s + 1.0);
                    e[i + 1] = g * r;
                    c = 1.0 / r;
                    s *= c;
                }
                g = d[i + 1] - p;
                r = (d[i] - g) * s + 2.0 * c * b;
                p = s * r;
                d[i + 1] = g + p;
                g = c * r - b;
                f = z[i + 1];
                z[i + 1] = s * z[i] + c * f;
                z[i] = c * z[i] - s * f;
            }
            d[l] -= p;
            e[l] = g;
            e[m] = 0.0;
        }                       /* while (1) */
      continue_l:
        ;
    }

    for (i = 0; i < nm1; ++i) {
        int k = i;
        int j;

        p = d[i];

        for (j = i + 1; j < n; ++j) {
            if (d[j] < p) {
                k = j;
                p = d[j];
            }
        }

        if (k != i) {
            d[k] = d[i];
            d[i] = p;
            p = z[i];
            z[i] = z[k];
            z[k] = p;
        }
    }
}

/* Borrowed from "sclibdbl.f" (SCPACK).
 * This routine computes the nodes t[] and weights w[] for Gauss-Jacobi
 * quadrature formulas. These are used when one wishes to approximate
 *
 *   integral (from a to b) f(x) w(x) dx
 *
 * by
 *
 *    n
 *   sum w[j] f(t[j])
 *   j=1
 *
 * (w(x) and w[j] have no connection with each other), where w(x) is the weight
 * function,
 *
 *   w(x) = (1-x)**alpha * (1+x)**beta
 *
 * on (-1, 1); alpha, beta > -1.
 *
 * @param n     The number of points used for the quadrature rule
 * @param alpha Parameter of the weight function
 * @param beta  Parameter of the weight function
 * @param b     Work array [n]
 * @param t     Nodes [n] (output)
 * @param w     Weights [n] (output)
 *
 * Subroutines required: class, imtql2
 *
 * Reference:
 *
 * The routine has been adapted from the more general routine gaussq by Golub
 * and Welsch. See Golub, G. H., and Welsch, J. H., "Calculation of Gaussian
 * quadrature rules", Mathematics of computation, 23 (April 1969), pp. 221-230.
 */
static void gaussj(int n, double alpha, double beta, double* b, double* t, double* w)
{
    double mu0 = class(n, alpha, beta, b, t);
    int i;

    w[0] = 1.0;
    imtql2(n, t, b, w);

    for (i = 0; i < n; ++i)
        w[i] = mu0 * w[i] * w[i];
}

/* Creates Schwarz-Christoffel transform structure.
 * @param n Number of prevertices
 * @param nq Number of nodes in Gauss-Jacobi quadrature formulas
 * @param betas Array of angle parameters [n]: betas[i] = alpha[i]/pi - 1, 
 *              where alpha[i] -- interior angle at i-th vertex
 * @return Pointer to Schwarz-Christoffel transform structure
 */
swcr* sc_create(int n, int nq, double betas[])
{
    swcr* sc;
    int i, ii;

    if (nq <= 0)
        return NULL;

    sc = malloc(sizeof(swcr));

    sc->n = n;
    sc->nq = nq;
    sc->betas = betas;
    sc->nodes = malloc((n + 1) * nq * sizeof(double));
    sc->weights = calloc((n + 1) * nq, sizeof(double));
    sc->work = malloc(nq * sizeof(double));
    sc->singular = 0;

    for (i = 0; i < n; ++i) {
        ii = nq * i;
        if (betas[i] > -1.0)
            gaussj(nq, 0.0, betas[i], sc->work, &sc->nodes[ii], &sc->weights[ii]);
    }
    ii = nq * n;
    gaussj(nq, 0.0, 0.0, sc->work, &sc->nodes[ii], &sc->weights[ii]);

    return sc;
}

/* Destroys Schwarz-Christoffel transform structure.
 * @param sc Structure to be destroyed
 */
void sc_destroy(swcr* sc)
{
    if (sc == NULL)
        return;
    free(sc->nodes);