C UMDT3D    SOURCE    FD218221  24/02/07    21:15:30     11834          
      subroutine umdt3d(E0,nu,souplesse66,dt3,umdt66,
     # a,b,x,ipzero,ngf,errg,iso)
      
c     calcul du tensseur d endommagement de traction directe
      
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      
      real*8 souplesse66(6,6)
      real*8 souplesse66d(6,6),raideur66d(6,6)
      real*8 dt3(3)
      integer ngf
      real*8 X(ngf),B(ngf),A(ngf,(ngf+1))
      integer errg,ipzero(ngf)
      real*8 umdt66(6,6)
      logical iso      
      real*8 dx
      real*8 sn3(3),d66(6,6)
      real*8 nu,E0
      

      if (iso) then
c           cas de la loi elastique isotrope
c           determination via solution type cf.
 
C           A. Sellier, G. Casaux-Ginestet, L. Buffo-Lacarrière, X. Bourbon,
C           Orthotropic damage coupled with localized crack reclosure processing. Part I: Constitutive laws,
C           Engineering Fracture Mechanics,
C           Volume 97,
C           2013,
C           Pages 148-167,
C           ISSN 0013-7944,
C           https://doi.org/10.1016/j.engfracmech.2012.10.012.

            do i=1,3
c               print*,dt3(i)
                 sn3(i)=1.d0/(1.d0-dt3(i))
             end do
c            calcul de la matrice d endommagement sans la borner et en
c            attenuant les effets de Poisson
             call b3_d66(nu,sn3,d66,.false.,.false.)
             do i=1,6
                 do j=1,6
                     if(i.eq.j) then
                         umdt66(i,j)=1.d0-d66(i,j)
                     else
                         umdt66(i,j)=-d66(i,j)
                     end if
                 end do
            end do
c           call afic66(umdt66)
        else
            print*,'verifier calcul de 1-Dt dans umdt3'
c           cas general : 1-d=(sd^-1)(s)       
c           print*,'ngf dans umdt3d',ngf
c           matrice de souplesse endommagee
c           seuls les termes de la diagonale sont affectes
            do i=1,6
             do j=1,6
              if(i.eq.j) then
               if(i.le.3) then
                  souplesse66d(i,j)=souplesse66(i,j)/(1.d0-dt3(i))
               else
                  call indce0(i,k,l)
                  dx=max(dt3(k),dt3(l))
                  souplesse66d(i,j)=souplesse66(i,j)/(1.d0-dx)
               end if
              else
               souplesse66d(i,j)=souplesse66(i,j)
              end if
             end do
            end do
c           matrice de la raideur endommagee 
c           obtenue par inversion de souplesse66d
            do i=1,6
c               chaque colonne est solution de ax=b        
                do j=1,6
                  do k=1,6 
                    a(j,k)=souplesse66d(j,k)  
                  end do
                  if(j.eq.i) then
                      b(j)=1.d0 
                  else
                      b(j)=0.d0  
                  end if
                end do
                call gaus3d(6,A,X,B,ngf,errg,ipzero)
                do j=1,6
                    raideur66d(j,i)=x(j)
                end do
           end do  
c          construction du tenseur umdt66
           do i=1,6
            do j=1,6
             umdt66(i,j)=0.d0
             do k=1,6            
               umdt66(i,j)=umdt66(i,j)+raideur66d(i,k)*souplesse66(k,j)
             end do
            end do
           end do
      end if
c      print*,'umdt dans umdt66'
c      call afic66(umdt66)
c     read*
      return
      end