|
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 00022 module Mesh_class 00023 00024 use DOUG_utils 00025 use globals 00026 use Graph_class 00027 use Polygon_class 00028 00029 implicit none 00030 00031 #include<doug_config.h> 00032 00033 integer, parameter :: D_FREEDOM_INNER = 0 00034 integer, parameter :: D_FREEDOM_INTERF = 1 00035 00036 !-------------------------------------------------------------- 00037 !! Mesh type 00038 !-------------------------------------------------------------- 00039 type Mesh 00040 00041 !! Mesh parameters: 00042 integer :: nell = -1 !< Number of elements in a mesh 00043 integer :: ngf = -1 !< Number of global freedoms 00044 integer :: mfrelt = -1 !< Max number of freedoms in element 00045 integer :: nsd = -1 !< Number of spacial dimensions 00046 integer :: nnode = -1 !< Number of nodes 00047 00048 ! Main mesh data: 00050 integer, dimension(:), pointer :: nfrelt 00051 00052 integer, dimension(:,:), pointer :: mhead 00053 00055 integer(kind=1), dimension(:), pointer :: freemask 00056 00057 integer, dimension(:), pointer :: freemap 00058 00060 real(kind=xyzk), dimension(:,:), pointer :: coords 00061 00063 integer :: lnnode = -1 00064 real(kind=xyzk), dimension(:,:), pointer :: lcoords 00065 integer, dimension(:), pointer :: lfreemap 00066 00067 00068 ! Auxiliary mesh data: 00069 ! Hash table 00070 00071 ! Try to implement it as linked list(?) The need for mesh 00072 ! assessment will disappear. 00073 ! It will be possible to add elements to the list in a natural way. 00074 ! 00075 ! (variables needed to be sent to slaves) 00076 ! Node hash table freedom<->element : hash[hashsize,2] 00077 ! hash(:,1) - freedoms, hash(:,2) - elements 00078 !! hash['freedom number', 'element number'] 00079 integer, dimension(:,:), pointer :: hash 00080 !! Size of hash table 00081 integer :: hashsize = -1 00082 !! auxiliary parameters 00083 integer, dimension(:), pointer :: hashlook 00084 integer :: numbins = 0 00085 integer :: binsize = 0 00086 integer :: hashentries = 0 00087 real(kind=rk) :: nperfree = 0.0_rk 00088 real(kind=rk) :: nperstd = 0.0_rk 00089 real(kind=rk) :: hscale = 0.0_rk 00090 00091 !! Partition 00092 logical :: parted = .false. 00093 integer :: nparts = -1 00094 !! Number of local freedoms (calculated locally) 00095 integer :: nlf = -1 00096 ! NB! change to(?) nlell - number of local elements : 00097 !! Number of elements in each partition: partnelems[nparts] 00098 integer, dimension(:), pointer :: partnelems 00099 00101 integer, dimension(:), pointer :: eptnmap 00102 !! Inner/interface mask for (local) freedoms : 00103 !! inner_interf_fmask[nlf] - '0' - inner, '1' - interf. 00104 integer, dimension(:), pointer :: inner_interf_fmask 00105 00106 ! MOVE IT FROM HERE => MAKE GLOBAL 00108 integer :: nnghbrs 00109 00110 integer, dimension(:), pointer :: nghbrs 00111 ! <= 00112 00115 integer, dimension(:), pointer :: nfreesend_map 00116 integer, dimension(:), pointer :: nghostsend_map 00117 ! Mappings 00119 integer, dimension(:), pointer :: gl_fmap 00120 00121 integer, dimension(:), pointer :: lg_fmap 00122 ! Data structures for assembled matrices case with its particluar 00123 ! overlap size: 00124 integer :: ntobsent !< #inner freedoms to be sent to neighbours 00125 integer :: ninonol !< #inner freedoms that are on overlap 00126 integer :: ninner !< #inner freedoms total (totally inner:nino) 00127 integer :: indepoutol !< points to the end of freedoms on outer 00128 !! overlap that does not get comm during Ax-op. 00129 !! ie. freedoms indepoutol+1:nlf get value 00130 !! through comm in Ax-op. 00131 00132 ! we are organising local freedoms as follows: 00133 00134 ! ('ol' - overlap) 00135 !1,2,...,M%ntobsent,...,M%ninonol,...,M%ninner,...,M%indepoutol,...,M%nlf| 00136 !<-feedoms4send -> |<-rest inol->|<-independ.>|<-indep.onoutol>|<receivd>| 00137 !<- inner overlap -> |<-freedoms->|<- outer overlap ->| 00138 !<- all inner freedoms ->| 00139 00140 ! if ol=0 then 00141 !<-feedoms4send -> |<-rest inol->|<-independ.>|<-indep.onoutol>|<receivd>| 00142 ! | / | / | 00143 ! | . . . . . | . . . . . . . | 00144 ! |/ |/ | 00145 ! actual freedoms | ghost freedoms | 00146 00149 type(indlist),dimension(:),pointer :: ax_recvidx,ax_sendidx 00150 type(indlist),dimension(:),pointer :: ol_outer !< overlap with each neighbour process in its region 00151 type(indlist),dimension(:),pointer :: ol_inner !< overlap with each neighbour process in my region 00152 type(indlist),dimension(:),pointer :: ol_solve !< overlap with each neighbour process in all regions (its,my,or third) 00153 00154 00155 !! Graph 00156 type(Graph) :: G 00157 00158 end type Mesh 00159 00160 ! Private methods 00161 private :: & 00162 Mesh_findNLF, & 00163 Mesh_findNghbrs, & 00164 Mesh_MapsAndNghbrs, & 00165 Mesh_assessHash, & 00166 hashfn 00167 00168 contains 00169 00170 00171 !---------------------------- 00172 !! Constructor 00173 !---------------------------- 00174 function Mesh_New() result(M) 00175 implicit none 00176 00177 type(Mesh) :: M 00178 00179 M%mhead => NULL() 00180 00181 M%nell = -1 00182 M%ngf = -1 00183 M%nsd = -1 00184 M%mfrelt = -1 00185 M%nnode = -1 00186 M%lnnode = -1 00187 M%nnghbrs = 0 00188 00189 nullify(M%nfrelt) 00190 nullify(M%mhead) 00191 nullify(M%coords) 00192 nullify(M%hash) 00193 nullify(M%hashlook) 00194 nullify(M%freemap) 00195 nullify(M%eptnmap) 00196 nullify(M%freemask) 00197 nullify(M%inner_interf_fmask) 00198 nullify(M%gl_fmap) 00199 nullify(M%lg_fmap) 00200 nullify(M%nfreesend_map) 00201 nullify(M%partnelems) 00202 nullify(M%nghbrs) 00203 nullify(M%lcoords) 00204 nullify(M%lfreemap) 00205 00206 M%G = Graph_New() 00207 00208 end function Mesh_New 00209 00210 00211 !----------------------------------------------------- 00212 !! Initialiser 00213 !----------------------------------------------------- 00214 subroutine Mesh_Init(M, nell, ngf, nsd, mfrelt, nnode) 00215 implicit none 00216 00217 integer, intent(in) :: nell, ngf, nsd, mfrelt, nnode 00218 type(Mesh), intent(in out) :: M ! Mesh (not filled, not allocated) 00219 00220 M%nell = nell 00221 M%ngf = ngf 00222 M%nsd = nsd 00223 M%mfrelt = mfrelt 00224 M%nnode = nnode 00225 M%lnnode = nnode 00226 00227 end subroutine Mesh_Init 00228 00229 00230 !------------------------------------------------------------- 00231 !! Constructor and initialiser 00232 !------------------------------------------------------------- 00233 function Mesh_newInit(nell, ngf, nsd, mfrelt, nnode) result(M) 00234 implicit none 00235 00236 integer, intent(in) :: nell, ngf, nsd, mfrelt, nnode 00237 type(Mesh) :: M ! Mesh (not filled, just allocated) 00238 00239 M = Mesh_New() 00240 call Mesh_Init(M, nell, ngf, nsd, mfrelt, nnode) 00241 00242 end function Mesh_newInit 00243 00244 !------------------------- 00245 !! Destructor 00246 !------------------------- 00247 subroutine Mesh_Destroy(M) 00248 implicit none 00249 00250 type(Mesh), intent(in out) :: M ! Mesh 00251 00252 if (associated(M%nfrelt)) deallocate(M%nfrelt) 00253 if (associated(M%mhead)) deallocate(M%mhead) 00254 if (associated(M%coords)) deallocate(M%coords) 00255 if (associated(M%hash)) deallocate(M%hash) 00256 if (associated(M%hashlook)) deallocate(M%hashlook) 00257 if (associated(M%freemap)) deallocate(M%freemap) 00258 if (associated(M%eptnmap)) deallocate(M%eptnmap) 00259 if (associated(M%freemask)) deallocate(M%freemask) 00260 if (associated(M%inner_interf_fmask)) deallocate(M%inner_interf_fmask) 00261 if (associated(M%gl_fmap)) deallocate(M%gl_fmap) 00262 if (associated(M%lg_fmap)) deallocate(M%lg_fmap) 00263 if (associated(M%nfreesend_map)) deallocate(M%nfreesend_map) 00264 if (associated(M%partnelems)) deallocate(M%partnelems) 00265 if (associated(M%nghbrs)) deallocate(M%nghbrs) 00266 if (.NOT.ismaster().AND.associated(M%lcoords)) deallocate(M%lcoords) ! for master it = coords 00267 if (associated(M%lfreemap)) deallocate(M%lfreemap) 00268 !if (associated(M%)) deallocate(M%) 00269 00270 M%nell = -1 00271 M%ngf = -1 00272 M%nsd = -1 00273 M%mfrelt = -1 00274 M%nnode = -1 00275 M%lnnode = -1 00276 00277 call Mesh_destroyGraph(M) 00278 00279 end subroutine Mesh_Destroy 00280 !======================================== 00281 ! 00282 ! File I/O 00283 ! 00284 !---------------------------------------- 00293 subroutine Mesh_initFromFile(M, fnInfo) 00294 use globals, only: mctls, sctls ! To correct ENTWIFE's inexact value of 00295 ! 'ngf' for multivariable problems 00296 implicit none 00297 00298 type(Mesh), intent(in out) :: M ! Mesh 00299 character*(*) :: fnInfo 00300 00301 integer(4) :: tmp(5) 00302 integer :: nell, ngf, nsd, mfrelt, nnode 00303 00304 00305 open(50, FILE=fnInfo, STATUS='OLD', FORM='UNFORMATTED') !XXX TODO Infofile 00306 read (50) tmp 00307 nell=tmp(1); ngf=tmp(2); nsd=tmp(3); mfrelt=tmp(4); nnode = tmp(5) 00308 close(50) 00309 00310 ! To correct ENTWIFE's inexact value of 'ngf' (which comes 00311 ! from 'info'-file) for multivariable problems 00312 if (sctls%number_of_blocks > 1) then 00313 ! Initialise with WRONG 'ngf' value at first, assuming all 00314 ! the others are correct 00315 call Mesh_Init(M, nell, ngf, nsd, mfrelt, nnode) 00316 call Mesh_readFileFreelists(M, trim(mctls%freedom_lists_file)) 00317 ngf = maxval(M%mhead) 00318 end if 00319 00320 call Mesh_Init(M, nell, ngf, nsd, mfrelt, nnode) 00321 00322 end subroutine Mesh_initFromFile 00323 00324 00325 !------------------------------------------------ 00328 subroutine Mesh_BuildSquare(M,n) 00329 implicit none 00330 type(Mesh), intent(in out) :: M ! Mesh 00331 integer, intent(in) :: n 00332 integer :: nell,ngf,nsd,mfrelt,nnode,i,j,k 00333 00334 ngf=n*n 00335 nell=0 ! just for now 00336 nsd=2 00337 mfrelt=4 00338 nnode=ngf 00339 call Mesh_Init(M, nell, ngf, nsd, mfrelt, nnode) 00340 allocate(M%coords(M%nsd,M%nnode),M%freemap(M%ngf)) ! needing the coords... 00341 M%freemap= (/ (i,i=1,ngf) /) 00342 do j=1,n 00343 do i=1,n 00344 k=(j-1)*n+i 00345 M%coords(1,k)=1.0_xyzk*(i-1)/(n-1) 00346 M%coords(2,k)=1.0_xyzk*(j-1)/(n-1) 00347 enddo 00348 enddo 00349 M%nlf=ngf 00350 end subroutine Mesh_BuildSquare 00351 00352 !------------------------------------------------ 00355 function Mesh_newInitFromFile(fnInfo) result(M) 00356 implicit none 00357 00358 character*(*) :: fnInfo 00359 type(Mesh) :: M ! Mesh 00360 00361 M = Mesh_New() 00362 call Mesh_initFromFile(M, fnInfo) 00363 00364 end function Mesh_newInitFromFile 00365 00366 00367 !-------------------------------- 00370 subroutine Mesh_readFromFile(M, & 00371 fnFreelists, & 00372 fnCoords, & 00373 fnFreemap, & 00374 fnFreemask) 00375 00376 use globals, only : stream,& 00377 sctls ! To correct ENTWIFE's inexact value of 'ngf' for 00378 ! multivariable problems 00379 00380 implicit none 00381 00382 type(Mesh), intent(in out) :: M ! Mesh 00383 character*(*), optional :: fnFreelists, fnCoords, fnFreemap, fnFreemask 00384 00385 logical :: key=.false. 00386 00387 00388 if (present(fnFreelists)) then 00389 if (sctls%number_of_blocks <= 1) then ! Otherwise freedoms lists were 00390 ! already read in to correct 00391 ! ENTWIFE's error 00392 call Mesh_readFileFreelists(M, fnFreelists) 00393 key = .true. 00394 end if 00395 end if 00396 if (present(fnCoords).and.& 00397 (sctls%plotting == D_PLOT_YES .or. sctls%levels == 2 )) then 00398 ! Load nodes coordinates only if we want to plot a mesh. 00399 ! NB! We would need them also for geometric multilevel methods at 00400 ! construction of coarse meshes. 00401 call Mesh_readFileCoords(M, fnCoords) 00402 key = .true. 00403 end if 00404 if (present(fnFreemap)) then 00405 call Mesh_readFileFreemap(M, fnFreemap) 00406 key = .true. 00407 end if 00408 if (present(fnFreemask)) then 00409 call Mesh_readFileFreemask(M, fnFreemask) 00410 key = .true. 00411 end if 00412 00413 if (.not.key) then 00414 call DOUG_abort('[Mesh_readFromFile] : no files specified for reading') 00415 end if 00416 00417 end subroutine Mesh_readFromFile 00418 00419 00420 !-------------------------------------------- 00425 subroutine Mesh_readFileFreelists(M, fnFreelists) 00426 00427 use globals, only : stream 00428 00429 implicit none 00430 00431 type(Mesh), intent(in out) :: M ! Mesh 00432 character*(*) :: fnFreelists 00433 00434 integer :: i 00435 integer(4) :: nfrelt, mhead(M%mfrelt) 00436 00437 ! Check for errors 00438 if (M%nell <= 0) & 00439 call DOUG_abort('[Mesh_readFileFreelists] : Mesh object must be'//& 00440 ' initialised first.', -1) 00441 00442 if (.not.associated(M%nfrelt)) allocate(M%nfrelt(M%nell)) 00443 if (.not.associated(M%mhead)) allocate(M%mhead(M%mfrelt,M%nell)) 00444 00445 ! Zero array 00446 M%mhead = 0 00447 00448 ! Read in element nodes data 00449 write(stream, FMT='(a)', advance='no') 'Reading in element nodes'//& 00450 ' numbering ... ' 00451 open(50, FILE=fnFreelists, STATUS='OLD', FORM='UNFORMATTED') !XXX TODO freedom_list_file 00452 do i = 1,M%nell 00453 read (50) nfrelt, mhead(1:nfrelt) 00454 M%nfrelt(i) = nfrelt 00455 M%mhead(1:M%nfrelt(i),i) = mhead(1:M%nfrelt(i)) 00456 if (maxval(M%mhead(:,i))>M%ngf) then 00457 write (*, *) i, M%ngf, maxval(M%mhead(:,i)) 00458 call DOUG_abort('[Mesh_readFileFreelists] - freedom index out of range', -1) 00459 end if 00460 enddo 00461 00462 close(50) 00463 write(stream, *) 'done' 00464 00465 end subroutine Mesh_readFileFreelists 00466 00467 00468 !------------------------------------------ 00472 subroutine Mesh_readFileCoords(M, fnCoords) 00473 00474 use globals, only : stream 00475 00476 implicit none 00477 00478 type(Mesh), intent(in out) :: M ! Mesh 00479 character*(*) :: fnCoords 00480 00481 if ((M%nsd <= 0).and.(M%nnode <= 0)) & 00482 call DOUG_abort('[Mesh_readFileCoords] : Mesh object must be'//& 00483 ' initialised first.',-1) 00484 00485 if (.not.associated(M%coords)) allocate(M%coords(M%nsd,M%nnode)) 00486 00487 ! Read in node coordinates 00488 write(stream, FMT='(a)', advance='no') 'Reading in nodes coordinates ... ' 00489 open(50, FILE=fnCoords, STATUS='OLD', FORM='UNFORMATTED') !XXX TODO coords_file 00490 read (50) M%coords 00491 close(50) 00492 write(stream, *) 'done' 00493 00494 end subroutine Mesh_readFileCoords 00495 00496 00497 !-------------------------------------------- 00501 subroutine Mesh_readFileFreemap(M, fnFreemap) 00502 00503 use globals, only : stream 00504 00505 implicit none 00506 00507 type(Mesh), intent(in out) :: M ! Mesh 00508 character*(*) :: fnFreemap 00509 integer(4) :: freemap(M%ngf) 00510 00511 if (M%ngf <= 0) & 00512 call DOUG_abort('[Mesh_readFileFreemap] : Mesh object must be'//& 00513 ' initialised first.',-1) 00514 00515 if (.not.associated(M%freemap)) allocate(M%freemap(M%ngf)) 00516 00517 write(stream, FMT='(a)', advance='no') 'Reading in freedoms'' map ... ' 00518 open(50, FILE=fnFreemap, STATUS='OLD', FORM='UNFORMATTED') !XXX TODO freemap_file 00519 read (50) freemap 00520 M%freemap = freemap 00521 close(50) 00522 write(stream, *) 'done' 00523 00524 end subroutine Mesh_readFileFreemap 00525 00526 00527 !-------------------------------------------------------- 00531 subroutine Mesh_readFileFreemask(M, fnFreemask) 00532 00533 use globals, only : stream 00534 00535 implicit none 00536 00537 type(Mesh), intent(in out) :: M ! Mesh 00538 character*(*) :: fnFreemask 00539 00540 if (M%ngf <= 0) & 00541 call DOUG_abort('[Mesh_readFileFreemask] : Mesh object must be'//& 00542 ' initialised first.',-1) 00543 00544 if (.not.associated(M%freemask)) allocate(M%freemask(M%ngf)) 00545 00546 write(stream, FMT='(a)', advance='no') 'Reading in freedoms'' mask ... ' 00547 open(50, FILE=fnFreemask, STATUS='OLD', FORM='UNFORMATTED') !XXX TODO freedom_mask_file 00548 read (50) M%freemask 00549 close(50) 00550 write(stream, *) 'done' 00551 00552 end subroutine Mesh_readFileFreemask 00553 !================================== 00554 ! 00555 !! Graph creation methods 00556 ! 00557 !---------------------------------- 00558 !! Builds mesh dual graph 00559 !---------------------------------- 00560 subroutine Mesh_buildGraphDual(M) 00561 00562 use globals, only : stream 00563 00564 implicit none 00565 00566 type(Mesh), intent(in out) :: M ! Mesh 00567 00568 integer :: binsize, numbins, hashsize 00569 real(kind=rk) :: nperfree, nperstd, hscale 00570 00571 integer, dimension(:), pointer :: xadj 00572 integer, dimension(:), pointer :: adjncy 00573 integer :: nedges 00574 00575 00576 call Mesh_buildHash(M) 00577 call Mesh_buildDualAdjncy(M, nedges, xadj, adjncy) 00578 00579 ! Create graph object representing mesh's dual graph 00580 M%G = Graph_newInit(M%nell, nedges, xadj, adjncy, D_GRAPH_DUAL) 00581 00582 ! Deallocate temporary arrays 00583 if (associated(xadj)) deallocate(xadj) 00584 if (associated(adjncy)) deallocate(adjncy) 00585 00586 end subroutine Mesh_buildGraphDual 00587 00588 00589 !----------------------------------- 00590 !! Destroy graph associated with mesh 00591 !----------------------------------- 00592 subroutine Mesh_destroyGraph(M) 00593 implicit none 00594 type(Mesh), intent(in out) :: M ! Mesh 00595 00596 call Graph_Destroy(M%G) 00597 end subroutine Mesh_destroyGraph 00598 00599 00600 !---------------------------- 00601 !! Assess mesh connectivity 00602 !---------------------------- 00603 subroutine Mesh_assessHash(M) 00604 00605 use globals, only : stream, D_MSGLVL 00606 00607 implicit none 00608 00609 type(Mesh), intent(in out) :: M ! Mesh 00610 real(kind=rk) :: tmp 00611 00612 integer :: i, j, nactivef 00613 integer, dimension(:), pointer :: freecount 00614 00615 00616 allocate(freecount(M%ngf)) 00617 freecount = 0 00618 00619 do i = 1,M%nell 00620 do j = 1,M%nfrelt(i) 00621 freecount(M%mhead(j,i)) = freecount(M%mhead(j,i)) + 1 00622 enddo 00623 enddo 00624 00625 nactivef = 0 00626 do i = 1,M%ngf 00627 if (freecount(i) > 0) then 00628 nactivef = nactivef + 1 00629 endif 00630 enddo 00631 00632 ! Average connection 00633 M%nperfree = sum(M%nfrelt) / nactivef 00634 00635 do i = 1,M%ngf 00636 if (freecount(i) > 0) then 00637 tmp = freecount(i) - M%nperfree 00638 M%nperstd = M%nperstd + tmp*tmp 00639 endif 00640 enddo 00641 00642 ! Standard deviation 00643 M%nperstd = sqrt(M%nperstd/(nactivef-1)) 00644 00645 00646 ! connectivity data 00647 !!$ M%numbins = int((M%nperfree+M%nperstd)*M%ngf/100.0 + 0.95) + 1 00648 M%numbins = int((M%nperfree+M%nperstd)*nactivef/100.0 + 0.95) + 1 00649 M%binsize = 101 00650 M%hashsize = M%numbins*M%binsize !- 1 00651 00652 M%hscale = 1.0*100/(M%nperfree+M%nperstd) 00653 00654 if (D_MSGLVL >= 1) then 00655 write(stream, *) 00656 write(stream, *) 'Statistics of mesh assessment:' 00657 write(stream, '(a,i9)') ' nactivef = ', nactivef 00658 write(stream, '(a,i9)') ' numbins = ', M%numbins 00659 write(stream, '(a,i9)') ' hashsize = ', M%hashsize 00660 write(stream, '(a,f10.5)') ' hscale = ', M%hscale 00661 write(stream, '(a,f10.5)') ' Average connection :', M%nperfree 00662 write(stream, '(a,f10.5)') ' Standard deviation :', M%nperstd 00663 call flush(stream) 00664 end if 00665 00666 deallocate(freecount) 00667 00668 end subroutine Mesh_assessHash 00669 00670 00671 !--------------------------- 00672 !! Builds node hash table 00673 !--------------------------- 00674 subroutine Mesh_buildHash(M) 00675 use globals, only : stream 00676 implicit none 00677 00678 type(Mesh), intent(in out) :: M ! Mesh 00679 00680 integer :: i, j, h, v1 00681 integer :: hashsize 00682 integer, dimension(:,:), allocatable :: hashtmp 00683 00684 00685 ! Asses mesh to build hash table of nodes 00686 call Mesh_assessHash(M) 00687 00688 allocate(M%hash(M%hashsize,2)) 00689 M%hash = 0 00690 00691 allocate(M%hashlook(M%numbins)) 00692 do i = 1,M%numbins 00693 M%hashlook(i) = M%binsize*(i-1) + 1 00694 enddo 00695 00696 ! Actually build hash table 00697 M%hashentries = 0 00698 do i = 1,M%nell 00699 ! loop over nodes in element 00700 do j = 1,M%nfrelt(i) 00701 v1 = M%mhead(j,i) 00702 ! check hash table full 00703 if (M%hashentries == M%hashsize) then 00704 call DOUG_abort('[Mesh_buildHash] : Node/element hash table'//& 00705 ' overrun', 52) 00706 endif 00707 ! store node appears in element no. i - 00708 ! access hash "fast lookup" table 00709 h = M%hashlook(int(v1/M%hscale)+1) 00710 00711 if (h > M%hashsize) then 00712 write(stream, *) 'hashsize = ', M%hashsize 00713 call flush(stream) 00714 call DOUG_abort('[Mesh_buildHash] : Over run hashtable',52) 00715 endif 00716 00717 990 continue 00718 do while (M%hash(h,1) > 0) 00719 !!$ h = mod(h,hashsize) + 1 00720 h = h + 1 00721 enddo 00722 if (h == M%hashsize) then 00723 h = 1 00724 goto 990 00725 endif 00726 00727 ! found space - store it 00728 M%hash(h,1) = v1 00729 M%hash(h,2) = i 00730 M%hashentries = M%hashentries + 1 00731 enddo 00732 enddo 00733 00734 end subroutine Mesh_buildHash 00735 00736 00737 !----------------------------------------------- 00738 !! Destroy hash table 00739 !----------------------------------------------- 00740 subroutine Mesh_destroyHash(M) 00741 implicit none 00742 00743 type(Mesh), intent(in out) :: M ! Mesh 00744 00745 if (associated(M%hash)) then 00746 deallocate(M%hash) 00747 M%hash => NULL() 00748 end if 00749 if (associated(M%hashlook)) then 00750 deallocate(M%hashlook) 00751 M%hashlook => NULL() 00752 end if 00753 end subroutine Mesh_destroyHash 00754 00755 !----------------------------------------------- 00756 !! Hash function for the edge hash table creation 00757 !----------------------------------------------- 00758 function hashfn(ii, jj, M) result(res) 00759 implicit none 00760 00761 integer, intent(in) :: ii, jj 00762 real(kind=rk), intent(in) :: M 00763 real(kind=rk) :: res 00764 00765 integer :: i, j, k 00766 00767 i = ii-1 00768 j = jj-1 00769 k = int(real(i)/10.0) 00770 00771 res = real(k)*int(M/10.0)+int(real(j-1)/10.0) 00772 end function hashfn 00773 00774 00775 !------------------------------------------------------- 00776 !! Build adjacency structure of the mesh's dual graph 00777 !------------------------------------------------------- 00778 subroutine Mesh_buildDualAdjncy(M, nedges, xadj, adjncy) 00779 use globals, only : stream, D_MSGLVL 00780 implicit none 00781 00782 type(Mesh), intent(in out) :: M ! Mesh 00783 integer, intent(out) :: nedges 00784 integer, dimension(:), pointer :: xadj 00785 integer, dimension(:), pointer :: adjncy 00786 00787 integer :: i, j, k, l, n 00788 real(kind=rk) :: edgeest, edgeband, edgepr 00789 integer :: edgebins, edgebsize, edgehsize, edgec, mxpfree 00790 integer :: h, je0, je, ke, he 00791 integer :: elemcnt, count 00792 logical :: set 00793 integer :: s, s1, n1, n2 00794 integer :: sadjncy 00795 integer, dimension(:,:), pointer :: edgehash 00796 integer, dimension(:), pointer :: edgelook 00797 integer, dimension(:), pointer :: elemuse 00798 00799 !call Mesh_buildEdgeHash(Msh, edgehash) 00800 00801 ! Guess of number of "faces" of the element + an additional bit 00802 ! for those elements with nfrelt<nsd 00803 ! Note that I have no idea what this does in the multi variable case (LS) 00804 00805 if (M%mfrelt /= 2) then 00806 j = 0 00807 do i = 1,M%nell 00808 if (M%nfrelt(i) < M%nsd) j = j+1 00809 enddo 00810 00811 ! CHANGES FOR 3D (LS) 00812 ! edgeest=(nperfree+2*nperstd)*ngf*nell/(2.0*nellg)+j*nperfree 00813 00814 if (M%nsd == 2) then 00815 edgeest = (M%nperfree+2*M%nperstd)*M%ngf*M%nell/(2.0*M%nell)+ & 00816 j*M%nperfree 00817 else 00818 edgeest = (M%nsd-1)*((M%nperfree+2*M%nperstd)*M%ngf*M%nell/ & 00819 (2.0*M%nell)+j*M%nperfree) 00820 endif 00821 00822 ! increase it a bit, make it too "exact" and we do too much searching 00823 edgeest = edgeest*1.3 00824 00825 edgebins = int(edgeest/100.0+0.9) 00826 ! scale factors needed for the hashfn 00827 edgeband = M%nell**(1/M%nsd) 00828 edgepr = hashfn(M%nell,M%nell,edgeband)/(edgebins-0.001) 00829 00830 edgebsize = 100 00831 00832 ! something is very wrong here - Temp Fix. - fix of which bit??? 00833 ! Does not work in 3D - I think!!!! 00834 edgehsize = edgebsize*edgebins-1 00835 00836 else !the case with assembled matrices: 00837 00838 if (M%nsd == 2) then 00839 edgeest = (M%nperfree+2*M%nperstd)*M%nell/1.5 00840 else 00841 edgeest = (M%nsd-1)*((M%nperfree+2*M%nperstd)*M%ngf*M%nell/ & 00842 (1.0*M%nell)) 00843 endif 00844 00845 ! increase it a bit, make it too "exact" and we do too much searching 00846 edgeest = edgeest*1.4 00847 00848 edgebins = int(edgeest/100.0+0.9) 00849 edgeband = M%nell**(1/M%nsd) 00850 edgepr = hashfn(M%nell,M%nell,edgeband)/(edgebins-0.001) 00851 edgebsize = 100 00852 edgehsize = edgebsize*edgebins-1 00853 endif 00854 00855 allocate(edgehash(2,edgehsize)) 00856 edgehash = 0 00857 edgehash(1,edgehsize) = -1 00858 edgehash(2,edgehsize-1) = -1 00859 00860 edgec = 2 00861 00862 allocate(edgelook(edgebins)) 00863 edgelook = 0 00864 00865 ! put in a pseudo random ordering? 00866 do i = 0,edgebins-1 00867 j = i*13 + 1 00868 if (j > edgebins) then 00869 j = mod(j,edgebins)+1 00870 endif 00871 do while (edgelook(j) > 0) 00872 j = j+1 00873 if (j > edgebins) then 00874 j = j-edgebins 00875 endif 00876 enddo 00877 edgelook(j) = i*edgebsize+1 00878 enddo 00879 00880 ! worst case of number of elements using a particular node? 00881 mxpfree = 100 00882 allocate(elemuse(mxpfree)) 00883 00884 00885 ! ** stage 2 ** 00886 ! build hash table of element "edges" 00887 00888 do i = 1,M%ngf 00889 ! find entries in hash table 00890 elemcnt = 0 00891 h = M%hashlook(int(i/M%hscale)+1) 00892 do while (M%hash(h,1) > 0) 00893 if (M%hash(h,1) == i) then 00894 elemcnt = elemcnt+1 00895 ! elements that the freedom 'i' belongs to 00896 elemuse(elemcnt) = M%hash(h,2) 00897 endif 00898 h = h + 1 00899 end do 00900 00901 ! hash add new entries (if any) 00902 do j = 1,elemcnt-1 00903 je0 = elemuse(j) 00904 do k = j+1,elemcnt 00905 ke = elemuse(k) 00906 00907 je = min(je0,ke) 00908 ke = max(je0,ke) 00909 00910 ! only proceed now IF either : 00911 ! (a) nfrelt(je) <= nsd 00912 ! (b) nfrelt(ke) <= nsd 00913 ! (c) >=nsd shared freedoms for je and ke 00914 set = .true. 00915 if (M%mfrelt == 2) goto 1001 00916 if ( (M%nfrelt(ke) >= M%nsd) .and. (M%nfrelt(je)>= M%nsd) ) then 00917 ! count shared freedoms between elements je and ke 00918 count = 0 00919 do l = 1,M%nfrelt(je) 00920 do n = 1,M%nfrelt(ke) 00921 if (M%mhead(l,je) == M%mhead(n,ke)) then 00922 count = count+1 00923 ! as soon as we exceed the condition then jump out 00924 if (count >= M%nsd) goto 1001 00925 endif 00926 enddo 00927 enddo 00928 if (count < M%nsd) set = .false. 00929 endif 00930 00931 1001 continue 00932 if (set) then 00933 he = edgelook(int(hashfn(je,ke,edgeband)/edgepr)+1) 00934 990 continue 00935 do while (edgehash(1,he) > 0) 00936 ! would prefer a "break" construct, goto's are messy 00937 if (edgehash(1,he) == je) then 00938 if (edgehash(2,he) == ke) goto 999 00939 endif 00940 he = he+1 00941 enddo 00942 ! check for reaching the end of the hash table 00943 if (he == edgehsize) then 00944 he = 1 00945 goto 990 00946 endif 00947 00948 edgehash(1,he) = je 00949 edgehash(2,he) = ke 00950 edgec = edgec + 1 00951 if (edgehsize == edgec) then 00952 write(stream, *) 'edgehsize = ', edgehsize 00953 call flush(stream) 00954 call DOUG_abort('[Mesh_buildDualAdjncy] : Edge hash'//& 00955 ' table has overrun', 52) 00956 endif 00957 999 continue 00958 endif 00959 enddo 00960 enddo 00961 enddo 00962 00963 if (D_MSGLVL > 1) then 00964 write(stream, *) 00965 write(stream, *) 'Statistics on edge hash table:' 00966 write(stream, '(a,i8)') ' Size of edge hash table : ', edgehsize 00967 write(stream, 5600) edgec-2,100*(edgec-2)/edgehsize 00968 call flush(stream) 00969 5600 format(' Number of edges found : ',i8,' (',i3,'% of hash table)') 00970 end if 00971 00972 00973 ! allocation for the adjacency data 00974 allocate(xadj(M%nell+1)) 00975 xadj = 0 00976 00977 ! scan edge list and create the necessary adjacency structure 00978 ! for the METIS call 00979 00980 ! pass 1 00981 do i = 1,edgehsize 00982 if (edgehash(1,i) > 0) then 00983 xadj(edgehash(1,i)) = xadj(edgehash(1,i)) + 1 00984 xadj(edgehash(2,i)) = xadj(edgehash(2,i)) + 1 00985 endif 00986 enddo 00987 00988 ! Allocate partition array 'eptnmap'. 00989 ! Will store result of partitioning made by 'Graph_Partition()' 00990 if (.not.associated(M%eptnmap)) allocate(M%eptnmap(M%nell)) 00991 M%eptnmap = 0 00992 00993 s = xadj(1) 00994 xadj(1) = 1 00995 M%eptnmap(1) = 1 00996 do i = 2,M%nell 00997 s1 = xadj(i) 00998 xadj(i) = xadj(i-1)+s 00999 M%eptnmap(i) = xadj(i) 01000 s = s1 01001 enddo 01002 xadj(M%nell+1) = xadj(M%nell) + s 01003 01004 ! size of 'adjncy' array - 2*(number of graph edges) 01005 sadjncy = xadj(M%nell+1) - 1 01006 01007 ! allocate array for the graph adjacency 01008 allocate(adjncy(sadjncy)) 01009 01010 ! pass 2 of the data 01011 do i = 1,edgehsize-1 01012 if (edgehash(1,i) /= 0) then 01013 n1 = edgehash(1,i) 01014 n2 = edgehash(2,i) 01015 adjncy(M%eptnmap(n1)) = n2 01016 M%eptnmap(n1) = M%eptnmap(n1)+1 01017 adjncy(M%eptnmap(n2)) = n1 01018 M%eptnmap(n2) = M%eptnmap(n2)+1 01019 endif 01020 enddo 01021 01022 ! number of edges of the graph * 2 01023 nedges = sadjncy/2 01024 01025 01026 ! Deallocate temporary arrays 01027 deallocate( & 01028 edgehash, & 01029 edgelook, & 01030 elemuse) 01031 01032 end subroutine Mesh_buildDualAdjncy 01033 !============================================================= 01034 ! 01035 !! Partitioning 01036 ! 01037 !------------------------------------------------------------- 01038 !! Partition mesh's dual graph 01039 !------------------------------------------------------------- 01040 subroutine Mesh_partitionDual(M, nparts, method, part_options) 01041 01042 use globals, only : stream, D_MSGLVL 01043 01044 implicit none 01045 01046 type(Mesh) :: M 01047 integer, intent(in) :: nparts, method 01048 integer, dimension(6), intent(in) :: part_options 01049 01050 integer :: i 01051 01052 call Graph_Partition(M%G, nparts, method, part_options) 01053 01054 if (D_MSGLVL > 0 ) & 01055 write(stream, FMT='(a,i5)') ' Number of edgecuts : ', M%G%edgecut 01056 if (D_DEBUGLVL > 4 ) then 01057 do i = 1,size(M%G%part,1) 01058 write(stream, FMT='(a,i5,a,i3)') 'M%G%part(', i, ') = ', M%G%part(i) 01059 end do 01060 end if 01061 01062 ! Save result in Mesh object 01063 M%parted = M%G%parted 01064 M%nparts = M%G%nparts 01065 M%eptnmap = M%G%part 01066 01067 end subroutine Mesh_partitionDual 01068 01069 01070 !------------------------------------------- 01071 !! Calculate number of elements in partitions 01072 !------------------------------------------- 01073 subroutine Mesh_calcElemsInParts(M) 01074 01075 type(Mesh), intent(in out) :: M 01076 integer :: p, el 01077 01078 allocate(M%partnelems(M%nparts)) 01079 M%partnelems = 0 01080 do el = 1,M%nell 01081 do p = 1,M%nparts 01082 if (M%eptnmap(el) == p) then 01083 M%partnelems(p) = M%partnelems(p) + 1 01084 end if 01085 end do 01086 end do 01087 end subroutine Mesh_calcElemsInParts 01088 01089 01090 !------------------------------------------------ 01091 !! Build maps: 01092 !! - two maps to map global degrees of freedoms to 01093 !! process's local ones and back 01094 !! - send map for freedoms 01095 !! Finds neighbours to given partition: 01096 !! - neighbours (calls local Mesh_findNghbrs) 01097 !------------------------------------------------ 01098 subroutine Mesh_MapsAndNghbrs(M) 01099 use globals, only: myrank 01100 implicit none 01101 01102 type(Mesh), intent(in out) :: M 01103 01104 integer :: gf, j, ptn, h, lf, interf_end 01105 integer :: nlf ! # of local freedoms 01106 integer :: nlf_interf ! interf freedoms counter 01107 integer :: nlf_inner ! inner freedoms counter 01108 integer :: ptncnt ! partition counter 01109 logical :: set 01110 01111 integer, dimension(M%nparts) :: ptnuse ! auxiliary: to which partition 01112 logical :: is_inner,is_mine 01113 ! the given freedom belogs to 01114 integer,dimension(:),pointer :: lg_fmap ! working-space 01115 01116 01117 if (M%hashsize <= 0) & 01118 call DOUG_abort('[Mesh_MapsAndNghbrs] : hash table was not built!',-1) 01119 if ((.not.associated(M%eptnmap)).or.(any(M%eptnmap == 0))) & 01120 call DOUG_abort('[Mesh_MapsAndNghbrs] : partitions map have to be'//& 01121 ' allocated and filled in!',-1) 01122 01123 allocate(lg_fmap(M%ngf)) 01124 lg_fmap = 0 01125 01126 allocate(M%nfreesend_map(M%nparts)) 01127 M%nfreesend_map = 0 01128 01129 nlf=0 01130 interf_end=0 01131 do gf = 1,M%ngf ! loop over global freedoms 01132 is_mine=.false. 01133 is_inner=.true. 01134 ptncnt = 0 01135 h = M%hashlook(int(gf/M%hscale)+1) 01136 do while (M%hash(h,1) > 0) 01137 if (M%hash(h,1) == gf) then 01138 ! Handle partition access to it 01139 set = .true. 01140 ptn = M%eptnmap(M%hash(h,2)) 01141 if (ptn == (myrank+1)) then 01142 is_mine=.true. 01143 else ! element 'gf'-node belongs to is NOT mine 01144 is_inner=.false. 01145 do j = 1,ptncnt 01146 if (ptnuse(j) == ptn) then 01147 set = .false. 01148 endif 01149 enddo 01150 if (set) then ! We are keeping the list of partitions the node 01151 ptncnt = ptncnt + 1 ! belongs to 01152 ptnuse(ptncnt) = ptn 01153 endif 01154 endif 01155 endif 01156 h = h + 1 01157 enddo ! do while 01158 if (is_mine) then 01159 nlf=nlf+1 01160 if (is_inner) then 01161 ! Book freedom 01162 lg_fmap(nlf) = gf 01163 else ! is interface freedom: (assigned -) 01164 interf_end=interf_end+1 01165 lg_fmap(nlf) = -gf 01166 ! preserve the list of partitions the node additionally belongs to: 01167 do j = 1,ptncnt 01168 M%nfreesend_map(ptnuse(j)) = M%nfreesend_map(ptnuse(j)) + 1 01169 enddo 01170 endif 01171 endif 01172 end do ! gf 01173 01174 allocate(M%gl_fmap(M%ngf),M%lg_fmap(nlf)) 01175 M%gl_fmap=0 01176 M%nlf=nlf 01177 allocate(M%inner_interf_fmask(M%nlf)) 01178 nlf_interf = 0 01179 nlf_inner = interf_end 01180 do lf=1,nlf 01181 gf=lg_fmap(lf) 01182 if (gf<0) then ! interface freedom 01183 nlf_interf = nlf_interf + 1 01184 M%lg_fmap(nlf_interf)=-gf 01185 M%gl_fmap(-gf)=nlf_interf 01186 else ! inner freedom 01187 nlf_inner = nlf_inner + 1 01188 M%lg_fmap(nlf_inner)=gf 01189 M%gl_fmap(gf)=nlf_inner 01190 endif 01191 enddo 01192 M%inner_interf_fmask(1:interf_end) = D_FREEDOM_INTERF 01193 M%inner_interf_fmask(interf_end+1:nlf) = D_FREEDOM_INNER 01194 01195 deallocate(lg_fmap) 01196 01197 if (D_MSGLVL > 4) then 01198 write(stream,'(/a,i6,a)') 'Global to local mapping: [M%gl_fmap]'//& 01199 ' :size [',M%ngf,']:' 01200 do j = 1,M%ngf 01201 !write(stream, '(a,i6,a,e22.15)') ' [',j,']=',M%gl_fmap(j) 01202 write(stream, '(a,i6,a,i16)') ' [',j,']=',M%gl_fmap(j) 01203 end do 01204 call flush(stream) 01205 write(stream,'(/a,i6,a)') 'Local to global mapping: [M%lg_fmap]'//& 01206 ' :size [',M%nlf,']:' 01207 do j = 1,M%nlf 01208 !write(stream, '(a,i6,a,e22.15)') ' [',j,']=',M%lg_fmap(j) 01209 write(stream, '(a,i6,a,i16)') ' [',j,']=',M%lg_fmap(j) 01210 end do 01211 call flush(stream) 01212 end if 01213 ! Now we are able very simply to calculate the number of 01214 ! our neighbours and their ranks 01215 call Mesh_findNghbrs(M) 01216 01217 end subroutine Mesh_MapsAndNghbrs 01218 01219 01220 !---------------------------- 01221 !! Finds my neighbours : ranks 01222 !---------------------------- 01223 subroutine Mesh_findNghbrs(M) 01224 01225 implicit none 01226 01227 type(Mesh), intent(in out) :: M 01228 01229 integer :: p, i 01230 01231 if (.not.associated(M%nfreesend_map)) & 01232 call DOUG_abort('[Mesh_findNghbrs] : "Number of freedoms to send"'//& 01233 ' map have to be built first!',-1) 01234 01235 M%nnghbrs = 0 01236 do p = 1,M%nparts 01237 if (M%nfreesend_map(p) /= 0) & 01238 M%nnghbrs = M%nnghbrs + 1 01239 end do 01240 01241 allocate(M%nghbrs(M%nnghbrs)) 01242 M%nghbrs = 0 01243 i = 0 01244 do p = 1,M%nparts 01245 if (M%nfreesend_map(p) /= 0) then 01246 i = i + 1 01247 M%nghbrs(i) = p - 1 01248 end if 01249 end do 01250 if (i /= M%nnghbrs) & 01251 call DOUG_abort('[Mesh_findNghbrs] : SEVERE : j /= M%nnghbrs ',-1) 01252 01253 write(stream,'(/a,i3,a)',advance='no') 'RANK<',myrank,'> my neighbours [' 01254 do i = 1,M%nnghbrs 01255 write(stream,'(i3)',advance='no') M%nghbrs(i) 01256 if (i /= M%nnghbrs) write(stream,'(a)',advance='no') ', ' 01257 end do 01258 write(stream,'(a)') ']' 01259 call flush(stream) 01260 01261 end subroutine Mesh_findNghbrs 01262 01263 01264 !------------------------------ 01265 !! Find number of local freedoms 01266 !! (local to sub-partition) 01267 !------------------------------ 01268 subroutine Mesh_findNLF(M) 01269 use globals, only: myrank 01270 implicit none 01271 01272 type(Mesh), intent(in out) :: M 01273 01274 integer :: gf, h 01275 logical :: f_booked 01276 01277 M%nlf = 0 01278 01279 if (associated(M%gl_fmap)) then ! we can do it very fast 01280 ! if we have 'Mesh%gl_fmap' built 01281 M%nlf = maxval(M%gl_fmap) 01282 else ! otherwise count them explicitly 01283 f_booked = .false. 01284 do gf = 1,M%ngf ! walk through global freedoms 01285 h = M%hashlook(int(gf/M%hscale)+1) 01286 do while (M%hash(h,1) > 0) 01287 if (M%hash(h,1) == gf) then ! found this freedom in hash table 01288 ! the element this freedom belongs to is in my partition 01289 if (M%eptnmap(M%hash(h,2)) == (myrank+1)) then 01290 if (.not.f_booked) then 01291 M%nlf = M%nlf + 1 01292 f_booked = .true. 01293 end if 01294 end if 01295 end if 01296 h = h + 1 01297 end do ! do while 01298 f_booked = .false. 01299 end do 01300 end if 01301 01302 end subroutine Mesh_findNLF 01303 01304 !----------------------------------------------------------- 01305 !! Builds: 01306 !! - global to local, local to global freedoms maps; 01307 !! - freedoms/elements send maps for local freedoms/elements; 01308 !! - inner/interface mask for local freedoms 01309 !! Allocates and fills in: 01310 !! - gl_fmap, lg_fmap; 01311 !! - nfreesend_map, nelemsend_map; 01312 !! - inner_interf_fmask 01313 !----------------------------------------------------------- 01314 subroutine Mesh_buldMapsNghbrsMasksNLF(M) 01315 implicit none 01316 01317 type(Mesh), intent(in out) :: M 01318 01319 ! global to local, local to global maps, 01320 ! neighbours,interface-inner decisions, 01321 ! changed it into ALL-IN-ONE subroutine: 01322 call Mesh_MapsAndNghbrs(M) 01323 01324 end subroutine Mesh_buldMapsNghbrsMasksNLF 01325 !================================ 01326 ! 01327 !! MPI wrappers 01328 ! 01329 !----------------------------- 01330 !! Master Send Mesh's info to slaves 01331 !----------------------------- 01332 subroutine Mesh_paramsMPIISEND(M) 01333 implicit none 01334 01335 type(Mesh), intent(in) :: M 01336 01337 integer, parameter :: bsize = 5 01338 integer, dimension(bsize) :: buf 01339 integer :: request 01340 integer :: i, ierr 01341 01342 if (.not.ismaster()) & 01343 call DOUG_abort('[Mesh_paramsMPIISEND] : Only master is allowed'//& 01344 ' to use me.',-1) 01345 01346 buf(1) = M%nell 01347 buf(2) = M%ngf 01348 buf(3) = M%mfrelt 01349 buf(4) = M%nsd 01350 buf(5) = M%nnode 01351 01352 do i = 1,numprocs 01353 call MPI_ISEND(buf, bsize, MPI_INTEGER, & 01354 i, D_TAG_MESH_INFO, MPI_COMM_WORLD, request, ierr) 01355 end do 01356 end subroutine Mesh_paramsMPIISEND 01357 01358 01359 !--------------------------------- 01360 !! Receive info for the Mesh object 01361 !--------------------------------- 01362 subroutine Mesh_paramsMPIRECV(M) 01363 implicit none 01364 01365 type(Mesh), intent(in out) :: M 01366 01367 integer, parameter :: bsize = 5 01368 integer, dimension(bsize) :: buf 01369 integer :: status(MPI_STATUS_SIZE) 01370 integer :: ierr 01371 01372 if (.not.isslave()) & 01373 call DOUG_abort('[Mesh_paramsMPIRECV] : Only slaves are allowed'//& 01374 ' to use me.',-1) 01375 01376 call MPI_RECV(buf, bsize, MPI_INTEGER, & 01377 D_MASTER, D_TAG_MESH_INFO, MPI_COMM_WORLD, status, ierr) 01378 01379 M%nell = buf(1) 01380 M%ngf = buf(2) 01381 M%mfrelt = buf(3) 01382 M%nsd = buf(4) 01383 M%nnode = buf(5) 01384 01385 end subroutine Mesh_paramsMPIRECV 01386 01387 01388 !----------------------------------- 01389 !! Master broadcasts Mesh's info data 01390 !----------------------------------- 01391 subroutine Mesh_paramsMPIBCAST(M) 01392 implicit none 01393 type(Mesh), intent(in out) :: M 01394 01395 integer, parameter :: bsize = 5 01396 integer, dimension(bsize) :: buf 01397 integer :: ierr 01398 01399 if (ismaster()) then 01400 buf(1) = M%nell 01401 buf(2) = M%ngf 01402 buf(3) = M%mfrelt 01403 buf(4) = M%nsd 01404 buf(5) = M%nnode 01405 end if 01406 01407 call MPI_BCAST(buf, bsize, MPI_INTEGER, & 01408 D_MASTER, MPI_COMM_WORLD, ierr) 01409 01410 if (isslave()) then 01411 M%nell = buf(1) 01412 M%ngf = buf(2) 01413 M%mfrelt = buf(3) 01414 M%nsd = buf(4) 01415 M%nnode = buf(5) 01416 end if 01417 01418 end subroutine Mesh_paramsMPIBCAST 01419 01420 01421 !------------------------------------------------------------------ 01422 !! Master sends out Mesh's data components to the others: 01423 !! nfrelt, mhead, freemap, coords, eptnmap, hash 01424 !------------------------------------------------------------------ 01425 subroutine Mesh_dataMPIISENDRECV(M, & 01426 nfrelt, & 01427 mhead, & 01428 freemap,& 01429 coords, & 01430 eptnmap, & 01431 hash, & 01432 freemask) 01433 implicit none 01434 01435 type(Mesh), intent(in out) :: M 01436 logical, intent(in), optional :: nfrelt, mhead, freemap 01437 logical, intent(in), optional :: coords, eptnmap, hash, freemask 01438 integer :: request, status(MPI_STATUS_SIZE) 01439 integer :: p, ierr 01440 integer, parameter :: D_TAG_MESH_NFRELT = 201 01441 integer, parameter :: D_TAG_MESH_MHEAD = 202 01442 integer, parameter :: D_TAG_MESH_FREEMAP = 203 01443 integer, parameter :: D_TAG_MESH_COORDS = 204 01444 integer, parameter :: D_TAG_MESH_EPTNMAP = 205 01445 integer, parameter :: D_TAG_MESH_HASH = 206 01446 integer, parameter :: D_TAG_MESH_FREEMASK = 207 01447 01448 ! Send/recv 'nfrelt' 01449 if (present(nfrelt)) then 01450 if (ismaster()) then 01451 do p = 1,numprocs-1 01452 call MPI_ISEND(M%nfrelt, M%nell, MPI_INTEGER, & 01453 p, D_TAG_MESH_NFRELT, MPI_COMM_WORLD, request, ierr) 01454 end do 01455 else 01456 if (.NOT.associated(M%nfrelt)) allocate(M%nfrelt(M%nell)) 01457 call MPI_RECV(M%nfrelt, M%nell, MPI_INTEGER, & 01458 D_MASTER, D_TAG_MESH_NFRELT, MPI_COMM_WORLD, status, ierr) 01459 end if 01460 end if 01461 01462 ! Send/recv 'mhead' 01463 if (present(mhead)) then 01464 if (ismaster()) then 01465 do p = 1,numprocs-1 01466 call MPI_ISEND(M%mhead, M%mfrelt*M%nell, MPI_INTEGER, & 01467 p, D_TAG_MESH_MHEAD, MPI_COMM_WORLD, request, ierr) 01468 end do 01469 else 01470 if (.NOT.associated(M%mhead)) allocate(M%mhead(M%mfrelt,M%nell)) 01471 call MPI_RECV(M%mhead, M%mfrelt*M%nell, MPI_INTEGER, & 01472 D_MASTER, D_TAG_MESH_MHEAD, MPI_COMM_WORLD, status, ierr) 01473 end if 01474 end if 01475 01476 ! Send/recv 'freemap' 01477 if (present(freemap)) then 01478 if (ismaster()) then 01479 do p = 1,numprocs-1 01480 call MPI_ISEND(M%freemap, M%ngf, MPI_INTEGER, & 01481 p, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, request, ierr) 01482 end do 01483 else 01484 if (.not.associated(M%freemap)) allocate(M%freemap(M%ngf)) 01485 call MPI_RECV(M%freemap, M%ngf, MPI_INTEGER, & 01486 D_MASTER, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, status, ierr) 01487 end if 01488 end if 01489 01490 ! Send/recv 'coords' 01491 if (present(coords)) then 01492 if (ismaster()) then 01493 do p = 1,numprocs-1 01494 call MPI_ISEND(M%coords, M%nnode*M%nsd, MPI_xyzkind, & 01495 p, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, request, ierr) 01496 end do 01497 else 01498 if (.not.associated(M%coords)) allocate(M%coords(M%nsd,M%nnode)) 01499 call MPI_RECV(M%coords, M%nnode*M%nsd, MPI_xyzkind, & 01500 D_MASTER, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, status, ierr) 01501 end if 01502 end if 01503 01504 ! Send/recv 'eptnmap' 01505 if (present(eptnmap)) then 01506 call MPI_BCAST(M%nparts, 1, MPI_INTEGER, & 01507 D_MASTER, MPI_COMM_WORLD, ierr) 01508 01509 if (ismaster()) then 01510 do p = 1,numprocs-1 01511 call MPI_ISEND(M%eptnmap, M%nell, MPI_INTEGER, & 01512 p, D_TAG_MESH_EPTNMAP, MPI_COMM_WORLD, request, ierr) 01513 end do 01514 else 01515 if (.not.associated(M%eptnmap)) allocate(M%eptnmap(M%nell)) 01516 call MPI_RECV(M%eptnmap, M%nell, MPI_INTEGER, & 01517 D_MASTER, D_TAG_MESH_EPTNMAP, MPI_COMM_WORLD, status, ierr) 01518 M%parted = .true. 01519 end if 01520 end if 01521 01522 ! Send/recv 'hash' 01523 if (present(hash)) then 01524 01525 ! First, broadcast size of a hash table to be sent later on 01526 call MPI_BCAST(M%hashsize, 1, MPI_INTEGER, & 01527 D_MASTER, MPI_COMM_WORLD, ierr) 01528 01529 ! Slaves allocate hash table 01530 if (isslave()) allocate(M%hash(M%hashsize,2)) 01531 01532 ! Send/recv hash table 01533 if (ismaster()) then 01534 do p = 1,numprocs-1 01535 call MPI_ISEND(M%hash, 2*M%hashsize, MPI_INTEGER, & 01536 p, D_TAG_MESH_HASH, MPI_COMM_WORLD, request, ierr) 01537 end do 01538 else 01539 call MPI_RECV(M%hash, 2*M%hashsize, MPI_INTEGER, & 01540 D_MASTER, D_TAG_MESH_HASH, MPI_COMM_WORLD, status, ierr) 01541 end if 01542 end if 01543 01544 01545 ! Send/recv 'freemask' 01546 if (present(freemask)) then 01547 ! Send/recv freedoms mask 01548 if (ismaster()) then 01549 do p = 1,numprocs-1 01550 call MPI_ISEND(M%freemask, M%ngf, MPI_BYTE, & 01551 p, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, request, ierr) 01552 end do 01553 else 01554 if (.not.associated(M%freemask)) allocate(M%freemask(M%ngf)) 01555 call MPI_RECV(M%freemask, M%ngf, MPI_BYTE, & 01556 D_MASTER, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, status, ierr) 01557 end if 01558 end if 01559 end subroutine Mesh_dataMPIISENDRECV 01560 !================================ 01561 ! 01562 !! Printing routines 01563 ! 01564 !-------------------------------- 01565 !! Prints out some data about mesh 01566 !-------------------------------- 01567 subroutine Mesh_printInfo(M) 01568 01569 use globals, only : stream 01570 01571 implicit none 01572 01573 type(Mesh), intent(in) :: M 01574 01575 write(stream, *) 01576 write(stream, *) 'Mesh object [info]:' 01577 write(stream, FMT='(A,I7)') ' # of elements (nell) = ',& 01578 M%nell 01579 write(stream, FMT='(A,I7)') ' # of global freedoms (ngf) = ',& 01580 M%ngf 01581 write(stream, FMT='(A,I7)') ' # of local freedoms (nlf) = ',& 01582 M%nlf 01583 write(stream, FMT='(A,I7)') ' # of spacial dimensions (nsd) = ',& 01584 M%nsd 01585 write(stream, FMT='(A,I7)') ' max # of freedoms per elem (mfrelt) = ',& 01586 M%mfrelt 01587 write(stream, FMT='(A,I7)') ' total # of nodes in a mesh (nnode) = ',& 01588 M%nnode 01589 write(stream, *) 01590 01591 end subroutine Mesh_printInfo 01592 01593 01594 !---------------------------------------------- 01595 !! Prints freedoms belonging to elements (mhead) 01596 !---------------------------------------------- 01597 subroutine Mesh_printElemFree(M) 01598 01599 use globals, only : stream 01600 01601 implicit none 01602 01603 type(Mesh), intent(in) :: M 01604 integer :: i, j, k, main 01605 01606 write(stream, *) 'Mesh object [element freedoms]:' 01607 01608 if (.not.associated(M%mhead)) then 01609 write(stream, *) ' object not initialised.' 01610 return 01611 end if 01612 01613 ! main part of 'mhead' 01614 main = M%nell - mod(M%nell,10) 01615 do k = 1,main,10 01616 do j = k,k+9 01617 if (j == k+9) then 01618 write(stream, FMT='(a,i5,a)', advance='no') '|',j,'|' 01619 else 01620 write(stream, FMT='(a,i5)', advance='no') '|',j 01621 end if 01622 end do 01623 write(stream, *) 01624 do i = 1,M%mfrelt 01625 write(stream, FMT='(10i6)') M%mhead(i,k:k+9) 01626 end do 01627 write(stream, *) 01628 end do 01629 01630 ! reminder of 'mhead' 01631 do j = main+1,M%nell 01632 if (j == M%nell) then 01633 write(stream, FMT='(a,i5,a)', advance='no') '|',j,'|' 01634 else 01635 write(stream, FMT='(a,i5)', advance='no') '|',j 01636 end if 01637 end do 01638 write(stream, *) 01639 do i = 1,M%mfrelt 01640 write(stream, FMT='(10i6)') M%mhead(i,main+1:M%nell) 01641 end do 01642 01643 end subroutine Mesh_printElemFree 01644 01645 !================================== 01646 ! 01647 !! I/O subroutines 01648 ! 01649 !---------------------------------- 01650 !! Print out integer vector 01651 !---------------------------------- 01652 subroutine Mesh_IVect_Print(x, str) 01653 implicit none 01654 integer, dimension(:), intent(in) :: x 01655 character*(*), optional, intent(in) :: str 01656 01657 integer :: i, n 01658 01659 n = size(x) 01660 if (present(str)) then 01661 write(stream,'(/a,i6,a)') str//' :size [',n,']:' 01662 else 01663 write(stream,'(/a,i6,a)') 'vector :size [',n,']:' 01664 end if 01665 do i = 1,n 01666 write(stream, '(a,i6,a,i10)') ' [',i,']=',x(i) 01667 end do 01668 call flush(stream) 01669 end subroutine Mesh_IVect_Print 01670 01671 01672 01673 end module Mesh_class
1.7.3-20110217