/* nag_mesh2d_smooth (d06cac) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. * * Mark 8 revised, 2004 * */ #include #include #include #include #include #include #define EDGE(I, J) edge[3*((J) -1)+(I) -1] #define CONN(I, J) conn[3*((J) -1)+(I) -1] #define COOR(I, J) coor[2*((J) -1)+(I) -1] int main(void) { /* Integer scalar and array declarations */ const Integer nvfix = 0; Integer exit_status = 0; Integer i, imax, imaxm1, ind, itrace, j, jmax, jmaxm1, k, me1, me2, me3, nedge, nelt, nqint, nv, reftk, lstate; Integer *conn = 0, *edge = 0, *numfix = 0, *state = 0; /* Character scalar and array declarations */ char pmesh[2]; /* NAG structures */ NagError fail; /* Double scalar and array declarations */ double delta, hx, hy, pi, dpi, r, rad, sk, theta, x1, x2, x3, y1, y2, y3; double one_draw[1]; double *coor = 0; /* Choose the base generator */ Nag_BaseRNG genid = Nag_Basic; Integer subid = 0; /* Set the seed */ Integer seed[] = { 1762543 }; Integer lseed = 1; /* Initialise the error structure */ INIT_FAIL(fail); /* Get the length of the state array */ lstate = -1; nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" nag_mesh2d_smooth (d06cac) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n] "); /* Read imax and jmax, the number of vertices */ /* in the x and y directions respectively. */ scanf("%ld", &imax); scanf("%ld", &jmax); scanf("%*[^\n] "); /* Read distortion percentage and calculate radius */ /* of distortion neighbourhood so that cross-over */ /* can only occur at 100% or greater. */ scanf("%lf", &delta); scanf("%*[^\n] "); nv = imax*jmax; imaxm1 = imax - 1; jmaxm1 = jmax - 1; nelt = 2*imaxm1*jmaxm1; nedge = 2*(imaxm1 + jmaxm1); /* Allocate memory */ if (!(coor = NAG_ALLOC(2*nv, double)) || !(conn = NAG_ALLOC(3*nelt, Integer)) || !(edge = NAG_ALLOC(3*nedge, Integer)) || !(state = NAG_ALLOC(lstate, Integer)) || !(numfix = NAG_ALLOC(1, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } scanf(" ' %1s '", pmesh); scanf("%*[^\n] "); hx = 1.0/(double) imaxm1; hy = 1.0/(double) jmaxm1; rad = 0.01*delta*(hx > hy?hy:hx)/2.0; pi = 4.0*atan(1.0); /* Initialise the generator to a repeatable sequence */ nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Generate a simple uniform mesh and then distort it */ /* randomly within the distortion neighbourhood of each */ /* node. */ ind = 0; for (j = 1; j <= jmax; ++j) { for (i = 1; i <= imax; ++i) { /* Generate one uniform variate between 0 and rad */ nag_rand_uniform(1, 0.0, rad, state, one_draw, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_uniform (g05sqc).\n%s\n", fail.message); exit_status = 1; goto END; } r = one_draw[0]; dpi = 2.0*pi; /* Generate one uniform variate between 0 and dpi */ nag_rand_uniform(1, 0.0, dpi, state, one_draw, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_uniform (g05sqc).\n%s\n", fail.message); exit_status = 1; goto END; } theta = one_draw[0]; if (i == 1 || i == imax || j == 1 || j == jmax) r = 0.0; k = (j-1)*imax + i; COOR(1, k) = (i-1)*hx + r *cos(theta); COOR(2, k) = (j-1)*hy + r *sin(theta); if (i < imax && j < jmax) { ++ind; CONN(1, ind) = k; CONN(2, ind) = k + 1; CONN(3, ind) = k + imax + 1; ++ind; CONN(1, ind) = k; CONN(2, ind) = k + imax + 1; CONN(3, ind) = k + imax; } } } if (pmesh[0] == 'N') { printf(" The complete distorted mesh characteristics\n"); printf(" nv =%6ld\n", nv); printf(" nelt =%6ld\n\n", nelt); } else if (pmesh[0] == 'Y') { /* Output the mesh to view it using the NAG Graphics Library */ printf(" %10ld%10ld\n", nv, nelt); for (i = 1; i <= nv; ++i) printf(" %15.6e %15.6e \n", COOR(1, i), COOR(2, i)); } else { printf("Problem with the printing option Y or N\n"); } reftk = 0; for (k = 1; k <= nelt; ++k) { me1 = CONN(1, k); me2 = CONN(2, k); me3 = CONN(3, k); x1 = COOR(1, me1); x2 = COOR(1, me2); x3 = COOR(1, me3); y1 = COOR(2, me1); y2 = COOR(2, me2); y3 = COOR(2, me3); sk = 0.5*((x2-x1)*(y3-y1) - (y2-y1)*(x3-x1)); if (sk < 0.0) { printf("Error the surface of the element is negative\n"); printf(" k = %6ld\n", k); printf(" sk = %15.6e\n", sk); exit_status = -1; goto END; } if (pmesh[0] == 'Y') printf(" %10ld%10ld%10ld%10ld\n", CONN(1, k), CONN(2, k), CONN(3, k), reftk); } /* Boundary edges */ ind = 0; for (i = 1; i <= imaxm1; ++i) { ++ind; EDGE(1, ind) = i; EDGE(2, ind) = i + 1; EDGE(3, ind) = 0; } for (i = 1; i <= jmaxm1; ++i) { ++ind; EDGE(1, ind) = i*imax; EDGE(2, ind) = (i+1)*imax; EDGE(3, ind) = 0; } for (i = 1; i <= (imax - 1); ++i) { ++ind; EDGE(1, ind) = imax*jmax - i + 1; EDGE(2, ind) = imax*jmax - i; EDGE(3, ind) = 0; } for (i = 1; i <= jmaxm1; ++i) { ++ind; EDGE(1, ind) = (jmax - i)*imax + 1; EDGE(2, ind) = (jmax - i - 1)*imax + 1; EDGE(3, ind) = 0; } itrace = 1; nqint = 10; /* Call the smoothing routine */ /* nag_mesh2d_smooth (d06cac). * Uses a barycentering technique to smooth a given mesh */ nag_mesh2d_smooth(nv, nelt, nedge, coor, edge, conn, nvfix, numfix, itrace, 0, nqint, &fail); if (fail.code == NE_NOERROR) { if (pmesh[0] == 'N') { printf("\n The complete smoothed mesh characteristics\n"); printf(" nv =%6ld\n", nv); printf(" nelt =%6ld\n", nelt); } else if (pmesh[0] == 'Y') { /* Output the mesh to view it using the NAG Graphics Library */ printf(" %10ld%10ld\n", nv, nelt); for (i = 1; i <= nv; ++i) printf(" %15.6e %15.6e \n", COOR(1, i), COOR(2, i)); reftk = 0; for (k = 1; k <= nelt; ++k) printf( " %10ld%10ld%10ld%10ld\n", CONN(1, k), CONN(2, k), CONN(3, k), reftk); } else { printf("Problem with the printing option Y or N\n"); exit_status = -1; goto END; } } else { printf("Error from nag_mesh2d_smooth (d06cac).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (coor) NAG_FREE(coor); if (conn) NAG_FREE(conn); if (edge) NAG_FREE(edge); if (numfix) NAG_FREE(numfix); if (state) NAG_FREE(state); return exit_status; }