Index: UPNODE.f =================================================================== RCS file: /cvs/cmiss/cm/source/fe23/UPNODE.f,v retrieving revision 1.452 diff -u -r1.452 UPNODE.f --- UPNODE.f 8 Mar 2006 00:43:21 -0000 1.452 +++ UPNODE.f 22 Jul 2006 02:41:28 -0000 @@ -1172,189 +1172,142 @@ C GMH 23/4/96 Loop for subset of nodes DO nolist=1,NPLIST(0) np=NPLIST(nolist) -C Find all lines that go through the node - DO njj=1,NJ_LOC(NJL_GEOM,0,nr) - nj=NJ_LOC(NJL_GEOM,njj,nr) - DX(nj)=0.0d0 - ENDDO !nj - NL_COUNT=0 - DO nl=1,NLT - ni=NPL(1,0,nl) - NP1=NPL(2,1,nl) - NP2=NPL(3,1,nl) -C new PJH 15Sep92 - IF((NP1.EQ.np.OR.NP2.EQ.np).AND. - ' ni.EQ.NO_DERIV) THEN !update specified derivative - NL_COUNT=NL_COUNT+1 - DO njj=1,NJ_LOC(NJL_GEOM,0,nr) - nj=NJ_LOC(NJL_GEOM,njj,nr) + DO nv=1,NVJP(nj,np) + CALL WRITES(IOOP,OP_STRING,ERROR,*9999) +C Initialise temp memory + DO njj=1,NJ_LOC(NJL_GEOM,0,nr) + nj=NJ_LOC(NJL_GEOM,njj,nr) + DX(nj)=0.0d0 + ENDDO !njj +C Find all lines that go through the node and version + NL_COUNT=0 + DO nl=1,NLT ! loop over all lines + ni=NPL(1,0,nl) + NP1=NPL(2,1,nl) + NP2=NPL(3,1,nl) + ! find the version of each node that is + ! connected to the line + nv1=NVJL(1,nj,nl); + nv2=NVJL(2,nj,nl); + + IF((((NP1.EQ.np).AND.(nv1.EQ.nv)).OR. + & ((NP2.EQ.np).AND.(nv2.EQ.nv))) + & .AND.ni.EQ.NO_DERIV) THEN + NL_COUNT=NL_COUNT+1 + DO njj=1,NJ_LOC(NJL_GEOM,0,nr) + nj=NJ_LOC(NJL_GEOM,njj,nr) C!!! NOTE: this will need updating for different versions - DX(nj)=DX(nj)+XP(1,1,nj,NP2)-XP(1,1,nj,NP1) + !DX(nj)=DX(nj)+XP(1,1,nj,NP2)-XP(1,1,nj,NP1) + +C Glenn Ramsey 15/7/06 Added version selection so that only the +C version relevant to the line is considered instead of +C taking the average of all versions as it was previously. + ! update the delta using only the relevant version + DX(nj)=DX(nj)+XP(1,nv2,nj,NP2)-XP(1,nv1,nj,NP1) + + IF(DOP) THEN + WRITE(OP_STRING,'(''nl '',I3,'' np '',I3,' + & //''' nv '',I1,I3,'':'',I1,' + & //'''<--->'',I3,'':'',I1)') + & nl,np,nv,NP1,nv1,NP2,nv2 + CALL WRITES(IODI,OP_STRING,ERROR,*9999) + WRITE(OP_STRING,'(''DX('',I1,'')='',F8.3)') + & nj,DX(nj) + CALL WRITES(IODI,OP_STRING,ERROR,*9999) + ENDIF + C!!! changed GT to GE here AAY 29Aug1995 - IF(ITYP10(nr).GT.1.AND.nj.EQ.2.AND. - ' (XP(1,1,nj,NP2).GE.3.d0*PI/2.d0).AND. - ' (XP(1,1,nj,NP1).LT.PI/2.d0).AND. - ' DX(nj).GT.PI) THEN !across 0->2pi discontinuity - DX(nj)=DX(nj)-2.d0*PI - ELSE IF(ITYP10(nr).GT.1.AND.nj.EQ.2.AND. - ' (XP(1,1,nj,NP2).LT.PI/2.d0).AND. - ' (XP(1,1,nj,NP1).GT.3.d0*PI/2.d0).AND. - ' DX(nj).LT.-PI)THEN !across 2pi->0 discont'y - DX(nj)=DX(nj)+2.d0*PI - ELSE IF(ITYP10(nr).GT.2.AND.nj.EQ.3.AND. - ' (XP(1,1,nj,NP2).GT.3.d0*PI/2.d0).AND. - ' (XP(1,1,nj,NP1).LT.PI/2.d0).AND. - ' DX(nj).GT.PI)THEN !across 0->2pi discontinuity - DX(nj)=DX(nj)-2.d0*PI - ELSE IF(ITYP10(nr).GT.1.AND.nj.EQ.3.AND. - ' (XP(1,1,nj,NP2).LT.PI/2.d0).AND. - ' (XP(1,1,nj,NP1).GT.3.d0*PI/2.d0).AND. - ' DX(nj).LT.-PI)THEN !across 2pi->0 discont'y - DX(nj)=DX(nj)+2.d0*PI - ENDIF - ENDDO !njj (nj) - ENDIF !node,deriv - ENDDO !nl - SUM=0.d0 - DO njj=1,NJ_LOC(NJL_GEOM,0,nr) - nj=NJ_LOC(NJL_GEOM,njj,nr) - SUM=SUM+DX(nj)*DX(nj) - ENDDO !nj - IF(ARCLENGTH) THEN - IF(DABS(SUM).GT.RDELTA) THEN - DS=DSQRT(SUM) + IF(ITYP10(nr).GT.1.AND.nj.EQ.2.AND. + ' (XP(1,1,nj,NP2).GE.3.d0*PI/2.d0).AND. + ' (XP(1,1,nj,NP1).LT.PI/2.d0).AND. + ' DX(nj).GT.PI) THEN !across 0->2pi discontinuity + DX(nj)=DX(nj)-2.d0*PI + ELSE IF(ITYP10(nr).GT.1.AND.nj.EQ.2.AND. + ' (XP(1,1,nj,NP2).LT.PI/2.d0).AND. + ' (XP(1,1,nj,NP1).GT.3.d0*PI/2.d0).AND. + ' DX(nj).LT.-PI) THEN !across 2pi->0 discont'y + DX(nj)=DX(nj)+2.d0*PI + ELSE IF(ITYP10(nr).GT.2.AND.nj.EQ.3.AND. + ' (XP(1,1,nj,NP2).GT.3.d0*PI/2.d0).AND. + ' (XP(1,1,nj,NP1).LT.PI/2.d0).AND. + ' DX(nj).GT.PI)THEN !across 0->2pi discontinuity + DX(nj)=DX(nj)-2.d0*PI + ELSE IF(ITYP10(nr).GT.1.AND.nj.EQ.3.AND. + ' (XP(1,1,nj,NP2).LT.PI/2.d0).AND. + ' (XP(1,1,nj,NP1).GT.3.d0*PI/2.d0).AND. + ' DX(nj).LT.-PI)THEN !across 2pi->0 discont'y + DX(nj)=DX(nj)+2.d0*PI + ENDIF + ENDDO !njj (nj) + ENDIF !node,version + ENDDO !nl + SUM=0.d0 + DO njj=1,NJ_LOC(NJL_GEOM,0,nr) + nj=NJ_LOC(NJL_GEOM,njj,nr) + SUM=SUM+DX(nj)*DX(nj) + ENDDO !nj + IF(ARCLENGTH) THEN + IF(DABS(SUM).GT.RDELTA) THEN + DS=DSQRT(SUM) C!!! KAT 6Dec00: This only seems designed for arc-length scale-factors. C CPB 1/7/96 Normalise DX - DO njj=1,NJ_LOC(NJL_GEOM,0,nr) - nj=NJ_LOC(NJL_GEOM,njj,nr) - DX(nj)=DX(nj)/DS - ENDDO !nj - ENDIF - ELSE !xi derivs - IF(NL_COUNT.GT.0) THEN -C Take an average of values from the lines - DO njj=1,NJ_LOC(NJL_GEOM,0,nr) - nj=NJ_LOC(NJL_GEOM,njj,nr) - DX(nj)=DX(nj)/NL_COUNT - ENDDO !nj - ENDIF !NL_COUNT>0 - ENDIF !ARCLENGTH + DO njj=1,NJ_LOC(NJL_GEOM,0,nr) + nj=NJ_LOC(NJL_GEOM,njj,nr) + DX(nj)=DX(nj)/DS + ENDDO !nj + ENDIF + ELSE !xi derivs + IF(NL_COUNT.GT.0) THEN +C Take an average of values from the lines + DO njj=1,NJ_LOC(NJL_GEOM,0,nr) + nj=NJ_LOC(NJL_GEOM,njj,nr) + DX(nj)=DX(nj)/NL_COUNT + ENDDO !nj + ENDIF !NL_COUNT>0 + ENDIF !ARCLENGTH C KAT 6Dec00: We probably could try to obtain nk from NPL (but C according to MPN 3Feb96 these are wrong) or from C inverting NUNK but we just use the simple NKNI. This C only needs to be done once in this subroutine but if it C were done using NUNK it would be done here. - nk=NKNI(NO_DERIV) - IF(NJTOT.EQ.NJT) THEN - DO njj=1,NJ_LOC(NJL_GEOM,0,nr) - nj=NJ_LOC(NJL_GEOM,njj,nr) + nk=NKNI(NO_DERIV) + IF(NJTOT.EQ.NJT) THEN + DO njj=1,NJ_LOC(NJL_GEOM,0,nr) + nj=NJ_LOC(NJL_GEOM,njj,nr) C CS 28/12/2000 adding loop over versions - DO nv=1,NVJP(nj,np) - IF (nv.GE.2) THEN - WRITE(OP_STRING,'('' >>WARNING: Updating' - ' //' higher versions too'')') - CALL WRITES(IOER,OP_STRING,ERROR,*9999) +C GR 15/7/06 loop no longer needed because we only want to update +C the relevant version, stored in nv. +! DO nv=1,NVJP(nj,np) +! IF (nv.GE.2) THEN +! WRITE(OP_STRING,'('' >>WARNING: Updating' +! ' //' higher versions too'')') +! CALL WRITES(IOOP,OP_STRING,ERROR,*9999) +! ENDIF +! XP(nk,nv,nj,np)=DX(nj) +! ENDDO + IF(DOP) THEN + WRITE(OP_STRING,'(''updating XP np '',I3,'' nv ' + & //''',I3,'' nk '',I1,' + & //''' DX('',I1,'')='',F8.3)') + & np,nv,nk,nj,DX(nj) + CALL WRITES(IODI,OP_STRING,ERROR,*9999) ENDIF XP(nk,nv,nj,np)=DX(nj) ENDDO - ENDDO - ELSE IF(NJTOT.EQ.1)THEN - DO nv=1,NVJP(nj,np) - IF (nv.GE.2) THEN - WRITE(OP_STRING,'('' >>WARNING: Updating' - ' //' higher versions too'')') - CALL WRITES(IOER,OP_STRING,ERROR,*9999) - ENDIF - XP(nk,1,nj_update,np)=DX(nj_update) - ENDDO - ENDIF !NJTOT - -C KAT 6Dec00: We have calculated an average derivative at np. We don't -C want to store this at NP1 and NP2 for every line that -C uses this node. -C DO nl=1,NLT -C ni=NPL(1,0,nl) -C NP1=NPL(2,1,nl) -C NP2=NPL(3,1,nl) -CC GMH 8/1/97 Update cmgui link -C CALL NODE_CHANGE(NP1,.FALSE.,ERROR,*9999) -CC GMH 8/1/97 Update cmgui link -C CALL NODE_CHANGE(NP2,.FALSE.,ERROR,*9999) -C IF(((NP1.EQ.np).OR.(NP2.EQ.np)).AND. -C ' (ni.EQ.NO_DERIV)) THEN !update specified derivative -CC find interpolation type and scaling type for line -CC KAT 6Dec00: This stuff either doesn't work or isn't used. -CC AVEARCLENGTH=.FALSE. -CC ARCLENGTH=.FALSE. -C NK1=0 -C NK2=0 -C IF(NJTOT.EQ.1) THEN !one geometric variable only -C IF(NPL(1,nj_update,nl).EQ.4) THEN !cubic Herm line -CC KAT 6Dec00: This stuff either doesn't work or isn't used. -CC DO nb=1,NBFT -CC IF(IBT(1,ni,nb).EQ.2.AND. -CC ' IBT(2,ni,nb).EQ.1) THEN -CC AVEARCLENGTH=(NBI(nb).EQ.6.OR.NBI(nb).EQ.7) -CC ARCLENGTH=(NBI(nb).EQ.5) -CC ENDIF -CC ENDDO !nb -C NK1=NPL(4,1,nl) -C NK2=NPL(5,1,nl) -C ENDIF -C ELSE !loop thru all geometric variables -C DO njj=1,NJ_LOC(NJL_GEOM,0,nr) -C nj=NJ_LOC(NJL_GEOM,njj,nr) -C IF(NPL(1,nj,nl).EQ.4)THEN !cubic Hermite line -CC KAT 6Dec00: This stuff either doesn't work or isn't used. -CC DO nb=1,NBFT -CC IF(IBT(1,ni,nb).EQ.2.AND.IBT(2,ni,nb).EQ.1) -CC ' THEN -CC !cubic Hermite -CC AVEARCLENGTH= -CC ' (NBI(nb).EQ.6.OR.NBI(nb).EQ.7) -CCC!!! KAT 6Dec00: What about ARCLENGTH? -CC ENDIF -CC ENDDO !nb -C NK1=NPL(4,1,nl) -C NK2=NPL(5,1,nl) -C ENDIF -C ENDDO !njj (nj) -C ENDIF !NJTOT - -CC KAT 6Dec00: This stuff either doesn't work or isn't used. -CC IF(AVEARCLENGTH.OR.ARCLENGTH) -CC ' THEN !arc length scaling factors -CC IF(DL(3,nl).EQ.0.d0) THEN !no lengths calculated -CCC CPB 1/7/96 ALREADY NORMALISED -CCC DS=DSQRT(SUM) !is strt line arc-length btw nodes -CCC update scaleing factors -CCC!!! KAT 6Dec00: DS is not set -CC DL(1,nl)=DS !is deriv scaling factor at LH end -CC DL(2,nl)=DS !is deriv scaling factor at RH end -CC DL(3,nl)=DS !is arc-length -CC ELSE -CC DS=DL(3,nl) !is arc-length between nodes -CC ENDIF -CC ELSE !assume angle change or whatever in DL(1,nl) -CC DS=DL(1,nl) -CC ENDIF -C IF(NJTOT.EQ.NJT.AND.NK1.GT.0.AND.NK2.GT.0)THEN -C DO njj=1,NJ_LOC(NJL_GEOM,0,nr) -C nj=NJ_LOC(NJL_GEOM,njj,nr) -CC XP(NK1,1,nj,NP1)=DX(nj)/DS -CC XP(NK2,1,nj,NP2)=DX(nj)/DS -C XP(NK1,1,nj,NP1)=DX(nj) -C XP(NK2,1,nj,NP2)=DX(nj) -C ENDDO -C ELSE IF(NJTOT.EQ.1.AND.NK1.GT.0.AND.NK2.GT.0)THEN -CC XP(NK1,1,nj_update,NP1)=DX(nj_update)/DS -CC XP(NK2,1,nj_update,NP2)=DX(nj_update)/DS -C XP(NK1,1,nj_update,NP1)=DX(nj_update) -C XP(NK2,1,nj_update,NP2)=DX(nj_update) -C ENDIF !NJTOT -C ENDIF ! ni.EQ.NO_DERIV -C ENDDO !nl + ELSE IF(NJTOT.EQ.1)THEN +! DO nv=1,NVJP(nj,np) +! IF (nv.GE.2) THEN +! WRITE(OP_STRING,'('' >>WARNING: Updating' +! ' //' higher versions too'')') +! CALL WRITES(IOOP,OP_STRING,ERROR,*9999) +! ENDIF +! XP(nk,1,nj_update,np)=DX(nj_update) +! ENDDO + XP(nk,nv,nj_update,np)=DX(nj_update) + ENDIF !NJTOT +! removed commented out code see CVS + ENDDO !nv ENDDO !np DO nb=1,NBFT CALL DLSE(IBT,IDO,nb,NEELEM,NLL,NNL,NPL,