Fortran Coder

查看: 1082|回复: 6

[求助] 求问 error#6366,matmul问题

[复制链接]

5

帖子

1

主题

0

精华

入门

F 币
38 元
贡献
16 点
发表于 2021-7-1 09:56:34 | 显示全部楼层 |阅读模式
[Fortran] 纯文本查看 复制代码
      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为什么会出错呢   



回复

使用道具 举报

5

帖子

1

主题

0

精华

入门

F 币
38 元
贡献
16 点
 楼主| 发表于 2021-7-1 09:58:25 | 显示全部楼层
是因为matmul只能赋给数组,而real(8)::temp2是一个变量吗

44

帖子

1

主题

1

精华

专家

Vim

F 币
466 元
贡献
233 点

规矩勋章

发表于 2021-7-1 10:39:38 | 显示全部楼层
本帖最后由 Transpose 于 2021-7-1 10:46 编辑

[Fortran] 纯文本查看 复制代码
   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] 纯文本查看 复制代码
temp2=dot_product(matmul(stress,P),stress)


5

帖子

1

主题

0

精华

入门

F 币
38 元
贡献
16 点
 楼主| 发表于 2021-7-1 10:56:34 | 显示全部楼层
仔细看了一下,发现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]       

148

帖子

2

主题

0

精华

宗师

F 币
1299 元
贡献
670 点

规矩勋章

发表于 2021-7-1 21:14:36 | 显示全部楼层
本帖最后由 风平老涡 于 2021-7-1 21:20 编辑
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))

5

帖子

1

主题

0

精华

入门

F 币
38 元
贡献
16 点
 楼主| 发表于 2021-7-3 16:29:44 | 显示全部楼层
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_ ...

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

5

帖子

1

主题

0

精华

入门

F 币
38 元
贡献
16 点
 楼主| 发表于 2021-7-3 16:30:19 | 显示全部楼层
风平老涡 发表于 2021-7-1 21:14
因为你的temp1和stress_V定义的都是二维数组,你可以定义temp2为二维数组temp2(1,1),用matmul;或者temp1 ...

谢谢大佬,讲的很清楚,我明白啦
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2021-9-24 10:37

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表