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