[Fortran] 纯文本查看 复制代码
#include "fintrf.h"
subroutine mexFunction(OutSum,OutVar,InSum,InVar) !函数接口名称必须为mexFunction,
!OutSum:输出参数个数
!OutVar:输出参数数组指针
!InSum:输入参数个数
!InVar:输入参数数组指针
!参数顺序不能随意更改
Integer InSum,OutSum
mwPointer InVar(*),OutVar(*) !mwPointer专门用于表示指针变量,这个不能随意用Integer代替
mwPointer mxGetPr, mxGetDimensions, mxCreateNumericArray !这个对返回指针函数的再次申明,
integer, parameter :: dp = 8
Integer , parameter :: myINT = SELECTED_INT_KIND(8)
Real(dp),Allocatable :: P(:,:,:,:), F1(:,:,:,:), F2(:,:,:,:), Gamma(:,:,:,:), rc(:,:,:,:), uind(:,:,:,:) ! input and output uind
Integer,Allocatable :: Dim_P(:), Dim_F(:), Dim_G(:), Dim_rc(:) ! input array dimension
Real(dp) :: n,co ! input n and cut_off distance
Real(dp) :: uindx,uindy,uindz,temuindx,temuindy,temuindz
Integer I1,J1,K1,I2,J2,K2,NBP,NTP,NSP,NBF,NTF,NSF,PD,FD,GD,RD
If(InSum/=7)Then
call mexErrMsgIdAndTxt('MATLAB:InputTooBig','输入参数个数必须为7个')
Return
EndIf
PD = mxGetNumberOfDimensions(P)
FD = mxGetNumberOfDimensions(F1)
GD = mxGetNumberOfDimensions(Gamma)
RD = mxGetNumberOfDimensions(rc)
Allocate(Dim_P(PD), Dim_F(FD), Dim_G(GD), Dim_rc(RD))
Dim_P = mxGetDimensions(InVar(1)) !获取第1个输入参数P的维度
Dim_F = mxGetDimensions(InVar(2)) !获取第2 (3)个输入参数F1 (F2)的维度
Dim_G = mxGetDimensions(InVar(4)) !获取第4个输入参数Gamma的维度
Dim_rc = mxGetDimensions(InVar(5)) !获取第5个输入参数rc的维度
Allocate(P(Dim_P(1),Dim_P(2),Dim_P(3),Dim_P(4)))
Allocate(F1(Dim_F(1),Dim_F(2),Dim_F(3),Dim_F(4)))
Allocate(F2(Dim_F(1),Dim_F(2),Dim_F(3),Dim_F(4)))
Allocate(Gamma(Dim_G(1),Dim_G(2),Dim_G(3),Dim_G(4)))
Allocate(rc(Dim_rc(1),Dim_rc(2),Dim_rc(3),Dim_rc(4)))
Allocate(Uind(Dim_P(1),Dim_P(2),Dim_P(3),Dim_P(4)))
Call mxCopyPtrToReal8(mxGetPr(InVar(1)),P,Dim_P(1)*Dim_P(2)*Dim_P(3)*Dim_P(4)) !将第1个参数数组赋值给P变量
Call mxCopyPtrToReal8(mxGetPr(InVar(2)),F1,Dim_F(1)*Dim_F(2)*Dim_F(3)*Dim_F(4)) !将第2个参数数组赋值给F1变量
Call mxCopyPtrToReal8(mxGetPr(InVar(3)),F2,Dim_F(1)*Dim_F(2)*Dim_F(3)*Dim_F(4)) !将第3个参数数组赋值给F2变量
Call mxCopyPtrToReal8(mxGetPr(InVar(4)),Gamma,Dim_G(1)*Dim_G(2)*Dim_G(3)*Dim_G(4)) !将第4个参数数组赋值给Gamma变量
Call mxCopyPtrToReal8(mxGetPr(InVar(5)),rc,Dim_rc(1)*Dim_rc(2)*Dim_rc(3)*Dim_rc(4)) !将第5个参数数组赋值给rc变量
Call mxCopyPtrToReal8(mxGetPr(InVar(6)),n,1)!将第6个整数变量赋值给n
Call mxCopyPtrToReal8(mxGetPr(InVar(7)),co,1)!将第7个整数变量赋值给co cut_off_distance
NBP = Dim_P(1)
NTP = Dim_P(3)
NSP = Dim_P(4)
NBF = Dim_F(1)
NTF = Dim_F(3)
NSF = Dim_F(4)
DO K1 = 1,NBP
DO J1 = 1,NTP
DO I1 = 1,NSP
temuindx = 0.0
temuindy = 0.0
temuindz = 0.0
DO K2 = 1,NBF
DO J2 = 1,NTF
DO I2 = 1,NSF
call Biot_p2p(P(I1,1,J1,K1),P(I1,2,J1,K1),P(I1,3,J1,K1),F1(I2,1,J2,K2),F1(I2,2,J2,K2),F1(I2,3,J2,K2),F2(I2,1,J2,K2),F2(I2,2,J2,K2),F2(I2,3,J2,K2),&
Gamma(I1,1,J1,K1),rc(I1,1,J1,K1),co,n,temuindx,temuindy,temuindz)
temuindx = temuindx + temuindx
temuindy = temuindy + temuindy
temuindy = temuindy + temuindy
ENDDO
ENDDO
ENDDO
Uind(I1,1,J1,K1) = temuindx
Uind(I1,2,J1,K1) = temuindy
Uind(I1,3,J1,K1) = temuindz
ENDDO
ENDDO
ENDDO
OutVar(1)=mxCreateNumericArray(PD,Dim_P,REAL*8,0)!给返回参数分配内存
Call mxCopyReal8ToPtr(Uind,mxGetPr(OutVar(1)),Dim_P(1)*Dim_P(2)*Dim_P(3)*Dim_P(4))!将返回参数赋值给分配的内存
DeAllocate(P,F1,F2,Gamma,rc)!释放临时分配的内存
Return
End SubRoutine
Subroutine Biot_p2p(px,py,pz,f1x,f1y,f1z,f2x,f2y,f2z,gamma,rc,co,n,temuindx,temuindy,temuindz)
implicit none
integer, parameter :: dp = 8
Real(dp),Intent(In):: px,py,pz,f1x,f1y,f1z,f2x,f2y,f2z,gamma
Real(dp),Intent(Out):: temuindx,temuindy,temuindz
Real(dp) :: ldx, ldy, ldz, pxx1, pxx2, pyy1, pyy2, pzz1, pzz2, r1, r2, ldr, len, r1dr2, r1tr2, cnu, den, ubar,cmn,rc
Real :: co,cutoff,cored,n
Real, parameter :: pi = acos(-1.0)
Real, parameter :: infinity = 999999.0
cutoff = co
cored = n
temuindx = 0.0
temuindy = 0.0
temuindz = 0.0
ldx = f2x-f1x
ldy = f2y-f1y
ldz = f2z-f1z
pxx1 = px-f1x
pxx2 = px-f2x
pyy1 = py-f1y
pyy2 = py-f2y
pzz1 = pz-f1z
pzz2 = pz-f2z
cmn = cored
r1 = sqrt(pxx1*pxx1+pyy1*pyy1+pzz1*pzz1)
r2 = sqrt(pxx2*pxx2+pyy2*pyy2+pzz2*pzz2)
if((r1<cutoff).AND.(r2<cutoff)) then
ldr = (ldx*pxx1 + ldy*pxx2 + ldz*pzz1)
ldr = ldr*ldr
len = ldx*ldx + ldy*ldy + ldz*ldz !!L^2
r1dr2 = pxx1*pxx2+pyy1*pyy2+pzz1*pzz2
r1tr2 = r1*r2
cnu = (r1*r1)-(ldr/len)
cnu = cnu*((rc**cmn+cnu**cmn)**(-1.0/cmn))
den = 4.0*pi*(r1tr2*(r1tr2 + r1dr2))
ubar = cnu*gamma*(r1+r2)
ubar = ubar/den
if (isnan(ubar)) then
temuindx = 0.0
temuindy = 0.0
temuindz = 0.0
else if (ubar.GE.infinity) then
temuindx = 0.0
temuindy = 0.0
temuindz = 0.0
else
temuindx = ubar*(pyy1*pzz2-pzz1*pyy2)
temuindy = ubar*(pzz1*pxx2-pxx1*pzz2)
temuindz = ubar*(pxx1*pyy2-pyy1*pxx2)
endif
endif
End SubRoutine Biot_p2p
for errors. Please consult the External Interfaces Guide for information
on debugging MEX-files.