DOUG 0.2

Mesh_plot.F90

Go to the documentation of this file.
00001 module Mesh_plot_mod
00002   use Mesh_class
00003   use Aggregate_mod
00004 
00005   implicit none
00006 
00007 
00008 #include<doug_config.h>
00009 
00010 contains
00011   !================================================
00012   !
00013   !! Plotting routines
00014   !
00015 
00017   subroutine Mesh_pl2D_mesh(Msh)
00018     type(Mesh), intent(inout) :: Msh
00019 
00020     call Mesh_pl2D_pointCloud(Msh,D_PLPLOT_INIT)
00021     ! Plots mesh's dual graph
00022     call Mesh_pl2D_plotGraphDual(Msh,D_PLPLOT_END)
00023     ! Mesh & its Dual Graph
00024     call Mesh_pl2D_plotMesh(Msh, D_PLPLOT_INIT)
00025     call Mesh_pl2D_plotGraphDual(Msh, D_PLPLOT_END)
00026   end subroutine Mesh_pl2D_mesh
00027 
00029   subroutine Mesh_pl2D_partitions(Msh)
00030     type(Mesh), intent(inout) :: Msh
00031 
00032     ! Draw colored partitoined graph
00033     call Mesh_pl2D_plotGraphParted(Msh)
00034 
00035     ! Plot partitions of the mesh
00036     ! NB: Check for multivariable case! TODO
00037     call Mesh_pl2D_Partition(Msh)
00038     ! Partition with Dual Graph upon it
00039     call Mesh_pl2D_Partition(Msh, D_PLPLOT_INIT)
00040     call Mesh_pl2D_plotGraphDual(Msh, D_PLPLOT_CONT)
00041     call Mesh_pl2D_pointCloud(Msh,D_PLPLOT_END)
00042   end subroutine Mesh_pl2D_partitions
00043 
00044   !------------------------------------------------
00045   !! Plot cloud field
00046   !------------------------------------------------
00047   subroutine Mesh_pl2D_pointCloud(M, INIT_CONT_END)
00048 
00049     use globals, only : stream
00050 
00051     implicit none
00052 
00053     type(Mesh), intent(in out)           :: M
00054     integer,    intent(in),    optional  :: INIT_CONT_END
00055 
00056     real(kind=xyzk), dimension(:), pointer :: xc, yc
00057     real(kind=xyzk)                        :: xmin, xmax, ymin, ymax
00058     integer                              :: n, i, j
00059     character*2                          :: buf2
00060     character*5                          :: buf5
00061 
00062 #ifdef D_WANT_PLPLOT_YES
00063 
00064     write(stream, *)
00065     write(stream, *) '[Mesh_pl2D_pointCloud] : Plotting 2D mesh cloud'//&
00066          ' of points.'
00067     if (M%nsd /= 2) then
00068        write(stream, *) '   Only for 2D meshes. Skipping...'
00069        return
00070     end if
00071 
00072     n = size(M%coords)
00073 
00074     ! We are not going to change array's 'coords' values, so just point on them
00075     ! Do it merely out of convenience for having short names for variables
00076     xc => M%coords(1,:)
00077     yc => M%coords(2,:)
00078 
00079     if (.not.present(INIT_CONT_END).or.&
00080          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00081 
00082        xmin = minval(xc)
00083        xmax = maxval(xc)
00084        ymin = minval(yc)
00085        ymax = maxval(yc)
00086 
00087        xmin = xmin - xmax/10.0
00088        xmax = xmax + xmax/10.0
00089        ymin = ymin - ymax/10.0
00090        ymax = ymax + ymax/10.0
00091 
00092        call plsdev("xwin")
00093        call plinit()
00094 
00095        call plenv (xmin, xmax, ymin, ymax, 0, 0);
00096 
00097        call plcol0(1) ! red
00098        write(buf5, '(i5)') M%nnode
00099        call pllab( '(x)', '(y)', 'Mesh : cloud of points ['//buf5//']' )
00100     end if
00101 
00102     call plcol0(15) ! white
00103     call plssym(0.0d0, 2.0d0)
00104     call plpoin(n/2, xc, yc, 1)
00105 
00106     if (M%nnode <= 100) then
00107        call plcol0(6) ! wheat
00108        call plssym(0.0d0, 5.0d0)
00109        do i = 1,n/2
00110           write(buf2, '(i2)') i
00111           call plptex(xc(i), yc(i), 0.0, 0.0, 0, buf2)
00112        end do
00113     end if
00114 
00115     if (.not.present(INIT_CONT_END).or.&
00116          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_END))) then
00117        call plend()
00118     end if
00119 
00120     nullify(xc, yc)
00121 
00122 #else
00123     write(stream, '(/a)') ' [Mesh_pl2D_pointCloud] : Compiled w/o plotting'//&
00124          ' support!'
00125 #endif
00126 
00127   end subroutine Mesh_pl2D_pointCloud
00128 
00129 
00130   !----------------------------------------------
00131   !! Plot mesh
00132   !----------------------------------------------
00133   subroutine Mesh_pl2D_plotMesh(M, INIT_CONT_END)
00134 
00135     use globals, only : stream
00136 
00137     implicit none
00138 
00139     type(Mesh), intent(in)            :: M
00140     integer,    intent(in), optional  :: INIT_CONT_END
00141 
00142     real(kind=xyzk), dimension(:), pointer :: xc, yc
00143     real(kind=xyzk)                        :: xmin, xmax, ymin, ymax
00144 
00145     type(Polygon)                       :: ep ! Coordinates of element/polygon
00146     integer                             :: npol ! Number of vertices in polygon
00147     integer                             :: n, e, i
00148     character*6                         :: buf61, buf62, buf63
00149     ! indexes to get element node coordinates from 'xc', 'yc'
00150     integer,   dimension(:), allocatable :: ind_coords
00151 
00152 #ifdef D_WANT_PLPLOT_YES
00153 
00154     write(stream, *)
00155     write(stream, *) '[Mesh_pl2D_plotMesh] : Plotting 2D mesh.'
00156     if (M%nsd /= 2) then
00157        write(stream, *) '   Only for 2D meshes. Skipping...'
00158        return
00159     end if
00160 
00161     n = size(M%coords)
00162 
00163     xc => M%coords(1,:)
00164     yc => M%coords(2,:)
00165 
00166     if (.not.present(INIT_CONT_END).or.&
00167          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00168 
00169        xmin = minval(xc)
00170        xmax = maxval(xc)
00171        ymin = minval(yc)
00172        ymax = maxval(yc)
00173 
00174        xmin = xmin - xmax/10.0
00175        xmax = xmax + xmax/10.0
00176        ymin = ymin - ymax/10.0
00177        ymax = ymax + ymax/10.0
00178 
00179        call plsdev("xwin")
00180        call plinit()
00181 
00182        call plenv (xmin, xmax, ymin, ymax, 0, 0);
00183 
00184        call plcol0(1) ! red
00185 
00186        write(buf61, '(i6)') M%nell
00187        write(buf62, '(i6)') M%ngf
00188        write(buf63, '(i6)') M%nnode
00189        call pllab( '(x)', '(y)', &
00190             'Mesh : nell='//buf61//'; ngf='//buf62//'; nnode='//buf63)
00191     end if
00192 
00193     allocate(ind_coords(M%mfrelt))
00194     ind_coords = 0
00195 
00196 
00197     ! Another algorithm (must be faster):
00198     ! 1. plot cloud of white points - all nodes in object
00199     ! 2. plot cloud of yellow points - only nodes from 'M%freemap' array
00200     ! 3. plot elements - yellow points connected by green lines
00201 
00202     ! Plot all nodes in white
00203     call plcol0(15) ! white
00204     call plssym(0.0d0, 2.0d0)
00205     call plpoin(n/2, xc, yc, 1)
00206 
00207     do e = 1,M%nell
00208 
00209        if (M%nfrelt(e) > M%nsd) then
00210 
00211           ! Get global node numbering to fetch nodes'
00212           ! global coordinate values later.
00213           npol = 0
00214           do i = 1,M%nfrelt(e)
00215              if (M%mhead(i,e) <= M%ngf) then
00216                 npol = npol + 1
00217                 ind_coords(npol) = M%freemap(M%mhead(i,e))
00218              end if
00219           end do
00220 
00221           if (npol == 3) then
00222              ! NB: Just to speed up plotting
00223              !
00224              ! Trivial 3-node (triangle) element case
00225              !
00226              ! There can be more than 3-nodes elements on the boundary
00227              ! (appeared due to Dirichlet BC), having only 3 nodes left
00228              ! on one line. So, it will fail in multi-variable case.
00229              ! Check it too! TODO
00230 
00231              ! Plot polygon
00232              call plcol0(7) ! grey
00233              ! "Close" 3-node elements to form polygons
00234              call plline(4, &
00235                   (/xc(ind_coords(1:3)),xc(ind_coords(1))/), &
00236                   (/yc(ind_coords(1:3)),yc(ind_coords(1))/) )
00237 
00238 
00239              ! Plot vertexes
00240              call plcol0(1) ! red
00241              call plssym(0.0d0, 2.0d0)
00242              call plpoin(M%nfrelt(e), &
00243                   xc(ind_coords(1:3)), yc(ind_coords(1:3)), 1)
00244 
00245           else if (npol > 3) then ! More than 3 is a little bit tricky
00246              ep = Polygon_New(npol)
00247 
00248              call Polygon_Init(ep, &
00249                   xc(ind_coords(1:npol)), yc(ind_coords(1:npol)) )
00250              call Polygon_sortVerts(ep)
00251 
00252              call Polygon_pl2D_Plot(ep, D_PLPLOT_CONT)
00253 
00254              call Polygon_Destroy(ep)
00255           end if
00256 
00257        end if
00258     end do
00259 
00260     if (.not.present(INIT_CONT_END).or.&
00261          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_END))) then
00262        call plend()
00263     end if
00264 
00265     nullify(xc, yc)
00266     deallocate(ind_coords)
00267 
00268 #else
00269     write(stream, '(/a)') ' [Mesh_pl2D_plotMesh] : Compiled w/o plotting'//&
00270          ' support!'
00271 #endif
00272 
00273   end subroutine Mesh_pl2D_plotMesh
00274 
00275 
00276 
00277   !----------------------------------------------
00278   !! Plot aggregates
00279   !----------------------------------------------
00280   subroutine Mesh_pl2D_plotAggregate(aggr,M,rowstart,colnrs,filename, &
00281     caggrnum,INIT_CONT_END)
00282     use globals !, only : stream
00283     implicit none
00284     type(Aggrs), intent(in)            :: aggr
00285     type(Mesh), intent(in)             :: M
00286     integer, dimension(:), pointer     :: aggrnum ! aggregate # for each node
00287     integer, dimension(:), pointer     :: rowstart,colnrs
00288     character*(*),intent(in)           :: filename
00289     integer,    pointer                :: neighood,nagrs,nisolated
00290     integer,dimension(:),pointer,optional :: caggrnum ! coarse aggr# for each aggr
00291     integer,    intent(in), optional   :: INIT_CONT_END
00292     real(kind=xyzk),dimension(:),pointer :: xc, yc
00293     real(kind=xyzk)                      :: xmin, xmax, ymin, ymax
00294     integer                            :: nc,i,j,jj,nr,c,ani,anj
00295     character*2                        :: buf2
00296     character*3                        :: buf3
00297     character*5                        :: buf5
00298     character*10                       :: buf10
00299 #ifdef D_WANT_PLPLOT_YES
00300     neighood  = aggr%radius
00301     nagrs     = aggr%nagr
00302     nisolated = aggr%nisolated
00303     aggrnum   = aggr%num
00304     write(stream, *)
00305     write(stream, *) '[Mesh_pl2D_plotAggregate] : Plotting 2D mesh.'
00306     if (M%nsd /= 2) then
00307        write(stream, *) '   Only for 2D meshes. Skipping...'
00308        return
00309     end if
00310     nc = size(M%coords)
00311     xc => M%coords(1,:)
00312     yc => M%coords(2,:)
00313     if (.not.present(INIT_CONT_END).or.&
00314          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00315        xmin = minval(xc)
00316        xmax = maxval(xc)
00317        ymin = minval(yc)
00318        ymax = maxval(yc)
00319        xmin = xmin - xmax/10.0
00320        xmax = xmax + xmax/10.0
00321        ymin = ymin - ymax/10.0
00322        ymax = ymax + ymax/10.0
00323        write(buf2, '(i2)') neighood
00324        write(buf10, '(e10.2)') sctls%strong1
00325        write(buf5, '(i5)') nagrs
00326        write(buf3, '(i3)') nisolated
00327        print *,trim(filename(1:index(filename,'.',.true.)-1))// &
00328              '_'//trim(buf2(index(buf2,' ',.true.)+1:))//'_'// &
00329               trim(buf10(index(buf10,' ',.true.)+1:))//'_'// &
00330               trim(buf5(index(buf5,' ',.true.)+1:))//'_'// &
00331               trim(buf3(index(buf3,' ',.true.)+1:))
00332        if (sctls%plotting==1) then
00333          call plsdev("xwin")
00334        else
00335          call plsdev("tk")
00336        endif
00337        call plinit()
00338        ! scale fontsize:
00339        call plschr(0.d0,0.5d0) ! d is crucial!
00340        call plenv (xmin, xmax, ymin, ymax, 0, 0);
00341        call plcol0(1) ! red
00342        call pllab( '(x)', '(y)', &
00343             trim(filename)// & !'; ngf='//buf61// &
00344              ' rad='//buf2//' thr='//buf10 &
00345               //' Na='//buf5//' isl='//buf3)
00346     end if
00347     call plssym(0.0d0, 2.0d0)
00348     call plpoin(nc/2, xc, yc, 1)
00349     do i=1,nc/2
00350       ! Plot vertices
00351       call plcol0(15) ! white
00352       call plssym(0.0d0, 2.0d0)
00353       call plpoin(1,xc(i),yc(i),1)
00354     enddo
00355     nr=size(rowstart)-1
00356     !print *,'nnnn nr=',nr
00357     do i=1,nr
00358       ani=aggrnum(i)
00359       if (ani/=0) then
00360         if (present(caggrnum)) then
00361           ani=caggrnum(ani)
00362         endif
00363         if (ani/=0) then
00364           c=1+modulo(ani,13)
00365           if (c==7) c=c+1
00366           call plcol0(c) ! cycling colors...
00367           do j=rowstart(i),rowstart(i+1)-1
00368             jj=colnrs(j)
00369             anj=aggrnum(jj)
00370             if (anj/=0) then
00371               if (present(caggrnum)) then
00372                 anj=caggrnum(anj)
00373               endif
00374               if (ani==anj.and.aggrnum(i)==aggrnum(jj)) then
00375                 call plline(2, &
00376                          (/xc(M%freemap(i)),xc(M%freemap(jj))/), &
00377                          (/yc(M%freemap(i)),yc(M%freemap(jj))/)  )
00378               endif
00379             endif
00380           enddo
00381         endif
00382       endif
00383     end do
00384     if (.not.present(INIT_CONT_END).or.&
00385          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_END))) then
00386        call plend()
00387     end if
00388     nullify(xc, yc)
00389 #else
00390     write(stream, '(/a)') ' [Mesh_pl2D_plotAggregate] : Compiled w/o plotting'//&
00391          ' support!'
00392 #endif
00393   end subroutine Mesh_pl2D_plotAggregate
00394 
00395 
00396   !---------------------------------------------------
00397   !! Plot mesh's dual graph
00398   !---------------------------------------------------
00399   subroutine Mesh_pl2D_plotGraphDual(M, INIT_CONT_END)
00400 
00401     use globals, only : stream
00402 
00403     implicit none
00404 
00405     type(Mesh), intent(in out)               :: M
00406     integer,    intent(in),    optional      :: INIT_CONT_END
00407 
00408     real(kind=xyzk), dimension(:), pointer   :: xc, yc ! all node coordinates
00409     real(kind=xyzk)                          :: xmin, xmax, ymin, ymax
00410     integer                                  :: vplotted
00411     integer                                  :: n, i, j
00412     character*5                              :: buf, buf51, buf52, buf53
00413 
00414     ! Elements' centres of mass
00415     real(kind=xyzk), dimension(:,:), allocatable :: elCentrMass
00416 
00417     integer                                  :: e, en_i ! elements counters
00418     type(Polygon)                            :: e_p
00419     type(Points2D)                           :: e_cntr
00420     ! auxiliary var. for 2-nodes element
00421     real(kind=xyzk)                            :: e_cntr_x, e_cntr_y
00422     ! indexes to get element node coordinates from 'xc', 'yc'
00423     integer,       dimension(:), allocatable :: ind_coords
00424     ! same for the neighbouring element
00425     integer,       dimension(:), allocatable :: ind_coords_en
00426     ! Counter for plotted graph nodes
00427     integer                                  :: graphNodesPlotted
00428     ! Accounting for drawn graph edges
00429     integer,       dimension(:), allocatable :: graphEdgesplotted
00430     logical                                  :: plotted=.false.
00431 
00432 #ifdef D_WANT_PLPLOT_YES
00433 
00434     write(stream, *)
00435     write(stream, *) '[Mesh_pl2D_plotGraphDual] : Plotting 2D mesh''s dual'//&
00436          ' graph.'
00437 
00438     if (M%nsd /= 2) then
00439        write(stream, *) '   Only for 2D meshes. Skipping...'
00440        return
00441     end if
00442 
00443     if (.not.associated(M%coords)) then
00444        write(stream, *) '   Mesh hasn''t nodes'' coordinates. Skipping...'
00445        return
00446     end if
00447 
00448 
00449     n = size(M%coords)
00450 
00451     xc => M%coords(1,:)
00452     yc => M%coords(2,:)
00453 
00454     if (.not.present(INIT_CONT_END).or.&
00455          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00456 
00457        xmin = minval(xc)
00458        xmax = maxval(xc)
00459        ymin = minval(yc)
00460        ymax = maxval(yc)
00461 
00462        xmin = xmin - xmax/10.0
00463        xmax = xmax + xmax/10.0
00464        ymin = ymin - ymax/10.0
00465        ymax = ymax + ymax/10.0
00466 
00467        call plsdev("xwin")
00468        call plinit()
00469 
00470        call plenv (xmin, xmax, ymin, ymax, 0, 0)
00471     end if
00472 
00473     ! Plot mesh nodes
00474     call plcol0(15) ! white
00475     call plssym(0.0d0, 2.0d0)
00476     call plpoin(n/2, xc, yc, 1)
00477 
00478     ! Coordinates' indexes
00479     allocate(ind_coords(M%mfrelt), ind_coords_en(M%mfrelt), &
00480          graphEdgesPlotted(size(M%G%adjncy)))
00481     ind_coords = 0
00482 
00483     graphEdgesPlotted = 0
00484     graphNodesPlotted = 0
00485 
00486 
00487     ! Find elements' centres of mass
00488     allocate(elCentrMass(M%nell,2))
00489     do e = 1,M%nell
00490 
00491        ! Get global node numbering to fetch nodes'
00492        ! global coordinate values later.
00493        do i = 1,M%nfrelt(e)
00494           ind_coords(i) = M%freemap(M%mhead(i,e))
00495        end do
00496 
00497        ! Calculate centre of mass of the element
00498        if (M%nfrelt(e) == 1) then ! 1-node boundary element
00499 
00500           ! Fake centre of 1-node element by simply assigning to
00501           ! it coordinates of the node itself
00502           elCentrMass(e,1) = xc(ind_coords(1))
00503           elCentrMass(e,2) = yc(ind_coords(1))
00504 
00505           ! Find my neighbours and calculate the centre of mass accordingly
00506           ! This is not necessary!
00507 
00508        else if (M%nfrelt(e) == 2) then ! two-nodes boundary element
00509 
00510           ! Fake centre of 2-node element by simply assigning to
00511           ! it coordinates of the middle point between element nodes
00512           elCentrMass(e,1) = xc(ind_coords(1)) + (xc(ind_coords(2)) - &
00513                xc(ind_coords(1))) / 2.0_xyzk
00514           elCentrMass(e,2) = yc(ind_coords(1)) + (yc(ind_coords(2)) - &
00515                yc(ind_coords(1))) / 2.0_xyzk
00516 
00517        else if (M%nfrelt(e) >= 3) then
00518 
00519           e_p = Polygon_New(M%nfrelt(e))
00520           call Polygon_Init(e_p, &
00521                xc(ind_coords(1:M%nfrelt(e))), yc(ind_coords(1:M%nfrelt(e))) )
00522 
00523           e_cntr = Points2D_New(1)
00524           e_cntr = Polygon_Centroid(e_p)
00525 
00526           elCentrMass(e,1) = e_cntr%x(1)
00527           elCentrMass(e,2) = e_cntr%y(1)
00528 
00529           call Polygon_Destroy(e_p)
00530           call Points2D_Destroy(e_cntr)
00531 
00532        end if
00533 
00534     end do
00535 
00536 
00537     ! Do main draw loop
00538     do e = 1,M%nell
00539 
00540        ! Plot point in the centre of mass of the element
00541        call plcol0(2) ! yellow
00542        call plssym(0.0d0, 2.0d0)
00543        call plpoin(1, elCentrMass(e,1), elCentrMass(e,2), 1)
00544        graphNodesPlotted = graphNodesPlotted + 1
00545 
00546        ! Draw edges with our neighbours
00547        do i = M%G%xadj(e),M%G%xadj(e+1)-1
00548 
00549           ! Index of neighbour element
00550           en_i = M%G%adjncy(i)
00551 
00552           ! Don't draw edges with nodes which has already drawn them with us
00553           do j = M%G%xadj(e),M%G%xadj(e+1)-1
00554              if ((M%G%adjncy(j) == en_i).and.(graphEdgesPlotted(j) == 1)) then
00555                 plotted = .true.
00556                 exit
00557              else
00558                 plotted = .false.
00559              end if
00560           end do
00561 
00562           if (.not.plotted) then
00563 
00564              ! Only if number of nodes more than two
00565              if (M%nfrelt(en_i) > M%nsd) then
00566 
00567                 ! Get global node numbering to fetch nodes'
00568                 ! global coordinate values later.
00569                 ! Inices of nodes for neighbouring element
00570                 do j = 1,M%nfrelt(en_i)
00571                    ind_coords_en(j) = M%freemap(M%mhead(j,en_i))
00572                 end do
00573 
00574                 ! Actually plot edge between centres of adjacent elements
00575                 call plcol0(9) ! Blue
00576                 call plline(2, &
00577                      (/elCentrMass(e,1),elCentrMass(en_i,1)/), &
00578                      (/elCentrMass(e,2),elCentrMass(en_i,2)/)  )
00579 
00580                 ! Find connection with the element which
00581                 ! is referencing us right now
00582                 ! and file us that we were plotted.
00583                 do j = M%G%xadj(en_i),M%G%xadj(en_i+1)-1
00584                    if (M%G%adjncy(j) == e) graphEdgesPlotted(j) = 1
00585                 end do
00586 
00587              else
00588 !!$                write(stream, '(a,i5,a,i5,a)') '  Adjacency for ', M%nfrelt(en_i),'-node(s) elem  [', en_i,&
00589 !!$                     ']. Not implemented yet. Skipping...'
00590              end if
00591           end if
00592        end do
00593 
00594        call Polygon_Destroy(e_p)
00595        call Points2D_Destroy(e_cntr)
00596 
00597        ind_coords = 0
00598     end do
00599 
00600     ! Plot figure title
00601     if (.not.present(INIT_CONT_END).or.&
00602          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00603        call plcol0(1) ! red
00604        write(buf51, '(i5)') M%G%nvtx
00605        write(buf52, '(i5)') graphNodesPlotted
00606        write(buf53, '(i5)') sum(graphEdgesPlotted)
00607        call pllab( '(x)', '(y)', 'Graph : '//buf51//&
00608             ' nodes. Plotted: nodes='//buf52//'; edges='//buf53//'.' )
00609     end if
00610 
00611     if (.not.present(INIT_CONT_END).or.&
00612          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_END))) then
00613        call plend()
00614     end if
00615 
00616     nullify(xc, yc)
00617     deallocate(ind_coords,  &
00618          ind_coords_en,     &
00619          graphEdgesPlotted, &
00620          elCentrMass)
00621 
00622 #else
00623     write(stream, '(/a)') ' [Mesh_pl2D_plotGraphDual] : Compiled w/o'//&
00624          ' plotting support!'
00625 #endif
00626 
00627   end subroutine Mesh_pl2D_plotGraphDual
00628 
00629 
00630   !-----------------------------------------------
00631   !! Plot 2D mesh partition
00632   !-----------------------------------------------
00633   subroutine Mesh_pl2D_Partition(M, INIT_CONT_END)
00634 
00635     use globals, only : stream
00636 
00637     implicit none
00638 
00639     type(Mesh), intent(in out)               :: M
00640     integer,    intent(in),    optional      :: INIT_CONT_END
00641 
00642     real(kind=xyzk), dimension(:), pointer   :: xc, yc
00643     real(kind=xyzk)                          :: xmin, xmax, ymin, ymax
00644     integer                                  :: n, e, i
00645     character*4                              :: buf4
00646 
00647     type(Polygon)                            :: ep
00648     real(kind=xyzk)                          :: ntmp_y
00649     ! indices to get element node coordinates from 'xc', 'yc'
00650     integer,       dimension(:), allocatable :: ind_coords
00651     integer                                  :: partcolor
00652 
00653 #ifdef D_WANT_PLPLOT_YES
00654 
00655     write(stream, *)
00656     write(stream, *) '[Mesh_pl2D_partns] : Plotting 2D mesh partition.'
00657 
00658     if (M%nsd /= 2) then
00659        write(stream, *) '   Only for 2D meshes. Skipping...'
00660        return
00661     end if
00662 
00663     if (M%parted.eqv.(.false.)) then
00664        write(stream, *) '   Mesh was not partitioned. Skipping...'
00665        return
00666     end if
00667 
00668 
00669     n = size(M%coords)
00670 
00671     xc => M%coords(1,:)
00672     yc => M%coords(2,:)
00673 
00674     if (.not.present(INIT_CONT_END).or.&
00675          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00676 
00677        xmin = minval(xc)
00678        xmax = maxval(xc)
00679        ymin = minval(yc)
00680        ymax = maxval(yc)
00681 
00682        xmin = xmin - xmax/10.0
00683        xmax = xmax + xmax/10.0
00684        ymin = ymin - ymax/10.0
00685        ymax = ymax + ymax/10.0
00686 
00687        call plsdev("xwin")
00688        call plinit()
00689 
00690        call plenv (xmin, xmax, ymin, ymax, 0, 0);
00691 
00692        call plcol0(1) ! red
00693        write(buf4, '(i4)') M%nparts
00694        call pllab( '(x)', '(y)', 'Mesh : '//buf4//' partitions' )
00695     end if
00696 
00697     ! Plot nodes
00698     call plcol0(15) ! white
00699     call plssym(0.0d0, 2.0d0)
00700     call plpoin(n/2, xc, yc, 1)
00701 
00702     allocate(ind_coords(M%mfrelt))
00703     ind_coords = 0
00704 
00705     do e = 1,M%nell
00706 
00707        if (M%nfrelt(e) > M%nsd) then
00708 
00709           ! Get global node numbering to fetch nodes'
00710           ! global coordinate values later.
00711           do i = 1,M%nfrelt(e)
00712              ind_coords(i) = M%freemap(M%mhead(i,e))
00713           end do
00714 
00715           if (M%nfrelt(e) == 3) then
00716              ! NB: Just to speed up plotting
00717              !
00718              ! Trivial 3-node (triangle) element case
00719              !
00720              ! There can be more than 3-nodes elements on the boundary
00721              ! with Dirichlet BC, having only 3 nodes left on one line.
00722              ! Check it too! TODO
00723 
00724              ! Choose colour according to partition number
00725              ! cycle colours : 1..15, 1..15,...
00726              call plcol0(1 + mod(M%eptnmap(e)-1,15))
00727              ! Plot filled polygon
00728              call plfill(M%nfrelt(e), &
00729                   xc(ind_coords(1:3)), yc(ind_coords(1:3)) )
00730 
00731              ! Plot polygon
00732              call plcol0(3) ! green
00733              ! "Close" 3-node elements to form polygons
00734              call plline(4, &
00735                   (/xc(ind_coords(1:3)),xc(ind_coords(1))/), &
00736                   (/yc(ind_coords(1:3)),yc(ind_coords(1))/) )
00737 
00738              ! Plot vertices
00739              call plcol0(2) ! yellow
00740              call plssym(0.0d0, 2.0d0)
00741              call plpoin(M%nfrelt(e), &
00742                   xc(ind_coords(1:3)), yc(ind_coords(1:3)), 1)
00743 
00744           else if (M%nfrelt(e) > 3) then ! More than 3 is a little bit tricky
00745              ep = Polygon_New(M%nfrelt(e))
00746 
00747              call Polygon_Init(ep, &
00748                   xc(ind_coords(1:M%nfrelt(e))), yc(ind_coords(1:M%nfrelt(e))))
00749              call Polygon_sortVerts(ep)
00750 
00751              ! Choose colour according to partition number
00752              ! cycle colours : 1..15, 1..15,...
00753              call plcol0(1 + mod(M%eptnmap(e)-1,15))
00754              call plfill(M%nfrelt(e), Polygon_getX(ep), Polygon_getY(ep) )
00755              call Polygon_pl2D_Plot(ep, D_PLPLOT_CONT)
00756 
00757              call Polygon_Destroy(ep)
00758           end if
00759 
00760        end if
00761 
00762     end do
00763 
00764     if (.not.present(INIT_CONT_END).or.&
00765          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_END))) then
00766        call plend()
00767     end if
00768 
00769     nullify(xc, yc)
00770     deallocate(ind_coords)
00771 
00772 #else
00773     write(stream, '(/a)') ' [Mesh_pl2D_Partition] : Compiled w/o plotting'//&
00774          ' support!'
00775 #endif
00776 
00777   end subroutine Mesh_pl2D_Partition
00778 
00779 
00780   !-----------------------------------------------------
00781   !! Plot coloured dual graph of partitioned mesh
00782   !-----------------------------------------------------
00783   subroutine Mesh_pl2D_plotGraphParted(M, INIT_CONT_END)
00784 
00785     use globals, only : stream
00786 
00787     implicit none
00788 
00789     type(Mesh), intent(in out)               :: M
00790     integer,    intent(in),    optional      :: INIT_CONT_END
00791 
00792     real(kind=xyzk), dimension(:), pointer   :: xc, yc ! all node coordinates
00793     real(kind=xyzk)                          :: xmin, xmax, ymin, ymax
00794     integer                                  :: n, i, j
00795     character*4                              :: buf4
00796 
00797     integer                                  :: e, en_i
00798     type(Polygon)                            :: e_p, en_p
00799     type(Points2D)                           :: e_cntr, en_cntr
00800     ! aux. var. for 2-nodes element
00801     real(kind=xyzk)                          :: e_cntr_x, e_cntr_y
00802     ! aux. var. for two-segment edges
00803     real(kind=xyzk)                          :: half_x, half_y
00804     ! indices to get element node coordinates from 'xc', 'yc'
00805     integer,       dimension(:), allocatable :: ind_coords
00806     ! same for the neighbouring element
00807     integer,       dimension(:), allocatable :: ind_coords_en
00808     ! Counter for plotted graph nodes
00809     integer                                  :: graphNodesPlotted
00810     ! Accounting for drawn graph edges
00811     integer,       dimension(:), allocatable :: graphEdgesplotted
00812     logical                                  :: plotted=.false.
00813     character*2                              :: buf2
00814 
00815 #ifdef D_WANT_PLPLOT_YES
00816 
00817     write(stream, *)
00818     write(stream, *) '[Mesh_pl2D_plotGraphParted] : Plotting 2D'//&
00819          ' partitioned graph.'
00820 
00821     if (M%nsd /= 2) then
00822        write(stream, *) '   Only for 2D meshes. Skipping...'
00823        return
00824     end if
00825 
00826     if (M%parted.eqv.(.false.)) then
00827        write(stream, *) '   Mesh was not partitioned. Skipping...'
00828        return
00829     end if
00830 
00831 
00832     n = size(M%coords)
00833 
00834     xc => M%coords(1,:)
00835     yc => M%coords(2,:)
00836 
00837     if (.not.present(INIT_CONT_END).or.&
00838          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00839 
00840        xmin = minval(xc)
00841        xmax = maxval(xc)
00842        ymin = minval(yc)
00843        ymax = maxval(yc)
00844 
00845        xmin = xmin - xmax/10.0
00846        xmax = xmax + xmax/10.0
00847        ymin = ymin - ymax/10.0
00848        ymax = ymax + ymax/10.0
00849 
00850        call plsdev("xwin")
00851        call plinit()
00852 
00853        call plenv (xmin, xmax, ymin, ymax, 0, 0);
00854 
00855        call plcol0(1) ! red
00856        write(buf4, '(i4)') M%nparts
00857        call pllab( '(x)', '(y)', 'Graph : '//buf4//' partitions' )
00858     end if
00859 
00860     ! Plot nodes
00861     call plcol0(15) ! white
00862     call plssym(0.0d0, 2.0d0)
00863     call plpoin(n/2, xc, yc, 1)
00864 
00865     ! Coordinates' indices
00866     allocate(ind_coords(M%mfrelt),           &
00867          ind_coords_en(M%mfrelt),            &
00868          graphEdgesPlotted(size(M%G%adjncy)) )
00869 
00870     ind_coords = 0
00871     graphEdgesPlotted = 0
00872     graphNodesPlotted = 0
00873 
00874 
00875     do e = 1,M%nell
00876 
00877        ! Get global node numbering to fetch nodes'
00878        ! global coordinate values later.
00879        do i = 1,M%nfrelt(e)
00880           ind_coords(i) = M%freemap(M%mhead(i,e))
00881        end do
00882 
00883        if (M%nfrelt(e) == 0) then
00884           cycle
00885 
00886        else if (M%nfrelt(e) == 1) then ! 1-node boundary element
00887 
00888           ! Fake centre of 1-node element by simply assigning to
00889           ! it coordinates of the node itself
00890           e_cntr = Points2D_newFill(&
00891                (/xc(ind_coords(1))/),(/yc(ind_coords(1))/))
00892 
00893        else if (M%nfrelt(e) == 2) then ! two-nodes boundary element
00894 
00895           ! Fake centre of 2-node element by simply assigning to
00896           ! it coordinates of the middle point between element nodes
00897           e_cntr_x = xc(ind_coords(1)) + (xc(ind_coords(2)) - &
00898                xc(ind_coords(1))) / 2.0_xyzk
00899           e_cntr_y = yc(ind_coords(1)) + (yc(ind_coords(2)) - &
00900                yc(ind_coords(1))) / 2.0_xyzk
00901           e_cntr = Points2D_newFill((/e_cntr_x/), (/e_cntr_y/))
00902 
00903        else if (M%nfrelt(e) >= 3) then
00904 
00905           e_p = Polygon_New(M%nfrelt(e))
00906           call Polygon_Init(e_p, &
00907                xc(ind_coords(1:M%nfrelt(e))), yc(ind_coords(1:M%nfrelt(e))) )
00908 
00909           e_cntr = Points2D_New(1)
00910           e_cntr = Polygon_Centroid(e_p)
00911 
00912        end if
00913 
00914        ! Plot point in the centre of mass of the element
00915        ! Choose colour according to partition number
00916        ! cycle colours : 1..15, 1..15,...
00917        call plcol0(1+mod(M%eptnmap(e)-1,15))
00918        call plssym(0.0d0, 2.0d0)
00919        call plpoin(1, e_cntr%x, e_cntr%y, 1)
00920        if (M%nnode <= 100) then
00921           call plssym(0.0d0, 5.0d0)
00922           write(buf2, '(i2)') e
00923           call plptex(e_cntr%x, e_cntr%y, 0.0, 0.0, 0, buf2)
00924        end if
00925        graphNodesPlotted = graphNodesPlotted + 1
00926 
00927        ! Find centres of masses of our neighbours
00928        do i = M%G%xadj(e),M%G%xadj(e+1)-1
00929 
00930           ! Index of neighbour element
00931           en_i = M%G%adjncy(i)
00932 
00933           ! Don't draw edges with nodes which has already drawn them with us
00934           do j = M%G%xadj(e),M%G%xadj(e+1)-1
00935              if ((M%G%adjncy(j) == en_i).and.(graphEdgesPlotted(j) == 1)) then
00936                 plotted = .true.
00937                 exit
00938              else
00939                 plotted = .false.
00940              end if
00941           end do
00942 
00943           if (.not.plotted) then
00944              ! Only if number of nodes more than two
00945              if (M%nfrelt(en_i) > M%nsd) then
00946 
00947                 ! Get global node numbering to fetch nodes'
00948                 ! global coordinate values later.
00949                 ! Inices of nodes for neighbouring element
00950                 do j = 1,M%nfrelt(en_i)
00951                    ind_coords_en(j) = M%freemap(M%mhead(j,en_i))
00952                 end do
00953                 ! Find its centre of mass
00954                 en_p = Polygon_New(M%nfrelt(en_i))
00955                 call Polygon_Init(en_p, xc(ind_coords_en(1:M%nfrelt(en_i))), &
00956                      yc(ind_coords_en(1:M%nfrelt(en_i))) )
00957                 en_cntr = Points2D_New(1)
00958                 en_cntr = Polygon_Centroid(en_p)
00959 
00960                 ! Actually plot edge between centres of adjacent elements
00961                 ! Choose colour according to partition number
00962                 ! cycle colours : 1..15, 1..15,...
00963                 call plcol0(1+mod(M%eptnmap(e)-1,15))
00964 
00965                 if (M%eptnmap(e) == M%eptnmap(en_i)) then
00966                    ! Elements belong to one partition
00967                    call plline(2, &
00968                         (/e_cntr%x,en_cntr%x/), (/e_cntr%y,en_cntr%y/))
00969                 else
00970                    ! Elements belong to different partitions
00971                    half_x = e_cntr%x(1) + (en_cntr%x(1) - e_cntr%x(1)) /2.0_xyzk
00972                    half_y = e_cntr%y(1) + (en_cntr%y(1) - e_cntr%y(1)) /2.0_xyzk
00973 
00974                    ! Draw two-segment coloured edge
00975                    call plcol0(1+mod(M%eptnmap(e)-1,15))
00976                    call plline(2, (/e_cntr%x,half_x/), (/e_cntr%y,half_y/))
00977                    call plcol0(1+mod(M%eptnmap(en_i)-1,15))
00978                    call plline(2, (/half_x,en_cntr%x/), (/half_y,en_cntr%y/))
00979                 end if
00980 
00981 
00982                 ! Find connection with the element is referencing us now
00983                 ! and file us that we were plotted.
00984                 !vplotted = vplotted + 1
00985                 do j = M%G%xadj(en_i),M%G%xadj(en_i+1)-1
00986                    if (M%G%adjncy(j) == e) graphEdgesPlotted(j) = 1
00987                 end do
00988 
00989                 call Polygon_Destroy(en_p)
00990                 call Points2D_Destroy(en_cntr)
00991              else
00992 !!$                write(stream, '(a,i5,a,i5,a)') '  Adjacency for ', M%nfrelt(en_i),'-node(s) elem  [', en_i,&
00993 !!$                     ']. Not implemented yet. Skipping...'
00994              end if
00995           end if
00996        end do
00997 
00998        call Polygon_Destroy(e_p)
00999        call Points2D_Destroy(e_cntr)
01000 
01001        ind_coords = 0
01002     end do
01003 
01004     ! Print out some statistics
01005     write(stream, '(a,i5,a,i5,a,i5)') '  Graph : ',M%G%nvtx, &
01006          ' nodes. Plotted: nodes=',graphNodesPlotted, &
01007          '; edges=',sum(graphEdgesPlotted)
01008 
01009     if (.not.present(INIT_CONT_END).or.&
01010          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_END))) then
01011        call plend()
01012     end if
01013 
01014     nullify(xc, yc)
01015     deallocate(ind_coords, ind_coords_en, graphEdgesPlotted)
01016 
01017 #else
01018     write(stream, '(/a)') ' [Mesh_pl2D_plotGraphParted] : Compiled w/o'//&
01019          ' plotting support!'
01020 #endif
01021 
01022   end subroutine Mesh_pl2D_plotGraphParted
01023 
01024 
01025 end module Mesh_plot_mod