Sources Esope | Cast3M

heflmi
C HEFLMI    SOURCE    CHAT      05/01/13    00:23:41     5004
 
      SUBROUTINE HESSFLMIEHE (SIG,VAR,AMAT,NDIMA)
      IMPLICIT INTEGER(I-N)
      integer ndima,m,i,j
      real*8  sig(3),var,amat(ndima,ndima)
      real*8  desv(3),aux1,aux2,alpha1,alpha2,alpha3,amataux(3,3)
      real*8  valor
      real*8            sigy0,kiso,siginf,velo,cpar,mpar
      common /miehdata/ sigy0,kiso,siginf,velo,cpar,mpar
      call zzero(amat,ndima*ndima)
      m=nint(mpar)
      call yieldmieheb(sig,valor)
      aux1=valor**(1-2*m)/2.D0*sqrt(2.D0/3.D0)
      aux2=aux1*cpar*3.D0**(2*m)/(2.D0**(2*m-1)+1.D0)*float(2*m-1)
      aux1=aux1*(1.D0-cpar)*float(2*m-1)
      call Desviador(sig,desv,3)
      if (m.gt.1) then
       alpha1=(desv(1)-desv(2))**(2*m-2)
       alpha2=(desv(2)-desv(3))**(2*m-2)
       alpha3=(desv(3)-desv(1))**(2*m-2)
      else
       alpha1=1.D0
       alpha2=1.D0
       alpha3=1.D0
      endif
      amat(1,1)=alpha1+alpha3
      amat(2,2)=alpha1+alpha2
      amat(3,3)=alpha2+alpha3
      amat(1,2)=-alpha1*aux1
      amat(1,3)=-alpha3*aux1
      amat(2,3)=-alpha2*aux1
      amat(2,1)=amat(1,2)
      amat(3,1)=amat(1,3)
      amat(3,2)=amat(2,3)
      if (m.gt.1) then
       do i=1,3
        amat(i,i)=aux1*amat(i,i)+aux2*desv(i)**(2*m-2)
       enddo
      else
       do i=1,3
        amat(i,i)=aux1*amat(i,i)+aux2
       enddo
      endif
c desviador de columnas
      do i=1,3
       call Desviador(amat(1,i),amataux(1,i),3)
      enddo
c trasponer amat
      aux1=amataux(1,2)
      amataux(1,2)=amataux(2,1)
      amataux(2,1)=aux1
      aux1=amataux(1,3)
      amataux(1,3)=amataux(3,1)
      amataux(3,1)=aux1
      aux1=amataux(2,3)
      amataux(2,3)=amataux(3,2)
      amataux(3,2)=aux1
c desviador de columnas
      do i=1,3
       call Desviador(amataux(1,i),amat(1,i),3)
      enddo
c no se traspone otra vez pq es simetrica
      call vflsigmiehe(sig,amataux(1,1))
      aux2=float(1-2*m)/valor*sqrt(3.D0/2.D0)
      do i=1,3
       do j=1,3
        amat(i,j)=amat(i,j)+amataux(i,1)*amataux(j,1)*aux2
       enddo
      enddo
      return
      end