Abaqus VUMAT error during compilation

Domv · ‎01-17-2024

I need help to solve the following problem, i understood there is an array already allocated, but i don't know how find it, can someone help me?

forrtl: severe (151): allocatable array is already allocated
Image PC Routine Line Source
libifcoremd.dll 00007FF97DE2CC0E Unknown Unknown Unknown
explicitU.dll 00007FF94DCE6F49 Unknown Unknown Unknown
explicitU.dll 00007FF94DCE2E1E Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92CFE2347 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92D08438D Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92D0ABD6A Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92C6F2874 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92C6CA58C Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92C5BB0E8 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E16B832 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E16033C Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E1A63AB Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E1A8F2F Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E19BA86 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E19E314 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E1B08F9 Unknown Unknown Unknown
package.exe 00007FF78275FD09 Unknown Unknown Unknown
package.exe 00007FF78275DD86 Unknown Unknown Unknown
package.exe 00007FF782761210 Unknown Unknown Unknown
KERNEL32.DLL 00007FF9B0AD244D Unknown Unknown Unknown
ntdll.dll 00007FF9B272DFB8 Unknown Unknown Unknown

@Steve_Lionel

andrew_4619 · ‎01-17-2024

well the error message explains itself. Why not post your VUMAT routine and suggestions on how to fix it might be possible then

Domv · ‎01-17-2024

subroutine vumat(
C Read only (unmodifiable)variables -
1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
2 stepTime, totalTime, dt, cmname, coordMp, charLength,
3 props, density, strainInc, relSpinInc,
4 tempOld, stretchOld, defgradOld, fieldOld,
5 stressOld, stateOld, enerInternOld, enerInelasOld,
6 tempNew, stretchNew, defgradNew, fieldNew,
C Write only (modifiable) variables -
7 stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
C
dimension props(nprops), density(nblock), coordMp(nblock,*),
1 charLength(nblock), strainInc(nblock,ndir+nshr),
2 relSpinInc(nblock,nshr), tempOld(nblock),
3 stretchOld(nblock,ndir+nshr),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
8 defgradNew(nblock,ndir+nshr+nshr),
9 fieldNew(nblock,nfieldv),
1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
2 enerInternNew(nblock), enerInelasNew(nblock)
C
dimension F_new(3,3), xU(3,3), U_1(3,3), rot(3,3),
3 C_new(3,3), a0(1,3), b0(1,3), a0_old(1,3), b0_old(1,3),
4 F_old(3,3), DWDI4(1,5), DWDI8(1,5), DWDI10(1,5),
5 vectI4(1,5),vectI8(1,5),vectI10(1,5),
6 vectI4_full(1,5),vectI8_full(1,5),vectI10_full(1,5),
7 S(3,3), sigma_co(3,3), sigma_new(3,3), vectW(1,12),
8 vect_Wtot(1,12), D(3,3), a0xa0(3,3),
3 b0xb0(3,3), a0xb0(3,3), b0xa0(3,3), aS1(3,3),
4 aS3(3,3), DW(1,3)

c real*8 a_a, a_b, b_a, b_b, d_a, d_b, pmin_a, pmin_b,
cc + c_a, c_b, s1, s2, s3, sT_a, sT_b, sC_a, sC_b, rho, xJU,
cc + aI4, aI6, aI8, aI10, WTa_bar, WTb_bar, kTC_a, kTC_b,
c + WTa, WCa, WTb, WCb, WS, WSTa, WSTb, kSCa, kSCb, W_tot

parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
3 four=4.d0,
2 third=one/three, half= .5d0, twoThirds=two/three,
3 threeHalfs = 1.5d0, fourth=one/four, six=6.d0)

character*80 cmname
C
do i = 1,nblock

C Read material properties

a_a=props(1)

a_b=a_a

b_a=props(2)

b_b=b_a

d_a=props(3)

d_b=d_a

pmin_a=props(4)

pmin_b=pmin_a

c_a=props(5)

c_b=c_a

s1=props(6)

s2=props(7)

s3=props(8)

sT_a=props(9)

sT_b=sT_a

sC_a=props(10)

sC_b=sC_a

c ----------

c Deformation gradiEnt before the increment

F_old(1,1)=defgradOld(i,1)
F_old(1,2)=defgradOld(i,4)
F_old(1,3)=defgradOld(i,9)
F_old(2,1)=defgradOld(i,7)
F_old(2,2)=defgradOld(i,2)
F_old(2,3)=defgradOld(i,5)
F_old(3,1)=defgradOld(i,6)
F_old(3,2)=defgradOld(i,8)
F_old(3,3)=defgradOld(i,3)

c Deformation gradiant after the increment

F_new(1,1)=defgradNew(i,1)
F_new(1,2)=defgradNew(i,4)
F_new(1,3)=defgradNew(i,9)
F_new(2,1)=defgradNew(i,7)
F_new(2,2)=defgradNew(i,2)
F_new(2,3)=defgradNew(i,5)
F_new(3,1)=defgradNew(i,6)
F_new(3,2)=defgradNew(i,8)
F_new(3,3)=defgradNew(i,3)

c Stretch tensor after the increment

xU(1,1)=stretchNew(i,1)
xU(2,2)=stretchNew(i,2)
xU(3,3)=stretchNew(i,3)
xU(1,2)=stretchNew(i,4)
xU(2,3)=stretchNew(i,5)
xU(3,1)=stretchNew(i,6)
xU(2,1)=xU(1,2)
xU(3,2)=xU(2,3)
xU(1,3)=xU(3,1)

c Inverse of Stretch

CALL CalDet(xU,xJU)
CALL CalInv(xU,xJU,U_1)

c We need to calculate the rotational matrix: rot=FU^-1
c
c Rotation Matrix

rot=matmul(F_new,U_1)

c Right cauchy stress tensor

C_new=matmul(TRANSPOSE(F_new),F_new)

c Inizial fiber direction

a0_old(1,1)=one
a0_old(1,2)=zero
a0_old(1,3)=zero

b0_old(1,1)=zero
b0_old(1,2)=one
b0_old(1,3)=zero

a0=matmul(a0_old,F_old)
b0=matmul(b0_old,F_old)

c Calcolo invariante I4

aI4=sum(matmul(a0,C_new)*a0)

c Calcolo invariante I6

aI6=sum(matmul(a0,C_new)*b0)

c Calcolo invariante I8

aI8=sum(matmul(b0,C_new)*b0)

c Calcolo invariante I10

aI10=acos(sum(a0*b0))-acos(aI6/sqrt(aI4*aI8))

C -----------------------

WTa_bar=a_a*(SQRT(aI4)-1)**three+b_a*(SQRT(aI4)-one)**two

WTb_bar=a_b*(SQRT(aI8)-one)**three+b_b*(SQRT(aI8)-one)**two

kTC_a=(one-pmin_a)*exp(c_a*(SQRT(aI8)-one))+pmin_a

kTC_b=(one-pmin_b)*exp(c_b*(SQRT(aI4)-one))+pmin_b

WTa=WTa_bar*kTC_a

WCa=d_a*(SQRT(aI4)-one)**two

WTb=WTb_bar*kTC_b

WCb=d_b*(SQRT(aI8)-one)**two

WS=half*s1*aI10**two+third*s2*aI10**two+fourth*s3*aI10**four

WSTa=third*aI10**three*sT_a*(aI4-one)**two

WSTb=third*aI10**three*sT_b*(aI8-one)**two

kSCa=exp(sC_a*(aI4-one)**three)

kSCb=exp(sC_b*(aI8-one)**three)

C Definizione vettore vectW

vectW(1,1)=WTa
vectW(1,2)=WTb
vectW(1,3)=WS
vectW(1,4)=WSTa
vectW(1,5)=WSTb
vectW(1,6)=WCb
vectW(1,7)=WS*kSCb
vectW(1,8)=WSTa*kSCb
vectW(1,9)=WCa
vectW(1,10)=WS*kSCa
vectW(1,11)=WSTb*kSCa
vectW(1,12)=WS*kSCa*kSCb

c Definizione vettori vectI4_full:

C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b

vectI4_full(1,1)=(a_a*(threehalfs*SQRT(aI4)-
3 three+threehalfs*one/SQRT(aI4))+b_a*(one-
1 one/SQRT(aI4)))*((one-pmin_a)*exp(c_a*(SQRT(aI8)-one))-one)

vectI4_full(1,2)=(a_b*(SQRT(aI8)-1)**three+
1 b_b*(SQRT(aI8)-one)**two)*(one-
3 pmin_b)*(c_b*(half*one/SQRT(aI4)))*exp(c_b*(SQRT(aI4)-one))

vectI4_full(1,3)=(twothirds*aI10**three*sT_a*(aI4-one))*(one+kSCb)

vectI4_full(1,4)=(one-one/SQRT(aI4))*d_a

vectI4_full(1,5)=(WS+WSTb+WS*kSCb)*(three*aI4**one-six*aI4+
1 three)*exp(sC_a*(aI4**three-three*aI4**two+three*aI4-one))

c Definizione vettori vectI8_full:

C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b

vectI8_full(1,1)=(a_a*(SQRT(aI4)-one)**three+b_a*(SQRT(aI4)-
1 one)**two)*(one-pmin_a)*(half*c_a/SQRT(aI8))*exp(c_b*(SQRT(aI8)-
2 one))

vectI8_full(1,2)=(a_b*(threehalfs*SQRT(aI8)-
2 three+three/(two*SQRT(aI8)))+b_b*(one-
1 one/SQRT(aI8)))*((one-pmin_b)*exp(c_b*(SQRT(aI4)-one)+pmin_b))

vectI8_full(1,3)=(twothirds*aI10**three*sT_b*(aI8-one))*(one+kSCa)

vectI8_full(1,4)=d_b*(one-one/SQRT(aI8))

vectI8_full(1,5)=(WS+WSTa+WS*kSCa)*((three*aI8**two-six*aI8+
2 three)*exp(sC_b*(aI8**three-three*aI8**two+three*aI4-one)))

c Definizione vettori vectI10_full:

C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b

vectI10_full(1,1)=(s1*aI10+s2*aI10**two+s3*aI10**three)*(one+kSCb+
1 kSCa+kSCa*kSCb)

vectI10_full(1,2)=(one+kSCb)*(aI10**two*sT_a*(SQRT(aI4)-one)**two)

vectI10_full(1,3)=(one+kSCa)*(aI10**two*sT_b*(SQRT(aI8)-one)**two)

C -----------------------
C IFCICLE W_tot=vectW(1,12)*vect_Wtot(1,12)
c
c vect_Wtot(,,,,,,,,,,,)
c W_tot=matmul(vectW(1,12)*vect_Wtot(1,12))
c
c vectI4=(,,,,)
c vectI8=(,,,,)
c vectI10=(,,)

c DWDI4=matmul(vectI4_full,transpose(vectI4))
c DWDI8=matmul(vectI8_full,transpose(vectI8))
c DWDI10=matmul(vectI10_full,transpose(vectI10))

C IF aI4.GT.one.AND.aI8.GT.one

IF(aI4.GT.one) THEN
IF(aI8.GT.one) THEN

c vect_Wtot=(0,1,1,1,1,0,0,0,0,0,0,0)
vect_Wtot(1,1)=zero
vect_Wtot(1,2)=one
vect_Wtot(1,3)=one
vect_Wtot(1,4)=one
vect_Wtot(1,5)=one
vect_Wtot(1,6)=zero
vect_Wtot(1,7)=zero
vect_Wtot(1,8)=zero
vect_Wtot(1,9)=zero
vect_Wtot(1,10)=zero
vect_Wtot(1,11)=zero
vect_Wtot(1,12)=zero

W_tot=sum(vectW*vect_Wtot)

c vectI4=(1,1,1,0,0)

vectI4(1,1)=one
vectI4(1,1)=one
vectI4(1,1)=one
vectI4(1,1)=zero
vectI4(1,1)=zero

c vectI8=(1,1,1,0,0)

vectI8(1,1)=one
vectI8(1,1)=one
vectI8(1,1)=one
vectI8(1,1)=zero
vectI8(1,1)=zero

c vectI10=(1,1,1)

vectI10(1,1)=one
vectI10(1,1)=one
vectI10(1,1)=one

DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8
DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10

END IF
END IF

IF(aI4.GT.one) THEN
IF(aI8.LT.one.OR.aI8.EQ.one) THEN

c vect_Wtot(1,0,0,0,0,1,1,1,0,0,0,0)

vect_Wtot(1,1)=one
vect_Wtot(1,2)=zero
vect_Wtot(1,3)=zero
vect_Wtot(1,4)=zero
vect_Wtot(1,5)=zero
vect_Wtot(1,6)=one
vect_Wtot(1,7)=one
vect_Wtot(1,8)=one
vect_Wtot(1,9)=zero
vect_Wtot(1,10)=zero
vect_Wtot(1,11)=zero
vect_Wtot(1,12)=zero

W_tot=sum(vectW*vect_Wtot)

c vectI4=(1,0,0,1,1)

vectI4(1,1)=one
vectI4(1,1)=zero
vectI4(1,1)=zero
vectI4(1,1)=one
vectI4(1,1)=one

c vectI8=(1,0,0,1,1)

vectI8(1,1)=one
vectI8(1,1)=zero
vectI8(1,1)=zero
vectI8(1,1)=one
vectI8(1,1)=one

c vectI10=(1,1,0)

vectI10(1,1)=one
vectI10(1,1)=one
vectI10(1,1)=zero

DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8
DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10

END IF
END IF

C IF aI4 LT and EQ one and one and aI8 GT one

IF(aI4.LT.one.OR.aI4.EQ.one) THEN

IF(aI8.GT.one) THEN

c vect_Wtot(0,1,0,0,0,0,0,0,1,1,1,0)

vect_Wtot(1,1)=zero
vect_Wtot(1,2)=one
vect_Wtot(1,3)=zero
vect_Wtot(1,4)=zero
vect_Wtot(1,5)=zero
vect_Wtot(1,6)=zero
vect_Wtot(1,7)=zero
vect_Wtot(1,8)=zero
vect_Wtot(1,9)=one
vect_Wtot(1,10)=one
vect_Wtot(1,11)=one
vect_Wtot(1,12)=zero

W_tot=sum(vectW*vect_Wtot)

c vectI4=(0,1,1,0,0)

vectI4(1,1)=zero
vectI4(1,1)=one
vectI4(1,1)=one
vectI4(1,1)=zero
vectI4(1,1)=zero

c vectI8=(0,1,0,1,1)

vectI8(1,1)=zero
vectI8(1,1)=one
vectI8(1,1)=zero
vectI8(1,1)=one
vectI8(1,1)=one

c vectI10=(1,0,1)

vectI10(1,1)=one
vectI10(1,1)=zero
vectI10(1,1)=one

DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8
DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10

END IF
END IF

c IF aI4 LT one and EQ one and aI8 LT one and EQ one

IF(aI4.LT.one.OR.aI4.EQ.one) THEN
IF(aI8.LT.one.OR.aI8.EQ.one) THEN

c vect_Wtot=(0,0,0,0,0,1,0,0,1,0,0,1)

vect_Wtot(1,1)=zero
vect_Wtot(1,2)=zero
vect_Wtot(1,3)=zero
vect_Wtot(1,4)=zero
vect_Wtot(1,5)=zero
vect_Wtot(1,6)=one
vect_Wtot(1,7)=zero
vect_Wtot(1,8)=zero
vect_Wtot(1,9)=one
vect_Wtot(1,10)=zero
vect_Wtot(1,11)=zero
vect_Wtot(1,12)=one

W_tot=sum(vectW*vect_Wtot)

c vectI4=(0,0,0,1,1)

vectI4(1,1)=zero
vectI4(1,1)=zero
vectI4(1,1)=zero
vectI4(1,1)=one
vectI4(1,1)=one

c vectI8=(0,0,0,1,1)

vectI8(1,1)=zero
vectI8(1,1)=zero
vectI8(1,1)=zero
vectI8(1,1)=one
vectI8(1,1)=one

c vectI10=(1,0,0)

vectI10(1,1)=one
vectI10(1,1)=zero
vectI10(1,1)=zero

DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8
DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10

END IF
END IF

C -----------------------

c tensor products

CALL KronProd(a0,a0,a0xa0)

CALL KronProd(b0,b0,b0xb0)

CALL KronProd(a0,b0,a0xb0)

CALL KronProd(b0,b0,b0xa0)

aS1=a0xa0*DW(1,1)+b0xb0*DW(1,2)
aS2=DW(1,3)*(one/SQRT(aI4*aI8-aI6**two))
aS3=aI6*(a0xa0/aI4+b0xb0/aI8)

S=two*aS1+aS2*(a0xb0 + b0xa0 - aS3)

sigma_co=one/aJ*matmul(F_new,matmul(S,transpose(F_new)))

sigma_new=matmul(matmul(transpose(rot),sigma_co),rot)

stressNew(i,1)=sigma_new(1,1)
stressNew(i,2)=sigma_new(2,2)
stressNew(i,3)=sigma_new(3,3)
stressNew(i,4)=sigma_new(1,2)
stressNew(i,5)=sigma_new(2,3)
stressNew(i,6)=sigma_new(3,1)

c stateNew(i,1)=W_tot/rho

c deallocate(C)

return
end do

end

c calcolo del determinante

subroutine CalDet(A,aJ)

implicit none
REAL*8 A(3,3), aJ

aJ = A(1,1) * A(2,2) * A(3,3) -
1 A(1,2) * A(2,1) * A(3,3) +
2 + A(1,2) * A(2,3) * A(3,1) +
1 A(1,3) * A(3,2) * A(2,1) -
2 A(1,3) * A(3,1) * A(2,2) -
3 A(2,3) * A(3,2) * A(1,1)

return
end

c calcolo della matrice inversa

subroutine CalInv(A,aJ,B)
include 'vaba_param.inc'
dimension A(3,3),B(3,3)

B(1,1)=A(2,2)*A(3,3)-A(2,3)*A(3,2)
B(1,2)=A(1,3)*A(3,2)-A(1,2)*A(3,3)
B(1,3)=A(1,2)*A(2,3)-A(1,3)*A(2,2)
B(2,1)=A(2,3)*A(3,1)-A(2,1)*A(3,3)
B(2,2)=A(1,1)*A(3,3)-A(1,3)*A(3,1)
B(2,3)=A(1,3)*A(2,1)-A(1,1)*A(2,3)
B(3,1)=A(2,1)*A(3,2)-A(2,2)*A(3,1)
B(3,2)=A(1,2)*A(3,1)-A(1,1)*A(3,2)
B(3,3)=A(1,1)*A(2,2)-A(1,2)*A(2,1)

B=B/aJ

return
end

c calcolo prodotto tensioriale

function KronProd(A,B) result(C)
IMPLICIT NONE
real, dimension (:,:), intent(in) :: A, B
real, dimension (:,:), allocatable :: C
integer :: i = 0, j = 0, k = 0, l = 0
integer :: m = 0, n = 0, p = 0, q = 0

allocate(C(size(A,1)*size(B,1),size(A,2)*size(B,2)))
C = 0

do i = 1,size(A,1)
do j = 1,size(A,2)
n=(i-1)*size(B,1) + 1
m=n+size(B,1) - 1
p=(j-1)*size(B,2) + 1
q=p+size(B,2) - 1
C(n:m,p:q) = A(i,j)*B
enddo
enddo

end function KronProd

Domv · ‎01-17-2024

     subroutine vumat(
C Read only (unmodifiable)variables -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2  stepTime, totalTime, dt, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     5  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only (modifiable) variables -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
C
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     1  charLength(nblock), strainInc(nblock,ndir+nshr),
     2  relSpinInc(nblock,nshr), tempOld(nblock),
     3  stretchOld(nblock,ndir+nshr),
     4  defgradOld(nblock,ndir+nshr+nshr),
     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(nblock),
     8  stretchNew(nblock,ndir+nshr),
     8  defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv),
     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock)
C
         dimension F_new(3,3), xU(3,3), U_1(3,3), rot(3,3),
     3  C_new(3,3), a0(1,3), b0(1,3), a0_old(1,3), b0_old(1,3),
     4  F_old(3,3), DWDI4(1,5), DWDI8(1,5), DWDI10(1,5),
     5  vectI4(1,5),vectI8(1,5),vectI10(1,5),
     6  vectI4_full(1,5),vectI8_full(1,5),vectI10_full(1,5),
     7  S(3,3), sigma_co(3,3), sigma_new(3,3), vectW(1,12), 
     8  vect_Wtot(1,12), D(3,3), a0xa0(3,3),
     3  b0xb0(3,3), a0xb0(3,3), b0xa0(3,3), aS1(3,3),
     4  aS3(3,3), DW(1,3)
     
         
c         real*8  a_a, a_b, b_a, b_b, d_a, d_b, pmin_a, pmin_b,
cc     +   c_a, c_b, s1, s2, s3, sT_a, sT_b, sC_a, sC_b, rho, xJU,
cc    +   aI4, aI6, aI8, aI10, WTa_bar, WTb_bar, kTC_a, kTC_b,
c     +   WTa, WCa, WTb, WCb, WS, WSTa, WSTb, kSCa, kSCb, W_tot   

         
         parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
     3 four=4.d0,
     2 third=one/three, half= .5d0, twoThirds=two/three, 
     3 threeHalfs = 1.5d0, fourth=one/four, six=6.d0)
      
      character*80 cmname
C
      do i = 1,nblock

C     Read material properties
      
      a_a=props(1)

      a_b=a_a

      b_a=props(2)

      b_b=b_a

      d_a=props(3)

      d_b=d_a

      pmin_a=props(4)

      pmin_b=pmin_a

      c_a=props(5)

      c_b=c_a

      s1=props(6)

      s2=props(7)

      s3=props(8)

      sT_a=props(9)

      sT_b=sT_a

      sC_a=props(10)

      sC_b=sC_a


c      ----------
      
c     Deformation gradiEnt before the increment
          
      F_old(1,1)=defgradOld(i,1)
      F_old(1,2)=defgradOld(i,4)
      F_old(1,3)=defgradOld(i,9)
      F_old(2,1)=defgradOld(i,7)
      F_old(2,2)=defgradOld(i,2)
      F_old(2,3)=defgradOld(i,5)
      F_old(3,1)=defgradOld(i,6)
      F_old(3,2)=defgradOld(i,8)
      F_old(3,3)=defgradOld(i,3)

c     Deformation gradiant after the increment      
      
      F_new(1,1)=defgradNew(i,1)
      F_new(1,2)=defgradNew(i,4)
      F_new(1,3)=defgradNew(i,9)
      F_new(2,1)=defgradNew(i,7)
      F_new(2,2)=defgradNew(i,2)
      F_new(2,3)=defgradNew(i,5)
      F_new(3,1)=defgradNew(i,6)
      F_new(3,2)=defgradNew(i,8)
      F_new(3,3)=defgradNew(i,3)
      
c     Stretch tensor after the increment      
      
      xU(1,1)=stretchNew(i,1)
      xU(2,2)=stretchNew(i,2)
      xU(3,3)=stretchNew(i,3)
      xU(1,2)=stretchNew(i,4)
      xU(2,3)=stretchNew(i,5)
      xU(3,1)=stretchNew(i,6)
      xU(2,1)=xU(1,2)
      xU(3,2)=xU(2,3)
      xU(1,3)=xU(3,1)
      
c     Inverse of Stretch
      
      CALL CalDet(xU,xJU)
      CALL CalInv(xU,xJU,U_1)

c     We need to calculate the rotational matrix: rot=FU^-1
c
c     Rotation Matrix

      rot=matmul(F_new,U_1)
      
c     Right cauchy stress tensor
      
      C_new=matmul(TRANSPOSE(F_new),F_new)
      
c     Inizial fiber direction
      
      a0_old(1,1)=one
      a0_old(1,2)=zero
      a0_old(1,3)=zero

      b0_old(1,1)=zero
      b0_old(1,2)=one
      b0_old(1,3)=zero
     
      a0=matmul(a0_old,F_old)
      b0=matmul(b0_old,F_old)
      

c     Calcolo invariante I4      
      
      aI4=sum(matmul(a0,C_new)*a0)
      
c     Calcolo invariante I6
   
      aI6=sum(matmul(a0,C_new)*b0)
      
c     Calcolo invariante I8
   
      aI8=sum(matmul(b0,C_new)*b0)
      
c     Calcolo invariante I10
      
      aI10=acos(sum(a0*b0))-acos(aI6/sqrt(aI4*aI8))
      
C     ----------------------- 
      
      WTa_bar=a_a*(SQRT(aI4)-1)**three+b_a*(SQRT(aI4)-one)**two

      WTb_bar=a_b*(SQRT(aI8)-one)**three+b_b*(SQRT(aI8)-one)**two

      kTC_a=(one-pmin_a)*exp(c_a*(SQRT(aI8)-one))+pmin_a

      kTC_b=(one-pmin_b)*exp(c_b*(SQRT(aI4)-one))+pmin_b

      WTa=WTa_bar*kTC_a

      WCa=d_a*(SQRT(aI4)-one)**two

      WTb=WTb_bar*kTC_b

      WCb=d_b*(SQRT(aI8)-one)**two 

      WS=half*s1*aI10**two+third*s2*aI10**two+fourth*s3*aI10**four

      WSTa=third*aI10**three*sT_a*(aI4-one)**two

      WSTb=third*aI10**three*sT_b*(aI8-one)**two

      kSCa=exp(sC_a*(aI4-one)**three)

      kSCb=exp(sC_b*(aI8-one)**three)
      
C     Definizione vettore vectW
       
      vectW(1,1)=WTa
      vectW(1,2)=WTb
      vectW(1,3)=WS
      vectW(1,4)=WSTa
      vectW(1,5)=WSTb
      vectW(1,6)=WCb
      vectW(1,7)=WS*kSCb
      vectW(1,8)=WSTa*kSCb
      vectW(1,9)=WCa
      vectW(1,10)=WS*kSCa
      vectW(1,11)=WSTb*kSCa
      vectW(1,12)=WS*kSCa*kSCb
      
c     Definizione vettori vectI4_full:    
      
C      a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b
      
      vectI4_full(1,1)=(a_a*(threehalfs*SQRT(aI4)-
     3 three+threehalfs*one/SQRT(aI4))+b_a*(one-
     1 one/SQRT(aI4)))*((one-pmin_a)*exp(c_a*(SQRT(aI8)-one))-one)
      
      vectI4_full(1,2)=(a_b*(SQRT(aI8)-1)**three+
     1 b_b*(SQRT(aI8)-one)**two)*(one-
     3 pmin_b)*(c_b*(half*one/SQRT(aI4)))*exp(c_b*(SQRT(aI4)-one))
      
      vectI4_full(1,3)=(twothirds*aI10**three*sT_a*(aI4-one))*(one+kSCb)
      
      vectI4_full(1,4)=(one-one/SQRT(aI4))*d_a
      
      vectI4_full(1,5)=(WS+WSTb+WS*kSCb)*(three*aI4**one-six*aI4+
     1 three)*exp(sC_a*(aI4**three-three*aI4**two+three*aI4-one))

c      Definizione vettori vectI8_full: 

C      a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b


      vectI8_full(1,1)=(a_a*(SQRT(aI4)-one)**three+b_a*(SQRT(aI4)-
     1 one)**two)*(one-pmin_a)*(half*c_a/SQRT(aI8))*exp(c_b*(SQRT(aI8)-
     2 one))
      
      vectI8_full(1,2)=(a_b*(threehalfs*SQRT(aI8)-
     2 three+three/(two*SQRT(aI8)))+b_b*(one-
     1 one/SQRT(aI8)))*((one-pmin_b)*exp(c_b*(SQRT(aI4)-one)+pmin_b))
      
      vectI8_full(1,3)=(twothirds*aI10**three*sT_b*(aI8-one))*(one+kSCa)
      
      vectI8_full(1,4)=d_b*(one-one/SQRT(aI8)) 
      
      vectI8_full(1,5)=(WS+WSTa+WS*kSCa)*((three*aI8**two-six*aI8+
     2 three)*exp(sC_b*(aI8**three-three*aI8**two+three*aI4-one)))
      
      
c     Definizione vettori  vectI10_full: 
      
C      a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b
      
      vectI10_full(1,1)=(s1*aI10+s2*aI10**two+s3*aI10**three)*(one+kSCb+
     1 kSCa+kSCa*kSCb)
      
      vectI10_full(1,2)=(one+kSCb)*(aI10**two*sT_a*(SQRT(aI4)-one)**two)
      
      vectI10_full(1,3)=(one+kSCa)*(aI10**two*sT_b*(SQRT(aI8)-one)**two)
      
C     -----------------------      
C     IFCICLE W_tot=vectW(1,12)*vect_Wtot(1,12) 
c      
c     vect_Wtot(,,,,,,,,,,,) 
c     W_tot=matmul(vectW(1,12)*vect_Wtot(1,12)) 
c      
c     vectI4=(,,,,) 
c     vectI8=(,,,,) 
c     vectI10=(,,)
      
c     DWDI4=matmul(vectI4_full,transpose(vectI4))
c     DWDI8=matmul(vectI8_full,transpose(vectI8))
c     DWDI10=matmul(vectI10_full,transpose(vectI10))

C    IF aI4.GT.one.AND.aI8.GT.one 
      
      IF(aI4.GT.one) THEN
          IF(aI8.GT.one) THEN
            
c      vect_Wtot=(0,1,1,1,1,0,0,0,0,0,0,0)
      vect_Wtot(1,1)=zero      
      vect_Wtot(1,2)=one
      vect_Wtot(1,3)=one
      vect_Wtot(1,4)=one
      vect_Wtot(1,5)=one
      vect_Wtot(1,6)=zero
      vect_Wtot(1,7)=zero
      vect_Wtot(1,8)=zero
      vect_Wtot(1,9)=zero
      vect_Wtot(1,10)=zero
      vect_Wtot(1,11)=zero
      vect_Wtot(1,12)=zero
      
      W_tot=sum(vectW*vect_Wtot)
      
c      vectI4=(1,1,1,0,0) 
      
      vectI4(1,1)=one
      vectI4(1,1)=one
      vectI4(1,1)=one
      vectI4(1,1)=zero
      vectI4(1,1)=zero
      
c      vectI8=(1,1,1,0,0) 
      
      vectI8(1,1)=one
      vectI8(1,1)=one
      vectI8(1,1)=one
      vectI8(1,1)=zero
      vectI8(1,1)=zero
      
c      vectI10=(1,1,1)
      
      vectI10(1,1)=one
      vectI10(1,1)=one
      vectI10(1,1)=one
      
      DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
      DW(1,2)=sum(vectI8_full*vectI8)  ! DWDI8
      DW(1,3)=sum(vectI10_full*vectI10)  ! DWDI10
      
          END IF
      END IF 
      
      IF(aI4.GT.one) THEN
          IF(aI8.LT.one.OR.aI8.EQ.one) THEN
              
c     vect_Wtot(1,0,0,0,0,1,1,1,0,0,0,0) 
              
      vect_Wtot(1,1)=one      
      vect_Wtot(1,2)=zero
      vect_Wtot(1,3)=zero
      vect_Wtot(1,4)=zero
      vect_Wtot(1,5)=zero
      vect_Wtot(1,6)=one
      vect_Wtot(1,7)=one
      vect_Wtot(1,8)=one
      vect_Wtot(1,9)=zero
      vect_Wtot(1,10)=zero
      vect_Wtot(1,11)=zero
      vect_Wtot(1,12)=zero        
              
      W_tot=sum(vectW*vect_Wtot) 
      
c      vectI4=(1,0,0,1,1) 
      
      vectI4(1,1)=one
      vectI4(1,1)=zero
      vectI4(1,1)=zero
      vectI4(1,1)=one
      vectI4(1,1)=one
      
c    vectI8=(1,0,0,1,1) 
      
      vectI8(1,1)=one
      vectI8(1,1)=zero
      vectI8(1,1)=zero
      vectI8(1,1)=one
      vectI8(1,1)=one
      
c      vectI10=(1,1,0)
      
      vectI10(1,1)=one
      vectI10(1,1)=one
      vectI10(1,1)=zero
      
      DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
      DW(1,2)=sum(vectI8_full*vectI8)  ! DWDI8
      DW(1,3)=sum(vectI10_full*vectI10)  ! DWDI10
      
          END IF
      END IF
      
C      IF aI4 LT and EQ one and one and aI8 GT one 

      IF(aI4.LT.one.OR.aI4.EQ.one) THEN
          
          IF(aI8.GT.one) THEN      

c      vect_Wtot(0,1,0,0,0,0,0,0,1,1,1,0) 

      vect_Wtot(1,1)=zero      
      vect_Wtot(1,2)=one
      vect_Wtot(1,3)=zero
      vect_Wtot(1,4)=zero
      vect_Wtot(1,5)=zero
      vect_Wtot(1,6)=zero
      vect_Wtot(1,7)=zero
      vect_Wtot(1,8)=zero
      vect_Wtot(1,9)=one
      vect_Wtot(1,10)=one
      vect_Wtot(1,11)=one
      vect_Wtot(1,12)=zero
      
      
      W_tot=sum(vectW*vect_Wtot) 
       
c      vectI4=(0,1,1,0,0) 
       
      vectI4(1,1)=zero
      vectI4(1,1)=one
      vectI4(1,1)=one
      vectI4(1,1)=zero
      vectI4(1,1)=zero
      
      
c      vectI8=(0,1,0,1,1) 
      
      vectI8(1,1)=zero
      vectI8(1,1)=one
      vectI8(1,1)=zero
      vectI8(1,1)=one
      vectI8(1,1)=one
      
c      vectI10=(1,0,1)
      
      vectI10(1,1)=one
      vectI10(1,1)=zero
      vectI10(1,1)=one
      
      DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
      DW(1,2)=sum(vectI8_full*vectI8)  ! DWDI8
      DW(1,3)=sum(vectI10_full*vectI10)  ! DWDI10

          END IF
      END IF      
      
c      IF aI4 LT one and EQ one and aI8 LT one and EQ one
      
      IF(aI4.LT.one.OR.aI4.EQ.one) THEN
          IF(aI8.LT.one.OR.aI8.EQ.one) THEN
              
c      vect_Wtot=(0,0,0,0,0,1,0,0,1,0,0,1) 
              
      vect_Wtot(1,1)=zero      
      vect_Wtot(1,2)=zero
      vect_Wtot(1,3)=zero
      vect_Wtot(1,4)=zero
      vect_Wtot(1,5)=zero
      vect_Wtot(1,6)=one
      vect_Wtot(1,7)=zero
      vect_Wtot(1,8)=zero
      vect_Wtot(1,9)=one
      vect_Wtot(1,10)=zero
      vect_Wtot(1,11)=zero
      vect_Wtot(1,12)=one 
              
      W_tot=sum(vectW*vect_Wtot) 
       
c      vectI4=(0,0,0,1,1) 
       
      vectI4(1,1)=zero
      vectI4(1,1)=zero
      vectI4(1,1)=zero
      vectI4(1,1)=one
      vectI4(1,1)=one
      
c      vectI8=(0,0,0,1,1) 
       
      vectI8(1,1)=zero
      vectI8(1,1)=zero
      vectI8(1,1)=zero
      vectI8(1,1)=one
      vectI8(1,1)=one
      
c      vectI10=(1,0,0)
      
      vectI10(1,1)=one
      vectI10(1,1)=zero
      vectI10(1,1)=zero
       
      DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
      DW(1,2)=sum(vectI8_full*vectI8)  ! DWDI8
      DW(1,3)=sum(vectI10_full*vectI10)  ! DWDI10
      
          END IF
      END IF      
      
C     ----------------------- 
      
c     tensor products  
      
      
      CALL KronProd(a0,a0,a0xa0)
      
      CALL KronProd(b0,b0,b0xb0)
        
      CALL KronProd(a0,b0,a0xb0)
            
      CALL KronProd(b0,b0,b0xa0)   
      

      aS1=a0xa0*DW(1,1)+b0xb0*DW(1,2)
      aS2=DW(1,3)*(one/SQRT(aI4*aI8-aI6**two))
      aS3=aI6*(a0xa0/aI4+b0xb0/aI8)
      
      S=two*aS1+aS2*(a0xb0 + b0xa0 - aS3)

      sigma_co=one/aJ*matmul(F_new,matmul(S,transpose(F_new)))
      
      sigma_new=matmul(matmul(transpose(rot),sigma_co),rot)
      
      stressNew(i,1)=sigma_new(1,1)
      stressNew(i,2)=sigma_new(2,2)
      stressNew(i,3)=sigma_new(3,3)
      stressNew(i,4)=sigma_new(1,2)
      stressNew(i,5)=sigma_new(2,3)
      stressNew(i,6)=sigma_new(3,1)
       
c      stateNew(i,1)=W_tot/rho
       
c      deallocate(C)
      
      
      return
      end do
      

      end  

c     calcolo del determinante      
      
      subroutine CalDet(A,aJ)

      implicit none 
      REAL*8 A(3,3), aJ

      aJ = A(1,1) * A(2,2) * A(3,3) -
     1 A(1,2) * A(2,1) * A(3,3) +
     2 + A(1,2) * A(2,3) * A(3,1) +
     1 A(1,3) * A(3,2) * A(2,1) -
     2 A(1,3) * A(3,1) * A(2,2) -
     3 A(2,3) * A(3,2) * A(1,1)  
      
      return 
      end 
      
c    calcolo della matrice inversa
      
      subroutine CalInv(A,aJ,B)
      include 'vaba_param.inc'
      dimension A(3,3),B(3,3)
      
      B(1,1)=A(2,2)*A(3,3)-A(2,3)*A(3,2)
      B(1,2)=A(1,3)*A(3,2)-A(1,2)*A(3,3)
      B(1,3)=A(1,2)*A(2,3)-A(1,3)*A(2,2)
      B(2,1)=A(2,3)*A(3,1)-A(2,1)*A(3,3)
      B(2,2)=A(1,1)*A(3,3)-A(1,3)*A(3,1)
      B(2,3)=A(1,3)*A(2,1)-A(1,1)*A(2,3)
      B(3,1)=A(2,1)*A(3,2)-A(2,2)*A(3,1)
      B(3,2)=A(1,2)*A(3,1)-A(1,1)*A(3,2)
      B(3,3)=A(1,1)*A(2,2)-A(1,2)*A(2,1)
      
      B=B/aJ
      
      return
      end
      
c     calcolo prodotto tensioriale
      
      function KronProd(A,B) result(C)
       IMPLICIT NONE
       real, dimension (:,:), intent(in)  :: A, B
       real, dimension (:,:), allocatable :: C
       integer :: i = 0, j = 0, k = 0, l = 0
       integer :: m = 0, n = 0, p = 0, q = 0

       allocate(C(size(A,1)*size(B,1),size(A,2)*size(B,2)))
       C = 0
  
       do i = 1,size(A,1)
        do j = 1,size(A,2)
         n=(i-1)*size(B,1) + 1
         m=n+size(B,1) - 1
         p=(j-1)*size(B,2) + 1
         q=p+size(B,2) - 1
         C(n:m,p:q) = A(i,j)*B
        enddo
       enddo
       
      end function KronProd

andrew_4619 · ‎01-17-2024

Well it appears that the message "forrtl: severe (151): allocatable array is already allocated" is not generated by the code you pasted, though it maybe your code sets things that cause other parts of the ABAQUS solver to crash. I see many bugs... For example

c      vectI4=(1,1,1,0,0) 
      
      vectI4(1,1)=one
      vectI4(1,1)=one
      vectI4(1,1)=one
      vectI4(1,1)=zero
      vectI4(1,1)=zero
      
c      vectI8=(1,1,1,0,0) 
      
      vectI8(1,1)=one
      vectI8(1,1)=one
      vectI8(1,1)=one
      vectI8(1,1)=zero
      vectI8(1,1)=zero
      
c      vectI10=(1,1,1)
      
      vectI10(1,1)=one
      vectI10(1,1)=one
      vectI10(1,1)=one

you are setting vectI4(1,1) five times as this array is declared (1,5) you really want:

c      vectI4=(1,1,1,0,0) 
      
      vectI4(1,1)=one
      vectI4(1,2)=one
      vectI4(1,3)=one
      vectI4(1,4)=zero
      vectI4(1,5)=zero
....

but is would make more sense to have

      vectI4(1,:) = [ one, one, one, zero, zero ]
      vectI8(1,:) = [ one, one, one, zero, zero ]   
      vectI10(1,1:3) =  one
c note vectI10(1,4) and (1,5) are not initialised is that intentional

Fortran has syntax for whole array and array slice operations.

Domv · ‎01-17-2024

I really don't know how i made this mistake, i've already correct and still receiving the same error, with abaqus message:

ERROR in job messaging system: Error in connection to analysis.

Error in job Job-1: Abaqus/Explicit Packager exited with an error - Please see the status file for possible error messages if the file exists.
Job Job-1 aborted due to errors.

andrew_4619 · ‎01-17-2024

There are I suspect more errors in the code I only had a quick look at it. You can get some Fortran help here but not much help with ABAQUS. Did you start by having a working sample VUMAT that ran OK? That should be the first test.