! This is a Fortran 90 program which solves any three dimensional linear elastostatic (without body forces) problem using constant boundary elements. ! This is a serial code ! 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 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 'outbem96.txt' and 'timebem96.txt'. The file 'outbem96.txt' contains the solution, i.e., displacements (along x, y, and z direction) or tractions at each of the boundary elements. The file 'timebem96.txt' shows the time taken to solve the problem, i.e., the Wall-clock time from the first call to DCLOCK() to the second call to DCLOCK(). ! 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) to learn the usage of the library routine 'dgesvx'. program bemconst USE IFPORT implicit none double precision :: start,finish double precision :: E,mu,pi 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(:,:) 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,lda=96*3,ldb=96*3,nrhs=1,n=96*3,ipiv(96*3),ldaf=96*3,ldx=96*3,iwork(96*3),istat,totalsize double precision :: af(96*3,96*3),work(4*96*3),r(96*3),c(96*3),x(96*3,1),rcond,ferr,berr character :: fact='E',trans='N',equed='n' ! Enter the value of Young's modulus and Poisson's ratio 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)) 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 ! For 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 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 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 ! DCLOCK() requires the module 'IFPORT' which is included in the Intel compiler start=DCLOCK() do i=1,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 do i=1,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 !Solve !soln=K\F ! Solve the system of equations call dgesvx(fact,trans,n,nrhs,K,lda,af,ldaf,ipiv,equed,r,c,F,ldb,x,ldx,rcond,ferr,berr,work,iwork,info) finish=DCLOCK() ! Display unknown dofs for each element open(102, file="outbem96.txt", status="unknown", action="write", form="formatted", IOSTAT=istat) write (102, *) x close(102) deallocate (nx,ny,nz,xbar,ybar,zbar,jacobian) deallocate (CPVITxxE,CPVITxyE,CPVITxzE,CPVITyxE,CPVITyyE,CPVITyzE,CPVITzxE,CPVITzyE,CPVITzzE,K,F) open(103, file="timebem96.txt", status="unknown", action="write", form="formatted", IOSTAT=istat) write (103, *) finish-start close(103) end program bemconst