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/PSPTHR.f | 503 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 503 insertions(+) create mode 100644 Dragon/src/PSPTHR.f (limited to 'Dragon/src/PSPTHR.f') diff --git a/Dragon/src/PSPTHR.f b/Dragon/src/PSPTHR.f new file mode 100644 index 0000000..9abdbd7 --- /dev/null +++ b/Dragon/src/PSPTHR.f @@ -0,0 +1,503 @@ +*DECK PSPTHR + SUBROUTINE PSPTHR(IPTRK ,ISPSP ,IPRINT,ICEL ,NDIM ,NFREG , + > MAXMSH,MXGREG,MAXPIN,KPSP ,COLREG,FACT , + > CELLPO,IDREG ,ITPIN ,DCMESH,DRAPIN, + > NBPTS ,POSTRI) +* +*----------------------------------------------------------------------- +* +*Purpose: +* To draw an hexagonal cell with triangular mesh according +* to its explicit position in the assembly. +* +*Copyright: +* Copyright (C) 2012 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): +* G. Marleau +* +*Parameters: input +* IPTRK pointer to the TRACKING data structure in +* update or creation mode. +* ISPSP POSTSCRIPT file index. +* IPRINT print level. +* ICEL cell number. +* NDIM problem dimensions. +* NFREG number of regions. +* MAXMSH maximum number of elements in MESH array. +* MXGREG maximum number of region for any geometry. +* MAXPIN maximum number of pins in a cell. +* KPSP PSP plot options. +* FACT scale factor for drawing. +* CELLPO global cell position in space. +* COLREG region color. +* +*Parameters: temporary storage +* IDREG local region identifier. +* ITPIN pin type identifier. +* DCMESH meshing vector for geometries. +* DRAPIN pin position identifier. +* NBPTS number of corners for regions in first sector. +* POSTRI positions of corners for regions in first sector. +* +*---------- +* + USE GANLIB + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + TYPE(C_PTR) IPTRK,ISPSP + INTEGER IPRINT,ICEL,NDIM,NFREG,MAXMSH,MXGREG,MAXPIN + INTEGER KPSP(7) + REAL COLREG(4,NFREG) + DOUBLE PRECISION FACT,CELLPO(2,2) + INTEGER IDREG(MXGREG),ITPIN(3,MAXPIN) + DOUBLE PRECISION DCMESH(-1:MAXMSH,4),DRAPIN(-1:4,MAXPIN) + INTEGER NBPTS(MXGREG) + DOUBLE PRECISION POSTRI(2,4,MXGREG) +*---- +* Local parameters +*---- + INTEGER IOUT + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,NAMSBR='PSPTHR') + INTEGER NSTATE + PARAMETER (NSTATE=40) + REAL WLINE + PARAMETER (WLINE=0.002) + DOUBLE PRECISION DZERO,DONE,DTWO,DHALF,DSQ3O2 + PARAMETER (DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0, + > DHALF=0.5D0,DSQ3O2=0.86602540378444D0) + INTEGER MAXDIM + PARAMETER (MAXDIM=4) +*---- +* Functions +*---- + DOUBLE PRECISION XDRCST,PI +*---- +* Local variables +*---- + INTEGER ILEV,ICONT,ICOL,KFS,KFR,KSS,KSR,NPTS,NSEG,NINT + INTEGER IEDIMC(NSTATE),IEDIMP(NSTATE) + INTEGER NX,MESHC(4),NREGC,NTPIN,NRTP + INTEGER IDIR,ILXY,IREG,IR,IY,IX,ILOC,IORDER(16) + INTEGER JR,NAR,ITREG,IAREG,JREG + REAL WLFAC,XYPOS(2,4),CENTER(2),RCIRC,RADANG(2,16) + CHARACTER NAMCEL*9,NAMREC*12 + REAL COLWHI(4),CENTEP(2),CENTED(2),CENTEB(2) + INTEGER IWCOL,IPIN,ISEG,NBRP,MESHP(4) + DOUBLE PRECISION ROTAX,COSDIR(3) + INTEGER ISECT,ITRI,KREG,NCPT + DOUBLE PRECISION SIDET,SIDEL,SIDEH,POSCXD,POSCYD,POSCX,POSCY, + > SIDEU,SIDEUH,DTMPX,DTMPY +*---- +* Data +*---- + CHARACTER CDIR(MAXDIM)*1 + SAVE CDIR + CHARACTER CLEV(2)*1 + SAVE CLEV + DATA CDIR /'X','Y','Z','R'/ + DATA CLEV /'C','P'/ +*---- +* Processing starts: +* print routine openning output header if required +* and initialize various parameters. +*---- + IF(IPRINT .GE. 100) THEN + WRITE(IOUT,6000) NAMSBR + ENDIF + IF(NDIM .GE. 3) CALL XABORT(NAMSBR// + >': PSP cannot treat 3D geometries') + PI=XDRCST('Pi',' ') + ILEV=1 + IWCOL=1 + COLWHI(:4)=1.0 +*---- +* PSP print control +*---- + WLFAC=1.0 + ICONT=KPSP(1) + ICOL=KPSP(2) + IF(KPSP(3) .EQ. 1) WLFAC=2.5 + KFS=KPSP(4) + KFR=KPSP(5) + KSS=KPSP(6) + KSR=KPSP(7) + NPTS=3 + NINT=16 +*---- +* Read cell information +*---- + WRITE(NAMCEL,'(A1,I8.8)') CLEV(ILEV),ICEL + NAMREC=NAMCEL//'DIM' + IEDIMC(:NSTATE)=0 + CALL LCMGET(IPTRK,NAMREC,IEDIMC) + NX=IEDIMC(3) + NAR=IEDIMC(2) + MESHC(1)=NX + MESHC(2)=NX + MESHC(3)=IEDIMC(5) + MESHC(4)=NAR + NREGC=IEDIMC(8) + NTPIN=IEDIMC(18) + NAMREC=NAMCEL//'RID' + CALL LCMGET(IPTRK,NAMREC,IDREG) + DO IDIR=1,4 + NAMREC=NAMCEL//'SM'//CDIR(IDIR) + IF(MESHC(IDIR) .GT. 0) THEN + CALL LCMGET(IPTRK,NAMREC,DCMESH(-1,IDIR)) + ENDIF + ENDDO + IF(NTPIN .GT .0) THEN + NAMREC=NAMCEL//'PIN' + CALL LCMGET(IPTRK,NAMREC,DRAPIN) + NAMREC=NAMCEL//'PNT' + CALL LCMGET(IPTRK,NAMREC,ITPIN) + ENDIF + CENTER(1)=REAL(FACT*(DCMESH(-1,1)+CELLPO(1,1))) + CENTER(2)=REAL(FACT*(DCMESH(-1,2)+CELLPO(2,1))) + CENTEB(1)=-CENTER(1) + CENTEB(2)=-CENTER(2) + NRTP=NX**2 + SIDEU=DCMESH(1,1)-DCMESH(0,1) + SIDET=SIDEU + IF(NX .GT. 1) THEN + SIDET=DCMESH(2,1)-DCMESH(1,1) + ENDIF + SIDEL=SIDET/DSQ3O2 + SIDEH=DHALF*SIDEL +*---- +* Find triangle corners +* 1- First sectors +* All crown +*---- + IREG=0 + POSCXD=DZERO + POSCYD=SIDEH + DO IX=1,NX + POSCX=POSCXD + POSCY=POSCYD +*---- +* right triangles +*---- + DO IR=1,IX-1 + IREG=IREG+1 + POSTRI(1,1,IREG)=POSCX+SIDET + POSTRI(2,1,IREG)=POSCY + POSTRI(1,2,IREG)=POSCXD + POSTRI(2,2,IREG)=POSCY+SIDEH + POSTRI(1,3,IREG)=POSCXD + POSTRI(2,3,IREG)=POSCY-SIDEH + NBPTS(IREG)=3 + POSCY=POSCY+SIDEL + ENDDO +*---- +* Loop over Left triangles on the line +*---- + POSCYD=POSCYD-SIDEH + POSCX=POSCXD + POSCY=POSCYD + DO IR=1,IX + IREG=IREG+1 + POSTRI(1,1,IREG)=POSCX + POSTRI(2,1,IREG)=POSCY + POSTRI(1,2,IREG)=POSCX+SIDET + POSTRI(2,2,IREG)=POSCY-SIDEH + POSTRI(1,3,IREG)=POSCX+SIDET + POSTRI(2,3,IREG)=POSCY+SIDEH + POSCY=POSCY+SIDEL + NBPTS(IREG)=3 + ENDDO + POSCXD=POSCXD+SIDET + ENDDO +*---- +* Replace right triangles in last crown by +* Isosceles trapezoid when required +*---- + IF(SIDEU .LT. SIDET) THEN + SIDEU=SIDET-SIDEU + SIDEUH=DHALF*SIDEU/DSQ3O2 + IX=NX + IREG=IREG-2*IX+1 +*---- +* Cut left triangle in last crown when required +*---- + DO IR=1,IX-1 + IREG=IREG+1 + POSTRI(1,4,IREG)=POSTRI(1,1,IREG)-SIDEU + POSTRI(2,4,IREG)=POSTRI(2,1,IREG)-SIDEUH + POSTRI(1,1,IREG)=POSTRI(1,1,IREG)-SIDEU + POSTRI(2,1,IREG)=POSTRI(2,1,IREG)+SIDEUH + NBPTS(IREG)=4 + ENDDO +*---- +* Loop over Left triangles on the line +*---- + DO IR=1,IX + IREG=IREG+1 + POSTRI(1,2,IREG)=POSTRI(1,2,IREG)-SIDEU + POSTRI(2,2,IREG)=POSTRI(2,2,IREG)+SIDEUH + POSTRI(1,3,IREG)=POSTRI(1,3,IREG)-SIDEU + POSTRI(2,3,IREG)=POSTRI(2,3,IREG)-SIDEUH + ENDDO + ENDIF +*---- +* General location of regions in first sector with respect to +* cell center are now known +* Draw all regions in first sector +*---- + IREG=0 + ISECT=1 + DO IR=1,NRTP + IREG=IREG+1 + ITREG=IREG*(NAR+1) + NCPT=NBPTS(IR) + KREG=ABS(IDREG(ITREG)) + DO ITRI=1,NCPT + XYPOS(1,ITRI)=REAL(FACT*(POSTRI(1,ITRI,IR)+CELLPO(1,1))) + XYPOS(2,ITRI)=REAL(FACT*(POSTRI(2,ITRI,IR)+CELLPO(2,1))) + ENDDO + IF(IPRINT .GE. 100) THEN + WRITE(IOUT,6010) ISECT,IR,IREG,ITREG,KREG,NCPT + WRITE(IOUT,6011) (XYPOS(1,ITRI),XYPOS(2,ITRI),ITRI=1,NCPT) + ENDIF + IF(KREG .NE. 0) THEN +*---- +* Color and trace result +*---- + CALL PSDREG(ISPSP,NCPT,XYPOS) + IF(ICOL. GT. 0) THEN + CALL PSFILL(ISPSP,ICOL,COLREG(1,KREG),KFS,KFR) + ENDIF + IF(ICONT.EQ.1) THEN + CALL PSSTRK(ISPSP,WLFAC*WLINE,KSS,KSR) + ENDIF + ENDIF + DO JR=NAR,1,-1 + IAREG=(NAR+1)*(IREG-1)+JR + JREG=ABS(IDREG(IAREG)) + IF(IPRINT .GE. 100) THEN + WRITE(IOUT,6020) JR,IAREG,JREG,FACT*DCMESH(JR,4) + ENDIF +*---- +* Annular region +*---- + IF(JREG .NE. 0) THEN + RCIRC=REAL(FACT*DCMESH(JR,4)) +*---- +* Move cursor to center of annulus +*---- + CALL PSMOVE(ISPSP,CENTER,-3) +*---- +* Color and trace result +*---- + CALL PSPRAI(NINT,NPTS,XYPOS,CENTER,RCIRC, + > NSEG,IORDER,RADANG) + CALL PSDRAI(ISPSP,NSEG,IORDER,CENTER,RADANG) + IF(ICOL. GT. 0) THEN + CALL PSFILL(ISPSP,ICOL,COLREG(1,JREG),KFS,KFR) + ENDIF + IF(ICONT.EQ.1) THEN + CALL PSSTRK(ISPSP,WLFAC*WLINE,KSS,KSR) + ENDIF +*---- +* Return cursor to original position +*---- + CALL PSMOVE(ISPSP,CENTEB,-3) + ENDIF + ENDDO + ENDDO +*---- +* 2) Rotate regions for remaining sectors +* and draw +*---- + DO ISECT=2,6 + DO IR=1,NRTP + IREG=IREG+1 + ITREG=IREG*(NAR+1) + NCPT=NBPTS(IR) + KREG=ABS(IDREG(ITREG)) + DO ITRI=1,NCPT +*---- +* rotation by pi/3 from previous sector +*---- + DTMPX=DHALF*POSTRI(1,ITRI,IR) + > -DSQ3O2*POSTRI(2,ITRI,IR) + DTMPY=DSQ3O2*POSTRI(1,ITRI,IR) + > +DHALF*POSTRI(2,ITRI,IR) + POSTRI(1,ITRI,IR)=DTMPX + POSTRI(2,ITRI,IR)=DTMPY + XYPOS(1,ITRI)=REAL(FACT*(POSTRI(1,ITRI,IR)+CELLPO(1,1))) + XYPOS(2,ITRI)=REAL(FACT*(POSTRI(2,ITRI,IR)+CELLPO(2,1))) + ENDDO + IF(IPRINT .GE. 100) THEN + WRITE(IOUT,6010) ISECT,IR,IREG,ITREG,KREG,NCPT + WRITE(IOUT,6011) (XYPOS(1,ITRI),XYPOS(2,ITRI),ITRI=1,NCPT) + ENDIF + IF(KREG .NE. 0) THEN +*---- +* Color and trace result +*---- + CALL PSDREG(ISPSP,NCPT,XYPOS) + IF(ICOL. GT. 0) THEN + CALL PSFILL(ISPSP,ICOL,COLREG(1,KREG),KFS,KFR) + ENDIF + IF(ICONT.EQ.1) THEN + CALL PSSTRK(ISPSP,WLFAC*WLINE,KSS,KSR) + ENDIF + ENDIF + DO JR=NAR,1,-1 + IAREG=(NAR+1)*(IREG-1)+JR + JREG=ABS(IDREG(IAREG)) + IF(IPRINT .GE. 100) THEN + WRITE(IOUT,6020) JR,IAREG,JREG,DCMESH(JR,4) + ENDIF +*---- +* Annular region +*---- + IF(JREG .NE. 0) THEN + RCIRC=REAL(FACT*DCMESH(JR,4)) +*---- +* Move cursor to center of annulus +*---- + CALL PSMOVE(ISPSP,CENTER,-3) +*---- +* Color and trace result +*---- + CALL PSPRAI(NINT,NPTS,XYPOS,CENTER,RCIRC, + > NSEG,IORDER,RADANG) + CALL PSDRAI(ISPSP,NSEG,IORDER,CENTER,RADANG) + IF(ICOL. GT. 0) THEN + CALL PSFILL(ISPSP,ICOL,COLREG(1,JREG),KFS,KFR) + ENDIF + IF(ICONT.EQ.1) THEN + CALL PSSTRK(ISPSP,WLFAC*WLINE,KSS,KSR) + ENDIF +*---- +* Return cursor to original position +*---- + CALL PSMOVE(ISPSP,CENTEB,-3) + ENDIF + ENDDO + ENDDO + ENDDO +*---- +* Pins +*---- + IF(NTPIN .GT. 0) THEN + NPTS=4 + ILEV=2 + MESHP(1)=1 + MESHP(2)=1 + MESHP(3)=1 + MESHP(4)=1 + CENTEP(1)=0.0 + CENTEP(2)=0.0 + DO IPIN=1,NTPIN +*---- +* Locate pin position +*---- + COSDIR(1)=DRAPIN(0,IPIN)*COS(DRAPIN(-1,IPIN)) + COSDIR(2)=DRAPIN(0,IPIN)*SIN(DRAPIN(-1,IPIN)) + CENTED(1)=CENTER(1)+REAL(FACT*COSDIR(1)) + CENTED(2)=CENTER(2)+REAL(FACT*COSDIR(2)) + RCIRC=REAL(FACT*DRAPIN(4,IPIN)) + ROTAX=PI/DTWO-DRAPIN(-1,IPIN) +*---- +* Move cursor to center of pin +*---- + CALL PSMOVE(ISPSP,CENTED,-3) +*---- +* Read pin information +*---- + WRITE(NAMCEL,'(A1,I8.8)') CLEV(ILEV),ITPIN(2,IPIN) + NAMREC=NAMCEL//'DIM' + IEDIMP(:NSTATE)=0 + CALL LCMGET(IPTRK,NAMREC,IEDIMP) + MESHP(1)=IEDIMP(3) + MESHP(2)=IEDIMP(4) + MESHP(3)=IEDIMP(5) + MESHP(4)=IEDIMP(2) + NBRP=MESHP(4) + NAMREC=NAMCEL//'RID' + CALL LCMGET(IPTRK,NAMREC,IDREG) + DO IDIR=1,4 + NAMREC=NAMCEL//'SM'//CDIR(IDIR) + IF(MESHP(IDIR) .GT. 0) THEN + CALL LCMGET(IPTRK,NAMREC,DCMESH(-1,IDIR)) + ENDIF + ENDDO + DO IY=MESHP(2),1,-1 + DO IX=MESHP(1),1,-1 + ILXY=((IY-1)*MESHP(1)+(IX-1))*NBRP +*---- +* Cartesian region +*---- + XYPOS(1,1)=REAL(FACT*(DCMESH(IX-1,1))) + XYPOS(2,1)=REAL(FACT*(DCMESH(IY-1,2))) + XYPOS(1,2)=REAL(FACT*(DCMESH(IX,1))) + XYPOS(2,2)=XYPOS(2,1) + XYPOS(1,3)=XYPOS(1,2) + XYPOS(2,3)=REAL(FACT*(DCMESH(IY,2))) + XYPOS(1,4)=XYPOS(1,1) + XYPOS(2,4)=XYPOS(2,3) + DO IR=MESHP(4),1,-1 + ILOC=ILXY+IR + IREG=ABS(IDREG(ILOC)) + RCIRC=REAL(FACT*DCMESH(IR,4)) +*---- +* Annular pin regions +*---- + IF(IREG .NE. 0) THEN + CALL PSPRAI(NINT,NPTS,XYPOS,CENTEP,RCIRC, + > NSEG,IORDER,RADANG) +*---- +* Rotate pins intersection points +*---- + DO ISEG=1,NSEG + RADANG(2,ISEG)=RADANG(2,ISEG)-REAL(ROTAX) + ENDDO +*---- +* Color and trace result +*---- + CALL PSDRAI(ISPSP,NSEG,IORDER,CENTEP,RADANG) + IF(ICOL. GT. 0) THEN + CALL PSFILL(ISPSP,ICOL,COLREG(1,IREG),KFS,KFR) + ENDIF + IF(ICONT.EQ.1) THEN + CALL PSSTRK(ISPSP,WLFAC*WLINE,KSS,KSR) + ENDIF + ENDIF + ENDDO + ENDDO + ENDDO +*---- +* Return cursor to original position +*---- + CENTED(1)=-CENTED(1) + CENTED(2)=-CENTED(2) + CALL PSMOVE(ISPSP,CENTED,-3) + ENDDO + ENDIF + IF(IPRINT .GE. 100) THEN + WRITE(IOUT,6001) NAMSBR + ENDIF + RETURN +*---- +* Output formats +*---- + 6000 FORMAT('(* Output from --',A6,'-- follows ') + 6001 FORMAT(' Output from --',A6,'-- completed *)') + 6010 FORMAT(' Sector ',I5,3X,'Triangle ',I5,5X,'T region ',I5,5X, + >'G region ',I5,5X, + >'Region number',I5,5X,'Number of corners ',I5) + 6011 FORMAT(2F20.7) + 6020 FORMAT(' Annulus ',I5,3X,'A Region ',I5,5X, + >'Region number',I5,5X,'radius ',F20.10) + END -- cgit v1.2.3