Civil Engineering Reference
In-Depth Information
CALL bcast_inputdata_p121(numpe,npes,nels,nxe,nze,nip,aa,bb,cc,
&
e,v,tol,limit)
CALL calc_nels_pp ; ndof=nod*nodof ; ntot=ndof
nn_temp = 0; neq_temp = 0 ; nye = nels/nxe/nze ; nle = nxe/5
nr = ((2*nxe+1)*(nze+1)+(nxe+1)*nze)*2+((2*nye-1)*nze+(nye-1)*nze)*2
&
+ (2*nye-1)*(nxe+1) + (nye-1)*nxe
loaded_freedoms = 3*nle*nle + 4*nle + 1
ALLOCATE (points(nip,ndim),dee(nst,nst),coord(nod,ndim),jac(ndim,ndim), &
der(ndim,nod),deriv(ndim,nod),eld_pp(ntot,nels_pp), &
bee(nst,ntot),km(ntot,ntot),eld(ntot),eps(nst),sigma(nst), &
g(ntot),pmul_pp(ntot,nels_pp),utemp_pp(ntot,nels_pp),num(nod), &
p_g_co_pp(nod,ndim,nels_pp),g_num_pp(nod,nels_pp),weights(nip),&
g_g_pp(ntot,nels_pp),no(loaded_freedoms),rest(nr,nodof+1),
&
val(loaded_freedoms),no_local_temp(loaded_freedoms),
&
diag_precon_tmp(ntot,nels_pp))
CALL cube_bc20(nxe,nye,nze,rest) ; CALL rearrange(rest)
CALL loading( nxe,nze,nle,no,val); val = - val * aa * bb / 12._iwp
CALL deemat(dee,e,v);CALL sample(element,points,weights); ielpe=iel_start
!---------- loop the elements for global cordinates etc ------------------
elements_0: DO iel = 1 , nels_pp
CALL geometry_20bxz(ielpe,nxe,nze,aa,bb,cc,coord,num)
CALL find_g3(num,g,rest) ; g_num_pp(:,iel) = num
p_g_co_pp(:,:,iel) = coord; g_g_pp(:,iel) = g ; ielpe=ielpe+1
i = MAXVAL(g); j = MAXVAL(num)
IF(i>neq_temp)neq_temp = i ; IF(j>nn_temp)nn_temp = j
END DO elements_0
neq = reduce(neq_temp); nn = reduce(nn_temp)
CALL calc_neq_pp ; CALL make_ggl(g_g_pp)
ALLOCATE(p_pp(neq_pp),r_pp(neq_pp),x_pp(neq_pp),xnew_pp(neq_pp),
&
u_pp(neq_pp),diag_precon_pp(neq_pp),d_pp(neq_pp))
r_pp = .0_iwp; p_pp = .0_iwp; x_pp = .0_iwp; xnew_pp = .0_iwp
diag_precon_pp=.0_iwp; diag_precon_tmp=.0_iwp
!------ element stiffness integration and build the preconditioner-------
iel=1;CALL geometry_20bxz(iel,nxe,nze,aa,bb,cc,coord,num);km=0.0_iwp
gauss_pts_1: 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)
km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
END DO gauss_pts_1
elements_1: DO iel = 1,nels_pp
DO k=1,ndof;diag_precon_tmp(k,iel)=diag_precon_tmp(k,iel) &
+km(k,k); END DO
END DO elements_1
CALL scatter(diag_precon_pp,diag_precon_tmp);DEALLOCATE(diag_precon_tmp)
IF(numpe==1)THEN
OPEN (11,FILE='p121.res',STATUS='REPLACE',ACTION='WRITE')
WRITE(11,'(A,I5,A)') "This job ran on ", npes ," processors"
WRITE(11,'(A)') "Global coordinates and node numbers "
DO i = 1 , nels_pp,nels_pp-1
WRITE(11,'(A,I8)')"Element ",i ; num = g_num_pp(:,i)
DO k = 1,nod;WRITE(11,'(A,I7,3E12.4)')
&
" Node ",num(k),p_g_co_pp(k,:,i); END DO
END DO
WRITE(11,'(A,3(I8,A))')"There are ",nn," nodes",nr,
&
" restrained and ", neq," equations"
WRITE(11,*) "Time after setup is
:", elap_time( ) - timest(1)
END IF
