GRIDMAN
grid managment library
grid2d.f
Go to the documentation of this file.
1 C> @file grid2D/grid2d.f
2 C> Basic subroutines, calculation of geometrical properties
3 C GRIDMAN, grid managment library. Author: Vladislav Kotov, v.kotov@fz-juelich.de
4 
5 ! Copyright (c) 2017 Forschungszentrum Juelich GmbH
6 ! Vladislav Kotov
7 !
8 ! This file is part of GRIDMAN.
9 !
10 ! GRIDMAN is free software: you can redistribute it and/or modify
11 ! it under the terms of the GNU General Public License as published by
12 ! the Free Software Foundation, either version 3 of the License, or
13 ! (at your option) any later version.
14 !
15 ! GRIDMAN is distributed in the hope that it will be useful,
16 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ! GNU General Public License for more details.
19 !
20 ! You should have received a copy of the GNU General Public License
21 ! along with GRIDMAN. If not, see <http://www.gnu.org/licenses/>.
22 
23 C*****************************************************************************************
24 C> Check correctness of the 2D grid object
25 C*****************************************************************************************
26 C>
27 C> - RES<100: see GRIDMAN_GRID_CHECK in grid.f
28 C> - RES=100: the grid is not of 2D-type
29 C> - RES=102: the chain of edges belonging to one cell is brocken
30 C> or no single chain
31 C> - RES=103: edges with zero length
32 C> - RES=111: intersecting edges
33 C> (only edges withing one cell are checked)
34 C> (edges which do not belong to any cell are not checked)
35 C
36  SUBROUTINE gridman_grid2d_check(GRID,RES,IERR)
37  USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
38  u gridman_dbg,gridman_indlist
39  USE gridman_lib,ONLY:gridman_grid_check,gridman_intersect2d,
40  u gridman_grid_cells,gridman_grid_points,
41  u gridman_indlist_deallocate
42  IMPLICIT NONE
43  INTRINSIC tiny,max,min
44 C> Grid object to be checked
45  TYPE(gridman_grid) :: GRID
46 C> Code of the result
47  INTEGER,INTENT(OUT) :: RES
48 C> Error code
49  INTEGER,INTENT(OUT) :: IERR
50 
51  INTEGER(GRIDMAN_SP) :: ICELL,IE,IE1,I,
52  i ipoint,ipoint1,ipoint2,ipoint0,
53  i ipoint11,ipoint12,ipoint21,ipoint22,
54  i itmp,iedge0,iedge,iedge1,iedge2
55  INTEGER :: IERR0
56  REAL(GRIDMAN_DP) :: X11,Y11,X12,Y12,X21,Y21,X22,Y22,XI,YI,
57  r x1,x2,y1,y2,s
58  TYPE(gridman_indlist) :: POINTS,CELLS
59  LOGICAL :: LPASS
60 
61  IF(gridman_dbg)
62  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_CHECK"
63  ierr=0
64  res=0
65 
66  IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
67  res=100
68  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CHECK: ",
69  w "the grid is not of 2D type"
70  WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
71  w grid%TYPE,grid%PDIM,grid%EDIM
72  RETURN
73  END IF
74 
75  CALL gridman_grid_check(grid,res,ierr)
76  IF(res.NE.0.OR.ierr.GT.0) RETURN
77 
78  CALL gridman_grid_cells(cells,grid,ierr)
79  IF(ierr.NE.0) GOTO 1000
80  CALL gridman_grid_points(points,grid,ierr)
81  IF(ierr.NE.0) GOTO 500
82 C
83 C CHECK THAT ALL CHAINS ARE CLOSED AND SINGLE
84 C
85  DO icell=1,grid%NCELLS
86 C TAKE FIRST EDGE OF THE CELL
87  iedge0=cells%IND(cells%IFIRST(icell))
88  IF(iedge0.LT.1.OR.iedge0.GT.grid%NEDGES) THEN
89  WRITE(gridman_unit,*)
90  w "ERROR in GRIDMAN_GRID2D_CHECK:",
91  w " internal error - wrong edge index"
92  WRITE(gridman_unit,*) " ICELL, IEDGE0 ",icell,iedge
93  GOTO 100
94  END IF
95  ipoint0=grid%POINTS(1,iedge0) !FIRST POINT IN THE CHAIN
96  ipoint=grid%POINTS(2,iedge0) !CURRENT POINT IN THE CHAIN
97 C NUMBER OF ITERATIONS IN THE CYCLE IS EQUAL TO THE NUMBER
98 C OF EDGES BELONGING TO THE CELL MINUS ONE
99  ipoint2=0
100  DO i=cells%IFIRST(icell)+1,cells%ILAST(icell)
101 C LOOK FOR THE NEXT EDGE IN THE CHAIN: EDGE CONENCTED TO
102 C THE CURRENT POINT IPOINT AND BELONGING TO THE CURRENT CELL ICELL
103  lpass=.false.
104  DO ie=points%IFIRST(ipoint),points%ILAST(ipoint)
105  iedge=points%IND(ie)
106  IF(iedge.NE.iedge0) THEN !IF NOT PREVIOUS EDGE
107  IF(iedge.LT.1.OR.iedge.GT.grid%NEDGES) THEN
108  WRITE(gridman_unit,*)
109  w "ERROR in GRIDMAN_GRID2D_CHECK:",
110  w " internal error - wrong edge index"
111  WRITE(gridman_unit,*) " ICELL, IEDGE ",icell,iedge
112  GOTO 100
113  END IF
114  IF(grid%CELLS(1,iedge).EQ.icell.OR.
115  f grid%CELLS(2,iedge).EQ.icell) THEN
116 C NEXT EDGE IS FOUND
117  ipoint1=grid%POINTS(1,iedge)
118  ipoint2=grid%POINTS(2,iedge)
119  IF(lpass) THEN
120  res=102
121  WRITE(gridman_unit,*) " GRIDMAN_GRID2D_CHECK: ",
122  w "more than two edges ",
123  w "connected to a point in one cell"
124  WRITE(gridman_unit,*) " ICELL, IPOINT1, IPOINT2 ",
125  w icell, ipoint1, ipoint2
126  WRITE(gridman_unit,*) " IEDGE0, IEDGE ", iedge0,iedge
127  GOTO 100
128  END IF
129  iedge2=iedge
130  lpass=.true.
131  IF(ipoint1.NE.ipoint) THEN
132  itmp=ipoint2
133  ipoint2=ipoint1
134  ipoint1=itmp
135  IF(ipoint1.NE.ipoint) THEN
136  ierr=400
137  WRITE(gridman_unit,*)
138  w "ERROR in GRIDMAN_GRID2D_CHECK:",
139  w " internal error - the found edge is wrong"
140  WRITE(gridman_unit,*)
141  w " ICELL, IEDGE, IPOINT, IPOINT1, IPOINT2 ",
142  w icell, iedge, ipoint, ipoint1, ipoint2
143  GOTO 100
144  END IF
145  END IF
146  END IF
147  END IF !IF(IEDGE.NE.IEDGE0)
148  END DO !DO IE=POINTS%IFIRST(IPOINT)+1,POINTS%ILAST(IPOINT)
149 C NOW: IPOINT1=IPOINT ------------ IPOINT2
150 C SHIFT IPOINT AND IEDGE0 TO THE NEXT POSITION
151  IF(ipoint2.EQ.0) THEN
152  res=102
153  WRITE(gridman_unit,*) " GRIDMAN_GRID2D_CHECK: ",
154  w "the chain is brocken"
155  WRITE(gridman_unit,*) " ICELL, IPOINT, IPOINT0 ",
156  w icell, ipoint, ipoint0
157  GOTO 100
158  END IF
159  ipoint=ipoint2
160  iedge0=iedge2
161  END DO !DO I=CELLS%IFIRST(ICELL)+1,CELLS%ILAST(ICELL)
162 C FIRST AND LAST POINT OF THE CHAIN MUST BE SAME
163  IF(ipoint0.NE.ipoint) THEN
164  res=102
165  WRITE(gridman_unit,*) " GRIDMAN_GRID2D_CHECK: ",
166  w "the chain is brocken"
167  WRITE(gridman_unit,*) " ICELL, IPOINT, IPOINT0 ",
168  w icell, ipoint, ipoint0
169  GOTO 100
170  END IF
171  END DO !DO ICELL=1,GRID%NCELLS
172 
173 C
174 C CHECK THAT THE EDGE LENGTH IS NOT ZERO
175 C
176  DO iedge=1,grid%NEDGES
177  ipoint1=grid%POINTS(1,iedge)
178  ipoint2=grid%POINTS(2,iedge)
179  x1=grid%X(1,ipoint1)
180  y1=grid%X(2,ipoint1)
181  x2=grid%X(1,ipoint2)
182  y2=grid%X(2,ipoint2)
183  s=(x2-x1)**2+(y2-y1)**2
184  IF(.NOT.s.GT.tiny(x1)) THEN
185  res=103
186  WRITE(gridman_unit,*) " GRIDMAN_GRID2D_CHECK: ",
187  w "edge with zero length"
188  WRITE(gridman_unit,*) " IEDGE ",iedge
189  WRITE(gridman_unit,*) " X1, Y1, X2, Y2 ",x1,y1,x2,y2
190  GOTO 100
191  END IF
192  END DO
193 
194 C
195 C CHECK FOR INTERSECTION BETWEEN EDGES WITHING ONE CELL
196 C
197  DO icell=1,grid%NCELLS
198  DO ie1=cells%IFIRST(icell),cells%ILAST(icell)
199  iedge1=cells%IND(ie1)
200  ipoint11=grid%POINTS(1,iedge1)
201  ipoint12=grid%POINTS(2,iedge1)
202  x11=grid%X(1,ipoint11)
203  y11=grid%X(2,ipoint11)
204  x12=grid%X(1,ipoint12)
205  y12=grid%X(2,ipoint12)
206  DO ie=ie1+1,cells%ILAST(icell)
207  iedge=cells%IND(ie)
208  ipoint21=grid%POINTS(1,iedge)
209  ipoint22=grid%POINTS(2,iedge)
210  IF(ipoint21.EQ.ipoint11.OR.ipoint21.EQ.ipoint12.OR.
211  f ipoint22.EQ.ipoint11.OR.ipoint22.EQ.ipoint12) cycle
212  x21=grid%X(1,ipoint21)
213  y21=grid%X(2,ipoint21)
214  x22=grid%X(1,ipoint22)
215  y22=grid%X(2,ipoint22)
216 C CHECK FOR INTERSECTION BETWEEN TWO INTERVALS
217  IF(gridman_intersect2d(x11,y11,x12,y12,
218  f x21,y21,x22,y22,xi,yi)) THEN
219  res=111
220  WRITE(gridman_unit,*) " GRIDMAN_GRID2D_CHECK: ",
221  w " intersecting edges"
222  WRITE(gridman_unit,*) " ICELL ",icell
223  WRITE(gridman_unit,*) " IEDGE1, IP1, IP2 ",
224  w iedge1,ipoint11,ipoint12
225  WRITE(gridman_unit,*) " X1, X2, Y1, Y2 ",
226  w x11,x12,y11,y12
227  WRITE(gridman_unit,*) " IEDGE, IP1, IP2 ",
228  w iedge,ipoint21,ipoint22
229  WRITE(gridman_unit,*) " X1, X2, Y1, Y2 ",
230  w x21,x22,y21,y22
231  WRITE(gridman_unit,*) " INTERSECTION: XI, YI ",xi,yi
232  GOTO 100
233  END IF
234  END DO !DO IE=IE1+1,CELLS%ILAST(ICELL)
235  END DO !DO IE1=CELLS%IFIRST(ICELL),CELLS%ILAST(ICELL)
236  END DO !DO ICELL=1,GRID%NCELL
237 
238  CALL gridman_indlist_deallocate(cells,ierr0)
239  CALL gridman_indlist_deallocate(points,ierr0)
240 
241  IF(gridman_dbg)
242  f WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CHECK finished"
243 
244  RETURN
245 
246  100 CALL gridman_indlist_deallocate(cells,ierr0)
247  CALL gridman_indlist_deallocate(points,ierr0)
248  RETURN
249 
250  500 CALL gridman_indlist_deallocate(cells,ierr0)
251  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CHECK terminated"
252  RETURN
253 
254  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CHECK terminated"
255 
256  END SUBROUTINE gridman_grid2d_check
257 
258 C*****************************************************************************************
259 C> Find closed chain of points which form each cell
260 C
261 C> In each element of the list first and last index points to the same point.
262 C> IERR=450 is returned of no closed chain is found.
263 C> WARNING: if the object CHAINS already exist it will be overwritten!
264 C*****************************************************************************************
265  SUBROUTINE gridman_grid2d_chains(GRID,CHAINS,IERR)
266  USE gridman,ONLY:gridman_grid,gridman_sp,gridman_unit,
267  u gridman_dbg,gridman_indlist
268  USE gridman_lib,ONLY:gridman_grid_cells,gridman_grid_points,
269  u gridman_indlist_allocate,
270  u gridman_indlist_deallocate
271  IMPLICIT NONE
272 C> Grid object for which the list of chains is generated
273  TYPE(gridman_grid) :: GRID
274 C> Resulting list of chains
275  TYPE(gridman_indlist) :: CHAINS
276 C> Error code
277  INTEGER,INTENT(OUT) :: IERR
278 
279  INTEGER :: IERR0
280  INTEGER(GRIDMAN_SP) :: ICELL,IEDGE0,IEDGE,I,IE,ITMP,
281  i ipoint,ipoint1,ipoint2,n1,n2
282  TYPE(gridman_indlist) :: CELLS,POINTS
283 
284  IF(gridman_dbg)
285  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_CHAINS"
286 
287  ierr=0
288 
289  IF(grid%NCELLS.LT.1) THEN
290  chains%N=0
291  chains%L=0
292  RETURN
293  END IF
294 
295  IF(grid%EDIM.NE.2) THEN
296  ierr=100
297  WRITE(gridman_unit,*) "Error in GRIDMAN_GRID2D_CHAINS: "
298  WRITE(gridman_unit,*) "This subroutine is applicable only if ",
299  w "edges are defined by exactly two points: EDIM=2"
300  WRITE(gridman_unit,*) " EDIM ",grid%EDIM
301  RETURN
302  END IF
303 
304  CALL gridman_grid_cells(cells,grid,ierr)
305  IF(ierr.NE.0) GOTO 1000
306  CALL gridman_grid_points(points,grid,ierr)
307  IF(ierr.NE.0) GOTO 1000
308  CALL gridman_indlist_allocate(chains,cells%N,cells%L+cells%N,ierr)
309  IF(ierr.NE.0) GOTO 1000
310 
311  chains%IFIRST(1)=1
312  DO icell=1,grid%NCELLS
313 C TAKE FIRST EDGE OF THE CELL
314  iedge0=cells%IND(cells%IFIRST(icell)) !PREVIOUS EDGE
315  IF(iedge0.LT.1.OR.iedge0.GT.grid%NEDGES) THEN
316  WRITE(gridman_unit,*)
317  w "ERROR in GRIDMAN_GRID2D_CHAINS:",
318  w " internal error - wrong edge index"
319  WRITE(gridman_unit,*) " ICELL, IEDGE0 ",icell,iedge
320  GOTO 100
321  END IF
322  ipoint1=grid%POINTS(1,iedge0) !FIRST POINT IN THE CHAIN
323  ipoint=grid%POINTS(2,iedge0) !CURRENT POINT IN THE CHAIN
324  chains%ILAST(icell)=chains%IFIRST(icell)
325  chains%IND( chains%ILAST(icell) )=ipoint1
326  chains%ILAST(icell)=chains%ILAST(icell)+1
327  chains%IND( chains%ILAST(icell) )=ipoint
328 C NUMBER OF ITERATIONS IN THE CYCLE IS EQUAL TO THE NUMBER
329 C OF EDGES BELONGING TO THE CELL MINUS ONE
330  DO i=cells%IFIRST(icell)+1,cells%ILAST(icell)
331 C LOOK FOR THE NEXT EDGE IN THE CHAIN: EDGE CONENCTED TO
332 C THE CURRENT POINT IPOINT AND BELONGING TO THE CURRENT CELL ICELL
333  DO ie=points%IFIRST(ipoint),points%ILAST(ipoint)
334  iedge=points%IND(ie)
335  IF(iedge.NE.iedge0) THEN !IF NOT PREVIOUS EDGE
336  IF(iedge.LT.1.OR.iedge.GT.grid%NEDGES) THEN
337  WRITE(gridman_unit,*)
338  w "ERROR in GRIDMAN_GRID2D_CHAINS:",
339  w " internal error - wrong edge index"
340  WRITE(gridman_unit,*) " ICELL, IEDGE ",icell,iedge
341  GOTO 100
342  END IF
343  IF(grid%CELLS(1,iedge).EQ.icell.OR.
344  f grid%CELLS(2,iedge).EQ.icell) THEN
345 C NEXT EDGE IS FOUND
346  ipoint1=grid%POINTS(1,iedge) !FIRST POINT IN THE CHAIN
347  ipoint2=grid%POINTS(2,iedge) !CURRENT POINT IN THE CHAIN
348  IF(ipoint1.NE.ipoint) THEN
349  itmp=ipoint2
350  ipoint2=ipoint1
351  ipoint1=itmp
352  IF(ipoint1.NE.ipoint) THEN
353  ierr=400
354  WRITE(gridman_unit,*)
355  w "ERROR in GRIDMAN_GRID2D_CHAINS:",
356  w " internal error - the found edge is wrong"
357  WRITE(gridman_unit,*)
358  w " ICELL, IEDGE, IPOINT, IPOINT1, IPOINT2 ",
359  w icell, iedge, ipoint, ipoint1, ipoint2
360  GOTO 100
361  END IF
362  END IF
363 C NOW IPOINT1=IPOINT ------------ IPOINT2
364 C ADD POINT2 TO THE CHAIN
365  chains%ILAST(icell)=chains%ILAST(icell)+1
366  chains%IND( chains%ILAST(icell) )=ipoint2
367 C SHIFT IPOINT AND IEDGE0 TO THE NEXT POSITION
368  ipoint=ipoint2
369  iedge0=iedge
370  EXIT
371  END IF
372  END IF !IF(IEDGE.NE.IEDGE0)
373  END DO !DO IE=POINTS%IFIRST(IPOINT)+1,POINTS%ILAST(IPOINT)
374  IF(icell.LT.grid%NCELLS)
375  f chains%IFIRST(icell+1)=chains%ILAST(icell)+1
376  END DO !DO I=CELLS%IFIRST(ICELL)+1,CELLS%ILAST(ICELL)
377 C FIRST AND LAST POINT OF THE CHAIN MUST BE SAME
378  ipoint1=chains%IND(chains%IFIRST(icell))
379  ipoint2=chains%IND(chains%ILAST(icell))
380  IF(ipoint1.NE.ipoint2) THEN
381  ierr=450
382  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CHAINS:" ,
383  w " internal error - no closed chain"
384  WRITE(gridman_unit,*) " ICELL, IPOINT, IPOINT0 ",
385  w icell, ipoint2, ipoint1
386  GOTO 100
387  END IF
388 C CHECK
389  n1=chains%ILAST(icell)-chains%IFIRST(icell)
390  n2=cells%ILAST(icell)-cells%IFIRST(icell)
391  IF(n1.NE.n2+1) THEN
392  ierr=400
393  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CHAINS:" ,
394  w " internal error - mismatch between CHAIN and CELLS"
395  WRITE(gridman_unit,*) " ICELL ",icell
396  WRITE(gridman_unit,*) " CHAIN, ILAST-IFIRST ", n1
397  WRITE(gridman_unit,*) " CELLS, ILAST-IFIRST+1",n2
398  GOTO 100
399  END IF
400  END DO !DO ICELL=1,GRID%NCELLS
401 
402  CALL gridman_indlist_deallocate(cells,ierr0)
403  CALL gridman_indlist_deallocate(points,ierr0)
404 
405  IF(gridman_dbg)
406  f WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CHAINS finished"
407  RETURN
408 
409  100 CALL gridman_indlist_deallocate(chains,ierr0)
410  CALL gridman_indlist_deallocate(cells,ierr0)
411  CALL gridman_indlist_deallocate(points,ierr0)
412  RETURN
413 
414  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CHAINS terminated"
415  CALL gridman_indlist_deallocate(chains,ierr0)
416  CALL gridman_indlist_deallocate(cells,ierr0)
417  CALL gridman_indlist_deallocate(points,ierr0)
418 
419  END SUBROUTINE gridman_grid2d_chains
420 
421 *****************************************************************************************
422 C> Calculate lengths of the cell edges
423 C*****************************************************************************************
424  SUBROUTINE gridman_grid2d_lengths(GRID,LEDGES,IERR)
425  USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
426  u gridman_dbg,gridman_check
427  USE gridman_lib,ONLY:gridman_grid2d_check
428  IMPLICIT NONE
429  INTRINSIC sqrt
430 C> Grid object for which the areas are calculated
431  TYPE(gridman_grid) :: GRID
432 C> Resulting array of lengths (has to be allocated)
433  REAL(GRIDMAN_DP) :: LEDGES(grid%nedges)
434 C> Error code
435  INTEGER,INTENT(OUT) :: IERR
436 
437  INTEGER(GRIDMAN_SP) :: IEDGE,IPOINT
438  INTEGER :: RES
439  REAL(GRIDMAN_DP) :: X1,X2,Y1,Y2
440 
441  IF(gridman_dbg)
442  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_LENGTHS"
443 
444  IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
445  ierr=100
446  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_LENGTHS: ",
447  w "the grid is not of 2D type"
448  WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
449  w grid%TYPE,grid%PDIM,grid%EDIM
450  RETURN
451  END IF
452 
453  IF(gridman_check) THEN
454  CALL gridman_grid2d_check(grid,res,ierr)
455  IF(res.NE.0.OR.ierr.GT.0) THEN
456  ierr=100
457  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_LENGTHS: ",
458  w "the grid is incorrect"
459  RETURN
460  END IF
461  END IF
462 
463  DO iedge=1,grid%NEDGES
464  ipoint=grid%POINTS(1,iedge)
465  x1=grid%X(1,ipoint)
466  y1=grid%X(2,ipoint)
467  ipoint=grid%POINTS(2,iedge)
468  x2=grid%X(1,ipoint)
469  y2=grid%X(2,ipoint)
470  ledges(iedge)=sqrt((x1-x2)**2+(y1-y2)**2)
471  END DO
472 
473  IF(gridman_dbg)
474  f WRITE(gridman_unit,*) "GRIDMAN_GRID2D_LENGTHS completed"
475 
476  END SUBROUTINE gridman_grid2d_lengths
477 
478 
479 C*****************************************************************************************
480 C> Calculate cross section area of the cells
481 C*****************************************************************************************
482  SUBROUTINE gridman_grid2d_crossect(GRID,SCELLS,IERR)
483  USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
484  u gridman_dbg,gridman_check,gridman_indlist
485  USE gridman_lib,ONLY:gridman_grid2d_check,gridman_grid2d_chains,
486  u gridman_indlist_deallocate
487  IMPLICIT NONE
488  INTRINSIC abs
489 C> Grid for which the areas are calculated
490  TYPE(gridman_grid) :: GRID
491 C> Resulting arrray with cell areas (has to be allocated)
492  REAL(GRIDMAN_DP) :: SCELLS(grid%ncells)
493 C> Error code
494  INTEGER,INTENT(OUT) :: IERR
495 
496  INTEGER(GRIDMAN_SP) :: ICELL,IP,IPOINT1,IPOINT2
497  INTEGER :: RES,IERR0
498  REAL(GRIDMAN_DP) :: X1,X2,Y1,Y2,S
499  TYPE(gridman_indlist) :: CHAINS
500 
501  IF(gridman_dbg)
502  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_CROSSECT"
503 
504  ierr=0
505 
506  IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
507  ierr=100
508  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CROSSECT: ",
509  w "the grid is not of 2D type"
510  WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
511  w grid%TYPE,grid%PDIM,grid%EDIM
512  RETURN
513  END IF
514 
515  IF(gridman_check) THEN
516  CALL gridman_grid2d_check(grid,res,ierr)
517  IF(res.NE.0.OR.ierr.GT.0) THEN
518  ierr=100
519  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CROSSECT: ",
520  w "the grid is incorrect"
521  RETURN
522  END IF
523  END IF
524 
525  CALL gridman_grid2d_chains(grid,chains,ierr)
526  IF(ierr.NE.0) GOTO 1000
527 
528  DO icell=1,grid%NCELLS
529  s=0.
530  DO ip=chains%IFIRST(icell),chains%ILAST(icell)-1
531  ipoint1=chains%IND(ip)
532  ipoint2=chains%IND(ip+1)
533  x1=grid%X(1,ipoint1)
534  y1=grid%X(2,ipoint1)
535  x2=grid%X(1,ipoint2)
536  y2=grid%X(2,ipoint2)
537  s=s+x1*y2-x2*y1
538  END DO
539  scells(icell)=0.5*abs(s)
540  END DO
541 
542  CALL gridman_indlist_deallocate(chains,ierr)
543 
544  IF(gridman_dbg)
545  f WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CROSSECT finished"
546  RETURN
547 
548  1000 CALL gridman_indlist_deallocate(chains,ierr0)
549  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CROSSECT terminated"
550 
551  END SUBROUTINE gridman_grid2d_crossect
552 
553 C*****************************************************************************************
554 C> Calculate coordinates of the cell centers
555 C>
556 C> Coordinates of the center of mass of planar polygons
557 C*****************************************************************************************
558  SUBROUTINE gridman_grid2d_center(GRID,XCN,IERR)
559  USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
560  u gridman_dbg,gridman_check,gridman_indlist,
561  u gridman_tol
562  USE gridman_lib,ONLY:gridman_grid2d_check,gridman_grid2d_chains,
563  u gridman_indlist_deallocate
564  IMPLICIT NONE
565  INTRINSIC abs,tiny
566 C> Grid object for which the coordinates are calculated
567  TYPE(gridman_grid) :: GRID
568 C> Resulting X (1st index=1) and Y (1st index=2) coordinates of the centers.
569 C> The array must be allocated
570  REAL(GRIDMAN_DP) :: XCN(2,grid%ncells)
571 C> Error code
572  INTEGER,INTENT(OUT) :: IERR
573 
574  INTEGER(GRIDMAN_SP) :: ICELL,IP,IPOINT,IPOINT1,IPOINT2,N
575  INTEGER :: RES,IERR0
576  REAL(GRIDMAN_DP) :: X1,X2,Y1,Y2,S,S2,XC,YC,XC2,YC2
577  TYPE(gridman_indlist) :: CHAINS
578 
579  IF(gridman_dbg)
580  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_CENTER"
581 
582  ierr=0
583 
584  IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
585  ierr=100
586  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CENTER: ",
587  w "the grid is not of 2D type"
588  WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
589  w grid%TYPE,grid%PDIM,grid%EDIM
590  RETURN
591  END IF
592 
593  IF(gridman_check) THEN
594  CALL gridman_grid2d_check(grid,res,ierr)
595  IF(res.NE.0.OR.ierr.GT.0) THEN
596  ierr=100
597  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CENTER: ",
598  w "the grid is incorrect"
599  RETURN
600  END IF
601  END IF
602 
603  CALL gridman_grid2d_chains(grid,chains,ierr)
604  IF(ierr.NE.0) GOTO 1000
605 
606  DO icell=1,grid%NCELLS
607  s=0.
608  xc=0.
609  yc=0.
610 C S2, XC2, YC2 ARE CALCULATED ONLY FOR DIAGNOSTIC
611  s2=0.
612  yc2=0.
613  xc2=0.
614  DO ip=chains%IFIRST(icell),chains%ILAST(icell)-1
615  ipoint1=chains%IND(ip)
616  ipoint2=chains%IND(ip+1)
617  x1=grid%X(1,ipoint1)
618  y1=grid%X(2,ipoint1)
619  x2=grid%X(1,ipoint2)
620  y2=grid%X(2,ipoint2)
621 C UPDATE THE SUMS
622  s=s+(x2+x1)*(y2-y1)
623  s2=s2+(x1*y2-x2*y1)
624  xc=xc+(x2+x1)*(x1*y2-x2*y1)
625  xc2=xc2+(x2**2+x1**2+x2*x1)*(y2-y1)
626  yc=yc+(y2+y1)*(x1*y2-x2*y1)
627  yc2=yc2+(y2**2+y1**2+y2*y1)*(x1-x2)
628  END DO
629 C CHECKS
630  IF(abs(s-s2).GT.gridman_tol*(abs(s)+abs(s2))) THEN
631  ierr=400
632  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CENTER: ",
633  w "the area calculation is incorrect"
634  WRITE(gridman_unit,*) " ICELL, S1, S2 ", icell, s, s2
635  GOTO 100
636  END IF
637  IF(abs(xc-xc2).GT.gridman_tol*(abs(xc)+abs(xc2))) THEN
638  ierr=400
639  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CENTER: ",
640  w "the center coordinate calculation is incorrect"
641  WRITE(gridman_unit,*) " ICELL, XC1, XC2 ",icell,xc,xc2
642  GOTO 100
643  END IF
644  IF(abs(yc-yc2).GT.gridman_tol*(abs(yc)+abs(yc2))) THEN
645  ierr=400
646  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CENTER: ",
647  w "the center coordinate calculation is incorrect"
648  WRITE(gridman_unit,*) " ICELL, YC1, YC2 ",icell,yc,yc2
649  GOTO 100
650  END IF
651  IF(abs(s).GT.tiny(s)) THEN
652 C UPDATE COORDINATES
653  xcn(1,icell)=xc/(3.*s)
654  xcn(2,icell)=yc/(3.*s)
655  ELSE
656 C SPECIAL CASE OF ZERO AREA
657  IF(gridman_dbg) WRITE(gridman_unit,*)
658  w "WARNING from GRIDMAN_GRID2D_CENTER: ",
659  w "cell has zero area, ICELL ", icell
660  xc=0
661  yc=0
662  DO ip=chains%IFIRST(icell),chains%ILAST(icell)-1
663  ipoint=chains%IND(ip)
664  xc=xc+grid%X(1,ipoint)
665  yc=yc+grid%X(2,ipoint)
666  END DO
667  n=chains%ILAST(icell)-chains%IFIRST(icell)
668  IF(n.GT.0) THEN
669  xcn(1,icell)=xc/REAL(n,gridman_dp)
670  xcn(2,icell)=yc/REAL(n,gridman_dp)
671  END IF
672  END IF
673  END DO !DO ICELL=1,GRID%NCELLS
674 
675  CALL gridman_indlist_deallocate(chains,ierr)
676 
677  IF(gridman_dbg)
678  f WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CENTER finished"
679  RETURN
680 
681  100 CALL gridman_indlist_deallocate(chains,ierr0)
682  RETURN
683 
684  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CENTER terminated"
685 
686  END SUBROUTINE gridman_grid2d_center
687 
688 C*****************************************************************************************
689 C> Calculate unit normal vectors to grid edges
690 C*****************************************************************************************
691  SUBROUTINE gridman_grid2d_norm(GRID,VN,IERR)
692  USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
693  u gridman_dbg,gridman_check,gridman_indlist,
694  u gridman_tol
695  USE gridman_lib,ONLY:gridman_grid2d_check,gridman_grid2d_chains,
696  u gridman_indlist_deallocate
697  IMPLICIT NONE
698  INTRINSIC abs,sign,tiny
699 C> Grid for which the normals are calculated
700  TYPE(gridman_grid) :: GRID
701 C> Resulting array with X (1st index=1) and Y (1st index=2) coordinates
702 C> of the normals (array has to be allocated)
703  REAL(GRIDMAN_DP) :: VN(2,grid%nedges)
704 C> Error code
705  INTEGER,INTENT(OUT) :: IERR
706 
707  INTEGER(GRIDMAN_SP) :: IEDGE,ICELL,IP,IPOINT1,IPOINT2,
708  i ipoint01,ipoint02
709  INTEGER :: IERR0,RES
710  REAL(GRIDMAN_DP) :: XN,YN,NORM,PROD,X1,Y1,X2,Y2,S
711  TYPE(gridman_indlist) :: CHAINS
712 
713  IF(gridman_dbg)
714  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_NORM"
715 
716  ierr=0
717 
718  IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
719  ierr=100
720  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_NORM: ",
721  w "the grid is not of 2D type"
722  WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
723  w grid%TYPE,grid%PDIM,grid%EDIM
724  RETURN
725  END IF
726 
727  IF(gridman_check) THEN
728  CALL gridman_grid2d_check(grid,res,ierr)
729  IF(res.NE.0.OR.ierr.GT.0) THEN
730  ierr=100
731  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_NORM: ",
732  w "the grid is incorrect"
733  RETURN
734  END IF
735  END IF
736 
737  CALL gridman_grid2d_chains(grid,chains,ierr)
738  IF(ierr.NE.0) GOTO 1000
739 
740  DO iedge=1,grid%NEDGES ! CYCLE OVER EDGES
741 C COORDINATES OF THE VECTOR WHICH POINTS FROM X1 TO X2
742  ipoint01=grid%POINTS(1,iedge)
743  ipoint02=grid%POINTS(2,iedge)
744  x1=grid%X(1,ipoint01)
745  y1=grid%X(2,ipoint01)
746  x2=grid%X(1,ipoint02)
747  y2=grid%X(2,ipoint02)
748 C COORDINATES OF THE NORMAL VECTOR: ROTATED
749 C CONTER-CLOCKWISE RELATIVE TO (X1,X2) X2 ^
750  xn=y1-y2 ! |
751  yn=x2-x1 ! <--- X1 |
752 C NOW CORRECT THE DIRECTION OF NORMAL, IF NECESSARY.
753 C> Convention about the direction of normal vectors
754 C> 1. the normal must point out of the cell CELLS(1,IEDGE)
755 C> 2. If CELLS(1,IEDGE)<0 then the normal must point out of the cell CELLS(2,IEDGE)
756 C> 3. IF both CELLS(:,IEDGE)==0 then the normal must point to the left
757 C> from the direction POINT(1,1,IEDGE) to POINT(1,2,IEDGE)
758 C>
759 C> The normal vectors, thus, always point outwards from the grid
760  icell=grid%CELLS(1,iedge)
761  IF(icell.LT.1) THEN
762  icell=grid%CELLS(2,iedge)
763  IF(icell.LT.1) THEN
764 C BOTH CELLS ARE ZERO
765  vn(1,iedge)=xn
766  vn(2,iedge)=yn
767  cycle
768  END IF
769 C SECOND CELL IS "NON-ZERO"
770  END IF
771 C SIGNED AREA OF THE CELL
772  s=0.
773  DO ip=chains%IFIRST(icell),chains%ILAST(icell)-1
774  ipoint1=chains%IND(ip)
775  ipoint2=chains%IND(ip+1)
776  x1=grid%X(1,ipoint1)
777  y1=grid%X(2,ipoint1)
778  x2=grid%X(1,ipoint2)
779  y2=grid%X(2,ipoint2)
780  s=s+x1*y2-x2*y1
781 C REVERT NORMAL IF DIRECTION OF IEDGE DIFFERS
782 C FROM THE DIRECTION OF CHAIN
783  IF(ipoint01.EQ.ipoint2.AND.ipoint02.EQ.ipoint1) THEN
784  xn=-xn
785  yn=-yn
786  END IF
787  END DO
788 C NORMALIZE
789  norm=sqrt(xn**2+yn**2)
790  IF(abs(norm).LE.10.*tiny(norm)) norm=1.0
791 C CHANGE SIGN IF SIGNED AREA S > 0
792  vn(1,iedge)=-xn*sign(1._gridman_dp,s)/norm
793  vn(2,iedge)=-yn*sign(1._gridman_dp,s)/norm
794 C CHECK
795  ipoint1=grid%POINTS(1,iedge)
796  ipoint2=grid%POINTS(2,iedge)
797  x1=grid%X(1,ipoint1)
798  y1=grid%X(2,ipoint1)
799  x2=grid%X(1,ipoint2)
800  y2=grid%X(2,ipoint2)
801 C SCALAR PRODUCT OF THE NORMAL AND EDGE VECTORS MUST BE ZERO
802  prod=(x2-x1)*vn(1,iedge)+(y2-y1)*vn(2,iedge)
803  IF(abs(prod).GT.gridman_tol*norm) THEN
804  ierr=400
805  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_NORM: ",
806  w "the calculation of normal vector is incorrect"
807  WRITE(gridman_unit,*) "IEDGE, S, X2-X1, Y2-Y1, XN, YN, ",
808  w iedge,s,x2-x1,y2-y1,vn(1,iedge),vn(2,iedge)
809  GOTO 100
810  END IF
811  END DO !DO IEDGE=1,GRID%NEDGES
812 
813  CALL gridman_indlist_deallocate(chains,ierr)
814 
815  IF(gridman_dbg)
816  w WRITE(gridman_unit,*) "GRIDMAN_GRID2D_NORM finished"
817  RETURN
818 
819  100 CALL gridman_indlist_deallocate(chains,ierr0)
820  RETURN
821 
822  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID2D_NORM terminated"
823 
824  END SUBROUTINE gridman_grid2d_norm
825 
826 C*****************************************************************************************
827 C> Find if cells are convex polygons or not
828 C*****************************************************************************************
829  SUBROUTINE gridman_grid2d_isconvex(GRID,ISCONVEX,IERR)
830  USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
831  u gridman_dbg,gridman_check,gridman_indlist
832  USE gridman_lib,ONLY:gridman_grid2d_check,gridman_grid2d_chains,
833  u gridman_indlist_deallocate
834  IMPLICIT NONE
835  INTRINSIC sign,abs,tiny,int
836 C> Grid which will be processed
837  TYPE(gridman_grid) :: GRID
838 C> Resulting array where .TRUE. stand for convex polygons and .FALS. otherwise
839 C> (has to be allocated)
840  LOGICAL :: ISCONVEX(grid%ncells)
841 C> Erorr code
842  INTEGER,INTENT(OUT) :: IERR
843 
844  INTEGER(GRIDMAN_SP) :: ICELL,IP,IPOINT1,IPOINT2,N,N0,NZ
845  INTEGER :: IERR0,RES
846  REAL(GRIDMAN_DP) :: X1,X2,Y1,Y2,VX,VY,VX0,VY0,V
847  TYPE(gridman_indlist) :: CHAINS
848 
849  IF(gridman_dbg)
850  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_ISCONVEX"
851 
852  ierr=0
853 
854  IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
855  ierr=100
856  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_ISCONVEX: ",
857  w "the grid is not of 2D type"
858  WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
859  w grid%TYPE,grid%PDIM,grid%EDIM
860  RETURN
861  END IF
862 
863  IF(gridman_check) THEN
864  CALL gridman_grid2d_check(grid,res,ierr)
865  IF(res.NE.0.OR.ierr.GT.0) THEN
866  ierr=100
867  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_ISCONVEX: ",
868  w "the grid is incorrect"
869  RETURN
870  END IF
871  END IF
872 
873  CALL gridman_grid2d_chains(grid,chains,ierr)
874  IF(ierr.NE.0) GOTO 1000
875 C
876 C> The algorithm assumes closed polygons without self intersections.
877 C> Such polygon is convex if all angles between two neighbour edges < 180 degrees
878 C
879  DO icell=1,grid%NCELLS
880  n=0
881  nz=0
882  ip=chains%ILAST(icell)-1
883  ipoint1=chains%IND(ip)
884  ip=chains%ILAST(icell)
885  ipoint2=chains%IND(ip)
886  x1=grid%X(1,ipoint1)
887  y1=grid%X(2,ipoint1)
888  x2=grid%X(1,ipoint2)
889  y2=grid%X(2,ipoint2)
890  vx0=x2-x1
891  vy0=y2-y1
892  DO ip=chains%IFIRST(icell),chains%ILAST(icell)-1
893  ipoint1=chains%IND(ip)
894  ipoint2=chains%IND(ip+1)
895  x1=grid%X(1,ipoint1)
896  y1=grid%X(2,ipoint1)
897  x2=grid%X(1,ipoint2)
898  y2=grid%X(2,ipoint2)
899  vx=x2-x1
900  vy=y2-y1
901 C THE COMBINATION VX0*VY-VX*VY0 HAS TO HAVE SAME SIGN FOR EACH PAIR OF EDGES (CORNER)
902  v=vx0*vy-vx*vy0
903  IF(abs(v).LT.10.*tiny(v)) THEN
904  nz=nz+1 !THE ANGLE IS EXACTLY 180 DEG
905  ELSE
906  n=n+int(sign(1._gridman_dp,v),gridman_sp)
907  END IF
908  vx0=vx
909  vy0=vy
910  END DO !DO IP=CHAINS%IFIRST(ICELL),CHAINS%ILAST(ICELL)-1
911 C TOTAL NUMBER OF CORNERS
912  n0=chains%ILAST(icell)-chains%IFIRST(icell)
913  n=abs(n)+nz
914  IF(abs(n).LT.n0) THEN
915 C SOME CORNERS ARE LARGER THAN 180 GRAD => THE POLYGON IS NOT CONVEX
916  isconvex(icell)=.false.
917  ELSE
918 C ALL CORNERS <= 180 GRAD => THE POLYGON IS CONVEX
919  isconvex(icell)=.true.
920  END IF
921  END DO !DO ICELL=1,GRID%NCELLS
922 
923  CALL gridman_indlist_deallocate(chains,ierr)
924 
925  IF(gridman_dbg)
926  f WRITE(gridman_unit,*) "GRIDMAN_GRID2D_ISCONVEX finished"
927  RETURN
928 
929  1000 CALL gridman_indlist_deallocate(chains,ierr0)
930  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_ISCONVEX terminated"
931 
932  END SUBROUTINE gridman_grid2d_isconvex
933 
934 *****************************************************************************************
935 C> Calculate cylindrical areas of the cell edges
936 C*****************************************************************************************
937 C>
938 C> Calculate areas of the cell edges asumming that the grid
939 C> defines cylindrical coordinates, with X(1) being the radius -
940 C> - distance from the axis of rotation. If the optional
941 C> argument ANGLE is not defined then the azimuthal angle
942 C> occupied by the edges is assumed to be 2*pi.
943  SUBROUTINE gridman_grid2d_cylareas(GRID,SEDGES,IERR,ANGLE)
944  USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
945  u gridman_dbg,gridman_check,gridman_pi
946  USE gridman_lib,ONLY:gridman_grid2d_check
947  IMPLICIT NONE
948  INTRINSIC sqrt,PRESENT
949 C> Grid object for which the areas are calculated
950  TYPE(gridman_grid) :: GRID
951 C> Resulting array of surface areas (has to be allocated)
952  REAL(GRIDMAN_DP) :: SEDGES(grid%nedges)
953 C> Error code
954  INTEGER,INTENT(OUT) :: IERR
955 C> Azimuthal angle occupied by the edges, in radian, default 2*pi
956  REAL(GRIDMAN_DP),OPTIONAL :: ANGLE
957 
958  INTEGER(GRIDMAN_SP) :: IEDGE,IPOINT
959  INTEGER :: RES
960  REAL(GRIDMAN_DP) :: X1,X2,Y1,Y2,FANG
961 
962  IF(gridman_dbg)
963  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_CYLAREAS"
964 
965  ierr=0
966 
967  IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
968  ierr=100
969  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CYLAREAS: ",
970  w "the grid is not of 2D type"
971  WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
972  w grid%TYPE,grid%PDIM,grid%EDIM
973  RETURN
974  END IF
975 
976  IF(gridman_check) THEN
977  CALL gridman_grid2d_check(grid,res,ierr)
978  IF(res.NE.0.OR.ierr.GT.0) THEN
979  ierr=100
980  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CYLAREAS: ",
981  w "the grid is incorrect"
982  RETURN
983  END IF
984  END IF
985 
986  IF(PRESENT(angle)) THEN
987  fang=0.5*angle
988  ELSE
989  fang=gridman_pi
990  END IF
991 
992  DO iedge=1,grid%NEDGES !CYCLE OVER EDGES
993  ipoint=grid%POINTS(1,iedge)
994  x1=grid%X(1,ipoint)
995  y1=grid%X(2,ipoint)
996  ipoint=grid%POINTS(2,iedge)
997  x2=grid%X(1,ipoint)
998  y2=grid%X(2,ipoint)
999 C SIDE AREA OF A CONICAL FRUSTUM
1000  sedges(iedge)=fang*(x1+x2)*sqrt((x1-x2)**2+(y1-y2)**2)
1001  END DO
1002 
1003  IF(gridman_dbg)
1004  f WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CYLAREAS finished"
1005 
1006  END SUBROUTINE gridman_grid2d_cylareas
1007 
1008 C*****************************************************************************************
1009 C> Calculate cylindrical cell volumes
1010 C*****************************************************************************************
1011 C>
1012 C> Calculate volumes of the cells asumming that the grid
1013 C> defines cylindrical coordinates, with X(1) being the radius -
1014 C> - distance from the axis of rotation. If the optional
1015 C> argument ANGLE is not defined then the azimuthal angle
1016 C> occupied by the cells is assumed to be 2*pi.
1017  SUBROUTINE gridman_grid2d_cylvolumes(GRID,VCELLS,IERR,ANGLE)
1018  USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
1019  u gridman_pi,gridman_dbg,gridman_check,gridman_indlist,
1020  u gridman_tol
1021  USE gridman_lib,ONLY:gridman_grid2d_check,gridman_grid2d_chains
1022  IMPLICIT NONE
1023  INTRINSIC abs,PRESENT
1024 C> Grid object for which the volumes are calculated
1025  TYPE(gridman_grid) :: GRID
1026 C> Resulting array of cell volumes (has to be allocated)
1027  REAL(GRIDMAN_DP) :: VCELLS(grid%ncells)
1028 C> Error code
1029  INTEGER,INTENT(OUT) :: IERR
1030 C> Azimuthal angle occupied by the cells, in radian, default 2*pi
1031  REAL(GRIDMAN_DP),OPTIONAL :: ANGLE
1032 
1033  INTEGER(GRIDMAN_SP) :: ICELL,IP,IPOINT1,IPOINT2
1034  INTEGER :: RES,IERR0
1035  REAL(GRIDMAN_DP) :: X1,X2,Y1,Y2,VOL,VOL2,FTOR
1036  TYPE(gridman_indlist) :: CHAINS
1037 C VOL2 IS NEEDED ONLY AS A DYAGNOSTIC
1038 
1039  IF(gridman_dbg)
1040  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_CYLVOLUMES"
1041 
1042  ierr=0
1043 
1044  IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
1045  ierr=100
1046  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CYLVOLUMES: ",
1047  w "the grid is not of 2D type"
1048  WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
1049  w grid%TYPE,grid%PDIM,grid%EDIM
1050  RETURN
1051  END IF
1052 
1053  IF(gridman_check) THEN
1054  CALL gridman_grid2d_check(grid,res,ierr)
1055  IF(res.NE.0.OR.ierr.GT.0) THEN
1056  ierr=100
1057  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CYLVOLUMES: ",
1058  w "the grid is incorrect"
1059  RETURN
1060  END IF
1061  END IF
1062 
1063  CALL gridman_grid2d_chains(grid,chains,ierr)
1064  IF(ierr.NE.0) GOTO 1000
1065 
1066 C FACTOR WHICH TAKES INTO ACCOUNT PRESCIBED TOROIDAL ANGLE OF THE GRID
1067 C DA/2PI*PI/3
1068  IF(PRESENT(angle)) THEN
1069  ftor=angle/6._gridman_dp
1070  ELSE
1071  ftor=gridman_pi/3._gridman_dp
1072  END IF
1073 
1074  DO icell=1,grid%NCELLS
1075  vol=0.
1076  vol2=0.
1077  DO ip=chains%IFIRST(icell),chains%ILAST(icell)-1
1078  ipoint1=chains%IND(ip)
1079  ipoint2=chains%IND(ip+1)
1080  x1=grid%X(1,ipoint1)
1081  y1=grid%X(2,ipoint1)
1082  x2=grid%X(1,ipoint2)
1083  y2=grid%X(2,ipoint2)
1084  vol=vol+(x2+x1)*(x1*y2-x2*y1) !UPDATE THE SUM
1085  vol2=vol2+(x2**2+x1**2+x2*x1)*(y2-y1)
1086  END DO
1087 C CHECK
1088  IF(abs(vol-vol2).GT.gridman_tol*(abs(vol)+abs(vol2))) THEN
1089  ierr=400
1090  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CYLVOLUMES: ",
1091  w "the volume calculation is incorrect"
1092  WRITE(gridman_unit,*)
1093  w " ICELL, VOL1, VOL2, TOL, TOL*(VOL+VOL2)",
1094  w icell,vol,vol2,gridman_tol,gridman_tol*(vol+vol2)
1095  GOTO 100
1096  END IF
1097  vcells(icell)=abs(vol)*ftor
1098  END DO !DO ICELL=1,GRID%NCELLS
1099 
1100  CALL gridman_indlist_deallocate(chains,ierr)
1101 
1102  IF(gridman_dbg)
1103  f WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CYLVOLUMES finished"
1104  RETURN
1105 
1106  100 CALL gridman_indlist_deallocate(chains,ierr0)
1107  RETURN
1108 
1109  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CYLVOLUMES terminated"
1110 
1111  END SUBROUTINE gridman_grid2d_cylvolumes
gridman_grid_points
subroutine gridman_grid_points(EDGES, GRID, IERR)
Create table of edges connected to each point.
Definition: grid2.f:361
gridman_grid_check
subroutine gridman_grid_check(GRID, RES, IERR)
Check consistency of the grid data.
Definition: grid1.f:289
gridman_intersect2d
logical function gridman_intersect2d(X11, Y11, X12, Y12, X21, Y21, X22, Y22, XI, YI)
Intersection of two intervals on a plane.
Definition: geom.f:167
gridman_indlist_allocate
subroutine gridman_indlist_allocate(INDLIST, N, L, IERR)
Allocate list of elements.
Definition: indlist.f:32
gridman::gridman_unit
integer, save, public gridman_unit
Index of the standard output unit.
Definition: gridman.f:116
gridman_grid2d_cylvolumes
subroutine gridman_grid2d_cylvolumes(GRID, VCELLS, IERR, ANGLE)
Calculate cylindrical cell volumes.
Definition: grid2d.f:1018
gridman::gridman_dp
integer, parameter, public gridman_dp
Kind parameter for real numbers.
Definition: gridman.f:93
gridman_lib
Explicit interfaces to GRIDMAN subroutines and functions.
Definition: gridman.f:251
gridman_grid2d_chains
subroutine gridman_grid2d_chains(GRID, CHAINS, IERR)
Find closed chain of points which form each cell.
Definition: grid2d.f:266
gridman_grid_cells
subroutine gridman_grid_cells(EDGES, GRID, IERR)
Create a list of edges which belong to each cell.
Definition: grid2.f:177
gridman_grid2d_check
subroutine gridman_grid2d_check(GRID, RES, IERR)
Check correctness of the 2D grid object.
Definition: grid2d.f:37
gridman_grid2d_cylareas
subroutine gridman_grid2d_cylareas(GRID, SEDGES, IERR, ANGLE)
Calculate cylindrical areas of the cell edges.
Definition: grid2d.f:944
gridman_indlist_deallocate
subroutine gridman_indlist_deallocate(INDLIST, IERR)
Deallocate list of indices.
Definition: indlist.f:106
gridman::gridman_grid
Data-type which describes a grid as a set of edges, methods in grid.f.
Definition: gridman.f:168
gridman_grid2d_center
subroutine gridman_grid2d_center(GRID, XCN, IERR)
Calculate coordinates of the cell centers.
Definition: grid2d.f:559
gridman_grid2d_isconvex
subroutine gridman_grid2d_isconvex(GRID, ISCONVEX, IERR)
Find if cells are convex polygons or not.
Definition: grid2d.f:830
gridman_grid2d_norm
subroutine gridman_grid2d_norm(GRID, VN, IERR)
Calculate unit normal vectors to grid edges.
Definition: grid2d.f:692
gridman::gridman_indlist
Data-type which describes lists of elements with variable number of indices for each element...
Definition: gridman.f:225
gridman
Definition of data types, global constants and variables.
Definition: gridman.f:83
gridman_grid2d_lengths
subroutine gridman_grid2d_lengths(GRID, LEDGES, IERR)
Calculate lengths of the cell edges.
Definition: grid2d.f:425
gridman_grid2d_crossect
subroutine gridman_grid2d_crossect(GRID, SCELLS, IERR)
Calculate cross section area of the cells.
Definition: grid2d.f:483
gridman::gridman_sp
integer, parameter, public gridman_sp
Kind parameter for integer numbers.
Definition: gridman.f:95
Generated on Thu Oct 5 2017 09:49:57 for GRIDMAN by   doxygen 1.8.9.1