/* nag_mesh2d_bound (d06bac) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #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; }; static double fbnd(Integer , double , double , Nag_Comm *); #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 COORUS(I,J) coorus[2*((J)-1)+(I)-1] int main(int argc, char* argv[]) { const Integer nus=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, npropa, nv, nvb, reftk, l; NagError fail; char pmesh[2]; double *coor=0, *coorch=0, *coorus=0, *rate=0, *weight=0; Integer *conn=0, *edge=0, *lcomp=0, *lined=0, *nlcomp=0; INIT_FAIL(fail); exit_status = 0; Vprintf(" d06bac Example Program Results\n\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); /* Initialise boundary mesh inputs: */ /* the number of line and of the characteristic points of */ /* the boundary mesh */ Vscanf("%ld%*[^\n] ", &nlines); /* Allocate memory */ if ( !(coor = NAG_ALLOC(2*nvmax, double)) || !(coorch = NAG_ALLOC(2*nlines, double)) || !(coorus = NAG_ALLOC(2*nus, 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)) ) { Vprintf("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) Vscanf("%lf", &COORCH(1,j)); Vscanf("%*[^\n] "); for (j = 1; j <= nlines; ++j) Vscanf("%lf", &COORCH(2,j)); Vscanf("%*[^\n] "); for (j = 1; j <= nus; ++j) Vscanf("%lf", &COORUS(1,j)); Vscanf("%*[^\n] "); for (j = 1; j <= nus; ++j) Vscanf("%lf", &COORUS(2,j)); Vscanf("%*[^\n] "); /* The lines of the boundary mesh */ for (j = 1; j <= nlines; ++j) { for (i = 1; i <= 4; ++i) Vscanf("%ld", &LINED(i,j)); Vscanf("%lf", &rate[j-1]); } Vscanf("%*[^\n] "); /* The number of connected components */ /* to the boundary and their information */ Vscanf("%ld%*[^\n] ", &ncomp); /* Allocate memory */ if (!(nlcomp = NAG_ALLOC(ncomp, Integer)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } j = 0; for (i = 0; i < ncomp; ++i) { Vscanf("%ld", &nlcomp[i]); Vscanf("%*[^\n] "); l = j + abs(nlcomp[i]); for (k = j; k < l; ++k) Vscanf("%ld", &lcomp[k]); Vscanf("%*[^\n] "); j += abs(nlcomp[i]); } Vscanf(" ' %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 */ d06bac(nlines, coorch, lined, fbnd, coorus, nus, rate, ncomp, nlcomp, lcomp, nvmax, nedmx, &nvb, coor, &nedge, edge, itrace, 0, &comm, &fail); if (fail.code == NE_NOERROR) { if (pmesh[0] == 'N') { Vprintf(" Boundary mesh characteristics\n"); Vprintf(" nvb =%6ld\n", nvb); Vprintf(" nedge =%6ld\n", nedge); } else if (pmesh[0] == 'Y') { /* Output the mesh to view it using the NAG Graphics Library */ Vprintf(" %10ld%10ld\n", nvb, nedge); for (i = 1; i <= nvb; ++i) Vprintf(" %4ld %12.6e %12.6e \n", i, COOR(1,i), COOR(2,i)); for (i = 1; i <= nedge; ++i) Vprintf(" %4ld%4ld%4ld%4ld\n", i, EDGE(1,i), EDGE(2,i), EDGE(3,i)); } else { Vprintf("Problem with the printing option Y or N\n"); exit_status = -1; goto END; } } else { Vprintf("Error from 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 */ d06abc(nvb, nvint, nvmax, nedge, edge, &nv, &nelt, coor, conn, weight, npropa, itrace, 0, &fail); if (fail.code == NE_NOERROR) { if (pmesh[0] == 'N') { Vprintf(" Complete mesh characteristics (Delaunay-Voronoi)\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 d06abc.\n%s\n", fail.message); exit_status = 1; goto END; } /* Call to the 2D Advancing front mesh generator */ d06acc(nvb, nvint, nvmax, nedge, edge, &nv, &nelt, coor, conn, weight, itrace, 0, &fail); if (fail.code == NE_NOERROR) { if (pmesh[0] == 'N') { Vprintf(" Complete mesh characteristics (Advancing Front)\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 d06acc.\n%s\n", fail.message); exit_status = 1; goto END; } END: if (coor) NAG_FREE(coor); if (coorch) NAG_FREE(coorch); if (coorus) NAG_FREE(coorus); 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 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; }