Sources Esope | Cast3M

fatig3
C FATIG3    SOURCE    KICH      22/09/22    21:15:03     11465          
      SUBROUTINE FATIG3(sig,p,raytau,nbrobl,nombre,icri,
     &alpha,beta,ycrit,s,ib,igau)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 centre(5),stab(3,3),d(3),e(3),v(3),xnorm(3,3),u(5)
      REAL*8 sig(nbrobl,nombre),p(nombre),raytau(nombre)
      REAL*8 s(nbrobl-1,nombre)
 
      SQ2 = DSQRT(2.d0)
      SQ6 = DSQRT(6.d0)
C
c     *********************************************************
c     appel des critere
c     *********************************************************
c
      if (icri.eq.3) then
C PAPADOPOULOS
c
*      call dvpoly(s,centre,elmin,nbrobl,nombre,ib,igau)
      call dvboul(s,centre,ray,nbrobl,nombre,ib,igau)
c
*      pmax = -1.0d10
*      pmin =  1.0d10
       pmax = p(1)
c
      ycri = -1.D10
      rayon = 0.D0
      do i=1,nombre
        stab(1,2) = 0.d0
        stab(2,1) = stab(1,2)
        stab(1,3) = 0.d0
        stab(3,1) = stab(1,3)
        if (p(i).gt.pmax) pmax = p(i)
*        if (p(i).lt.pmin) pmin = p(i)
        stab(1,1) = (s(1,i) - centre(1))/SQ2
     .             +(s(2,i) - centre(2))/SQ6
        stab(2,2) =-(s(1,i) - centre(1))/SQ2
     .             +(s(2,i) - centre(2))/SQ6
        stab(3,3) = - stab(1,1) - stab(2,2)
       if(nbrobl.gt.4) then
        stab(1,2) = s(5,i) - centre(5)
        stab(2,1) = s(5,i) - centre(5)
        stab(1,3) = s(4,i) - centre(4)
        stab(3,1) = s(4,i) - centre(4)
       endif
        stab(2,3) = s(3,i) - centre(3)
        stab(3,2) = s(3,i) - centre(3)
        raytau(i)  = dsqrt(stab(1,1)**2+stab(2,2)**2+stab(3,3)**2+
     .   2.*(stab(1,2)**2+stab(2,3)**2+stab(1,3)**2))/SQ2
 
       ycri0 = (raytau(i) + alpha*p(i) - beta)/(beta-alpha*p(i))
*      write(6,*) 'PP',i,p(i),raytau(i),ycri
       if (raytau(i).gt.rayon) rayon=raytau(i)
       if (ycri0.gt.ycri) ycri=ycri0
      enddo
c      ycrit0 = (ray + alpha*pmax - beta)/(beta-alpha*pmax)
      ycrit = (rayon + alpha*pmax - beta)/(beta-alpha*pmax)
*      write(6,*) 'PPF',pmax,ray,rayon,ycrit,ycri
c
      return
      endif
c
      if(icri.eq.2) then
C DANG VAN
c
      critere = -1.d10
c
      call dvboul(s,centre,ray,nbrobl,nombre,ib,igau)
*      if (ib.eq.1.and.igau.eq.1) then
*      write(6,*) 'ct',((centre(jj)/1.e6),jj= 1,5)
*      write(6,*) (ray/1.e6),nbrobl,nombre
*      endif
c
      do i=1,nombre
        stab(1,2) = 0.d0
        stab(2,1) = stab(1,2)
        stab(1,3) = 0.d0
        stab(3,1) = stab(1,3)
        stab(1,1) = (s(1,i) - centre(1))/SQ2
     .             +(s(2,i) - centre(2))/SQ6
        stab(2,2) =-(s(1,i) - centre(1))/SQ2
     .             +(s(2,i) - centre(2))/SQ6
        stab(3,3) = - stab(1,1) - stab(2,2)
       if(nbrobl.gt.4) then
        stab(1,2) = s(5,i) - centre(5)
        stab(2,1) = s(5,i) - centre(5)
        stab(1,3) = s(4,i) - centre(4)
        stab(3,1) = s(4,i) - centre(4)
       endif
        stab(2,3) = s(3,i) - centre(3)
        stab(3,2) = s(3,i) - centre(3)
        call dvdiag(stab,d)
c
c       Classement des valeurs propres
c
        d3 = dabs(d(1)-d(2))
        d2 = dabs(d(1)-d(3))
        d1 = dabs(d(3)-d(2))
        if (d1.lt.1.d-04) d(2) = d(3)
        if (d2.lt.1.d-04) d(1) = d(3)
        if (d3.lt.1.d-04) d(2) = d(1)
c
        n1 = 1
        n2 = 2
        n3 = 3
        if (d(2).gt.d(1).and.d(3).gt.d(2)) then
           n1 = 3
           n2 = 2
           n3 = 1
        endif
        if (d(2).gt.d(3).and.d(3).gt.d(1)) then
           n1 = 2
           n2 = 3
           n3 = 1
        endif
        if (d(3).gt.d(1).and.d(1).gt.d(2)) then
           n1 = 3
           n2 = 1
           n3 = 2
        endif
        if (d(1).gt.d(3).and.d(3).gt.d(2)) then
           n1 = 1
           n2 = 3
           n3 = 2
        endif
        if (d(2).gt.d(1).and.d(1).gt.d(3)) then
           n1 = 2
           n2 = 1
           n3 = 3
        endif
c
        raytau(i) = dabs(d(n1)-d(n3))/2.
        crit = raytau(i) + alpha*p(i) - beta
        crit = crit / beta
c
c
        if (crit.gt.critere) then
           critere    = crit
*           critere    = crit / (beta-alpha*p(i))
           v(1)       = d(n1)
           v(2)       = d(n2)
           v(3)       = d(n3)
           xnorm(1,1) = stab(1,n1)
           xnorm(2,1) = stab(2,n1)
           xnorm(3,1) = stab(3,n1)
           xnorm(1,2) = stab(1,n2)
           xnorm(2,2) = stab(2,n2)
           xnorm(3,2) = stab(3,n2)
           xnorm(1,3) = stab(1,n3)
           xnorm(2,3) = stab(2,n3)
           xnorm(3,3) = stab(3,n3)
        endif
      enddo
c
      ycrit = critere
c
      return
      endif
 
      if(icri.eq.4) then
C SINES
c
      dmax = 0.d0
      do i=1,nombre-1
        do j=i+1,nombre
          dist = (s(1,i)-s(1,j))**2+(s(2,i)-s(2,j))**2
     .      +2.d0*((s(3,i)-s(3,j))**2)
          if (nbrobl.gt.4) then
           dist = dist + 2.d0*((s(4,i)-s(4,j))**2
     .          +(s(5,i)-s(5,j))**2)
          endif
          dist = dsqrt(dist)
          if (dist.gt.dmax) then
           dmax=dist
          endif
        enddo
      enddo
      dmax = dmax/2.d0/SQ2
c
      pmoyenne = 0.d0
      do i=1,nombre
        pmoyenne = pmoyenne+p(i)
      enddo
      pmoyenne = pmoyenne/float(nombre)
      raytau(1) = pmoyenne
      raytau(2) = dmax
      ycrit = (dmax+alpha*pmoyenne-beta)/(beta-alpha*pmoyenne)
      return
      endif
 
      if(icri.eq.5) then
C CROSS
c
      dmax = 0.d0
      do i=1,nombre-1
        do j=i+1,nombre
          dist = (s(1,i)-s(1,j))**2+(s(2,i)-s(2,j))**2
     .      +2.d0*((s(3,i)-s(3,j))**2)
          if (nbrobl.gt.4) then
           dist = dist + 2.d0*((s(4,i)-s(4,j))**2
     .          +(s(5,i)-s(5,j))**2)
          endif
          dist = dsqrt(dist)
          if (dist.gt.dmax) then
            dmax=dist
          endif
        enddo
      enddo
      dmax = dmax/2.d0/SQ2
c
*      pmax = -1.0d10
      pmax = p(1)
      do i=1,nombre
        if (p(i).gt.pmax) pmax=p(i)
      enddo
      raytau(1) = pmax
      raytau(2) = dmax
      ycrit = (dmax + alpha*pmax - beta)/(beta - alpha*pmax)
      return
      endif
 
      if(icri.eq.6) then
C DC
c
      dmax  = 0.d0
      dmaxp = 0.d0
c
c     ***********************
c     plus grand diametre ...
c     ***********************
      do i=1,nombre-1
        do j=i+1,nombre
          dist = (s(1,i)-s(1,j))**2+(s(2,i)-s(2,j))**2
     .      +2.d0*((s(3,i)-s(3,j))**2)
          if (nbrobl.gt.4) then
           dist = dist + 2.d0*((s(4,i)-s(4,j))**2
     .          +(s(5,i)-s(5,j))**2)
          endif
          dist = dsqrt(dist)
          if (dist.gt.dmax) then
            dmax=dist
            ic = i
            jc = j
          endif
        enddo
      enddo
      if(dmax.eq.0.) then
        dmaxp = 0.
      else
c       ***********************
c       deuxieme diametre ...
c       ***********************
        do k=1,nbrobl-1
          u(k) = (s(k,ic)-s(k,jc))/dmax
        enddo
        do k=1,nombre
         prod = u(1)*s(1,k)+u(2)*s(2,k)
     .     +2.d0*u(3)*s(3,k)
          if (nbrobl.gt.4) then
           prod = prod + 2.d0*(u(4)*s(4,k)+u(5)*s(5,k))
          endif         
         do kj=1,nbrobl-1
           s(kj,k) = s(kj,k)-prod*u(kj)
         enddo
        enddo
        do i=1,nombre-1
          do j=i+1,nombre
            dist = (s(1,i)-s(1,j))**2+(s(2,i)-s(2,j))**2
     .        +2.*((s(3,i)-s(3,j))**2)
          if (nbrobl.gt.4) then
           dist = dist + 2.d0*((s(4,i)-s(4,j))**2
     .          +(s(5,i)-s(5,j))**2)
          endif         
            dist = dsqrt(dist)
            if (dist.gt.dmaxp) dmaxp = dist
          enddo
        enddo
      endif
c
*      pmax = -1.0d10
      pmax = p(1)
      do i=1,nombre
        if (p(i).gt.pmax) pmax=p(i)
      enddo
c      write(6,*) 'dc',pmax,dmax,dmaxp
      a = dsqrt((dmax/2.)**2+(dmaxp/2.)**2)/SQ2
      raytau(1) = pmax
      raytau(2) = a
      ycrit = (a + alpha*pmax - beta)/(beta - alpha*pmax)
      return
      endif
 
      if(ICRI.eq.7) then
C VON MISES
c
      vonmax = 0.d0
      do i=1,nombre
          if (nbrobl.gt.4) then
        von = dsqrt(s(1,i)**2+s(2,i)**2
     .        +2.*(s(3,i)**2+s(4,i)**2+s(5,i)**2))/SQ2
          else
        von = dsqrt(s(1,i)**2+s(2,i)**2
     .        +2.*(s(3,i)**2))/SQ2
          endif
        if (von.gt.vonmax) then
         vonmax = von
        endif
      enddo
c
      ycrit = vonmax
      return
      endif
 
c
 900  format(9x,'1 = dvk',/9x,'2 = papado',/9X,'3 = cross',
     .   /9X,'4 = sines',/9X,'5 = dc',/9X,'6 = vm'/)
c
c
c
      end