/* nag_mesh2d_smooth (d06cac) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 7, 2001.
*
* Mark 8 revised, 2004
*
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nagd06.h>
#include <nagg05.h>
#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");
fflush(stdout);
/* 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
*/
fflush(stdout);
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:
NAG_FREE(coor);
NAG_FREE(conn);
NAG_FREE(edge);
NAG_FREE(numfix);
NAG_FREE(state);
return exit_status;
}