AMSS_NCKU_source/prolongrestrict_cell.f90

@@ -1933,10 +1933,16 @@
   integer,dimension(3) :: lbc,ubc,lbf,ubf,lbp,ubp,lbpc,ubpc
 ! when if=1 -> ic=0, this is different to vertex center grid 
   real*8, dimension(-2:extc(1),-2:extc(2),-2:extc(3))   :: funcc
-  integer,dimension(3) :: cxI
-  integer :: i,j,k,ii,jj,kk
-  real*8, dimension(6,6) :: tmp2
-  real*8, dimension(6) :: tmp1
+  integer,dimension(3) :: cxI
+  integer :: i,j,k,ii,jj,kk
+  real*8, dimension(6,6) :: tmp2
+  real*8, dimension(6) :: tmp1
+  integer, dimension(extf(1)) :: cix
+  integer, dimension(extf(2)) :: ciy
+  integer, dimension(extf(3)) :: ciz
+  logical, dimension(extf(1)) :: evenx
+  logical, dimension(extf(2)) :: eveny
+  logical, dimension(extf(3)) :: evenz
 
   real*8, parameter :: C1=7.7d1/8.192d3,C2=-6.93d2/8.192d3,C3=3.465d3/4.096d3
   real*8, parameter :: C6=6.3d1/8.192d3,C5=-4.95d2/8.192d3,C4=1.155d3/4.096d3
@@ -2003,10 +2009,10 @@
   kmini=lbpc(3)-lbc(3) + 1
   kmaxi=ubpc(3)-lbc(3) + 1
 
-  if(imino.lt.1.or.jmino.lt.1.or.kmino.lt.1.or.&
-     imini.lt.1.or.jmini.lt.1.or.kmini.lt.1.or.&
-     imaxo.gt.extf(1).or.jmaxo.gt.extf(2).or.kmaxo.gt.extf(3).or.&
-     imaxi.gt.extc(1)-2.or.jmaxi.gt.extc(2)-2.or.kmaxi.gt.extc(3)-2)then
+  if(imino.lt.1.or.jmino.lt.1.or.kmino.lt.1.or.&
+     imini.lt.1.or.jmini.lt.1.or.kmini.lt.1.or.&
+     imaxo.gt.extf(1).or.jmaxo.gt.extf(2).or.kmaxo.gt.extf(3).or.&
+     imaxi.gt.extc(1)-2.or.jmaxi.gt.extc(2)-2.or.kmaxi.gt.extc(3)-2)then
           write(*,*)"error in prolongation for"
           write(*,*)"from"
           write(*,*)llbc,uubc
@@ -2017,33 +2023,40 @@
           write(*,*)"want"
           write(*,*)llbp,uubp
           write(*,*)lbp,ubp,lbpc,ubpc
-          return
-  endif
-
-  call symmetry_bd(3,extc,func,funcc,SoA)
-     
-!~~~~~~> prolongation start...
-  do k = kmino,kmaxo
-   do j = jmino,jmaxo
-    do i = imino,imaxo
-       cxI(1) = i
-       cxI(2) = j
-       cxI(3) = k
-! change to coarse level reference
-!|---*--- ---*--- ---*--- ---*--- ---*--- ---*--- ---*--- ---*---| 
-!|=======x===============x===============x===============x=======|
-       cxI = (cxI+lbf-1)/2
-! change to array index      
-       cxI = cxI - lbc + 1
-
-       if(any(cxI+3 > extc)) write(*,*)"error in prolong"
-       ii=i+lbf(1)-1
-       jj=j+lbf(2)-1
-       kk=k+lbf(3)-1
-#if 0
-       if(ii/2*2==ii)then
-         if(jj/2*2==jj)then
-           if(kk/2*2==kk)then
+          return
+  endif
+
+  do i = imino,imaxo
+     ii = i + lbf(1) - 1
+     cix(i) = ii/2 - lbc(1) + 1
+     evenx(i) = (ii/2*2 == ii)
+  enddo
+  do j = jmino,jmaxo
+     jj = j + lbf(2) - 1
+     ciy(j) = jj/2 - lbc(2) + 1
+     eveny(j) = (jj/2*2 == jj)
+  enddo
+  do k = kmino,kmaxo
+     kk = k + lbf(3) - 1
+     ciz(k) = kk/2 - lbc(3) + 1
+     evenz(k) = (kk/2*2 == kk)
+  enddo
+
+  call symmetry_bd(3,extc,func,funcc,SoA)
+     
+!~~~~~~> prolongation start...
+  do k = kmino,kmaxo
+   do j = jmino,jmaxo
+    do i = imino,imaxo
+       cxI(1) = cix(i)
+       cxI(2) = ciy(j)
+       cxI(3) = ciz(k)
+
+       if(any(cxI+3 > extc)) write(*,*)"error in prolong"
+#if 0
+       if(ii/2*2==ii)then
+         if(jj/2*2==jj)then
+           if(kk/2*2==kk)then
              tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
                    C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
                    C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
@@ -2126,11 +2139,11 @@
            endif
          endif
        endif
-#else 
-       if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
+#else 
+       if(evenz(k))then
+             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
+                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
+                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
                    C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
                    C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
                    C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
@@ -2139,21 +2152,21 @@
                    C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
                    C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
                    C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-       endif
-
-       if(jj/2*2==jj)then
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-       else
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-       endif
-
-       if(ii/2*2==ii)then
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-       else
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-       endif
+                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
+                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
+       endif
+
+       if(eveny(j))then
+             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
+       else
+             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
+       endif
+
+       if(evenx(i))then
+             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
+       else
+             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
+       endif
 #endif
     enddo
    enddo