Civil Engineering Reference
In-Depth Information
Program 12.3 Three dimensional Laplacian flow. Compare Program 7.5.
PROGRAM p123
!------------------------------------------------------------------------
! Program 7.5 three dimensional analysis of Laplace's equation
! using 8-node brick elements, preconditioned conjugate gradient solver
! only integrate one element , diagonal preconditioner diag_precon
! parallel version ; central loaded or fixed freedom ; box_bc
!------------------------------------------------------------------------
USE new_library; USE geometry_lib; USE precision; USE utility
USE mp_module ; USE timing ; USE global_variables1
USE gather_scatter6; IMPLICIT NONE
! ndof, nels, neq , ntot are now global variables - not declared
INTEGER::nxe,nye,nze,nn,nr,nip,nodof=1,nod=8, nres, is , it ,
&
i,j,k,ndim=3,iters,limit,iel,num_no,no_index_start,
&
neq_temp,nn_temp , loaded_freedoms, fixed_freedoms
REAL(iwp)::aa,bb,cc,kx,ky,kz,det,tol,up,alpha,beta,q,penalty=1.e20_iwp
CHARACTER(LEN=15)::element= 'hexahedron'; LOGICAL :: converged
!-------------------------- dynamic arrays---------------------------------
REAL(iwp),ALLOCATABLE :: points(:,:),kc(:,:),coord(:,:), weights(:), &
p_g_co_pp(:,:,:), jac(:,:), der(:,:), deriv(:,:),&
col(:,:),row(:,:),kcx(:,:),kcy(:,:),kcz(:,:),
&
diag_precon_pp(:),p_pp(:),r_pp(:),x_pp(:),
&
xnew_pp(:),u_pp(:),pmul_pp(:,:),utemp_pp(:,:),
&
d_pp(:),diag_precon_tmp(:,:),val(:),val_f(:),
&
store_pp(:),eld(:)
INTEGER, ALLOCATABLE :: rest(:,:),g(:),num(:),g_num_pp(:,:),g_g_pp(:,:), &
no(:),no_f(:),no_local_temp(:),no_local_temp_f(:),no_local(:)
!-------------------------input and initialisation-------------------------
timest(1) = elap_time( )
;
CALL find_pe_procs(numpe,npes)
IF (numpe==npes) THEN
OPEN (10,FILE='p123.dat',STATUS= 'OLD',ACTION='READ')
READ (10,*) nels,nxe,nze,nip,aa,bb,cc,kx,ky,kz, tol,limit ,
&
loaded_freedoms,fixed_freedoms
END IF
CALL bcast_inputdata_p123(numpe,npes,nels,nxe,nze,nip,aa,bb,cc,kx,ky,kz, &
tol,limit,loaded_freedoms,fixed_freedoms)
CALL calc_nels_pp ; ndof=nod*nodof ; ntot=ndof ; nye = nels/nxe/nze
neq_temp = 0; nn_temp = 0 ; nr=(nxe+1)*(nye+1)+(nxe+1)*nze+nye*nze
ALLOCATE (points(nip,ndim),coord(nod,ndim),jac(ndim,ndim),kc(ntot,ntot), &
der(ndim,nod),deriv(ndim,nod),rest(nr,nodof+1),kcx(ntot,ntot), &
g(ntot),pmul_pp(ntot,nels_pp),utemp_pp(ntot,nels_pp),col(ntot,1), &
p_g_co_pp(nod,ndim,nels_pp),g_num_pp(nod,nels_pp),weights(nip),
&
num(nod),g_g_pp(ntot,nels_pp),no(1),kcy(ntot,ntot),val(1),
&
no_local_temp(1),row(1,ntot),diag_precon_tmp(ntot,nels_pp),
&
val_f(1),eld(ntot),no_f(1),no_local_temp_f(1),kcz(ntot,ntot))
CALL box_bc8(nxe,nye,nze,rest); ielpe=iel_start; nres=nxe*(nze-1)+1
IF(loaded_freedoms>0) THEN; no = nres; val = 10._iwp; END IF
IF(fixed_freedoms>0)THEN; no_f=nres; val_f = 100._iwp ; END IF
CALL sample(element,points,weights); CALL rearrange_2(rest)
!---------- loop the elements for global cordinates etc -------------------
elements_0: DO iel = 1 , nels_pp
CALL geometry_8bxz(ielpe,nxe,nze,aa,bb,cc,coord,num)
CALL find_g4(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
