* D06CAF Example Program Text * Mark 20 Release. NAG Copyright 2001. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NBEDMX, NVMAX, NELTMAX, NVFIXMX, LNUME, LIWORK, + LRWORK, MSTATE, LSEED PARAMETER (NBEDMX=100,NVMAX=400,NELTMAX=2*NVMAX-1, + NVFIXMX=20,LNUME=3*NELTMAX, + LIWORK=2*NVMAX+5*NELTMAX+LNUME, + LRWORK=2*NVMAX+NELTMAX,MSTATE=200,LSEED=1) * .. Local Scalars .. DOUBLE PRECISION DELTA, HX, HY, PI2, R, RAD, SK, THETA, X1, X2, + X3, Y1, Y2, Y3 INTEGER GENID, I, IFAIL, IMAX, IND, ITRACE, J, JMAX, K, + LENSD, LSTATE, ME1, ME2, ME3, NEDGE, NELT, NQINT, + NV, NVFIX, REFTK, SUBID CHARACTER PMESH * .. Local Arrays .. DOUBLE PRECISION COOR(2,NVMAX), RWORK(LRWORK) INTEGER CONN(3,NELTMAX), EDGE(3,NBEDMX), IWORK(LIWORK), + NUMFIX(NVFIXMX), SEED(LSEED), STATE(MSTATE) * .. External Functions .. DOUBLE PRECISION X01AAF EXTERNAL X01AAF * .. External Subroutines .. EXTERNAL D06CAF, G05KFF, G05SQF * .. Intrinsic Functions .. INTRINSIC COS, DBLE, MIN, SIN * .. Executable Statements .. CONTINUE * LSTATE = MSTATE LENSD = LSEED * WRITE (NOUT,*) 'D06CAF Example Program Results' WRITE (NOUT,*) * * Skip heading in data file * READ (NIN,*) * * Read IMAX and JMAX, the number of vertices * in the x and y directions respectively. * READ (NIN,*) IMAX, JMAX * * Read distortion percentage and calculate radius * of distortion neighbourhood so that cross-over * can only occur at 100% or greater. * READ (NIN,*) DELTA * NV = IMAX*JMAX IF (NV.GT.NVMAX) THEN WRITE (NOUT,99999) NV, NVMAX GO TO 220 END IF * READ (NIN,*) PMESH * HX = 1.0D0/DBLE(IMAX-1) HY = 1.0D0/DBLE(JMAX-1) RAD = 0.005D0*DELTA*MIN(HX,HY) PI2 = 2.0D0*X01AAF(PI2) IND = 0 * * Initialise the seed for the random number generator SEED(1) = 1762541 * GENID and SUBID identify the base generator GENID = 1 SUBID = 1 * Initialise the generator to a repeatable sequence IFAIL = 1 CALL G05KFF(GENID,SUBID,SEED,LENSD,STATE,LSTATE,IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,*) WRITE (NOUT,99993) ' ** G05KFF returned with IFAIL = ', IFAIL GO TO 220 END IF * Generate two sets of uniform random variates IFAIL = 0 CALL G05SQF(NV,0.0D0,RAD,STATE,RWORK,IFAIL) CALL G05SQF(NV,0.0D0,PI2,STATE,RWORK(NV+1),IFAIL) * * Generate a simple uniform mesh and then distort it * randomly within the distortion neighbourhood of each * node. * K = 0 DO 40 J = 1, JMAX DO 20 I = 1, IMAX * K = K + 1 R = RWORK(K) THETA = RWORK(NV+K) * IF (I.EQ.1 .OR. I.EQ.IMAX .OR. J.EQ.1 .OR. J.EQ.JMAX) + R = 0.D0 * COOR(1,K) = DBLE(I-1)*HX + R*COS(THETA) COOR(2,K) = DBLE(J-1)*HY + R*SIN(THETA) * IF (I.LT.IMAX .AND. J.LT.JMAX) THEN IND = IND + 1 CONN(1,IND) = K CONN(2,IND) = K + 1 CONN(3,IND) = K + IMAX + 1 IND = IND + 1 CONN(1,IND) = K CONN(2,IND) = K + IMAX + 1 CONN(3,IND) = K + IMAX END IF 20 CONTINUE 40 CONTINUE * NELT = IND * IF (PMESH.EQ.'N') THEN WRITE (NOUT,*) 'The complete distorted mesh characteristics' WRITE (NOUT,99998) 'Number of vertices =', NV WRITE (NOUT,99998) 'Number of elements =', NELT ELSE IF (PMESH.EQ.'Y') THEN * * Output the mesh to view it using the NAG Graphics Library * WRITE (NOUT,99997) NV, NELT DO 60 I = 1, NV WRITE (NOUT,99996) COOR(1,I), COOR(2,I) 60 CONTINUE ELSE WRITE (NOUT,*) 'Problem: the printing option must be Y or N' GO TO 220 END IF * REFTK = 0 DO 80 K = 1, NELT 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 = ((X2-X1)*(Y3-Y1)-(Y2-Y1)*(X3-X1))/2.D0 IF (SK.LT.0.D0) THEN WRITE (NOUT,*) + 'Error: the surface of the element is negative' WRITE (NOUT,99998) 'element number = ', K WRITE (NOUT,99994) 'element surface = ', SK GO TO 220 END IF IF (PMESH.EQ.'Y') WRITE (NOUT,99995) CONN(1,K), CONN(2,K), + CONN(3,K), REFTK 80 CONTINUE * * Boundary edges * NEDGE = 0 DO 100 I = 1, IMAX - 1 NEDGE = NEDGE + 1 EDGE(1,NEDGE) = I EDGE(2,NEDGE) = I + 1 EDGE(3,NEDGE) = 0 100 CONTINUE * DO 120 I = 1, JMAX - 1 NEDGE = NEDGE + 1 EDGE(1,NEDGE) = I*IMAX EDGE(2,NEDGE) = (I+1)*IMAX EDGE(3,NEDGE) = 0 120 CONTINUE * DO 140 I = 1, IMAX - 1 NEDGE = NEDGE + 1 EDGE(1,NEDGE) = IMAX*JMAX - I + 1 EDGE(2,NEDGE) = IMAX*JMAX - I EDGE(3,NEDGE) = 0 140 CONTINUE * DO 160 I = 1, JMAX - 1 NEDGE = NEDGE + 1 EDGE(1,NEDGE) = (JMAX-I)*IMAX + 1 EDGE(2,NEDGE) = (JMAX-I-1)*IMAX + 1 EDGE(3,NEDGE) = 0 160 CONTINUE * NVFIX = 0 NUMFIX(1) = 0 ITRACE = 1 NQINT = 10 IFAIL = 1 * * Call the smoothing routine * CALL D06CAF(NV,NELT,NEDGE,COOR,EDGE,CONN,NVFIX,NUMFIX,ITRACE, + NQINT,IWORK,LIWORK,RWORK,LRWORK,IFAIL) * IF (IFAIL.NE.0) THEN WRITE (NOUT,*) WRITE (NOUT,99993) ' ** D06CAF returned with IFAIL = ', IFAIL GO TO 220 ELSE IF (PMESH.EQ.'N') THEN WRITE (NOUT,*) WRITE (NOUT,*) 'The complete smoothed mesh characteristics' WRITE (NOUT,99998) 'Number of vertices =', NV WRITE (NOUT,99998) 'Number of elements =', NELT ELSE IF (PMESH.EQ.'Y') THEN * * Output the mesh to view it using the NAG Graphics Library * WRITE (NOUT,99997) NV, NELT DO 180 I = 1, NV WRITE (NOUT,99996) COOR(1,I), COOR(2,I) 180 CONTINUE * REFTK = 0 DO 200 K = 1, NELT WRITE (NOUT,99995) CONN(1,K), CONN(2,K), CONN(3,K), REFTK 200 CONTINUE END IF * 220 CONTINUE * 99999 FORMAT (1X,'Problem: the number of vertices, NV is ',I6,',',/' ', + ' the maximum allowed here is, NVMAX =',I6,'.') 99998 FORMAT (1X,A,I6) 99997 FORMAT (1X,2I10) 99996 FORMAT (2(2X,E12.6)) 99995 FORMAT (1X,4I10) 99994 FORMAT (1X,A,E12.6) 99993 FORMAT (1X,A,I5) END