From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/LDRASS.f | 467 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 467 insertions(+) create mode 100644 Dragon/src/LDRASS.f (limited to 'Dragon/src/LDRASS.f') diff --git a/Dragon/src/LDRASS.f b/Dragon/src/LDRASS.f new file mode 100644 index 0000000..609d187 --- /dev/null +++ b/Dragon/src/LDRASS.f @@ -0,0 +1,467 @@ +*DECK LDRASS + LOGICAL FUNCTION LDRASS(IPGEOM,IPRT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Reads the geometry on LCM and check compatibility for cell assemblies. +* +*Copyright: +* Copyright (C) 2002 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): R. Roy +* +*Parameters: input +* IPGEOM pointer to the geometry LCM object (L_GEOM). +* IPRT print flag (iprt=0: no print). +* +*Parameters: output +* LDRASS checking flag: =.true. if everything was 'ok' with assembly; +* =.false. if nothing was checked. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPGEOM + INTEGER IPRT +*---- +* LOCAL VARIABLES +*---- + PARAMETER ( IOUT=6, NLCM=26, NSTATE=40, NIXS=12, NIST=1 ) + CHARACTER*12 LCMNM(NLCM), GEONAM, TEXT12 + INTEGER LNLCM(NLCM),INVLCM(NIXS),INVSTA(NIST),ISTATE(NSTATE), + > MAXGRI(3),NCODE(6) + LOGICAL LL1,LL2,LCELL,LDRCEL,EMPTY,LCM +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: KEYTYP,KMERGE,KTURN,KEYCEL, + > KCHECK +*---- +* DATA STATEMENTS +*---- + DATA INVLCM/ 6, 7, 8, 9, 10, 12, 16, 17, 18, + > 20, 21, 22 / + DATA INVSTA/11 / + DATA LCMNM / 'MIX', 'MESHX', 'MESHY', 'MESHZ', 'RADIUS', + > 'SIDE', 'SPLITX', 'SPLITY', 'SPLITZ', 'SPLITR', + > 'CELL', 'COORD', 'MERGE', 'TURN','CLUSTER', + > 'NPIN', 'RPIN', 'APIN', 'BIHET', 'POURCE', + > 'PROCEL', 'IHEX', 'STATE', 'NCODE', 'ZCODE', + > 'ICODE'/ +* +* I: CELL TYPE +* J: CELL TURN (=J.EQ.1.OR.J.EQ.2) + IADR( IAXIS, I1,J1, I2,J2)= IAXIS + > + 3*(J1-1+2*(J2-1+2*(I1-1+NTYPES*(I2-1)))) +* + CALL LCMINF(IPGEOM,GEONAM,TEXT12,EMPTY,ILONG,LCM) + LDRASS= .TRUE. + IF( IPRT.GT.1 )THEN + WRITE(IOUT,'(/20H CHECKING ASSEMBLY: ,A12)') GEONAM + ENDIF + DO 10 ILCM= 1, NLCM + CALL LCMLEN(IPGEOM, LCMNM(ILCM), LNLCM(ILCM), ITPLCM ) + 10 CONTINUE +*---- +* ELIMINATES OPTIONS NOT CHECKED BY THE ROUTINE +*---- + DO 20 IIXS= 1, NIXS + IF( LNLCM(INVLCM(IIXS)).NE.0 )THEN + LDRASS= .FALSE. + GO TO 9999 + ENDIF + 20 CONTINUE + CALL LCMLEN(IPGEOM,'STATE-VECTOR',ILEN,ITPLCM) + IF( ITPLCM.NE.1 )THEN + LDRASS= .FALSE. + CALL XABORT( 'LDRASS: STATE VECTOR IS NOT AN INTEGER') + ENDIF + IF( ILEN.NE.NSTATE )THEN + LDRASS= .FALSE. + CALL XABORT( 'LDRASS: GEOMETRY HAS INVALID STATE VECTOR') + ENDIF + ISTATE(:NSTATE)=0 + CALL LCMGET(IPGEOM,'STATE-VECTOR',ISTATE) +*---- +* ELIMINATES THE INVALID OPTIONS +*---- + DO 30 IIST= 1, NIST + IF( ISTATE(INVSTA(IIST)).NE.0 )THEN + LDRASS= .FALSE. + GO TO 9999 + ENDIF + 30 CONTINUE +* + ITYPE= ISTATE(1) +*---- +* FIX DIMENSION OF CELL ASSEMBLY +*---- + IF( ITYPE.EQ.5 )THEN + NDIM= 2 + ELSEIF( ITYPE.EQ.7 )THEN + NDIM= 3 + ELSE + LDRASS= .FALSE. + GO TO 9999 + ENDIF + LREG= ISTATE(6) + NBMIX= ISTATE(7) + ISUB1= ISTATE(8) + ISUB2= ISTATE(9) +* + IF( ISUB1.NE.0 )THEN + MAXGRI(1)= MAX(1,ISTATE(3)) + MAXGRI(2)= MAX(1,ISTATE(4)) + MAXGRI(3)= MAX(1,ISTATE(5)) + NBLOCK= MAXGRI(1)*MAXGRI(2)*MAXGRI(3) + IF( NBLOCK.EQ.1 )THEN +* +* JUST ONE CELL + LDRASS=.FALSE. + GO TO 9999 + ENDIF +* +* MANY CELLS + NTYPES= ISUB2 + IF( IPRT.GT.1 )THEN + WRITE(IOUT,'(6H ,I1,13H-D ASSEMBLY: ,3I4)') + > NDIM, (MAXGRI(I),I=1,NDIM) + ENDIF + ELSE +* +* JUST ONE CELL + LDRASS=.FALSE. + GO TO 9999 + ENDIF +*---- +* RECOVERS BOUNDARY CONDITIONS. +*---- + CALL LCMLEN(IPGEOM,'NCODE', ILEN, ITPLCM) + IF( ITPLCM.NE.1 ) + > CALL XABORT('LDRASS: THE NCODE BLOCK IS '// + > 'NOT ADEQUATELY DEFINED') + IF( ILEN.NE.6 ) + > CALL XABORT('LDRASS: THE NCODE BLOCK HAS '// + > 'INCORRECT DIMENSION') + CALL LCMGET(IPGEOM,'NCODE',NCODE) + IF( NDIM.EQ.2.AND.(NCODE(5).NE.0.OR.NCODE(6).NE.0) ) + > CALL XABORT('LDRASS: 3-D NCODE VALUES FOR A 2-D ASSEMBLY...') + NOCELL= NBLOCK + NDIAG=0 + DO 40 IAL= 1, 2*NDIM + IF( NCODE(IAL).EQ.0 ) + > CALL XABORT('LDRASS: A BOUNDARY CONDITION IS MISSING.') + IF( NCODE(IAL).EQ.3 )THEN + NDIAG=NDIAG+1 + ENDIF + 40 CONTINUE +* + LL1= .FALSE. + LL2= .FALSE. + MXDIAG=0 + NOC1=0 + NOC2=0 + NOCO=0 + IF( NDIAG.GT.0 )THEN + IF( NDIAG.NE.2 ) + > CALL XABORT('LDRASS: NO MORE THAN 2 DIAGONAL CONDITIONS') + IF( MAXGRI(1).NE.MAXGRI(2)) + > CALL XABORT('LDRASS: LX=LY WITH A DIAGONAL SYMMETRY.') + MXDIAG= MAXGRI(1) + NOC1=((MXDIAG+1)*MXDIAG)/2 + NOC2=MXDIAG*MXDIAG + LL1=((NCODE(2).EQ.3).AND.(NCODE(3).EQ.3)) + LL2=((NCODE(1).EQ.3).AND.(NCODE(4).EQ.3)) + IF( LL1 )THEN + NCODE(2)= NCODE(1) + NCODE(3)= NCODE(4) + ELSEIF( LL2 )THEN + NCODE(1)= NCODE(2) + NCODE(4)= NCODE(3) + ELSE + CALL XABORT('LDRASS: THE DIAGONAL CONDITIONS '// + > 'X+: DIAG Y-: DIAG OR '// + > 'X-: DIAG Y+: DIAG ARE THE ONLY PERMITTED.') + ENDIF + NOCO=NOC1 + IF(NOC2*MAXGRI(3) .EQ. LREG) NOCO=NOC2 + NOCELL=NOCO*MAXGRI(3) + ENDIF + IF( NOCELL.GT.LREG )THEN + CALL XABORT('LDRASS: # OF CELLS OF ASSEMBLY TOO LARGE...') + ELSEIF( NOCELL.LT.LREG )THEN + CALL XABORT('LDRASS: # OF CELLS OF ASSEMBLY TOO SMALL...') + ENDIF +*---- +* CHECK IF GENERATING CELL VECTOR IS 'OK' +*---- + CALL LCMLEN(IPGEOM,'MIX', KCELL, ITPLCM) + IF( ITPLCM.NE.1 ) + > CALL XABORT('LDRASS: THE MIX BLOCK WITH CELLS IS '// + > 'NOT ADEQUATELY DEFINED') + IF( KCELL.GT.NOCELL )THEN + CALL XABORT('LDRASS: THE ASSEMBLY HAS TOO MANY CELLS...') + ELSEIF( KCELL.LT.NOCELL )THEN + CALL XABORT('LDRASS: THE ASSEMBLY HAS NOT ENOUGH CELLS...') + ENDIF + ALLOCATE(KEYTYP(NBLOCK)) + CALL LCMGET(IPGEOM,'MIX', KEYTYP) + DO 41 I=1,KCELL + KEYTYP(I)=-KEYTYP(I) + 41 CONTINUE +*---- +* CHECK IF 'MERGE' ARE CORRECTLY DEFINED +*---- + CALL LCMLEN(IPGEOM,'MERGE', KCELL, ITPLCM) + IF( KCELL.NE.0 )THEN + IF( ITPLCM.NE.1 ) + > CALL XABORT('LDRASS: THE MERGE BLOCK IS '// + > 'NOT ADEQUATELY DEFINED') + IF( KCELL.GT.NOCELL )THEN + CALL XABORT('LDRASS: THE ASSEMBLY HAS TOO MANY MERGE...') + ELSEIF( KCELL.LT.NOCELL )THEN + CALL XABORT('LDRASS: THE ASSEMBLY HAS NOT ENOUGH MERGE...') + ENDIF + ALLOCATE(KMERGE(NOCELL)) + CALL LCMGET(IPGEOM,'MERGE', KMERGE) + DO 46 IM= 1, NOCELL + IDKT= KMERGE(IM) + IF( IDKT.GT.NOCELL ) + > CALL XABORT( 'LDRASS: MERGE NUMBER > TOTAL # OF CELLS') + DO 45 JM= IM+1, NOCELL + JDKT= KMERGE(JM) + IF( JDKT.EQ.IDKT )THEN + IF( KEYTYP(IDKT).NE.KEYTYP(JDKT) )THEN + CALL XABORT( 'LDRASS: MERGE NUMBERS ARE NOT '// + > 'CONSISTENT WITH GEOMETRIC DEFINITION') + ENDIF + ENDIF + 45 CONTINUE + 46 CONTINUE + DEALLOCATE(KMERGE) + ENDIF + ALLOCATE(KTURN(NBLOCK)) + DO 47 IT= 1, NOCELL + KTURN(IT)= 1 + 47 CONTINUE + CALL LCMLEN(IPGEOM,'TURN', KCELL, ITPLCM) + IF( KCELL.NE.0 )THEN + IF( ITPLCM.NE.1 ) + > CALL XABORT('LDRASS: THE MERGE BLOCK IS '// + > 'NOT ADEQUATELY DEFINED') + IF( KCELL.GT.NOCELL )THEN + CALL XABORT('LDRASS: THE ASSEMBLY HAS TOO MANY TURN...') + ELSEIF( KCELL.LT.NOCELL )THEN + CALL XABORT('LDRASS: THE ASSEMBLY HAS NOT ENOUGH TURN...') + ENDIF + CALL LCMGET(IPGEOM,'TURN', KTURN) + ENDIF +*---- +* EFFECTIVE MOD2 FOR TURN IN CELLS +*---- + DO 48 IT= 1, NOCELL + KTURN(IT)= MOD( KTURN(IT)+1,2 )+1 + 48 CONTINUE +*---- +* CHECK IF CELL NAMES ARE 'OK' +*---- + CALL LCMLEN(IPGEOM,'CELL', KTYPES, ITPLCM) + IF( ITPLCM.NE.3 ) + > CALL XABORT('LDRASS: THE CELL NAMES ARE NOT STORED IN' + > //' CHARACTER*4') + IF( KTYPES.GT.3*NTYPES )THEN + CALL XABORT('LDRASS: THE ASSEMBLY HAS TOO MANY CELL NAMES') + ELSEIF( KTYPES.LT.3*NTYPES )THEN + CALL XABORT('LDRASS: THE ASSEMBLY HAS NOT ENOUGH CELLS NAMES') + ENDIF + ALLOCATE(KEYCEL(3*NTYPES),KCHECK(12*NTYPES*NTYPES)) + DO 50 IT= 1,12*NTYPES*NTYPES + KCHECK(IT)= 0 + 50 CONTINUE + CALL LCMGET(IPGEOM,'CELL', KEYCEL) +*---- +* FILL UP "KEYTYP" ARRAY IN THE CASE OF DIAGONAL SYMMETRY +*---- + IF( LL1 )THEN + IF(NOCO .EQ. NOC1) THEN +*---- +* LOCATE DIAGONAL ELEMENTS IN THEIR RESPECTIVE PLANES +* WHILE UNFOLDING +*---- + K=LREG + DO 72 IZ=MAXGRI(3),1,-1 + IOFF=(IZ-1)*MXDIAG*MXDIAG + DO 71 IY=MXDIAG,1,-1 + DO 60 IX=MXDIAG,IY+1,-1 + KEYTYP(IOFF+(IY-1)*MXDIAG+IX)=KEYTYP(IOFF+(IX-1)*MXDIAG+IY) + KTURN(IOFF+(IY-1)*MXDIAG+IX)=KTURN(IOFF+(IX-1)*MXDIAG+IY) + 60 CONTINUE + DO 70 IX=IY,1,-1 + KEYTYP(IOFF+(IY-1)*MXDIAG+IX)=KEYTYP(K) + KTURN(IOFF+(IY-1)*MXDIAG+IX)=KTURN(K) + K=K-1 + 70 CONTINUE + 71 CONTINUE + 72 CONTINUE + DO 77 IZ=MAXGRI(3),1,-1 + IOFF=(IZ-1)*MXDIAG*MXDIAG + DO 76 IY=1,MXDIAG + DO 75 IX=IY+1,MXDIAG + KTURN(IOFF+(IY-1)*MXDIAG+IX)= + > MOD(KTURN(IOFF+(IY-1)*MXDIAG+IX),2) + 1 + 75 CONTINUE + 76 CONTINUE + 77 CONTINUE + IF (K.NE.0) + > CALL XABORT( 'LDRASS: UNABLE TO UNFOLD '// + > 'X+: DIAG Y-: DIAG ASSEMBLY...') + ENDIF + ELSEIF( LL2 )THEN + IF(NOCO .EQ. NOC1) THEN +*---- +* LOCATE DIAGONAL ELEMENTS IN THEIR RESPECTIVE PLANES +*---- + K=LREG + DO 82 IZ=MAXGRI(3),1,-1 + IOFF=(IZ-1)*MXDIAG*MXDIAG + DO 81 IY=MXDIAG,1,-1 + DO 80 IX=MXDIAG,IY,-1 + KEYTYP(IOFF+(IY-1)*MXDIAG+IX)=KEYTYP(K) + KTURN(IOFF+(IY-1)*MXDIAG+IX)=KTURN(K) + K=K-1 + 80 CONTINUE + 81 CONTINUE + 82 CONTINUE +*---- +* UNFOLD DIAGONAL ELEMENTS FOR EACH PLANE +*---- + DO 92 IZ=1,MAXGRI(3) + IOFF=(IZ-1)*MXDIAG*MXDIAG + DO 91 IY=1,MXDIAG + DO 90 IX=1,IY-1 + KEYTYP(IOFF+(IY-1)*MXDIAG+IX)=KEYTYP(IOFF+(IX-1)*MXDIAG+IY) + KTURN(IOFF+(IY-1)*MXDIAG+IX)=KTURN(IOFF+(IX-1)*MXDIAG+IY) + 90 CONTINUE + 91 CONTINUE + 92 CONTINUE + IF (K.NE.0) + > CALL XABORT( 'LDRASS: UNABLE TO UNFOLD '// + > 'X-: DIAG Y+: DIAG ASSEMBLY...') + DO 97 IZ=MAXGRI(3),1,-1 + IOFF=(IZ-1)*MXDIAG*MXDIAG + DO 96 IY=1,MXDIAG + DO 95 IX=1,IY-1 + KTURN(IOFF+(IY-1)*MXDIAG+IX)= + > MOD(KTURN(IOFF+(IY-1)*MXDIAG+IX),2) + 1 + 95 CONTINUE + 96 CONTINUE + 97 CONTINUE + ENDIF + ENDIF + IF(IPRT .GE. 10) THEN + WRITE(IOUT,6100) + DO 600 IZ=MAXGRI(3),1,-1 + DO 601 IY=MAXGRI(2),1,-1 + IOFF=(IZ-1)*MAXGRI(2)*MAXGRI(1)+(IY-1)*MAXGRI(1) + WRITE(IOUT,6110) (KEYTYP(IOFF+IX),IX=1,MAXGRI(1)) + 601 CONTINUE + WRITE(IOUT,6111) + 600 CONTINUE + WRITE(IOUT,6101) + DO 610 IZ=MAXGRI(3),1,-1 + DO 611 IY=MAXGRI(2),1,-1 + IOFF=(IZ-1)*MAXGRI(2)*MAXGRI(1)+(IY-1)*MAXGRI(1) + WRITE(IOUT,6110) (KTURN(IOFF+IX),IX=1,MAXGRI(1)) + 611 CONTINUE + WRITE(IOUT,6111) + 610 CONTINUE + ENDIF +*---- +* TRANSLATION B.C.: CHECK BEGIN-TO-END CONNEXIONS +*---- + DO 100 IC= 1, NDIM + IF( NCODE(2*IC-1).EQ.4 )THEN + IF( NCODE(2*IC).NE.4 ) + > CALL XABORT( 'LDRASS: TRANSLATION B.C. IS NOT WELL DEFINED') + ENDIF + 100 CONTINUE +*---- +* CHECK CELL INTERFACES +*---- + IOFF1= 0 + DO 112 IZ= 1, MAXGRI(3) + DO 111 IY= 1, MAXGRI(2) + DO 110 IX= 1, MAXGRI(1) + IOFF1= IOFF1+1 + IT1= KEYTYP(IOFF1) + JT1= KTURN(IOFF1) + IF(IX.NE.MAXGRI(1).OR.(IX.EQ.MAXGRI(1).AND.NCODE(1).EQ.4))THEN + IF(IX.NE.MAXGRI(1))THEN + IOFF2= IOFF1 + 1 + ELSE + IOFF2= IOFF1 + 1 - MAXGRI(1) + ENDIF + IT2= KEYTYP(IOFF2) + JT2= KTURN(IOFF2) + IF( KCHECK(IADR(1,IT1,JT1,IT2,JT2)).EQ.0 )THEN + LCELL= LDRCEL(IPGEOM, IT1,JT1, IT2,JT2, KEYCEL, + > NTYPES, 1, NDIM, IPRT) + KCHECK(IADR(1,IT1,JT1,IT2,JT2))= 1 + ENDIF + ENDIF + IF(IY.NE.MAXGRI(2).OR.(IY.EQ.MAXGRI(2).AND.NCODE(3).EQ.4))THEN + IF( IY.NE.MAXGRI(2) )THEN + IOFF2= IOFF1 + MAXGRI(1) + ELSE + IOFF2= IOFF1 +(1-MAXGRI(2)) * MAXGRI(1) + ENDIF + IT2= KEYTYP(IOFF2) + JT2= KTURN(IOFF2) + IF( KCHECK(IADR(2,IT1,JT1,IT2,JT2)).EQ.0 )THEN + LCELL= LDRCEL(IPGEOM, IT1,JT1, IT2,JT2, KEYCEL, + > NTYPES, 2, NDIM, IPRT) + KCHECK(IADR(2,IT1,JT1,IT2,JT2))= 1 + ENDIF + ENDIF + IF(IZ.NE.MAXGRI(3).OR.(IZ.EQ.MAXGRI(3).AND.NCODE(5).EQ.4))THEN + IF( IZ.NE.MAXGRI(3) )THEN + IOFF2= IOFF1 + MAXGRI(1)*MAXGRI(2) + ELSE + IOFF2= IOFF1 +(1-MAXGRI(3)) * MAXGRI(1) * MAXGRI(2) + ENDIF + IT2= KEYTYP(IOFF2) + JT2= KTURN(IOFF2) + IF( KCHECK(IADR(3,IT1,JT1,IT2,JT2)).EQ.0 )THEN + LCELL= LDRCEL(IPGEOM, IT1,JT1, IT2,JT2, KEYCEL, + > NTYPES, 3, NDIM, IPRT) + KCHECK(IADR(3,IT1,JT1,IT2,JT2))= 1 + ENDIF + ENDIF + 110 CONTINUE + 111 CONTINUE + 112 CONTINUE + DEALLOCATE(KCHECK,KEYCEL,KTURN,KEYTYP) + 9999 IF( IPRT.GT.0 )THEN + IF( LDRASS )THEN + WRITE(IOUT,'(1H ,A12,25H ASSEMBLY IS **OK** )') GEONAM + ELSE + WRITE(IOUT,'(1H ,A12,25H WAS NOT ASSEMBLY CHECKED )') GEONAM + ENDIF + ENDIF + RETURN +*---- +* FORMATS +*---- + 6100 FORMAT(/' Assembly by cell number' ) + 6101 FORMAT(/' Assembly turns' ) + 6110 FORMAT(40(1X,I5)) + 6111 FORMAT(' ') + END -- cgit v1.2.3