C
C One-purpose code to compare the finite-difference and analytical C solutions in the twisted-crystal model. C C The program reads the output files TCFDX='tc-fdx.dat' and C TCFDY='tc-fdx.dat' of Matlab scripts 'tc-fdx.m' and 'tc-fdy.m', C calculated for relatively simple and well-defined initial C conditions by a finite difference method. C The program reads the output file TCSTRESS='stress.out' of program C 'tcgreen.for', containing the frequency-dependent initial stress C corresponding to the analytical solution (initial displacement C is constant). C The program then converts the given finite-difference solution to C the finite-difference solution with the initial conditions C corresponding to the analytical solution, and writes it into C output file TCFD='tc-fd.out'. C The program also reads analytical solution GREENTC='tc-anal.out' C calculated by program 'tcgreen.for', and writes the relative C difference between finite-difference solution TCFD='tc-fd.out' and C analytical solution GREENTC='tc-anal.out' into output file C TCCOMP='tc-fderr.out'. C C Date: 2003, May 26 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C CHARACTER*80 FILSEP,FILFDX,FILFDY,FILSTR,FILFD,FILTC,FILERR PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,LU5=5,LU6=6) C READ(*,*) FILSEP CALL RSEP1(1,FILSEP) CALL RSEP3T('TCFDX' ,FILFDX,'tc-fdx.dat') CALL RSEP3T('TCFDY' ,FILFDY,'tc-fdy.dat') CALL RSEP3T('TCSTRESS',FILSTR,'stress.out') CALL RSEP3T('TCFD' ,FILFD ,'tc-fd.out') CALL RSEP3T('GREENTC' ,FILTC ,'tc-anal.out') CALL RSEP3T('TCCOMP' ,FILERR,'tc-fderr.out') C OPEN(LU1,FILE=FILFDX) OPEN(LU2,FILE=FILFDY) OPEN(LU3,FILE=FILSTR) OPEN(LU4,FILE=FILFD ) OPEN(LU5,FILE=FILTC ) OPEN(LU6,FILE=FILERR) C 10 CONTINUE READ(LU1,*,END=90) F1,G11R,G11Y,G21R,G21Y READ(LU2,*,END=90) F2,G12R,G12Y,G22R,G22Y READ(LU3,*,END=90) F3,S11R,S11Y,S12R,S12Y,S22R,S22Y READ(LU5,*,END=90) F3,U11R,U11Y,U21R,U21Y,U31R,U31Y, * U12R,U12Y,U22R,U22Y IF(F1.NE.F3) THEN WRITE(*,*) F1,F3 STOP 1 END IF IF(F2.NE.F3) THEN WRITE(*,*) F2,F3 STOP 2 END IF V11R=G11R+G11Y*S11R+G12Y*S12R V11Y= G11Y*S11Y+G12Y*S12Y V21R=G21R+G21Y*S11R+G22Y*S12R V21Y= G21Y*S11Y+G22Y*S12Y V12R=G12R+G11Y*S12R+G12Y*S22R V12Y= G11Y*S12Y+G12Y*S22Y V22R=G22R+G21Y*S12R+G22Y*S22R V22Y= G21Y*S12Y+G22Y*S22Y WRITE(LU4,'(9(F9.6,A))') * F3,' ',V11R,' ',V11Y,' ',V21R,' ',V21Y,' 0 0 ' * ,V12R,' ',V12Y,' ',V22R,' ',V22Y,' 0 0 0 0 0 0 0 0' DIF=(V11R-U11R)**2+(V21R-U21R)**2+(V12R-U12R)**2+(V22R-U22R)**2 * +(V11R-U11R)**2+(V21R-U21R)**2+(V12R-U12R)**2+(V22R-U22R)**2 DIF=SQRT(0.5*DIF) WRITE(LU6,'(9(F9.6,1X))') F3,DIF GO TO 10 C 90 CONTINUE CLOSE(LU1) CLOSE(LU2) CLOSE(LU3) CLOSE(LU4) CLOSE(LU5) CLOSE(LU6) STOP END C C======================================================================= C INCLUDE 'sep.for' INCLUDE 'error.for' INCLUDE 'length.for' C C======================================================================= C