DOUG 0.2

CreateCoarse.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 
00131 
00132 module CreateCoarseGrid
00133 
00134 use RealKind
00135 
00138     real(kind=xyzk), parameter :: meanpow=1.0_xyzk
00139 
00140     private :: CreateHangingNodes, CreateCoarseMesh, CreateCoarseFreemap
00141 contains
00142 
00144 subroutine CreateCoarse(M,C)
00145         use RealKind
00146         use CoarseGrid_class
00147         use globals 
00148         
00149         implicit none
00150 
00152         type(CoarseGrid), intent(inout) :: C
00153 
00154         type(Mesh), intent(in) :: M
00155 
00156         ! Make a choice between different possibilities
00157         selectcase(mctls%center_type)
00158             case (COARSE_CENTER_GEOM) ! 1
00159                 call CreateCoarseMesh(M,C,chooseGeometricCenter)
00160             case (COARSE_CENTER_MEAN) ! 2
00161                 call CreateCoarseMesh(M,C,chooseMeanCenter)
00162             case (COARSE_CENTER_MERID) ! 3
00163                 call CreateCoarseMesh(M,C,chooseMeridianCenter) 
00164             case default
00165                 call DOUG_abort("No such center choosing method")
00166         endselect
00167 
00168         ! Create the freemap
00169         call CreateCoarseFreemap(C,M)
00170 end subroutine CreateCoarse
00171 
00173 subroutine CreateCoarseMesh(M, C, choosecenter)
00174         use RealKind
00175         use CoarseGrid_class
00176         use globals
00177         
00178         implicit none
00179 
00181         type(CoarseGrid), intent(inout) :: C
00182 
00184         type(Mesh), intent(in) :: M
00185 
00186         interface
00187 
00188             subroutine ChooseCenter(pt,cpt,pts,elmap,el,refels,flags,minv,maxv)
00189                 use RealKind
00190                 use CoarseGrid_class
00191                 real(kind=xyzk), intent(out) :: pt(:) !< the new center to output
00192                 real(kind=xyzk), intent(in) :: cpt(:) !< coordinates of the prev cn
00193                 real(kind=xyzk), intent(in) :: pts(:,:) !< points array
00194                 integer, intent(in) :: elmap(:) !< indices to pts
00195                 integer, intent(in) :: el !< index of the element being divided
00196                 type(RefinedElem), intent(in) :: refels(:) !< refined elements 
00197                 integer, intent(in) :: flags !< flags saying which dir to div
00198                 real(kind=xyzk), intent(in) :: minv(:), maxv(:)  !< elem bounds
00199            end subroutine ChooseCenter
00200         end interface
00201         
00202         real(kind=xyzk) :: minv(M%nsd), maxv(M%nsd), ln(M%nsd)
00203         
00204         real(kind=xyzk) :: buf, multbuf
00205         integer :: i, j, k, nd, nb, refpt, hangpt, par
00206         integer ::  el, flags, pcnt
00207         
00208         integer :: counts(2**M%nsd), cnt
00209         real(kind=xyzk) :: h1(M%nsd), pt(M%nsd), pt2(M%nsd)
00210         real(kind=xyzk), pointer :: ncoords(:,:)
00211         type(BHeap) :: queue
00212 
00213         type(RefinedElem), pointer :: cur, new
00214 
00215         logical :: inuse(M%nnode) ! A filter for weeding unused fine nodes
00216         integer, allocatable :: nsame(:) ! A pointer array for refined els
00217         integer, allocatable :: ctiremap(:) ! For remaping coarse nodes
00218 
00219         ! Mark which nodes have freedoms attached
00220         inuse=.false.
00221         do i=1,M%nell
00222             if (M%eptnmap(i)<=M%nparts .and. M%eptnmap(i)>0) then 
00223             do j=1,M%nfrelt(i)
00224                 inuse(M%freemap(M%mhead(j, i)))=.true.
00225             enddo
00226             endif
00227         enddo
00228 
00229         ! Allocate memory for the basic structures 
00230         call CoarseGrid_allocate(C,M%nsd)
00231 
00232         !***********************************************************
00233         ! Find the minimum and maximum coordinates
00234         !***********************************************************
00235         do i=1,M%nsd
00236             ! Find the extremes
00237             C%minvg(i)=minval(M%coords(i,:),MASK=inuse)
00238             C%maxvg(i)=maxval(M%coords(i,:),MASK=inuse)
00239 
00240             ! And stretch them a little outward
00241             ln(i)=(C%maxvg(i)-C%minvg(i))
00242             C%minvg(i)=C%minvg(i)-ln(i)*0.005_xyzk
00243             C%maxvg(i)=C%maxvg(i)+ln(i)*0.005_xyzk
00244             ln(i)=ln(i)*1.01_xyzk
00245         enddo
00246  
00247         !***********************************************************
00248         ! Find the appropriate size for the initial grid elements
00249         !***********************************************************
00250 
00251         ! Find how many to use in each direction
00252         ! Should give quite good approximations to squares for el. shapes
00253         ! (Currently, C%nc(i) means number of elements in that direction)
00254         ! Uses formulae derived by solving prod(C%nc(:))=mctls%maxcie
00255         ! and C%nc(i)/C%nc(j) = ln(i)/ln(j), 1<=i,j<=M%nsd
00256         ! and then truncating C%nc(i) thus obtaining whole number ratios.
00257         ! Averaging might be better, but maxcie could not be enforced then
00258         if (M%nsd==2) then
00259             C%nc(1)=anint((mctls%maxcie*(ln(1)/ln(2)))**(0.5_xyzk)) 
00260             C%nc(2)=anint((mctls%maxcie*(ln(2)/ln(1)))**(0.5_xyzk))
00261         else if (M%nsd==3) then
00262             C%nc(1)=anint((mctls%maxcie*ln(1)*ln(1)/(ln(2)*ln(3)))**(0.33_xyzk))
00263             C%nc(2)=anint((mctls%maxcie*ln(2)*ln(2)/(ln(1)*ln(3)))**(0.33_xyzk))
00264             C%nc(3)=anint((mctls%maxcie*ln(3)*ln(3)/(ln(1)*ln(2)))**(0.33_xyzk))
00265         endif
00266         C%nc=max(C%nc,1)
00267         C%h0=ln/C%nc
00268         C%elnum=product(C%nc)
00269 
00270         if (sctls%verbose>3) &
00271                 write (stream,*) "Initial coarse grid has dimensions ", &
00272                                 C%nc," with a total of ",C%elnum," elements"
00273 
00274         ! Change nc from number of elements to number of nodes
00275         C%nc=C%nc+1
00276         C%ncti=product(C%nc)
00277 
00278         !***********************************************************
00279         ! Create the actual initial grid
00280         !***********************************************************
00281 
00282         ! Calcluate the number of refined elements and set nct
00283         C%refnum=max(mctls%maxnd-C%ncti,0)
00284         C%nct=max(mctls%maxnd,C%ncti)
00285  
00286         ! Allocate memory for the grid nodes and elements
00287         call CoarseGrid_allocate(C,M%nsd,nnode=M%nnode, &
00288                         coords=.true.,els=.true.,refels=.true.)
00289 
00290         ! Create initial coarse grid nodes
00291         nd=1
00292         if (M%nsd==2) then ! 2D
00293             do i=0,C%nc(1)-1; do j=0,C%nc(2)-1
00294                 C%coords(1,nd)=i*C%h0(1)+C%minvg(1)
00295                 C%coords(2,nd)=j*C%h0(2)+C%minvg(2)
00296                 nd=nd+1
00297             enddo; enddo
00298         else ! 3D
00299             do i=0,C%nc(1)-1; do j=0,C%nc(2)-1; do k=0,C%nc(3)-1
00300                 C%coords(1,nd)=i*C%h0(1)+C%minvg(1)
00301                 C%coords(2,nd)=j*C%h0(2)+C%minvg(2)
00302                 C%coords(3,nd)=k*C%h0(3)+C%minvg(3)
00303                 nd=nd+1
00304             enddo; enddo; enddo
00305         endif
00306 
00307         ! Create initial coarse grid elements
00308         el=1
00309         if (M%nsd==2) then ! 2D
00310             do i=1,C%nc(1)-1; do j=1,C%nc(2)-1
00311                 C%els(el)%nfs=0; C%els(el)%nref=0; 
00312                 C%els(el)%rbeg=-1
00313                 
00314                 allocate(C%els(el)%n(4)) ! 2**2=4
00315                 ! indices of the nodes that bound this element
00316                 nb=(i-1)*C%nc(2)+j
00317                 ! It is sensible to keep them sorted
00318                 C%els(el)%n(1)=nb
00319                 C%els(el)%n(2)=nb+1
00320                 C%els(el)%n(3)=nb+C%nc(2)
00321                 C%els(el)%n(4)=nb+C%nc(2)+1
00322                 el=el+1
00323             enddo; enddo
00324         else ! 3D
00325             do i=1,C%nc(1)-1; do j=1,C%nc(2)-1; do k=1,C%nc(3)-1
00326                 C%els(el)%nfs=0; C%els(el)%nref=0; 
00327                 C%els(el)%rbeg=-1
00328                 allocate(C%els(el)%n(8)) ! 2**3=8
00329 
00330                 ! indices of the nodes that bound this element
00331                 nb=((i-1)*C%nc(2)+(j-1))*C%nc(3)+k
00332 
00333                 ! It is sensible to keep them in this sorted order
00334                 C%els(el)%n(1)=nb
00335                 C%els(el)%n(2)=nb+1
00336                 C%els(el)%n(3)=nb+C%nc(3)
00337                 C%els(el)%n(4)=nb+C%nc(3)+1
00338                 C%els(el)%n(5)=nb+C%nc(3)*C%nc(2)
00339                 C%els(el)%n(6)=nb+C%nc(3)*C%nc(2)+1
00340                 C%els(el)%n(7)=nb+C%nc(3)*(C%nc(2)+1)
00341                 C%els(el)%n(8)=nb+C%nc(3)*(C%nc(2)+1)+1
00342                 el=el+1
00343             enddo; enddo; enddo
00344         endif
00345         
00346         !***********************************************************
00347         ! Create the elmap for quick retrieval of fine nodes
00348         !***********************************************************
00349 
00350         ! Allocate nodes into the elements
00351         do nd=1,M%nnode
00352             if (inuse(nd)) then
00353                 ! Get the coarse element containing it
00354                 el=getelem(M%coords(:,nd),C%minvg,C%h0,C%nc)
00355                 C%els(el)%nfs=C%els(el)%nfs+1
00356             endif
00357         enddo
00358 
00359         ! Counting sort the nodes by elements they belong to
00360 
00361         ! Since we already have the counts, just find the positions
00362         C%els(1)%lbeg=1
00363         do i=2, C%elnum
00364             C%els(i)%lbeg=C%els(i-1)%lbeg+C%els(i-1)%nfs
00365         enddo
00366 
00367         ! Set the elements into elmap
00368         do nd=1,M%nnode
00369             if (inuse(nd)) then
00370                 ! Get the coarse element containing it
00371                 el=getelem(M%coords(:,nd),C%minvg,C%h0,C%nc)
00372                 C%elmap(C%els(el)%lbeg)=nd
00373                 ! And move the writer head for freedoms in this element by one
00374                 C%els(el)%lbeg=C%els(el)%lbeg+1
00375             endif
00376         enddo
00377 
00378         ! Init the priority queue (MaxHeap)
00379         queue = BHeap_new()
00380         call BHeap_init(queue,C%elnum+C%refnum)
00381 
00382         ! Recalculate the counts and fill the priority queue
00383         C%els(1)%lbeg=1; buf=sqrt(real(M%nsd,kind=xyzk))
00384         if (C%els(1)%nfs>0) call BHeap_insert(queue,nint(buf*C%els(1)%nfs),1)
00385         do i=2, C%elnum
00386            C%els(i)%lbeg=C%els(i-1)%lbeg+C%els(i-1)%nfs
00387            if (C%els(i)%nfs>0) call BHeap_insert(queue,nint(buf*C%els(i)%nfs),i)
00388         enddo
00389         
00390         !***********************************************************
00391         ! Refine the initial elements
00392         !***********************************************************
00393 
00394         ! Allocate the array nsame for a linking same level refinements
00395         allocate(nsame(C%refnum))
00396        
00397         ! Set the initial mlvl to 1
00398         C%mlvl=1
00399        
00400         ! Refine the elements while we can
00401         refpt=1; hangpt=C%nct
00402         do while (C%ncti+refpt<=hangpt)!(refpt<=C%refnum)
00403             ! Get the least balanced element number and its imbalance
00404             el=BHeap_maxi(queue)
00405             pcnt=BHeap_maxv(queue)
00406             
00407             ! Remove the element from the queue
00408             call BHeap_delmax(queue)
00409 
00410             ! If the largest imbalance is smaller than cutoff, stop refinement
00411             if (pcnt<=2*mctls%cutbal) exit
00412 
00413             ! Otherwise refine
00414             if (el<=C%elnum) then ! We need to refine an initial grid element
00415                 ! Set new to point to the element being created
00416                 new=>C%refels(refpt)
00417                 
00418                 ! Create the refined element
00419                 new%level=1
00420                 new%node=C%ncti+refpt
00421                 new%parent=-el
00422                 new%next=-1
00423 
00424                 ! Set beg and end to cover the initial grid element
00425                 new%lbeg=C%els(el)%lbeg
00426                 new%lend=C%els(el)%lbeg+C%els(el)%nfs-1
00427                 new%lstop=new%lend
00428                 nsame(refpt)=0
00429 
00430                 ! Create the geometric center for passing to chooser
00431                 pt2=C%coords(:,C%els(el)%n(1))+C%h0/2.0_xyzk
00432            
00433                 ! Choose the appropriate center for the element
00434                 call ChooseCenter(pt, pt2, &
00435                         M%coords,C%elmap(new%lbeg:new%lend), &
00436                         new%parent,C%refels,-1, &
00437                         C%coords(:,C%els(el)%n(1)), & ! mins
00438                         C%coords(:,C%els(el)%n(2**M%nsd))) ! maxs
00439      
00440                 ! Set it as center
00441                 C%coords(:, new%node ) = pt
00442                         
00443                 ! Modify the initial grid element as needed
00444                 !C%els(el)%ncs=1
00445                 C%els(el)%rbeg=refpt
00446 
00447                 ! And put this new element into the queue
00448                 call BHeap_insert(queue, &
00449                         C%els(el)%nfs,C%elnum+refpt)
00450 
00451                 ! Try to create the associated hanging nodes
00452                 nullify(new%hnds)
00453                 if (mctls%hanging_nodes) &
00454                 call CreateHangingNodes(refpt,hangpt,M%nsd,nsame, &
00455                         C%coords(:,C%els(el)%n(1)), & ! mins
00456                         C%coords(:,C%els(el)%n(2**M%nsd)),C) 
00457 
00458                 refpt=refpt+1
00459             else ! We need to refine an already refined element
00460                
00461                 ! Get the element to subdivide and the initial el it belongs to
00462                 el=el-C%elnum
00463                 par=C%refels(el)%parent
00464  
00465                 ! Set new/cur to point to the el being created/divided
00466                 new=>C%refels(refpt); cur=>C%refels(el)
00467                 
00468                 ! Create the framework for the new refined element
00469                 new%level=cur%level+1
00470                 new%node=C%ncti+refpt
00471                 new%parent=el
00472             
00473                 ! Get the center and reach of the divisible element
00474                 pt=C%coords(:,cur%node)
00475                 call getRefBounds(el,C,M%nsd,minv,maxv)
00476              
00477                 ! Count the surrounding nodes to see which to subdivide
00478                 counts=0
00479 
00480                 do i=cur%lbeg,cur%lend
00481                     pt2=M%coords(:,C%elmap(i))-pt
00482 
00483                     ! Determine which region the node belongs to
00484                     flags=getDir(pt2,M%nsd)                    
00485                     counts(flags)=counts(flags)+1
00486                 enddo
00487 
00488                 ! Find which area has the most nodes and how many that is
00489                 flags=maxloc(counts,dim=1)
00490                 cnt=counts(flags)
00491 
00492                 ! If the new node would be finer than allowed, end refinement
00493                 if (cnt<=mctls%cutbal) exit
00494 
00495                 ! Rearrange the parents fine node list
00496                 !  and take a part of it for the new node
00497                 ! Use the method used inside QuickSort implementations
00498                 i=cur%lbeg; k=cur%lend
00499                 do
00500                     ! Find one in the first part that is in our region
00501                     do i=i,k  
00502                     if (getDir(M%coords(:,C%elmap(i))-pt,M%nsd)==flags ) exit
00503                     enddo
00504 
00505                     ! Find one in the second part that is not
00506                     do k=k,i,-1  
00507                     if (getDir(M%coords(:,C%elmap(k))-pt,M%nsd)/=flags ) exit
00508                     enddo
00509 
00510                     ! If the ends meet, stop
00511                     if (i>=k) exit ! Note i=k+1 is also possible
00512 
00513                     ! Otherwise switch them
00514                     j=C%elmap(k); C%elmap(k)=C%elmap(i); C%elmap(i)=j
00515                 enddo
00516 
00517                 ! Use the end part for this nodes list
00518                 new%lbeg=i
00519                 new%lend=cur%lend
00520                 new%lstop=new%lend
00521  
00522                 ! Add the element to the refined linked list after its parent
00523                 new%next=cur%next
00524                 cur%next=refpt
00525 
00526                 ! Add the element to the nsame linked list
00527                 nsame(refpt)=0
00528                 if (new%next>0) then
00529                 if (C%refels(new%next)%level==new%level) then
00530                         nsame(refpt)=new%next
00531                 endif;  endif
00532                 
00533                 ! Cut the tail away from its previous owner
00534                 cur%lend=i-1
00535 
00536                 ! Choose the appropriate center for the element
00537                 pt2=M%coords(:,C%elmap(new%lbeg))
00538                 call adjustBounds(pt,pt2,M%nsd,minv,maxv)
00539                 call ChooseCenter(pt, C%coords(:,cur%node), &
00540                         M%coords,C%elmap(new%lbeg:new%lend), &
00541                         el,C%refels,flags,minv,maxv)
00542 
00543                 ! Set it as the center
00544                 C%coords(:,new%node)=pt
00545 
00546                 ! Put this new element into the queue
00547                 call BHeap_insert(queue, &
00548                         cnt,C%elnum+refpt)
00549                 
00550                 ! And put the old element we subdivided back in as well
00551                 call BHeap_insert(queue, &
00552                         pcnt-cnt,C%elnum+el)
00553                 
00554                 ! Update the max lvl if neccessary
00555                 if (C%mlvl<new%level) C%mlvl=new%level
00556 
00557                 ! And try to create the associated hanging nodes
00558                 nullify(new%hnds)
00559                 if (mctls%hanging_nodes) &
00560                 call CreateHangingNodes(refpt,hangpt,M%nsd,nsame,minv,maxv,C)
00561        
00562                 ! Increase the pointer
00563                 refpt=refpt+1
00564             endif
00565         enddo
00566         
00567         ! Change the counts to reflect reality
00568         C%refnum=refpt-1; C%nhn=C%nct-hangpt; C%nct=C%ncti+C%refnum+C%nhn
00569 
00570         ! Free the memory
00571         call BHeap_destroy(queue); deallocate(nsame)
00572 
00573         !************************************************************
00574         ! Create a new coords array
00575         !************************************************************
00576 
00577         ! Locate the unused initial grid nodes by marking all the used ones
00578         allocate(ctiremap(0:C%ncti)); ctiremap=0
00579         do i=1,C%elnum
00580         if (C%els(i)%nfs/=0) then
00581             do k=1,2**M%nsd
00582                 ctiremap(C%els(i)%n(k))=1
00583             enddo
00584         endif
00585         enddo
00586 
00587         ! Find the number of unused grid nodes
00588         pcnt=count(ctiremap==0)-1; C%nct=C%nct-pcnt
00589 
00590         if (sctls%verbose>3) &
00591                 write (stream,*) "Removing ",pcnt," empty coarse elements"
00592 
00593         allocate(ncoords(M%nsd,C%nct))
00594 
00595         if (pcnt>0) then
00596 
00597             ! Create the remap and fill the initial node part of ncoords
00598             do i=1,C%ncti
00599             if (ctiremap(i)==1) then
00600                 ctiremap(i)=ctiremap(i-1)+1
00601                 ncoords(:,ctiremap(i))=C%coords(:,i)
00602             else
00603                 ctiremap(i)=ctiremap(i-1)
00604             endif
00605             enddo
00606 
00607             ! Change the initial element count accordingly
00608             C%ncti=C%ncti-pcnt
00609 
00610             ! Walk the elements and change their indices as required
00611             do i=1,C%elnum
00612             if (C%els(i)%nfs>0) then
00613                 C%els(i)%n=ctiremap(C%els(i)%n(:))
00614             else
00615                 C%els(i)%n=0 ! To make it explicitly seem useless
00616             endif
00617             enddo
00618        else
00619             ! Simply copy initial nodes over
00620             ncoords(:,1:C%ncti)=C%coords(:,1:C%ncti)
00621         endif
00622 
00623         ! Deallocate ctiremap
00624         deallocate(ctiremap)
00625 
00626         ! Move the hanging nodes over as they are
00627         if(C%nhn>0) then
00628            ncoords(:,C%nct-C%nhn+1:C%nct)=C%coords(:,hangpt+1:size(C%coords))
00629         end if
00630 
00631         ! Walk the hanging nodes through if needed
00634         if (hangpt/=C%ncti+C%refnum) then
00635         pcnt=hangpt-(C%ncti+C%refnum)
00636         do i=1,C%refnum
00637             if (associated(C%refels(i)%hnds)) then
00638                 C%refels(i)%hnds=max(0,C%refels(i)%hnds-pcnt)
00639             endif
00640         enddo
00641         endif
00642  
00643         ! Walk the initial grid elements and add their refined elems
00644         k=C%ncti+1
00645         do i=1, C%elnum
00646            C%els(i)%nref=k
00647            
00648            ! Walk the refined elements belonging to it
00649            j=C%els(i)%rbeg
00650            do while (j/=-1)
00651                ncoords(:,k)=C%coords(:,C%refels(j)%node)
00652                C%refels(j)%node=k
00653                k=k+1; j=C%refels(j)%next
00654            enddo
00655 
00656            C%els(i)%nref=k-C%els(i)%nref
00657         enddo
00658 
00659         ! Deallocate old coords
00660         deallocate(C%coords)
00661         
00662         ! Put in the new ones
00663         C%coords=>ncoords
00664  
00665         ! Print a small report
00666         if (sctls%verbose>3) then
00667                 write (stream,*) "Global coarse mesh has: "
00668                 write (stream,*) "Grid nodes:   ",C%ncti
00669                 write (stream,*) "Refinements:  ",C%refnum
00670                 write (stream,*) "Hanging nodes:",C%nhn
00671         endif
00672 
00673 end subroutine CreateCoarseMesh
00674 
00675 subroutine CreateHangingNodes(refpt,coordpt,nsd,nsame,minv,maxv,C)
00676     use CoarseGrid_class
00677     implicit none
00678     integer, intent(in) :: refpt
00679     integer, intent(inout) :: coordpt
00680     integer, intent(in) :: nsd
00681     integer, intent(in) :: nsame(:)
00682     real(kind=xyzk),intent(in) :: minv(:),maxv(:)
00683     type(CoarseGrid), intent(inout) :: C
00684 
00685     integer :: i,j,k,nb
00686     type(RefinedElem),pointer :: new, cur
00687     real(kind=xyzk),pointer :: pt(:)
00688     real(kind=xyzk) :: pt2(nsd)
00689 
00690     new=>C%refels(refpt); pt=>C%coords(:,new%node); k=1
00691     do i=1,nsd
00692         !************************************
00693         ! Positive direction
00694         !************************************
00695 
00696         ! Check if we have the room                    
00697         if (C%ncti+refpt>=coordpt) exit
00698 
00699         ! Locate the neighbour in positive direction
00700         nb=getNeighbourEl(refpt,k,nsd,nsame,C)
00701         !write(stream,*) i,"> ",refpt,nb
00702                     
00703         ! Create the new nodes coordinates
00704         j=-1
00705         if (nb>0) then ! use the average of the two centers
00706         if (C%coords(i,C%refels(nb)%node)<maxv(i)) write(stream,*) "SCREAM!!!",i
00707 
00708             if (C%refels(nb)%level==new%level)then ! assuming we can
00709                pt2=(C%coords(:,C%refels(nb)%node)+pt)/2.0_xyzk; j=nb
00710             endif
00711         else if (nb==0) then
00712             j=0; pt2=pt ! and use the current element center
00713         endif                     
00714         pt2(i)=maxv(i) ! put it on the dividing line ( moreorless)
00715 
00716                     
00717         if (j>=0) then ! we do actually need to add one
00718             ! Add the node to coordinates list
00719             C%coords(:,coordpt)=pt2
00720 
00721             ! Allocate memory for the adressing
00722             if (.not.associated(new%hnds)) then
00723                 allocate(new%hnds(2*nsd)); new%hnds=0
00724             endif
00725 
00726             ! And add this node to the map
00727             new%hnds(i)=coordpt
00728 
00729             ! Same stuff for the other node aswell
00730             if (j>0) then
00731                 cur=>C%refels(j)
00732                 if (.not.associated(cur%hnds)) then
00733                     allocate(cur%hnds(2*nsd)); cur%hnds=0
00734                 endif
00735 
00736                 cur%hnds(nsd+i)=coordpt
00737             endif
00738 
00739             ! Decrease the coordinate adress
00740             coordpt=coordpt-1
00741         endif
00742 
00743         !************************************
00744         ! Negative direction - completely analogous
00745         !************************************
00746 
00747         ! Check if we have the room                    
00748         if (C%ncti+refpt>=coordpt) exit
00749 
00750         ! Locate the neighbour in negative direction
00751         nb=getNeighbourEl(refpt,-k,nsd,nsame,C)
00752                     
00753         ! Create the new nodes coordinates
00754         j=-1
00755         if (nb>0) then ! use the average of the two centers
00756             if (C%refels(nb)%level==new%level)then ! assuming we can
00757                pt2=(C%coords(:,C%refels(nb)%node)+pt)/2.0_xyzk; j=nb
00758             endif
00759         else if (nb==0) then
00760             j=0; pt2=pt ! and use the current element center
00761         endif                    
00762         pt2(i)=minv(i) ! put it on the dividing line ( moreorless)
00763                     
00764         if (j>=0) then ! we do actually need to add one
00765             ! Add the node to coordinates list
00766             C%coords(:,coordpt)=pt2
00767 
00768             ! Allocate memory for the adressing
00769             if (.not.associated(new%hnds)) then
00770                 allocate(new%hnds(2*nsd)); new%hnds=0
00771             endif
00772 
00773             ! And add this node to the map
00774             new%hnds(nsd+i)=coordpt
00775 
00776             ! Same stuff for the other node aswell
00777             if (j>0) then
00778                 cur=>C%refels(j)
00779                 if (.not.associated(cur%hnds)) then
00780                     allocate(cur%hnds(2*nsd)); cur%hnds=0
00781                 endif
00782                 cur%hnds(i)=coordpt
00783             endif
00784 
00785             ! Decrease the coordinate adress
00786             coordpt=coordpt-1
00787 
00788         endif
00789 
00790         k=2*k
00791     enddo
00792 
00793 end subroutine CreateHangingNodes
00794 
00796 subroutine CreateCoarseFreemap(C,M)
00797         use RealKind
00798         use CoarseGrid_class
00799         
00800         implicit none
00801 
00802         type(CoarseGrid), intent(inout) :: C !< The Coarse Grid to generate it for.
00803 
00804         type(Mesh), intent(in) :: M
00805 
00806         integer :: i
00807 
00808         ! As there is currently one freedom for each node..
00809         C%ngfc=C%nct; C%nlfc=C%ngfc
00810 
00811         call CoarseGrid_allocate(C,cfreemap=.true.)
00812        
00813         do i=1, C%ngfc
00814             C%cfreemap(i)=i
00815         enddo
00816 
00817 end subroutine CreateCoarseFreemap
00818 
00819 ! Some choices for choosing the center
00820 
00822 subroutine ChooseGeometricCenter(pt,cpt,pts,elmap,el,refels,flags,minv,maxv)
00823     use RealKind
00824     use CoarseGrid_class
00825 
00826     implicit none
00827     
00828     real(kind=xyzk), intent(out) :: pt(:) !< the new center to output
00829     real(kind=xyzk), intent(in) :: cpt(:) !< coordinates of the prev cn
00830     real(kind=xyzk), intent(in) :: pts(:,:) !< points array
00831     integer, intent(in) :: elmap(:) !< indices to pts
00832     integer, intent(in) :: el !< index of the element being divided
00833     type(RefinedElem), intent(in) :: refels(:) !< refined elements 
00834     integer, intent(in) :: flags !< flags saying which dir to div
00835     real(kind=xyzk), intent(in) :: minv(:),maxv(:) !< elem bounds
00836 
00837     real(kind=xyzk) :: h(size(minv))
00838     integer :: i, fl
00839 
00840     ! Simply use the geometric center
00841     pt=(maxv+minv)/(2.0_xyzk)
00842 
00843 end subroutine ChooseGeometricCenter
00844 
00846 subroutine ChooseMeanCenter(pt,cpt,pts,elmap,el,refels,flags,minv,maxv)
00847     use RealKind
00848     use CoarseGrid_class
00849 
00850     implicit none
00851 
00852     real(kind=xyzk), intent(out) :: pt(:) ! the new center to output
00853     real(kind=xyzk), intent(in) :: cpt(:) ! coordinates of the prev cn
00854     real(kind=xyzk), intent(in) :: pts(:,:) ! points array
00855     integer, intent(in) :: elmap(:) ! indices to pts
00856     integer, intent(in) :: el ! index of the element being divided
00857     type(RefinedElem), intent(in) :: refels(:) ! refined elements 
00858     integer, intent(in) :: flags ! flags saying which dir to div
00859     real(kind=xyzk), intent(in) :: minv(:),maxv(:)
00860 
00861     integer :: i,j
00862     real(kind=xyzk) :: bval
00863 
00864     if (meanpow==0.0_xyzk) then ! Geometric mean
00865         pt=1.0_xyzk
00866         
00867         ! Multiply all the coordinates together
00868         do i=1,size(elmap)
00869             pt=pt*pts(:,elmap(i))  
00870         enddo
00871 
00872         ! Take the appropriate root
00873         pt=pt**(1.0_xyzk/size(elmap))
00874     else if (meanpow==1.0_xyzk) then ! we can avoid costly exponentations
00875         pt=0.0_xyzk
00876 
00877         ! Find the sum
00878         do i=1,size(elmap)
00879             pt=pt+pts(:,elmap(i))
00880         enddo
00881 
00882         ! And divide by count
00883         pt=pt/size(elmap)
00884     else ! Some other mean
00885         pt=0.0_xyzk
00886 
00887         ! Add the coordinates together
00888         do i=1,size(elmap)
00889             pt=pt+(pts(:,elmap(i))**meanpow)
00890         enddo
00891 
00892         ! Get their average
00893         pt=(pt/size(elmap))**(1.0_xyzk/meanpow)
00894 
00895     endif
00896 
00897     ! Degeneracy avoidance (until a better way of doing this is found)
00898     do i=1,size(pt)
00899         bval=0.25_xyzk*(maxv(i)-minv(i))
00900         if (pt(i)>maxv(i)-bval) then; pt(i)=maxv(i)-bval
00901         elseif (pt(i)<minv(i)+bval) then; pt(i)=minv(i)+bval
00902         endif
00903     enddo
00904  
00905 end subroutine ChooseMeanCenter
00906 
00911 subroutine ChooseMeridianCenter(pt,cpt,pts,elmap,el,refels,flags,minv,maxv)
00912     use RealKind
00913     use CoarseGrid_class
00914 
00915     implicit none
00916 
00917     real(kind=xyzk), intent(out) :: pt(:) !< the new center to output
00918     real(kind=xyzk), intent(in) :: cpt(:) !< coordinates of the prev cn
00919     real(kind=xyzk), intent(in) :: pts(:,:) !< points array
00920     integer, intent(in) :: elmap(:) !< indices to pts
00921     integer, intent(in) :: el !< index of the element being divided
00922     type(RefinedElem), intent(in) :: refels(:) !< refined elements 
00923     integer, intent(in) :: flags !< flags saying which dir to div
00924     real(kind=xyzk), intent(in) :: minv(:),maxv(:)
00925 
00926     integer :: i,k,j,c,frml,tol
00927     real(kind=xyzk) :: bval
00928     integer :: sz, get
00929     integer :: ar(size(elmap))
00930 
00931     sz=size(elmap); get=sz/2
00932 
00933     ar=elmap
00934 
00935     do c=1,size(pt)
00936     
00937     frml=1; tol=sz
00938 
00939     if (sz==0) call DOUG_abort("Coarse grid too fine, adjust parameters")
00940 
00941     do while (frml<tol)
00942         ! Choose a pivot at random
00943         call random_number(bval)
00944         j=floor(bval*(tol-frml+1))+frml
00945         bval=pts(c,ar(j))
00946         
00947         ! Move larger to one side and smaller to the other
00948         i=frml; k=tol
00949         do
00950             ! Find one in the first part that should be in second
00951             do i=i,k  
00952             if ( bval<=pts(c,ar(i)) ) exit
00953             enddo
00954 
00955             ! Find one in the second part that should be in first
00956             do k=k,i,-1  
00957             if ( bval>pts(c,ar(k)) ) exit
00958             enddo
00959 
00960             ! If the ends meet, stop
00961             if (i>=k) exit ! Note i=k+1 is also possible
00962 
00963             ! Otherwise switch them
00964             j=ar(k); ar(k)=ar(i); ar(i)=j
00965         enddo
00966    
00967         ! Set the new bounds
00968         if (get==i) then
00969             exit
00970         else if (get>i) then
00971             do frml=i,get
00972                if (pts(c,ar(frml))/=pts(c,ar(frml+1))) exit
00973             enddo
00974         else
00975             tol=i-1
00976         endif
00977     enddo
00978     
00979     ! Set it as the result coordinate
00980     pt(c)=pts(c,ar(frml))
00981     
00982     enddo
00983     
00984     ! Degeneracy avoidance (until a better way of doing this is found)
00985     do i=1,size(pt)
00986         bval=0.25_xyzk*(maxv(i)-minv(i))
00987         if (pt(i)>maxv(i)-bval) then; pt(i)=maxv(i)-bval
00988         elseif (pt(i)<minv(i)+bval) then; pt(i)=minv(i)+bval
00989         endif
00990     enddo
00991     
00992 end subroutine ChooseMeridianCenter
00993 
00994 end module CreateCoarseGrid