求问 error#6366，matmul问题

marlene

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

      module K_Hill48Module
        implicit none;

        ! 定义函数重载接口，后缀为dp为双精度函数，后缀为sp为单精度函数
        interface K_getEquStress
            module procedure K_getEquStress_dp;
            module procedure K_getEquStress_sp;
        end interface

        interface K_getDelF_DelSigma
            module procedure K_getDelF_DelSigma_dp;
            module procedure K_getDelF_DelSigma_sp;
        end interface

        interface K_getDelF_DelSigma_fast
            module procedure K_getDelF_DelSigma_fast_dp;
            module procedure K_getDelF_DelSigma_fast_sp;
        end interface

      
        

      contains


      
      
      ! Stress为行向量，3个分量
      real(8) function K_getP(param_r0, param_r45, param_r90)
        implicit none;
        real(8) :: param_r0;
        real(8) :: param_r45;
        real(8) :: param_r90;
        real(8) :: P12;
        real(8) :: P22;
        real(8) :: P66;
        dimension K_getP(3, 3);

        P12 = param_r0 / (1.d0 + param_r90);
        P22 = param_r0 *(1.d0 + param_r90) / (param_r90 
     &                               * (1.d0 + param_r0));
        P66 = (1.d0 +2.d0 * param_r45) * (param_r0 + param_r90)  
     &              / ((param_r0 + 1.d0)*param_r90);
         ! reshape前的数组构造器必须是列优先的
         K_getP = reshape((/  1.d0,           -P12,      0d0,
     &                        -P12,         P22,      0d0,
     &                        0d0,          0d0,      P66/)
     &                      , (/3, 3/));
       
      end function





      ! 求Hill48等效应力（双精度）
      ! Stress为行向量，3个分量
      real(8) function K_getEquStress_dp(param_r0, param_r45, param_r90
     &                                   ,Stress)
        implicit none;
        real(8) :: param_r0;
        real(8) :: param_r45;
        real(8) :: param_r90;
        real(8) :: P(3, 3);
        real(8) :: Stress_T(1, 3);
        real(8) :: Stress_V(3, 1);
        real(8) :: Stress(3);
        real(8) :: temp1(1, 3);
        real(8) :: temp2;
        
        P = K_getP(param_r0, param_r45, param_r90);
         Stress_T = reshape((/Stress(1), Stress(2), Stress(3)/)
     &                          , (/1, 3/));
         Stress_V = reshape((/Stress(1), Stress(2), Stress(3)/)
     &                          , (/3, 1/));
        temp1 = matmul(Stress_T, P);
        temp2 = matmul(temp1, Stress_V);
        K_getEquStress_dp = sqrt(temp2 / 2.d0);
        
      end function


      ! 求Hill48等效应力（单精度）
      ! Stress为行向量，3个分量
      real function K_getEquStress_sp(param_r0_sp, param_r45_sp,
     &                                param_r90_sp, Stress_sp)
        implicit none;
        real :: param_r0_sp;
        real :: param_r45_sp;
        real :: param_r90_sp;
        real :: Stress_sp(3);
        real(8) :: Stress_dp(3);
        real(8) :: equStress_dp;
        real(8) :: param_r0_dp;
        real(8) :: param_r45_dp;
        real(8) :: param_r90_dp;
        
        param_r0_dp = param_r0_sp;
        param_r45_dp = param_r45_sp;
        param_r90_dp = param_r90_sp;
        Stress_dp = Stress_sp;
        equStress_dp = K_getEquStress_dp(param_r0_dp, param_r45_dp
     &                              , param_r90_dp, Stress_dp);
        K_getEquStress_sp = equStress_dp;
      end function

      ! 求屈服函数f对应力Sigma(k)的偏导数df/dSigma(k); k = 1, 2, 3（双精度）
      ! Stress为行向量，3个分量
      ! df/dSigma也存储为行向量
      real(8) function K_getDelF_DelSigma_dp(param_r0, param_r45,
     &                                    param_r90, Stress)
        implicit none;
        real(8) :: param_r0;
        real(8) :: param_r45;
        real(8) :: param_r90;
        real(8) :: Stress_V(3, 1);
        real(8) :: Stress(3);
        real(8) :: P(3, 3);
        real(8) :: equStress;
        real(8) :: temp;
        dimension K_getDelF_DelSigma_dp(3,1);

        P = K_getP(param_r0, param_r45, param_r90);
        Stress_V = reshape((/Stress(1), Stress(2), Stress(3)/)
     &                          , (/3, 1/));
        equStress = K_getEquStress_dp(param_r0, param_r45, param_r90,
     &                                       Stress);
        temp = 2.d0 * equStress;
        K_getDelF_DelSigma_dp = matmul(P, Stress_V) / temp;
      end function

      ! 求屈服函数f对应力Sigma(k)的偏导数df/dSigma(k); k = 1, 2, 3（单精度）
      ! Stress为行向量，3个分量
      ! df/dSigma也存储为行向量
      real function K_getDelF_DelSigma_sp(param_r0_sp, param_r45_sp,
     &                                    param_r90_sp, Stress_sp)
        implicit none;
        real :: param_r0_sp;
        real :: param_r45_sp;
        real :: param_r90_sp;
        real :: Stress_sp(3);
        real(8) :: Stress_dp(3);
        real(8) :: param_r0_dp;
        real(8) :: param_r45_dp;
        real(8) :: param_r90_dp;
        real(8) :: DelF_DelSigma_dp(3,1);
        dimension K_getDelF_DelSigma_sp(3,1);
        param_r0_dp = param_r0_sp;
        param_r45_dp = param_r45_sp;
        param_r90_dp = param_r90_sp;

        Stress_dp = Stress_sp;

        DelF_DelSigma_dp = K_getDelF_DelSigma_dp(param_r0_dp, 
     &        param_r45_dp, param_r90_dp, Stress_dp);
        K_getDelF_DelSigma_sp = DelF_DelSigma_dp;
      end function


      ! 快速求屈服函数f对应力Sigma(k)的偏导数df/dSigma(k); k = 1, 2, 3（双精度）
      ! 给定equStress，减少计算量
      ! Stress为行向量，3个分量
      ! df/dSigma也存储为行向量
      real(8) function K_getDelF_DelSigma_fast_dp(param_r0, param_r45, 
     &                                  param_r90, Stress, equStress)
        implicit none;
        real(8) :: param_r0;
        real(8) :: param_r45;
        real(8) :: param_r90;
        real(8) :: Stress_V(3, 1);
        real(8) :: Stress(3);
        real(8) :: P(3, 3);
        real(8) :: equStress;
        real(8) :: temp;
        dimension K_getDelF_DelSigma_fast_dp(3,1);

        P = K_getP(param_r0, param_r45, param_r90);
        Stress_V = reshape((/Stress(1), Stress(2), Stress(3)/)
     &                          , (/3, 1/));
        temp = 2.d0 * equStress;
        
        K_getDelF_DelSigma_fast_dp = matmul(P, Stress_V) / temp;
      end function


      ! 快速求屈服函数f对应力Sigma(k)的偏导数df/dSigma(k); k = 1, 2, 3（单精度）
      ! 给定equStress，减少计算量
      ! Stress采用Voigt标记，为行向量，3个分量
      ! df/dSigma也存储为行向量
      real function K_getDelF_DelSigma_fast_sp(param_r0_sp, 
     &         param_r45_sp, param_r90_sp  ,Stress_sp, equStress_sp)
        implicit none;
        real :: param_r0_sp;
        real :: param_r45_sp;
        real :: param_r90_sp;
        real :: Stress_sp(3);
        real :: equStress_sp;
        
        real(8) :: Stress_dp(3);
        real(8) :: param_r0_dp;
        real(8) :: param_r45_dp;
        real(8) :: param_r90_dp;
        real(8) :: equStress_dp;
        real(8) :: DelF_DelSigma_dp(3,1);
        dimension K_getDelF_DelSigma_fast_sp(3,1);
        param_r0_dp = param_r0_sp;
        param_r45_dp = param_r45_sp;
        param_r90_dp = param_r90_sp;
        Stress_dp = Stress_sp;
        equStress_dp = equStress_sp;

        DelF_DelSigma_dp = K_getDelF_DelSigma_fast_dp(param_r0_dp,
     &      param_r45_dp, param_r90_dp, Stress_dp, equStress_dp);
        K_getDelF_DelSigma_fast_sp = DelF_DelSigma_dp;
      end function


     


      end module K_Hill48Module

错误    1     error #6366: The shapes of the array expressions do not conform.   [TEMP2]    求问大神们，在算等效应力时temp2为什么会出错呢   

marlene

是因为matmul只能赋给数组，而real（8）：：temp2是一个变量吗

Transpose

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

   real(8) :: P(3, 3);
   real(8) :: Stress_T(1, 3);
   real(8) :: Stress_V(3, 1);
   real(8) :: Stress(3);
   real(8) :: temp1(1, 3);
   real(8) :: temp2;
   temp1 = matmul(Stress_T, P);
   temp2 = matmul(temp1, Stress_V);

维度不匹配
(1 ,3) * (3, 1) =(1,1)
是一个二维数组，不是一个数

在使用matmul时，如果输入一维数数组，Fortran会自动判断
对于c=matmul(a,b)
如果a是一维数组，b是二维数组，会自动解析为 a(1,n) b(n,m) 返回一维数组 c(m)
如果a是二维数组，b是一维数组，会自动解析为 a(m,n) b(n,1) 返回一维数组 c(m)
不需要定义为 a(1,3)这种数组
最后的是两个向量乘法，使用dot_product即可

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

temp2=dot_product(matmul(stress,P),stress)

marlene

仔细看了一下，发现matmul必须有一个是二维数组，我这个是两个一维数组，需要dot_product,但是用了dot_product,temp2 = dot_product(temp1, Stress_V);
error #6372: One-dimensional array-valued arguments are required in this context.   [DOT_PRODUCT]       

风平老涡

marlene 发表于 2021-7-1 10:56
仔细看了一下，发现matmul必须有一个是二维数组，我这个是两个一维数组，需要dot_product,但是用了dot_prod ...

因为你的temp1和stress_V定义的都是二维数组，你可以定义temp2为二维数组temp2(1,1)，用matmul；或者temp1和stress_V赋值给二个一维数组，用dot_product；又或者用temp2=dot_product(temp1(1,:), stress_v(:, 1))

marlene

Transpose 发表于 2021-7-1 10:39
[mw_shl_code=fortran,true]   real(8) :: P(3, 3);
   real(8) :: Stress_T(1, 3);
   real(8) :: Stress_ ...

谢谢大佬，讲的很清楚，我明白啦

marlene

风平老涡 发表于 2021-7-1 21:14
因为你的temp1和stress_V定义的都是二维数组，你可以定义temp2为二维数组temp2(1,1)，用matmul；或者temp1 ...

谢谢大佬，讲的很清楚，我明白啦