/* nag_mesh_dim2_smooth_bary (d06cac) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.3, 2022.
*
*
*/
#include <math.h>
#include <nag.h>
#include <stdio.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;
/* Initialize the error structure */
INIT_FAIL(fail);
/* Get the length of the state array */
lstate = -1;
nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf(" nag_mesh_dim2_smooth_bary (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("%" NAG_IFMT "", &imax);
scanf("%" NAG_IFMT "", &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);
/* Initialize the generator to a repeatable sequence */
nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_init_repeat (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_dist_uniform(1, 0.0, rad, state, one_draw, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_dist_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_dist_uniform(1, 0.0, dpi, state, one_draw, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_dist_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 =%6" NAG_IFMT "\n", nv);
printf(" nelt =%6" NAG_IFMT "\n\n", nelt);
} else if (pmesh[0] == 'Y') {
/* Output the mesh to view it using the NAG Graphics Library */
printf(" %10" NAG_IFMT "%10" NAG_IFMT "\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 = %6" NAG_IFMT "\n", k);
printf(" sk = %15.6e\n", sk);
exit_status = -1;
goto END;
}
if (pmesh[0] == 'Y')
printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "\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_mesh_dim2_smooth_bary (d06cac).
* Uses a barycentering technique to smooth a given mesh
*/
fflush(stdout);
nag_mesh_dim2_smooth_bary(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 =%6" NAG_IFMT "\n", nv);
printf(" nelt =%6" NAG_IFMT "\n", nelt);
} else if (pmesh[0] == 'Y') {
/* Output the mesh to view it using the NAG Graphics Library */
printf(" %10" NAG_IFMT "%10" NAG_IFMT "\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(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT
"\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_mesh_dim2_smooth_bary (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;
}