/* nag_mesh2d_front (d06acc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include /* Structure to allow data to be passed onto */ /* the d06bac user-supplied function fbnd */ struct user { /* details of the double NACA0012 and the circle around it */ double x0, y0, x1, y1, radius; }; 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(void) { const Integer nus=1, nvmax=2000, nedmx=200, nvint=40; struct user geom_Naca; double dnvint, radius, x0, x1, y0, y1; Integer exit_status, i, itrace, j, k, l, ncomp, nedge, nelt, nlines, nv, nvb, nvint2, reftk; 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; NagError fail; Nag_Comm comm; INIT_FAIL(fail); exit_status = 0; Vprintf(" d06acc Example Program Results\n\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); /* Initialise boundary mesh inputs: the number of lines and */ /* the number of characteristic points of the boundary mesh */ Vscanf("%ld", &nlines); Vscanf("%*[^\n] "); /* 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(nvint, 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 double NACA0012 and the circle around it */ 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] "); /* 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", &ncomp); Vscanf("%*[^\n] "); /* 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 */ x0 = 1.5; y0 = 0.0; radius = 4.5; x1 = 0.8; y1 = -0.3; comm.p = (Pointer)&geom_Naca; geom_Naca.x0 = x0; geom_Naca.y0 = y0; geom_Naca.radius = radius; geom_Naca.x1 = x1; geom_Naca.y1 = y1; itrace = 0; /* Call to the 2D 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; /* Generation of interior vertices */ /* for the wake of the first NACA */ nvint2 = nvint/2; dnvint = 5.0/(nvint2 + 1.0); for (i = 1; i <= nvint2; ++i) { reftk = nvb + i; COOR(1, reftk) = i*dnvint + 1.0; COOR(2, reftk) = 0.0; weight[i-1] = 0.05; } /* for the wake of the second one */ dnvint = 4.19/(nvint2 + 1.0); for (i = nvint2+1; i <= nvint; ++i) { reftk = nvb + i; COOR(1, reftk) = (i - nvint2)*dnvint + 1.8; COOR(2, reftk) = -0.3; weight[i-1] = 0.05; } /* 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\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; double c, radius, x0, x1, y0, y1; struct user *geom_Naca = (struct user *)pcomm->p; x0 = geom_Naca->x0; y0 = geom_Naca->y0; radius = geom_Naca->radius; x1 = geom_Naca->x1; y1 = geom_Naca->y1; ret_val = 0.0; switch (i) { case 1: /* upper NACA0012 wing beginning at the origin */ c = 1.008930411365; ret_val = 0.6*(0.2969*sqrt(c*x) - 0.126*(c*x) - 0.3516*pow(c*x,2.0) + 0.2843*pow(c*x,3.0) - 0.1015*pow(c*x,4.0)) - c*y; break; case 2: /* lower NACA0012 wing beginning at the origin */ c = 1.008930411365; ret_val = 0.6*(0.2969*sqrt(c*x) - 0.126*(c*x) - 0.3516*pow(c*x,2.0) + 0.2843*pow(c*x,3.0) - 0.1015*pow(c*x,4.0)) + c*y; break; case 3: /* the circle around the double NACA */ ret_val = (x-x0)*(x-x0) + (y-y0)*(y-y0) - radius*radius; break; case 4: /* upper NACA0012 wing beginning at (X1;Y1) */ c = 1.008930411365; ret_val = 0.6*(0.2969*sqrt(c*(x-x1)) - 0.126*c*(x-x1) - 0.3516*pow(c*(x-x1),2.0) + 0.2843*pow(c*(x-x1),3.0) - 0.1015*pow(c*(x-x1),4.0)) - c*(y-y1); break; case 5: /* lower NACA0012 wing beginning at (X1;Y1) */ c = 1.008930411365; ret_val = 0.6*(0.2969*sqrt(c*(x-x1)) - 0.126*(c*(x-x1)) - 0.3516*pow(c*(x-x1),2.0) + 0.2843*pow(c*(x-x1),3.0) - 0.1015*pow(c*(x-x1),4.0)) + c*(y-y1); break; } return ret_val; }