Next: DC3D.F Up: Listing of key Previous: ELAS3D.F
c *********************** thermal3d.f ***************************
c BACKGROUND
c Program adjusts dimensions of unit cell,
c [(1 + macrostrain) times dimension],
c in response to phases that have a non-zero eigenstrain and arbitrary
c elastic moduli tensors.
c All six macrostrains can adjust their values (3-d program), and are
c stored in the last two positions in the displacement vector u,
c as listed below. Periodic boundaries are maintained.
c In the comments below, (USER) means that this is a section of code
c that the user might have to change for his particular problem.
c Therefore the user is encouraged to search for this string.
c PROBLEM AND VARIABLE DEFINITION
c The problem being solved is the minimization of the elastic energy
c 1/2 uAu + bu + C + Tu + Y, where b and C are also functions of the
c macrostrains.
c The small array zcon computes the thermal strain energy associated
c with macrostrains (C term), T is the thermal energy term linear in the
c displacements (built from ss), b is the regular energy term linear in the
c b is the regular energy term linear in the
c displacements, u is the displacements including the macrostrains,
c gb is the energy gradient vector, h,Ah are auxiliary vectors,
c dk is the single pixel stiffness matrix, pix is the phase
c identification vector, and ib is the
c integer matrix for mapping labels from the 1-27 nearest neighbor
c labelling to the 1-d system labelling.
c The array prob(i) contains the volume fractions of the i'th phase,
c strxx, etc. are the six independent (Voigt notation) volume
c averaged stresses, sxx, etc. are the six independent (Voigt notation)
c volume averaged strains (not counting the thermal strains).
c The variable cmod(i,6,6) gives the elastic moduli tensor
c of the i'th phase, eigen(i,6) gives the six independent elements
c of the eigenstrain tensor for the i'th phase (Voigt notation)
c and dk(i,8,3,8,3) is the stiffness matrix of the i'th
c phase. The parameter nphase gives the number of phases being considered
c in the problem, and is set by the user.
c DIMENSIONS
c The main arrays of the problem, u, gb, h, Ah, b, and T, are dimensioned
c as (nx*ny*nz)+2, which is the number of nodal displacements plus two for
c the macrostrains.
c Currently the program assumes the number of different phases is
c 100, since phasemod and eigen (the moduli and eigenstrains for each phase)
c and dk are dimensioned to have at most 100 different kinds. This
c is easily changed, by changing the dimension of these three variables
c throughout the program. The parameter nphase gives the number of phases
c All major arrays are passed to subroutines via simple common statements.
c NOTE ON USE OF PROGRAM
c Program is set up to allow the macrostrains,
c which control the overall size of the system, to be dynamic
c variables, which are adjusted in order to minimize the overall
c energy. That means that if there are no eigenstrains specified
c for any of the phases, the overall strain will always relax to
c zero. If it is desired to simply apply a fixed strain, with no
c eigenstrains, then in subroutines Energy and Dembx, one must
c zero out the elements of gb (in Energy and in Dembx) that
c correspond to the macrostrains. This is easily done.
c This will fix the gradients of the macrostrains to always to be
c zero, so that they will not change, so the applied strain (initial
c values of the macrostrains) will remain fixed.
c STRONGLY SUGGESTED: READ MANUAL BEFORE USING PROGRAM!!!
c (USER) Change these dimensions and in other subroutines at same
c time. For example, search and replace all occurrences throughout
c the program of "(8002" by "(64000", to go from a 20 x 20 x 20
c system to a 40 x 40 x 40 system.
real u(8002,3),gb(8002,3),b(8002,3)
real h(8002,3),Ah(8002,3),T(8002,3)
real C,dk(100,8,3,8,3)
real cmod(100,6,6),ss(100,8,3),eigen(100,6)
real zcon(2,3,2,3),pk(6,8,3)
real phasemod(100,2),prob(100)
integer in(27),jn(27),kn(27)
integer*4 ib(8002,27)
integer*2 pix(8002)
common/list1/strxx,stryy,strzz,strxz,stryz,strxy
common/list2/h,Ah
common/list3/ib
common/list4/pix
common/list5/dk,b,C,zcon,Y
common/list6/u
common/list7/gb
common/list8/cmod,T,eigen
common/list10/phasemod,nphase,ss
common/list11/sxx,syy,szz,sxz,syz,sxy
c (USER) Unit 9 is the microstructure input file, unit 7
c is the results output file.
open (9,file='microstructure.dat')
open (7,file='outputfile.out')
c (USER) nx,ny,nz are the size of the lattice
nx=20
ny=20
nz=20
c ns=total number of sites
ns=nx*ny*nz
write(7,9010) nx,ny,nz,ns
9010 format(' nx= ',i4,' ny= ',i4,' nz= ',i4,' ns = ',i8)
c Add two more entries in the displacement vector for the 6 macrostrains,
c u(ns+1,1) = exx,u(ns+1,2) = eyy, u(ns+1,3) = ezz,
c u(nss,1) = exz, u(nss,2) = eyz, u(nss,3) = exy
nss=ns+2
c (USER) nphase is the number of phases being considered in the problem.
c The values of pix(m) will run from 1 to nphase.
nphase=2
c (USER) gtest is the stopping criterion, compared to gg=gb*gb.
c If gtest=abc*ns, when gg < gtest, the rms value per pixel
c of gb is less than sqrt(abc)
gtest=1.e-20*(nx*ny*nz)
write(7,*) 'relaxation criterion gtest = ',gtest
c (USER)
c The parameter phasemod(i,j) is the bulk (i,1) and shear (i,2) moduli of
c the i'th phase. These can be
c input in terms of Young's modulus E (i,1) and Poisson's ratio nu (i,2).
c The program, in the do 1144 loop, changes them to bulk and shear
c moduli, using relations for isotropic elastic moduli.
c For anisotropic moduli tensors, one can directly input the whole tensor
c cmod in subroutine femat, and skip this part.
c If you wish to input in terms of bulk (i,1) and shear (i,2) moduli,
c then simply comment out do 1144 loop.
phasemod(1,1)=1.0
phasemod(1,2)=0.2
phasemod(2,1)=1.0
phasemod(2,2)=0.2
do 1144 i=1,nphase
save=phasemod(i,1)
phasemod(i,1)=phasemod(i,1)/3./(1.-2.*phasemod(i,2))
phasemod(i,2)=save/2./(1.+phasemod(i,2))
1144 continue
c (USER) input eigenstrains for each phase
c (1=xx, 2=yy, 3=zz, 4=xz, 5=yz, 6=xy).
eigen(1,1)=0.
eigen(1,2)=0.
eigen(1,3)=0.
eigen(1,4)=0.
eigen(1,5)=0.
eigen(1,6)=0.
eigen(2,1)=0.1
eigen(2,2)=0.
eigen(2,3)=0.
eigen(2,4)=0.
eigen(2,5)=0.
eigen(2,6)=0.
c Construct the 27 neighbor table, ib(m,n)
c First construct the 27 neighbor table in terms of delta i, delta j,
c and delta k information (see Table 3 in manual)
in(1)=0
in(2)=1
in(3)=1
in(4)=1
in(5)=0
in(6)=-1
in(7)=-1
in(8)=-1
jn(1)=1
jn(2)=1
jn(3)=0
jn(4)=-1
jn(5)=-1
jn(6)=-1
jn(7)=0
jn(8)=1
do 555 n=1,8
kn(n)=0
kn(n+8)=-1
kn(n+16)=1
in(n+8)=in(n)
in(n+16)=in(n)
jn(n+8)=jn(n)
jn(n+16)=jn(n)
555 continue
in(25)=0
in(26)=0
in(27)=0
jn(25)=0
jn(26)=0
jn(27)=0
kn(25)=-1
kn(26)=1
kn(27)=0
c Now construct neighbor table according to 1-d labels
c Matrix ib(m,n) gives the 1-d label of the n'th neighbor (n=1,27) of
c the node labelled m.
nxy=nx*ny
do 1020 k=1,nz
do 1020 j=1,ny
do 1020 i=1,nx
m=nxy*(k-1)+nx*(j-1)+i
do 1004 n=1,27
i1=i+in(n)
j1=j+jn(n)
k1=k+kn(n)
if(i1.lt.1) i1=i1+nx
if(i1.gt.nx) i1=i1-nx
if(j1.lt.1) j1=j1+ny
if(j1.gt.ny) j1=j1-ny
if(k1.lt.1) k1=k1+nz
if(k1.gt.nz) k1=k1-nz
m1=nxy*(k1-1)+nx*(j1-1)+i1
ib(m,n)=m1
1004 continue
1020 continue
c Compute the average stress and strain, as well as the macrostrains (overall
c system size and shape) in each microstructure.
c (USER) npoints is the number of microstructures to use.
npoints=1
do 8000 micro=1,npoints
c Read in a microstructure in subroutine ppixel, and set up pix(m)
c with the appropriate phase assignments.
call ppixel(nx,ny,nz,ns,nphase)
c Count and output the volume fractions of the different phases
call assig(ns,nphase,prob)
do 8050 i=1,nphase
write(7,9065) i,prob(i)
9065 format(' Volume fraction of phase ',i3,' is ',f10.8)
8050 continue
c output elastic moduli (bulk and shear) for each phase
write(7,*) ' Phase Moduli'
do 111 i=1,nphase
write(7,9020) i,phasemod(i,1),phasemod(i,2)
9020 format(' Phase ',i3,' bulk = ',f12.6,' shear = ',f12.6)
111 continue
c output thermal strains for each phase
write(7,*) ' Thermal Strains'
do 119 i=1,nphase
write(7,9029) i,eigen(i,1),eigen(i,2),eigen(i,3)
write(7,9029) i,eigen(i,4),eigen(i,5),eigen(i,6)
9029 format('Phase ',i3,' ',3f6.2)
119 continue
c (USER) Set inital macrostrains of computational cell
u(ns+1,1)=0.0
u(ns+1,2)=0.0
u(ns+1,3)=0.0
u(nss,1)=0.0
u(nss,2)=0.0
u(nss,3)=0.0
c Apply homogeneous macroscopic strain as the initial condition
c to displacement variables
do 1050 k=1,nz
do 1050 j=1,ny
do 1050 i=1,nx
m=nxy*(k-1)+nx*(j-1)+i
x=float(i-1)
y=float(j-1)
z=float(k-1)
u(m,1)=x*u(ns+1,1)+y*u(nss,3)+z*u(nss,1)
u(m,2)=x*u(nss,3)+y*u(ns+1,2)+z*u(nss,2)
u(m,3)=x*u(nss,1)+y*u(nss,2)+z*u(ns+1,3)
1050 continue
c Set up the finite element stiffness matrices,the constant, C,
c the vector, b, required for the energy. b and C depend on the macrostrains.
c When they are updated, the values of b and C are updated too via
c calling subroutine femat.
c Only compute the thermal strain terms the first time femat is called,
c (iskip=0) as they are unaffected by later changes (iskip=1) in
c displacements and macrostrains.
c Compute initial value of gradient gb and gg=gb*gb.
iskip=0
call femat(nx,ny,nz,ns,iskip)
call energy(nx,ny,nz,ns,utot)
gg=0.0
do 100 m3=1,3
do 100 m=1,nss
gg=gg+gb(m,m3)*gb(m,m3)
100 continue
write(7,9042) utot,gg
9042 format(' energy = ',e15.8,' gg= ',e15.8)
call flush(7)
c Relaxation loop
c (USER) kmax is the maximum number of times that dembx will be called,
c with ldemb conjugate gradient steps performed during each call.
c The total number of conjugate gradient steps allowed for a given elastic
c computation is kmax*ldemb.
kmax=40
ldemb=50
ltot=0
do 5000 kkk=1,kmax
c Call dembx to implement conjugate gradient routine
write(7,*) 'Going into dembx, call no. ',kkk
call dembx(nx,ny,nz,ns,Lstep,gg,gtest,ldemb,kkk)
ltot=ltot+Lstep
c Call energy to compute energy after dembx call. If gg < gtest, this
c will be the final energy. If gg is still larger than gtest, then this
c will give an intermediate energy with which to check how the
c relaxation process is coming along. The call to energy does not
c change the gradient or the value of gg.
c Need to first call femat to update the vector b, as the value of the
c components of b depend on the macrostrains.
iskip=1
call femat(nx,ny,nz,ns,iskip)
call energy(nx,ny,nz,ns,utot)
write(7,9043) utot,gg,ltot
9043 format(' energy = ',e15.8,' gg= ',e15.8,' ltot = ',i6)
call flush(7)
c If relaxation process is finished, jump out of loop
if(gg.lt.gtest) goto 444
c Output stresses, strains, and macrostrains as an additional aid in judging
c how well the relaxation process is proceeding.
call stress(nx,ny,nz,ns)
write(7,*) ' stresses: xx,yy,zz,xz,yz,xy'
write(7,*) strxx,stryy,strzz,strxz,stryz,strxy
write(7,*) ' strains: xx,yy,zz,xz,yz,xy'
write(7,*) sxx,syy,szz,sxz,syz,sxy
write(7,*) ' macrostrains in same order'
write(7,*) u(ns+1,1),u(ns+1,2),u(ns+1,3)
write(7,*) u(nss,1),u(nss,2),u(nss,3)
write(7,*) 'avg = ',(u(ns+1,1)+u(ns+1,2)+u(ns+1,3))/3.
5000 continue
444 call stress(nx,ny,nz,ns)
write(7,*) ' stresses: xx,yy,zz,xz,yz,xy'
write(7,*) strxx,stryy,strzz,strxz,stryz,strxy
write(7,*) ' strains: xx,yy,zz,xz,yz,xy'
write(7,*) sxx,syy,szz,sxz,syz,sxy
write(7,*) ' macrostrains in same order'
write(7,*) u(ns+1,1),u(ns+1,2),u(ns+1,3)
write(7,*) u(nss,1),u(nss,2),u(nss,3)
write(7,*) 'avg = ',(u(ns+1,1)+u(ns+1,2)+u(ns+1,3))/3.
8000 continue
end
c Subroutine sets up the stiffness matrices, the linear term in the
c regular displacements, b, and the constant term, C, which come from
c the periodic boundary conditions, the term linear in the displacments,
c T, that comes from the thermal strains, and the constant term Y.
subroutine femat(nx,ny,nz,ns,iskip)
real u(8002,3),b(8002,3),T(8002,3)
real dk(100,8,3,8,3),phasemod(100,2),dndx(8),dndy(8),dndz(8)
real g(3,3,3),econ,ck(6,6),cmu(6,6),cmod(100,6,6)
real es(6,8,3),zcon(2,3,2,3),ss(100,8,3)
real eigen(100,6),delta(8,3)
integer is(8),iskip
integer*4 ib(8002,27)
integer*2 pix(8002)
common/list3/ib
common/list4/pix
common/list5/dk,b,C,zcon,Y
common/list6/u
common/list8/cmod,T,eigen
common/list10/phasemod,nphase,ss
nxy=nx*ny
c Generate dk, zcon, T, and Y on first pass. After that they are
c constant, snce they are independent of the macrostrains. Only b gets
c upgraded as the macrostrains change.
c Line number 1221 is the routine for b.
if(iskip.eq.1) goto 1221
c initialize stiffness matrices
do 40 m=1,nphase
do 40 l=1,3
do 40 k=1,3
do 40 j=1,8
do 40 i=1,8
dk(m,i,k,j,l)=0.0
40 continue
c initialize zcon matrix (gives C term for arbitrary macrostrains)
do 42 i=1,2
do 42 j=1,2
do 42 mi=1,3
do 42 mj=1,3
zcon(i,mi,j,mj)=0.0
42 continue
c (USER) An anisotropic elastic moduli tensor could be input at this point,
c bypassing this part, which assumes isotropic elasticity, so that there
c are only two independent numbers making up the elastic moduli tensor,
c the bulk modulus K and the shear modulus G.
c Set up elastic moduli matrices for each kind of element
c ck and cmu are the bulk modulus and shear modulus matrices, which
c need to multiplied by the actual bulk and shear moduli in each phase.
ck(1,1)=1.0
ck(1,2)=1.0
ck(1,3)=1.0
ck(1,4)=0.0
ck(1,5)=0.0
ck(1,6)=0.0
ck(2,1)=1.0
ck(2,2)=1.0
ck(2,3)=1.0
ck(2,4)=0.0
ck(2,5)=0.0
ck(2,6)=0.0
ck(3,1)=1.0
ck(3,2)=1.0
ck(3,3)=1.0
ck(3,4)=0.0
ck(3,5)=0.0
ck(3,6)=0.0
ck(4,1)=0.0
ck(4,2)=0.0
ck(4,3)=0.0
ck(4,4)=0.0
ck(4,5)=0.0
ck(4,6)=0.0
ck(5,1)=0.0
ck(5,2)=0.0
ck(5,3)=0.0
ck(5,4)=0.0
ck(5,5)=0.0
ck(5,6)=0.0
ck(6,1)=0.0
ck(6,2)=0.0
ck(6,3)=0.0
ck(6,4)=0.0
ck(6,5)=0.0
ck(6,6)=0.0
cmu(1,1)=4.0/3.0
cmu(1,2)=-2.0/3.0
cmu(1,3)=-2.0/3.0
cmu(1,4)=0.0
cmu(1,5)=0.0
cmu(1,6)=0.0
cmu(2,1)=-2.0/3.0
cmu(2,2)=4.0/3.0
cmu(2,3)=-2.0/3.0
cmu(2,4)=0.0
cmu(2,5)=0.0
cmu(2,6)=0.0
cmu(3,1)=-2.0/3.0
cmu(3,2)=-2.0/3.0
cmu(3,3)=4.0/3.0
cmu(3,4)=0.0
cmu(3,5)=0.0
cmu(3,6)=0.0
cmu(4,1)=0.0
cmu(4,2)=0.0
cmu(4,3)=0.0
cmu(4,4)=1.0
cmu(4,5)=0.0
cmu(4,6)=0.0
cmu(5,1)=0.0
cmu(5,2)=0.0
cmu(5,3)=0.0
cmu(5,4)=0.0
cmu(5,5)=1.0
cmu(5,6)=0.0
cmu(6,1)=0.0
cmu(6,2)=0.0
cmu(6,3)=0.0
cmu(6,4)=0.0
cmu(6,5)=0.0
cmu(6,6)=1.0
do 31 k=1,nphase
do 21 j=1,6
do 11 i=1,6
cmod(k,i,j)=phasemod(k,1)*ck(i,j)+phasemod(k,2)*cmu(i,j)
11 continue
21 continue
31 continue
c Set up Simpson's integration rule weight vector
do 30 k=1,3
do 30 j=1,3
do 30 i=1,3
nm=0
if(i.eq.2) nm=nm+1
if(j.eq.2) nm=nm+1
if(k.eq.2) nm=nm+1
g(i,j,k)=4.0**nm
30 continue
c Loop over the nphase kinds of pixels and
c Simpson's rule quadrature points in order to compute the stiffness
c matrices. Stiffness matrices of trilinear finite elements are quadratic
c in x, y, and z, so that Simpson's rule quadrature gives exact results.
do 4000 ijk=1,nphase
do 3000 k=1,3
do 3000 j=1,3
do 3000 i=1,3
x=float(i-1)/2.0
y=float(j-1)/2.0
z=float(k-1)/2.0
c dndx means the negative derivative with respect to x, of the shape
c matrix N (see manual, Sec. 2.2), dndy and dndz are similar.
dndx(1)=-(1.0-y)*(1.0-z)
dndx(2)=(1.0-y)*(1.0-z)
dndx(3)=y*(1.0-z)
dndx(4)=-y*(1.0-z)
dndx(5)=-(1.0-y)*z
dndx(6)=(1.0-y)*z
dndx(7)=y*z
dndx(8)=-y*z
dndy(1)=-(1.0-x)*(1.0-z)
dndy(2)=-x*(1.0-z)
dndy(3)=x*(1.0-z)
dndy(4)=(1.0-x)*(1.0-z)
dndy(5)=-(1.0-x)*z
dndy(6)=-x*z
dndy(7)=x*z
dndy(8)=(1.0-x)*z
dndz(1)=-(1.0-x)*(1.0-y)
dndz(2)=-x*(1.0-y)
dndz(3)=-x*y
dndz(4)=-(1.0-x)*y
dndz(5)=(1.0-x)*(1.0-y)
dndz(6)=x*(1.0-y)
dndz(7)=x*y
dndz(8)=(1.0-x)*y
c now build strain matrix
do 2799 n1=1,6
do 2799 n2=1,8
do 2799 n3=1,3
es(n1,n2,n3)=0.0
2799 continue
do 2797 n=1,8
es(1,n,1)=dndx(n)
es(2,n,2)=dndy(n)
es(3,n,3)=dndz(n)
es(4,n,1)=dndz(n)
es(4,n,3)=dndx(n)
es(5,n,2)=dndz(n)
es(5,n,3)=dndy(n)
es(6,n,1)=dndy(n)
es(6,n,2)=dndx(n)
2797 continue
c now do matrix multiply to determine value at (x,y,z), multiply by
c proper weight, and sum into dk, the stiffness matrix
do 900 mm=1,3
do 900 nn=1,3
do 900 ii=1,8
do 900 jj=1,8
c define sum over strain matrices and elastic moduli matrix for
c stiffness matrix
sum=0.0
do 890 kk=1,6
do 890 ll=1,6
sum=sum+es(kk,ii,mm)*cmod(ijk,kk,ll)*es(ll,jj,nn)
890 continue
dk(ijk,ii,mm,jj,nn)=dk(ijk,ii,mm,jj,nn)+g(i,j,k)*sum/216.
900 continue
3000 continue
4000 continue
c Now compute the ss matrices, which give the thermal strain terms
c for the i'th phase, single pixel.
dndx(1)=-0.25
dndx(2)=0.25
dndx(3)=0.25
dndx(4)=-0.25
dndx(5)=-0.25
dndx(6)=0.25
dndx(7)=0.25
dndx(8)=-0.25
dndy(1)=-0.25
dndy(2)=-0.25
dndy(3)=0.25
dndy(4)=0.25
dndy(5)=-0.25
dndy(6)=-0.25
dndy(7)=0.25
dndy(8)=0.25
dndz(1)=-0.25
dndz(2)=-0.25
dndz(3)=-0.25
dndz(4)=-0.25
dndz(5)=0.25
dndz(6)=0.25
dndz(7)=0.25
dndz(8)=0.25
c now build average strain matrix
do 3799 n1=1,6
do 3799 n2=1,8
do 3799 n3=1,3
es(n1,n2,n3)=0.0
3799 continue
do 3797 n=1,8
es(1,n,1)=dndx(n)
es(2,n,2)=dndy(n)
es(3,n,3)=dndz(n)
es(4,n,1)=dndz(n)
es(4,n,3)=dndx(n)
es(5,n,2)=dndz(n)
es(5,n,3)=dndy(n)
es(6,n,1)=dndy(n)
es(6,n,2)=dndx(n)
3797 continue
do 3598 mmm=1,nphase
do 3798 nn=1,3
do 3798 mm=1,8
sum=0.0
do 3698 nm=1,6
do 3698 n=1,6
sum=sum+cmod(mmm,n,nm)*es(n,mm,nn)*eigen(mmm,nm)
3698 continue
ss(mmm,mm,nn)=sum
3798 continue
3598 continue
c now call subroutine const to generate zcon
c zcon is a (2,3) x (2,3) matrix
call const(dk,ns,zcon,nx,ny,nz)
c Now set up linear term, T, for thermal energy. It does not depend
c on the macrostrains or displacements, so there is no need to update it
c as the macrostrains change. T is built up out of the ss matrices.
nss=ns+2
do 6066 m3=1,3
do 6066 m=1,nss
T(m,m3)=0.0
6066 continue
c For all cases, the correspondence between 1-8 finite element node
c labels and the 1-27 neighbor labels is (see Table 4 in manual):
c 1:ib(m,27), 2:ib(m,3),3:ib(m,2),4:ib(m,1)
c 5:ib(m,26),6:ib(m,19),7:ib(m,18),8:ib(m,17)
is(1)=27
is(2)=3
is(3)=2
is(4)=1
is(5)=26
is(6)=19
is(7)=18
is(8)=17
c Do all points, but no macrostrain terms
c note: factor of 2 on linear thermal term is cancelled
c by factor of 1/2 out in front of total energy term
do 6601 k=1,nz
do 6601 j=1,ny
do 6601 i=1,nx
m=nxy*(k-1)+nx*(j-1)+i
do 6600 mm=1,8
do 6600 nn=1,3
T(ib(m,is(mm)),nn)=T(ib(m,is(mm)),nn)-ss(pix(m),mm,nn)
6600 continue
6601 continue
c now need to pick up and sum in all terms multiplying macrostrains
do 7788 ipp=1,2
do 7788 jpp=1,3
exx=0.0
eyy=0.0
ezz=0.0
exz=0.0
eyz=0.0
exy=0.0
if(ipp.eq.1.and.jpp.eq.1) exx=1.0
if(ipp.eq.1.and.jpp.eq.2) eyy=1.0
if(ipp.eq.1.and.jpp.eq.3) ezz=1.0
if(ipp.eq.2.and.jpp.eq.1) exz=1.0
if(ipp.eq.2.and.jpp.eq.2) eyz=1.0
if(ipp.eq.2.and.jpp.eq.3) exy=1.0
c x=nx face
do 6001 i3=1,3
do 6001 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.3.or.i8.eq.6.or.i8.eq.7) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
6001 continue
do 6000 j=1,ny-1
do 6000 k=1,nz-1
m=nxy*(k-1)+j*nx
do 6900 nn=1,3
do 6900 mm=1,8
T(ns+ipp,jpp)=T(ns+ipp,jpp)-ss(pix(m),mm,nn)*delta(mm,nn)
6900 continue
6000 continue
c y=ny face
do 6011 i3=1,3
do 6011 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.3.or.i8.eq.4.or.i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
6011 continue
do 6010 i=1,nx-1
do 6010 k=1,nz-1
m=nxy*(k-1)+nx*(ny-1)+i
do 6901 nn=1,3
do 6901 mm=1,8
T(ns+ipp,jpp)=T(ns+ipp,jpp)-ss(pix(m),mm,nn)*delta(mm,nn)
6901 continue
6010 continue
c z=nz face
do 6021 i3=1,3
do 6021 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.5.or.i8.eq.6.or.i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
6021 continue
do 6020 i=1,nx-1
do 6020 j=1,ny-1
m=nxy*(nz-1)+nx*(j-1)+i
do 6902 nn=1,3
do 6902 mm=1,8
T(ns+ipp,jpp)=T(ns+ipp,jpp)-ss(pix(m),mm,nn)*delta(mm,nn)
6902 continue
6020 continue
c x=nx y=ny edge
do 6031 i3=1,3
do 6031 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.6) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.4.or.i8.eq.8) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.3.or.i8.eq.7) then
delta(i8,1)=exy*ny+exx*nx
delta(i8,2)=eyy*ny+exy*nx
delta(i8,3)=eyz*ny+exz*nx
end if
6031 continue
do 6030 k=1,nz-1
m=nxy*k
do 6903 nn=1,3
do 6903 mm=1,8
T(ns+ipp,jpp)=T(ns+ipp,jpp)-ss(pix(m),mm,nn)*delta(mm,nn)
6903 continue
6030 continue
c x=nx z=nz edge
do 6041 i3=1,3
do 6041 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.3) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.5.or.i8.eq.8) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.6.or.i8.eq.7) then
delta(i8,1)=exz*nz+exx*nx
delta(i8,2)=eyz*nz+exy*nx
delta(i8,3)=ezz*nz+exz*nx
end if
6041 continue
do 6040 j=1,ny-1
m=nxy*(nz-1)+nx*(j-1)+nx
do 6904 nn=1,3
do 6904 mm=1,8
T(ns+ipp,jpp)=T(ns+ipp,jpp)-ss(pix(m),mm,nn)*delta(mm,nn)
6904 continue
6040 continue
c y=ny z=nz edge
do 6051 i3=1,3
do 6051 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.5.or.i8.eq.6) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.3.or.i8.eq.4) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exy*ny+exz*nz
delta(i8,2)=eyy*ny+eyz*nz
delta(i8,3)=eyz*ny+ezz*nz
end if
6051 continue
do 6050 i=1,nx-1
m=nxy*(nz-1)+nx*(ny-1)+i
do 6905 nn=1,3
do 6905 mm=1,8
T(ns+ipp,jpp)=T(ns+ipp,jpp)-ss(pix(m),mm,nn)*delta(mm,nn)
6905 continue
6050 continue
c x=nx y=ny z=nz corner
do 6061 i3=1,3
do 6061 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.4) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.5) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.8) then
delta(i8,1)=exy*ny+exz*nz
delta(i8,2)=eyy*ny+eyz*nz
delta(i8,3)=eyz*ny+ezz*nz
end if
if(i8.eq.6) then
delta(i8,1)=exx*nx+exz*nz
delta(i8,2)=exy*nx+eyz*nz
delta(i8,3)=exz*nx+ezz*nz
end if
if(i8.eq.3) then
delta(i8,1)=exx*nx+exy*ny
delta(i8,2)=exy*nx+eyy*ny
delta(i8,3)=exz*nx+eyz*ny
end if
if(i8.eq.7) then
delta(i8,1)=exx*nx+exy*ny+exz*nz
delta(i8,2)=exy*nx+eyy*ny+eyz*nz
delta(i8,3)=exz*nx+eyz*ny+ezz*nz
end if
6061 continue
m=nx*ny*nz
do 6906 mm=1,8
do 6906 nn=1,3
T(ns+ipp,jpp)=T(ns+ipp,jpp)-ss(pix(m),mm,nn)*delta(mm,nn)
6906 continue
7788 continue
c now compute Y, the 0.5(eigen)Cij(eigen) energy, doesn't ever change
c with macrostrain or displacements
Y=0.0
do 8811 m=1,ns
do 8811 n=1,6
do 8811 nn=1,6
Y=Y+0.5*eigen(pix(m),n)*cmod(pix(m),n,nn)*eigen(pix(m),nn)
8811 continue
c Following needs to be run after every change in macrostrain
c when energy is recomputed.
1221 continue
c Use auxiliary variables (exx, etc.) instead of u() variable, for
c convenience, and to make the following code easier to read.
exx=u(ns+1,1)
eyy=u(ns+1,2)
ezz=u(ns+1,3)
exz=u(nss,1)
eyz=u(nss,2)
exy=u(nss,3)
c Now set up vector for linear term that comes from periodic boundary
c conditions. Notation and conventions same as for T term.
c This is done using the stiffness matrices, and the periodic terms
c in the macrostrains. It is easier to set up b this way than to
c analytically write out all the terms involved.
do 5000 m3=1,3
do 5000 m=1,ns
b(m,m3)=0.0
5000 continue
c For all cases, the correspondence between 1-8 finite element node
c labels and the 1-27 neighbor labels is (see Table 4 in manual):
c 1:ib(m,27), 2:ib(m,3),3:ib(m,2),4:ib(m,1)
c 5:ib(m,26),6:ib(m,19),7:ib(m,18),8:ib(m,17)
is(1)=27
is(2)=3
is(3)=2
is(4)=1
is(5)=26
is(6)=19
is(7)=18
is(8)=17
C=0.0
c x=nx face
do 2001 i3=1,3
do 2001 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.3.or.i8.eq.6.or.i8.eq.7) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
2001 continue
do 2000 j=1,ny-1
do 2000 k=1,nz-1
m=nxy*(k-1)+j*nx
do 1900 nn=1,3
do 1900 mm=1,8
sum=0.0
do 1899 m3=1,3
do 1899 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
1899 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1900 continue
2000 continue
c y=ny face
do 2011 i3=1,3
do 2011 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.3.or.i8.eq.4.or.i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
2011 continue
do 2010 i=1,nx-1
do 2010 k=1,nz-1
m=nxy*(k-1)+nx*(ny-1)+i
do 1901 nn=1,3
do 1901 mm=1,8
sum=0.0
do 2099 m3=1,3
do 2099 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2099 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1901 continue
2010 continue
c z=nz face
do 2021 i3=1,3
do 2021 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.5.or.i8.eq.6.or.i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
2021 continue
do 2020 i=1,nx-1
do 2020 j=1,ny-1
m=nxy*(nz-1)+nx*(j-1)+i
do 1902 nn=1,3
do 1902 mm=1,8
sum=0.0
do 2019 m3=1,3
do 2019 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2019 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1902 continue
2020 continue
c x=nx y=ny edge
do 2031 i3=1,3
do 2031 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.6) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.4.or.i8.eq.8) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.3.or.i8.eq.7) then
delta(i8,1)=exy*ny+exx*nx
delta(i8,2)=eyy*ny+exy*nx
delta(i8,3)=eyz*ny+exz*nx
end if
2031 continue
do 2030 k=1,nz-1
m=nxy*k
do 1903 nn=1,3
do 1903 mm=1,8
sum=0.0
do 2029 m3=1,3
do 2029 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2029 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1903 continue
2030 continue
c x=nx z=nz edge
do 2041 i3=1,3
do 2041 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.3) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.5.or.i8.eq.8) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.6.or.i8.eq.7) then
delta(i8,1)=exz*nz+exx*nx
delta(i8,2)=eyz*nz+exy*nx
delta(i8,3)=ezz*nz+exz*nx
end if
2041 continue
do 2040 j=1,ny-1
m=nxy*(nz-1)+nx*(j-1)+nx
do 1904 nn=1,3
do 1904 mm=1,8
sum=0.0
do 2039 m3=1,3
do 2039 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2039 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1904 continue
2040 continue
c y=ny z=nz edge
do 2051 i3=1,3
do 2051 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.5.or.i8.eq.6) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.3.or.i8.eq.4) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exy*ny+exz*nz
delta(i8,2)=eyy*ny+eyz*nz
delta(i8,3)=eyz*ny+ezz*nz
end if
2051 continue
do 2050 i=1,nx-1
m=nxy*(nz-1)+nx*(ny-1)+i
do 1905 nn=1,3
do 1905 mm=1,8
sum=0.0
do 2049 m3=1,3
do 2049 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2049 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1905 continue
2050 continue
c x=nx y=ny z=nz corner
do 2061 i3=1,3
do 2061 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.4) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.5) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.8) then
delta(i8,1)=exy*ny+exz*nz
delta(i8,2)=eyy*ny+eyz*nz
delta(i8,3)=eyz*ny+ezz*nz
end if
if(i8.eq.6) then
delta(i8,1)=exx*nx+exz*nz
delta(i8,2)=exy*nx+eyz*nz
delta(i8,3)=exz*nx+ezz*nz
end if
if(i8.eq.3) then
delta(i8,1)=exx*nx+exy*ny
delta(i8,2)=exy*nx+eyy*ny
delta(i8,3)=exz*nx+eyz*ny
end if
if(i8.eq.7) then
delta(i8,1)=exx*nx+exy*ny+exz*nz
delta(i8,2)=exy*nx+eyy*ny+eyz*nz
delta(i8,3)=exz*nx+eyz*ny+ezz*nz
end if
2061 continue
m=nx*ny*nz
do 1906 nn=1,3
do 1906 mm=1,8
sum=0.0
do 2059 m3=1,3
do 2059 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2059 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1906 continue
return
end
c Subroutine that computes derivatives of the b-vector with respect
c to the macrostrains. Since b is linear in the macrostrains, the
c derivative with respect to any one of them can be computed simply
c by letting that macrostrain, within the subroutine, be equal to one,
c and all the other macrostrains to be zero.
c Very similar to 1221 loop in femat for b.
subroutine bgrad(nx,ny,nz,ns,exx,eyy,ezz,exz,eyz,exy)
real b(8002,3)
real dk(100,8,3,8,3),delta(8,3),zcon(2,3,2,3)
integer is(8)
integer*4 ib(8002,27)
integer*2 pix(8002)
common/list3/ib
common/list4/pix
common/list5/dk,b,C,zcon,Y
nxy=nx*ny
c exx, eyy, ezz, exz, eyz, exy are the artificial macrostrains used
c to get the gradient terms (appropriate combinations of 1's and 0's).
c Set up vector for linear term
do 5000 m3=1,3
do 5000 m=1,ns
b(m,m3)=0.0
5000 continue
is(1)=27
is(2)=3
is(3)=2
is(4)=1
is(5)=26
is(6)=19
is(7)=18
is(8)=17
c x=nx face
do 2001 i3=1,3
do 2001 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.3.or.i8.eq.6.or.i8.eq.7) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
2001 continue
do 2000 j=1,ny-1
do 2000 k=1,nz-1
m=nxy*(k-1)+j*nx
do 1900 nn=1,3
do 1900 mm=1,8
sum=0.0
do 1899 m3=1,3
do 1899 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
1899 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1900 continue
2000 continue
c y=ny face
do 2011 i3=1,3
do 2011 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.3.or.i8.eq.4.or.i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
2011 continue
do 2010 i=1,nx-1
do 2010 k=1,nz-1
m=nxy*(k-1)+nx*(ny-1)+i
do 1901 nn=1,3
do 1901 mm=1,8
sum=0.0
do 2099 m3=1,3
do 2099 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2099 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1901 continue
2010 continue
c z=nz face
do 2021 i3=1,3
do 2021 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.5.or.i8.eq.6.or.i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
2021 continue
do 2020 i=1,nx-1
do 2020 j=1,ny-1
m=nxy*(nz-1)+nx*(j-1)+i
do 1902 nn=1,3
do 1902 mm=1,8
sum=0.0
do 2019 m3=1,3
do 2019 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2019 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1902 continue
2020 continue
c x=nx y=ny edge
do 2031 i3=1,3
do 2031 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.6) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.4.or.i8.eq.8) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.3.or.i8.eq.7) then
delta(i8,1)=exy*ny+exx*nx
delta(i8,2)=eyy*ny+exy*nx
delta(i8,3)=eyz*ny+exz*nx
end if
2031 continue
do 2030 k=1,nz-1
m=nxy*k
do 1903 nn=1,3
do 1903 mm=1,8
sum=0.0
do 2029 m3=1,3
do 2029 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2029 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1903 continue
2030 continue
c x=nx z=nz edge
do 2041 i3=1,3
do 2041 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.3) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.5.or.i8.eq.8) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.6.or.i8.eq.7) then
delta(i8,1)=exz*nz+exx*nx
delta(i8,2)=eyz*nz+exy*nx
delta(i8,3)=ezz*nz+exz*nx
end if
2041 continue
do 2040 j=1,ny-1
m=nxy*(nz-1)+nx*(j-1)+nx
do 1904 nn=1,3
do 1904 mm=1,8
sum=0.0
do 2039 m3=1,3
do 2039 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2039 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1904 continue
2040 continue
c y=ny z=nz edge
do 2051 i3=1,3
do 2051 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.5.or.i8.eq.6) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.3.or.i8.eq.4) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exy*ny+exz*nz
delta(i8,2)=eyy*ny+eyz*nz
delta(i8,3)=eyz*ny+ezz*nz
end if
2051 continue
do 2050 i=1,nx-1
m=nxy*(nz-1)+nx*(ny-1)+i
do 1905 nn=1,3
do 1905 mm=1,8
sum=0.0
do 2049 m3=1,3
do 2049 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2049 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1905 continue
2050 continue
c x=nx y=ny z=nz corner
do 2061 i3=1,3
do 2061 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.4) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.5) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.8) then
delta(i8,1)=exy*ny+exz*nz
delta(i8,2)=eyy*ny+eyz*nz
delta(i8,3)=eyz*ny+ezz*nz
end if
if(i8.eq.6) then
delta(i8,1)=exx*nx+exz*nz
delta(i8,2)=exy*nx+eyz*nz
delta(i8,3)=exz*nx+ezz*nz
end if
if(i8.eq.3) then
delta(i8,1)=exx*nx+exy*ny
delta(i8,2)=exy*nx+eyy*ny
delta(i8,3)=exz*nx+eyz*ny
end if
if(i8.eq.7) then
delta(i8,1)=exx*nx+exy*ny+exz*nz
delta(i8,2)=exy*nx+eyy*ny+eyz*nz
delta(i8,3)=exz*nx+eyz*ny+ezz*nz
end if
2061 continue
m=nx*ny*nz
do 1906 nn=1,3
do 1906 mm=1,8
sum=0.0
do 2059 m3=1,3
do 2059 m8=1,8
sum=sum+delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)
2059 continue
b(ib(m,is(mm)),nn)=b(ib(m,is(mm)),nn)+sum
1906 continue
return
end
c Subroutine computes the quadratic term in the macrostrains, that comes
c from the periodic boundary conditions, and sets it up as a
c (2,3) x (2,3) matrix that couples to the six macrostrains
subroutine const(dk,ns,zcon,nx,ny,nz)
real dk(100,8,3,8,3),zcon(2,3,2,3),delta(8,3)
real pp(6,6),s(6,6)
integer*2 pix(8002)
common/list4/pix
c routine to set up 6 x 6 matrix for energy term involving macro-strains
c only, pulled out of femat
c Idea is to evaluate the quadratic term in the macrostrains repeatedly
c for choices of strain like
c exx=1, exy=1, all others = 0, build up 21 choices, then recombine
c to get matrix elements by themselves
nxy=nx*ny
nss=ns+2
do 1111 i=1,6
do 1111 j=1,6
s(i,j)=0.0
pp(i,j)=0.0
1111 continue
do 5000 ii=1,6
do 5000 jj=ii,6
econ=0.0
exx=0.0
eyy=0.0
ezz=0.0
exz=0.0
eyz=0.0
exy=0.0
if(ii.eq.1.and.jj.eq.1) exx=1.0
if(ii.eq.2.and.jj.eq.2) eyy=1.0
if(ii.eq.3.and.jj.eq.3) ezz=1.0
if(ii.eq.4.and.jj.eq.4) exz=1.0
if(ii.eq.5.and.jj.eq.5) eyz=1.0
if(ii.eq.6.and.jj.eq.6) exy=1.0
if(ii.eq.1.and.jj.eq.2) then
exx=1.0
eyy=1.0
end if
if(ii.eq.1.and.jj.eq.3) then
exx=1.0
ezz=1.0
end if
if(ii.eq.1.and.jj.eq.4) then
exx=1.0
exz=1.0
end if
if(ii.eq.1.and.jj.eq.5) then
exx=1.0
eyz=1.0
end if
if(ii.eq.1.and.jj.eq.6) then
exx=1.0
exy=1.0
end if
if(ii.eq.2.and.jj.eq.3) then
eyy=1.0
ezz=1.0
end if
if(ii.eq.2.and.jj.eq.4) then
eyy=1.0
exz=1.0
end if
if(ii.eq.2.and.jj.eq.5) then
eyy=1.0
eyz=1.0
end if
if(ii.eq.2.and.jj.eq.6) then
eyy=1.0
exy=1.0
end if
if(ii.eq.3.and.jj.eq.4) then
ezz=1.0
exz=1.0
end if
if(ii.eq.3.and.jj.eq.5) then
ezz=1.0
eyz=1.0
end if
if(ii.eq.3.and.jj.eq.6) then
ezz=1.0
exy=1.0
end if
if(ii.eq.4.and.jj.eq.5) then
exz=1.0
eyz=1.0
end if
if(ii.eq.4.and.jj.eq.6) then
exz=1.0
exy=1.0
end if
if(ii.eq.5.and.jj.eq.6) then
eyz=1.0
exy=1.0
end if
c x=nx face
do 2001 i3=1,3
do 2001 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.3.or.i8.eq.6.or.i8.eq.7) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
2001 continue
do 2000 j=1,ny-1
do 2000 k=1,nz-1
m=nxy*(k-1)+j*nx
do 1900 nn=1,3
do 1900 mm=1,8
do 1899 m3=1,3
do 1899 m8=1,8
econ=econ+0.5*delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)*delta(mm,nn)
1899 continue
1900 continue
2000 continue
c y=ny face
do 2011 i3=1,3
do 2011 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.3.or.i8.eq.4.or.i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
2011 continue
do 2010 i=1,nx-1
do 2010 k=1,nz-1
m=nxy*(k-1)+nx*(ny-1)+i
do 1901 nn=1,3
do 1901 mm=1,8
do 2099 m3=1,3
do 2099 m8=1,8
econ=econ+0.5*delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)*delta(mm,nn)
2099 continue
1901 continue
2010 continue
c z=nz face
do 2021 i3=1,3
do 2021 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.5.or.i8.eq.6.or.i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
2021 continue
do 2020 i=1,nx-1
do 2020 j=1,ny-1
m=nxy*(nz-1)+nx*(j-1)+i
do 1902 nn=1,3
do 1902 mm=1,8
do 2019 m3=1,3
do 2019 m8=1,8
econ=econ+0.5*delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)*delta(mm,nn)
2019 continue
1902 continue
2020 continue
c x=nx y=ny edge
do 2031 i3=1,3
do 2031 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.6) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.4.or.i8.eq.8) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.3.or.i8.eq.7) then
delta(i8,1)=exy*ny+exx*nx
delta(i8,2)=eyy*ny+exy*nx
delta(i8,3)=eyz*ny+exz*nx
end if
2031 continue
do 2030 k=1,nz-1
m=nxy*k
do 1903 nn=1,3
do 1903 mm=1,8
do 2029 m3=1,3
do 2029 m8=1,8
econ=econ+0.5*delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)*delta(mm,nn)
2029 continue
1903 continue
2030 continue
c x=nx z=nz edge
do 2041 i3=1,3
do 2041 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2.or.i8.eq.3) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.5.or.i8.eq.8) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.6.or.i8.eq.7) then
delta(i8,1)=exz*nz+exx*nx
delta(i8,2)=eyz*nz+exy*nx
delta(i8,3)=ezz*nz+exz*nx
end if
2041 continue
do 2040 j=1,ny-1
m=nxy*(nz-1)+nx*(j-1)+nx
do 1904 nn=1,3
do 1904 mm=1,8
do 2039 m3=1,3
do 2039 m8=1,8
econ=econ+0.5*delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)*delta(mm,nn)
2039 continue
1904 continue
2040 continue
c y=ny z=nz edge
do 2051 i3=1,3
do 2051 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.5.or.i8.eq.6) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.3.or.i8.eq.4) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.7.or.i8.eq.8) then
delta(i8,1)=exy*ny+exz*nz
delta(i8,2)=eyy*ny+eyz*nz
delta(i8,3)=eyz*ny+ezz*nz
end if
2051 continue
do 2050 i=1,nx-1
m=nxy*(nz-1)+nx*(ny-1)+i
do 1905 nn=1,3
do 1905 mm=1,8
do 2049 m3=1,3
do 2049 m8=1,8
econ=econ+0.5*delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)*delta(mm,nn)
2049 continue
1905 continue
2050 continue
c x=nx y=ny z=nz corner
do 2061 i3=1,3
do 2061 i8=1,8
delta(i8,i3)=0.0
if(i8.eq.2) then
delta(i8,1)=exx*nx
delta(i8,2)=exy*nx
delta(i8,3)=exz*nx
end if
if(i8.eq.4) then
delta(i8,1)=exy*ny
delta(i8,2)=eyy*ny
delta(i8,3)=eyz*ny
end if
if(i8.eq.5) then
delta(i8,1)=exz*nz
delta(i8,2)=eyz*nz
delta(i8,3)=ezz*nz
end if
if(i8.eq.8) then
delta(i8,1)=exy*ny+exz*nz
delta(i8,2)=eyy*ny+eyz*nz
delta(i8,3)=eyz*ny+ezz*nz
end if
if(i8.eq.6) then
delta(i8,1)=exx*nx+exz*nz
delta(i8,2)=exy*nx+eyz*nz
delta(i8,3)=exz*nx+ezz*nz
end if
if(i8.eq.3) then
delta(i8,1)=exx*nx+exy*ny
delta(i8,2)=exy*nx+eyy*ny
delta(i8,3)=exz*nx+eyz*ny
end if
if(i8.eq.7) then
delta(i8,1)=exx*nx+exy*ny+exz*nz
delta(i8,2)=exy*nx+eyy*ny+eyz*nz
delta(i8,3)=exz*nx+eyz*ny+ezz*nz
end if
2061 continue
m=nx*ny*nz
do 1906 nn=1,3
do 1906 mm=1,8
do 2059 m3=1,3
do 2059 m8=1,8
econ=econ+0.5*delta(m8,m3)*dk(pix(m),m8,m3,mm,nn)*delta(mm,nn)
2059 continue
1906 continue
pp(ii,jj)=econ*2.
5000 continue
do 6000 i=1,6
do 6000 j=i,6
if(i.eq.j) s(i,j)=pp(i,j)
if(i.ne.j) then
s(i,j)=pp(i,j)-pp(i,i)-pp(j,j)
end if
6000 continue
do 7000 i=1,6
do 7000 j=1,6
pp(i,j)=0.5*(s(i,j)+s(j,i))
7000 continue
c now map pp(i,j) into zcon(2,3,2,3)
do 7200 i=1,2
do 7200 j=1,2
do 7200 mi=1,3
do 7200 mj=1,3
if(i.eq.1) ii=i+mi-1
if(i.eq.2) ii=i+mi+1
if(j.eq.1) jj=j+mj-1
if(j.eq.2) jj=j+mj+1
zcon(i,mi,j,mj)=pp(ii,jj)
write(7,1) i,mi,j,mj,zcon(i,mi,j,mj)
1 format(4i5,f20.10)
7200 continue
return
end
c Subroutine computes the total energy, utot, and the gradient, gb,
c for the regular displacements as well as for the macrostrains
subroutine energy(nx,ny,nz,ns,utot)
real u(8002,3),gb(8002,3)
real b(8002,3),T(8002,3)
real cmod(100,6,6),pk(6,8,3),eigen(100,6)
real dk(100,8,3,8,3),zcon(2,3,2,3),C
integer*4 ib(8002,27)
integer*2 pix(8002)
common/list3/ib
common/list4/pix
common/list5/dk,b,C,zcon,Y
common/list6/u
common/list7/gb
common/list8/cmod,T,eigen
nss=ns+2
c Do global matrix multiply via small stiffness matrices, Ah = A * h
c The long statement below correctly brings in all the terms from
c the global matrix A using only the small stiffness matrices.
do 2090 m3=1,3
do 2090 m=1,nss
gb(m,m3)=0.0
2090 continue
do 3000 j=1,3
do 3000 n=1,3
do 3000 m=1,ns
gb(m,j)=gb(m,j)+u(ib(m,1),n)*( dk(pix(ib(m,27)),1,j,4,n)
&+dk(pix(ib(m,7)),2,j,3,n)
&+dk(pix(ib(m,25)),5,j,8,n)+dk(pix(ib(m,15)),6,j,7,n) )+
&u(ib(m,2),n)*( dk(pix(ib(m,27)),1,j,3,n)
&+dk(pix(ib(m,25)),5,j,7,n) )+
&u(ib(m,3),n)*( dk(pix(ib(m,27)),1,j,2,n)+dk(pix(ib(m,5)),4,j,3,n)+
&dk(pix(ib(m,13)),8,j,7,n)+dk(pix(ib(m,25)),5,j,6,n) )+
&u(ib(m,4),n)*( dk(pix(ib(m,5)),4,j,2,n)
&+dk(pix(ib(m,13)),8,j,6,n) )+
&u(ib(m,5),n)*( dk(pix(ib(m,6)),3,j,2,n)+dk(pix(ib(m,5)),4,j,1,n)+
&dk(pix(ib(m,14)),7,j,6,n)+dk(pix(ib(m,13)),8,j,5,n) )+
&u(ib(m,6),n)*( dk(pix(ib(m,6)),3,j,1,n)
&+dk(pix(ib(m,14)),7,j,5,n) )+
&u(ib(m,7),n)*( dk(pix(ib(m,6)),3,j,4,n)+dk(pix(ib(m,7)),2,j,1,n)+
&dk(pix(ib(m,14)),7,j,8,n)+dk(pix(ib(m,15)),6,j,5,n) )+
&u(ib(m,8),n)*( dk(pix(ib(m,7)),2,j,4,n)
&+dk(pix(ib(m,15)),6,j,8,n) )+
&u(ib(m,9),n)*( dk(pix(ib(m,25)),5,j,4,n)
&+dk(pix(ib(m,15)),6,j,3,n) )+
&u(ib(m,10),n)*( dk(pix(ib(m,25)),5,j,3,n) )+
&u(ib(m,11),n)*( dk(pix(ib(m,13)),8,j,3,n)
&+dk(pix(ib(m,25)),5,j,2,n) )+
&u(ib(m,12),n)*( dk(pix(ib(m,13)),8,j,2,n) )+
&u(ib(m,13),n)*( dk(pix(ib(m,13)),8,j,1,n)
&+dk(pix(ib(m,14)),7,j,2,n) )+
&u(ib(m,14),n)*( dk(pix(ib(m,14)),7,j,1,n) )+
&u(ib(m,15),n)*( dk(pix(ib(m,14)),7,j,4,n)
&+dk(pix(ib(m,15)),6,j,1,n) )+
&u(ib(m,16),n)*( dk(pix(ib(m,15)),6,j,4,n) )+
&u(ib(m,17),n)*( dk(pix(ib(m,27)),1,j,8,n)
&+dk(pix(ib(m,7)),2,j,7,n) )+
&u(ib(m,18),n)*( dk(pix(ib(m,27)),1,j,7,n) )+
&u(ib(m,19),n)*( dk(pix(ib(m,27)),1,j,6,n)
&+dk(pix(ib(m,5)),4,j,7,n) )+
&u(ib(m,20),n)*( dk(pix(ib(m,5)),4,j,6,n) )+
&u(ib(m,21),n)*( dk(pix(ib(m,5)),4,j,5,n)
&+dk(pix(ib(m,6)),3,j,6,n) )+
&u(ib(m,22),n)*( dk(pix(ib(m,6)),3,j,5,n) )+
&u(ib(m,23),n)*( dk(pix(ib(m,6)),3,j,8,n)
&+dk(pix(ib(m,7)),2,j,5,n) )+
&u(ib(m,24),n)*( dk(pix(ib(m,7)),2,j,8,n) )+
&u(ib(m,25),n)*( dk(pix(ib(m,14)),7,j,3,n)
&+dk(pix(ib(m,13)),8,j,4,n)+
&dk(pix(ib(m,15)),6,j,2,n)+dk(pix(ib(m,25)),5,j,1,n) )+
&u(ib(m,26),n)*( dk(pix(ib(m,6)),3,j,7,n)
&+dk(pix(ib(m,5)),4,j,8,n)+
&dk(pix(ib(m,27)),1,j,5,n)+dk(pix(ib(m,7)),2,j,6,n) )+
&u(ib(m,27),n)*( dk(pix(ib(m,27)),1,j,1,n)
&+dk(pix(ib(m,7)),2,j,2,n)+
&dk(pix(ib(m,6)),3,j,3,n)+dk(pix(ib(m,5)),4,j,4,n)
&+dk(pix(ib(m,25)),5,j,5,n)+
&dk(pix(ib(m,15)),6,j,6,n)+dk(pix(ib(m,14)),7,j,7,n)+
&dk(pix(ib(m,13)),8,j,8,n) )
3000 continue
utot=0.0
gtot=0.0
g2=0.0
g1=0.0
do 3100 m3=1,3
do 3100 m=1,ns
utot=utot+0.5*u(m,m3)*gb(m,m3)+b(m,m3)*u(m,m3)
c this is gradient of energy with respect to normal displacements
gb(m,m3)=gb(m,m3)+b(m,m3)
3100 continue
c compute "constant" macrostrain energy term
C=0.0
do 7200 i=1,2
do 7200 j=1,2
do 7200 mi=1,3
do 7200 mj=1,3
C=C+0.5*u(ns+i,mi)*zcon(i,mi,j,mj)*u(ns+j,mj)
7200 continue
utot=utot+C
c now add in constant term from thermal energy, Y
utot=utot+Y
c now add in linear term in thermal energy
do 7171 m3=1,3
do 7171 m=1,nss
utot=utot+T(m,m3)*u(m,m3)
7171 continue
c now compute gradient with respect to macrostrains
c put in piece from first derivative of zcon quadratic term
do 7300 i=1,2
do 7300 mi=1,3
sum=0.0
do 7250 j=1,2
do 7250 mj=1,3
sum=sum+zcon(i,mi,j,mj)*u(ns+j,mj)
7250 continue
gb(ns+i,mi)=sum
7300 continue
c add in piece of gradient, for displacements as well as macrostrains,
c that come from linear term in thermal energy
do 3150 m3=1,3
do 3150 m=1,nss
gb(m,m3)=gb(m,m3)+T(m,m3)
3150 continue
c now generate part that comes from b . u term
c do by calling b generation with appropriate macrostrain set to 1 to
c get that partial derivative, just use bgrad (taken from femat),
c skip dk and zcon part
do 8100 ii=1,6
exx=0.0
eyy=0.0
ezz=0.0
exz=0.0
eyz=0.0
exy=0.0
if(ii.eq.1) exx=1.0
if(ii.eq.2) eyy=1.0
if(ii.eq.3) ezz=1.0
if(ii.eq.4) exz=1.0
if(ii.eq.5) eyz=1.0
if(ii.eq.6) exy=1.0
call bgrad(nx,ny,nz,ns,exx,eyy,ezz,exz,eyz,exy)
sum=0.0
do 8200 m3=1,3
do 8200 m=1,ns
sum=sum+u(m,m3)*b(m,m3)
8200 continue
if(ii.eq.1) gb(ns+1,1)=gb(ns+1,1)+sum
if(ii.eq.2) gb(ns+1,2)=gb(ns+1,2)+sum
if(ii.eq.3) gb(ns+1,3)=gb(ns+1,3)+sum
if(ii.eq.4) gb(ns+2,1)=gb(ns+2,1)+sum
if(ii.eq.5) gb(ns+2,2)=gb(ns+2,2)+sum
if(ii.eq.6) gb(ns+2,3)=gb(ns+2,3)+sum
8100 continue
return
end
c Subroutine that carries out the conjugate gradient relaxation process
subroutine dembx(nx,ny,nz,ns,Lstep,gg,gtest,ldemb,kkk)
real u(8002,3),gb(8002,3),b(8002,3)
real h(8002,3),Ah(8002,3)
real dk(100,8,3,8,3),zcon(2,3,2,3)
real lambda,gamma
integer*4 ib(8002,27)
integer*2 pix(8002)
common/list2/h,Ah
common/list3/ib
common/list4/pix
common/list5/dk,b,C,zcon,Y
common/list6/u
common/list7/gb
nss=ns+2
c Initialize the conjugate direction vector on first call to dembx only.
c For calls to dembx after the first, we want to continue using the value
c of h determined in the previous call. Of course, if npoints
c is greater than 1, then this initialization step will be run each time
c a new microstructure is used, as kkk will be reset to 1 every time
c the counter micro is increased.
if(kkk.eq.1) then
do 500 m3=1,3
do 500 m=1,nss
h(m,m3)=gb(m,m3)
500 continue
end if
c Lstep counts the number of conjugate gradient steps taken
c in each call to dembx
Lstep=0
do 800 ijk=1,ldemb
Lstep=Lstep+1
do 290 m3=1,3
do 290 m=1,nss
Ah(m,m3)=0.0
290 continue
c Do global matrix multiply via small stiffness matrices, Ah = A * h
c The long statement below correctly brings in all the terms from
c the global matrix A using only the small stiffness matrices.
do 400 j=1,3
do 400 n=1,3
do 400 m=1,ns
Ah(m,j)=Ah(m,j)+h(ib(m,1),n)*( dk(pix(ib(m,27)),1,j,4,n)
&+dk(pix(ib(m,7)),2,j,3,n)
&+dk(pix(ib(m,25)),5,j,8,n)+dk(pix(ib(m,15)),6,j,7,n) )+
&h(ib(m,2),n)*( dk(pix(ib(m,27)),1,j,3,n)
&+dk(pix(ib(m,25)),5,j,7,n) )+
&h(ib(m,3),n)*( dk(pix(ib(m,27)),1,j,2,n)+dk(pix(ib(m,5)),4,j,3,n)+
&dk(pix(ib(m,13)),8,j,7,n)+dk(pix(ib(m,25)),5,j,6,n) )+
&h(ib(m,4),n)*( dk(pix(ib(m,5)),4,j,2,n)
&+dk(pix(ib(m,13)),8,j,6,n) )+
&h(ib(m,5),n)*( dk(pix(ib(m,6)),3,j,2,n)+dk(pix(ib(m,5)),4,j,1,n)+
&dk(pix(ib(m,14)),7,j,6,n)+dk(pix(ib(m,13)),8,j,5,n) )+
&h(ib(m,6),n)*( dk(pix(ib(m,6)),3,j,1,n)
&+dk(pix(ib(m,14)),7,j,5,n) )+
&h(ib(m,7),n)*( dk(pix(ib(m,6)),3,j,4,n)+dk(pix(ib(m,7)),2,j,1,n)+
&dk(pix(ib(m,14)),7,j,8,n)+dk(pix(ib(m,15)),6,j,5,n) )+
&h(ib(m,8),n)*( dk(pix(ib(m,7)),2,j,4,n)
&+dk(pix(ib(m,15)),6,j,8,n) )+
&h(ib(m,9),n)*( dk(pix(ib(m,25)),5,j,4,n)
&+dk(pix(ib(m,15)),6,j,3,n) )+
&h(ib(m,10),n)*( dk(pix(ib(m,25)),5,j,3,n) )+
&h(ib(m,11),n)*( dk(pix(ib(m,13)),8,j,3,n)
&+dk(pix(ib(m,25)),5,j,2,n) )+
&h(ib(m,12),n)*( dk(pix(ib(m,13)),8,j,2,n) )+
&h(ib(m,13),n)*( dk(pix(ib(m,13)),8,j,1,n)
&+dk(pix(ib(m,14)),7,j,2,n) )+
&h(ib(m,14),n)*( dk(pix(ib(m,14)),7,j,1,n) )+
&h(ib(m,15),n)*( dk(pix(ib(m,14)),7,j,4,n)
&+dk(pix(ib(m,15)),6,j,1,n) )+
&h(ib(m,16),n)*( dk(pix(ib(m,15)),6,j,4,n) )+
&h(ib(m,17),n)*( dk(pix(ib(m,27)),1,j,8,n)
&+dk(pix(ib(m,7)),2,j,7,n) )+
&h(ib(m,18),n)*( dk(pix(ib(m,27)),1,j,7,n) )+
&h(ib(m,19),n)*( dk(pix(ib(m,27)),1,j,6,n)
&+dk(pix(ib(m,5)),4,j,7,n) )+
&h(ib(m,20),n)*( dk(pix(ib(m,5)),4,j,6,n) )+
&h(ib(m,21),n)*( dk(pix(ib(m,5)),4,j,5,n)
&+dk(pix(ib(m,6)),3,j,6,n) )+
&h(ib(m,22),n)*( dk(pix(ib(m,6)),3,j,5,n) )+
&h(ib(m,23),n)*( dk(pix(ib(m,6)),3,j,8,n)
&+dk(pix(ib(m,7)),2,j,5,n) )+
&h(ib(m,24),n)*( dk(pix(ib(m,7)),2,j,8,n) )+
&h(ib(m,25),n)*( dk(pix(ib(m,14)),7,j,3,n)
&+dk(pix(ib(m,13)),8,j,4,n)+
&dk(pix(ib(m,15)),6,j,2,n)+dk(pix(ib(m,25)),5,j,1,n) )+
&h(ib(m,26),n)*( dk(pix(ib(m,6)),3,j,7,n)
&+dk(pix(ib(m,5)),4,j,8,n)+
&dk(pix(ib(m,27)),1,j,5,n)+dk(pix(ib(m,7)),2,j,6,n) )+
&h(ib(m,27),n)*( dk(pix(ib(m,27)),1,j,1,n)
&+dk(pix(ib(m,7)),2,j,2,n)+
&dk(pix(ib(m,6)),3,j,3,n)+dk(pix(ib(m,5)),4,j,4,n)
&+dk(pix(ib(m,25)),5,j,5,n)+
&dk(pix(ib(m,15)),6,j,6,n)+dk(pix(ib(m,14)),7,j,7,n)+
&dk(pix(ib(m,13)),8,j,8,n) )
400 continue
c The above accurately gives the second derivative matrix with respect
c to nodal displacements, but fails to give the 2nd derivative terms that
c include the macrostrains [ du d(strain) and d(strain)d(strain) ].
c Use repeated calls to bgrad to generate mixed 2nd derivatives terms,
c plus use zcon in order to correct the matrix multiply and correctly bring
c in macrostrain terms (see manual, Sec. 2.4).
do 8100 ii=1,6
e11=0.0
e22=0.0
e33=0.0
e13=0.0
e23=0.0
e12=0.0
if(ii.eq.1) e11=1.0
if(ii.eq.2) e22=1.0
if(ii.eq.3) e33=1.0
if(ii.eq.4) e13=1.0
if(ii.eq.5) e23=1.0
if(ii.eq.6) e12=1.0
call bgrad(nx,ny,nz,ns,e11,e22,e33,e13,e23,e12)
c now fill in terms from matrix multiply
c right hand sides, 1 to ns
do 3333 m=1,ns
do 3333 m1=1,3
if(ii.eq.1) Ah(m,m1)=Ah(m,m1)+b(m,m1)*h(ns+1,1)
if(ii.eq.2) Ah(m,m1)=Ah(m,m1)+b(m,m1)*h(ns+1,2)
if(ii.eq.3) Ah(m,m1)=Ah(m,m1)+b(m,m1)*h(ns+1,3)
if(ii.eq.4) Ah(m,m1)=Ah(m,m1)+b(m,m1)*h(nss,1)
if(ii.eq.5) Ah(m,m1)=Ah(m,m1)+b(m,m1)*h(nss,2)
if(ii.eq.6) Ah(m,m1)=Ah(m,m1)+b(m,m1)*h(nss,3)
3333 continue
c now do across bottom, 1 to ns
do 3334 m=1,ns
if(ii.eq.1) Ah(ns+1,1)=Ah(ns+1,1)+b(m,1)*h(m,1)+
+b(m,2)*h(m,2)+b(m,3)*h(m,3)
if(ii.eq.2) Ah(ns+1,2)=Ah(ns+1,2)+b(m,1)*h(m,1)+
+b(m,2)*h(m,2)+b(m,3)*h(m,3)
if(ii.eq.3) Ah(ns+1,3)=Ah(ns+1,3)+b(m,1)*h(m,1)+
+b(m,2)*h(m,2)+b(m,3)*h(m,3)
if(ii.eq.4) Ah(nss,1)=Ah(nss,1)+b(m,1)*h(m,1)+
+b(m,2)*h(m,2)+b(m,3)*h(m,3)
if(ii.eq.5) Ah(nss,2)=Ah(nss,2)+b(m,1)*h(m,1)+
+b(m,2)*h(m,2)+b(m,3)*h(m,3)
if(ii.eq.6) Ah(nss,3)=Ah(nss,3)+b(m,1)*h(m,1)+
+b(m,2)*h(m,2)+b(m,3)*h(m,3)
3334 continue
c now do righthand corner terms, ns+1 to nss
do 3335 m=1,2
do 3335 m1=1,3
if(ii.eq.1) Ah(ns+1,1)=Ah(ns+1,1)+zcon(1,1,m,m1)*h(ns+m,m1)
if(ii.eq.2) Ah(ns+1,2)=Ah(ns+1,2)+zcon(1,2,m,m1)*h(ns+m,m1)
if(ii.eq.3) Ah(ns+1,3)=Ah(ns+1,3)+zcon(1,3,m,m1)*h(ns+m,m1)
if(ii.eq.4) Ah(nss,1)=Ah(nss,1)+zcon(2,1,m,m1)*h(ns+m,m1)
if(ii.eq.5) Ah(nss,2)=Ah(nss,2)+zcon(2,2,m,m1)*h(ns+m,m1)
if(ii.eq.6) Ah(nss,3)=Ah(nss,3)+zcon(2,3,m,m1)*h(ns+m,m1)
3335 continue
8100 continue
hAh=0.0
do 530 m3=1,3
do 530 m=1,nss
hAh=hAh+h(m,m3)*Ah(m,m3)
530 continue
lambda=gg/hAh
do 540 m3=1,3
do 540 m=1,nss
u(m,m3)=u(m,m3)-lambda*h(m,m3)
gb(m,m3)=gb(m,m3)-lambda*Ah(m,m3)
540 continue
gglast=gg
gg=0.0
do 550 m3=1,3
do 550 m=1,nss
gg=gg+gb(m,m3)*gb(m,m3)
550 continue
if(gg.lt.gtest) goto 1000
gamma=gg/gglast
do 570 m3=1,3
do 570 m=1,nss
h(m,m3)=gb(m,m3)+gamma*h(m,m3)
570 continue
800 continue
1000 continue
return
end
c Subroutine that computes the six average stresses and six average strains
subroutine stress(nx,ny,nz,ns)
real u(8002,3),uu(8,3)
real T(8002,3),eigen(100,6)
real dndx(8),dndy(8),dndz(8),es(6,8,3),cmod(100,6,6)
integer*4 ib(8002,27)
integer*2 pix(8002)
common/list1/strxx,stryy,strzz,strxz,stryz,strxy
common/list3/ib
common/list4/pix
common/list6/u
common/list8/cmod,T,eigen
common/list11/sxx,syy,szz,sxz,syz,sxy
nxy=nx*ny
nss=ns+2
exx=u(ns+1,1)
eyy=u(ns+1,2)
ezz=u(ns+1,3)
exz=u(nss,1)
eyz=u(nss,2)
exy=u(nss,3)
c set up single pixel strain matrix
dndx(1)=-0.25
dndx(2)=0.25
dndx(3)=0.25
dndx(4)=-0.25
dndx(5)=-0.25
dndx(6)=0.25
dndx(7)=0.25
dndx(8)=-0.25
dndy(1)=-0.25
dndy(2)=-0.25
dndy(3)=0.25
dndy(4)=0.25
dndy(5)=-0.25
dndy(6)=-0.25
dndy(7)=0.25
dndy(8)=0.25
dndz(1)=-0.25
dndz(2)=-0.25
dndz(3)=-0.25
dndz(4)=-0.25
dndz(5)=0.25
dndz(6)=0.25
dndz(7)=0.25
dndz(8)=0.25
c Build average strain matrix, follows code in femat, but for average
c strain over a pixel, not the strain at a point
do 2799 n1=1,6
do 2799 n2=1,8
do 2799 n3=1,3
es(n1,n2,n3)=0.0
2799 continue
do 2797 n=1,8
es(1,n,1)=dndx(n)
es(2,n,2)=dndy(n)
es(3,n,3)=dndz(n)
es(4,n,1)=dndz(n)
es(4,n,3)=dndx(n)
es(5,n,2)=dndz(n)
es(5,n,3)=dndy(n)
es(6,n,1)=dndy(n)
es(6,n,2)=dndx(n)
2797 continue
c now compute average stresses and strains in each pixel
sxx=0.0
syy=0.0
szz=0.0
sxz=0.0
syz=0.0
sxy=0.0
strxx=0.0
stryy=0.0
strzz=0.0
strxz=0.0
stryz=0.0
strxy=0.0
do 470 k=1,nz
do 470 j=1,ny
do 470 i=1,nx
m=(k-1)*nxy+(j-1)*nx+i
c load in elements of 8-vector using pd. bd. conds.
do 9898 mm=1,3
uu(1,mm)=u(m,mm)
uu(2,mm)=u(ib(m,3),mm)
uu(3,mm)=u(ib(m,2),mm)
uu(4,mm)=u(ib(m,1),mm)
uu(5,mm)=u(ib(m,26),mm)
uu(6,mm)=u(ib(m,19),mm)
uu(7,mm)=u(ib(m,18),mm)
uu(8,mm)=u(ib(m,17),mm)
9898 continue
c Correct for periodic boundary conditions, some displacements are wrong
c for a pixel on a periodic boundary. Since they come from an opposite
c face, need to put in applied strain to correct them.
if(i.eq.nx) then
uu(2,1)=uu(2,1)+exx*nx
uu(2,2)=uu(2,2)+exy*nx
uu(2,3)=uu(2,3)+exz*nx
uu(3,1)=uu(3,1)+exx*nx
uu(3,2)=uu(3,2)+exy*nx
uu(3,3)=uu(3,3)+exz*nx
uu(6,1)=uu(6,1)+exx*nx
uu(6,2)=uu(6,2)+exy*nx
uu(6,3)=uu(6,3)+exz*nx
uu(7,1)=uu(7,1)+exx*nx
uu(7,2)=uu(7,2)+exy*nx
uu(7,3)=uu(7,3)+exz*nx
end if
if(j.eq.ny) then
uu(3,1)=uu(3,1)+exy*ny
uu(3,2)=uu(3,2)+eyy*ny
uu(3,3)=uu(3,3)+eyz*ny
uu(4,1)=uu(4,1)+exy*ny
uu(4,2)=uu(4,2)+eyy*ny
uu(4,3)=uu(4,3)+eyz*ny
uu(7,1)=uu(7,1)+exy*ny
uu(7,2)=uu(7,2)+eyy*ny
uu(7,3)=uu(7,3)+eyz*ny
uu(8,1)=uu(8,1)+exy*ny
uu(8,2)=uu(8,2)+eyy*ny
uu(8,3)=uu(8,3)+eyz*ny
end if
if(k.eq.nz) then
uu(5,1)=uu(5,1)+exz*nz
uu(5,2)=uu(5,2)+eyz*nz
uu(5,3)=uu(5,3)+ezz*nz
uu(6,1)=uu(6,1)+exz*nz
uu(6,2)=uu(6,2)+eyz*nz
uu(6,3)=uu(6,3)+ezz*nz
uu(7,1)=uu(7,1)+exz*nz
uu(7,2)=uu(7,2)+eyz*nz
uu(7,3)=uu(7,3)+ezz*nz
uu(8,1)=uu(8,1)+exz*nz
uu(8,2)=uu(8,2)+eyz*nz
uu(8,3)=uu(8,3)+ezz*nz
end if
c stresses and strains in a pixel
str11=0.0
str22=0.0
str33=0.0
str13=0.0
str23=0.0
str12=0.0
s11=0.0
s22=0.0
s33=0.0
s13=0.0
s23=0.0
s12=0.0
c********compute average stress and strain tensor in each pixel*************
c First put thermal strain-induced stresses into stress tensor
do 465 n=1,6
str11=str11-cmod(pix(m),1,n)*eigen(pix(m),n)
str22=str22-cmod(pix(m),2,n)*eigen(pix(m),n)
str33=str33-cmod(pix(m),3,n)*eigen(pix(m),n)
str13=str13-cmod(pix(m),4,n)*eigen(pix(m),n)
str23=str23-cmod(pix(m),5,n)*eigen(pix(m),n)
str12=str12-cmod(pix(m),6,n)*eigen(pix(m),n)
465 continue
do 466 n3=1,3
do 466 n8=1,8
c compute non-thermal strains in each pixel
s11=s11+es(1,n8,n3)*uu(n8,n3)
s22=s22+es(2,n8,n3)*uu(n8,n3)
s33=s33+es(3,n8,n3)*uu(n8,n3)
s13=s13+es(4,n8,n3)*uu(n8,n3)
s23=s23+es(5,n8,n3)*uu(n8,n3)
s12=s12+es(6,n8,n3)*uu(n8,n3)
do 466 n=1,6
c compute stresses in each pixel that include both non-thermal
c and thermal strains
str11=str11+cmod(pix(m),1,n)*es(n,n8,n3)*uu(n8,n3)
str22=str22+cmod(pix(m),2,n)*es(n,n8,n3)*uu(n8,n3)
str33=str33+cmod(pix(m),3,n)*es(n,n8,n3)*uu(n8,n3)
str13=str13+cmod(pix(m),4,n)*es(n,n8,n3)*uu(n8,n3)
str23=str23+cmod(pix(m),5,n)*es(n,n8,n3)*uu(n8,n3)
str12=str12+cmod(pix(m),6,n)*es(n,n8,n3)*uu(n8,n3)
466 continue
c Sum local stresses and strains into global stresses and strains
strxx=strxx+str11
stryy=stryy+str22
strzz=strzz+str33
strxz=strxz+str13
stryz=stryz+str23
strxy=strxy+str12
sxx=sxx+s11
syy=syy+s22
szz=szz+s33
sxz=sxz+s13
syz=syz+s23
sxy=sxy+s12
470 continue
c Volume average global stresses and strains
strxx=strxx/float(ns)
stryy=stryy/float(ns)
strzz=strzz/float(ns)
strxz=strxz/float(ns)
stryz=stryz/float(ns)
strxy=strxy/float(ns)
sxx=sxx/float(ns)
syy=syy/float(ns)
szz=szz/float(ns)
sxz=sxz/float(ns)
syz=syz/float(ns)
sxy=sxy/float(ns)
return
end
c Subroutine to count volume fractions of various phases
subroutine assig(ns,nphase,prob)
integer*2 pix(8002)
real prob(100)
common/list4/pix
do 999 i=1,nphase
prob(i)=0.0
999 continue
do 1000 m=1,ns
do 1000 i=1,nphase
if(pix(m).eq.i) then
prob(i)=prob(i)+1
end if
1000 continue
do 998 i=1,nphase
prob(i)=prob(i)/float(ns)
998 continue
return
end
c Subroutine to set up image of microstructure
subroutine ppixel(nx,ny,nz,ns,nphase)
integer*2 pix(8002)
integer*4 ib(8002,27)
common/list3/ib
common/list4/pix
c (USER) If you want to set up a test image inside the program, instead
c of reading it in from a file, this should be done inside this subroutine.
nxy=nx*ny
do 200 k=1,nz
do 200 j=1,ny
do 200 i=1,nx
m=nxy*(k-1)+nx*(j-1)+i
read(9,*) pix(m)
200 continue
c Check for wrong phase labels--less than 1 or greater than nphase
do 500 m=1,ns
if(pix(m).lt.1) then
write(7,*) 'Phase label in pix < 1--error at ',m
end if
if(pix(m).gt.nphase) then
write(7,*) 'Phase label in pix > nphase--error at ',m
end if
500 continue
return
end