! This is a Fortran 90 program which solves any three dimensional linear elastostatic (without body forces) problem using constant boundary elements. ! This is a parallelized code (uses BLACS and ScaLAPACK) ! Strongly singular integrals are handled using rigid body modes ! Weakly singular integrals are not handled by any special methods. The program uses 16 by 16 (256 in total) integration points over each of the boundary elements; high number of numerical integration points take care of weak singularity. Program can be easily modified to use 4 by 4, 8 by 8, or 32 by 32 integration points (locations and weights of the integration points for these cases are listed in this code too). Test runs showed that using 4 by 4 integration does not give very accure results whereas 8 by 8 integration is fine. But use 16 by 16 integration to ensure that one does not get erroneous results because of inaccurate integration. There is not much improvement in accuracy by using more than 16 by 16 integration points. ! This code has been tested using Mvapich2-1.8 and Intel Fortran compiler (Version 2011.5.220). Intel Math Kernel Library was also essential. ! This program uses a geometry that is described by 96 triangular boundary elements. However, the program can be very easily modified to solve any problem that uses a geometry that is described by any number of surface triangles. ! Input is through a file named 'inbem96.txt'. The file lists the x, y, and z coordinates of the three vertices of boundary elements that together describe the geometry. ! Output is through files named 'outbem96p.txt' and 'timebem96p.txt'. The file 'outbem96p.txt' contains the solution, i.e., displacements (along x, y, and z direction) or tractions at each of the boundary elements. The file 'timebem96p.txt' shows the time taken to solve the problem, i.e., the Wall-clock time from the first call to MPI_WTIME() to the second call to MPI_WTIME(). Outputs are also written to the screen. ! In this program, element numbers from 65 to 80 are fixed. Element numbers from 1 to 64 are subjected to zero traction, and element numbers from 81 to 96 are subjected to a traction '4' in the y direction. Hence, in total, 16 elements are subjected to displacement boundary condition, and 80 elements are subjected to traction boundary condition. When using this program to solve another problem, suitable modifications have to be made to incorporate appropriate boundary conditions. ! It may help to refer to "Intel Math Kernel Library Reference Manual" (from Intel) and "Parallel ESSL Guide and Reference" (from IBM) to learn the usage of the library routines, and also to learn some theory (like Block-Cyclic Data Distribution) in some cases. program bemconstp implicit none integer :: mypnum,nprocs,icontxt,NUMROC,ia,ja,ib,jb,INDXL2G double precision :: start,finish,time double precision :: E,mu,pi,MPI_WTIME double precision :: dispbcelements(16,4), fbcelements(80,4), xyzofelements(96*3,4),xyzofelements0(4,96*3) double precision,allocatable :: nx(:),ny(:),nz(:),xbar(:),ybar(:),zbar(:),jacobian(:) double precision :: tonedim(16),vonedim(16),weightonedim(16),t(256),v(256),u(256),weight(256),xx(256),y(256),z(256) double precision,allocatable :: K(:,:),F(:,:),Kp(:,:),Fp(:,:) integer :: ndex,vindex,tindex,h,i,j,nelement,sizedispbcelements,sizefbcelements,ndispbcelements,nfbcelements double precision :: x1,x2,x3,y1,y2,y3,z1,z2,z3,d,alpha,beta,gama,sigma,rr,drbydx,drbydy,drbydz,costheta,G,C1,C2,C3,CC double precision :: ITxxE,ITxyE,ITxzE,ITyxE,ITyyE,ITyzE,ITzxE,ITzyE,ITzzE,IUxxE,IUxyE,IUxzE,IUyxE,IUyyE,IUyzE,IUzxE,IUzyE double precision :: IUzzE,Txx,Txy,Txz,Tyx,Tyy,Tyz,Tzx,Tzy,Tzz,Uxx,Uxy,Uxz,Uyx,Uyy,Uyz,Uzx,Uzy,Uzz double precision,allocatable :: CPVITxxE(:),CPVITxyE(:),CPVITxzE(:),CPVITyxE(:),CPVITyyE(:),CPVITyzE(:) double precision,allocatable :: CPVITzxE(:),CPVITzyE(:),CPVITzzE(:) integer :: info,nrhs=1,n=96*3,istat,totalsize,descK(9),descF(9),block double precision :: X(96*3,1),ra(96*3,96*3),ca(96*3,96*3) ! 'nprocsr' is the total number of rows in the process grid. 'nprocsc' is the total number of columns in the process grid. (nprocsr * nprocsc) is equal to the total number of processes ('nprocs') in the parallel program. 'mbK' is the row block size and 'nbK' is the column block size for the 'characteristic' matrix. 'mbF' is the row block size and 'nbF' is the column block size for the right hand size. 'mbK' should be equal to 'mbF', and 'nbF' should be equal to '1'. integer :: mK=96*3,nK=96*3,mbK=16,nbK=16,mF=96*3,nF=1,mbF=16,nbF=1,lldK,lldF,lldKc,lldFc,nprocsr=2,nprocsc=2,myrow,mycol integer,allocatable :: ipiv(:) call blacs_pinfo(mypnum,nprocs) call blacs_get(0,0,icontxt) call blacs_gridinit(icontxt,'R',nprocsr,nprocsc) call blacs_pcoord(icontxt, mypnum, myrow, mycol) ! Enter values of Young's modulus 'E' and Poisson's ratio 'mu' E=200000.0 mu=0.33 do i=1,16 dispbcelements(i,1)=i+64 dispbcelements(i,2)=0 dispbcelements(i,3)=0 dispbcelements(i,4)=0 enddo do i=1,64 fbcelements(i,1)=i fbcelements(i,2)=0 fbcelements(i,3)=0 fbcelements(i,4)=0 enddo do i=65,80 fbcelements(i,1)=i+16 fbcelements(i,2)=0 fbcelements(i,3)=4 fbcelements(i,4)=0 enddo open(101, file="inbem96.txt", status="unknown", action="read", form="formatted", IOSTAT=istat) read (101, *) xyzofelements0 close(101) do i=1,96*3 do j=1,4 xyzofelements(i,j)=xyzofelements0(j,i) enddo enddo pi=3.141592653589793 G=E/(2*(1+mu)) CC=1/(16*pi*G*(1-mu)) C1=3-4*mu C2=1/(8*pi*(1-mu)) C3=1-2*mu totalsize=size(xyzofelements,1) nelement=totalsize/3 sizedispbcelements=size(dispbcelements,1) ndispbcelements=sizedispbcelements sizefbcelements=size(fbcelements,1) nfbcelements=sizefbcelements allocate (nx(nelement)) allocate (ny(nelement)) allocate (nz(nelement)) allocate (xbar(nelement)) allocate (ybar(nelement)) allocate (zbar(nelement)) allocate (jacobian(nelement)) ! For using 4 by 4 integration ! tonedim=(/-0.861136311594 -0.339981043585 0.339981043585 0.861136311594/) ! vonedim=(/-0.861136311594 -0.339981043585 0.339981043585 0.861136311594/) ! weightonedim=(/0.347854845137 0.652145154863 0.652145154863 0.347854845137/) ! For using 8 by 8 integration ! tonedim=(/-0.960289856498 -0.796666477414 -0.525532409916 -0.183434642496 0.183434642496 0.525532409916 0.796666477414 0.960289856498/) ! vonedim=(/-0.960289856498 -0.796666477414 -0.525532409916 -0.183434642496 0.183434642496 0.525532409916 0.796666477414 0.960289856498/) ! weightonedim=(/0.10122853629 0.222381034453 0.313706645878 0.362683783378 0.362683783378 0.313706645878 0.222381034453 0.10122853629/) ! 16 by 16 integration tonedim=(/-0.989400934992, -0.944575023073, -0.865631202388, -0.755404408355, -0.617876244403, -0.458016777657, -0.281603550779, & -0.0950125098376, 0.0950125098376, 0.281603550779, 0.458016777657, 0.617876244403, 0.755404408355, 0.865631202388, & 0.944575023073, 0.989400934992/) vonedim=(/-0.989400934992, -0.944575023073, -0.865631202388, -0.755404408355, -0.617876244403, -0.458016777657, -0.281603550779, & -0.0950125098376, 0.0950125098376, 0.281603550779, 0.458016777657, 0.617876244403, 0.755404408355, 0.865631202388, & 0.944575023073, 0.989400934992/) weightonedim=(/0.027152459411, 0.0622535239372, 0.0951585116838, 0.124628971256, 0.149595988817, 0.169156519395, 0.182603415045, & 0.189450610455, 0.189450610455, 0.182603415045, 0.169156519395, 0.149595988817, 0.124628971256, 0.0951585116838, & 0.0622535239372, 0.027152459411/) ! For using 32 by 32 integration ! tonedim=(/-0.997263861849 -0.985611511545 -0.964762255588 -0.934906075938 -0.896321155766 -0.849367613733 -0.794483795968 -0.73218211874 -0.66304426693 -0.587715757241 -0.506899908932 -0.421351276131 -0.331868602282 -0.239287362252 -0.144471961583 -0.0483076656877 0.0483076656877 0.144471961583 0.239287362252 0.331868602282 0.421351276131 0.506899908932 0.587715757241 0.66304426693 0.73218211874 0.794483795968 0.849367613733 0.896321155766 0.934906075938 0.964762255588 0.985611511545 0.997263861849/) ! vonedim=(/-0.997263861849 -0.985611511545 -0.964762255588 -0.934906075938 -0.896321155766 -0.849367613733 -0.794483795968 -0.73218211874 -0.66304426693 -0.587715757241 -0.506899908932 -0.421351276131 -0.331868602282 -0.239287362252 -0.144471961583 -0.0483076656877 0.0483076656877 0.144471961583 0.239287362252 0.331868602282 0.421351276131 0.506899908932 0.587715757241 0.66304426693 0.73218211874 0.794483795968 0.849367613733 0.896321155766 0.934906075938 0.964762255588 0.985611511545 0.997263861849/) ! weightonedim=(/0.00701814576495 0.0162774265831 0.0253910098329 0.0342745478477 0.0428359896785 0.0509978738117 0.0586839394615 0.0658220603578 0.0723456094297 0.078193695762 0.083311711103 0.0876518688047 0.0911736454878 0.0938441590423 0.0956384754512 0.0965398415811 0.0965398415811 0.0956384754512 0.0938441590423 0.0911736454878 0.0876518688047 0.083311711103 0.078193695762 0.0723456094297 0.0658220603578 0.0586839394615 0.0509978738117 0.0428359896785 0.0342745478477 0.0253910098329 0.0162774265831 0.00701814576495/) ndex=1 do tindex=1,16 do vindex=1,16 t(ndex)=tonedim(tindex) v(ndex)=vonedim(vindex) weight(ndex)=weightonedim(tindex)*weightonedim(vindex) t(ndex)=0.5*t(ndex)+0.5 v(ndex)=0.5*v(ndex)+0.5 ndex=ndex+1 enddo enddo u=t*(1-v) allocate (CPVITxxE(nelement)) allocate (CPVITxyE(nelement)) allocate (CPVITxzE(nelement)) allocate (CPVITyxE(nelement)) allocate (CPVITyyE(nelement)) allocate (CPVITyzE(nelement)) allocate (CPVITzxE(nelement)) allocate (CPVITzyE(nelement)) allocate (CPVITzzE(nelement)) allocate (K(3*nelement,3*nelement)) allocate (F(3*nelement,1)) CPVITxxE=0 CPVITxyE=0 CPVITxzE=0 CPVITyxE=0 CPVITyyE=0 CPVITyzE=0 CPVITzxE=0 CPVITzyE=0 CPVITzzE=0 block=ceiling(real(nelement)/real(nprocs)) call blacs_barrier(icontxt,'A') do i=1,nelement x1=xyzofelements(3*i-2,2) y1=xyzofelements(3*i-2,3) z1=xyzofelements(3*i-2,4) x2=xyzofelements(3*i-1,2) y2=xyzofelements(3*i-1,3) z2=xyzofelements(3*i-1,4) x3=xyzofelements(3*i,2) y3=xyzofelements(3*i,3) z3=xyzofelements(3*i,4) d=(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1))**2+((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1))**2+((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))**2)**0.5 nx(i)=((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1))/d ny(i)=((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1))/d nz(i)=((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))/d alpha=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2) beta=sqrt((x2-x3)**2+(y2-y3)**2+(z2-z3)**2) gama=sqrt((x3-x1)**2+(y3-y1)**2+(z3-z1)**2) sigma=(alpha+beta+gama)/2 jacobian(i)=2*sqrt(sigma*(sigma-alpha)*(sigma-beta)*(sigma-gama)) xbar(i)=(x2-x1)*0.25+(x3-x1)*0.5+x1 ybar(i)=(y2-y1)*0.25+(y3-y1)*0.5+y1 zbar(i)=(z2-z1)*0.25+(z3-z1)*0.5+z1 enddo ! Upto this line, this program runs on every process (processor) ! The following 'if' statement runs on process '0' only call blacs_barrier(icontxt,'A') if (mypnum==0) then start=MPI_WTIME() endif call blacs_barrier(icontxt,'A') ! The following 'do' loop is executed on all processes but for different values of i's (embarrasingly parallel) do i=block*mypnum+1,min(block*mypnum+block,nelement) do j=1,nfbcelements x1=xyzofelements(3*int(fbcelements(j,1))-2,2) y1=xyzofelements(3*int(fbcelements(j,1))-2,3) z1=xyzofelements(3*int(fbcelements(j,1))-2,4) x2=xyzofelements(3*int(fbcelements(j,1))-1,2) y2=xyzofelements(3*int(fbcelements(j,1))-1,3) z2=xyzofelements(3*int(fbcelements(j,1))-1,4) x3=xyzofelements(3*int(fbcelements(j,1)),2) y3=xyzofelements(3*int(fbcelements(j,1)),3) z3=xyzofelements(3*int(fbcelements(j,1)),4) xx=(x2-x1)*u+(x3-x1)*v+x1 y=(y2-y1)*u+(y3-y1)*v+y1 z=(z2-z1)*u+(z3-z1)*v+z1 ITxxE=0 ITxyE=0 ITxzE=0 ITyxE=0 ITyyE=0 ITyzE=0 ITzxE=0 ITzyE=0 ITzzE=0 IUxxE=0 IUxyE=0 IUxzE=0 IUyxE=0 IUyyE=0 IUyzE=0 IUzxE=0 IUzyE=0 IUzzE=0 do h=1,256 rr=sqrt((xx(h)-xbar(i))**2+(y(h)-ybar(i))**2+(z(h)-zbar(i))**2) drbydx=(xx(h)-xbar(i))/rr drbydy=(y(h)-ybar(i))/rr drbydz=(z(h)-zbar(i))/rr costheta=(1/rr)*((xx(h)-xbar(i))*nx(int(fbcelements(j,1)))+(y(h)-ybar(i))*ny(int(fbcelements(j,1)))+ & (z(h)-zbar(i))*nz(int(fbcelements(j,1)))) Txx=(-C2/rr**2)*((C3+3*drbydx**2)*costheta) Tyy=(-C2/rr**2)*((C3+3*drbydy**2)*costheta) Tzz=(-C2/rr**2)*((C3+3*drbydz**2)*costheta) Txy=(-C2/rr**2)*(3*drbydx*drbydy*costheta-C3*(ny(int(fbcelements(j,1)))*drbydx-nx(int(fbcelements(j,1)))*drbydy)) Tyx=(-C2/rr**2)*(3*drbydy*drbydx*costheta-C3*(nx(int(fbcelements(j,1)))*drbydy-ny(int(fbcelements(j,1)))*drbydx)) Tyz=(-C2/rr**2)*(3*drbydy*drbydz*costheta-C3*(nz(int(fbcelements(j,1)))*drbydy-ny(int(fbcelements(j,1)))*drbydz)) Tzy=(-C2/rr**2)*(3*drbydz*drbydy*costheta-C3*(ny(int(fbcelements(j,1)))*drbydz-nz(int(fbcelements(j,1)))*drbydy)) Txz=(-C2/rr**2)*(3*drbydx*drbydz*costheta-C3*(nz(int(fbcelements(j,1)))*drbydx-nx(int(fbcelements(j,1)))*drbydz)) Tzx=(-C2/rr**2)*(3*drbydz*drbydx*costheta-C3*(nx(int(fbcelements(j,1)))*drbydz-nz(int(fbcelements(j,1)))*drbydx)) ITxxE=ITxxE+Txx*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) ITxyE=ITxyE+Txy*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) ITxzE=ITxzE+Txz*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) ITyxE=ITyxE+Tyx*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) ITyyE=ITyyE+Tyy*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) ITyzE=ITyzE+Tyz*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) ITzxE=ITzxE+Tzx*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) ITzyE=ITzyE+Tzy*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) ITzzE=ITzzE+Tzz*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) Uxx=(CC/rr)*(C1+drbydx**2) Uyy=(CC/rr)*(C1+drbydy**2) Uzz=(CC/rr)*(C1+drbydz**2) Uxy=(CC/rr)*drbydx*drbydy Uyx=Uxy Uyz=(CC/rr)*drbydy*drbydz Uzy=Uyz Uxz=(CC/rr)*drbydx*drbydz Uzx=Uxz IUxxE=IUxxE+Uxx*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) IUxyE=IUxyE+Uxy*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) IUxzE=IUxzE+Uxz*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) IUyxE=IUyxE+Uyx*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) IUyyE=IUyyE+Uyy*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) IUyzE=IUyzE+Uyz*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) IUzxE=IUzxE+Uzx*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) IUzyE=IUzyE+Uzy*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) IUzzE=IUzzE+Uzz*(1-v(h))*jacobian(int(fbcelements(j,1)))*weight(h) enddo ITxxE=0.5*ITxxE ITxyE=0.5*ITxyE ITxzE=0.5*ITxzE ITyxE=0.5*ITyxE ITyyE=0.5*ITyyE ITyzE=0.5*ITyzE ITzxE=0.5*ITzxE ITzyE=0.5*ITzyE ITzzE=0.5*ITzzE IUxxE=0.5*IUxxE IUxyE=0.5*IUxyE IUxzE=0.5*IUxzE IUyxE=0.5*IUyxE IUyyE=0.5*IUyyE IUyzE=0.5*IUyzE IUzxE=0.5*IUzxE IUzyE=0.5*IUzyE IUzzE=0.5*IUzzE CPVITxxE(i)=CPVITxxE(i)-ITxxE CPVITxyE(i)=CPVITxyE(i)-ITxyE CPVITxzE(i)=CPVITxzE(i)-ITxzE CPVITyxE(i)=CPVITyxE(i)-ITyxE CPVITyyE(i)=CPVITyyE(i)-ITyyE CPVITyzE(i)=CPVITyzE(i)-ITyzE CPVITzxE(i)=CPVITzxE(i)-ITzxE CPVITzyE(i)=CPVITzyE(i)-ITzyE CPVITzzE(i)=CPVITzzE(i)-ITzzE if (i==int(fbcelements(j,1))) then CPVITxxE(i)=CPVITxxE(i)+ITxxE CPVITxyE(i)=CPVITxyE(i)+ITxyE CPVITxzE(i)=CPVITxzE(i)+ITxzE CPVITyxE(i)=CPVITyxE(i)+ITyxE CPVITyyE(i)=CPVITyyE(i)+ITyyE CPVITyzE(i)=CPVITyzE(i)+ITyzE CPVITzxE(i)=CPVITzxE(i)+ITzxE CPVITzyE(i)=CPVITzyE(i)+ITzyE CPVITzzE(i)=CPVITzzE(i)+ITzzE endif K(3*i-2,3*int(fbcelements(j,1))-2)=ITxxE K(3*i-2,3*int(fbcelements(j,1))-1)=ITxyE K(3*i-2,3*int(fbcelements(j,1)))=ITxzE K(3*i-1,3*int(fbcelements(j,1))-2)=ITyxE K(3*i-1,3*int(fbcelements(j,1))-1)=ITyyE K(3*i-1,3*int(fbcelements(j,1)))=ITyzE K(3*i,3*int(fbcelements(j,1))-2)=ITzxE K(3*i,3*int(fbcelements(j,1))-1)=ITzyE K(3*i,3*int(fbcelements(j,1)))=ITzzE F(3*i-2,1)=F(3*i-2,1)+fbcelements(j,2)*IUxxE+fbcelements(j,3)*IUxyE+fbcelements(j,4)*IUxzE F(3*i-1,1)=F(3*i-1,1)+fbcelements(j,2)*IUyxE+fbcelements(j,3)*IUyyE+fbcelements(j,4)*IUyzE F(3*i,1)=F(3*i,1)+fbcelements(j,2)*IUzxE+fbcelements(j,3)*IUzyE+fbcelements(j,4)*IUzzE enddo do j=1,ndispbcelements x1=xyzofelements(3*int(dispbcelements(j,1))-2,2) y1=xyzofelements(3*int(dispbcelements(j,1))-2,3) z1=xyzofelements(3*int(dispbcelements(j,1))-2,4) x2=xyzofelements(3*int(dispbcelements(j,1))-1,2) y2=xyzofelements(3*int(dispbcelements(j,1))-1,3) z2=xyzofelements(3*int(dispbcelements(j,1))-1,4) x3=xyzofelements(3*int(dispbcelements(j,1)),2) y3=xyzofelements(3*int(dispbcelements(j,1)),3) z3=xyzofelements(3*int(dispbcelements(j,1)),4) xx=(x2-x1)*u+(x3-x1)*v+x1 y=(y2-y1)*u+(y3-y1)*v+y1 z=(z2-z1)*u+(z3-z1)*v+z1 ITxxE=0 ITxyE=0 ITxzE=0 ITyxE=0 ITyyE=0 ITyzE=0 ITzxE=0 ITzyE=0 ITzzE=0 IUxxE=0 IUxyE=0 IUxzE=0 IUyxE=0 IUyyE=0 IUyzE=0 IUzxE=0 IUzyE=0 IUzzE=0 do h=1,256 rr=sqrt((xx(h)-xbar(i))**2+(y(h)-ybar(i))**2+(z(h)-zbar(i))**2) drbydx=(xx(h)-xbar(i))/rr drbydy=(y(h)-ybar(i))/rr drbydz=(z(h)-zbar(i))/rr costheta=(1/rr)*((xx(h)-xbar(i))*nx(int(dispbcelements(j,1)))+(y(h)-ybar(i))*ny(int(dispbcelements(j,1)))+ & (z(h)-zbar(i))*nz(int(dispbcelements(j,1)))) Txx=(-C2/rr**2)*((C3+3*drbydx**2)*costheta) Tyy=(-C2/rr**2)*((C3+3*drbydy**2)*costheta) Tzz=(-C2/rr**2)*((C3+3*drbydz**2)*costheta) Txy=(-C2/rr**2)*(3*drbydx*drbydy*costheta-C3*(ny(int(dispbcelements(j,1)))*drbydx-nx(int(dispbcelements(j,1)))*drbydy)) Tyx=(-C2/rr**2)*(3*drbydy*drbydx*costheta-C3*(nx(int(dispbcelements(j,1)))*drbydy-ny(int(dispbcelements(j,1)))*drbydx)) Tyz=(-C2/rr**2)*(3*drbydy*drbydz*costheta-C3*(nz(int(dispbcelements(j,1)))*drbydy-ny(int(dispbcelements(j,1)))*drbydz)) Tzy=(-C2/rr**2)*(3*drbydz*drbydy*costheta-C3*(ny(int(dispbcelements(j,1)))*drbydz-nz(int(dispbcelements(j,1)))*drbydy)) Txz=(-C2/rr**2)*(3*drbydx*drbydz*costheta-C3*(nz(int(dispbcelements(j,1)))*drbydx-nx(int(dispbcelements(j,1)))*drbydz)) Tzx=(-C2/rr**2)*(3*drbydz*drbydx*costheta-C3*(nx(int(dispbcelements(j,1)))*drbydz & -nz(int(dispbcelements(j,1)))*drbydx)) ITxxE=ITxxE+Txx*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) ITxyE=ITxyE+Txy*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) ITxzE=ITxzE+Txz*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) ITyxE=ITyxE+Tyx*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) ITyyE=ITyyE+Tyy*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) ITyzE=ITyzE+Tyz*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) ITzxE=ITzxE+Tzx*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) ITzyE=ITzyE+Tzy*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) ITzzE=ITzzE+Tzz*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) Uxx=(CC/rr)*(C1+drbydx**2) Uyy=(CC/rr)*(C1+drbydy**2) Uzz=(CC/rr)*(C1+drbydz**2) Uxy=(CC/rr)*drbydx*drbydy Uyx=Uxy Uyz=(CC/rr)*drbydy*drbydz Uzy=Uyz Uxz=(CC/rr)*drbydx*drbydz Uzx=Uxz IUxxE=IUxxE+Uxx*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) IUxyE=IUxyE+Uxy*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) IUxzE=IUxzE+Uxz*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) IUyxE=IUyxE+Uyx*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) IUyyE=IUyyE+Uyy*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) IUyzE=IUyzE+Uyz*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) IUzxE=IUzxE+Uzx*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) IUzyE=IUzyE+Uzy*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) IUzzE=IUzzE+Uzz*(1-v(h))*jacobian(int(dispbcelements(j,1)))*weight(h) enddo ITxxE=0.5*ITxxE ITxyE=0.5*ITxyE ITxzE=0.5*ITxzE ITyxE=0.5*ITyxE ITyyE=0.5*ITyyE ITyzE=0.5*ITyzE ITzxE=0.5*ITzxE ITzyE=0.5*ITzyE ITzzE=0.5*ITzzE IUxxE=0.5*IUxxE IUxyE=0.5*IUxyE IUxzE=0.5*IUxzE IUyxE=0.5*IUyxE IUyyE=0.5*IUyyE IUyzE=0.5*IUyzE IUzxE=0.5*IUzxE IUzyE=0.5*IUzyE IUzzE=0.5*IUzzE CPVITxxE(i)=CPVITxxE(i)-ITxxE CPVITxyE(i)=CPVITxyE(i)-ITxyE CPVITxzE(i)=CPVITxzE(i)-ITxzE CPVITyxE(i)=CPVITyxE(i)-ITyxE CPVITyyE(i)=CPVITyyE(i)-ITyyE CPVITyzE(i)=CPVITyzE(i)-ITyzE CPVITzxE(i)=CPVITzxE(i)-ITzxE CPVITzyE(i)=CPVITzyE(i)-ITzyE CPVITzzE(i)=CPVITzzE(i)-ITzzE if (i==int(fbcelements(j,1))) then CPVITxxE(i)=CPVITxxE(i)+ITxxE CPVITxyE(i)=CPVITxyE(i)+ITxyE CPVITxzE(i)=CPVITxzE(i)+ITxzE CPVITyxE(i)=CPVITyxE(i)+ITyxE CPVITyyE(i)=CPVITyyE(i)+ITyyE CPVITyzE(i)=CPVITyzE(i)+ITyzE CPVITzxE(i)=CPVITzxE(i)+ITzxE CPVITzyE(i)=CPVITzyE(i)+ITzyE CPVITzzE(i)=CPVITzzE(i)+ITzzE endif K(3*i-2,3*int(dispbcelements(j,1))-2)=IUxxE K(3*i-2,3*int(dispbcelements(j,1))-1)=IUxyE K(3*i-2,3*int(dispbcelements(j,1)))=IUxzE K(3*i-1,3*int(dispbcelements(j,1))-2)=IUyxE K(3*i-1,3*int(dispbcelements(j,1))-1)=IUyyE K(3*i-1,3*int(dispbcelements(j,1)))=IUyzE K(3*i,3*int(dispbcelements(j,1))-2)=IUzxE K(3*i,3*int(dispbcelements(j,1))-1)=IUzyE K(3*i,3*int(dispbcelements(j,1)))=IUzzE F(3*i-2,1)=F(3*i-2,1)+dispbcelements(j,2)*ITxxE+dispbcelements(j,3)*ITxyE+dispbcelements(j,4)*ITxzE F(3*i-1,1)=F(3*i-1,1)+dispbcelements(j,2)*ITyxE+dispbcelements(j,3)*ITyyE+dispbcelements(j,4)*ITyzE F(3*i,1)=F(3*i,1)+dispbcelements(j,2)*ITzxE+dispbcelements(j,3)*ITzyE+dispbcelements(j,4)*ITzzE if (i==int(dispbcelements(j,1))) then F(3*i-2,1)=F(3*i-2,1)-dispbcelements(j,2)*ITxxE-dispbcelements(j,3)*ITxyE-dispbcelements(j,4)*ITxzE F(3*i-1,1)=F(3*i-1,1)-dispbcelements(j,2)*ITyxE-dispbcelements(j,3)*ITyyE-dispbcelements(j,4)*ITyzE F(3*i,1)=F(3*i,1)-dispbcelements(j,2)*ITzxE-dispbcelements(j,3)*ITzyE-dispbcelements(j,4)*ITzzE endif enddo enddo call blacs_barrier(icontxt,'A') ! The following 'do' loop is executed on all processes but for different values of i's (embarrasingly parallel) do i=block*mypnum+1,min(block*mypnum+block,nelement) do j=1,nfbcelements if (i==int(fbcelements(j,1))) then K(3*i-2,3*int(fbcelements(j,1))-2)=CPVITxxE(i) K(3*i-2,3*int(fbcelements(j,1))-1)=CPVITxyE(i) K(3*i-2,3*int(fbcelements(j,1)))=CPVITxzE(i) K(3*i-1,3*int(fbcelements(j,1))-2)=CPVITyxE(i) K(3*i-1,3*int(fbcelements(j,1))-1)=CPVITyyE(i) K(3*i-1,3*int(fbcelements(j,1)))=CPVITyzE(i) K(3*i,3*int(fbcelements(j,1))-2)=CPVITzxE(i) K(3*i,3*int(fbcelements(j,1))-1)=CPVITzyE(i) K(3*i,3*int(fbcelements(j,1)))=CPVITzzE(i) endif enddo do j=1,ndispbcelements if (i==int(dispbcelements(j,1))) then F(3*i-2,1)=F(3*i-2,1)+dispbcelements(j,2)*CPVITxxE(i)+dispbcelements(j,3)*CPVITxyE(i)+dispbcelements(j,4)*CPVITxzE(i) F(3*i-1,1)=F(3*i-1,1)+dispbcelements(j,2)*CPVITyxE(i)+dispbcelements(j,3)*CPVITyyE(i)+dispbcelements(j,4)*CPVITyzE(i) F(3*i,1)=F(3*i,1)+dispbcelements(j,2)*CPVITzxE(i)+dispbcelements(j,3)*CPVITzyE(i)+dispbcelements(j,4)*CPVITzzE(i) endif enddo enddo call blacs_barrier(icontxt,'A') ! 'Assemble' 'K' and 'F' from pieces from different processes. 'Assembled' 'K' and 'F' are stored on every process. call dgamx2d(icontxt,'A',' ',3*nelement,3*nelement,K,3*nelement,ra,ca,-1,-1,-1) call dgamx2d(icontxt,'A',' ',3*nelement,1,F,3*nelement,ra,ca,-1,-1,-1) call blacs_barrier(icontxt,'A') ! Determine the sizes of local arrays on each of the processes. Following four lines execute on all of the processes (separately). lldK=NUMROC(3*nelement, mbK, myrow, 0, nprocsr) lldF=NUMROC(3*nelement, mbF, myrow, 0, nprocsr) lldKc=NUMROC(3*nelement, nbK, mycol, 0, nprocsc) lldFc=1 call blacs_barrier(icontxt,'A') ! Allocate local arrays on each of the processes. Following four lines execute on all of the processes (separately). allocate (Kp(lldK,lldKc)) allocate (Fp(lldF,1)) allocate (ipiv(lldK*lldKc+mbK)) call blacs_barrier(icontxt,'A') ! The following two 'do' loops construct 'Kp' and 'Fp' on each of the processes, separately. Effectively, they Block-Cyclically 'distribute' the matrices 'K' and 'F' over all the processes in the process grid. do i=lldK,1,-1 ia=INDXL2G(i,mbK,myrow,0,nprocsr) do j=lldKc,1,-1 ja=INDXL2G(j,nbK,mycol,0,nprocsc) Kp(i,j)=K(ia,ja) enddo enddo do i=lldF,1,-1 ib=INDXL2G(i,mbF,myrow,0,nprocsr) Fp(i,1)=F(ib,1) enddo call blacs_barrier(icontxt,'A') CALL DESCINIT(descK,mK,nK,mbK,nbK,0,0,icontxt,lldK,info) CALL DESCINIT(descF,mF,nF,mbF,nbF,0,0,icontxt,lldF,info) call blacs_barrier(icontxt,'A') ! Solve the system of equations on multiple processors. The matrix 'K' itself is the 'submatrix' in 'pdgesv' call pdgesv(3*nelement,1,Kp,1,1,descK,ipiv,Fp,1,1,descF,info) call blacs_barrier(icontxt,'A') ! 'Reverse' the Block-Cyclic Distribution of the solution ('do' loop executes on all the processes, separately) do i=1,lldF X(INDXL2G(i,mbF,myrow,0,nprocsr),1)=Fp(i,1) enddo call blacs_barrier(icontxt,'A') ! 'Assemble' 'X' from pieces from different processes. 'Assembled' 'X' is stored in process '0' only. 'X' is the solution. call dgamx2d(icontxt,'A',' ',3*nelement,1,X,3*nelement,ra,ca,-1,0,0) call blacs_barrier(icontxt,'A') ! The following 'if' statement executes on process '0' only if (mypnum==0) then finish=MPI_WTIME() time=finish-start open(102, file="outbem96p.txt", status="unknown", action="write", form="formatted", IOSTAT=istat) write (102, *) X close(102) open(103, file="timebem96p.txt", status="unknown", action="write", form="formatted", IOSTAT=istat) write (103, *) time close(103) write (*, *) X write (*, *) time endif call blacs_barrier(icontxt,'A') ! The following two lines execute on each of the processes, separately deallocate (nx,ny,nz,xbar,ybar,zbar,jacobian) deallocate (CPVITxxE,CPVITxyE,CPVITxzE,CPVITyxE,CPVITyyE,CPVITyzE,CPVITzxE,CPVITzyE,CPVITzzE,K,F,Kp,Fp) call blacs_barrier(icontxt,'A') call blacs_exit(0) end program bemconstp