Civil Engineering Reference
In-Depth Information
!-------------------
solve the simultaneous equations by pcg -----------
cjiters = 0
conjugate_gradients: DO
cjiters = cjiters + 1 ; u_pp = .0_iwp
; pmul_pp = .0_iwp
CALL gather(p_pp,pmul_pp)
elements_3 : DO iel=1,nels_pp
utemp_pp(:,iel)=MATMUL(storkm_pp(:,:,iel),pmul_pp(:,iel))
END DO elements_3 ; CALL scatter(u_pp,utemp_pp)
!-------------------------------pcg process ------------------------------
up =DOT_PRODUCT_P(r_pp,d_pp); alpha=up/DOT_PRODUCT_P(p_pp,u_pp)
xnew_pp = x_pp + p_pp* alpha; r_pp = r_pp - u_pp*alpha
d_pp = diag_precon_pp*r_pp ;beta = DOT_PRODUCT_P(r_pp,d_pp)/up
p_pp = d_pp + p_pp * beta ; cj_converged = .TRUE.
CALL checon_par(xnew_pp,x_pp,cjtol,cj_converged,neq_pp)
IF(cj_converged.or.cjiters==cjits) EXIT
END DO conjugate_gradients
cjtot = cjtot + cjiters
!---------------------------- end of pcg process ------------------------
loads_pp = xnew_pp ; pmul_pp = .0_iwp
!----------------------- check plastic convergence -------------------
CALL checon_par(loads_pp,oldis_pp,plastol,plastic_converged,neq_pp)
IF(plasiters==1)plastic_converged=.FALSE.
IF(plastic_converged.OR.plasiters==plasits)bdylds_pp=.0_iwp
CALL gather(loads_pp,pmul_pp) ; utemp_pp = .0_iwp
!------------------------ go round the Gauss Points ---------------------
elements_4: DO iel = 1 , nels_pp
bload=.0_iwp; coord =p_g_co_pp( : , :, iel) ; eld = pmul_pp(:,iel)
gauss_points_2 : DO i = 1 , nip
CALL shape_der ( der,points,i); jac=MATMUL(der,coord)
det = determinant(jac);CALL invert(jac); deriv= MATMUL(jac,der)
CALL beemat (bee,deriv); eps=MATMUL(bee,eld)
eps = eps - evpt_pp(: ,i ,iel) ; sigma = MATMUL(dee,eps)
stress = sigma+tensor_pp(: , i, iel)
CALL invar(stress,sigm,dsbar,lode_theta)
!---------------------- check whether yield is violated ----------------
CALL mocouf( phi,c,sigm,dsbar,lode_theta,f)
IF(plastic_converged.OR.plasiters==plasits) THEN
devp=stress
ELSE
IF(f>=.0) THEN ; CALL mocouq(psi,dsbar,lode_theta,dq1,dq2,dq3)
CALL formm(stress,m1,m2,m3) ; flow=f*(m1*dq1+m2*dq2+m3*dq3)
erate=MATMUL(flow,stress) ; evp=erate*dt
evpt_pp(:,i,iel)=evpt_pp(:,i,iel)+evp ; devp=MATMUL(dee,evp)
END IF; END IF
IF(f>=.0) THEN
eload=MATMUL(TRANSPOSE(bee),devp);bload=bload+eload*det*weights(i)
END IF
IF(plastic_converged.OR.plasiters==plasits)THEN
!---------------------- update the Gauss Point stresses -----------------
tensor_pp(: , i , iel) = stress
END IF
END DO gauss_points_2
!-------------compute the total bodyloads vector -------------------------
utemp_pp(:,iel) = utemp_pp(:,iel) + bload
END DO elements_4 ; CALL scatter(bdylds_pp,utemp_pp)
IF(plastic_converged.OR.plasiters==plasits)EXIT
END DO plastic_iterations
totd_pp=totd_pp+loads_pp
