/* nag_mesh2d_smooth (d06cac) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #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(int argc, char* argv[]) { const Integer nvfix=0; const double zero=0.0; double delta, hx, hy, pi, dpi, r, rad, sk, theta, x1, x2, x3, y1, y2, y3; Integer exit_status, i, imax, imaxm1, ind, itrace, j, jmax, jmaxm1, k, me1, me2, me3, nedge, nelt, nqint, nv, reftk; char pmesh[2]; double *coor=0; Integer *conn=0, *edge=0, *numfix=0; NagError fail; INIT_FAIL(fail); exit_status = 0; Vprintf(" d06cac Example Program Results\n\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); /* Read imax and jmax, the number of vertices */ /* in the x and y directions respectively. */ Vscanf("%ld", &imax); Vscanf("%ld", &jmax); Vscanf("%*[^\n] "); /* Read distortion percentage and calculate radius */ /* of distortion neighbourhood so that cross-over */ /* can only occur at 100% or greater. */ Vscanf("%lf", &delta); Vscanf("%*[^\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)) || !(numfix = NAG_ALLOC(1, Integer)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } Vscanf(" ' %1s '", pmesh); Vscanf("%*[^\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); g05cbc(zero); /* 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) { r = g05dac(zero, rad); dpi = 2.0*pi; theta = g05dac(zero, dpi); 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') { Vprintf(" The complete distorted mesh characteristics\n"); Vprintf(" nv =%6ld\n", nv); Vprintf(" nelt =%6ld\n", nelt); } else if (pmesh[0] == 'Y') { /* Output the mesh to view it using the NAG Graphics Library */ Vprintf(" %10ld%10ld\n", nv, nelt); for (i = 1; i <= nv; ++i) Vprintf(" %12.6e %12.6e \n", COOR(1,i), COOR(2,i)); } else { Vprintf("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) { Vprintf("Error the surface of the element is negative\n"); Vprintf(" k = %6ld\n", k); Vprintf(" sk = %12.6e\n", sk); exit_status = -1; goto END; } if (pmesh[0] == 'Y') Vprintf(" %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 */ d06cac(nv, nelt, nedge, coor, edge, conn, nvfix, numfix, itrace, 0, nqint, &fail); if (fail.code == NE_NOERROR) { if (pmesh[0] == 'N') { Vprintf(" The complete smoothed mesh characteristics\n"); Vprintf(" nv =%6ld\n", nv); Vprintf(" nelt =%6ld\n", nelt); } else if (pmesh[0] == 'Y') { /* Output the mesh to view it using the NAG Graphics Library */ Vprintf(" %10ld%10ld\n", nv, nelt); for (i = 1; i <= nv; ++i) Vprintf(" %12.6e %12.6e \n", COOR(1,i), COOR(2,i)); reftk = 0; for (k = 1; k <= nelt; ++k) Vprintf(" %10ld%10ld%10ld%10ld\n", CONN(1,k), CONN(2,k), CONN(3,k), reftk); } else { Vprintf("Problem with the printing option Y or N\n"); exit_status = -1; goto END; } } else { Vprintf("Error from 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); return exit_status; }