DOUG 0.2

CoarseGrid.F90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00035 module CoarseGrid_class
00036   
00037     use Mesh_class
00038     use SpMtx_class
00039     use RealKind
00040     use BinaryHeap
00041 
00042     implicit none
00043 
00044 #include<doug_config.h>
00045 
00047     integer, parameter :: COARSE_CENTER_GEOM  = 1
00048     integer, parameter :: COARSE_CENTER_MEAN  = 2
00049     integer, parameter :: COARSE_CENTER_MERID = 3
00050 
00052     integer, parameter :: COARSE_INTERPOLATION_MULTLIN  = 1
00053     integer, parameter :: COARSE_INTERPOLATION_INVDIST  = 2
00054     integer, parameter :: COARSE_INTERPOLATION_KRIGING  = 3
00055 
00063     type CoarseGridElem
00064         integer :: nfs !< number of fine mesh nodes within 
00065         integer, dimension(:), pointer :: n !< nodes in the corners
00066         integer :: nref !< number of refined elems within
00067         integer :: rbeg !< beginning of the contained refined element list
00068         integer :: lbeg !< beginning of its nodes in elmap
00069     end type
00070 
00071     type RefinedElem
00072         integer :: level !< the level of subdivision it is on
00073         integer :: node !< index of the node this element corresponds to
00074         integer :: next !< as they are held in a linked list
00075         integer :: parent !< the direct parent it belongs to
00076                           !< if it is an initial grid element, it is negative
00077         ! The next three are into CoarseGrid%elmap
00078         integer :: lbeg, lend !< indices to list indicating nodes belonging to el
00079         integer :: lstop !< index to indicate the end of nodes within this el
00080         integer,pointer :: hnds(:) !< indices to hanging nodes (into coords)
00081     end type
00082 
00084     type CoarseGrid
00085         integer :: ncti=-1 !< number of nodes in the initial grid mesh
00086         integer :: nct=-1 !< number of nodes in the mesh (total)
00087         integer :: elnum=-1 !< number of elements in the initial grid
00088         integer :: refnum=-1 !< number of refined elements
00089         integer :: nhn=-1 !< number of hanging nodes
00090         integer :: ngfc=-1 !< number of global coarse freedoms
00091         integer :: nlfc=-1 !< number of local coarse freedoms
00092         integer :: mlvl=-1 !< maximum level of refinement (global value)
00093 
00095         real(kind=xyzk), dimension(:,:), pointer :: coords
00096 
00097         real(kind=xyzk), dimension(:), pointer :: h0
00098 
00099         real(kind=xyzk), dimension(:), pointer :: minvg, maxvg
00100 
00101         integer, dimension(:), pointer :: nc
00102 
00103         type(CoarseGridElem), dimension(:), pointer :: els
00104 
00105         type(RefinedElem), dimension(:), pointer :: refels 
00106 
00108         integer, dimension(:), pointer :: elmap
00109         
00111         integer, dimension(:), pointer :: cfreemap
00112 
00114         integer, dimension(:), pointer :: gl_fmap
00115 
00116         integer, dimension(:), pointer :: lg_fmap
00117 
00118     end type
00119 
00120 contains
00121 
00122   function CoarseGrid_new() result(C)
00123     type(CoarseGrid) :: C
00124 
00125     C%ncti=-1
00126     C%nct=-1
00127     C%elnum=-1
00128     C%refnum=-1
00129     C%nhn=-1
00130     C%ngfc=-1
00131     C%nlfc=-1
00132     C%mlvl=-1
00133 
00134     C%coords => NULL()
00135     C%h0 => NULL()
00136     C%minvg => NULL()
00137     C%maxvg => NULL()
00138     C%nc => NULL()
00139     C%els => NULL()
00140     C%refels => NULL()
00141     C%elmap => NULL()
00142     C%cfreemap => NULL()
00143     C%gl_fmap => NULL()
00144     C%lg_fmap => NULL()
00145     
00146   end function CoarseGrid_new
00147 
00149     subroutine CoarseGrid_allocate(C, nsd, nnode, coords, els, &
00150                                         refels,cfreemap, local)
00151         implicit none
00152         !! The CoarseGrid to destroy
00153         type(CoarseGrid), intent(inout) :: C
00154         !! Dimension and number of nodes in the fine mesh
00155         integer, intent(in), optional :: nsd, nnode ! Not contained in C
00156         !! Flags saying what to allocate
00157         logical, intent(in), optional :: coords, els, refels, cfreemap, local
00158 
00159         if (present(nsd) .and. nsd>0) then
00160             if (.not.associated(C%h0)) allocate(C%h0(nsd))
00161             if (.not.associated(C%minvg)) allocate(C%minvg(nsd))
00162             if (.not.associated(C%maxvg)) allocate(C%maxvg(nsd))
00163             if (.not.associated(C%nc)) allocate(C%nc(nsd))
00164             
00165             if (present(coords) .and. C%nct>0 .and. &
00166                 .not.associated(C%coords) ) then
00167                 allocate(C%coords(nsd,C%nct))
00168             endif
00169         endif
00170 
00171         if (present(els) .and. C%elnum>0 .and. &
00172             .not.associated(C%els) ) then
00173             allocate(C%els(C%elnum))
00174         endif 
00175 
00176         if (present(refels) .and. C%refnum>0 .and. &
00177             .not.associated(C%refels)) then
00178             allocate(C%refels(C%refnum))
00179         endif
00180 
00181         if (present(nnode) .and. .not.associated(C%elmap)) then
00182             allocate(C%elmap(nnode))
00183         endif
00184 
00185         if (present(cfreemap) .and. .not.associated(C%cfreemap)) then
00186             if (C%nlfc>0) then
00187                 allocate(C%cfreemap(C%nlfc))
00188             elseif (C%ngfc>0) then
00189                 allocate(C%cfreemap(C%ngfc))
00190             endif
00191         endif
00192 
00193         if (present(local)) then
00194             if (C%nlfc>0 .and. .not.associated(C%lg_fmap)) then
00195                 allocate(C%lg_fmap(C%nlfc))
00196             endif
00197             if (C%ngfc>0 .and. .not.associated(C%gl_fmap)) then
00198                 allocate(C%gl_fmap(C%ngfc))
00199             endif
00200         endif
00201 
00202     end subroutine CoarseGrid_allocate
00203 
00205     subroutine CoarseGrid_Destroy(C)
00206         implicit none
00208         type(CoarseGrid), intent(inout) :: C
00209 
00210         integer :: i
00211 
00212         if (associated(C%h0)) deallocate(C%h0)
00213         if (associated(C%coords)) deallocate(C%coords)
00214         if (associated(C%minvg)) deallocate(C%minvg)
00215         if (associated(C%maxvg)) deallocate(C%maxvg)
00216         if (associated(C%nc)) deallocate(C%nc)
00217         if (associated(C%els)) then
00218             do i=1,size(C%els,dim=1)
00219                 if (associated(C%els(i)%n)) deallocate(C%els(i)%n)
00220             enddo
00221             deallocate(C%els)
00222         endif
00223         if (associated(C%refels)) then
00224             do i=1,C%refnum
00225                 if (associated(C%refels(i)%hnds)) deallocate(C%refels(i)%hnds)
00226             enddo
00227             deallocate(C%refels)
00228         endif
00229         if (associated(C%elmap)) deallocate(C%elmap)
00230         if (associated(C%cfreemap)) deallocate(C%cfreemap)
00231         if (associated(C%lg_fmap)) deallocate(C%lg_fmap)
00232         if (associated(C%gl_fmap)) deallocate(C%gl_fmap)
00233 
00234     end subroutine CoarseGrid_Destroy
00235 
00236     subroutine CoarseGrid_pl2D_plotMesh(C, INIT_CONT_END)
00237         use globals, only: stream
00238         implicit none
00239 
00240         type(CoarseGrid), intent(in)      :: C
00241         integer,    intent(in), optional  :: INIT_CONT_END
00242 
00243         real(kind=xyzk), dimension(:), pointer :: xc, yc
00244         real(kind=xyzk) :: xmin, xmax, ymin, ymax
00245 
00246         integer :: e, i, k, lvl
00247         character*6 :: buf61, buf62, buf63
00248         
00249         integer,parameter :: nsd=2, tnsd=4
00250 
00251         type(CoarseGridElem), pointer :: el
00252         type(RefinedElem), pointer :: rel
00253 
00254         real(kind=xyzk) :: mins(nsd,0:C%mlvl), maxs(nsd,0:C%mlvl)
00255         real(kind=xyzk), pointer :: ct(:), parct(:), pt(:)
00256         
00257 #ifdef D_WANT_PLPLOT_YES
00258 
00259         write(stream, *)
00260         write(stream, *) '[CoarseGrid_pl2D_plotMesh] : Plotting 2D coarse mesh.'
00261     
00262         !if (nsd /= 2) then
00263         !    write(stream, *) '   Only for 2D meshes. Skipping...'
00264         !    return
00265         !end if
00266 
00267         !tnsd=2**nsd
00268         
00269         xc => C%coords(1,:)
00270         yc => C%coords(2,:)
00271 
00272         if (.not.present(INIT_CONT_END).or.&
00273             (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00274 
00275              xmin = minval(xc)
00276              xmax = maxval(xc)
00277              ymin = minval(yc)
00278              ymax = maxval(yc)
00279 
00280              xmin = xmin - (xmax-xmin)/10.0
00281              xmax = xmax + (xmax-xmin)/10.0
00282              ymin = ymin - (ymax-ymin)/10.0
00283              ymax = ymax + (ymax-ymin)/10.0
00284 
00285              call plsdev("xwin")
00286              call plinit()
00287 
00288              call plenv (xmin, xmax, ymin, ymax, 0, 0);
00289 
00290              call plcol0(1) ! red
00291 
00292              write(buf61, '(i6)') C%elnum
00293              write(buf62, '(i6)') C%refnum
00294              write(buf63, '(i6)') C%mlvl
00295              call pllab( '(x)', '(y)', &
00296                'Coarse : elnum='//buf61//'; refnum='//buf62//'; mlvl='//buf63)
00297         end if
00298 
00299         call plcol0(7) ! grey
00300         call plssym(0.0d0, 2.0d0)
00301         !call plpoin(C%nct, xc, yc, 1)
00302 
00303         ! Draw Elements
00304         do e=1,C%elnum
00305             el=>C%els(e)
00306             
00307             ! Dont try to draw discarded elements
00308             if (el%nfs==0) cycle
00309             
00310             ! Get the initial coarse element bounds into mins/maxs
00311             mins(:,0)=C%coords(:,el%n(1)); maxs(:,0)=mins(:,0)
00312             do i=2,tnsd
00313                 do k=1,nsd
00314                     if (C%coords(k,el%n(i))<mins(k,0)) then
00315                         mins(k,0)=C%coords(k,el%n(i))
00316                     else if (C%coords(k,el%n(i))>maxs(k,0)) then
00317                         maxs(k,0)=C%coords(k,el%n(i))
00318                     endif
00319                 enddo
00320             enddo
00321             call plcol0(11) ! grey
00322 
00323             ! Draw the initial grid box
00324             call pljoin(mins(1,0),maxs(2,0),maxs(1,0),maxs(2,0)) ! top
00325             call pljoin(mins(1,0),mins(2,0),maxs(1,0),mins(2,0)) ! bottom
00326             call pljoin(mins(1,0),maxs(2,0),mins(1,0),mins(2,0)) ! left
00327             call pljoin(maxs(1,0),maxs(2,0),maxs(1,0),mins(2,0)) ! right
00328 
00329             ! Walk the refined elements witin
00330             i=el%rbeg
00331             if (i/=-1) then
00332                 call plcol0(11) ! grey
00333 
00334                 ! Deal with the first refinement
00335                 ct=>C%coords(:,C%refels(i)%node) ! the center
00336                 call pljoin(mins(1,0),ct(2),maxs(1,0),ct(2)) ! left-right
00337                 call pljoin(ct(1),mins(2,0),ct(1),maxs(2,0)) ! top-bottom
00338 
00339                 ! Move on to the other refinements
00340                 i=C%refels(i)%next
00341                 do while (i/=-1)
00342                     rel=>C%refels(i)
00343                 
00344                     lvl=rel%level-1
00345                     ct=>C%coords(:,rel%node) ! center
00346                     parct=>C%coords(:,C%refels(rel%parent)%node) ! parent center
00347                
00348                     ! Adjust mins and maxs as needed
00349                     do k=1,nsd
00350                         if (parct(k)>ct(k)) then
00351                             maxs(k,lvl)=parct(k)
00352                             mins(k,lvl)=mins(k,lvl-1)
00353                         else
00354                             maxs(k,lvl)=maxs(k,lvl-1)
00355                             mins(k,lvl)=parct(k)
00356                         endif
00357                     enddo
00358 
00359                     ! Draw the division
00360                     call pljoin(mins(1,lvl),ct(2),maxs(1,lvl),ct(2)) ! l-r
00361                     call pljoin(ct(1),mins(2,lvl),ct(1),maxs(2,lvl)) ! t-b
00362 
00363                     !if (associated(rel%hnds)) &
00364                     !    write(stream,*) ACHAR(rel%node+32),ACHAR(rel%hnds+32)
00365 
00366                
00367                     ! Move on
00368                     i=rel%next
00369                 enddo
00370             endif
00371         enddo
00372        
00373         call plcol0(12)
00374         call plssym(0.0d0,5.0d0) 
00375         call plpoin(C%ncti+C%refnum,C%coords(1,:),C%coords(2,:),1)
00376         call plpoin(C%nhn,C%coords(1,C%nct-C%nhn+1:C%nct), &
00377                           C%coords(2,C%nct-C%nhn+1:C%nct),1)
00378 
00379 !        ! Debugging
00380 !        call plcol0(11)
00381 !        call plssym(0.0d0,1.0d0) 
00382 !        if (C%ncti>0) then
00383 !        do i=1,C%nct
00384 !                call plpoin(1,C%coords(1,i),C%coords(2,i),32+i)
00385 !        enddo; endif
00386 
00387 
00388         if (.not.present(INIT_CONT_END).or.&
00389            (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_END))) then
00390             call plend()
00391         end if
00392 
00393 #else
00394         write(stream, '(/a)') ' [CoarseGrid_pl2D_plotMesh] : Compiled w/o'//&
00395                 ' plotting support!'
00396 #endif
00397 
00398     end subroutine CoarseGrid_pl2D_plotMesh
00399 
00400 
00403     function getelem(coords,mins,h,nc) result (ind)
00404         use RealKind
00405         implicit none
00406         !! Coordinates of the point being mapped
00407         real(kind=xyzk), dimension(:), intent(in) :: coords, mins, h
00408         !! The maximum values of dimensions
00409         integer, dimension(:), intent(in) :: nc
00410 
00411         real(kind=xyzk) :: rx(size(coords)) 
00412         integer :: i,ind
00413 
00414         ! Scale to grid
00415         rx=(coords-mins)/h
00416         
00417         ! And determine the element it belongs to
00418         ind=aint(rx(1))
00419         do i=2,size(rx,dim=1)
00420             ind=ind*(nc(i)-1) + (aint(rx(i)))
00421         enddo
00422 
00423         ind=ind+1        
00424     end function getelem
00425             
00427     subroutine getRefBounds(refi,C,nsd,minv,maxv)
00428         integer, intent(in) :: refi
00429         type(CoarseGrid), intent(in) :: C
00430         integer, intent(in) :: nsd
00431         real(kind=xyzk), intent(out) :: minv(:), maxv(:)        
00432 
00433         real(kind=xyzk),pointer :: ct(:)
00434         type(RefinedElem),pointer :: ref
00435         type(CoarseGridElem),pointer :: el
00436         integer :: i,k
00437         logical :: mi(nsd), ma(nsd)
00438 
00439         ref=>C%refels(refi); mi=.false.; ma=.false.
00440 
00441         ! Get the initial element "center point" into bounds
00442         minv(:)=C%coords(:,ref%node); maxv(:)=minv(:)
00443 
00444         ! Get the bounds provided by the refined tree
00445         do while (ref%parent>0)
00446            ref=>C%refels(ref%parent)
00447            ct=>C%coords(:,ref%node)
00448 
00449            ! Adjust mins and maxs as needed
00450            do k=1,nsd
00451               if (.not.ma(k) .and. ct(k)>maxv(k)) then
00452                 maxv(k)=ct(k); ma(k)=.true.
00453               else if (.not.mi(k) .and. ct(k)<=minv(k)) then
00454                 minv(k)=ct(k); mi(k)=.true.
00455               endif
00456            enddo
00457 
00458         enddo
00459 
00460         el=>C%els(-ref%parent)
00461         ct=>C%coords(:,C%refels(refi)%node)
00462 
00463         ! Fill in the gaps the refined tree left
00464         do i=1,2**nsd
00465             do k=1,nsd
00466                 if (.not.mi(k) .and. C%coords(k,el%n(i))<minv(k)) then
00467                     minv(k)=C%coords(k,el%n(i)); mi(k)=.true.
00468                 else if (.not.ma(k) .and. C%coords(k,el%n(i))>maxv(k)) then
00469                     maxv(k)=C%coords(k,el%n(i)); ma(k)=.true.
00470                 endif
00471             enddo
00472         enddo
00473     end subroutine getRefBounds
00474 
00477     subroutine adjustBounds(ct,pt,nsd,minv,maxv)
00478         integer, intent(in) :: nsd
00479         real(kind=xyzk),intent(in) :: ct(nsd), pt(nsd)
00480         real(kind=xyzk),intent(inout) :: minv(nsd),maxv(nsd)
00481 
00482         integer :: i
00483 
00484         do i=1,nsd
00485             if (ct(i)<=pt(i)) then
00486                 minv(i)=ct(i)
00487             else
00488                 maxv(i)=ct(i)
00489             endif   
00490         enddo
00491     end subroutine adjustBounds
00492 
00498     function getDir(ds,nsd) result (n)
00499         real(kind=xyzk), intent(in) :: ds(:)
00500         integer, intent(in) :: nsd
00501         integer :: n
00502         
00503          n=1
00504          if (ds(1)<0) n=n+1
00505          if (ds(2)<0) n=n+2
00506          if (nsd==3) then
00507               if (ds(3)<0) n=n+4
00508          endif
00509     end function getDir
00510 
00514     function getNextElem(eli,dir,nsd,C) result (nel)
00515         integer,intent(in) :: eli
00516         integer,intent(in) :: dir
00517         integer,intent(in) :: nsd
00518         type(CoarseGrid), intent(in) :: C
00519 
00520         integer :: nel 
00521 
00522         if (nsd==2) then
00523             if (dir==-1) then 
00524                 if ( eli<=C%nc(2)-1 ) then; nel=0
00525                 else; nel=eli-(C%nc(2)-1); endif
00526             else if (dir==1) then
00527                 if ( eli>C%elnum-(C%nc(2)-1) ) then; nel=0
00528                 else; nel=eli+(C%nc(2)-1); endif
00529             else if (dir==-2) then
00530                 if ( mod( eli-1 , C%nc(2)-1 ) == 0 ) then; nel=0 ! before first
00531                 else; nel=eli-1; endif
00532             else if (dir==2) then
00533                 if ( mod( eli , C%nc(2)-1 ) == 0 ) then; nel=0  ! after last
00534                 else; nel=eli+1; endif
00535             endif
00536         else if (nsd==3) then
00537             if (dir==-1) then 
00538                 if ( eli<=(C%nc(2)-1)*(C%nc(3)-1) ) then; nel=0
00539                 else; nel=eli-(C%nc(2)-1)*(C%nc(3)-1); endif
00540             else if (dir==1) then
00541                 if ( eli>C%elnum-(C%nc(2)-1)*(C%nc(3)-1) ) then; nel=0
00542                 else; nel=eli+(C%nc(2)-1)*(C%nc(3)-1); endif
00543             else if (dir==-2) then
00544                 if ( mod((eli-1)/(C%nc(3)-1),C%nc(2)-1 ) == 0 ) then
00545                      nel=0 ! before first
00546                 else; nel=eli-(C%nc(3)-1); endif
00547             else if (dir==2) then
00548                 if ( mod((eli-1)/(C%nc(3)-1)+1,C%nc(2)-1 ) == 0 ) then
00549                      nel=0  ! after last
00550                 else; nel=eli+(C%nc(3)-1); endif 
00551             else if (dir==-4) then
00552                 if ( mod( eli-1 , C%nc(3)-1 ) == 0 ) then; nel=0 ! before first
00553                 else; nel=eli-1; endif
00554             else if (dir==4) then
00555                 if ( mod( eli , C%nc(3)-1 ) == 0 ) then; nel=0  ! after last
00556                 else; nel=eli+1; endif
00557             endif
00558  
00559         endif
00560 
00561     end function getNextElem
00562         
00565     function  getNeighbourEl(el,dir,nsd,nsame,C) result (nb)
00566         integer, intent(in) :: el !< the element whose neighbour we want
00567         integer, intent(in) :: dir !< the direction to get the neighbour from
00568         integer, intent(in) :: nsd !< num of dimensions
00569         integer, intent(in) :: nsame(:) !< next refinements of same level
00570         type(CoarseGrid),intent(in) :: C !< the coarse grid itself
00571         integer :: nb !< The neighbours index in refels 
00572 
00573         integer :: stack(C%mlvl)
00574         integer :: i,j,k
00575         integer :: flag
00576         type(RefinedElem),pointer :: p
00577         real(kind=xyzk) :: pt(nsd)
00578 
00579         ! Move rootward in the tree trying to find a place to take the step
00580         j=el; nb=0; stack(C%refels(j)%level)=0
00581         do 
00582            if (C%refels(j)%parent<0) then
00583                 nb=-getNextElem(-C%refels(j)%parent,dir,nsd,C)
00584                 if (nb/=0) then
00585                     if (C%els(-nb)%nfs==0) then; nb=0; ! "overlook" empty els
00586                     else if (C%els(-nb)%rbeg>0) then; nb=C%els(-nb)%rbeg
00587                     endif
00588                 endif
00589                 exit
00590            else
00591                 p=>C%refels(C%refels(j)%parent)
00592                 pt=C%coords(:,C%refels(j)%node)-C%coords(:,p%node)
00593                 stack(p%level)=getDir(pt,nsd)
00594                 j=C%refels(j)%parent
00595 
00596                 ! Find the direction to which to try to move back to
00597                 flag=stack(p%level)-1
00598                 if (dir>0) then
00599                     if (iand(flag,dir)==dir) flag=flag-dir
00600                 else
00601                     if (iand(flag,-dir)==0) flag=flag-dir
00602                 endif
00603                 flag=flag+1
00604                 ! If no such move is possible, go on upwards
00605                 if (flag==stack(p%level)) then
00606                     stack(p%level)=stack(p%level)+dir; cycle
00607                 endif
00608 
00609                 ! otherwise, try all the children of this parent
00610                 ! and hope one matches
00611                 k=p%next
00612                 do while (k>0)
00613                     pt=C%coords(:,C%refels(k)%node)-C%coords(:,p%node)
00614                     if (flag==getDir(pt,nsd)) then ! it did
00615                         nb=k; exit ! so use it
00616                     endif
00617                     k=nsame(k)
00618                 enddo
00619 
00620                 ! if it didnt match, current parent is its neighbour
00621                 if (nb==0) nb=j;
00622 
00623                 exit;
00624            endif
00625         enddo
00626 
00627         ! Move back leafward as much as is possible and sufficient
00628         if (nb>0 .and. nb/=j) then
00629         do
00630            p=>C%refels(nb); flag=stack(p%level)
00631            if (flag==0) exit ! stop if at the same level as we started
00632 
00633            ! See if there is a child of this parent in the desired direction
00634            k=p%next
00635            if (k<=0) exit
00636            if (C%refels(k)%level==p%level+1) then
00637            do while (k>0)
00638               pt=C%coords(:,C%refels(k)%node)-C%coords(:,p%node)
00639               if (flag==getDir(pt,nsd)) then ! there was
00640                   nb=k; flag=0; exit ! so use it
00641               endif
00642               k=nsame(k)
00643            enddo
00644            endif
00645 
00646            if (flag/=0) exit
00647         enddo      
00648         endif
00649 
00650     end function getNeighbourEl
00651 
00652 end module CoarseGrid_class