C D2ENDP    SOURCE    CB215821  16/04/21    21:16:13     8920

      SUBROUTINE DER2_ENERELAS_DPRAL(X,E,NDIMX,NMODEL)
      IMPLICIT INTEGER(I-N)
      integer nmodel,i,ndimx,j
      real*8 x(3),e(ndimx,ndimx),aux0,aux1,aux2,aux3,aux(3,3)
      real*8                  you,xnu
      common /elasdata/       you,xnu
      real*8                  kap,xmu,xmu1,alp1,xmu2,alp2,xmu3,alp3
      logical                 ogden
      common /elasdata_ogden/ kap,xmu,xmu1,alp1,xmu2,alp2,xmu3,alp3,
     .                        ogden

      call zzero(e,ndimx*ndimx)
      if (ogden) goto 10
c Modelo cuadratico
      aux0=you/((1.D0+xnu)*(1.D0-2.D0*xnu))
      aux1=(1.D0-xnu)*aux0
      aux2=xnu*aux0
      aux3=you/(2.D0*(1.D0+xnu))
      do i=1,3
        E(i,i)=aux1
      enddo
      E(1,2)=aux2
      E(1,3)=aux2
      E(2,1)=aux2
      E(2,3)=aux2
      E(3,1)=aux2
      E(3,2)=aux2
      if (ndimx.eq.4) E(4,4)=1.D0
c hacer desviador para VMT
      if (nmodel.eq.8) goto 20
      return
c Modelo ogden
10    continue
      do i=1,3
       aux(i,1)=EXP(x(i))
      enddo
      e(1,1)=xmu*(xmu1*alp1*aux(1,1)**alp1+
     .            xmu2*alp2*aux(1,1)**alp2+
     .            xmu3*alp3*aux(1,1)**alp3)
      e(2,2)=xmu*(xmu1*alp1*aux(2,1)**alp1+
     .            xmu2*alp2*aux(2,1)**alp2+
     .            xmu3*alp3*aux(2,1)**alp3)
      e(3,3)=xmu*(xmu1*alp1*aux(3,1)**alp1+
     .            xmu2*alp2*aux(3,1)**alp2+
     .            xmu3*alp3*aux(3,1)**alp3)
c     desviador de columnas
20    continue
      do i=1,3
       call Desviador(e(1,i),aux(1,i),3)
      enddo
c     trasponer amat
      aux1    =aux(1,2)
      aux(1,2)=aux(2,1)
      aux(2,1)=aux1
      aux1    =aux(1,3)
      aux(1,3)=aux(3,1)
      aux(3,1)=aux1
      aux1    =aux(2,3)
      aux(2,3)=aux(3,2)
      aux(3,2)=aux1
c     desviador de columnas
      do i=1,3
       call Desviador(aux(1,i),e(1,i),3)
      enddo
      goto 100
c Calculadas con derivacion, dan igual!
      call zzero(e,ndimx*ndimx)
      call der_enerelas_dpral(x,aux(1,1),nmodel)
      do j=1,3
        do i=1,3
          aux(i,2)=x(i)
        enddo
        aux(j,3)=max(x(j)*1.D-6,1.D-8)
        aux(j,2)=x(j)+aux(j,3)
        call der_enerelas_dpral(aux(1,2),e(1,j),nmodel)
        do i=1,3
         e(i,j)=(e(i,j)-aux(i,1))/aux(j,3)
        enddo
      enddo
c
100   continue
      if (ndimx.eq.4) E(4,4)=1.D0
      return
      end