PROGRAM ELEMENTI FINITI QUADRATI ISOPARAMETRICI
      PARAMETER(NODES=9)
      PARAMETER(NELEMS=4)
      PARAMETER(IFMAX=3)
      PARAMETER(ICNMAX=3)
      PARAMETER(MBAND=10)
C     DIMENSIONAMENTI VARIABILI DELLA STRUTTURA
      DIMENSION XY(2,NODES),NVERT(4,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
     $ICNXY(2,ICNMAX),CNXY(2,ICNMAX),
     $STRUTK(2*NODES,MBAND),FORCE(2*NODES),UV(2*NODES),SIG(32,NELEMS)
C     APERTURA DEI FILE
      DIMENSION D(3,3),B(3,8),DELEL(8),IPOINT(8),AK(8,8),BDELEL(3,1),
     $AKTEMP(8,8),SIGEL(3),C(4),E(4),BT(8,3),DB(3,8),AJ(2,2),AJINV(2,2)
      WRITE(*,*)' '
      WRITE(*,*)' PROGRAMMA FEM PER STRUTTURE IN TP O DP'
      WRITE(*,*)' --------------------------------------'
      WRITE(*,*)' '
c      open(1,file='mesh_iso4_comp.dat', STATUS='OLD')
c      open(2,file='output_iso4_comp.3.txt', STATUS='UNKNOWN')
c      open(1,file='mesh_iso4_flex.dat', STATUS='OLD', STATUS='UNKNOWN')
c      open(2,file='output_iso4_flex.3.txt', STATUS='UNKNOWN')
      open(1,file='mesh_iso4_verifica.dat', STATUS='OLD')
      open(2,file='output_iso4_verifica.3.txt', STATUS='UNKNOWN')      
      WRITE(*,*)' '
      WRITE(*,*)'DATI DI INPUT:'
      WRITE(*,*)' '
C     LETTURA DEI FILE DI DATI
      CALL READIN(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,FXY,
     $ICNMAX,ICNXY,CNXY)
      WRITE(*,*)' '
      WRITE(*,*)'INIZIO DELLA PROCEDURA DI CALCOLO:'
      WRITE(*,*)' '
!      CALL WIDTHBAND(NELEMS,NVERT,MBAND)       <--- Š commentato perchŠ fortran 77 non accetta dimensioni di matrici variabili. Quindi occorre definire MABAND come PARAMETER!
      CALL CALC(NODES,NELEMS,NVERT,IFMAX,IFXY,FXY,ICNMAX,ICNXY,CNXY,YM,
     $PR,ITDP,XY,MBAND,FORCE)
      WRITE(*,*)' '
      WRITE(*,*)'IL PROGRAMMA HA COMPLETATO I CALCOLI:'
      WRITE(*,*)' '
      WRITE(*,*)' I RISULTATI SI TROVANO NEL FILE: RISULTATI_ISOQ.OUT'
      WRITE(*,*)' '
c      READ(*,*)
      STOP
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE READIN:
C
      SUBROUTINE READIN(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,
     $FXY,ICNMAX,ICNXY,CNXY)
      DIMENSION XY(2,NODES),NVERT(4,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
     $ICNXY(2,ICNMAX),CNXY(2,ICNMAX)
C     Lettura di YM, PR e ITDP:
      READ(1,*)YM,PR,ITDP
      WRITE(*,*)'YM,PR,ITDP=',YM,PR,ITDP
      WRITE(2,*)'YM,PR,ITDP=',YM,PR,ITDP
C     Lettura delle coordinate nodali (x,y):
      WRITE(*,*)'COORDINATE NODALI'
      WRITE(2,*)'COORDINATE NODALI'
      DO 10,I10=1,NODES
      READ(1,*)NODE,XY(1,I10),XY(2,I10)
      IF(I10.NE.NODE)THEN
      WRITE(*,*)'ERRORE DI INDIRIZZO DEL NODO',NODE
      WRITE(*,*)'DOVREBBE ESSERE',I10
      END IF
      WRITE(*,*)NODE,XY(1,I10),XY(2,I10)
      WRITE(2,*)NODE,XY(1,I10),XY(2,I10)
   10 CONTINUE
C     Lettura dei vertici degli elementi quadrati isoparametrici:
      WRITE(*,*)'NOMI DEI VERTICI DELL*ELEMENTO'
      WRITE(2,*)'NOMI DEI VERTICI DELL*ELEMENTO'
      DO 20,I20=1,NELEMS
      READ(1,*)NELEM,NVERT(1,I20),NVERT(2,I20),NVERT(3,I20),NVERT(4,I20)
      IF(I20.NE.NELEM)THEN
      WRITE(*,*)'ERRORE DI INDIRIZZO DELL*ELEMENTO',NELEM
      WRITE(*,*)'DOVREBBE ESSERE',I20
      END IF
      WRITE(*,*)NELEM,NVERT(1,I20),NVERT(2,I20),NVERT(3,I20),
     $NVERT(4,I20)
      WRITE(2,*)NELEM,NVERT(1,I20),NVERT(2,I20),NVERT(3,I20),
     $NVERT(4,I20)
   20 CONTINUE
C     Lettura delle forze nodali:
      WRITE(*,*)'FORZE NODALI'
      WRITE(2,*)'FORZE NODALI'
      DO 30,I30=1,IFMAX
      READ(1,*)IFXY(I30),FXY(1,I30),FXY(2,I30)
      WRITE(*,*)IFXY(I30),FXY(1,I30),FXY(2,I30)
      WRITE(2,*)IFXY(I30),FXY(1,I30),FXY(2,I30)
   30 CONTINUE
C     Lettura delle condizioni di vincolamento:
      WRITE(*,*)'VINCOLAMENTI'
      WRITE(2,*)'VINCOLAMENTI'
      DO 40,I40=1,ICNMAX
      READ(1,*)ICNXY(1,I40),ICNXY(2,I40),CNXY(1,I40),CNXY(2,I40)
      WRITE(*,*)ICNXY(1,I40),ICNXY(2,I40),CNXY(1,I40),CNXY(2,I40)
      WRITE(2,*)ICNXY(1,I40),ICNXY(2,I40),CNXY(1,I40),CNXY(2,I40)
   40 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE WIDTHBAND:
C
!      SUBROUTINE WIDTHBAND(NELEMS,NVERT,MBAND)
!      DIMENSION NVERT(4,NELEMS),IPOINT(8)
!      MBAND=0.                             !Metto MBAND=1 perchŠ Š una
!      DO 10,I10=1,NELEMS                  !variabile intera, e nei PC
!      IPOINT(1)=NVERT(1,I10)*2-1          !non esiste lo ZERO assoluto.
!      IPOINT(2)=NVERT(1,I10)*2
!      IPOINT(3)=NVERT(2,I10)*2-1
!      IPOINT(4)=NVERT(2,I10)*2
!      IPOINT(5)=NVERT(3,I10)*2-1
!      IPOINT(6)=NVERT(3,I10)*2
!      IPOINT(7)=NVERT(4,I10)*2-1
!      IPOINT(8)=NVERT(4,I10)*2
!      IPONTMAX=MAX(IPOINT(2),IPOINT(4),IPOINT(6),IPOINT(8))
!      IPOINTMIN=MIN(IPOINT(1),IPOINT(3),IPOINT(5),IPOINT(7))
!      MBANDTEMP=2*(IPOINTMAX-IPOINTMIN+1)
!      IF (MBANDTEMP.GT.MBAND) THEN
!      MBAND=MBANDTEMP
!      END IF
!   10 CONTINUE
!      WRITE(*,*)'LA LARGHEZZA DI BANDA DELLA MATRICE DI RIGIDEZZA'
!      WRITE(*,*)'GLOBALE DELLA STRUTTURA VALE: ',MBAND
!      WRITE(2,*)'LA LARGHEZZA DI BANDA DELLA MATRICE DI RIGIDEZZA'
!      WRITE(2,*)'GLOBALE DELLA STRUTTURA VALE: ',MBAND
!      RETURN
!      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUROUTINE CLEAR:
C
      SUBROUTINE CLEAR(K1,K2,A)
      DIMENSION A(K1,K2)
      DO 10,I10=1,K1
      DO 20,I20=1,K2
      A(I10,I20)=0.
   20 CONTINUE
   10 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE DMAT:
C
      SUBROUTINE DMAT(YM,PR,ITDP,D)
      DIMENSION D(3,3)
      IF(ITDP.EQ.0)GO TO 100
      IF(ITDP.EQ.1)GO TO 200
C     Per Stato di Tensione Piana:
  100 CONTINUE
      CTP=YM/(1-PR**2)
      G=YM/(2*(1+PR))
C     Struttura della matrice D per Stato di Tensione Piana:
C                     | 1*CTP    PR*CTP   0 |
C                D =  | PR*CTP   1*CTP    0 |
C                     | 0          0      G |
      D(1,1)=1.
      D(1,2)=PR
      D(2,1)=PR
      D(2,2)=1.
      DO 10,I10=1,2
      DO 20,I20=1,2
      D(I10,I20)=D(I10,I20)*CTP
   20 CONTINUE
   10 CONTINUE
      D(1,3)=0.
      D(2,3)=0.
      D(3,1)=0.
      D(3,2)=0.
      D(3,3)=G
      RETURN
C     Per Stato di Deformazione Piana:
  200 CONTINUE
      CDP=YM*(1-PR)/((1+PR)*(1-2*PR))
      CPR=PR/(1-PR)
      G=YM/((1+PR)*2)
C     Struttura della matrice D per Stato di Deformazione Piana:
C                     | 1*CDP     CDP*CPR   0 |
C                D =  | CDP*CPR    1*CDP    0 |
C                     | 0            0      G |
      D(1,1)=1.
      D(1,2)=CPR
      D(2,1)=CPR
      D(2,2)=1.
      DO 30,I30=1,2
      DO 40,I40=1,2
      D(I30,I40)=D(I30,I40)*CDP
   40 CONTINUE
   30 CONTINUE
      D(1,3)=0.
      D(2,3)=0.
      D(3,1)=0.
      D(3,2)=0.
      D(3,3)=G
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE KISOQ:
C
      SUBROUTINE KISOQ(YM,PR,X1,Y1,X2,Y2,X3,Y3,X4,Y4,D,B,AK)
C     DETERMINAZIONE DELLA MATRICE DI RIGIDEZZA ELK PER
C     L*ELEMENTO ISOPARAMETRICO QUADRILATERO
      DIMENSION D(3,3),B(3,8),BT(8,3),DB(3,8),AKTEMP(8,8),AK(8,8)
      CALL DMAT(YM,PR,ITDP,D)
      CALL CLEAR(8,8,AK)
C     INTEGRAZIONE GAUSSIANA A QUATTRO PUNTI.
C     I PESI DI GAUSS SONO TUTTI UNITARI
C     DETERMINAZIONE DI B PER CSI=-0.57735, ETA=-0.57735
      CSI=-0.57735
      ETA=-0.57735
      CALL BISOQ(X1,Y1,X2,Y2,X3,Y3,X4,Y4,CSI,ETA,B,DETJ)
      CALL TRSPMAT(3,8,B,BT)
      CALL PRODMAT(3,3,8,D,B,DB)
C     AKTEMP=MATRICE DI RIGIDEZZA K TEMPORANEA
      CALL PRODMAT(8,3,8,BT,DB,AKTEMP)
      CALL PRDMTSC(AKTEMP,DETJ)
C     AGGIORNAMENTO DELLA MATRICE DI RIGIDEZZA K
      CALL AGGIORN(AKTEMP,AK)
C     DETERMINAZIONE DI B PER CSI=-0.57735,ETA=+0.57735
      CSI=-0.57735
      ETA=+0.57735
      CALL BISOQ(X1,Y1,X2,Y2,X3,Y3,X4,Y4,CSI,ETA,B,DETJ)
      CALL TRSPMAT(3,8,B,BT)
      CALL PRODMAT(3,3,8,D,B,DB)
C     AKTEMP=MATRICE DI RIGIDEZZA K TEMPORANEA
      CALL PRODMAT(8,3,8,BT,DB,AKTEMP)
      CALL PRDMTSC(AKTEMP,DETJ)
C     AGGIORNAMENTO DELLA MATRICE DI RIGIDEZZA K
      CALL AGGIORN(AKTEMP,AK)
C     DETERMINAZIONE DI B PER CSI=+0.57735,ETA=+0.57735
      CSI=+0.57735
      ETA=+0.57735
      CALL BISOQ(X1,Y1,X2,Y2,X3,Y3,X4,Y4,CSI,ETA,B,DETJ)
      CALL TRSPMAT(3,8,B,BT)
      CALL PRODMAT(3,3,8,D,B,DB)
C     AKTEMP= MATRICE DI RIGIDEZZA K TEMPORANEA
      CALL PRODMAT(8,3,8,BT,DB,AKTEMP)
      CALL PRDMTSC(AKTEMP,DETJ)
C     AGGIORNAMENTO DELLA MATRICE DI RIGIDEZZA K
      CALL AGGIORN(AKTEMP,AK)
C     DETERMINAZIONE DI B PER CSI=+0.57735,ETA=-0.57735
      CSI=+0.57735
      ETA=-0.57735
      CALL BISOQ(X1,Y1,X2,Y2,X3,Y3,X4,Y4,CSI,ETA,B,DETJ)
      CALL TRSPMAT(3,8,B,BT)
      CALL PRODMAT(3,3,8,D,B,DB)
C     AKTEMP=MATRICE DI RIGIDEZZA K TEMPORANEA
      CALL PRODMAT(8,3,8,BT,DB,AKTEMP)
      CALL PRDMTSC(AKTEMP,DETJ)
C     AGGIORNAMENTO DELLA MATRICE DI RIGIDEZZA K
      CALL AGGIORN(AKTEMP,AK)
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE DIFFUV:
C
      SUBROUTINE DIFFUV(U1,V1,U2,V2,U3,V3,U4,V4,CSI,ETA,
     $DUDCSI,DVDCSI,DUDETA,DVDETA)
C     CALCOLO DI DUDCSI,DUDETA,DVDCSI,DVDETA
      DUDCSI=-(1-ETA)*U1/4+(1-ETA)*U2/4+(1+ETA)*U3/4-
     $(1+ETA)*U4/4
      DUDETA=-(1-CSI)*U1/4-(1+CSI)*U2/4+(1+CSI)*U3/4+
     $(1-CSI)*U4/4
      DVDCSI=-(1-ETA)*V1/4+(1-ETA)*V2/4+(1+ETA)*V3/4-
     $(1+ETA)*V4/4
      DVDETA=-(1-CSI)*V1/4-(1+CSI)*V2/4+(1+CSI)*V3/4+
     $(1-CSI)*V4/4
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE BISOQ:
C
      SUBROUTINE BISOQ(X1,Y1,X2,Y2,X3,Y3,X4,Y4,CSI,ETA,B,DETJ)
C     RIEMPIMENTO MATRICE B(3,8) PER
C     L*ELEMENTO ISOPARAMETRICO QUADRILATERO
C     AJ=MATRICE J
C     AJINV=J^(-1)
C     B=MATRICE B(3,8)
      DIMENSION AJ(2,2),AJINV(2,2),B(3,8)
      AJ(1,1)=-(1-ETA)*X1/4+(1-ETA)*X2/4+(1+ETA)*X3/4-(1+ETA)*X4/4
      AJ(1,2)=-(1-ETA)*Y1/4+(1-ETA)*Y2/4+(1+ETA)*Y3/4-(1+ETA)*Y4/4
      AJ(2,1)=-(1-CSI)*X1/4-(1+CSI)*X2/4+(1+CSI)*X3/4+(1-CSI)*X4/4
      AJ(2,2)=-(1-CSI)*Y1/4-(1+CSI)*Y2/4+(1+CSI)*Y3/4+(1-CSI)*Y4/4
      DETJ=AJ(1,1)*AJ(2,2)-AJ(2,1)*AJ(1,2)
      AJINV(1,1)=AJ(2,2)/DETJ
      AJINV(1,2)=-AJ(1,2)/DETJ
      AJINV(2,1)=-AJ(2,1)/DETJ
      AJINV(2,2)=AJ(1,1)/DETJ
C     CASO U1=1 E GLI ALTRI SPOSTAMENTI U E V NULLI
      CALL DIFFUV(1.,0.,0.,0.,0.,0.,0.,0.,CSI,ETA,DUDCSI,DVDCSI,
     $DUDETA,DVDETA)
      DUDX=AJINV(1,1)*DUDCSI+AJINV(1,2)*DUDETA
      DUDY=AJINV(2,1)*DUDCSI+AJINV(2,2)*DUDETA
      DVDX=AJINV(1,1)*DVDCSI+AJINV(1,2)*DVDETA
C     DVDY=AJINV(2,1)*DVDCSI+AJINV(2,2)*DVDETA
      B(1,1)=DUDX
C     B(2,1)=DVDY
      B(2,1)=0.
      B(3,1)=DUDY+DVDX
C     CASO V1=1 E GLI ALTRI SPOSTAMENTI U E V NULLI
      CALL DIFFUV(0.,1.,0.,0.,0.,0.,0.,0.,CSI,ETA,
     $DUDCSI,DVDCSI,DUDETA,DVDETA)
C     DUDX=AJINV(1,1)*DUDCSI+AJINV(1,2)*DUDETA
      DUDY=AJINV(2,1)*DUDCSI+AJINV(2,2)*DUDETA
      DVDX=AJINV(1,1)*DVDCSI+AJINV(1,2)*DVDETA
      DVDY=AJINV(2,1)*DVDCSI+AJINV(2,2)*DVDETA
C     B(1,2)=DUDX
      B(1,2)=0.
      B(2,2)=DVDY
      B(3,2)=DUDY+DVDX
C     CASO U2=1 E GLI ALTRI SPOSTAMENTI U E V NULLI
      CALL DIFFUV(0.,0.,1.,0.,0.,0.,0.,0.,CSI,ETA,
     $DUDCSI,DVDCSI,DUDETA,DVDETA)
      DUDX=AJINV(1,1)*DUDCSI+AJINV(1,2)*DUDETA
      DUDY=AJINV(2,1)*DUDCSI+AJINV(2,2)*DUDETA
      DVDX=AJINV(1,1)*DVDCSI+AJINV(1,2)*DVDETA
C     DVDY=AJINV(2,1)*DVDCSI+AJINV(2,2)*DVDETA
      B(1,3)=DUDX
C     B(2,3)=DVDY
      B(2,3)=0.
      B(3,3)=DUDY+DVDX
C     CASO V2=1 E GLI ALTRI SPOSTAMENTI U E V NULLI
      CALL DIFFUV(0.,0.,0.,1.,0.,0.,0.,0.,CSI,ETA,
     $DUDCSI,DVDCSI,DUDETA,DVDETA)
C     DUDX=AJINV(1,1)*DUDCSI+AJINV(1,2)*DUDETA
      DUDY=AJINV(2,1)*DUDCSI+AJINV(2,2)*DUDETA
      DVDX=AJINV(1,1)*DVDCSI+AJINV(1,2)*DVDETA
      DVDY=AJINV(2,1)*DVDCSI+AJINV(2,2)*DVDETA
C     B(1,4)=DUDX
      B(1,4)=0.
      B(2,4)=DVDY
      B(3,4)=DUDY+DVDX
C     CASO U3=1 E GLI ALTRI SPOSTAMENTI U E V NULLI
      CALL DIFFUV(0.,0.,0.,0.,1.,0.,0.,0.,CSI,ETA,
     $DUDCSI,DVDCSI,DUDETA,DVDETA)
      DUDX=AJINV(1,1)*DUDCSI+AJINV(1,2)*DUDETA
      DUDY=AJINV(2,1)*DUDCSI+AJINV(2,2)*DUDETA
      DVDX=AJINV(1,1)*DVDCSI+AJINV(1,2)*DVDETA
C     DVDY=AJINV(2,1)*DVDCSI+AJINV(2,2)*DVDETA
      B(1,5)=DUDX
C     B(2,5)=DVDY
      B(2,5)=0.
      B(3,5)=DUDY+DVDX
C     CASO V3=1 E GLI ALTRI SPOSTAMENTI U E V NULLI
      CALL DIFFUV(0.,0.,0.,0.,0.,1.,0.,0.,CSI,ETA,
     $DUDCSI,DVDCSI,DUDETA,DVDETA)
C     DUDX=AJINV(1,1)*DUDCSI+AJINV(1,2)*DUDETA
      DUDY=AJINV(2,1)*DUDCSI+AJINV(2,2)*DUDETA
      DVDX=AJINV(1,1)*DVDCSI+AJINV(1,2)*DVDETA
      DVDY=AJINV(2,1)*DVDCSI+AJINV(2,2)*DVDETA
C     B(1,6)=DUDX
      B(1,6)=0.
      B(2,6)=DVDY
      B(3,6)=DUDY+DVDX
C     CASO U4=1 E GLI ALTRI SPOSTAMENTI U E V NULLI
      CALL DIFFUV(0.,0.,0.,0.,0.,0.,1.,0.,CSI,ETA,
     $DUDCSI,DVDCSI,DUDETA,DVDETA)
      DUDX=AJINV(1,1)*DUDCSI+AJINV(1,2)*DUDETA
      DUDY=AJINV(2,1)*DUDCSI+AJINV(2,2)*DUDETA
      DVDX=AJINV(1,1)*DVDCSI+AJINV(1,2)*DVDETA
C     DVDY=AJINV(2,1)*DVDCSI+AJINV(2,2)*DVDETA
      B(1,7)=DUDX
C     B(2,7)=DVDY
      B(2,7)=0.
      B(3,7)=DUDY+DVDX
C     CASO V4=1 E GLI ALTRI SPOSTAMENTI U E V NULLI
      CALL DIFFUV(0.,0.,0.,0.,0.,0.,0.,1.,CSI,ETA,
     $DUDCSI,DVDCSI,DUDETA,DVDETA)
C     DUDX=AJINV(1,1)*DUDCSI+AJINV(1,2)*DUDETA
      DUDY=AJINV(2,1)*DUDCSI+AJINV(2,2)*DUDETA
      DVDX=AJINV(1,1)*DVDCSI+AJINV(1,2)*DVDETA
      DVDY=AJINV(2,1)*DVDCSI+AJINV(2,2)*DVDETA
C     B(1,8)=DUDX
      B(1,8)=0.
      B(2,8)=DVDY
      B(3,8)=DUDY+DVDX
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE TRSPMAT: Calcolo della matrice trasposta AT di una
C     matrice A.
C
      SUBROUTINE TRSPMAT(K1,K2,A,AT)
      DIMENSION A(K1,K2),AT(K2,K1)
      DO 10,I10=1,K1
      DO 20,I20=1,K2
      AT(I20,I10)=A(I10,I20)
   20 CONTINUE
   10 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE PRODMAT: Esecuzione del prodotto AB di due matrici
C     A e B.
C
      SUBROUTINE PRODMAT(K1,K2,K3,A,B,AB)
      DIMENSION A(K1,K2),B(K2,K3),AB(K1,K3)
      DO 10,I10=1,K1
      DO 30,I30=1,K3
      AB(I10,I30)=0.   !Inizializzazione a zero della matrice prodotto.
      DO 20,I20=1,K2
      AB(I10,I30)=AB(I10,I30)+A(I10,I20)*B(I20,I30)   !Prodotto.
   20 CONTINUE
   30 CONTINUE
   10 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE PRODMATSC:
C
      SUBROUTINE PRDMTSC(AKTEMP,DETJ)
C     PRODOTTO TRA MATRICE E SCALARE
      DIMENSION AKTEMP(8,8)
      DO 10,I10=1,8
      DO 20,I20=1,8
      AKTEMP(I10,I20)=AKTEMP(I10,I20)*DETJ
   20 CONTINUE
   10 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE AGGIORN:
C
      SUBROUTINE AGGIORN(AKTEMP,AK)
C     AGGIORNAMENTO DELLA MATRICE DI RIGIDEZZA
      DIMENSION AKTEMP(8,8),AK(8,8)
      DO 10,I10=1,8
      DO 20,I20=1,8
      AK(I10,I20)=AK(I10,I20)+AKTEMP(I10,I20)
   20 CONTINUE
   10 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE POINTVT: Costruzione del vettore puntatore IPOINT
C     necessario ad effettuare l'assemblaggio della matrice ELK di
C     rigidezza di ciascun elemento finito triangolare nella matrice
C     STRUTK di rigidezza globale della struttura.
C
      SUBROUTINE POINTVT(NELEMS,NELEM,NVERT,IPOINT)
      DIMENSION NVERT(4,NELEMS),IPOINT(8)
      IPOINT(1)=NVERT(1,NELEM)*2-1
      IPOINT(2)=NVERT(1,NELEM)*2
      IPOINT(3)=NVERT(2,NELEM)*2-1
      IPOINT(4)=NVERT(2,NELEM)*2
      IPOINT(5)=NVERT(3,NELEM)*2-1
      IPOINT(6)=NVERT(3,NELEM)*2
      IPOINT(7)=NVERT(4,NELEM)*2-1
      IPOINT(8)=NVERT(4,NELEM)*2
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE ASSEMBL:
C
      SUBROUTINE ASSEMBL(AK,NODES,IPOINT,MBAND,STRUTK)
      DIMENSION AK(8,8),IPOINT(8),STRUTK(2*NODES,MBAND)
      DO 10,I10=1,8
      DO 20,I20=1,8
      IND10=IPOINT(I10)
      IND20=IPOINT(I20)
      IND30=IPOINT(I20)-IPOINT(I10)+1
      IF(IND20.GE.IND10)THEN
      STRUTK(IND10,IND30)=STRUTK(IND10,IND30)+AK(I10,I20)
      END IF
   20 CONTINUE
   10 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE FORCES: Riempimento del vettore globale delle forze
C     nodali FORCE.
C
      SUBROUTINE FORCES(IFMAX,IFXY,FXY,NODES,FORCE)
      DIMENSION IFXY(IFMAX),FXY(2,IFMAX),FORCE(2*NODES)
      DO 10,I10=1,IFMAX
      INDX=IFXY(I10)*2-1
      INDY=IFXY(I10)*2
      FORCE(INDX)=FXY(1,I10)
      FORCE(INDY)=FXY(2,I10)
   10 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE CNSTNG:
C
      SUBROUTINE CNSTNG(ICNMAX,ICNXY,CNXY,NODES,FORCE,MBAND,STRUTK)
C     SUBROUTINE DI VINCOLAMENTO
C     ICNMAX= NUMERO MASSIMO DI NODI VINCOLATI
C     ICNXY= MATRICE CON IL NOME DEI NODI VINCOLATI ED IL ÍIPO DI VINCOLO
C     1: VINCOLATO LUNGO X
C     2: VINCOLATO LUNGO Y
C     3: VINCOLATO LUNGO X ED Y (CERNIERA)
C     CNXY== MATRICE CON GLI SPOSTAMENTÍ ASSEGNATI LUNGO X E Y
C     FORCE= VETTORE TERMINI NOTI
C     STRUTK= MATRICE DI RIGIDEZZA DELLA STRUTTURA
      DIMENSION ICNXY(2,ICNMAX),CNXY(2,ICNMAX),STRUTK(2*NODES,MBAND),
     $FORCE(2*NODES)
      DO 10,I10=1,ICNMAX
      IF(ICNXY(2,I10).EQ.1) GOTO 100
      IF(ICNXY(2,I10).EQ.2) GOTO 200
      IF(ICNXY(2,I10).EQ.3) GOTO 300
      WRITE(*,*)'ERRORE NEL VINCOLAMENTO DEL NODO',I10
      WRITE(*,*)'L*INDICE DI VINCOLAMENTO DEVE ESSERE 1, 2 o 3'
      STOP
  100 CONTINUE
      DO 110,I110=1,MBAND
      FORCE(2*ICNXY(1,I10)-1+I110-1)=FORCE(2*ICNXY(1,I10)-1+I110-1)-
     $STRUTK(2*ICNXY(1,I10)-1,I110)*CNXY(1,I10)
      STRUTK(2*ICNXY(1,I10)-1,I110)=0.
  110 CONTINUE
      IF (2*ICNXY(1,I10)-1.LE.MBAND) THEN
      DO 120,I120=1,2*ICNXY(1,I10)-1
      FORCE(2*ICNXY(1,I10)-1-I120+1)=FORCE(2*ICNXY(1,I10)-1-I120+1)-
     $STRUTK(2*ICNXY(1,I10)-1-I120+1,I120)*CNXY(1,I10)
      STRUTK(2*ICNXY(1,I10)-1-I120+1,I120)=0.
  120 CONTINUE
      ELSE
      DO 130,I130=1,MBAND
      FORCE(2*ICNXY(1,I10)-1-I130+1)=FORCE(2*ICNXY(1,I10)-1-I130+1)-
     $STRUTK(2*ICNXY(1,I10)-1-I130+1,I130)*CNXY(1,I10)
      STRUTK(2*ICNXY(1,I10)-1-I130+1,I130)=0.
  130 CONTINUE
      ENDIF
      STRUTK(2*ICNXY(1,I10)-1,1)=1.
      FORCE(2*ICNXY(1,I10)-1)=CNXY(1,I10)
      GOTO 10
  200 CONTINUE
      DO 210,I210=1,MBAND
      FORCE(2*ICNXY(1,I10)+I210-1)=FORCE(2*ICNXY(1,I10)+I210-1)-
     $STRUTK(2*ICNXY(1,I10),I210)*CNXY(2,I10)
      STRUTK(2*ICNXY(1,I10),I210)=0.
  210 CONTINUE
      IF(2*ICNXY(1,I10).LE.MBAND)THEN
      DO 220,I220=1,2*ICNXY(1,I10)
      FORCE(2*ICNXY(1,I10)-I220+1)=FORCE(2*ICNXY(1,I10)-I220+1)-
     $STRUTK(2*ICNXY(1,I10)-I220+1,I220)*CNXY(2,I10)
      STRUTK(2*ICNXY(1,I10)-I220+1,I220)=0.
  220 CONTINUE
      ELSE
      DO 230,I230=1,MBAND
      FORCE(2*ICNXY(1,I10)-I230+1)=FORCE(2*ICNXY(1,I10)-I230+1)-
     $STRUTK(2*ICNXY(1,I10)-I230+1,I230)*CNXY(2,I10)
      STRUTK(2*ICNXY(1,I10)-I230+1,I230)=0.
  230 CONTINUE
      ENDIF
      STRUTK(2*ICNXY(1,I10),1)=1.
      FORCE(2*ICNXY(1,I10))=CNXY(2,I10)
      GOTO 10
  300 CONTINUE
      DO 310,I310=1,MBAND
      FORCE(2*ICNXY(1,I10)-1+I310-1)=FORCE(2*ICNXY(1,I10)-1+I310-1)-
     $STRUTK(2*ICNXY(1,I10)-1,I310)*CNXY(1,I10)
      FORCE(2*ICNXY(1,I10)+I310-1)=FORCE(2*ICNXY(1,I10)+I310-1)-
     $STRUTK(2*ICNXY(1,I10),I310)*CNXY(2,I10)
      STRUTK(2*ICNXY(1,I10)-1,I310)=0.
      STRUTK(2*ICNXY(1,I10),I310)=0.
  310 CONTINUE
      IF (2*ICNXY(1,I10)-1.LE.MBAND) THEN
      DO 320,I320=1,2*ICNXY(1,I10)-1
      FORCE(2*ICNXY(1,I10)-1-I320+1)=FORCE(2*ICNXY(1,I10)-1-I320+1)-
     $STRUTK(2*ICNXY(1,I10)-1-I320+1,I320)*CNXY(1,I10)
      STRUTK(2*ICNXY(1,I10)-1-I320+1,I320)=0.
  320 CONTINUE
      ELSE
      DO 330,I330=1,MBAND
      FORCE(2*ICNXY(1,I10)-1-I330+1)=FORCE(2*ICNXY(1,I10)-1-I330+1)-
     $STRUTK(2*ICNXY(1,I10)-1-I330+1,I330)*CNXY(1,I10)
      STRUTK(2*ICNXY(1,I10)-1-I330+1,I330)=0.
  330 CONTINUE
      ENDIF
      STRUTK(2*ICNXY(1,I10)-1,1)=1.
      FORCE(2*ICNXY(1,I10)-1)=CNXY(1,I10)
      IF (2*ICNXY(1,I10).LE.MBAND) THEN
      DO 420,I420=1,2*ICNXY(1,I10)
      FORCE(2*ICNXY(1,I10)-I420+1)=FORCE(2*ICNXY(1,I10)-I420+1)-
     $STRUTK(2*ICNXY(1,I10)-I420+1,I420)*CNXY(2,I10)
      STRUTK(2*ICNXY(1,I10)-I420+1,I420)=0.
  420 CONTINUE
      ELSE
      DO 430,I430=1,MBAND
      FORCE(2*ICNXY(1,I10)-I430+1)=FORCE(2*ICNXY(1,I10)-I430+1)-
     $STRUTK(2*ICNXY(1,I10)-I430+1,I430)*CNXY(2,I10)
      STRUTK(2*ICNXY(1,I10)-I430+1,I430)=0.
  430 CONTINUE
      ENDIF
      STRUTK(2*ICNXY(1,I10),1)=1.
      FORCE(2*ICNXY(1,I10))=CNXY(2,I10)
   10 CONTINUE
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE COOK:
C
      SUBROUTINE COOK2(NSIZE,MBAND,S,R)
      DIMENSION S(NSIZE,MBAND),R(NSIZE)
      !RIDUZIONE IN FORMA TRIANGOLARE SUPERIORE
      DO 110,IBASE=1,NSIZE-1                    !CICLO PRINCIPALE PER AZZERARE I TERMINI SOTTODIAGONALI DELLA COLONNA "IBASE"
      DO 120,ISCORRI=1,MIN(MBAND-1,NSIZE-IBASE) !CICLO SECONDARIO SULLE RIGHE IBASE+1, IBASE+2, IBASE+3...
C                                              !FINO AD ESAURIMENTO BANDA IN IBASE+MBAND-1 O ESAURIMENTO MATRICE IN IBASE + (NSIZE-IBASE)
      COEFF=-S(IBASE,1+ISCORRI)/S(IBASE,1)
      DO 130,JSCORRI=1,MBAND-ISCORRI            !CICLO TERZIARIO SUGLI ELEMENTI DELLA RIGA ISCORRI, DA 1 A MBAND-1, POI DA 1 A MBAND-2,
C                                               !... INFINE DA 1 A MBAND-(MBAND-1)=1
      S(IBASE+ISCORRI,JSCORRI)=S(IBASE+ISCORRI,JSCORRI)+
     $COEFF*S(IBASE,ISCORRI+JSCORRI)
C                                               !L'ELEMENTO DELLA RIGA IBASE ASSOCIATO È NELLA BANDATA SCOSTATO DI ISCORRI COLONNE
  130 CONTINUE
      R(IBASE+ISCORRI)=R(IBASE+ISCORRI)+COEFF*R(IBASE)   !NON SCORDIAMOCI DI AGGIORNARE IL TERMINE NOTO.
  120 CONTINUE
  110 CONTINUE
C     !BACKSUBSTITUTION
      DO 210,IBASE=NSIZE,1,-1      !SCORRO A ROVESCIO IL VETTORE DELLE INCOGNITE, SOVRASCRITTO VIA VIA AL TERMINE NOTO R
      SOMMA = 0.
      DO 220,ISCORRI=1,MIN(NSIZE-IBASE,MBAND-1) !SCORRO SUI TERMINI GIA' RISOLTI
      SOMMA=SOMMA+R(IBASE+ISCORRI)*S(IBASE,1+ISCORRI)
  220 CONTINUE
       R(IBASE) = (R(IBASE) - SOMMA) / S(IBASE,1)
  210 CONTINUE
      END


        SUBROUTINE COOK(NSIZE,MBAND,S,R)

      DIMENSION S(NSIZE,MBAND),R(NSIZE)

700       DO 790 N=1,NSIZE

      DO 780 L=2,MBAND

      IF(S(N,L).EQ.0)GO TO 780

      I=N+L-1

      C=S(N,L)/S(N,1)

      J=0

      DO 750 K=L,MBAND

      J=J+1

      IF(J.GT.MBAND) THEN

      PRINT*,'AUMENTA MBAND!!!!'

      STOP

      ENDIF

750       S(I,J)=S(I,J)-C*S(N,K)

      S(N,L)=C

780       CONTINUE

790       CONTINUE

800       DO 830 N=1,NSIZE

      DO 820 L=2,MBAND

      IF (S(N,L).EQ. 0) GO TO 820

      I=N+L-1

      R(I)=R(I)-S(N,L)*R(N)

820       CONTINUE

830       R(N)=R(N)/S(N,1)

      DO 860 M=2,NSIZE

      N=NSIZE+1-M

      DO 850 L=2,MBAND

      IF (S(N,L).EQ.0) GO TO 850

      K=N+L-1

      R(N)=R(N)-S(N,L)*R(K)

850       CONTINUE

860       CONTINUE

        RETURN

      END

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE TENSMAT:
C
      SUBROUTINE TENSMAT(PR,ITDP,NODES,NELEMS,XY,NVERT,B,D,UV,SIG)
      DIMENSION XY(2,NODES),NVERT(4,NELEMS),B(3,8),D(3,3),UV(2*NODES),
     1DELEL(8),BDELEL(3,1),SIGEL(3),SIG(32,NELEMS),C(4),E(4)
      SMAXMISES=0.
      DO 10,I10=1,NELEMS
      X1=XY(1,NVERT(1,I10))
      Y1=XY(2,NVERT(1,I10))
      X2=XY(1,NVERT(2,I10))
      Y2=XY(2,NVERT(2,I10))
      X3=XY(1,NVERT(3,I10))
      Y3=XY(2,NVERT(3,I10))
      X4=XY(1,NVERT(4,I10))
      Y4=XY(2,NVERT(4,I10))
      DELEL(1)=UV(2*NVERT(1,I10)-1)
      DELEL(2)=UV(2*NVERT(1,I10))
      DELEL(3)=UV(2*NVERT(2,I10)-1)
      DELEL(4)=UV(2*NVERT(2,I10))
      DELEL(5)=UV(2*NVERT(3,I10)-1)
      DELEL(6)=UV(2*NVERT(3,I10))
      DELEL(7)=UV(2*NVERT(4,I10)-1)
      DELEL(8)=UV(2*NVERT(4,I10))
C     PUNTI DI GAUSS IN CUI CALCOLARE LE TENSIONI IN SENSO ANTIORARIO
C     CSI
      C(1)=-0.57735
      C(2)=+0.57735
      C(3)=+0.57735
      C(4)=-0.57735
C     ETA
      E(1)=-0.57735
      E(2)=-0.57735
      E(3)=+0.57735
      E(4)=+0.57735
C     INDICE DI GAUSS
      IG=0
      DO 20,I20=1,4
      CSI=C(I20)
      ETA=E(I20)
      CALL BISOQ(X1,Y1,X2,Y2,X3,Y3,X4,Y4,CSI,ETA,B,DETJ)
      CALL PRODMAT(3,8,1,B,DELEL,BDELEL)
      CALL PRODMAT(3,3,1,D,BDELEL,SIGEL)
      SIG(1+IG,I10)=SIGEL(1) !SIGMA_X
      SIG(2+IG,I10)=SIGEL(2) !SIGMA_Y
      IF (ITDP.EQ.1) THEN      !SIGMA_Z PER I DUE CASI...
      SIG(3+IG,I10)=PR*(SIGEL(1)+SIGEL(2)) !DEFORMAZIONE PIANA
      ELSE
      SIG(3+IG,I10)=0. !TENSIONE PIANA
      END IF
      SIG(4+IG,I10)=SIGEL(3)   !T_XY
      SX=SIG(1+IG,I10)    !
      SY=SIG(2+IG,I10)    ! RI-NOMINAZIONE PER FACILITARE LE FORMULE
      SZ=SIG(3+IG,I10)    ! SUCCESSIVE...
      TXY=SIG(4+IG,I10)   !
      SIG(5+IG,I10)=(SX+SY)/2+SQRT(((SX-SY)**2)/4+TXY**2)  ! TENSIONE PRINCIPALE LUNGO X (SIGMA_1)
      SIG(6+IG,I10)=(SX+SY)/2-SQRT(((SX-SY)**2)/4+TXY**2)  ! TENSIONE PRINCIPALE LUNGO Y (SIGMA_2)
      SIG(7+IG,I10)=SZ                                     ! TENSIONE PRINCIPALE LUNGO Z (SIGMA_3)
      SIG(8+IG,I10)=SQRT(SX**2+SY**2+SZ**2-(SX*SY+SY*SZ+SX*SZ)+3*TXY**2)  ! TENSIONE IDEALE SECONDO VON MISES
      XGG=((1-CSI)*(1-ETA)*X1/4)+((1+CSI)*(1-ETA)*X2/4)+
     $((1+CSI)*(1+ETA)*X3/4)+((1-CSI)*(1+ETA)*X4/4)
      YGG=((1-CSI)*(1-ETA)*Y1/4)+((1+CSI)*(1-ETA)*Y2/4)+                                 !  '(A6,A6,10A15)'
     $((1+CSI)*(1+ETA)*Y3/4)+((1-CSI)*(1+ETA)*Y4/4)                                      !   '(I6,I6,10E15.5)'
      WRITE(2,*)
      WRITE(2,*)'TENSIONI NEGLI ELEMENTI'
      WRITE(2,*)'ELEMENTO','   P.INTEGR.','            X',
     $'                   Y','                   SIGMA_x',
     $'           SIGMA_y','             TAU_xy',
     $'          VON MISES'
      WRITE(2,*)I10,'   ',I20,'        ',XGG,YGG,SIG(1+IG,I10),'  ',
     $SIG(2+IG,I10),SIG(4+IG,I10),'  ',SIG(8+IG,I10)
      IF (SIG(8+IG,I10).GT.SMAXMISES) THEN
      SMAXMISES=SIG(8+IG,I10)
      NMAXMISES=I10
      ENDIF
      IG=IG+8  !CAMBIA GLI INDICI PER IL VERTICE SUCCESSIVO...8x4=32 RIGHE DI SIG
   20 CONTINUE
   10 CONTINUE
      WRITE(2,*)
      WRITE(2,*)'LA TENSIONE IDEALE MASSIMA SECONDO VON MISES VALE:',
     1SMAXMISES,'E SI HA NELL''ELEMENTO ',NMAXMISES
      WRITE(*,*)'LA TENSIONE IDEALE MASSIMA SECONDO VON MISES VALE:',
     1SMAXMISES,'E SI HA NELL''ELEMENTO ',NMAXMISES
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     SUBROUTINE CALC:
C
      SUBROUTINE CALC(NODES,NELEMS,NVERT,IFMAX,IFXY,FXY,ICNMAX,
     $ICNXY,CNXY,YM,PR,ITDP,XY,MBAND,FORCE)
C     DIMENSIONAMENTI VARIABILI DELLA STRUTTURA
      DIMENSION XY(2,NODES),NVERT(4,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
     $ICNXY(2,ICNMAX),CNXY(2,ICNMAX),STRUTK(2*NODES,MBAND),
     $FORCE(2*NODES)
C     DIMENSIONAMENTI FISSI
      DIMENSION D(3,3),B(3,8),AK(8,8),IPOINT(8),SIG(32,NELEMS)
C     INIZIALIZZAZIONE DELLE MATRICI
      PRINT*,' INIZIALIZZAZIONE DELLE MATRICI'
      CALL CLEAR(2*NODES,MBAND,STRUTK)
      CALL CLEAR(2*NODES,1,FORCE)
      CALL CLEAR(3,8,B)
!     SCRITTURA DELLA MATRICE D CONTENENTE IL LEGAME COSTITUTIVO
      WRITE(*,*)' ELEBORAZIONE DEL LEGAME COSTITUTIVO'
      CALL DMAT(YM,PR,ITDP,D)
!     DO PRINCIPALE DELLA SCRITTURA DELLA MATRICE AK E ASSEMBLAGGIO NELLA
!     STUTK
      WRITE(*,*)'ASSEMBLAGGIO DELLA MATRICE GLOBALE '
      WRITE(*,*)'  DI RIGIDEZZA DELLA STRUTTURA'
      DO 10,I10=1,NELEMS
      XI=XY(1,NVERT(1,I10))
      YI=XY(2,NVERT(1,I10))
      XJ=XY(1,NVERT(2,I10))
      YJ=XY(2,NVERT(2,I10))
      XK=XY(1,NVERT(3,I10))
      YK=XY(2,NVERT(3,I10))
      XL=XY(1,NVERT(4,I10))
      YL=XY(2,NVERT(4,I10))
      CALL KISOQ(YM,PR,XI,YI,XJ,YJ,XK,YK,XL,YL,D,B,AK)
      CALL POINTVT(NELEMS,I10,NVERT,IPOINT)
      CALL ASSEMBL(AK,NODES,IPOINT,MBAND,STRUTK)
   10 CONTINUE
!     CARICAMENTO DELLA STRUTTURA
      WRITE(*,*)' CARICAMENTO DELLA STRUTTURA'
      CALL FORCES(IFMAX,IFXY,FXY,NODES,FORCE)
!     VINCOLAMENTO DELLA STRUTTURA
      WRITE(*,*)' VINCOLAMENTO DEI NODI'
      CALL CNSTNG(ICNMAX,ICNXY,CNXY,NODES,FORCE,MBAND,STRUTK)
!     SOLUZIONE DEL PROBLEMA DEGLI SPOSTAMENTI NODALI
      WRITE(*,*)' RISOLUZIONE DEL SISTEMA DI EQUAZIONI AGLI SPOSTAMENTI'
      write(2,*) "========================================="
      do i=1,2*nodes
      write(2,*) (strutk(i,j),j=1,MBAND)
      enddo
            write(2,*) (force(i),i=1,2*nodes)
      write(2,*) "========================================="
      CALL COOK(2*NODES,MBAND,STRUTK,FORCE)
      WRITE(2,*)' '
      WRITE(2,*)'SPOSTAMENTI NODALI'
      WRITE(2,'(A6,2A20)')'NODO:','DELTA_X:','DELTA_Y:'
      DO 20,I20=1,NODES
      WRITE(2,'(I6,2G20.5)')I20,FORCE(I20*2-1),FORCE(I20*2)
   20 CONTINUE
C     POST PROCESSING
      WRITE(*,*)' CALCOLO DELLE TENSIONI'
      CALL TENSMAT(PR,ITDP,NODES,NELEMS,XY,NVERT,B,D,FORCE,SIG)
      RETURN
      END