MODULE bie IMPLICIT NONE REAL(KIND=8),PARAMETER :: PI=3.141592653589793238462643383279502884197 CONTAINS SUBROUTINE abmatrel(n,nsum,nmax,kcode,pos,ub,vb,g,h,a,r) !This subroutine rearranges the matrices G and H and produces the ! matrices A and R !One-dimensional array specifying the type of boundary conditions at node I (I= 1,2,...,N) and taking the values: !KCODE(I)=1 when u and v are prescribed !KCODE(I)=2 when tx and ty are prescribed !KCODE(I)=3 when u and ty are prescribed !KCODE(1)=4 when tx and v are prescribed INTEGER(KIND=4),INTENT(IN) :: nmax,nsum INTEGER(KIND=4),DIMENSION(2),INTENT(IN) :: n INTEGER(KIND=4),DIMENSION(2,nmax),INTENT(IN) :: kcode,pos REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: ub,vb REAL(KIND=8),DIMENSION(2,2*nmax,2*nmax),INTENT(IN) :: g,h REAL(KIND=8),DIMENSION(2*nsum,2*nsum),INTENT(OUT) :: a REAL(KIND=8),DIMENSION(2*nsum),INTENT(OUT) :: r INTEGER(KIND=4) :: ii,jj,kk,l,m !Reorder the columns the system of equations and store them in A r=0.0 a=0.0 m=0 DO 70 ii=1,2 IF (ii.gt.1) THEN m=2*n(ii-1) END IF DO 70 kk=1,n(ii) l=pos(ii,kk) SELECT CASE(kcode(ii,kk)) CASE(1) DO 35 jj=1,2*n(ii) a(m+jj,2*l-1)=-g(ii,jj,2*kk-1) a(m+jj,2*l)=-g(ii,jj,2*kk) r(m+jj)=r(m+jj)-h(ii,jj,2*kk-1)*ub(ii,kk)-h(ii,jj,2*kk)*vb(ii,kk) 35 CONTINUE CASE(2) DO 45 jj=1,2*n(ii) a(m+jj,2*l-1)=h(ii,jj,2*kk-1) a(m+jj,2*l)=h(ii,jj,2*kk) r(m+jj)=r(m+jj)+g(ii,jj,2*kk-1)*ub(ii,kk)+g(ii,jj,2*kk)*vb(ii,kk) 45 CONTINUE CASE(3) DO 55 jj=1,2*n(ii) a(m+jj,2*l-1)=-g(ii,jj,2*kk-1) a(m+jj,2*l)=h(ii,jj,2*kk) r(m+jj)=r(m+jj)-h(ii,jj,2*kk-1)*ub(ii,kk)+g(ii,jj,2*kk)*vb(ii,kk) 55 CONTINUE CASE(4) DO 65 jj=1,2*n(ii) a(m+jj,2*l-1)=h(ii,jj,2*kk-1) a(m+jj,2*l)=-g(ii,jj,2*kk) r(m+jj)=r(m+jj)+g(ii,jj,2*kk-1)*ub(ii,kk)-h(ii,jj,2*kk)*vb(ii,kk) 65 CONTINUE CASE(12) DO 66 jj=1,2*n(ii) IF (ii.eq.1) THEN a(m+jj,2*l-1)=-g(ii,jj,2*kk-1) a(m+jj,2*l)=-g(ii,jj,2*kk) a(m+jj,2*(l+1)-1)=h(ii,jj,2*kk-1) a(m+jj,2*(l+1))=h(ii,jj,2*kk) ELSE a(m+jj,2*l-1)=g(ii,jj,2*kk-1) a(m+jj,2*l)=g(ii,jj,2*kk) a(m+jj,2*(l+1)-1)=h(ii,jj,2*kk-1) a(m+jj,2*(l+1))=h(ii,jj,2*kk) END IF 66 CONTINUE END SELECT 70 CONTINUE RETURN END FUNCTION fs(isub,xp,yp,xq,yq,pn) !The fundamental solution s_ij INTEGER(KIND=4),INTENT(IN) :: isub REAL(KIND=8),INTENT(IN) :: xp,yp,xq,yq REAL(KIND=8),DIMENSION(2),INTENT(IN) :: pn REAL(KIND=8) :: ra,dx,dy REAL(KIND=8) :: a1,a2,cosa,sina REAL(KIND=8),DIMENSION(3,3) :: fs a1=-(1+pn(isub))/(4*PI) a2=(1-pn(isub))/(1+pn(isub)) dx=xp-xq dy=yp-yq ra=DSQRT(dx**2+dy**2) cosa=-dx/ra sina=-dy/ra fs(1,1)=(a1/ra)*(a2*cosa+2*cosa**3) fs(1,2)=(a1/ra)*(-a2*sina+2*cosa**2*sina) fs(2,1)=(a1/ra)*(-a2*cosa+2*cosa*sina**2) fs(2,2)=(a1/ra)*(a2*sina+2*sina**3) fs(3,1)=(a1/ra)*(a2*sina+2*cosa**2*sina) fs(3,2)=(a1/ra)*(a2*cosa+2*cosa*sina**2) fs(3,3)=0.0 fs(1,3)=0.0 fs(2,3)=0.0 RETURN END FUNCTION fsb(isub,enx,eny,xp,yp,xq,yq,gl,pn) !The fundamental solution s_ij INTEGER(KIND=4),INTENT(IN) :: isub REAL(KIND=8),INTENT(IN) :: xp,yp,xq,yq,enx,eny REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8) :: ra,rn,dx,dy REAL(KIND=8) :: a1,a3,cosa,sina REAL(KIND=8),DIMENSION(3,3) :: fsb a1=-(1+pn(isub))/(4*PI) a3=-2*gl(isub)*a1 dx=xp-xq dy=yp-yq ra=DSQRT(dx**2+dy**2) cosa=-dx/ra sina=-dy/ra rn=cosa*enx+sina*eny fsb(1,1)=(a3/ra**2)*(2*cosa*rn*(1-4*cosa**2)+(2*cosa**2+1)*enx) fsb(1,2)=(a3/ra**2)*(-8*cosa**2*sina*rn+2*cosa*sina*enx+eny) fsb(2,1)=(a3/ra**2)*(-8*sina**2*cosa*rn+2*cosa*sina*eny+enx) fsb(2,2)=(a3/ra**2)*(2*sina*rn*(1-4*sina**2)+(2*sina**2+1)*eny) fsb(3,1)=fsb(1,2) fsb(3,2)=fsb(2,1) fsb(3,3)=0.0 fsb(1,3)=0.0 fsb(2,3)=0.0 RETURN END FUNCTION ft(isub,enx,eny,xp,yp,xq,yq,pn) !The fundamental solution Tij INTEGER(KIND=4),INTENT(IN) :: isub REAL(KIND=8),INTENT(IN) :: xp,yp,xq,yq,enx,eny REAL(KIND=8),DIMENSION(2),INTENT(IN) :: pn REAL(KIND=8) :: ra,dx,dy,rn,rt REAL(KIND=8) :: a1,a2,cosa,sina REAL(KIND=8), DIMENSION(3,3) :: ft a1=-(1+pn(isub))/(4*PI) a2=(1-pn(isub))/(1+pn(isub)) dx=xp-xq dy=yp-yq ra=DSQRT(dx**2+dy**2) cosa=dx/ra sina=dy/ra rn=cosa*enx+sina*eny rt=sina*enx-cosa*eny ft(1,1)=(a1/ra)*(a2+2*cosa**2)*rn ft(1,2)=(a1/ra)*(2*cosa*sina*rn-a2*rt) ft(2,1)=(a1/ra)*(2*cosa*sina*rn+a2*rt) ft(2,2)=(a1/ra)*(a2+2*sina**2)*rn ft(1,3)=0.0 ft(2,3)=0.0 ft(3,1)=0.0 ft(3,2)=0.0 ft(3,3)=0.0 RETURN END FUNCTION fu(isub,xp,yp,xq,yq,gl,pn) !The fundamental solution Uij INTEGER(KIND=4),INTENT(IN) :: isub REAL(KIND=8),INTENT(IN) :: xp,yp,xq,yq REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8) :: pig,ra,dx,dy REAL(KIND=8) :: cosa,sina REAL(KIND=8),DIMENSION(3,3) :: fu pig=1/(8*PI*gl(isub)) dx=xp-xq dy=yp-yq ra=DSQRT(dx**2+dy**2) cosa=dx/ra sina=dy/ra fu(1,1)=-pig*((3-pn(isub))*DLOG(ra)-(1+pn(isub))*cosa**2) fu(1,2)=pig*(1+pn(isub))*cosa*sina fu(2,2)=-pig*((3-pn(isub))*DLOG(ra)-(1+pn(isub))*sina**2) fu(2,1)=0.0 fu(1,3)=0.0 fu(2,3)=0.0 fu(3,1)=0.0 fu(3,2)=0.0 fu(3,3)=0.0 RETURN END SUBROUTINE gmatrel(n,nmax,x1,y1,x2,y2,xm,ym,gl,pn,g) !This subroutine computes the elements of the G matrix INTEGER(KIND=4),INTENT(IN) :: nmax INTEGER(KIND=4),DIMENSION(2),INTENT(IN) :: n REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,y1,x2,y2,xm,ym REAL(KIND=8),DIMENSION(2,2*nmax,2*nmax),INTENT(OUT) :: g INTEGER(KIND=4) :: ii,jj,kk,isub,iq,ip REAL(KIND=8),DIMENSION(3,3) :: res !compute the elements of matrix G DO 20 ii=1,2 isub=ii DO 20 jj=1,n(ii) iq=jj DO 20 kk=1,n(ii) ip=kk CALL rlintg1(isub,ip,iq,nmax,x1,y1,x2,y2,xm,ym,gl,pn,res) g(ii,2*jj-1,2*kk-1)=res(1,1) g(ii,2*jj-1,2*kk)=res(1,2) g(ii,2*jj,2*kk-1)=res(1,2) g(ii,2*jj,2*kk)=res(2,2) 20 CONTINUE RETURN END SUBROUTINE hmatrel(n,nmax,x1,y1,x2,y2,xm,ym,gl,pn,h) !This subroutine computes the elements of the H matrix INTEGER(KIND=4),INTENT(IN) :: nmax INTEGER(KIND=4),DIMENSION(2),INTENT(IN) :: n REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,y1,x2,y2,xm,ym REAL(KIND=8),DIMENSION(2,2*nmax,2*nmax),INTENT(OUT) :: h INTEGER(KIND=4) :: ii,jj,kk,ip,iq,isub REAL(KIND=8),DIMENSION(3,3) :: res !compute the elements of matrix H DO 40 ii=1,2 isub=ii DO 40 jj=1,n(ii) iq=jj DO 40 kk=1,n(ii) ip=kk CALL rlinth1(isub,ip,iq,nmax,x1,y1,x2,y2,xm,ym,gl,pn,res) h(ii,2*jj-1,2*kk-1)=res(1,1) h(ii,2*jj-1,2*kk)=res(2,1) h(ii,2*jj,2*kk-1)=res(1,2) h(ii,2*jj,2*kk)=res(2,2) 40 CONTINUE RETURN END SUBROUTINE inttij(int) !The singular integral of the fundamental solution Tij: REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: int int(1,1)=0.5 int(1,2)=0.0 int(2,1)=0.0 int(2,2)=0.5 int(1,3)=0.0 int(2,3)=0.0 int(3,1)=0.0 int(3,2)=0.0 int(3,3)=0.0 RETURN END SUBROUTINE intuij(isub,ip,nmax,x1,y1,x2,y2,gl,pn,int) !The singular integral of the fundamental solution Uij: INTEGER(KIND=4),INTENT(IN) :: isub,ip,nmax REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,y1,x2,y2 REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: int REAL(KIND=8) :: spig,dx,dy REAL(KIND=8) :: sl1,cosa,sina dx=x2(isub,ip)-x1(isub,ip) dy=y2(isub,ip)-y1(isub,ip) sl1=DSQRT(dx**2+dy**2) cosa=dx/sl1 sina=dy/sl1 spig=sl1/(8*PI*gl(isub)) int(1,1)=-spig*((3-pn(isub))*(DLOG(sl1/2)-1)-(1+pn(isub))*cosa**2) int(1,2)=spig*(1+pn(isub))*cosa*sina int(2,2)=-spig*((3-pn(isub))*(DLOG(sl1/2)-1)-(1+pn(isub))*sina**2) int(1,3)=0.0 int(2,1)=0.0 int(2,3)=0.0 int(3,1)=0.0 int(3,2)=0.0 int(3,3)=0.0 RETURN END SUBROUTINE quadraturea(isub,iq,nmax,xm,ym,gl,pn,f,a0,b0,int) INTEGER(KIND=4),INTENT(IN) :: isub,iq,f,nmax REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn,a0,b0 REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: xm,ym REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: int INTEGER(KIND=4) :: ii REAL(KIND=8) :: ax,ay,bx,by,sl1 REAL(KIND=8) :: xc,yc,enx,eny REAL(KIND=8),DIMENSION(6) :: xi,wg REAL(KIND=8),DIMENSION(3,3) :: func !DATA xi/-0.99877101, -0.99353017, -0.98412458, -0.97059159, -0.9529877 , & !-0.93138669, -0.90587914, -0.87657202, -0.84358826, -0.8070662 , & !-0.76715903, -0.72403413, -0.67787238, -0.6288674 , -0.57722473, & !-0.52316097, -0.4669029 , -0.40868648, -0.34875589, -0.28736249, & !-0.22476379, -0.16122236, -0.0970047 , -0.03238017, 0.03238017, & ! 0.0970047 , 0.16122236, 0.22476379, 0.28736249, 0.34875589, & ! 0.40868648, 0.4669029 , 0.52316097, 0.57722473, 0.6288674 , & ! 0.67787238, 0.72403413, 0.76715903, 0.8070662 , 0.84358826, & ! 0.87657202, 0.90587914, 0.93138669, 0.9529877 , 0.97059159, & ! 0.98412458, 0.99353017, 0.99877101/ !DATA wg/ 0.00315335, 0.00732755, 0.01147723, 0.01557932, 0.01961616,& !0.02357076, 0.02742651, 0.03116723, 0.03477722, 0.03824135,& !0.04154508, 0.04467456, 0.04761666, 0.05035904, 0.05289019,& !0.0551995 , 0.05727729, 0.05911484, 0.06070444, 0.06203942,& !0.06311419, 0.06392424, 0.06446616, 0.0647377 , 0.0647377 ,& !0.06446616, 0.06392424, 0.06311419, 0.06203942, 0.06070444,& !0.05911484, 0.05727729, 0.0551995 , 0.05289019, 0.05035904,& !0.04761666, 0.04467456, 0.04154508, 0.03824135, 0.03477722,& !0.03116723, 0.02742651, 0.02357076, 0.01961616, 0.01557932,& !0.01147723, 0.00732755, 0.00315335/ DATA xi/-0.9324695142031520278123,-0.6612093864662645136614,-0.2386191860831969086305, & 0.2386191860831969086305,0.6612093864662645136614,0.9324695142031520278123/ DATA wg/0.1713244923791703450403,0.36076157304813860757,0.46791393457269104739, & 0.46791393457269104739,0.36076157304813860757,0.1713244923791703450403/ !DATA xi/-0.86113631,-0.33998104,0.33998104,0.86113631/ !DATA wg/0.34785485,0.65214515,0.65214515,0.34785485/ ax=(b0(1)-a0(1))/2 ay=(b0(2)-a0(2))/2 bx=(b0(1)+a0(1))/2 by=(b0(2)+a0(2))/2 !Compute the Jacobian sl1=DSQRT(ax**2+ay**2) !Compute the line integral int=0.0 DO 30 ii=1,6 xc=ax*xi(ii)+bx yc=ay*xi(ii)+by SELECT CASE(f) CASE(1) func=fu(isub,xc,yc,xm(isub,iq),ym(isub,iq),gl,pn) CASE(2) enx=ay/sl1 eny=-ax/sl1 func=ft(isub,enx,eny,xc,yc,xm(isub,iq),ym(isub,iq),pn) CASE(3) enx=ay/sl1 eny=-ax/sl1 func=fs(isub,xc,yc,xm(isub,iq),ym(isub,iq),pn) CASE(4) enx=ay/sl1 eny=-ax/sl1 func=fsb(isub,enx,eny,xc,yc,xm(isub,iq),ym(isub,iq),gl,pn) END SELECT int=int+func*wg(ii)*sl1 30 CONTINUE RETURN END SUBROUTINE quadratureb(isub,iq,nmax,xm,ym,gl,pn,f,a0,b0,int) INTEGER(KIND=4),INTENT(IN) :: isub,iq,f,nmax REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn,a0,b0 REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: xm,ym REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: int INTEGER(KIND=4) :: ii REAL(KIND=8) :: ax,ay,bx,by,sl1 REAL(KIND=8) :: xc,yc,enx,eny REAL(KIND=8),DIMENSION(48) :: wg REAL(KIND=8),DIMENSION(3,3) :: func REAL(KIND=8),DIMENSION(2) :: vect_ab,vect_ay REAL(KIND=8) :: abs_ab,abs_ay,abs_ac REAL(KIND=8) :: c,r0,alpha REAL(KIND=8) :: sl2,sl2minus1,sl2plus1 REAL(KIND=8) :: eta1,eta2,eta1a,eta1b,eta2a,eta2b REAL(KIND=8) :: eta,etaa,etab REAL(KIND=8) :: xi,xi1,xi2 REAL(KIND=8),DIMENSION(48) :: zeta REAL(KIND=8) :: xc1,yc1,xc2,yc2 REAL(KIND=8),DIMENSION(3,3) :: func1,func2 REAL(KIND=8),DIMENSION(3,3) :: integral1,integral2 DATA zeta/-0.99877101, -0.99353017, -0.98412458, -0.97059159, -0.9529877 , & -0.93138669, -0.90587914, -0.87657202, -0.84358826, -0.8070662 , & -0.76715903, -0.72403413, -0.67787238, -0.6288674 , -0.57722473, & -0.52316097, -0.4669029 , -0.40868648, -0.34875589, -0.28736249, & -0.22476379, -0.16122236, -0.0970047 , -0.03238017, 0.03238017, & 0.0970047 , 0.16122236, 0.22476379, 0.28736249, 0.34875589, & 0.40868648, 0.4669029 , 0.52316097, 0.57722473, 0.6288674 , & 0.67787238, 0.72403413, 0.76715903, 0.8070662 , 0.84358826, & 0.87657202, 0.90587914, 0.93138669, 0.9529877 , 0.97059159, & 0.98412458, 0.99353017, 0.99877101/ DATA wg/ 0.00315335, 0.00732755, 0.01147723, 0.01557932, 0.01961616,& 0.02357076, 0.02742651, 0.03116723, 0.03477722, 0.03824135,& 0.04154508, 0.04467456, 0.04761666, 0.05035904, 0.05289019,& 0.0551995 , 0.05727729, 0.05911484, 0.06070444, 0.06203942,& 0.06311419, 0.06392424, 0.06446616, 0.0647377 , 0.0647377 ,& 0.06446616, 0.06392424, 0.06311419, 0.06203942, 0.06070444,& 0.05911484, 0.05727729, 0.0551995 , 0.05289019, 0.05035904,& 0.04761666, 0.04467456, 0.04154508, 0.03824135, 0.03477722,& 0.03116723, 0.02742651, 0.02357076, 0.01961616, 0.01557932,& 0.01147723, 0.00732755, 0.00315335/ !DATA zeta/-0.9324695142031520278123,-0.6612093864662645136614,-0.2386191860831969086305, & ! 0.2386191860831969086305,0.6612093864662645136614,0.9324695142031520278123/ !DATA wg/0.1713244923791703450403,0.36076157304813860757,0.46791393457269104739, & ! 0.46791393457269104739,0.36076157304813860757,0.1713244923791703450403/ !DATA zeta/-0.86113631,-0.33998104,0.33998104,0.86113631/ !DATA wg/0.34785485,0.65214515,0.65214515,0.34785485/ vect_ab(1)=b0(1)-a0(1) vect_ab(2)=b0(2)-a0(2) vect_ay(1)=xm(isub,iq)-a0(1) vect_ay(2)=ym(isub,iq)-a0(2) abs_ab=DSQRT(vect_ab(1)**2+vect_ab(2)**2) abs_ay=DSQRT(vect_ay(1)**2+vect_ay(2)**2) abs_ac=1/abs_ab*ABS(vect_ab(1)*vect_ay(1)+vect_ab(2)*vect_ay(2)) c=-1.+2.*abs_ac/abs_ab r0=DSQRT(abs_ay**2-abs_ac**2) ax=vect_ab(1)/2. ay=vect_ab(2)/2. bx=(b0(1)+a0(1))/2. by=(b0(2)+a0(2))/2. !jacobian x->xi sl1=DSQRT(ax**2+ay**2) alpha=r0/sl1 !jacobian xi->eta v bodech +-1. sl2plus1=DSQRT(alpha**2+(1.-c)**2) sl2minus1=DSQRT(alpha**2+(-1.-c)**2) IF (ABS(c).le.1.) THEN !jacobian eta->zeta eta1a=.5*(DLOG(alpha)-DLOG(sl2minus1-(1.+c))) eta1b=.5*(DLOG(alpha)+DLOG(sl2minus1-(1.+c))) eta2a=.5*(DLOG(sl2plus1+(1.-c))-DLOG(alpha)) eta2b=.5*(DLOG(sl2plus1+(1.-c))+DLOG(alpha)) !vypocet krivkoveho integralu integral1=0.0 integral2=0.0 DO 80 ii=1,48 eta1=eta1a*zeta(ii)+eta1b eta2=eta2a*zeta(ii)+eta2b xi1=.5*(DEXP(eta1)-alpha**2*DEXP(-eta1))+c xi2=.5*(DEXP(eta2)-alpha**2*DEXP(-eta2))+c xc1=ax*xi1+bx yc1=ay*xi1+by xc2=ax*xi2+bx yc2=ay*xi2+by SELECT CASE(f) CASE(1) func1=fu(isub,xc1,yc1,xm(isub,iq),ym(isub,iq),gl,pn) func2=fu(isub,xc2,yc2,xm(isub,iq),ym(isub,iq),gl,pn) CASE(2) enx=ay/sl1 eny=-ax/sl1 func1=ft(isub,enx,eny,xc1,yc1,xm(isub,iq),ym(isub,iq),pn) func2=ft(isub,enx,eny,xc2,yc2,xm(isub,iq),ym(isub,iq),pn) CASE(3) enx=ay/sl1 eny=-ax/sl1 func1=fs(isub,xc1,yc1,xm(isub,iq),ym(isub,iq),pn) func2=fs(isub,xc2,yc2,xm(isub,iq),ym(isub,iq),pn) CASE(4) enx=ay/sl1 eny=-ax/sl1 func1=fsb(isub,enx,eny,xc1,yc1,xm(isub,iq),ym(isub,iq),gl,pn) func2=fsb(isub,enx,eny,xc2,yc2,xm(isub,iq),ym(isub,iq),gl,pn) END SELECT sl2=DSQRT(alpha**2+(xi1-c)**2) integral1=integral1+func1*wg(ii)*sl1*sl2*eta1a sl2=DSQRT(alpha**2+(xi2-c)**2) integral2=integral2+func2*wg(ii)*sl1*sl2*eta2a 80 CONTINUE int=integral1+integral2 ELSE !jacobian eta->zeta etaa=.5*(DLOG(sl2plus1+(1.-c))-DLOG(sl2minus1-(1.+c))) etab=.5*(DLOG(sl2plus1+(1.-c))+DLOG(sl2minus1-(1.+c))) !vypocet krivkoveho integralu int=0.0 DO 90 ii=1,48 eta=etaa*zeta(ii)+etab xi=.5*(DEXP(eta)-alpha**2*DEXP(-eta))+c xc=ax*xi+bx yc=ay*xi+by SELECT CASE(f) CASE(1) func=fu(isub,xc,yc,xm(isub,iq),ym(isub,iq),gl,pn) CASE(2) enx=ay/sl1 eny=-ax/sl1 func=ft(isub,enx,eny,xc,yc,xm(isub,iq),ym(isub,iq),pn) CASE(3) enx=ay/sl1 eny=-ax/sl1 func=fs(isub,xc,yc,xm(isub,iq),ym(isub,iq),pn) CASE(4) enx=ay/sl1 eny=-ax/sl1 func=fsb(isub,enx,eny,xc,yc,xm(isub,iq),ym(isub,iq),gl,pn) END SELECT sl2=DSQRT(alpha**2+(xi-c)**2) int=int+func*wg(ii)*sl1*sl2*etaa 90 CONTINUE END IF RETURN END SUBROUTINE reorderel(n,nsum,nmax,kcode,pos,ub,vb,r,u,v,tx,ty) !This subroutine rearranges the arrays UB and VB !in such a way that all boundary values of u and v !are stored in UB and VB, while all boundary !values of tx and ty in TXB and TYB INTEGER(KIND=4),INTENT(IN) :: nmax,nsum INTEGER(KIND=4),DIMENSION(2),INTENT(IN) :: n INTEGER(KIND=4),DIMENSION(2,nmax),INTENT(IN) :: kcode,pos REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: ub,vb REAL(KIND=8),DIMENSION(4*nsum),INTENT(IN) :: r REAL(KIND=8),DIMENSION(2,nmax),INTENT(OUT) :: u,v,tx,ty INTEGER(KIND=4) :: ii,jj,l DO 50 ii=1,2 DO 50 jj=1,n(ii) l=pos(ii,jj) SELECT CASE(kcode(ii,jj)) CASE(1) u(ii,jj)=ub(ii,jj) v(ii,jj)=vb(ii,jj) tx(ii,jj)=r(2*l-1) ty(ii,jj)=r(2*l) CASE(2) u(ii,jj)=r(2*l-1) v(ii,jj)=r(2*l) tx(ii,jj)=ub(ii,jj) ty(ii,jj)=vb(ii,jj) CASE(3) u(ii,jj)=ub(ii,jj) v(ii,jj)=r(2*l) ty(ii,jj)=vb(ii,jj) tx(ii,jj)=r(2*l-1) CASE(4) u(ii,jj)=r(2*l-1) v(ii,jj)=vb(ii,jj) tx(ii,jj)=ub(ii,jj) ty(ii,jj)=r(2*l) CASE(12) IF (ii.eq.1) THEN tx(ii,jj)=r(2*l-1) ty(ii,jj)=r(2*l) u(ii,jj)=r(2*(l+1)-1) v(ii,jj)=r(2*(l+1)) ELSE tx(ii,jj)=-r(2*l-1) ty(ii,jj)=-r(2*l) u(ii,jj)=r(2*(l+1)-1) v(ii,jj)=r(2*(l+1)) END IF END SELECT 50 CONTINUE RETURN END SUBROUTINE rlintg1(isub,ip,iq,nmax,x1,y1,x2,y2,xm,ym,gl,pn,res) !This subroutsine computes the elements of the matrix G INTEGER(KIND=4),INTENT(IN) :: isub,ip,iq,nmax REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,x2,y1,y2,xm,ym REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: res INTEGER(KIND=4) :: f REAL(KIND=8),DIMENSION(2) :: a0,b0 f=1 a0(1)=x1(isub,ip) a0(2)=y1(isub,ip) b0(1)=x2(isub,ip) b0(2)=y2(isub,ip) IF (ip.eq.iq) THEN CALL intuij(isub,ip,nmax,x1,y1,x2,y2,gl,pn,res) ELSE CALL quadraturea(isub,iq,nmax,xm,ym,gl,pn,f,a0,b0,res) END IF RETURN END SUBROUTINE rlintg2(isub,ip,iq,nin,xin,yin,nmax,x1,y1,x2,y2,gl,pn,res) !This subroutsine computes the elements of the matrix G INTEGER(KIND=4),INTENT(IN) :: isub,ip,iq,nin,nmax REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,x2,y1,y2 REAL(KIND=8),DIMENSION(2,nin),INTENT(IN) :: xin,yin REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: res REAL(KIND=8) :: eps1,eps2,eps3,l0 INTEGER(KIND=4) :: f REAL(KIND=8),DIMENSION(2) :: a0,b0,c0 f=1 a0(1)=x1(isub,ip) a0(2)=y1(isub,ip) b0(1)=x2(isub,ip) b0(2)=y2(isub,ip) c0(1)=(a0(1)+b0(1))/2. c0(2)=(a0(2)+b0(2))/2. l0=DSQRT((b0(1)-a0(1))**2+(b0(2)-a0(2))**2)/2. eps1=DSQRT((a0(1)-xin(isub,iq))**2+(a0(2)-yin(isub,iq))**2) eps2=DSQRT((b0(1)-xin(isub,iq))**2+(b0(2)-yin(isub,iq))**2) eps3=DSQRT((c0(1)-xin(isub,iq))**2+(c0(2)-yin(isub,iq))**2) IF ((eps1.lt.l0).or.(eps2.lt.l0).or.(eps3.lt.l0)) THEN CALL quadratureb(isub,iq,nin,xin,yin,gl,pn,f,a0,b0,res) ELSE CALL quadraturea(isub,iq,nin,xin,yin,gl,pn,f,a0,b0,res) ENDIF RETURN END SUBROUTINE rlinth1(isub,ip,iq,nmax,x1,y1,x2,y2,xm,ym,gl,pn,res) !This subroutsine computes the elements of the matrix H INTEGER(KIND=4),INTENT(IN) :: isub,ip,iq,nmax REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,x2,y1,y2,xm,ym REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: res INTEGER(KIND=4) :: f REAL(KIND=8),DIMENSION(2) :: a0,b0 f=2 a0(1)=x1(isub,ip) a0(2)=y1(isub,ip) b0(1)=x2(isub,ip) b0(2)=y2(isub,ip) IF (ip.eq.iq) THEN CALL inttij(res) ELSE CALL quadraturea(isub,iq,nmax,xm,ym,gl,pn,f,a0,b0,res) END IF RETURN END SUBROUTINE rlinth2(isub,ip,iq,nin,xin,yin,nmax,x1,y1,x2,y2,gl,pn,res) !This subroutsine computes the elements of the matrix G INTEGER(KIND=4),INTENT(IN) :: isub,ip,iq,nin,nmax REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,x2,y1,y2 REAL(KIND=8),DIMENSION(2,nin),INTENT(IN) :: xin,yin REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: res REAL(KIND=8) :: l0,eps1,eps2,eps3 INTEGER(KIND=4) :: f REAL(KIND=8),DIMENSION(2) :: a0,b0,c0 f=2 a0(1)=x1(isub,ip) a0(2)=y1(isub,ip) b0(1)=x2(isub,ip) b0(2)=y2(isub,ip) c0(1)=(a0(1)+b0(1))/2. c0(2)=(a0(2)+b0(2))/2. l0=DSQRT((b0(1)-a0(1))**2+(b0(2)-a0(2))**2)/2. eps1=DSQRT((a0(1)-xin(isub,iq))**2+(a0(2)-yin(isub,iq))**2) eps2=DSQRT((b0(1)-xin(isub,iq))**2+(b0(2)-yin(isub,iq))**2) eps3=DSQRT((c0(1)-xin(isub,iq))**2+(c0(2)-yin(isub,iq))**2) IF ((eps1.lt.l0).or.(eps2.lt.l0).or.(eps3.lt.l0)) THEN CALL quadratureb(isub,iq,nin,xin,yin,gl,pn,f,a0,b0,res) ELSE CALL quadraturea(isub,iq,nin,xin,yin,gl,pn,f,a0,b0,res) ENDIF RETURN END SUBROUTINE slinth2(isub,ip,iq,nin,xin,yin,nmax,x1,y1,x2,y2,gl,pn,res,resb) !This subroutsine computes the elements of the matrix G and H INTEGER(KIND=4),INTENT(IN) :: isub,ip,iq,nin,nmax REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,x2,y1,y2 REAL(KIND=8),DIMENSION(2,nin),INTENT(IN) :: xin,yin REAL(KIND=8),DIMENSION(3,3),INTENT(OUT) :: res,resb INTEGER(KIND=4) :: f REAL(KIND=8),DIMENSION(2) :: a0,b0,c0 REAL(KIND=8) :: l0,eps1,eps2,eps3 a0(1)=x1(isub,ip) a0(2)=y1(isub,ip) b0(1)=x2(isub,ip) b0(2)=y2(isub,ip) c0(1)=(a0(1)+b0(1))/2. c0(2)=(a0(2)+b0(2))/2. l0=DSQRT((b0(1)-a0(1))**2+(b0(2)-a0(2))**2)/2. eps1=DSQRT((a0(1)-xin(isub,iq))**2+(a0(2)-yin(isub,iq))**2) eps2=DSQRT((b0(1)-xin(isub,iq))**2+(b0(2)-yin(isub,iq))**2) eps3=DSQRT((c0(1)-xin(isub,iq))**2+(c0(2)-yin(isub,iq))**2) IF ((eps1.lt.l0).or.(eps2.lt.l0).or.(eps3.lt.l0)) THEN f=3 CALL quadratureb(isub,iq,nin,xin,yin,gl,pn,f,a0,b0,res) f=4 CALL quadratureb(isub,iq,nin,xin,yin,gl,pn,f,a0,b0,resb) ELSE f=3 CALL quadraturea(isub,iq,nin,xin,yin,gl,pn,f,a0,b0,res) f=4 CALL quadraturea(isub,iq,nin,xin,yin,gl,pn,f,a0,b0,resb) ENDIF res(1,3)=0.0 res(2,3)=0.0 res(3,3)=0.0 resb(1,3)=0.0 resb(2,3)=0.0 resb(3,3)=0.0 RETURN END SUBROUTINE uvinter(isub,nin,xin,yin,n,nmax,x1,y1,x2,y2,gl,pn,ub,vb,txb,tyb,uin,vin) !This subroutine computes the values of the displacements u and v !at the internal points INTEGER(KIND=4),INTENT(IN) :: isub,nin,nmax INTEGER(KIND=4),DIMENSION(2),INTENT(IN) :: n REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(nin),INTENT(IN) :: xin,yin REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,y1,x2,y2,ub,vb,txb,tyb REAL(KIND=8),DIMENSION(nin),INTENT(OUT) :: uin,vin INTEGER(KIND=4) :: ii,jj,iq,ip REAL(KIND=8),DIMENSION(2,nin) :: xinm,yinm REAL(KIND=8),DIMENSION(3,3) :: gin,hin !Compute the values of u and v at the internal points uin=0.0 vin=0.0 xinm=0.0 yinm=0.0 xinm(isub,:)=xin yinm(isub,:)=yin DO 10 ii=1,nin iq=ii DO 10 jj=1,n(isub) ip=jj CALL rlintg2(isub,ip,iq,nin,xinm,yinm,nmax,x1,y1,x2,y2,gl,pn,gin) CALL rlinth2(isub,ip,iq,nin,xinm,yinm,nmax,x1,y1,x2,y2,gl,pn,hin) uin(ii)=uin(ii)-hin(1,1)*ub(isub,jj)-hin(2,1)*vb(isub,jj)+gin(1,1)*txb(isub,jj)+gin(1,2)*tyb(isub,jj) vin(ii)=vin(ii)-hin(1,2)*ub(isub,jj)-hin(2,2)*vb(isub,jj)+gin(1,2)*txb(isub,jj)+gin(2,2)*tyb(isub,jj) 10 CONTINUE RETURN END SUBROUTINE stressin(isub,nin,xin,yin,n,nmax,x1,y1,x2,y2,gl,pn,ub,vb,txb,tyb,sxin,syin,sxyin) !This subroutine computes the stresses at the internal points INTEGER(KIND=4),INTENT(IN) :: isub,nin,nmax INTEGER(KIND=4),DIMENSION(2),INTENT(IN) :: n REAL(KIND=8),DIMENSION(2),INTENT(IN) :: gl,pn REAL(KIND=8),DIMENSION(nin),INTENT(IN) :: xin,yin REAL(KIND=8),DIMENSION(2,nmax),INTENT(IN) :: x1,y1,x2,y2,ub,vb,txb,tyb REAL(KIND=8),DIMENSION(nin),INTENT(OUT) :: sxin,syin,sxyin INTEGER(KIND=4) :: ii,jj,ip,iq REAL(KIND=8),DIMENSION(2,nin) :: xinm,yinm REAL(KIND=8),DIMENSION(3,3) :: res,resb res=0.0 resb=0.0 sxin=0.0 syin=0.0 sxyin=0.0 xinm=0.0 yinm=0.0 xinm(isub,:)=xin yinm(isub,:)=yin DO 10 ii=1,nin iq=ii DO 10 jj=1,n(isub) ip=jj CALL slinth2(isub,ip,iq,nin,xinm,yinm,nmax,x1,y1,x2,y2,gl,pn,res,resb) sxin(ii)=sxin(ii)+res(1,1)*txb(isub,jj)+res(1,2)*tyb(isub,jj)-resb(1,1)*ub(isub,jj)-resb(1,2)*vb(isub,jj) syin(ii)=syin(ii)+res(2,1)*txb(isub,jj)+res(2,2)*tyb(isub,jj)-resb(2,1)*ub(isub,jj)-resb(2,2)*vb(isub,jj) sxyin(ii)=sxyin(ii)+res(3,1)*txb(isub,jj)+res(3,2)*tyb(isub,jj)-resb(3,1)*ub(isub,jj)-resb(3,2)*vb(isub,jj) 10 CONTINUE RETURN END END MODULE bie