DOUG 0.2

Mesh.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 
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