|
DOUG 0.2
|
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
1.7.3-20110217