cfl3d的重叠网格处理程序maggie真是奇芭，难怪自己都不发展了。

 

1 它是基于cellcenter来挖洞和找donnor的，因此它并非基于原始网格处理，而是将原始的顶点网格（vc）的cellcenter连接为一个新的网格（cc），因此新的cc网格是比原始vc网格缩水了一圈。这一点在如下代码中体现

      subroutine setup 
c 
c 此处省略若干行
c 
c     read input data and initialize
c 
      call initia 
c
      do 11 mesh = 1,nmesh
c 注意 m?max是cellsize ， +1后成为point size
      jdg = mjmax(mesh)+1
      kdg = mkmax(mesh)+1
      ldg = mlmax(mesh)+1
c 
c     read in grid
c
      if(ip3dgrd.eq.0) then
c       (cfl3d format)
        call rgrid( mesh,x,y,z,jdg,kdg,ldg )
      else
c       (plot3d/tlns3d format)
        call rp3d( mesh,x,y,z,jdg,kdg,ldg,ialph,nmesh )
      end if
c 
      npnts = jdg*kdg*ldg 
c 
c 此处省略若干行

c     check for branch cuts in grid
c
      if(mesh.eq.1) then
        write(6,*)
        write(6,*) '  *** beginning checks to identify branch',
     .  ' cuts ***'
      end if
      write(6,*)
      call branch(mesh,jdg,kdg,ldg,jimage,kimage,limage,
     .                 x,y,z)
c
c     transform grid to cell centered grid. 
c 
      call cellcen(x,y,z,jdg,kdg,ldg)
c 
c 此处省略若干行

关键在于cellcen的调用

      subroutine cellcen(x,y,z,jd,kd,ld)
c
c***************************************************************** 
c     Purpose: find the cell centers of the grid 
c***************************************************************** 
c
      include 'mag1.h'
c 
      dimension x(jd,kd,ld),y(jd,kd,ld),z(jd,kd,ld)
      common/temp1/   x1(jdim,kdim,ldim),y1(jdim,kdim,ldim),
     .                z1(jdim,kdim,ldim) 
c 
      jmax = jd
      kmax = kd
      lmax = ld
c 
      do 10 l=1,lmax-1 
      do 20 k=1,kmax-1 
      do 30 j=1,jmax-1
      x1(j,k,l)=(x(j,k,l)+x(j+1,k,l)+x(j,k+1,l)+x(j+1,k+1,l)+
     .           x(j,k,l+1)+x(j,k+1,l+1)+x(j+1,k,l+1)+ 
     .           x(j+1,k+1,l+1))/8.0 
      y1(j,k,l)=(y(j,k,l)+y(j+1,k,l)+y(j,k+1,l)+y(j+1,k+1,l)+
     .           y(j,k,l+1)+y(j,k+1,l+1)+y(j+1,k,l+1)+ 
     .           y(j+1,k+1,l+1))/8.0 
      z1(j,k,l)=(z(j,k,l)+z(j+1,k,l)+z(j,k+1,l)+z(j+1,k+1,l)+
     .           z(j,k,l+1)+z(j,k+1,l+1)+z(j+1,k,l+1)+ 
     .           z(j+1,k+1,l+1))/8.0 
30    continue
20    continue 
10    continue 
c
      do 140 l=1,lmax-1
      do 150 k=1,kmax-1 
      do 160 j=1,jmax-1
      x(j,k,l) = x1(j,k,l) 
      y(j,k,l) = y1(j,k,l) 
      z(j,k,l) = z1(j,k,l) 
160   continue 
150   continue
140   continue 
c
      return 
      end

hole是整体驱动

03658       subroutine hole
03659 c 
03660 c***********************************************************************
03661 c     Purpose: construct hole in mesh m introduced by the presence of 
03662 c     mesh m1 and compute the interpolation stencils for mesh m fringe 
03663 c     points surrounding the hole. 
03664 c 
03665 c     iblank = 0 for points interior to hole
03666 c     iblank < 0 for points surrounding hole (fringe points)
03667 c     iblank = 1 for regular field points
03668 c
03669 c     if ihplt(nh) > 0  plot3d files for the initial hole and normals are
03670 c     output. hole (grid file) is in hole_grd.nh, components of normal vector 
03671 c     (q file) are in hole_nrm.nh. plot hole points with "points" option
03672 c     in plot3d. plot normals using vector velocity. nh is the hole number.
03673 c***********************************************************************

intp是基于块操作的搜索

05174       subroutine intpt( itr,jimage,kimage,limage,jd,kd,
05175      .                   ld,m1,xm1,ym1,zm1,jd1,kd1,ld1,m,i1)
05176 c
05177 c*********************************************************************** 
05178 c     Purpose: find interpolation coefficients for interpolation to fringe 
05179 c     or boundary points of mesh m from mesh m1
05180 c***********************************************************************

search2给定点寻dono

07167       subroutine search2(i,m,m1,jd1,kd1,ld1,xm1,ym1,zm1,
07168      .                   jimage,kimage,limage,xp,yp,zp,iok,icall)     
07169 c
07170 c******************************************************************
07171 c     Purpose: search over mesh m1 to find the cell surrounding the 
07172 c     mesh m point xp,yp,zp and compute the required interpolation data
07173 c******************************************************************

其中对于给定点寻donor的关键在于topol，xe和xe2只是由最近cc点编号构造六面体进行Newton迭代找到局部坐标

07424       subroutine topol(m1,jdim,kdim,ldim,x,y,z,jimage,kimage,limage,
07425      .                 xp,yp,zp,xie,eta,zeta,jp,kp,lp,itmax,limit,iok,
07426      .                 js,je,ks,ke,ls,le,idsmin)
07427 c
07428 c***********************************************************************
07429 c     Purpose: search mesh m1 for  the mesh m point with 
07430 c     coordinates xp,yp,zp; determine corresponding xie,eta,zeta
07431 c     在m1网格中寻找 网格m的点(xp,yp,zp)  计算得到 (xie,eta,zeta)
07432 c***********************************************************************

 

在寻找点donor时，先在粗网格上用最小距离比较得到最近点

cc网格没有将vc网格的边界保留，因此当在边界附近插值时必然遇到问题，这时插值会被退化为零阶插值或者外推插值。

如下是hole函数片段

c       initialize block pointer array nblkpt
c       初始化数组nblkpt
c       ...nblkpt(i) = 0 if interp. stencil for fringe point i not yet found
c                        如果点i的插值模板还未被发现
c       ...nblkpt(i) = -(m1+mdim)  if interp. stencil for fringe point i 
c                                  has been found in mesh m1
c                                  如果点i的插值模板在网格m1中被发现
        do 112 i=i1,itotal
        nblkpt(i) = 0 
112     continue
c
c       search over all meshes in search list to find interpolation
c       stencils for fringe points of mesh m associated with hole n
c       遍历所有网格来搜索hole(n)的每个fringe点的插值模板
        do 12 nn = 1,nserch
c       网格m的第n个hole的fringe点搜索目的网格集合是lhole(m,n,:)
c       而nserch=mhole(m,n) 是网格m的第n个hole的fringe点搜索目的网格的个数
        m1  = lhole( m,n,nn ) 
c       获取m1网格的cellsize
        jd1 = mjmax(m1)
        kd1 = mkmax(m1)
        ld1 = mlmax(m1)
c       get a copy of cell-center grid for mesh m1 from the file 
c       temp_cen.m, along with branch-cut arrays
        call getgrd( m1,x,y,z,jimage,kimage,limage,jd1,kd1,ld1 ) 
c       find interpolation coefficients for interpolation to fringe 
c       or boundary points of mesh m from mesh m1        
        call intpt( -1,jimage,kimage,limage,jd,kd,ld, 
     .  m1,x,y,z,jd1,kd1,ld1,   m,i1 ) 
c       reject stencils that contain hole/fringe pts.
          call getibl( m1,iblank,jd1,kd1,ld1 )
c         对网格m的第n个hole的所有itotal个fringe点依次循环
          do 1114 i=i1,itotal
          if(nblkpt(i) .eq. -(m1+mdim)) then
c            这是 intpt 函数搜索到的donnor点         
             j     = jn(i)
             k     = kn(i)
             l     = ln(i)
c            根据donor点创建插值模板 (j:jp1, k:kp1, l:lp1)
             jp1   = min( j+1,jd1 )
             kp1   = min( k+1,kd1 )
             lp1   = min( l+1,ld1 )
c            得到插值模板的八个顶点的1D编号
             ii1 = j   + (k-1)*jd1   + (l-1)*jd1*kd1
             ii2 = jp1 + (k-1)*jd1   + (l-1)*jd1*kd1
             ii3 = jp1 + (kp1-1)*jd1 + (l-1)*jd1*kd1
             ii4 = j   + (kp1-1)*jd1 + (l-1)*jd1*kd1
             ii5 = j   + (k-1)*jd1   + (lp1-1)*jd1*kd1
             ii6 = jp1 + (k-1)*jd1   + (lp1-1)*jd1*kd1
             ii7 = jp1 + (kp1-1)*jd1 + (lp1-1)*jd1*kd1
             ii8 = j   + (kp1-1)*jd1 + (lp1-1)*jd1*kd1
c            检查插值模板是否含有fringe/solid点，如果有则拒绝，保持 nblkpt(i) = 0
             if(iblank(ii1) .le. 0 .or.
     .          iblank(ii2) .le. 0 .or.
     .          iblank(ii3) .le. 0 .or.
     .          iblank(ii4) .le. 0 .or.
     .          iblank(ii5) .le. 0 .or.
     .          iblank(ii6) .le. 0 .or.
     .          iblank(ii7) .le. 0 .or.
     .          iblank(ii8) .le. 0) then
                  nblkpt(i) = 0
                  write(88,*)'in hole, rejecting pt. i = ',i,
     .            ' mesh ',m
             end if
          end if
 1114     continue
12      continue
c 
c       check to see if interpolation stencils were not found for any 
c       mesh m fringe points. for such "orphan" points, use either the nearest
c       point from one of the meshes in the search list, with xi=eta=zeta=0 
c       (zeroth order interpolation), or, use extrapolation from the nearest 
c       point.  extrapolation is distinguished from interpolation in that an
c       interpolation stencil has 0 < xie,eta,zeta < 1 (to within a tolerance
c       epsc) while for extrapolation, either xie, eta or zeta is < 0 or > 1
c       use of the nearest neighbor or extrapolation should occur only from 
c       boundaries of mesh m1; otherwise search routine has failed in some way
c       统计插值模板无法构造的fringe点（orphan）个数
c       orphan点采用0阶插值  或者外插 
c       0阶插和外插仅当 在m1网格的边界处出现
        notok = 0
        do 113 i=i1,itotal
        if(nblkpt(i).eq.0) notok = notok + 1
113     continue
c
        if(notok.gt.0) then
          iflg = 0
          call orphan(m,n,nserch,i1,iorph,iflg)
        end if    

 

 

cc网格只考虑了单块网格的C/O构型中出现的cut面，而不考虑对接面，因此在对接面附近插值，也会退化为零阶插值或者外推插值，而无法利用对接网格的信息。

如下是topol函数片段

07462 c     find nearest point (with indicies jp,kp,lp) to xp,yp,zp; start with
07463 c     minimum distance search if idsmin > 0
07464 c     idsmin>0   找距离(xp,yp,zp)最近的点(jp,kp,lp)
07465 c     search only on coarser version of  mesh 
07466 c     在间隔为4的粗网格上搜索
07467       jskip = 4
07468       kskip = 4
07469       lskip = 4
07470       if(je-js.lt.4) jskip = 2
07471       if(ke-ks.lt.4) kskip = 2
07472       if(le-ls.lt.4) lskip = 2
07473 c     通过最小距离比较得到 (jp,kp,lp)
07474       if(idsmin.gt.0) call dsmin(jdim,kdim,ldim,x,y,z,
07475      .xp,yp,zp,jp,kp,lp,js,je,ks,ke,ls,le,jskip,kskip,lskip,dmin)
07476 c
07477 c     keep within bounds 1 < j,k,l < j/k/ldim-1 for interp/extrap stencils
07478 c     m1网格size是 1 : ?dim-1
07479 c     说明最近点的范围是  1 : ?dim-1
07480 c     在xe函数中，六面体由 (jp:jp+1,kp:kp+1, lp:lp+1)构造
07481 c     因此六面体是由cellcenter构成。如果frigne在与face邻接的cell中且靠近face，
07482 c     那么导致外插
07483       jp = min0( jp , jdim-1 )
07484       kp = min0( kp , kdim-1 )
07485       lp = min0( lp , ldim-1 )
07486       jp = max0( 1 , jp )
07487       kp = max0( 1 , kp )
07488       lp = max0( 1 , lp )
07489 c
07490       do 5555 intern=1,itmax
07491 c
07492       jfroz(intern) = jp
07493       kfroz(intern) = kp
07494       lfroz(intern) = lp
07495 c
07496       call trace(3,intern,idum2,idum3,idum4,dum1,dum2,dum3)
07497 c
07498 c     find local xie,eta,zeta via Newton iteraton in current target 
07499 c     cell jp,kp,lp
07500 c
07501       call trace(4,jp,kp,lp,m1,dum1,dum2,dum3)
07502 c     这里就是搜索donor
07503       call xe(jdim,kdim,ldim,x,y,z,jp,kp,lp,xp,yp,zp, 
07504      .              xie,eta,zeta,imiss)
07505 c
07506       call trace(5,idum1,idum2,idum3,idum4,xie,eta,zeta)
07507 c
07508 c     current target cell correct if imiss = 0
07509 c
07510       if (imiss.eq.0) go to 5556
07511 c
07512 c     update current guess for target cell based on result of Newton
07513 c     iteration, with max allowable change set by limit
07514 c
07515       if (xie.ge.0) jinc = abs(xie)
07516       if (xie.lt.0) jinc = abs(xie-1)
07517       if (eta.ge.0) kinc = abs(eta)
07518       if (eta.lt.0) kinc = abs(eta-1)
07519       if (zeta.ge.0) linc = abs(zeta)
07520       if (zeta.lt.0) linc = abs(zeta-1)
07521 c
07522       jinc = min0( jinc , limit )
07523       kinc = min0( kinc , limit )
07524       linc = min0( linc , limit )
07525 
07526 c     根据Newton迭代得到的xie eta zeta与区间[0,1]之间的关系
07527 c     确定下一个搜索点
07528       if (xie.gt.1.0) then
07529          jp = jp + jinc  
07530       else if (xie.lt.0.) then
07531          jp = jp - jinc  
07532       end if
07533       if (eta.gt.1.0) then
07534          kp = kp + kinc
07535       else if (eta.lt.0.) then
07536          kp = kp - kinc 
07537       end if
07538       if (zeta.gt.1.0) then
07539          lp = lp + linc
07540       else if (zeta.lt.0.) then
07541          lp = lp - linc 
07542       end if
07543 c
07544       xieg  = float(jp) 
07545       etag  = float(kp) 
07546       zetag = float(lp) 
07547 c
07548 c     keep within bounds of mesh m1
07549 c
07550       jp = min0( jp , jdim-1 )
07551       kp = min0( kp , kdim-1 )
07552       lp = min0( lp , ldim-1 )
07553       jp = max0( 1 , jp )
07554       kp = max0( 1 , kp )
07555       lp = max0( 1 , lp )
07556 c
07557 c     account for any branch cuts
07558 c     brach and image  是什么  
07559       jpc = jp
07560       kpc = kp
07561       lpc = lp
07562 c     xieg etag zetag 超出当前网格，需要修正
07563 c     image是针对C/O网格的cut来修正，但是无法处理多块对接情况
07564 c     因此新的jpc kpc lpc还是在m1网格内
07565       if(xieg .lt. 1. .or. xieg .gt.jdim-1) then
07566         jpc = jimage(jp,kp,lp)
07567         kpc = kimage(jp,kp,lp)
07568         lpc = limage(jp,kp,lp)
07569       end if
07570       if(etag .lt. 1. .or. etag .gt.kdim-1) then
07571         jpc = jimage(jp,kp,lp)
07572         kpc = kimage(jp,kp,lp)
07573         lpc = limage(jp,kp,lp)
07574       end if
07575       if(zetag.lt. 1. .or. zetag.gt.ldim-1) then
07576         jpc = jimage(jp,kp,lp)
07577         kpc = kimage(jp,kp,lp)
07578         lpc = limage(jp,kp,lp)
07579       end if
07580 c
07581       if(jpc.ne.jp .or. kpc.ne.kp .or. lpc.ne.lp) then
07582         call trace(23,m1,jpc,kpc,lpc,dum1,dum2,dum3)
07583         call trace(24,idum1,jp,kp,lp,dum1,dum2,dum3)
07584       end if
07585 c     更新jkl值
07586       jp  = jpc
07587       kp  = kpc
07588       lp  = lpc

 

因此可以总结说，maggie搜索的donor是尺寸为(1:jdim-1, 1:kdim-1, 1:ldim-1)的cc网格中的单元  (jp:jp+1, kp:kp+1, lp:lp+1)

 

 

3 在插值时，cfl3d首先构造一个基于cellcenter的网格来存储cc变量，但是这个网格却包含了边界信息，即它构造了cellcenter切面在边界面上的交点（即有关的face center和edge center，以及cornor）上的变量值。这个网格对应的数组的size是(1:jdim+1, 1:kdim+1, 1:ldim+1)，其中cc对应的编号是 (2:jdim, 2:kdim, 2:ldim)。那么maggie搜索到的cc单元ijk编号需要全部+1才能cfl3d的数据一致。具体见cfl3d的intrbc

      call augmntq(q,jdim,kdim,idim,nbl,ldim,qj0,qk0,qi0,qq,
     .             bcj,bck,bci,nou,bou,nbuf,ibufdim,icorr)
c     3D 情况
      if (i2d.eq.0) then
         do 12 n=1,ldim
         do 12 l=lsta,lend
c
c        set up interpolation coefficients 
         s1 = qq(jjig(l)+1, kkig(l)+1, iiig(l)+1, n)
         s2 = qq(jjig(l)+1, kkig(l)+1, iiig(l)+2, n)
         s3 = qq(jjig(l)+2, kkig(l)+1, iiig(l)+2, n)
         s4 = qq(jjig(l)+2, kkig(l)+1, iiig(l)+1, n)
         s5 = qq(jjig(l)+1, kkig(l)+2, iiig(l)+1, n)
         s6 = qq(jjig(l)+1, kkig(l)+2, iiig(l)+2, n)
         s7 = qq(jjig(l)+2, kkig(l)+2, iiig(l)+2, n)
         s8 = qq(jjig(l)+2, kkig(l)+2, iiig(l)+1, n)
c 三线性插值
         a1 =  s1 
         a2 = -s1+s2 
         a3 = -s1+s4 
         a4 = -s1+s5 
         a5 =  s1-s2+s3-s4
         a6 =  s1-s2-s5+s6
         a7 =  s1-s4-s5+s8
         a8 = -s1+s2-s3+s4+s5-s6+s7-s8 
c
c        interpolate and store in qb array 
c
         qb(l,n,iset) = a1 + a2*dxintg(l) 
     .                + a3*dyintg(l)
     .                + a4*dzintg(l)
     .                + a5*dxintg(l)*dyintg(l)
     .                + a6*dxintg(l)*dzintg(l)
     .                + a7*dyintg(l)*dzintg(l)
     .                + a8*dxintg(l)*dyintg(l)*dzintg(l) 
   12    continue
      else

 

可见其中jjig，kkig和llig全部被+1

 

最后参考augmntq可以看到它是如何构造基于cc扩展网格的。从中可以看出，虽然网格扩展到了边界上，但是边界上的点的变量是按照等距假设来外推的。事实上，边界上的单元只有半个单位。

      subroutine augmntq(q,jdim,kdim,idim,nbl,ldim,qj0,qk0,qi0,qq,
     .                   bcj,bck,bci,nou,bou,nbuf,ibufdim,ighost)
c
c     $Id: augmntq.F,v 1.2 2001/05/25 20:00:01 biedron Exp $
c
c***********************************************************************
c     Purpose:  Create an "augmented" q array (qq) of dimensions 
c     jdim+1 x kdim+1 x idim+1, in which the min/max indicies
c     contain cell-face center data and the interior points contain
c     cell-center data at the standard cfl3d cell center locations.
c***********************************************************************
c
#   ifdef CMPLX
      implicit complex(a-h,o-z)
#   endif
c
      character*120 bou(ibufdim,nbuf)
c
      dimension nou(nbuf)
      dimension q(jdim,kdim,idim,ldim),qi0(jdim,kdim,ldim,4),
     .          qj0(kdim,idim-1,ldim,4),qk0(jdim,idim-1,ldim,4)
      dimension qq(jdim+1,kdim+1,idim+1,ldim)
      dimension bcj(kdim,idim-1,2),bck(jdim,idim-1,2),bci(jdim,kdim,2)
c
      common /twod/ i2d
c
c     load interior q data into qq array
c     qq(2:dim ) = q(1:dim-1)
      do l=1,ldim
         do i=1,idim-1
            ii = i+1
            do k=1,kdim-1
               kk = k+ 1
               do j=1,jdim-1
                  jj = j+1
                  qq(jj,kk,ii,l) = q(j,k,i,l)
               end do
            end do
         end do
      end do
c
c     fill in block faces (excluding for now the edges and corners)
c     if ghost-cell data is not available, use duplicate interior
c     data to fill in the block face info.
c
c     note: eddy viscosity and turbulence quantities are always 
c     stored at cell centers, so must override the bcj/k/i arrays

c   /=5 是湍流  
      visflg = 1.
      if (ldim .ne. 5) visflg = 0.
c   
      if (ighost .ne. 0) then
c
c        calculate and store cell-face center data into qq
c
         do l=1,ldim
            do i=1,idim-1
               ii = i+1
               do k=1,kdim-1
                  kk = k+1
                  aa = 1. + bcj(k,i,1)*visflg
                  bb = 1. - bcj(k,i,1)*visflg
                  cc = 1. + bcj(k,i,2)*visflg
                  dd = 1. - bcj(k,i,2)*visflg
                  qq(1,kk,ii,l)      = (aa*qj0(k,i,l,1) +
     .                                  bb*q(1,k,i,l))*0.5
                  qq(jdim+1,kk,ii,l) = (cc*qj0(k,i,l,3) +
     .                                  dd*q(jdim-1,k,i,l))*0.5     
               end do
            end do
            do i=1,idim-1
               ii = i+1
               do j=1,jdim-1
                  jj = j+1
                  aa = 1. + bck(j,i,1)*visflg
                  bb = 1. - bck(j,i,1)*visflg
                  cc = 1. + bck(j,i,2)*visflg
                  dd = 1. - bck(j,i,2)*visflg
                  qq(jj,1,ii,l)      = (aa*qk0(j,i,l,1) +
     .                                  bb*q(j,1,i,l))*0.5
                  qq(jj,kdim+1,ii,l) = (cc*qk0(j,i,l,3) +
     .                                  dd*q(j,kdim-1,i,l))*0.5
               end do
            end do
            do k=1,kdim-1
               kk = k+1
               do j=1,jdim-1
                  jj = j+1
                  aa = 1. + bci(j,k,1)*visflg
                  bb = 1. - bci(j,k,1)*visflg
                  cc = 1. + bci(j,k,2)*visflg
                  dd = 1. - bci(j,k,2)*visflg
                  qq(jj,kk,1,l)      = (aa*qi0(j,k,l,1) + 
     .                                  bb*q(j,k,1,l))*0.5
                  qq(jj,kk,idim+1,l) = (cc*qi0(j,k,l,3) +
     .                                  dd*q(j,k,idim-1,l))*0.5
               end do
            end do
         end do
c
      else
c
c        load duplicate interior data in block faces (1st order approx)
c
         do l=1,ldim
            do i=1,idim-1
               ii = i+1
               do k=1,kdim-1
                  kk = k+1
                  qq(1,kk,ii,l)      = qq(2,kk,ii,l)
                  qq(jdim+1,kk,ii,l) = qq(jdim,kk,ii,l)
               end do
            end do
            do i=1,idim-1
               ii = i+1
               do j=1,jdim-1
                  jj = j+1
                  qq(jj,1,ii,l)      = qq(jj,2,ii,l)
                  qq(jj,kdim+1,ii,l) = qq(jj,kdim,ii,l)
               end do
            end do
            do k=1,kdim-1
               kk = k+1
               do j=1,jdim-1
                  jj = j+1
                  qq(jj,kk,1,l)      = qq(jj,kk,2,l)
                  qq(jj,kk,idim+1,l) = qq(jj,kk,idim,l)
               end do
            end do
         end do
      end if 
c
c     fill in edge values and corner values via extrapolation/averaging
c
      itwo = 2
      if (i2d.ne.0) itwo = 1
c
      do l=1,ldim
c
c        fill in j-k edge values via extapolation in both j and k
c        directions and averaging
c        直接外推再平均，当作均匀网格处理
         do i=2,idim
            do k=1,kdim+1,kdim
               k1 = 1
               k2 = 2
               if (k.eq.kdim+1) then
                  k1 = -1
                  k2 = -2
               end if
               do j=1,jdim+1,jdim
                  j1 = 1
                  j2 = 2
                  if (j.eq.jdim+1) then
                     j1 = -1
                     j2 = -2
                  end if
                  dqj = qq(j+j1,k,i,l) - qq(j+j2,k,i,l)
                  qj  = qq(j+j1,k,i,l) + dqj
                  dqk = qq(j,k+k1,i,l) - qq(j,k+k2,i,l)
                  qk  = qq(j,k+k1,i,l) + dqk
                  qq(j,k,i,l) = (qj + qk)*0.5
               end do
            end do
         end do
c
c        fill in i-j edge values via extapolation in both i and j
c        directions and averaging
c
         do i=1,idim+1,idim
            i1 = 1
            i2 = itwo
            if (i.eq.idim+1) then
               i1 = -1
               i2 = -itwo
            end if
            do k=2,kdim
               do j=1,jdim+1,jdim
                  j1 = 1
                  j2 = 2
                  if (j.eq.jdim+1)  then
                     j1 = -1
                     j2 = -2
                  end if
                  dqj = qq(j+j1,k,i,l) - qq(j+j2,k,i,l)
                  qj  = qq(j+j1,k,i,l) + dqj
                  dqi = qq(j,k,i+i1,l) - qq(j,k,i+i2,l)
                  qi  = qq(j,k,i+i1,l) + dqi
                  qq(j,k,i,l) = (qj + qi)*0.5
               end do
            end do
         end do
c
c        fill in i-k edge values via extapolation in both i and k
c        directions and averaging
c
         do i=1,idim+1,idim
            i1 = 1
            i2 = itwo
            if (i.eq.idim+1) then
               i1 = -1
               i2 = -itwo
            end if
            do k=1,kdim+1,kdim
               k1 = 1
               k2 = 2
               if (k.eq.kdim+1) then
                  k1 = -1
                  k2 = -2
               end if
               do j=2,jdim
                  dqi = qq(j,k,i+i1,l) - qq(j,k,i+i2,l)
                  qi  = qq(j,k,i+i1,l) + dqi
                  dqk = qq(j,k+k1,i,l) - qq(j,k+k2,i,l)
                  qk  = qq(j,k+k1,i,l) + dqk
                  qq(j,k,i,l) = (qi + qk)*0.5
               end do
            end do
         end do
c
c        fill in corner values via extapolation in all 3 directions
c        and averaging
c
         do k=1,kdim+1,kdim
            k1 = 1
            k2 = 2
            if (k.eq.kdim+1) then
               k1 = -1
               k2 = -2
            end if
            do j=1,jdim+1,jdim
               j1 = 1
               j2 = 2
               if (j.eq.jdim+1) then
                  j1 = -1
                  j2 = -2
               end if
               do i=1,idim+1,idim
                  i1 = 1
                  i2 = itwo
                  if (i.eq.idim+1) then
                     i1 = -1
                     i2 = -itwo
                  end if
                  dqi = qq(j,k,i+i1,l) - qq(j,k,i+i2,l)
                  qi  = qq(j,k,i+i1,l) + dqi
                  dqj = qq(j+j1,k,i,l) - qq(j+j2,k,i,l)
                  qj  = qq(j+j1,k,i,l) + dqj
                  dqk = qq(j,k+k1,i,l) - qq(j,k+k2,i,l)
                  qk  = qq(j,k+k1,i,l) + dqk
                  qq(j,k,i,l) = (qi + qj + qk)/3.
               end do
            end do
         end do
      end do
c
      return
      end

 

 

5 cfl3d读入ovrlp.bin的过程

overlay.bin文件的格式

每块数据结构相同

 

read(21) ibpnts,intpts,iipnts, iieptr(ig),iisptr(ig)

! ibpnts，receptor单元个数

! intpts(1:4)，XINTOUT_to_ovlp设置其为 ibpnts, 0, 0, 0

! iipnts，donnoer单元个数

! iieptr，donnor所包含的receptor单元的序号end标号。cfl3d并不读这个数。

! iisptr，donnor所包含的receptor单元的序号begin标号。cfl3d并不读这个数。

 cfl3d并不读最后两个数

 

      lsta = lig(nbl)
      lend = lsta+iipnts-1
      read(21)(jji,kki,iii,dyint,dzint,dxint,         l=lsta,lend)

读取当前块donnor单元的编号和包含的receptor单元中心在其内部的相对坐标 xi eta zeta

intrbc中利用jjig/kkig/llib获取donorstencil 和 在donor stencil中的自然坐标
 s1 = qq(jjig(l)+1, kkig(l)+1, iiig(l)+1, n)

对应的插值结果存在qb(l,n,iset) 

这里l是donor的一维id

      lsta = lbg(nbl)
      lend = lsta+ibpnts-1
      read(21)(jjbg,kkbg,iibg,ibcg,l=lsta,lend)

读取当前块中recptor单元的ijk编号和所在的donor的编号ibc

c     在xupdt中
c         lend = lsta+ibpntsg(nbl,1)-1
c         do 11 l=lsta,lend;    do 10 ll=1,ldim
c         q(jjbg(l),kkbg(l),iibg(l),ll) = qb(ibcg(l),ll,iset)
c     因为在intrbc中，qb是按照donor的排列顺序存储的
c     而这里是按照receptor的顺序设置q，因此需要ibcg(l)就是 receptor(l)的donor的id

最后

read(21)(((blank,j=1,jdim1),k=1,kdim1),i=1,idim1)

读取块中所有单元的iblank

 

 

6 为了做得比maggie好，我们需要确保在边界附近能够实现三线性插值而不是零阶插值或者外推，同时需要利用多块网格对接信息。不过我们不想改动cfl3d，因此一个设想是：

• 无匹配

• 无匹配