Civil Engineering Reference
In-Depth Information
READ(10,*)x_coords,y_coords
READ(10,*)dtim,nstep,theta,npri,nres,ntime
CALL sample(element,points,weights); globma=zero; gc=one
!---------create and store element and global lumped mass matrices--------
elements_1: DO iel=1,nels
CALL geom_rect(element,iel,x_coords,y_coords,coord,num,dir)
kay=zero; DO i=1,ndim; kay(i,i)=prop(i,etype(iel)); END DO
g_num(:,iel)=num; g_coord(:,num)=TRANSPOSE(coord); kc=zero; mm=zero
gauss_pts: DO i=1,nip
CALL shape_der(der,points,i); CALL shape_fun(fun,points,i)
jac=MATMUL(der,coord); det=determinant(jac); CALL invert(jac)
deriv=MATMUL(jac,der); if(type_2d=='axisymmetric')gc=MATMUL(fun,coord)
kc=kc+MATMUL(MATMUL(TRANSPOSE(deriv),kay),deriv)*det*weights(i)*gc(1)
CALL cross_product(fun,fun,ntn); mm=mm+ntn*det*weights(i)*gc(1)
END DO gauss_pts; store_kc(:,:,iel)=kc
DO i=1,nod; globma(num(i))=globma(num(i))+SUM(mm(i,:)); END DO
globma(0)=zero
END DO elements_1; CALL mesh(g_coord,g_num,12)
!-----------------------recover element A and B matrices------------------
elements_2: DO iel=1,nels
num=g_num(:,iel); kc=-store_kc(:,:,iel)*(one-theta)*dtim*pt5
mm=store_kc(:,:,iel)*theta*dtim*pt5
DO i=1,nod
mm(i,i)=mm(i,i)+globma(num(i)); kc(i,i)=kc(i,i)+globma(num(i))
END DO; CALL invert(mm); mm=MATMUL(mm,kc); store_kc(:,:,iel)=mm
END DO elements_2
!-----------------------specify initial and boundary values---------------
READ(10,*)loads(1:); loads(0)=zero; READ(10,*)fixed_freedoms
IF(fixed_freedoms/=0)THEN
ALLOCATE(node(fixed_freedoms),value(fixed_freedoms))
READ(10,*)(node(i),value(i),i=1,fixed_freedoms)
END IF
!-----------------------time stepping loop--------------------------------
WRITE(11,'(/a,i3,a)')"
Time
Pressure (node",nres,")"
WRITE(11,'(2e12.4)')0.0,loads(nres)
timesteps: DO j=1,nstep
time=j*dtim
!-----------------------first pass (1 to nels)----------------------------
elements_3: DO iel=1,nels
num=g_num(:,iel); mm=store_kc(:,:,iel)
loads(num)=MATMUL(mm,loads(num)); loads(0)=zero; loads(node)=value
END DO elements_3
!-----------------------second pass (nels to 1)---------------------------
elements_4: DO iel=nels,1,-1
num=g_num(:,iel); mm=store_kc(:,:,iel)
loads(num)=MATMUL(mm,loads(num)); loads(0)=zero; loads(node)=value
END DO elements_4
IF(nod==4.AND.j==ntime)THEN; READ(10,*)nci
CALL contour(loads,g_coord,g_num,nci,13); END IF
IF(j/npri*npri==j)WRITE(11,'(2e12.4)')time,loads(nres)
END DO timesteps
STOP
END PROGRAM p85
New dynamic real arrays:
store kc
stores element
kc
matrices
