Civil Engineering Reference
In-Depth Information
CALL shape_der(derf,points,i) ;jac=MATMUL(derf,coordf)
det=determinant(jac) ; CALL invert(jac)
derivf=MATMUL(jac,derf) ; rowf(1,:) = derivf(1,:)
c12 = c12 + MATMUL(funny,rowf)*det*weights(i)/rho
rowf(1,:) = derivf(2,:)
c32 = c32 + MATMUL(funny,rowf)*det*weights(i)/rho
rowf(1,:) = derivf(3,:)
c42 = c42 + MATMUL(funny,rowf)*det*weights(i)/rho
c21 = c21 + MATMUL(funnyf,row1)*det*weights(i)
c23 = c23 + MATMUL(funnyf,row2)*det*weights(i)
c24 = c24 + MATMUL(funnyf,row3)*det*weights(i)
END DO gauss_points_1
CALL formupvw(ke,c11,c12,c21,c23,c32,c24,c42); storke_pp(:,:,iel)=ke
END DO elements_2 ; diag_tmp = .0_iwp
elements_2a: DO iel=1,nels_pp
DO k=1,ntot; diag_tmp(k,iel) = diag_tmp(k,iel) + ke(k,k) ; END DO
END DO elements_2a ; CALL scatter(diag_pp,diag_tmp)
!----------- prescribed values of velocity and pressure ------------------
DO i=1,num_no ; k = no_local(i) - ieq_start +1
diag_pp(k)=diag_pp(k)+penalty
b_pp(k)=diag_pp(k)*val(no_index_start+i-1);store_pp(k)= diag_pp(k)
END DO
!--- solve the simultaneous equations element by element using BiCGSTAB---
!-------------------
initialisation phase
-------------------------
IF(iters==1) x_pp = x0 ; pmul_pp = .0_iwp
y_pp = x_pp ; y1_pp = .0_iwp ; CALL gather(y_pp,pmul_pp)
elements_3 : DO iel = 1 , nels_pp
ke = storke_pp( : , : , iel)
utemp_pp(:,iel)= MATMUL(ke,pmul_pp(:,iel))
END DO elements_3 ; CALL scatter(y1_pp,utemp_pp)
DO i=1,num_no ; k=no_local(i)-ieq_start+1
y1_pp(k) = y_pp(k) * store_pp(k)
END DO
y_pp=y1_pp; rt_pp = b_pp - y_pp
r_pp=.0_iwp;r_pp(:,1)=rt_pp ;u_pp=.0_iwp; gama=1.0_iwp; omega=1.0_iwp
k = 0;norm_r=norm_p(rt_pp);r0_norm=norm_r; error=1.0_iwp;cjiters = 0
!------------------- bicgstab(ell) iterations ------------------------
bicg_iterations : DO
cjiters = cjiters + 1 ; cj_converged = error < cjtol
IF(cjiters==cjits.OR. cj_converged) EXIT
gama = - omega*gama ; y_pp = r_pp(:,1)
DOj=1,ell
rho1 = DOT_PRODUCT_P(rt_pp,y_pp) ; beta = rho1/gama
u_pp(:,1:j)=r_pp(:,1:j)-beta * u_pp(:,1:j); pmul_pp=.0_iwp
y_pp = u_pp(:,j); y1_pp = .0_iwp; CALL gather(y_pp,pmul_pp)
elements_4 : DO iel = 1 , nels_pp
ke = storke_pp(: , : , iel)
utemp_pp(:,iel)=MATMUL(ke,pmul_pp(:,iel))
END DO elements_4 ; CALL scatter(y1_pp,utemp_pp)
DO i=1,num_no ; l=no_local(i) -ieq_start +1
y1_pp(l) = y_pp(l) * store_pp(l)
END DO
y_pp=y1_pp; u_pp(:,j+1) = y_pp
gama = DOT_PRODUCT_P(rt_pp,y_pp) ; alpha = rho1/gama
x_pp=x_pp+ alpha * u_pp(:,1)
r_pp(:,1:j)= r_pp(:,1:j)-alpha*u_pp(:,2:j+1);pmul_pp=.0_iwp
y_pp = r_pp(:,j); y1_pp = .0_iwp; CALL gather(y_pp,pmul_pp)
elements_5 : DO iel = 1 , nels_pp
