calcbody.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	calcbody (nef, body, ipobody, ibody, xbody, coel, vel, lakon, nactdohinv)

Function/Subroutine Documentation

calcbody()

subroutine calcbody	(	integer	nef,
		real*8, dimension(0:3,*)	body,
		integer, dimension(2,*)	ipobody,
		integer, dimension(3,*)	ibody,
		real*8, dimension(7,*)	xbody,
		real*8, dimension(3,*)	coel,
		real*8, dimension(nef,0:7)	vel,
		character*8, dimension(*)	lakon,
		integer, dimension(*)	nactdohinv
	)

!
!     calculation of the actual body force in each element. The body
!     force is the sum of gravity and centrifugal/Coriolis forces
!
      implicit none
!
      character*8 lakon(*)
!
      integer i,j,nef,index,ipobody(2,*),ibody(3,*),nactdohinv(*)
!
      real*8 om,body(0:3,*),p1(3),p2(3),xbody(7,*),omcor,q(3),coel(3,*),
     &  vel(nef,0:7),const,corio(3)
!
      do i=1,nef
         om=0.d0
         do j=1,3
            body(j,i)=0.d0
         enddo
!
         index=nactdohinv(i)
!
         do
            j=ipobody(1,index)
            if(j.eq.0) exit
            if(ibody(1,j).eq.1) then
               om=xbody(1,j)
               p1(1)=xbody(2,j)
               p1(2)=xbody(3,j)
               p1(3)=xbody(4,j)
               p2(1)=xbody(5,j)
               p2(2)=xbody(6,j)
               p2(3)=xbody(7,j)
!     
!     assigning gravity forces
!     
            elseif(ibody(1,j).eq.2) then
               body(1,i)=body(1,i)+xbody(1,j)*xbody(2,j)
               body(2,i)=body(2,i)+xbody(1,j)*xbody(3,j)
               body(3,i)=body(3,i)+xbody(1,j)*xbody(4,j)
            endif
            index=ipobody(2,index)
            if(index.eq.0) exit
         enddo
!
!        adding the centrifugal/Coriolis force (if any) to the
!        gravity loading
!
         if(om.gt.0.d0) then
            omcor=2.d0*dsqrt(om)
            do j=1,3
               q(j)=coel(j,i)-p1(j)
            enddo
            const=q(1)*p2(1)+q(2)*p2(2)+q(3)*p2(3)
!     
!     Coriolis forces
!     
            corio(1)=vel(i,2)*p2(3)-vel(i,3)*p2(2)
            corio(2)=vel(i,3)*p2(1)-vel(i,1)*p2(3)
            corio(3)=vel(i,1)*p2(2)-vel(i,2)*p2(1)
!     
!     inclusion of the centrifugal force into the body force
!     
            do j=1,3
               body(j,i)=body(j,i)+(q(j)-const*p2(j))*om+
     &              corio(j)*omcor
            enddo
         endif
!
      enddo
!     
      return