c----New grid generator with iz being the slowest varying index--- c c----Formulation with center prescription for displacement---- use portlib parameter (nx=6,ny=6,nz=3,nn=nx*ny*nz,nd=3*nn) real umat(nd,nd),rhs(nd,6) real xx(nn),yy(nn),zz(nn) real uu(3,3),imped(10) integer ipiv(nd) data imped/10*0./ open(1,file='imped.dat',status='unknown') itime1=time() c-----set Gaussian Quadrature weight for later integration---- call weights delta=2./nx c----Generate grid, obtain xx, yy, zz. The coordinate of the center c-------point of the elemental volumes.---- ip=0 do 100 iz=1,nz z=(2*iz-1.)*delta/2. do 100 iy=1,ny y=-1.+(2*iy-1.)*delta/2. do 100 ix=1,nx x=-1.+(2*ix-1.)*delta/2. ip=ip+1 xx(ip)=x yy(ip)=y zz(ip)=z 100 continue if(ip.ne.nn)write (*,*) 'Grid generation error.' c----setup up [U], assume the elemental volumes are cubical---- c----The half width is 1---- do 200 i=1,nn is=3*(i-1) xp=xx(i) yp=yy(i) zp=zz(i) c----set right hand side---- do iq1=1,3 do iq2=1,6 rhs(is+iq1,iq2)=0.0 enddo enddo rhs(is+1,1)=1. rhs(is+2,2)=1. rhs(is+3,3)=1. rhs(is+2,4)=-zp rhs(is+3,4)=yp rhs(is+1,5)=zp rhs(is+3,5)=-xp rhs(is+1,6)=-yp rhs(is+2,6)=xp do 200 j=1,nn io=3*(j-1) xo=xx(j) yo=yy(j) zo=zz(j) if(i.eq.j)then x=0. y=0. z=zo c=zp call intdiag(x,y,z,c,delta,uu) else x=xo-xp y=yo-yp z=zo c=zp if(x.eq.0.0 .and. y.eq.0.0)then call intuu4(x,y,z,c,delta,uu) else call intuu(x,y,z,c,delta,uu) endif endif do ii=1,3 do jj=1,3 umat(is+ii,io+jj)=uu(ii,jj) enddo enddo 200 continue call sgesv(nd,6,umat,nd,ipiv,rhs,nd,info) write (*,*) 'Exit condition:',info itime2=time()-itime1 iq1=itime2/60 iq2=mod(itime2,60) write (*,*) iq1,' min',iq2, 'sec' dv=delta**3 do i=1,nn is=3*(i-1) xp=xx(i) yp=yy(i) zp=zz(i) imped(1)=imped(1)+rhs(is+1,1)*dv imped(2)=imped(2)+rhs(is+2,2)*dv imped(3)=imped(3)+rhs(is+3,3)*dv imped(4)=imped(4)+(-zp*rhs(is+2,4)+yp*rhs(is+3,4))*dv imped(5)=imped(5)+(zp*rhs(is+1,5)-xp*rhs(is+3,5))*dv imped(6)=imped(6)+(-yp*rhs(is+1,6)+xp*rhs(is+3,6))*dv imped(7)=imped(7)+rhs(is+1,5)*dv imped(8)=imped(8)+rhs(is+2,4)*dv imped(9)=imped(9)+(-zp*rhs(is+2,2)+yp*rhs(is+3,2))*dv imped(10)=imped(10)+(zp*rhs(is+1,1)-xp*rhs(is+3,1))*dv enddo write (*,'(5f10.4)') imped write (1,'(i3,i3,i3)') nx,ny,nz write (1,'(5f10.4)') imped end subroutine intdiag(x,y,z,c,delta,uu) real uu(3,3),uh(3,3),ui(3,3) common/gauss4/w(4),e(4) a=delta/2. a3=a**3 do 1 i=1,3 do 1 j=1,3 1 uu(i,j)=0. do 100 ix=1,4 do 100 iy=1,4 do 100 iz=1,4 xx=a*e(ix)+x yy=a*e(iy)+y zz=a*e(iz)+z call getuu(xx,yy,zz,c,uh) call aehlove(xx,yy,zz-c,ui) do 2 i=1,3 do 2 j=1,3 2 uu(i,j)=uu(i,j)+a3*w(ix)*w(iy)*w(iz)*(uh(i,j)-ui(i,j)) 100 continue call istatic(x,y,z-c,delta,ui) do 3 i=1,3 do 3 j=1,3 3 uu(i,j)=uu(i,j)+ui(i,j) return end subroutine intuu(x,y,z,c,delta,uu) real uu(3,3),u(3,3) common/gauss3/w(3),e(3) a=delta/2. a3=a**3 do 1 i=1,3 do 1 j=1,3 1 uu(i,j)=0. do 100 ix=1,3 do 100 iy=1,3 do 100 iz=1,3 xx=a*e(ix)+x yy=a*e(iy)+y zz=a*e(iz)+z call getuu(xx,yy,zz,c,u) do 2 i=1,3 do 2 j=1,3 2 uu(i,j)=uu(i,j)+a3*w(ix)*w(iy)*w(iz)*u(i,j) 100 continue return end subroutine intuu4(x,y,z,c,delta,uu) real uu(3,3),u(3,3) common/gauss4/w(4),e(4) a=delta/2. a3=a**3 do 1 i=1,3 do 1 j=1,3 1 uu(i,j)=0. do 100 ix=1,4 do 100 iy=1,4 do 100 iz=1,4 xx=a*e(ix)+x yy=a*e(iy)+y zz=a*e(iz)+z call getuu(xx,yy,zz,c,u) do 2 i=1,3 do 2 j=1,3 2 uu(i,j)=uu(i,j)+a3*w(ix)*w(iy)*w(iz)*u(i,j) 100 continue return end subroutine getuu(x,y,z,c,uu) real gf(5),uu(3,3) c-----The order of the Green's Functions: rz, zz, rr, psi-r, zr----- x2=x*x y2=y*y r2=x2+y2 r=sqrt(r2) r3=r2*r call mindlin(r,z,c,gf) uu(1,1)=(gf(3)*x2-gf(4)*y2)/r3 uu(1,2)=(gf(3)+gf(4))*x*y/r3 uu(1,3)=gf(5)*x/r2 uu(2,1)=(gf(3)+gf(4))*x*y/r3 uu(2,2)=(gf(3)*y2-gf(4)*x2)/r3 uu(2,3)=gf(5)*y/r2 uu(3,1)=gf(1)*x/r2 uu(3,2)=gf(1)*y/r2 uu(3,3)=gf(2)/r return end subroutine mindlin(r,z,c,gf) real gf(5) real mu pi=3.1415926535 theta=0. c----mu is the Poisson's Ratio---- c----gf is to be multiplied by 1/(G*r)---- mu=1./3. c1=(1.-mu) c2=(1.-2.*mu) c3=(3.-4.*mu) r1=sqrt(r**2+(z-c)**2) r2=sqrt(r**2+(z+c)**2) t1=(z-c)/r1**3 t2=c3*(z-c)/r2**3 t3=-4.*c1*c2/(r2*(r2+z+c)) t4=6.*c*z*(z+c)/r2**5 u=r*(t1+t2+t3+t4)/(16.*pi*c1) t1=c3/r1 t2=(8.*c1*c1-c3)/r2 t3=(z-c)**2/r1**3 t4=(c3*(z+c)**2-2.*c*z)/r2**3 t5=(6.*c*z*(z+c)**2)/r2**5 w=(t1+t2+t3+t4+t5)/(16.*pi*c1) x=r*cos(theta) t1=c3/r1+1./r2+x**2/r1**3+c3*x**2/r2**3 t2=(2.*c*z/r2**3)*(1.-3.*x**2/r2**2) t3=(4.*c1*c2/(r2+z+c))*(1.-x**2/(r2*(r2+z+c))) uh=(t1+t2+t3)/(16.*pi*c1) t1=(z-c)/r1**3+c3*(z-c)/r2**3 t2=-6.*c*z*(z+c)/r2**5 t3=4.*c1*c2/(r2*(r2+z+c)) wh=(t1+t2+t3)*x/(16.*pi*c1) t1=1./r1**3 + c3/r2**3 - 6.*c*z/r2**5 t2=-4.*c1*c2/(r2*(r2+z+c)**2) vh=(t1+t2)*r**2/(16.*pi*c1) gf(1)=u*r gf(2)=w*r gf(3)=uh*r gf(4)=(vh-uh)*r gf(5)=wh*r return end subroutine istatic(x,y,z,delta,ww) real ww(3,3) a=delta/2. pi=3.1415926535 gamma=0.5 g2=gamma**2 c1=(1+g2)/(8.*pi) c2=(1-g2)/(8.*pi) ww(1,1)=c1*s0(x,y,z,a)+c2*sxx(x,y,z,a) ww(1,2)=c2*sxy(x,y,z,a) ww(1,3)=c2*sxz(x,y,z,a) ww(2,1)=c2*sxy(x,y,z,a) ww(2,2)=c1*s0(x,y,z,a)+c2*syy(x,y,z,a) ww(2,3)=c2*syz(x,y,z,a) ww(3,1)=c2*sxz(x,y,z,a) ww(3,2)=c2*syz(x,y,z,a) ww(3,3)=c1*s0(x,y,z,a)+c2*szz(x,y,z,a) return end subroutine aehlove(x,y,z,uu) real uu(3,3),u3(3) call vlove(y,z,x,u3) uu(1,1)=u3(3) uu(1,2)=u3(1) uu(1,3)=u3(2) call vlove(z,x,y,u3) uu(2,1)=u3(2) uu(2,2)=u3(3) uu(2,3)=u3(1) call vlove(x,y,z,u3) uu(3,1)=u3(1) uu(3,2)=u3(2) uu(3,3)=u3(3) return end subroutine vlove(x,y,z,u) real u(3) real gamma pi=3.1415926535 gamma=0.5 g2=gamma**2 c1=(1+g2)/(8.*pi) c2=(1-g2)/(8.*pi) r2=x*x+y*y+z*z r=sqrt(r2) r3=r*r2 u(1)=c2*z*x/r3 u(2)=c2*z*y/r3 u(3)=c1/r+c2*z*z/r3 return end subroutine weights common/gauss2/w2(2),e2(2) common/gauss3/w3(3),e3(3) common/gauss4/w4(4),e4(4) w2(1)=1. w2(2)=1. e2(1)=-1./sqrt(3.) e2(2)=-e2(1) w3(1)=5./9. w3(2)=8./9. w3(3)=w3(1) e3(1)=-sqrt(0.6) e3(2)=0. e3(3)=-e3(1) w4(1)=0.3478548 w4(2)=0.6521432 w4(3)=w4(2) w4(4)=w4(1) e4(1)=-0.8611363 e4(2)=-0.3399810 e4(3)=-e4(2) e4(4)=-e4(1) return end real function s0(x0,y0,z0,a) real x(2),y(2),z(2) x(1)=-a-x0 x(2)=a-x0 y(1)=-a-y0 y(2)=a-y0 z(1)=-a-z0 z(2)=a-z0 sum1=0. do 100 i=1,2 do 100 j=1,2 do 100 k=1,2 s=(-1.)**(i+j+k) sum1=sum1+s*( q1(y(j),z(k),x(i))+q1(x(i),z(k),y(j)) # +q1(x(i),y(j),z(k)) ) 100 continue sum2=0. do 200 i=1,2 do 200 j=1,2 s=(-1.)**(i+j) sum2=sum2+s*( q2(x(i),y(j),z(1),z(2))+q2(y(i),z(j),x(1),x(2)) # +q2(x(i),z(j),y(1),y(2)) ) 200 continue s0=sum1+sum2 return end real function szz(x0,y0,z0,a) real x(2),y(2),z(2) x(1)=-a-x0 x(2)=a-x0 y(1)=-a-y0 y(2)=a-y0 z(1)=-a-z0 z(2)=a-z0 sum1=0. do 100 i=1,2 do 100 j=1,2 do 100 k=1,2 s=(-1.)**(i+j+k) sum1=sum1+s*( q1(y(j),z(k),x(i))+q1(x(i),z(k),y(j)) # -q1(x(i),y(j),z(k)) ) 100 continue sum2=0. do 200 i=1,2 do 200 j=1,2 s=(-1.)**(i+j) sum2=sum2+s*q2(x(i),y(j),z(1),z(2)) 200 continue szz=sum1+sum2 return end real function sxz(x0,y0,z0,a) real x(2),y(2),z(2) x(1)=-a-x0 x(2)=a-x0 y(1)=-a-y0 y(2)=a-y0 z(1)=-a-z0 z(2)=a-z0 sum1=0. do 100 i=1,2 do 100 j=1,2 do 100 k=1,2 s=(-1.)**(i+j+k) t1=x(i)**2+z(k)**2 t2=sqrt(y(j)**2+t1) sum1=sum1-0.5*s*( y(j)*t2+t1*alog(y(j)+t2) ) 100 continue sxz=sum1 return end real function syz(x0,y0,z0,a) real x(2),y(2),z(2) x(1)=-a-x0 x(2)=a-x0 y(1)=-a-y0 y(2)=a-y0 z(1)=-a-z0 z(2)=a-z0 sum1=0. do 100 i=1,2 do 100 j=1,2 do 100 k=1,2 s=(-1.)**(i+j+k) t1=y(j)**2+z(k)**2 t2=sqrt(x(i)**2+t1) sum1=sum1-0.5*s*( x(i)*t2+t1*alog(x(i)+t2) ) 100 continue syz=sum1 return end real function sxx(z0,x0,y0,a) real x(2),y(2),z(2) x(1)=-a-x0 x(2)=a-x0 y(1)=-a-y0 y(2)=a-y0 z(1)=-a-z0 z(2)=a-z0 sum1=0. do 100 i=1,2 do 100 j=1,2 do 100 k=1,2 s=(-1.)**(i+j+k) sum1=sum1+s*( q1(y(j),z(k),x(i))+q1(x(i),z(k),y(j)) # -q1(x(i),y(j),z(k)) ) 100 continue sum2=0. do 200 i=1,2 do 200 j=1,2 s=(-1.)**(i+j) sum2=sum2+s*q2(x(i),y(j),z(1),z(2)) 200 continue sxx=sum1+sum2 return end real function syy(y0,z0,x0,a) real x(2),y(2),z(2) x(1)=-a-x0 x(2)=a-x0 y(1)=-a-y0 y(2)=a-y0 z(1)=-a-z0 z(2)=a-z0 sum1=0. do 100 i=1,2 do 100 j=1,2 do 100 k=1,2 s=(-1.)**(i+j+k) sum1=sum1+s*( q1(y(j),z(k),x(i))+q1(x(i),z(k),y(j)) # -q1(x(i),y(j),z(k)) ) 100 continue sum2=0. do 200 i=1,2 do 200 j=1,2 s=(-1.)**(i+j) sum2=sum2+s*q2(x(i),y(j),z(1),z(2)) 200 continue syy=sum1+sum2 return end real function sxy(x0,y0,z0,a) real x(2),y(2),z(2) x(1)=-a-x0 x(2)=a-x0 y(1)=-a-y0 y(2)=a-y0 z(1)=-a-z0 z(2)=a-z0 sum1=0. do 100 i=1,2 do 100 j=1,2 do 100 k=1,2 s=(-1.)**(i+j+k) t1=x(i)**2+y(j)**2 t2=sqrt(z(k)**2+t1) sum1=sum1-0.5*s*( z(k)*t2+t1*alog(z(k)+t2) ) 100 continue sxy=sum1 return end real function q1(a,b,c) t1=a*b t2=c*sqrt(a*a+b*b+c*c) if(t2.eq.0.0.and.t1.gt.0.0)then q1=3.141593/2. else if(t2.eq.0.0.and.t1.lt.0.0)then q1=-3.141593/2. else q1=atan(t1/t2) endif q1=-0.5*c*c*q1 return end real function q2(a,b,c1,c2) common/int1/h t2=c2+sqrt(a*a+b*b+c2*c2) t1=c1+sqrt(a*a+b*b+c1*c1) q2=a*b*alog(t2/t1) return end C*********************************************************************** SUBROUTINE MATINV(A,M,B,N) REAL A(M,M),B(M,N) C======================================================================= C SOLUTION OF REAL SYSTEM OF LINEAR EQUATIONS. C======================================================================= IMA=M-1 DO 100 I=1,IMA AMAX=ABS(A(I,I)) JMAX=I I1=I+1 DO 20 J=I1,M AJI=ABS(A(J,I)) IF(AMAX-AJI)10,20,20 10 AMAX=AJI JMAX=J 20 CONTINUE IF(AMAX)30,140,30 30 IF(I-JMAX)40,70,40 40 DO 50 K=I,M AT=A(I,K) A(I,K)=A(JMAX,K) 50 A(JMAX,K)=AT DO 60 K=1,N AT=B(I,K) B(I,K)=B(JMAX,K) 60 B(JMAX,K)=AT 70 DO 100 J=I1,M FAC=A(J,I)/A(I1-1,I1-1) DO 80 K=I1,M 80 A(J,K)=A(J,K)-FAC*A(I,K) IF(N.LE.0)GO TO 10 DO 90 K=1,N 90 B(J,K)=B(J,K)-FAC*B(I,K) 100 CONTINUE C---------------------------------------TRIANGULAR MATRIX READY--------- IF(ABS(A(M,M)).EQ.0.)GO TO 140 DO 110 K=1,N 110 B(M,K)=B(M,K)/A(M,M) DO 130 I=1,IMA J=M-I K1=J+1 DO 120 K=K1,M DO 120 L=1,N 120 B(J,L)=B(J,L)-A(J,K)*B(K,L) DO 130 L=1,N 130 B(J,L)=B(J,L)/A(J,J) RETURN 140 CONTINUE 150 WRITE(*,1000) 1000 FORMAT(1X,36HDETERMINANT EQUAL TO ZERO IN MATINV.) RETURN END