/* nag_mesh2d_join (d06dbc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 7, 2001.
* Mark 7b revised, 2004.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd06.h>
#ifdef __cplusplus
extern "C" {
#endif
static double NAG_CALL fbnd(Integer, double, double, Nag_Comm *);
#ifdef __cplusplus
}
#endif
static int ex1(void);
static int ex2(void);
#define EDGE1(I, J) edge1[3*((J) -1)+(I) -1]
#define EDGE2(I, J) edge2[3*((J) -1)+(I) -1]
#define EDGE3(I, J) edge3[3*((J) -1)+(I) -1]
#define EDGE4(I, J) edge4[3*((J) -1)+(I) -1]
#define EDGE5(I, J) edge5[3*((J) -1)+(I) -1]
#define CONN1(I, J) conn1[3*((J) -1)+(I) -1]
#define CONN2(I, J) conn2[3*((J) -1)+(I) -1]
#define CONN3(I, J) conn3[3*((J) -1)+(I) -1]
#define CONN4(I, J) conn4[3*((J) -1)+(I) -1]
#define CONN5(I, J) conn5[3*((J) -1)+(I) -1]
#define COOR1(I, J) coor1[2*((J) -1)+(I) -1]
#define COOR2(I, J) coor2[2*((J) -1)+(I) -1]
#define COOR3(I, J) coor3[2*((J) -1)+(I) -1]
#define COOR4(I, J) coor4[2*((J) -1)+(I) -1]
#define COOR5(I, J) coor5[2*((J) -1)+(I) -1]
#define TRANS(I, J) trans[6*((J) -1)+(I) -1]
#define LINED(I, J) lined[4*((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)
{
Integer exit_status_ex1 = 0;
Integer exit_status_ex2 = 0;
printf("nag_mesh2d_join (d06dbc) Example Program Results\n");
exit_status_ex1 = ex1();
exit_status_ex2 = ex2();
return (exit_status_ex1 == 0 && exit_status_ex2 == 0 )? 0 : 1;
}
int ex1(void)
{
const Integer nvmax = 900, nedmx = 200, neltmx = 2*nvmax+5, ntrans = 1,
mode = 0;
double eps;
Integer exit_status, i, imax, itrace, itrans, jmax, jtrans, k, ktrans;
Integer nedge1, nedge2, nedge3, nelt1, nelt2, nelt3, nv1, nv2, nv3,
reftk;
Integer imaxm1, jmaxm1, ind;
char pmesh[2];
double *coor1 = 0, *coor2 = 0, *coor3 = 0, *trans = 0;
Integer *conn1 = 0, *conn2 = 0, *conn3 = 0, *edge1 = 0, *edge2 = 0;
Integer *edge3 = 0, *itype = 0, *reft1 = 0, *reft2 = 0, *reft3 = 0;
NagError fail;
INIT_FAIL(fail);
exit_status = 0;
printf("\nExample 1\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%*[^\n] ");
/* Read the mesh: coordinates and connectivity of the 1st domain */
scanf("%ld", &nv1);
scanf("%ld", &nelt1);
scanf("%*[^\n] ");
nv2 = nv1;
nelt2 = nelt1;
imax = 20;
jmax = imax;
imaxm1 = imax - 1;
jmaxm1 = jmax - 1;
nedge1 = 2*(imaxm1 + jmaxm1);
nedge2 = nedge1;
/* Allocate memory */
if (!(coor1 = NAG_ALLOC(2*nv1, double)) ||
!(coor2 = NAG_ALLOC(2*nv2, double)) ||
!(coor3 = NAG_ALLOC(2*nvmax, double)) ||
!(trans = NAG_ALLOC(6*ntrans, double)) ||
!(conn1 = NAG_ALLOC(3*nelt1, Integer)) ||
!(conn2 = NAG_ALLOC(3*nelt2, Integer)) ||
!(conn3 = NAG_ALLOC(3*neltmx, Integer)) ||
!(edge1 = NAG_ALLOC(3*nedge1, Integer)) ||
!(edge2 = NAG_ALLOC(3*nedge2, Integer)) ||
!(edge3 = NAG_ALLOC(3*nedmx, Integer)) ||
!(itype = NAG_ALLOC(ntrans, Integer)) ||
!(reft1 = NAG_ALLOC(nelt1, Integer)) ||
!(reft2 = NAG_ALLOC(nelt2, Integer)) ||
!(reft3 = NAG_ALLOC(neltmx, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= nv1; ++i)
{
scanf("%lf", &COOR1(1, i));
scanf("%lf", &COOR1(2, i));
scanf("%*[^\n] ");
}
for (k = 1; k <= nelt1; ++k)
{
scanf("%ld", &CONN1(1, k));
scanf("%ld", &CONN1(2, k));
scanf("%ld", &CONN1(3, k));
scanf("%ld", &reftk);
scanf("%*[^\n] ");
reft1[k-1] = 1;
reft2[k-1] = 2;
}
scanf(" ' %1s '", pmesh);
scanf("%*[^\n] ");
/* the edges of the boundary */
ind = 0;
for (i = 1; i <= imaxm1; ++i)
{
++ind;
EDGE1(1, ind) = i;
EDGE1(2, ind) = i + 1;
EDGE1(3, ind) = 1;
}
for (i = 1; i <= jmaxm1; ++i)
{
++ind;
EDGE1(1, ind) = i*imax;
EDGE1(2, ind) = (i+1)*imax;
EDGE1(3, ind) = 1;
}
for (i = 1; i <= imaxm1; ++i)
{
++ind;
EDGE1(1, ind) = imax*jmax - i + 1;
EDGE1(2, ind) = imax*jmax - i;
EDGE1(3, ind) = 1;
}
for (i = 1; i <= jmaxm1; ++i)
{
++ind;
EDGE1(1, ind) = (jmax - i)*imax + 1;
EDGE1(2, ind) = (jmax - i - 1)*imax + 1;
EDGE1(3, ind) = 1;
}
for (ktrans = 0; ktrans < 2; ++ktrans)
{
/* Translation of the 1st domain to obtain the 2nd domain */
/* KTRANS = 0 leading to a domain overlapping */
/* KTRANS = 1 leading to a domain partition */
if (ktrans == 0)
{
itrans = imax - 5;
jtrans = jmax - 3;
}
else
{
itrans = imax - 1;
jtrans = 0;
}
itype[0] = 1;
TRANS(1, 1) = (double) itrans/(imax - 1.0);
TRANS(2, 1) = (double) jtrans/(jmax - 1.0);
itrace = 0;
/* nag_mesh2d_trans (d06dac).
* Generates a mesh resulting from an affine transformation
* of a given mesh
*/
nag_mesh2d_trans(mode, nv2, nedge2, nelt2, ntrans, itype, trans, coor1,
edge1, conn1, coor2, edge2, conn2, itrace, 0,
&fail);
if (fail.code == NE_NOERROR)
{
for (i = 1; i <= nedge2; ++i) EDGE2(3, i) = 2;
/* Call to the restitching driver */
itrace = 0;
eps = 0.01;
/* nag_mesh2d_join (d06dbc).
* Joins together two given adjacent (possibly overlapping)
* meshes
*/
nag_mesh2d_join(eps, nv1, nelt1, nedge1, coor1, edge1, conn1, reft1,
nv2, nelt2, nedge2, coor2, edge2, conn2, reft2, &nv3,
&nelt3, &nedge3, coor3, edge3, conn3, reft3, itrace,
0, &fail);
if (fail.code == NE_NOERROR)
{
if (pmesh[0] == 'N')
{
if (ktrans == 0)
{
printf("The restitched mesh characteristics\n");
printf("in the overlapping case\n");
}
else
{
printf("in the partition case\n");
}
printf(" nv =%6ld\n", nv3);
printf(" nelt =%6ld\n", nelt3);
printf(" nedge =%6ld\n", nedge3);
}
else if (pmesh[0] == 'Y')
{
/* Output the mesh to view it using the
NAG Graphics Library */
printf(" %10ld%10ld%10ld\n", nv3,
nelt3, nedge3);
for (i = 1; i <= nv3; ++i)
printf(" %15.6e %15.6e\n",
COOR3(1, i), COOR3(2, i));
for (k = 1; k <= nelt3; ++k)
printf(" %10ld%10ld%10ld"
"%10ld\n", CONN3(1, k), CONN3(2, k),
CONN3(3, k), reft3[k - 1]);
for (k = 1; k <= nedge3; ++k)
printf(" %10ld%10ld%10ld\n",
EDGE3(1, k), EDGE3(2, k), EDGE3(3, k));
}
else
{
printf("Problem with the printing option Y or N\n");
exit_status = -1;
goto END;
}
}
else
{
printf("Error from nag_mesh2d_join (d06dbc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
}
else
{
printf("Error from nag_mesh2d_trans (d06dac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
}
END:
NAG_FREE(coor1);
NAG_FREE(coor2);
NAG_FREE(coor3);
NAG_FREE(trans);
NAG_FREE(conn1);
NAG_FREE(conn2);
NAG_FREE(conn3);
NAG_FREE(edge1);
NAG_FREE(edge2);
NAG_FREE(edge3);
NAG_FREE(itype);
NAG_FREE(reft1);
NAG_FREE(reft2);
NAG_FREE(reft3);
return exit_status;
}
int ex2(void)
{
const Integer nvmax = 900, nedmx = 200, neltmx = 2*nvmax+5,
ntrans = 1, nus = 0, nvint = 0, nvfix = 0;
static double ruser[1] = {-1.0};
double eps;
Integer exit_status, i, itrace, j, k, l, ncomp, nedge1, nedge2, nedge3;
Integer nedge4, nedge5, nelt1, nelt2, nelt3, nelt4, nelt5, nlines;
Integer npropa, nqint, nv1, nv2, nv3, nv4, nv5, nvb1, nvb2, mode;
char pmesh[2];
double *coor1 = 0, *coor2 = 0, *coor3 = 0, *coor4 = 0, *coor5 = 0;
double *coorch = 0, *coorus = 0, *rate = 0, *trans = 0, *weight = 0;
Integer *conn1 = 0, *conn2 = 0, *conn3 = 0, *conn4 = 0, *conn5 = 0;
Integer *edge1 = 0, *edge2 = 0, *edge3 = 0, *edge4 = 0, *edge5 = 0;
Integer *itype = 0, *lcomp = 0, *lined = 0, *nlcomp = 0, *numfix = 0;
Integer *reft1 = 0, *reft2 = 0, *reft3 = 0, *reft4 = 0, *reft5 = 0;
NagError fail;
Nag_Comm comm;
INIT_FAIL(fail);
/* For communication with user-supplied functions: */
comm.user = ruser;
exit_status = 0;
printf("\nExample 2\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Build the mesh of the 1st domain */
/* Initialise boundary mesh inputs: */
/* the number of line and of the characteristic points of */
/* the boundary mesh */
scanf("%ld", &nlines);
scanf("%*[^\n] ");
/* Allocate memory */
if (!(coor1 = NAG_ALLOC(2*nvmax, double)) ||
!(coor2 = NAG_ALLOC(2*nvmax, double)) ||
!(coor3 = NAG_ALLOC(2*nvmax, double)) ||
!(coor4 = NAG_ALLOC(2*nvmax, double)) ||
!(coor5 = NAG_ALLOC(2*nvmax, double)) ||
!(coorch = NAG_ALLOC(2*nlines, double)) ||
!(coorus = NAG_ALLOC(1, double)) ||
!(rate = NAG_ALLOC(nlines, double)) ||
!(trans = NAG_ALLOC(6*ntrans, double)) ||
!(weight = NAG_ALLOC(1, double)) ||
!(conn1 = NAG_ALLOC(3*neltmx, Integer)) ||
!(conn2 = NAG_ALLOC(3*neltmx, Integer)) ||
!(conn3 = NAG_ALLOC(3*neltmx, Integer)) ||
!(conn4 = NAG_ALLOC(3*neltmx, Integer)) ||
!(conn5 = NAG_ALLOC(3*neltmx, Integer)) ||
!(edge1 = NAG_ALLOC(3*nedmx, Integer)) ||
!(edge2 = NAG_ALLOC(3*nedmx, Integer)) ||
!(edge3 = NAG_ALLOC(3*nedmx, Integer)) ||
!(edge4 = NAG_ALLOC(3*nedmx, Integer)) ||
!(edge5 = NAG_ALLOC(3*nedmx, Integer)) ||
!(itype = NAG_ALLOC(ntrans, Integer)) ||
!(lcomp = NAG_ALLOC(nlines, Integer)) ||
!(lined = NAG_ALLOC(4*nlines, Integer)) ||
!(numfix = NAG_ALLOC(1, Integer)) ||
!(reft1 = NAG_ALLOC(2*nvmax+5, Integer)) ||
!(reft2 = NAG_ALLOC(2*nvmax+5, Integer)) ||
!(reft3 = NAG_ALLOC(2*nvmax+5, Integer)) ||
!(reft4 = NAG_ALLOC(2*nvmax+5, Integer)) ||
!(reft5 = NAG_ALLOC(2*nvmax+5, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Characteristic points of the boundary geometry */
for (j = 1; j <= nlines; ++j)
{
scanf("%lf", &COORCH(1, j));
}
scanf("%*[^\n] ");
for (j = 1; j <= nlines; ++j)
{
scanf("%lf", &COORCH(2, j));
}
scanf("%*[^\n] ");
/* The lines of the boundary mesh */
for (j = 1; j <= nlines; ++j)
{
for (i = 1; i <= 4; ++i) scanf("%ld", &LINED(i, j));
scanf("%lf", &rate[j - 1]);
}
scanf("%*[^\n] ");
/* The number of connected components */
/* on the boundary and their data */
scanf("%ld", &ncomp);
scanf("%*[^\n] ");
/* Allocate memory */
if (!(nlcomp = NAG_ALLOC(ncomp, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
j = 1;
for (i = 1; i <= ncomp; ++i)
{
scanf("%ld", &nlcomp[i - 1]);
scanf("%*[^\n] ");
l = j + abs(nlcomp[i - 1]) - 1;
for (k = j; k <= l; ++k) scanf("%ld", &lcomp[k - 1]);
scanf("%*[^\n] ");
j += abs(nlcomp[i - 1]);
}
itrace = 0;
/* Call to the 2D boundary mesh generator */
/* nag_mesh2d_bound (d06bac).
* Generates a boundary mesh
*/
nag_mesh2d_bound(nlines, coorch, lined, fbnd, coorus, nus, rate, ncomp,
nlcomp, lcomp, nvmax, nedmx, &nvb1, coor1, &nedge1,
edge1, itrace, 0, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_mesh2d_bound (d06bac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Generate mesh using Delaunay-Voronoi method */
/* Initialise mesh control parameters */
itrace = 0;
npropa = 1;
/* nag_mesh2d_delaunay (d06abc).
* Generates a two-dimensional mesh using a Delaunay-Voronoi
* process
*/
nag_mesh2d_delaunay(nvb1, nvint, nvmax, nedge1, edge1, &nv1, &nelt1, coor1,
conn1, weight, npropa, itrace, 0, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_mesh2d_delaunay (d06abc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
for (k = 1; k <= nelt1; ++k) reft1[k - 1] = 1;
/* Call the smoothing routine */
nqint = 10;
/* nag_mesh2d_smooth (d06cac).
* Uses a barycentering technique to smooth a given mesh
*/
nag_mesh2d_smooth(nv1, nelt1, nedge1, coor1, edge1, conn1, nvfix, numfix,
itrace, 0, nqint, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_mesh2d_smooth (d06cac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Build the mesh of the 2nd domain */
scanf("%ld", &nlines);
scanf("%*[^\n] ");
/* Characteristic points of the boundary geometry */
for (j = 1; j <= nlines; ++j) scanf("%lf", &COORCH(1, j));
scanf("%*[^\n] ");
for (j = 1; j <= nlines; ++j) scanf("%lf", &COORCH(2, j));
scanf("%*[^\n] ");
/* The lines of the boundary mesh */
for (j = 1; j <= nlines; ++j)
{
for (i = 1; i <= 4; ++i) scanf("%ld", &LINED(i, j));
scanf("%lf", &rate[j - 1]);
}
scanf("%*[^\n] ");
/* The number of connected components */
/* to the boundary and their data */
scanf("%ld", &ncomp);
scanf("%*[^\n] ");
j = 1;
for (i = 1; i <= ncomp; ++i)
{
scanf("%ld", &nlcomp[i - 1]);
scanf("%*[^\n] ");
for (k = j; k <= j+abs(nlcomp[i-1])-1; ++k)
scanf("%ld", &lcomp[k - 1]);
scanf("%*[^\n] ");
j += abs(nlcomp[i-1]);
}
scanf(" ' %1s '", pmesh);
scanf("%*[^\n] ");
itrace = 0;
/* Call to the 2D boundary mesh generator */
/* nag_mesh2d_bound (d06bac), see above. */
nag_mesh2d_bound(nlines, coorch, lined, fbnd, coorus, nus, rate, ncomp,
nlcomp, lcomp, nvmax, nedmx, &nvb2, coor2, &nedge2, edge2,
itrace, 0, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_mesh2d_bound (d06bac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Generate mesh using the advancing front method */
itrace = 0;
/* nag_mesh2d_front (d06acc).
* Generates a two-dimensional mesh using an Advancing-front
* method
*/
nag_mesh2d_front(nvb2, nvint, nvmax, nedge2, edge2, &nv2, &nelt2, coor2,
conn2, weight, itrace, 0, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_mesh2d_front (d06acc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
for (k = 1; k <= nelt2; ++k) reft2[k - 1] = 2;
/* Rotation of the 2nd domain mesh */
/* to produce the 3rd mesh domain */
itype[0] = 3;
TRANS(1, 1) = 6.0;
TRANS(2, 1) = -1.0;
TRANS(3, 1) = -90.0;
itrace = 0;
mode = 0;
/* nag_mesh2d_trans (d06dac), see above. */
nag_mesh2d_trans(mode, nv2, nedge2, nelt2, ntrans, itype, trans, coor2,
edge2, conn2, coor3, edge3, conn3, itrace, 0, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_mesh2d_trans (d06dac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
nv3 = nv2;
nelt3 = nelt2;
nedge3 = nedge2;
for (k = 1; k <= nelt3; ++k) reft3[k - 1] = 3;
/* Restitching meshes 1 and 2 to form mesh 4 */
eps = 0.001;
itrace = 0;
/* nag_mesh2d_join (d06dbc), see above. */
nag_mesh2d_join(eps, nv1, nelt1, nedge1, coor1, edge1, conn1, reft1, nv2,
nelt2, nedge2, coor2, edge2, conn2, reft2, &nv4, &nelt4,
&nedge4, coor4, edge4, conn4, reft4, itrace, 0, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_mesh2d_join (d06dbc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Restitching meshes 3 and 4 to form mesh 5 */
itrace = 0;
/* nag_mesh2d_join (d06dbc), see above. */
nag_mesh2d_join(eps, nv4, nelt4, nedge4, coor4, edge4, conn4, reft4, nv3,
nelt3, nedge3, coor3, edge3, conn3, reft3, &nv5, &nelt5,
&nedge5, coor5, edge5, conn5, reft5, itrace, 0, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_mesh2d_join (d06dbc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
if (pmesh[0] == 'N')
{
printf("The restitched mesh characteristics\n");
printf(" nv =%6ld\n", nv5);
printf(" nelt =%6ld\n", nelt5);
printf(" nedge =%6ld\n", nedge5);
}
else if (pmesh[0] == 'Y')
{
/* Output the mesh to view it using the NAG Graphics Library */
printf(" %10ld%10ld%10ld\n", nv5, nelt5, nedge5);
for (i = 1; i <= nv5; ++i)
printf(" %15.6e %15.6e\n",
COOR5(1, i), COOR5(2, i));
for (k = 1; k <= nelt5; ++k)
printf("%10ld%10ld%10ld%10ld\n",
CONN5(1, k), CONN5(2, k), CONN5(3, k), reft5[k - 1]);
for (k = 1; k <= nedge5; ++k)
printf(" %10ld%10ld%10ld\n",
EDGE5(1, k), EDGE5(2, k), EDGE5(3, k));
}
else
{
printf("Problem with the printing option Y or N\n");
}
END:
NAG_FREE(coor1);
NAG_FREE(coor2);
NAG_FREE(coor3);
NAG_FREE(coor4);
NAG_FREE(coor5);
NAG_FREE(coorch);
NAG_FREE(coorus);
NAG_FREE(rate);
NAG_FREE(trans);
NAG_FREE(weight);
NAG_FREE(conn1);
NAG_FREE(conn2);
NAG_FREE(conn3);
NAG_FREE(conn4);
NAG_FREE(conn5);
NAG_FREE(edge1);
NAG_FREE(edge2);
NAG_FREE(edge3);
NAG_FREE(edge4);
NAG_FREE(edge5);
NAG_FREE(itype);
NAG_FREE(lcomp);
NAG_FREE(lined);
NAG_FREE(nlcomp);
NAG_FREE(numfix);
NAG_FREE(reft1);
NAG_FREE(reft2);
NAG_FREE(reft3);
NAG_FREE(reft4);
NAG_FREE(reft5);
return exit_status;
}
double NAG_CALL fbnd(Integer i, double x, double y, Nag_Comm *pcomm)
{
double radius2, x0, y0, ret_val;
if (pcomm->user[0] == -1.0)
{
printf("(User-supplied callback fbnd, first invocation.)\n");
pcomm->user[0] = 0.0;
}
ret_val = 0.0;
switch (i)
{
case 1:
/* inner circle */
x0 = 0.0;
y0 = 0.0;
radius2 = 1.0;
ret_val = (x-x0)*(x-x0) + (y-y0)*(y-y0) - radius2;
break;
case 2:
/* outer circle */
x0 = 0.0;
y0 = 0.0;
radius2 = 5.0;
ret_val = (x-x0)*(x-x0) + (y-y0)*(y-y0) - radius2;
break;
}
return ret_val;
}