Demo entry 6659135

lpw

   

Submitted by lpw1 on Nov 09, 2017 at 11:04
Language: Fortran. Code size: 10.8 kB.

    program My_Code

!!!!!!!!!!!!!!!!!!!!!!!!!!!  Variable  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none

	  integer:: nels,nn,nsource,iel,i,j,k,bn,nd,icd,ii,jj,jel,nband,mi,ma,          &
	            idiff,ncoef,i1,i2,i3,neq,ival,kk

   	  real(8):: pi,alpha0,x1,x2,x3,y1,y2,y3,area,length,sm1,sm2,sm3,dis1,dis2,      &
   	            rbound,xalpha,temp1,temp2,absorb0,absorb1,diffuse0,diffuse1,        &
   	            loads_intensity=1.0d0

      real(8),allocatable::g_coord(:,:),coord(:,:),loads(:),absorb(:),diffuse(:),   &
                           factord(:),aa(:),bb(:),cc(:),yl1(:),yl2(:),kp_de(:),     &
                           bp_de(:),loads0(:)

	  integer,allocatable::g_num(:,:),num(:),bnode(:,:),mm(:),detector(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!! Initial !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


      open (11, file = 'coord.dat' ,  status = 'old' , action = 'read' )
      open (12, file = 'num.dat' ,    status = 'old' , action = 'read' )
      open (13, file = 'boundary.dat',status = 'old' , action = 'read' )                                                
      open (14, file = 'detector.dat',status = 'old' , action = 'read' )      
      open (15, file = 'mesh1924.dat',status = 'old' , action = 'read' )

      open (31, file = 'exact_mua.dat' ,  status='replace',action = 'write')                                                
      open (32, file = 'exact_mus.dat' ,  status='replace',action = 'write')                                                

      open (33, file = 'result_phi.dat' ,  status='replace',action = 'write')                                                
      open (34, file = 'result_phi_log.dat' ,  status='replace',action = 'write')                                                
      open (35, file = 'result_boundary.dat' ,  status='replace',action = 'write')                                                

!	  read(10,*) nn,nels,nsource,bn,nd                                
!	  read(10,*) absorb0,diffuse0,alpha0     
!	  read(10,*) absorb1,diffuse1     
        
      nn=1924
      nels=3686
      nsource=2
      bn=64
      nd=64
      absorb0=0.025
      absorb1=0.05
      diffuse0=2.0
      diffuse1=5.0
      alpha0=0.476
    	                          
      pi=dacos(-1.0d0)

      write(*,*) "FEA preprocessing...."

      allocate(g_coord(2,nn),coord(3,2),g_num(3,nels),num(3),bnode(2,bn),absorb(nn),   &
               diffuse(nn),loads(nn),mm(nn),aa(3),bb(3),cc(3),yl1(nd),yl2(nd),         &
               factord(nn),detector(nd),loads0(nn))     

	  do i=1,bn
	     read(13,*) bnode(1,i),bnode(2,i),bnode(3,i)
	  end do

	  do i=1,nd
	     read(14,*) detector(i)
	  end do

!!!! global coordinates information:

      do i=1,nn
	     read(11,*) g_coord(1,i), g_coord(2,i)
      end do
	  g_coord=g_coord*15.0d0/25.0d0
	
!!!! global elements information:

      do iel=1,nels
         read(12,*) g_num(1,iel),g_num(2,iel),g_num(3,iel)
      end do

!!!! optical parameters:

      do i=1,nn
	     absorb(i)=absorb0
		 diffuse(i)=diffuse0
         x1=g_coord(1,i);y1=g_coord(2,i)
		 dis1=dsqrt(((x1-0.0d0)**2.+(y1-7.5d0)**2.))
		 if(dis1<=2.0d0) then
		    absorb(i)=absorb1
		    diffuse(i)=diffuse1
         endif
	     factord(i)=1.0d0/(3.0d0*(absorb(i)+diffuse(i)))
      end do

      write(31,*)" TITLE = ""Example: 2D Circle DOT Finite-Element Data"" "
	  write(31,*)" VARIABLES = ""X"", ""Y"", ""P"" "
	  write(31,*)" ZONE N=634, E=1202, F=FEPOINT, ET=TRIANGLE"
      do i=1,nn
	     write(31,*) g_coord(1,i),g_coord(2,i),absorb(i)
	  end do
      do iel=1,nels
	     write(31,*) g_num(1,iel),g_num(2,iel),g_num(3,iel)
	  end do
      write(32,*)" TITLE = ""Example: 2D Circle DOT Finite-Element Data"" "
	  write(32,*)" VARIABLES = ""X"", ""Y"", ""P"" "
	  write(32,*)" ZONE N=634, E=1202, F=FEPOINT, ET=TRIANGLE"
      do i=1,nn
	     write(32,*) g_coord(1,i),g_coord(2,i),diffuse(i)
	  end do
      do iel=1,nels
	     write(32,*) g_num(1,iel),g_num(2,iel),g_num(3,iel)
	  end do
      
!!! calculate the half band width

	  nband=0
      do iel=1,nels
         num=g_num(:,iel)			  			                       
	     mi=nn
		 ma=0
		 do i=1,3
		    ii=num(i)
			if(ii<mi) mi=ii
			if(ii>ma) ma=ii
		 end do
		 idiff=ma-mi
		 if(idiff>nband) nband=idiff
	  end do
	  write(*,*) "Half bandwidth is", nband

!!! use 2d storage method

      neq=nn
	  ncoef=neq*(nband+1)
      allocate(kp_de(ncoef),bp_de(ncoef))

      write(*,'(2(a,i5))')                                                    &
              "There are ",neq,"  equations and the half-bandwidth is ",nband

      write(*,*) "Now start calculating global stiffness matrix...."
                 
!!!! Calculating the global stiffness matrix: 

      kp_de=0.0d0;bp_de=0.0d0
      do iel=1,nels
         num=g_num(:,iel)			  			                       
         coord=transpose(g_coord(:,num)) 
         x1=coord(1,1); y1=coord(1,2)
         x2=coord(2,1); y2=coord(2,2)
         x3=coord(3,1); y3=coord(3,2)      
         aa(1)=x2*y3-x3*y2; aa(2)=x3*y1-x1*y3; aa(3)=x1*y2-x2*y1
         bb(1)=y2-y3; bb(2)=y3-y1; bb(3)=y1-y2
         cc(1)=-x2+x3; cc(2)=-x3+x1; cc(3)=-x1+x2
         area=(aa(1)+aa(2)+aa(3))*0.5d0
	     if(area<0.d0) write(*,*) "area error"
		 do i=1,3
		    ii=num(i)
		    do j=1,3
			   jj=num(j)
			   if(i==j) then
			      sm3=area/6.0d0
			   else
			      sm3=area/12.0d0
			   endif
			   icd=jj-ii+1
			   if(icd-1>=0) then
			      ival=neq*(icd-1)+ii
				  do k=1,3
				     if(i==j.and.j==k) then
					    sm1=area/10.0d0
					 elseif(i/=j.and.j/=k.and.k/=i) then
					    sm1=area/60.0d0
					 else
					    sm1=area/30.0d0
					 endif
				     kk=num(k)
				     kp_de(ival)=kp_de(ival)+(bb(i)*bb(j)+cc(i)*cc(j))/         &
				                 (12.0d0*area)*factord(kk)+sm1*absorb(kk)
				     bp_de(ival)=bp_de(ival)+sm3/3.0d0
				  end do
               endif
			end do
		 end do
	  end do

	  do jel=1,bn
	     i1=bnode(1,jel)
	     i2=bnode(2,jel)
         rbound=dsqrt(g_coord(1,i1)**2.d0+g_coord(2,i1)**2.d0)
	     length=dsqrt((g_coord(1,i1)-g_coord(1,i2))**2.d0+                      &
	                  (g_coord(2,i1)-g_coord(2,i2))**2.d0)
	     do i=1,2
		    ii=bnode(i,jel)
		    do j=1,2
			   jj=bnode(j,jel)
			   if(i==j) then
			      sm2=length/3.d0
			   else
			      sm2=length/6.d0
			   endif
			   icd=jj-ii+1
			   if(icd-1>=0) then
			      ival=neq*(icd-1)+ii
				  kp_de(ival)=kp_de(ival)+sm2*alpha0
			   endif
			end do
		 end do
	  end do

!!! factorise left hand side
    
	  call banred(kp_de,neq,ncoef)

!!!! source type:

      loads0=0.0d0
      loads0(nsource)=loads_intensity
	  call linmul(bp_de,loads0,loads,ncoef,neq)

!!!!!!!!!!!!!!!!!!!!!! Calculating the matrix equation !!!!!!!!!!!!!!!!!!!!!!!!!!

      write(*,*) "Start to solve the diffuse equations..."
      call bacsub(kp_de,loads,ncoef,neq)

      write(*,*) "Start to output the DE results..."

      write(33,*)" TITLE = ""Example: 2D Circle DOT Finite-Element Data"" "
	  write(33,*)" VARIABLES = ""X"", ""Y"", ""P"" "
	  write(33,*)" ZONE N=1924, E=3686, F=FEPOINT, ET=TRIANGLE"
      do i=1,nn
	     write(33,*) g_coord(1,i),g_coord(2,i),loads(i)
	  end do
      do iel=1,nels
	     write(33,*) g_num(1,iel),g_num(2,iel),g_num(3,iel)
	  end do

      write(34,*)" TITLE = ""Example: 2D Circle DOT Finite-Element Data"" "
	  write(34,*)" VARIABLES = ""X"", ""Y"", ""P"" "
	  write(34,*)" ZONE N=634, E=1202, F=FEPOINT, ET=TRIANGLE"
      do i=1,nn
	     write(34,*) g_coord(1,i),g_coord(2,i),log10(loads(i))
	  end do
      do iel=1,nels
	     write(34,*) g_num(1,iel),g_num(2,iel),g_num(3,iel)
	  end do

!!! detector results output:  

     YL1=0.d0;YL2=0.d0
     do i=1,nd
	      x1=g_coord(1,detector(i))
	      y1=g_coord(2,detector(i))
		  xalpha=atan2(y1,x1)
          if(xalpha>0.) YL1(i)=xalpha*180.d0/pi
		  if(xalpha<0.) YL1(i)=xalpha*180.d0/pi+360.d0
		  if((x1==0.).and.(y1>0.)) YL1(i)=90.d0
		  if((x1==0.).and.(y1<0.)) YL1(i)=270.d0
          if((y1==0.).and.(x1>0.)) YL1(i)=0.d0
		  if((y1==0.).and.(x1<0.)) YL1(i)=180.d0
	      YL2(i)=loads(detector(i))
	 end do			    
     do i=nd-1,1,-1
	    do j=1,i
        if(YL1(j)>YL1(j+1)) then
	       temp1=YL1(j)
		   YL1(j)=YL1(j+1)
		   YL1(j+1)=temp1
	       temp2=YL2(j)
		   YL2(j)=YL2(j+1)
		   YL2(j+1)=temp2
        end if
	    end do
	 end do
	 do i=1,nd
	    write(35,*) YL1(i),YL2(i)
	 end do

 end program My_Code

!!!!!!!!!!!! algorithm for solving the matrix equation !!!!!!!!!!!!!!!!!
         
subroutine banred(bk,n,ncoef)
! gaussian reduction on a vector stored as an upper triangle
implicit none
integer::n,ncoef
integer::i,il1,kbl,j,ij,nkb,m,ni,nj,iw 
real(8)::bk(ncoef)
real(8)::sum

iw = ubound(bk,1)/n-1
       do i=2,n
          il1=i-1;kbl=il1+iw+1
          if(kbl-n>0)kbl=n
          do j=i,kbl
             ij=(j-i)*n+i;sum=bk(ij);nkb=j-iw
             if(nkb<=0)nkb=1
             if(nkb-il1<=0)then
                do m=nkb,il1
                   ni=(i-m)*n+m ; nj=(j-m)*n+m
                   sum=sum-bk(ni)*bk(nj)/bk(m) 
                end do
             end if
             bk(ij)=sum
           end do
       end do
return
end subroutine banred   


subroutine bacsub(bk,loads,ncoef,neq)
! performs the complete gaussian backsubstitution
implicit none
integer::nkb,k,i,jn,jj,i1,iw,ncoef,neq
real(8)::sum
real(8)::bk(ncoef)
real(8)::loads(neq)
iw=ubound(bk,1)/neq-1
loads(1)=loads(1)/bk(1)
   do i=2,neq
      sum=loads(i);i1=i-1 ; nkb=i-iw
      if(nkb<=0)nkb=1
      do k=nkb,i1
         jn=(i-k)*neq+k;sum=sum-bk(jn)*loads(k)
      end do
      loads(i)=sum/bk(i)
   end do
   do jj=2,neq
      i=neq-jj+1;sum=.0;i1=i+1;nkb=i+iw
      if(nkb-neq>0)nkb=neq
      do k=i1,nkb
           jn=(k-i)*neq+i  ; sum=sum+bk(jn)*loads(k)
      end do
      loads(i)=loads(i)-sum/bk(i)
   end do
return
end subroutine bacsub                                                       

subroutine linmul(bk,disps,loads,ncoef,neq)
! matrix-vector multiply for symmetric matrix bk
! stored in upper triangular form as a vector

implicit none

integer::i,j,iw,ncoef,neq
real(8)::x 
real(8)::bk(ncoef),disps(neq),loads(neq)

iw=ubound(bk,1)/neq-1 
    do i=1,neq
       x=0.0d0
       do j=1,iw+1
          if(i+j<=neq+1)    x=x+bk(neq*(j-1)+i)*disps(i+j-1)
       end do
       do j=2,iw+1
          if(i-j+1>=1)    x=x+bk((neq-1)*(j-1)+i)*disps(i-j+1)
       end do    
       loads(i)=x
    end do
  return
end subroutine linmul

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).