Télécharger la plaquette Cast3M


fr en

Télécharger mtcldd.eso

Retour à la liste

Numérotation des lignes : 

mtcldd

  1. C MTCLDD    SOURCE    PV        22/04/19    16:18:07     11344          
  2. CCC
  3. C **********************************************************************
  4. CCC
  5.       SUBROUTINE MTC_LS_DPRAL_DENSI(EKMAT,NDIMK,X,NDIMX,LAM,
  6.      .                              xdensi,xmat,
  7.      .                              nmodel,nescri,ues,nnumer,deltax)
  8.  
  9.       IMPLICIT INTEGER(I-N)
  10.       IMPLICIT REAL*8 (a-h,o-z)
  11.       integer ndimk,ndimx,nmodel,nescri,ues,nnumer
  12.       integer i,j,k,ndims,ndimv
  13.       real*8  EKmat(ndimk,ndimk),x(ndimx),lam,deltax
  14.       real*8  amat(16),aamat(16),eemat(16),vecm(4),vecnaux(4),vecn(4),
  15.      .        gmat(16),Avecm(4),vecntA(4),kmat(4,4),g,void(1),sig(3)
  16.       real*8  xdensi,dxdensi,xdensi2,xmat(*),resu,resuaux
  17.       real*8  Mmat(4,3),kmat2(4,4),kmat3(4,4),ekmat2(3,3),ekmat3(3,3),
  18.      .        vecmaux(4)
  19.  
  20.       call zzero(amat ,16)
  21.       call zzero(aamat,16)
  22.       call zzero(eemat,16)
  23.       call zzero(gmat ,16)
  24.       call zzero(kmat ,16)
  25.       call zzero(vecm,4)
  26.       call zzero(vecn,4)
  27.       call zzero(vecnaux,4)
  28.       call zzero(vecn,4)
  29.       call zzero(Avecm,4)
  30.       call zzero(vecntA,4)
  31.       call zzero(sig,3)
  32.  
  33.       call zzero(mmat,12)
  34.       call zzero(kmat2,16)
  35.       call zzero(kmat3,16)
  36.       call zzero(ekmat2,9)
  37.       call zzero(ekmat3,9)
  38.       call zzero(vecmaux,4)
  39.  
  40.       void(1)=0.D0
  41.       ndimv=1
  42.       ndims=3
  43.       call der_enerelas_dpral(x,sig,nmodel)
  44. c m,n,A
  45.       if (ndimx.eq.3) then
  46.        call yieldd(sig,ndims,void,ndimv,resu,nmodel)
  47.        call vflsig(sig,ndims,void,ndimv,vecm,nmodel)
  48.        call vyisig(sig,ndims,void,ndimv,vecnaux,nmodel)
  49.        call HessFlsig(sig,ndims,void,ndimv,amat,ndimx,nmodel)
  50.       else if (ndimx.eq.4) then
  51.        call yieldd(sig,ndims,x(ndimx),ndimv,resu,nmodel)
  52.        call vflsig(sig,ndims,x(ndimx),ndimv,vecm,nmodel)
  53.        call vflvar(sig,ndims,x(ndimx),ndimv,vecm(ndimx),nmodel)
  54.        call vyisig(sig,ndims,x(ndimx),ndimv,vecnaux,nmodel)
  55.        call vyivar(sig,ndims,x(ndimx),ndimv,vecnaux(ndimx),nmodel)
  56.        call HessFlsig(sig,ndims,x(ndimx),ndimv,amat,ndimx,nmodel)
  57.       endif
  58. c E=d2_ener (ampliada de 3 a ndimx con 1 en la diagonal)
  59.       call der2_enerelas_dpral(x,eemat,ndimx,nmodel)
  60. c n^T=n^T*E
  61. c AA=A*E
  62.       do i=1,ndimx
  63.        vecn(i)=0.D0
  64.        do j=1,ndimx
  65.          vecn(i)=vecn(i)+vecnaux(j)*eemat(j+(i-1)*ndimx)
  66.          aamat(i+(j-1)*ndimx)=0.D0
  67.          do k=1,ndimx
  68.            aamat(i+(j-1)*ndimx)=aamat(i+(j-1)*ndimx)+
  69.      .                         amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx)
  70.          enddo
  71.        enddo
  72.       enddo
  73.       call zzero(amat,16 )
  74. c A=I+l*AA
  75.       do i=1,ndimx
  76.        amat(i+(i-1)*ndimx)=1.D0
  77.        do j=1,ndimx
  78.          amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  79.      .                       lam*aamat(i+(j-1)*ndimx)
  80.        enddo
  81.       enddo
  82. c G = A-1
  83.       call DescLU(Amat,ndimx)
  84.       call LUinv(Amat,Gmat,ndimx)
  85. c g = nt*A-1*m
  86.       g=0.D0
  87.       do i=1,ndimx
  88.        do j=1,ndimx
  89.           g=g+vecn(i)*Gmat(i+ndimx*(j-1))*vecm(j)
  90.        enddo
  91.       enddo
  92. c G*m; nt*G
  93.       do i=1,ndimx
  94.        do j=1,ndimx
  95.          Avecm(i)=Avecm(i)+Gmat(i+ndimx*(j-1))*vecm(j)
  96.          vecntA(i)=vecntA(i)+vecn(j)*Gmat(j+ndimx*(i-1))
  97.        enddo
  98.       enddo
  99. c K = G-(G*m * nt*G)/g
  100.       do i=1,ndimx
  101.        do j=1,ndimx
  102.           Kmat(i,j)=Gmat((j-1)*ndimx+i)-Avecm(i)*vecntA(j)/g
  103.        enddo
  104.       enddo
  105. CCCC
  106. C densi
  107. CCCC
  108. c Vecm y f perturbados
  109.  
  110. c      write(399,*)
  111. c      dxrel = 0.1
  112. c      do i = 1,10
  113. c      dx = xdensi * (dxrel**i)
  114. c      xdensi2 = xdensi + dx
  115. c      call carac_mate_rhmc_densi(XMAT,xdensi2)
  116. c      call yieldd(sig,ndims,void,ndimv,resuaux,nmodel)
  117. c      call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel)
  118. c      write(399,*)dxrel**i,(resuaux-resu)/dx
  119. c      write(399,*)'  ',((vecmaux(j)-vecm(j))/dx,j=1,3)
  120. c      call carac_mate_rhmc_densi(XMAT,xdensi)
  121. c      enddo
  122. c      do i = 1,10
  123. c      dx = - xdensi * (dxrel**i)
  124. c      xdensi2 = xdensi + dx
  125. c      call carac_mate_rhmc_densi(XMAT,xdensi2)
  126. c      call yieldd(sig,ndims,void,ndimv,resuaux,nmodel)
  127. c      call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel)
  128. c      write(399,*)-dxrel**i,(resuaux-resu)/dx
  129. c      write(399,*)'  ',((vecmaux(j)-vecm(j))/dx,j=1,3)
  130. c      call carac_mate_rhmc_densi(XMAT,xdensi)
  131. c      enddo
  132. CC
  133.       dxrel = 0.1D0
  134.       dx = (1.D0 + xdensi) * (dxrel**6)
  135.       xdensi2 = xdensi + dx
  136. c      call carac_mate_rhmc_densi(XMAT,xdensi2)
  137.       call carac_mate_densi(XMAT,xdensi2,nmodel)
  138.       if (ndimx.eq.3) then
  139.        call yieldd(sig,ndims,void,ndimv,resuaux,nmodel)
  140.        call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel)
  141.       else if (ndimx.eq.4) then
  142.        call yieldd(sig,ndims,x(ndimx),ndimv,resuaux,nmodel)
  143.        call vflsig(sig,ndims,x(ndimx),ndimv,vecmaux,nmodel)
  144.        call vflvar(sig,ndims,x(ndimx),ndimv,vecmaux(ndimx),nmodel)
  145.       endif
  146.       xdensi2 = xdensi - dx
  147. c      call carac_mate_rhmc_densi(XMAT,xdensi2)
  148.       call carac_mate_densi(XMAT,xdensi2,nmodel)
  149.       if (ndimx.eq.3) then
  150.        call yieldd(sig,ndims,void,ndimv,resu,nmodel)
  151.        call vflsig(sig,ndims,void,ndimv,vecm,nmodel)
  152.       else if (ndimx.eq.4) then
  153.        call yieldd(sig,ndims,x(ndimx),ndimv,resu,nmodel)
  154.        call vflsig(sig,ndims,x(ndimx),ndimv,vecm,nmodel)
  155.        call vflvar(sig,ndims,x(ndimx),ndimv,vecm(ndimx),nmodel)
  156.       endif
  157.       dx = dx * 2.D0
  158. c      call carac_mate_rhmc_densi(XMAT,xdensi)
  159.       call carac_mate_densi(XMAT,xdensi,nmodel)
  160. c Mmat = derivada de vecm (3 veces)
  161.       do i=1,ndimx
  162.           Mmat(i,1)=(vecmaux(i)-vecm(i))/dx
  163.           Mmat(i,2)=Mmat(i,1)
  164.           Mmat(i,3)=Mmat(i,1)
  165.       enddo
  166. c Kmat2 = Kmat * Mmat
  167. c Kmat3 = Avecm/g (3 veces)
  168.       do i=1,ndimx
  169.        do j=1,3
  170.          do k=1,ndimx
  171.           Kmat2(i,j)=Kmat2(i,j)+Kmat(i,k)*Mmat(k,j)
  172.          enddo
  173.          Kmat3(i,j)=Avecm(i)/g
  174.        enddo
  175.       enddo
  176. c auxiliares
  177.       aux1= xdensi*lam
  178.       aux2= xdensi*(resuaux-resu)/dx
  179. c EKmat  = E*Kmat ; EKmat2 = E*Kmat2 ; EKmat3 = E*Kmat3
  180. c EKmat  = EKmat + aux1*EKmat2 + aux2*EKmat3
  181.       call zzero(ekmat,ndimk**2)
  182.       do i=1,3
  183.        do j=1,3
  184.          do k=1,ndimx
  185.           ekmat(i,j) =ekmat(i,j)+
  186.      .                eemat(i+(k-1)*ndimx)*kmat(k,j)
  187.           ekmat2(i,j)=ekmat2(i,j)+
  188.      .                eemat(i+(k-1)*ndimx)*kmat2(k,j)
  189.           ekmat3(i,j)=ekmat3(i,j)+
  190.      .                eemat(i+(k-1)*ndimx)*kmat3(k,j)
  191.          enddo
  192.          ekmat(i,j)=ekmat(i,j)+aux1*ekmat2(i,j)+aux2*ekmat3(i,j)
  193.        enddo
  194.       enddo
  195.  
  196.       if (nescri.eq.1)
  197.      .   write(ues,999)' EKmat_dp ',((ekmat(i,j),j=1,3),i=1,3)
  198. 999   format(2x,A12,/,3(3(1x,E12.6),/))
  199.       return
  200.       end
  201.  
  202.  
  203.  
  204.  
  205.  
  206.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales