/* nag_mesh2d_bound (d06bac) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include /* Structure to allow data to be passed into */ /* the user-supplied function fbnd */ struct user { /* details of the ellipse containing the NAG logo */ double xa, xb, x0, y0; }; #ifdef __cplusplus extern "C" { #endif static double NAG_CALL fbnd(Integer, double, double, Nag_Comm *); #ifdef __cplusplus } #endif #define EDGE(I, J) edge[3*((J) -1)+(I) -1] #define LINED(I, J) lined[4*((J) -1)+(I) -1] #define CONN(I, J) conn[3*((J) -1)+(I) -1] #define COOR(I, J) coor[2*((J) -1)+(I) -1] #define COORCH(I, J) coorch[2*((J) -1)+(I) -1] #define CRUS(I, J) crus[2*((J) -1)+(I) -1] int main(int argc, char *argv[]) { FILE *fpin, *fpout; char *outfile = 0; const Integer sdcrus = 4, nvmax = 1000, nedmx = 300, nvint = 0; struct user ellipse; Nag_Comm comm; double x0, xa, xb, xmax, xmin, y0, ymax, ymin; Integer exit_status, i, itrace, j, k, ncomp, nedge, nelt, nlines; Integer npropa, nv, nvb, reftk, l; char pmesh[2]; double *coor = 0, *coorch = 0, *crus = 0, *rate = 0, *weight = 0; Integer *conn = 0, *edge = 0, *lcomp = 0, *lined = 0, *nlcomp = 0; NagError fail; INIT_FAIL(fail); /* Check for command-line IO options */ fpin = nag_example_file_io(argc, argv, "-data", NULL); fpout = nag_example_file_io(argc, argv, "-results", NULL); (void) nag_example_file_io(argc, argv, "-nag_write", &outfile); exit_status = 0; fprintf(fpout, " nag_mesh2d_bound (d06bac) Example Program Results\n\n"); /* Skip heading in data file */ fscanf(fpin, "%*[^\n] "); /* Initialise boundary mesh inputs: */ /* the number of line and of the characteristic points of */ /* the boundary mesh */ fscanf(fpin, "%ld%*[^\n] ", &nlines); /* Allocate memory */ if (!(coor = NAG_ALLOC(2*nvmax, double)) || !(coorch = NAG_ALLOC(2*nlines, double)) || !(crus = NAG_ALLOC(2*sdcrus, double)) || !(rate = NAG_ALLOC(nlines, double)) || !(weight = NAG_ALLOC(1, double)) || !(conn = NAG_ALLOC(3*(2*nvmax+5), Integer)) || !(edge = NAG_ALLOC(3*nedmx, Integer)) || !(lined = NAG_ALLOC(4*nlines, Integer)) || !(lcomp = NAG_ALLOC(nlines, Integer))) { fprintf(fpout, "Allocation failure\n"); exit_status = -1; goto END; } /* The ellipse boundary which envelops */ /* the NAG Logo, the N, the A and the G */ for (j = 1; j <= nlines; ++j) fscanf(fpin, "%lf", &COORCH(1, j)); fscanf(fpin, "%*[^\n] "); for (j = 1; j <= nlines; ++j) fscanf(fpin, "%lf", &COORCH(2, j)); fscanf(fpin, "%*[^\n] "); for (j = 1; j <= sdcrus; ++j) fscanf(fpin, "%lf", &CRUS(1, j)); fscanf(fpin, "%*[^\n] "); for (j = 1; j <= sdcrus; ++j) fscanf(fpin, "%lf", &CRUS(2, j)); fscanf(fpin, "%*[^\n] "); /* The lines of the boundary mesh */ for (j = 1; j <= nlines; ++j) { for (i = 1; i <= 4; ++i) fscanf(fpin, "%ld", &LINED(i, j)); fscanf(fpin, "%lf", &rate[j-1]); } fscanf(fpin, "%*[^\n] "); /* The number of connected components */ /* to the boundary and their information */ fscanf(fpin, "%ld%*[^\n] ", &ncomp); /* Allocate memory */ if (!(nlcomp = NAG_ALLOC(ncomp, Integer))) { fprintf(fpout, "Allocation failure\n"); exit_status = -1; goto END; } j = 0; for (i = 0; i < ncomp; ++i) { fscanf(fpin, "%ld", &nlcomp[i]); fscanf(fpin, "%*[^\n] "); l = j + abs(nlcomp[i]); for (k = j; k < l; ++k) fscanf(fpin, "%ld", &lcomp[k]); fscanf(fpin, "%*[^\n] "); j += abs(nlcomp[i]); } fscanf(fpin, " ' %1s '%*[^\n] ", pmesh); /* Data passed to the user-supplied function */ xmin = COORCH(1, 4); xmax = COORCH(1, 2); ymin = COORCH(2, 1); ymax = COORCH(2, 3); xa = (xmax-xmin)/2.0; xb = (ymax-ymin)/2.0; x0 = (xmin+xmax)/2.0; y0 = (ymin+ymax)/2.0; comm.p = (Pointer)&ellipse; ellipse.xa = xa; ellipse.xb = xb; ellipse.x0 = x0; ellipse.y0 = y0; itrace = -1; /* Call to the boundary mesh generator */ /* nag_mesh2d_bound (d06bac). * Generates a boundary mesh */ if (outfile) fclose(fpout); nag_mesh2d_bound(nlines, coorch, lined, fbnd, crus, sdcrus, rate, ncomp, nlcomp, lcomp, nvmax, nedmx, &nvb, coor, &nedge, edge, itrace, outfile, &comm, &fail); if (outfile && !(fpout = fopen(outfile, "a"))) { exit_status = 2; goto END; } if (fail.code == NE_NOERROR) { if (pmesh[0] == 'N') { fprintf(fpout, " Boundary mesh characteristics\n"); fprintf(fpout, " nvb =%6ld\n", nvb); fprintf(fpout, " nedge =%6ld\n", nedge); } else if (pmesh[0] == 'Y') { /* Output the mesh to view it using the NAG Graphics Library */ fprintf(fpout, " %10ld%10ld\n", nvb, nedge); for (i = 1; i <= nvb; ++i) fprintf(fpout, " %4ld %15.6e %15.6e \n", i, COOR(1, i), COOR(2, i)); for (i = 1; i <= nedge; ++i) fprintf(fpout, " %4ld%4ld%4ld%4ld\n", i, EDGE(1, i), EDGE(2, i), EDGE(3, i)); } else { fprintf(fpout, "Problem with the printing option Y or N\n"); exit_status = -1; goto END; } } else { fprintf(fpout, "Error from nag_mesh2d_bound (d06bac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Initialise mesh control parameters */ itrace = 0; npropa = 1; /* Call to the 2D Delaunay-Voronoi mesh generator */ /* nag_mesh2d_delaunay (d06abc). * Generates a two-dimensional mesh using a Delaunay-Voronoi * process */ if (outfile) fclose(fpout); nag_mesh2d_delaunay(nvb, nvint, nvmax, nedge, edge, &nv, &nelt, coor, conn, weight, npropa, itrace, outfile, &fail); if (outfile && !(fpout = fopen(outfile, "a"))) { exit_status = 2; goto END; } if (fail.code == NE_NOERROR) { if (pmesh[0] == 'N') { fprintf(fpout, " Complete mesh characteristics (Delaunay-Voronoi)\n"); fprintf(fpout, " nv =%6ld\n", nv); fprintf(fpout, " nelt =%6ld\n", nelt); } else if (pmesh[0] == 'Y') { /* Output the mesh to view it using the NAG Graphics Library */ fprintf(fpout, " %10ld%10ld\n", nv, nelt); for (i = 1; i <= nv; ++i) fprintf(fpout, " %15.6e %15.6e \n", COOR(1, i), COOR(2, i)); reftk = 0; for (k = 1; k <= nelt; ++k) fprintf(fpout, " %10ld%10ld%10ld%10ld\n", CONN(1, k), CONN(2, k), CONN(3, k), reftk); } else { fprintf(fpout, "Problem with the printing option Y or N\n"); exit_status = -1; goto END; } } else { fprintf(fpout, "Error from nag_mesh2d_delaunay (d06abc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Call to the 2D Advancing front mesh generator */ /* nag_mesh2d_front (d06acc). * Generates a two-dimensional mesh using an Advancing-front * method */ if (outfile) fclose(fpout); nag_mesh2d_front(nvb, nvint, nvmax, nedge, edge, &nv, &nelt, coor, conn, weight, itrace, outfile, &fail); if (outfile && !(fpout = fopen(outfile, "a"))) { exit_status = 2; goto END; } if (fail.code == NE_NOERROR) { if (pmesh[0] == 'N') { fprintf(fpout, " Complete mesh characteristics (Advancing Front)\n"); fprintf(fpout, " nv =%6ld\n", nv); fprintf(fpout, " nelt =%6ld\n", nelt); } else if (pmesh[0] == 'Y') { /* Output the mesh to view it using the NAG Graphics Library */ fprintf(fpout, " %10ld%10ld\n", nv, nelt); for (i = 1; i <= nv; ++i) fprintf(fpout, " %15.6e %15.6e \n", COOR(1, i), COOR(2, i)); reftk = 0; for (k = 1; k <= nelt; ++k) fprintf(fpout, " %10ld%10ld%10ld%10ld\n", CONN(1, k), CONN(2, k), CONN(3, k), reftk); } else { fprintf(fpout, "Problem with the printing option Y or N\n"); exit_status = -1; goto END; } } else { fprintf(fpout, "Error from nag_mesh2d_front (d06acc).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (fpin != stdin) fclose(fpin); if (fpout != stdout) fclose(fpout); if (coor) NAG_FREE(coor); if (coorch) NAG_FREE(coorch); if (crus) NAG_FREE(crus); if (rate) NAG_FREE(rate); if (weight) NAG_FREE(weight); if (conn) NAG_FREE(conn); if (edge) NAG_FREE(edge); if (lcomp) NAG_FREE(lcomp); if (lined) NAG_FREE(lined); if (nlcomp) NAG_FREE(nlcomp); return exit_status; } static double NAG_CALL fbnd(Integer i, double x, double y, Nag_Comm *pcomm) { double ret_val, d1, d2; double radius2, x0, xa, xb, y0; struct user *ellipse = (struct user *) pcomm->p; xa = ellipse->xa; xb = ellipse->xb; x0 = ellipse->x0; y0 = ellipse->y0; ret_val = 0.0; switch (i) { case 1: /* line 1,2,3, and 4: ellipse centred in (X0,Y0) with */ /* XA and XB as coefficients */ d1 = (x - x0)/xa; d2 = (y - y0)/xb; ret_val = d1*d1 + d2*d2 - 1.0; break; case 2: /* line 24, 27, 33 and 38 are a circle centred in (X0,Y0) */ /* with radius SQRT(RADIUS2) */ x0 = 20.5; y0 = 4.0; radius2 = 4.25; d1 = x - x0; d2 = y - y0; ret_val = d1*d1 + d2*d2 - radius2; break; case 3: x0 = 17.0; y0 = 8.5; radius2 = 5.0; d1 = x - x0; d2 = y - y0; ret_val = d1*d1 + d2*d2 - radius2; break; case 4: x0 = 17.0; y0 = 8.5; radius2 = 5.0; d1 = x - x0; d2 = y - y0; ret_val = d1*d1 + d2*d2 - radius2; break; case 5: x0 = 19.5; y0 = 4.0; radius2 = 1.25; d1 = x - x0; d2 = y - y0; ret_val = d1*d1 + d2*d2 - radius2; break; default: break; } return ret_val; }