c programma classico isoparametrico 4 nodi
c note negative
c - solo isoparametrico 4 nodi, non estendibile né adattabile
c - solo carichi concentrati, no carichi distribuiti 
c - solo vincolamenti in direzione x,y, no cambio di coordinate o servolink
c - utilizza massicciamente la dichiarazione implicita:
c   le variabili il cui nome inizia per [ijklmn] sono integer,
c   le altre real (singola precisione)
c
c note positive:
c - è già pronto.
c - funziona.
c

c234567========================================================= MAIN ==
      program iso4

c   parametri di dimensionamento:
c   - numero di nodi della struttura
      parameter (nodes=9)
c   - numero di elementi della struttura
      parameter (nelems=4)
c   - numero di nodi caricati e di nodi vincolati 
      parameter(ifmax=3,icnmax=3)

c   larghezza di banda
      parameter(mband=10)

C     dimensionamento delle matrici e dei vettori
      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)
C     dimensionamenti fissi
      dimension d(3,3),b(3,8),db(3,8),ak(8,8),ipoint(8),
     +deltael(8),sigmael(3)

c     punti di campionamento delle tensioni in coordinate xi (1a riga) 
c      ed eta (seconda riga)
      parameter(nsmppts=4)
      dimension  smppts(2,nsmppts)
c     definizione coordinate (per colonne)
      data smppts/ -0.57735,-0.57735, 
     $             -0.57735, 0.57735, 
     $              0.57735, 0.57735, 
     $              0.57735,-0.57735/
     
     
C     inizio del programma
c      open(1,file='mesh_iso4_comp.dat')
c      open(2,file='output_iso4_comp.4.txt')
c      open(1,file='mesh_iso4_flex.dat')
c      open(2,file='output_iso4_flex.4.txt')
      open(1,file='mesh_iso4_verifica.dat')
      open(2,file='output_iso4_verifica.4.txt')

c     leggo da file (unit=1) i dati del problema
      call readin(ym,pr,itdp,nodes,xy,nelems,nvert,ifmax,ifxy,fxy,icnmax
     +,icnxy,cnxy)

c     scrivo su file (unit=2) i dati letti, per verifica
      call echo(ym,pr,nodes,xy,nelems,nvert,ifmax,ifxy,fxy,icnmax,icnxy,
     +cnxy)

c     azzero la matrice di rigidezza e il vettore delle forze
      call clear(2*nodes,mband,strutk)
      call clear(2*nodes,1,force)

c     azzero la matrice b di trasformazione spostamenti -> deformazioni
c      per il singolo elemento (quello sotto esame in un dato punto del 
c      codice)
      call clear(3,8,b)

c     definisco la matrice di legame costitutivo, comune a tutti gli 
c      elementi
      call dmat(ym,pr,itdp,d)
      
c     ciclo di assemblaggio sugli elementi, ripeto per ogni elemento
      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))
      call kisoq(x1,y1,x2,y2,x3,y3,x4,y4,d,b,ym,pr,itdp,ak)
      call pointvt(nelems,i10,nvert,ipoint)
c TODO: togliere...
      write(2,*) (ipoint(i100),i100=1,8)
c
      call assembl(ak,nodes,ipoint,mband,strutk)
   10 continue

c a questo punto ho la matrice di rigidezza di struttura
c completamente assemblata; mancano solo i vincoli

c     costruisco il vettore dei termini noti, a partire dai carichi 
c      nodali
      call forces(ifmax,ifxy,fxy,nodes,force)
      
c     applico i vincoli alla struttura
      call cnstng(icnmax,icnxy,cnxy,nodes,force,mband,strutk)

c a questo punto ho una matrice di rigidezza della struttura
c assemblata e vincolata, ho un vettore dei termini noti
c assemblato (ho inserito i carichi esterni) e con termini
c legati allo spostamento imposto

      
c     per controllo, visualizzo in output la matrice di rigidezza della
c      struttura
c     print *,'Matrice di rigidezza'
c     do 11,i11=1,2*nodes
c     print *,(strutk(i11,i12),i12=1,mband)
c  11 continue

c     richiamo il solutore per matrici in forma bandata
      call cook(2*nodes,mband,strutk,force)
c     tale solutore lavora "in place", ovvero stocca dati temporanei
c      sulla matrice di sistema (dovrei farne una copia per non perderla
c      ) e scrive la soluzione sopra il vettore delle forze

c     ricopio tale soluzione sul vettore uv
      do 13,i13=1,2*nodes
      uv(i13)=force(i13)
   13 continue

c    scrivo in output gli spostamenti nodali
      write (2,*)'Spostamenti nodali'
      write (2,*)'Numero nodo, ux, uy'
      do 20,i20=1,nodes
      write(2,*) I20,uv(I20*2-1),uv(I20*2)
   20 continue

C     postprocessor per calcolare le tensioni negli elementi

c     intestazione output
      write (2,*)'Tensioni dell elemento'
      write (2,*)'Numero dell elemento, numero punto integrazione,'
      write (2,*)' Coord X, Coord Y, sx, sy, txy, Von Mises'

c     definisco punti di campionamento delle tensioni su elemento


c     ciclo sugli elementi
      do 30,I30=1,nelems

c     estraggo coordinate nodali
      x1=xy(1,nvert(1,i30))
      y1=xy(2,nvert(1,i30))
      x2=xy(1,nvert(2,i30))
      y2=xy(2,nvert(2,i30))
      x3=xy(1,nvert(3,i30))
      y3=xy(2,nvert(3,i30))
      x4=xy(1,nvert(4,i30))
      y4=xy(2,nvert(4,i30))

c     estraggo spostamenti nodali
      deltael(1)=uv(2*nvert(1,i30)-1)
      deltael(2)=uv(2*nvert(1,i30))
      deltael(3)=uv(2*nvert(2,i30)-1)
      deltael(4)=uv(2*nvert(2,i30))
      deltael(5)=uv(2*nvert(3,i30)-1)
      deltael(6)=uv(2*nvert(3,i30))
      deltael(7)=uv(2*nvert(4,i30)-1)
      deltael(8)=uv(2*nvert(4,i30))

c     procedo a calcolare coordinate e valori tensionali ai punti
c     di integrazione

      do 40,i40=1,nsmppts
      csi=smppts(1,i40)
      eta=smppts(2,i40)
      call bisoq(x1,y1,x2,y2,x3,y3,x4,y4,csi,eta,b,detj)
      call evaluv(x1,y1,x2,y2,x3,y3,x4,y4,csi,eta,x,y)

      call prodmat(3,3,8,d,b,db)
      call prodmat (3,8,1,db,deltael,sigmael)
      sigid1=sqrt(sigmael(1)**2+sigmael(2)**2-sigmael(1)*sigmael(2)+
     +3*sigmael(3)**2)
     
      write (2,*) I30,I40,x,y,sigmael (1),sigmael(2),sigmael(3),sigid1

   40 continue
   30 continue

      stop
      end

c234567================================================== SUB ASSEMBL ==
      subroutine assembl(ak,nodes,ipoint,mband,strutk)
C     Asseblaggio per la matrice bandata con larghezza di banda mband

      dimension ak(8,8),strutk(2*nodes,mband),ipoint(8)
      do 10,i10=1,8
      do 20,i20=1,8
      ind1=ipoint(i10)
      ind2=ipoint(i20)
C     Si valuta quanto l'elemento considerato e' distante dalla diagonale per determinare quale colonna occupa
      ind3=ipoint(i20)-ipoint(i10)+1
      if(ind3.gt.mband) then
      write(2,*)'aumentare mband!!!'
      stop
      endif
C     Questo if permette di andare ad assemblare l'elemento soltanto se esso si trova sopra la diagonale principale
      if(ind2.ge.ind1) then
      strutk(ind1,ind3)=strutk(ind1,ind3)+ak(i10,i20)
      endif
   20 continue
   10 continue
      return
      end

c234567===================================================== SUB COOK ==
c     SUBROUTINE SOLUTORE PER MATRICE BANDATA
       
      subroutine COOK(NSIZE,MBAND,S,r)
      real S(NSIZE,MBAND),r(NSIZE)

c     riduzione in forma triangolare superiore

      do ibase = 1,NSIZE-1                     
c                                              ciclo principale per azzerare 
c                                              i termini sottodiagonali della
c                                              colonna "ibase"
       do iscorri = 1,min(MBAND-1,NSIZE-ibase) 
c                                              ciclo secondario sulle righe 
c                                              ibase+1, ibase+2, ibase+3... 
c                                              fino ad esaurimento banda in
c                                              ibase+MBAND-1 o esaurimento 
c                                              matrice in ibase + (NSIZE-ibase)
       coeff = -S(ibase,1+iscorri)/S(ibase,1)
       do jscorri = 1,MBAND-iscorri           
c                                              ciclo terziario sugli elementi
c                                              della riga iscorri, da 1 a 
c                                              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 
c                                              è nella bandata scostato di iscorri colonne
        enddo
        r(ibase+iscorri) = r(ibase+iscorri)+coeff*r(ibase)   
c                                              non scordiamoci 
c                                              di aggiornare il termine noto.
       enddo
      enddo

c     backsubstitution

      do ibase = NSIZE,1,-1          
c                                     scorro a rovescio il vettore delle incognite, 
c                                     sovrascritto via via al termine noto r
       somma = 0.0
       do iscorri = 1,min(NSIZE-ibase,MBAND-1) !scorro sui termini gia' risolti
        somma = somma + r(ibase+iscorri)*S(ibase,1+iscorri)
       enddo
       r(ibase) = (r(ibase) - somma) / S(ibase,1)
      enddo
      end

c234567==================================================== SUB BISOQ ==
      subroutine bisoq(x1,y1,x2,y2,x3,y3,x4,y4,csi,eta,b,detj)
      dimension aj(2,2), ajinv(2,2),b(3,8)
c     calcolo matrice jacobiana trasposta
      call diffuv(x1,y1,x2,y2,x3,y3,x4,y4,csi,eta,aj(1,1),aj(1,2),
     +aj(2,1),aj(2,2))
c     calcolo matrice inversa della jacobiana trasposta
      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 con u1=1 e gli altri spostamenti 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     il termine dvdy viene omesso perché non dipende da u1 e quindi
C     in questo caso risulta nullo
      b(1,1)=dudx
C     il termine b(2,1) è nullo perché in questo caso deve essere uguale
C     a dvdy
      b(2,1)=0
      b(3,1)=dudy+dvdx
C     caso con v1=1 e tutti gli altri spostamenti nulli
      call diffuv(0.,1.,0.,0.,0.,0.,0.,0.,csi,eta,dudcsi,dvdcsi,dudeta,
     +dvdeta)
C     il termine dudx è nullo perché non dipende da v1
      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
      b(1,2)=0
C     dudx infatti è nullo
      b(2,2)=dvdy
      b(3,2)=dudy+dvdx
C     caso u2=1 e gli altri spostamenti 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
      b(1,3)=dudx
      b(2,3)=0
      b(3,3)=dudy+dvdx
C     caso con v2=1 e gli altri spostamenti nulli
      call diffuv(0.,0.,0.,1.,0.,0.,0.,0.,csi,eta,dudcsi,dvdcsi,dudeta,
     +dvdeta)
      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
      b(1,4)=0
      b(2,4)=dvdy
      b(3,4)=dudy+dvdx
C     caso con u3=1
      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
      b(1,5)=dudx
      b(2,5)=0
      b(3,5)=dudy+dvdx
C     caso con v3=1 e gli altri spostamenti nulli
      call diffuv(0.,0.,0.,0.,0.,1.,0.,0.,csi,eta,dudcsi,dvdcsi,dudeta,
     +dvdeta)
      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
      b(1,6)=0
      b(2,6)=dvdy
      b(3,6)=dudy+dvdx
C     caso con u4=1
      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
      b(1,7)=dudx
      b(2,7)=0
      b(3,7)=dudy+dvdx
C     caso con v4=1 e gli altri spostamenti nulli
      call diffuv(0.,0.,0.,0.,0.,0.,0.,1.,csi,eta,dudcsi,dvdcsi,dudeta,
     +dvdeta)
      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
      b(1,8)=0
      b(2,8)=dvdy
      b(3,8)=dudy+dvdx
      return
      end
c234567=================================================== SUB EVALUV ==
      subroutine evaluv(u1,v1,u2,v2,u3,v3,u4,v4,csi,eta,u,v)
C     campionamento di u,v in funzione dei valori nodali, eta e csi
      u=0.25*(1-csi)*(1-eta)*u1+0.25*(1+csi)*(1-eta)*
     +u2+0.25*(1+csi)*(1+eta)*u3+0.25*(1-csi)*(1+eta)*u4
      v=0.25*(1-csi)*(1-eta)*v1+0.25*(1+csi)*(1-eta)*
     +v2+0.25*(1+csi)*(1+eta)*v3+0.25*(1-csi)*(1+eta)*v4
      return
      end
c234567=================================================== SUB DIFFUV ==
      subroutine diffuv(u1,v1,u2,v2,u3,v3,u4,v4,csi,eta,dudcsi,dvdcsi,
     +dudeta,dvdeta)
C     calcolo delle derivate rispetto a eta e csi
      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
c234567==================================================== SUB KISOQ ==
      subroutine kisoq(x1,y1,x2,y2,x3,y3,x4,y4,d,b,ym,pr,itdp,ak)
C     determinazione della matrice di rigidezza dell'elemento
      dimension d(3,3),b(3,8),bt(8,3),db(3,8),aktemp(8,8),ak(8,8)

c     numero e definizione punti di gauss
      parameter(ngp=2)
      dimension gpc(ngp) , gpw(ngp)
c     assegno valori a posizione e pesi
      data gpc/-0.57735,0.57735/, gpw/1.,1./
      
c     matrice di legame costitutivo
      call dmat(ym,pr,itdp,d)

c     azzero la matrice di rigidezza di elemento
      call clear(8,8,ak)

C     integrazione gaussiana a 4 punti
c     scorro sui punti di gauss lungo l'asse xi, a posiz. xi=gpc(i10)
      do 10,i10=1,ngp
c      scorro sui punti di gauss lungo l'asse eta, a posiz eta=gpc(i20)
       do 20,i20=1,ngp
C     calcolo dei contributi relativi ai vari punti di integrazione
        call bisoq(x1,y1,x2,y2,x3,y3,x4,y4,gpc(i10),gpc(i20),b,detj)
        call btmat(b,bt)
        call prodmat(3,3,8,d,b,db)
        call prodmat(8,3,8,bt,db,aktemp)
C     accumulo dei contribuiti relativi ai vari punti di integrazione
        call accmat(8,8,gpw(i10)*gpw(i20)*detj,aktemp,ak)
   20  continue
   10 continue
   
      return
      end

c234567===================================================== SUB DMAT ==
      subroutine dmat(ym,pr,itdp,d)
C     d= matrice [3x3], sigma= d*epsilon
      dimension d(3,3)
      if(itdp.EQ.0) goto 100
      if (itdp.eq.1) goto 200
100   continue
      ctp=ym/(1-pr**2)
      g=ym/((1+pr)*2)
      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
200   continue
      cdp =ym*(1-pr)/((1+pr)*(1-2*pr))
      cpr=pr/(1-pr)
      g=ym/((1+pr)*2)
      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
      
c234567==================================================== SUB CLEAR ==
      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


c234567==================================================== SUB BTMAT ==
      subroutine btmat(b,bt)
C Trasposta di matrice A[k1,k2]=AT[k2,k1]
      dimension b(3,8),bt(8,3)
      do 10,I10=1,3
      do 20,I20=1,8
      bt(I20,I10)=b(I10,I20)
   20 continue
   10 continue
      return
      end

c234567================================================== SUB PRODMAT ==
      subroutine prodmat(n1,n2,n3,a,b,ab)
C Prodotto matrici A[k1,k2]*B[k2,k3]=AB[k1,k3]
      dimension a(n1,n2),b(n2,n3),ab(n1,n3)
      do 10,i10=1,n1
      do 30,i30=1,n3
      ab(i10,i30)=0.
      do 20,i20=1,n2
      ab(i10,i30)=ab(i10,i30)+a(i10,i20)*b(i20,i30)
   20 continue
   30 continue
   10 continue
      return
      end

c234567================================================== SUB AGGIORN ==
      subroutine accmat(nrow,ncol,scalare,contributo,accumulo)
      dimension contributo(nrow,ncol),accumulo(nrow,ncol)
      do 10,I10=1,nrow
      do 20,I20=1,ncol
      accumulo(I10,I20)=accumulo(i10,I20)+contributo(i10,I20)*scalare
   20 continue
   10 continue
      return
      end
      

c234567===================================================== SUB ECHO ==
      subroutine echo(ym,pr,nodes,xy,nelems,nvert,ifmax,ifxy,fxy,icnmax,
     +icnxy,cnxy)
C     controlla se il programma legge correttamente gli input pag 114
      dimension xy(2,nodes),nvert(4,nelems),ifxy(ifmax),fxy(2,ifmax),
     +icnxy(2,icnmax),cnxy(2,icnmax)
      write(*,*)'modulo di Young=',ym
      write(*,*)'coefficiente di Poisson=',pr
      write(*,*)'nodo:     cordx:       cordy:'
      do 10,I10=1,nodes
      write(*,*) I10,xy(1,I10),xy(2,I10)
   10 continue
      write(*,*)'elemento: v1:    v2:    v3:   v4:'
      do 20,I20=1,nelems
      write(*,*) I20,nvert(1,I20),nvert(2,I20),nvert(3,I20),nvert(4,I20)
   20 continue
      write(*,*)'nodo:  fx:    fy:'
      do 30,I30=1,ifmax
      write(*,*) ifxy(I30),fxy(1,I30),fxy(2,I30)
   30 continue
      write(*,*)'nodo: vinc:   spostx:    sposty:'
      do 40,I40=1,icnmax
      write(*,*) icnxy(1,I40),icnxy(2,I40),cnxy(1,I40),cnxy(2,I40)
   40 continue
      return
      end
      
c234567=================================================== SUB IPOINT ==
      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

c234567=================================================== SUB FORCES ==
      subroutine forces(ifmax,ifxy,fxy,nodes,force)
C     Subroutine di riempimento delle forze nodali
C     ifmax=numero massimo di nodi caricati
C     ifxy=vettore con il nome dei nodi caricati
C     fxy=matrice con le componenti delle forze nodali lungo x e y
C     force=vettore termini noti
      dimension ifxy(ifmax),fxy(2,ifmax),force(2*nodes)
      do 10,i10=1,ifmax
C     creo gli indici globali
      indx=ifxy(i10)*2-1
      indy=ifxy(i10)*2
C     creo il vettore dei termini noti
      force(indx)=fxy(1,i10)
      force(indy)=fxy(2,i10)
   10 continue
      return
      end

c234567=================================================== SUB CNSTNG ==
      subroutine cnstng(icnmax,icnxy,cnxy,nodes,force,mband,strutk)
C     Subroutine di vincolamento
C     icnmax=numero massimo di nodi vincolati
C     incxy=matrice con il nome dei nodi vincolati ed il tipo 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 spostamenti assegnati lungo x e y
C     force=vettore termini noti
C     strutk=matrice di rigidezza della struttura
      dimension icnxy(2,icnmax),cnxy(2,icnmax)
      dimension strutk(2*nodes,mband),force(2*nodes)

      do 10,i10=1,icnmax
      if     (icnxy(2,i10).eq.1) then
       call vincola(2*nodes,mband,strutk,force,
     $                                    icnxy(1,i10)*2-1, cnxy(1,i10))
      else if(icnxy(2,i10).eq.2) then
       call vincola(2*nodes,mband,strutk,force,
     $                                    icnxy(1,i10)*2  , cnxy(2,i10))
      else if(icnxy(2,i10).eq.3) then
       call vincola(2*nodes,mband,strutk,force,
     $                                    icnxy(1,i10)*2-1, cnxy(1,i10))
       call vincola(2*nodes,mband,strutk,force,
     $                                    icnxy(1,i10)*2  , cnxy(2,i10))
      else
       write(0,*)'errore su vincolamento'
       stop
      end if
   10 continue
   
      return
      end

c234567================================================== SUB VINCOLA ==
c      SUBROUTINE VINCOLAMENTO PER MATRICE BANDATA
       subroutine VINCOLA(NSIZE,MBAND,Aband,b,idof,val)
       real Aband(NSIZE,MBAND),b(NSIZE)

c      correzione termine diagonale A(idof,idof)==Aband(idof,1) e b(idof)
       Aband(idof,1)=1.0
       b(idof)=val

c      spostamento termini sopradiagonali di A (*val) a termine noto
c      e annullamento degli stessi
       do isd=1,min(MBAND-1,idof-1) !evito di sforare la banda e il limite superiore della matrice
         b(idof-isd) = b(idof-isd)-Aband(idof-isd,1+isd)*val
         Aband(idof-isd,1+isd)=0.0
       enddo

c      spostamento termini sottodiagonali di A (*val) a termine noto
c      e annullamento degli stessi
       do isd=1,min(MBAND-1,NSIZE-idof) !evito di sforare la banda e il limite inferiore della matrice
         b(idof+isd) = b(idof+isd)-Aband(idof,1+isd)*val
         Aband(idof,1+isd)=0.0
       enddo

       end subroutine

c234567=================================================== SUB READIN ==
      subroutine readin (ym,pr,itdp,nodes,xy,nelems,nvert,ifmax,ifxy,
     +fxy,icnmax,icnxy,cnxy)
C     itdp=indice di scelta tra tensione piana (itdp=0) e deformazione
C     piana (itdp=1)
      dimension xy(2,nodes), nvert(4,nelems),ifxy(ifmax),fxy(2,ifmax),
     +icnxy(2,icnmax),cnxy(2,icnmax)
C     legge da file le caratteristiche meccaniche, le coordinate nodali
C     e se l elemento è in tensione o deformazione piana
      read(1,*) ym,pr,itdp
      do 10,I10=1,nodes
C     legge le coordinate dei nodi
      read(1,*)node,xy(1,I10),xy(2,I10)
      if(I10.NE.node)then
      write(0,*) 'ERR: nodo', node , 'dovrebbe essere', I10
      stop
      end if
   10 continue
      write(2,*) 'nome 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(0,*) 'ERR: elemento', nelem, 'dovrebbe essere', I20
      stop
      end if
   20 continue
      do 30,I30=1,ifmax
      read(1,*)ifxy(I30),fxy(1,I30),fxy(2,I30)
      write(2,*)ifxy(I30),fxy(1,I30),fxy(2,I30)
   30 continue
      do 40, I40=1,icnmax
      read(1,*)icnxy(1,I40),icnxy(2,I40),cnxy(1,I40),cnxy(2,I40)
   40 continue
      return
      end