GRIDMAN
grid managment library
test_cut.f
Go to the documentation of this file.
1 C> @file tests/test_cut.f
2 C> Unit test for GRIRDMAN_GRID2D_CUT
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  PROGRAM test_cut
24  USE gridman
25  USE gridman_lib
26  IMPLICIT NONE
27 
28  TYPE(gridman_grid) :: grid,cutgrid
29  INTEGER :: ierr,res
30  INTEGER(GRIDMAN_SP),PARAMETER :: np1=9,np2=18,np3=13
31  REAL(GRIDMAN_DP) :: xp(6),yp(6),
32  r xp1(np1),yp1(np1),xp2(np2),yp2(np2),
33  r xp3(np3),yp3(np3)
34  LOGICAL :: lin
35 
36  INTERFACE
37  SUBROUTINE test_mappings(GRID0,GRID1,XP,YP,NP,IERR)
38  USE gridman
39  TYPE(gridman_grid) :: grid0,grid1
40  INTEGER(GRIDMAN_SP),INTENT(IN) :: np
41  REAL(GRIDMAN_DP) :: xp(np),yp(np)
42  INTEGER,INTENT(OUT) :: ierr
43  END SUBROUTINE test_mappings
44  END INTERFACE
45 
46 C TEST GRIDMAN_POINT_IN_POLYGON
47 C POLYGON 1
48  CALL test_polygon1(xp1,yp1,np1)
49  lin=gridman_point_in_polygon(xp1,yp1,np1,
50  . 7._gridman_dp,2._gridman_dp)
51  IF(.NOT.lin) THEN
52  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP1, (7,2)"
53  stop "TEST_CUT TERMINATED"
54  END IF
55  lin=gridman_point_in_polygon(xp1,yp1,np1,
56  . 10._gridman_dp,1.1_gridman_dp)
57  IF(.NOT.lin) THEN
58  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP1, (10,1.1)"
59  stop "TEST_CUT TERMINATED"
60  END IF
61  lin=gridman_point_in_polygon(xp1,yp1,np1,
62  . 6.1_gridman_dp,3._gridman_dp)
63  IF(.NOT.lin) THEN
64  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP1, (6.1,3)"
65  stop "TEST_CUT TERMINATED"
66  END IF
67 
68  lin=gridman_point_in_polygon(xp1,yp1,np1,
69  . 6._gridman_dp,6._gridman_dp)
70  IF(lin) THEN
71  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP1, (6,6)"
72  stop "TEST_CUT TERMINATED"
73  END IF
74  lin=gridman_point_in_polygon(xp1,yp1,np1,
75  . 10._gridman_dp,5._gridman_dp)
76  IF(lin) THEN
77  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP1, (10,5)"
78  stop "TEST_CUT TERMINATED"
79  END IF
80  lin=gridman_point_in_polygon(xp1,yp1,np1,
81  . 7._gridman_dp,-0.51_gridman_dp)
82  IF(lin) THEN
83  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP1, (7,-0.51)"
84  stop "TEST_CUT TERMINATED"
85  END IF
86 
87 C POLYGON 2
88  CALL test_polygon2(xp2,yp2,np2)
89  lin=gridman_point_in_polygon(xp2,yp2,np2,
90  . 15._gridman_dp,2._gridman_dp)
91  IF(.NOT.lin) THEN
92  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (15,2)"
93  stop "TEST_CUT TERMINATED"
94  END IF
95  lin=gridman_point_in_polygon(xp2,yp2,np2,
96  . 9.001_gridman_dp,-5.999_gridman_dp)
97  IF(.NOT.lin) THEN
98  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for ",
99  w "XP2, (9.001,-5.999)"
100  stop "TEST_CUT TERMINATED"
101  END IF
102  lin=gridman_point_in_polygon(xp2,yp2,np2,
103  . 18._gridman_dp,-1._gridman_dp)
104  IF(.NOT.lin) THEN
105  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (18,-1)"
106  stop "TEST_CUT TERMINATED"
107  END IF
108  lin=gridman_point_in_polygon(xp2,yp2,np2,
109  . 17._gridman_dp,7._gridman_dp)
110  IF(.NOT.lin) THEN
111  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (17,7)"
112  stop "TEST_CUT TERMINATED"
113  END IF
114  lin=gridman_point_in_polygon(xp2,yp2,np2,
115  . 16._gridman_dp,3.5_gridman_dp)
116  IF(.NOT.lin) THEN
117  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (16,3.5)"
118  stop "TEST_CUT TERMINATED"
119  END IF
120 
121  lin=gridman_point_in_polygon(xp2,yp2,np2,
122  . 15.5_gridman_dp,0.5_gridman_dp)
123  IF(lin) THEN
124  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (15.5,0.5)"
125  stop "TEST_CUT TERMINATED"
126  END IF
127  lin=gridman_point_in_polygon(xp2,yp2,np2,
128  . 19.0001_gridman_dp,-1._gridman_dp)
129  IF(lin) THEN
130  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for ",
131  w "XP2 (19.0001,-1)"
132  stop "TEST_CUT TERMINATED"
133  END IF
134  lin=gridman_point_in_polygon(xp2,yp2,np2,
135  . 12._gridman_dp,5.5_gridman_dp)
136  IF(lin) THEN
137  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (12,5.5)"
138  stop "TEST_CUT TERMINATED"
139  END IF
140  lin=gridman_point_in_polygon(xp2,yp2,np2,
141  . 14._gridman_dp,3._gridman_dp)
142  IF(lin) THEN
143  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (14,3)"
144  stop "TEST_CUT TERMINATED"
145  END IF
146  lin=gridman_point_in_polygon(xp2,yp2,np2,
147  . 14._gridman_dp,9.01_gridman_dp)
148  IF(lin) THEN
149  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (14,9.01)"
150  stop "TEST_CUT TERMINATED"
151  END IF
152  lin=gridman_point_in_polygon(xp2,yp2,np2,
153  . 14._gridman_dp,-6.51_gridman_dp)
154  IF(lin) THEN
155  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (14,-6.51)"
156  stop "TEST_CUT TERMINATED"
157  END IF
158  lin=gridman_point_in_polygon(xp2,yp2,np2,
159  . 12._gridman_dp,2._gridman_dp)
160  IF(lin) THEN
161  WRITE(*,*) "GRIDMAN_POINT_IN_POLYGON failed for XP2, (12,2)"
162  stop "TEST_CUT TERMINATED"
163  END IF
164 
165  WRITE(*,*) "test_cut: test 0 passed"
166 
167 C TEST GRIRDMAN_GRID2D_CUT
168  CALL grid_example1(grid,ierr)
169  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
170 
171 C TAKE WHOLE GRID
172  xp(1)=1.
173  yp(1)=-5.
174  xp(2)=1.
175  yp(2)=5.
176  xp(3)=13.
177  yp(3)=5.
178  xp(4)=13.
179  yp(4)=-5.
180  xp(5)=xp(1)
181  yp(5)=yp(1)
182  CALL gridman_grid2d_cut(cutgrid,grid,xp,yp,
183  c 5_gridman_sp,.false.,ierr)
184  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
185 
186  CALL gridman_grid_compare(grid,cutgrid,res,ierr)
187  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
188  IF(res.NE.7) THEN
189  WRITE(*,*) "Grids are different"
190  stop "TEST_GRID TERMINATED"
191  END IF
192 
193  CALL test_mappings(grid,cutgrid,xp,yp,5_gridman_sp,ierr)
194  IF(ierr.NE.0) THEN
195  WRITE(*,*) "Wrong indices in test 1.1"
196  stop "TEST_GRID TERMINATED"
197  END IF
198 
199  WRITE(*,*) "test_cut: test 1.1 passed"
200 
201 C TAKE SINGLE POINT
202  xp(1)=6.
203  yp(1)=1.
204  xp(2)=6.
205  yp(2)=3.
206  xp(3)=8.
207  yp(3)=3.
208  xp(4)=8.
209  yp(4)=1.
210  xp(5)=xp(1)
211  yp(5)=yp(1)
212  CALL gridman_grid2d_cut(cutgrid,grid,xp,yp,
213  c 5_gridman_sp,.false.,ierr)
214  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
215 
216  CALL gridman_grid_write(cutgrid,'cutA.grd',ierr)
217  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
218 
219  IF(cutgrid%NCELLS.NE.4.OR.cutgrid%NEDGES.NE.12.OR.
220  f cutgrid%NPOINTS.NE.9.OR.
221  f cutgrid%CELLINDEX(1)%NELEMENTS.NE.4.OR.
222  f cutgrid%EDGEINDEX(1)%NELEMENTS.NE.4) THEN
223  WRITE(*,*) "Unexpected resulting grid (0)"
224  stop "TEST_GRID TERMINATED"
225  END IF
226 
227  CALL test_mappings(grid,cutgrid,xp,yp,5_gridman_sp,ierr)
228  IF(ierr.NE.0) THEN
229  WRITE(*,*) "Wrong indices in test 1.2"
230  stop "TEST_GRID TERMINATED"
231  END IF
232 
233  WRITE(*,*) "test_cut: test 1.2 passed"
234 
235 C EXCLUDE SINGLE POINT
236  CALL gridman_grid2d_cut(cutgrid,grid,xp,yp,
237  c 5_gridman_sp,.true.,ierr)
238  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
239 
240  CALL gridman_grid_write(cutgrid,'cutAI.grd',ierr)
241  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
242 
243  IF(cutgrid%NCELLS.NE.10.OR.cutgrid%NEDGES.NE.33.OR.
244  f cutgrid%NPOINTS.NE.23.OR.
245  f cutgrid%CELLINDEX(1)%NELEMENTS.NE.10.OR.
246  f cutgrid%EDGEINDEX(1)%NELEMENTS.NE.25) THEN
247  WRITE(*,*) "Unexpected resulting grid (0-I)"
248  stop "TEST_GRID TERMINATED"
249  END IF
250 
251  CALL test_mappings(grid,cutgrid,xp,yp,5_gridman_sp,ierr)
252  IF(ierr.NE.0) THEN
253  WRITE(*,*) "Wrong indices in test 1.3"
254  stop "TEST_GRID TERMINATED"
255  END IF
256 
257  WRITE(*,*) "test_cut: test 1.3 passed"
258 
259 C EMPTY OUTPUT
260  xp=xp+20.
261  CALL gridman_grid2d_cut(cutgrid,grid,xp,yp,
262  c 5_gridman_sp,.false.,ierr)
263  IF(ierr.NE.100) THEN
264  WRITE(*,*) "Unexpected output, expected value 100, IERR ",ierr
265  stop "TEST_CUT TERMINATED"
266  END IF
267 
268 C POLYGON WITH INTERSECTING EDGES, MUST PRODUCE ERROR
269  gridman_check=.true.
270  xp(1)=6.
271  yp(1)=1.
272  xp(3)=6.
273  yp(3)=3.
274  xp(2)=8.
275  yp(2)=3.
276  xp(4)=8.
277  yp(4)=1.
278  xp(5)=xp(1)
279  yp(5)=yp(1)
280  CALL gridman_grid2d_cut(cutgrid,grid,xp,yp,
281  c 5_gridman_sp,.false.,ierr)
282  IF(ierr.NE.100) THEN
283  WRITE(*,*) "Unexpected output, expected value 100, IERR ",ierr
284  stop "TEST_GRID TERMINATED"
285  END IF
286 
287 C VERTEX ON THE EDGE AND EDGE ON THE VERTEX, MUST PRODUCE WARNING MESSAGE
288  gridman_check=.true.
289  xp(1)=6.
290  yp(1)=1.
291  xp(2)=5.
292  yp(2)=4.
293  xp(3)=8.
294  yp(3)=3.
295  xp(4)=8.
296  yp(4)=0.
297  xp(5)=8.
298  yp(5)=-1.
299  xp(6)=xp(1)
300  yp(6)=yp(1)
301  CALL gridman_grid2d_cut(cutgrid,grid,xp,yp,
302  c 6_gridman_sp,.false.,ierr)
303  IF(ierr.NE.400) THEN
304  WRITE(*,*) "Unexpected output, expected value 400, IERR ",ierr
305  stop "TEST_GRID TERMINATED"
306  END IF
307 
308  WRITE(*,*) "test_cut: test 1.4 passed"
309 
310 C TAKE FIRST POLYGON
311  CALL gridman_grid2d_cut(cutgrid,grid,xp1,yp1,np1,.false.,ierr)
312  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
313 
314  CALL gridman_grid_write(cutgrid,'cutB.grd',ierr)
315  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
316 
317  IF(cutgrid%NCELLS.NE.5.OR.cutgrid%NEDGES.NE.19.OR.
318  f cutgrid%NPOINTS.NE.15.OR.
319  f cutgrid%CELLINDEX(1)%NELEMENTS.NE.5.OR.
320  f cutgrid%EDGEINDEX(1)%NELEMENTS.NE.10) THEN
321  WRITE(*,*) "Unexpected resulting grid (1)"
322  stop "TEST_GRID TERMINATED"
323  END IF
324 
325  CALL test_mappings(grid,cutgrid,xp,yp,6_gridman_sp,ierr)
326  IF(ierr.NE.0) THEN
327  WRITE(*,*) "Wrong indices in test 2.1"
328  stop "TEST_GRID TERMINATED"
329  END IF
330 
331  WRITE(*,*) "test_cut: test 2.1 passed"
332 
333 C EXCLUDE FIRST POLYGON
334  CALL gridman_grid2d_cut(cutgrid,grid,xp1,yp1,np1,.true.,ierr)
335  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
336 
337  CALL gridman_grid_write(cutgrid,'cutBI.grd',ierr)
338  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
339 
340  IF(cutgrid%NCELLS.NE.10.OR.cutgrid%NEDGES.NE.31.OR.
341  f cutgrid%NPOINTS.NE.23.OR.
342  f cutgrid%CELLINDEX(1)%NELEMENTS.NE.10.OR.
343  f cutgrid%EDGEINDEX(1)%NELEMENTS.NE.22) THEN
344  WRITE(*,*) "Unexpected resulting grid (1-I)"
345  stop "TEST_GRID TERMINATED"
346  END IF
347 
348  CALL test_mappings(grid,cutgrid,xp1,yp1,np1,ierr)
349  IF(ierr.NE.0) THEN
350  WRITE(*,*) "Wrong indices in test 2.2"
351  stop "TEST_GRID TERMINATED"
352  END IF
353 
354  WRITE(*,*) "test_cut: test 2.2 passed"
355 
356 C TAKE SECOND POLYGON
357 C TAKE ANOTHER GRID
358  CALL grid_example4(grid,ierr)
359  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
360 
361  CALL gridman_grid2d_cut(cutgrid,grid,xp2,yp2,np2,.false.,ierr)
362  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
363 
364  CALL gridman_grid_write(cutgrid,'cutC.grd',ierr)
365  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
366 
367  IF(cutgrid%NCELLS.NE.9.OR.cutgrid%NEDGES.NE.34.OR.
368  f cutgrid%NPOINTS.NE.31.OR.
369  f cutgrid%CELLINDEX(1)%NELEMENTS.NE.9.OR.
370  f cutgrid%EDGEINDEX(1)%NELEMENTS.NE.17) THEN
371  WRITE(*,*) "Unexpected resulting grid (2)"
372  stop "TEST_GRID TERMINATED"
373  END IF
374 
375  CALL test_mappings(grid,cutgrid,xp2,yp2,np2,ierr)
376  IF(ierr.NE.0) THEN
377  WRITE(*,*) "Wrong indices in test 3.1"
378  stop "TEST_GRID TERMINATED"
379  END IF
380 
381  WRITE(*,*) "test_cut: test 3.1 passed"
382 
383 C EXCLUDE SECOND POLYGON
384  CALL gridman_grid2d_cut(cutgrid,grid,xp2,yp2,np2,.true.,ierr)
385  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
386 
387  CALL gridman_grid_write(cutgrid,'cutCI.grd',ierr)
388  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
389 
390  IF(cutgrid%NCELLS.NE.7.OR.cutgrid%NEDGES.NE.42.OR.
391  f cutgrid%NPOINTS.NE.36.OR.
392  f cutgrid%CELLINDEX(1)%NELEMENTS.NE.7.OR.
393  f cutgrid%EDGEINDEX(1)%NELEMENTS.NE.25) THEN
394  WRITE(*,*) "Unexpected resulting grid (4)"
395  stop "TEST_GRID TERMINATED"
396  END IF
397 
398  CALL test_mappings(grid,cutgrid,xp2,yp2,np2,ierr)
399  IF(ierr.NE.0) THEN
400  WRITE(*,*) "Wrong indices in test 3.2"
401  stop "TEST_GRID TERMINATED"
402  END IF
403 
404  WRITE(*,*) "test_cut: test 3.2 passed"
405 
406 C ONE MORE POLYGON
407  CALL test_polygon3(xp3,yp3,np3)
408 C RETURN TO FIRST GRID
409  CALL grid_example1(grid,ierr)
410  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
411 
412  CALL gridman_grid2d_cut(cutgrid,grid,xp3,yp3,np3,.false.,ierr)
413  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
414 
415  CALL gridman_grid_write(cutgrid,'cutD.grd',ierr)
416  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
417 
418  IF(cutgrid%NCELLS.NE.11.OR.cutgrid%NEDGES.NE.39.OR.
419  f cutgrid%NPOINTS.NE.30.OR.
420  f cutgrid%CELLINDEX(1)%NELEMENTS.NE.11.OR.
421  f cutgrid%EDGEINDEX(1)%NELEMENTS.NE.19) THEN
422  WRITE(*,*) "Unexpected resulting grid (5)"
423  stop "TEST_GRID TERMINATED"
424  END IF
425 
426  CALL test_mappings(grid,cutgrid,xp3,yp3,np3,ierr)
427  IF(ierr.NE.0) THEN
428  WRITE(*,*) "Wrong indices in test 4.1"
429  stop "TEST_GRID TERMINATED"
430  END IF
431 
432  WRITE(*,*) "test_cut: test 4.1 passed"
433 
434 C EXCLUDE
435  CALL gridman_grid2d_cut(cutgrid,grid,xp3,yp3,np3,.true.,ierr)
436  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
437 
438  CALL gridman_grid_write(cutgrid,'cutDI.grd',ierr)
439  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
440 
441  IF(cutgrid%NCELLS.NE.16.OR.cutgrid%NEDGES.NE.48.OR.
442  f cutgrid%NPOINTS.NE.36.OR.
443  f cutgrid%CELLINDEX(1)%NELEMENTS.NE.16.OR.
444  f cutgrid%EDGEINDEX(1)%NELEMENTS.NE.28) THEN
445  WRITE(*,*) "Unexpected resulting grid (6)"
446  stop "TEST_GRID TERMINATED"
447  END IF
448 
449  CALL test_mappings(grid,cutgrid,xp3,yp3,np3,ierr)
450  IF(ierr.NE.0) THEN
451  WRITE(*,*) "Wrong indices in test 4.2"
452  stop "TEST_GRID TERMINATED"
453  END IF
454 
455  WRITE(*,*) "test_cut: test 4.2 passed"
456 
457 C DEALLOCATE
458  CALL gridman_grid_deallocate(grid,ierr)
459  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
460  CALL gridman_grid_deallocate(cutgrid,ierr)
461  IF(ierr.NE.0) stop "TEST_CUT TERMINATED"
462 
463  WRITE(*,*) "TEST_CUT COMPLETED"
464 
465  END PROGRAM test_cut
466 C
467 C
468 C
469  SUBROUTINE test_polygon1(X,Y,N)
470  USE gridman
471  IMPLICIT NONE
472  INTEGER(GRIDMAN_SP),INTENT(IN) :: n
473  REAL(GRIDMAN_DP),INTENT(OUT) :: x(n),y(n)
474 
475  x(1)=6.
476  y(1)=1.
477  x(2)=6.
478  y(2)=5.
479  x(3)=8.
480  y(3)=6.
481  x(4)=8.
482  y(4)=3.
483  x(5)=11.
484  y(5)=3.
485  x(6)=11.
486  y(6)=1.
487  x(7)=8.
488  y(7)=0.5
489  x(8)=7.
490  y(8)=-0.5
491  x(9)=x(1)
492  y(9)=y(1)
493 
494  END SUBROUTINE test_polygon1
495 
496 C
497 C
498  SUBROUTINE test_polygon2(X,Y,N)
499  USE gridman
500  IMPLICIT NONE
501  INTEGER(GRIDMAN_SP),INTENT(IN) :: n
502  REAL(GRIDMAN_DP),INTENT(OUT) :: x(n),y(n)
503  REAL(GRIDMAN_DP) :: x0(n-1),y0(n-1)
504 
505  x0(1)=9.
506  y0(1)=-6.
507  x0(2)=9.
508  y0(2)=-2.5
509  x0(3)=17.
510  y0(3)=-2.5
511  x0(4)=17.
512  y0(4)=-1.
513  x0(5)=15.
514  y0(5)=-1.
515  x0(6)=16.
516  y0(6)=1.
517  x0(7)=14.
518  y0(7)=1.
519  x0(8)=13.
520  y0(8)=2.
521  x0(9)=15.
522  y0(9)=2.
523  x0(10)=15.
524  y0(10)=3.5
525  x0(11)=14.
526  y0(11)=4.5
527  x0(12)=13.
528  y0(12)=6.
529  x0(13)=12.
530  y0(13)=6.
531  x0(14)=12.
532  y0(14)=7.
533  x0(15)=12.
534  y0(15)=9.
535  x0(16)=19.
536  y0(16)=9.
537  x0(17)=19.
538  y0(17)=-7.
539 
540  x(1:n-1)=x0
541  y(1:n-1)=y0
542 
543  x(n)=x(1)
544  y(n)=y(1)
545 
546  END SUBROUTINE test_polygon2
547 
548 C
549 C
550  SUBROUTINE test_polygon3(X,Y,N)
551  USE gridman
552  IMPLICIT NONE
553  INTEGER(GRIDMAN_SP),INTENT(IN) :: n
554  REAL(GRIDMAN_DP),INTENT(OUT) :: x(n),y(n)
555  REAL(GRIDMAN_DP) :: x0(n-1),y0(n-1)
556 
557  x0(1)=8.
558  y0(1)=3.
559  x0(2)=9.
560  y0(2)=-5.
561  x0(3)=9.
562  y0(3)=-6.
563  x0(4)=4.5
564  y0(4)=-6.
565  x0(5)=4.5
566  y0(5)=-5.
567  x0(6)=4.5
568  y0(6)=5.
569  x0(7)=5.5
570  y0(7)=5.
571  x0(8)=5.5
572  y0(8)=3.
573  x0(9)=6.
574  y0(9)=5.
575  x0(10)=6.5
576  y0(10)=3.
577  x0(11)=6.5
578  y0(11)=5.
579  x0(12)=7.
580  y0(12)=-5.
581 
582  x(1:n-1)=x0
583  y(1:n-1)=y0
584 
585  x(n)=x(1)
586  y(n)=y(1)
587 
588  END SUBROUTINE test_polygon3
589 
590 C TEST CONSISTENCY OF THE INDEX MAPPINGS
591  SUBROUTINE test_mappings(GRID0,GRID1,XP,YP,NP,IERR)
592  USE gridman
593  USE gridman_lib
594  IMPLICIT NONE
595  INTRINSIC sqrt,abs
596  TYPE(gridman_grid) :: grid0,grid1
597  INTEGER(GRIDMAN_SP),INTENT(IN) :: np
598  REAL(GRIDMAN_DP) :: xp(np),yp(np)
599  INTEGER,INTENT(OUT) :: ierr
600  INTEGER(GRIDMAN_SP) :: iedge1,iedge0,icell,imin,imax,ie,
601  i icell1_1,icell1_2,icell0_1,icell0_2,
602  i ip,iseg,iedge
603  INTEGER(GRIDMAN_SP),ALLOCATABLE :: indcell(:)
604  INTEGER :: res
605  REAL(GRIDMAN_DP) :: ls,le
606  REAL(GRIDMAN_DP),ALLOCATABLE :: dl(:)
607 
608  ierr=0
609 
610  IF(((grid1%NEDGEINDEX.NE.grid0%NEDGEINDEX).AND.
611  f (grid1%NEDGEINDEX.NE.grid0%NEDGEINDEX+1)).OR.
612  f grid1%NCELLINDEX.NE.grid0%NCELLINDEX) THEN
613  WRITE(*,*) "ERROR: wrong number of indices"
614  WRITE(*,*) " NEDGEINDEX0, NEDGEINDEX1 ",
615  w grid0%NEDGEINDEX,grid1%NEDGEINDEX
616  WRITE(*,*) " NCELLINDEX0, NCELLINDEX1 ",
617  w grid0%NCELLINDEX,grid1%NCELLINDEX
618  ierr=1
619  RETURN
620  END IF
621 
622  CALL gridman_index_elmap(indcell,imin,imax,
623  c grid1%CELLINDEX(1),ierr)
624  IF(ierr.NE.0) RETURN
625 
626 C Consistency check between edge and cell index mappings
627  DO iedge1=1,grid1%EDGEINDEX(1)%NELEMENTS
628  IF(iedge1.NE.grid1%EDGEINDEX(1)%INDEXES(0,iedge1)) THEN
629  WRITE(*,*) "ERROR: wrong edge index"
630  WRITE(*,*) " IEDGE1, INDEXES(0,IEDGE1) ",
631  w iedge1, grid1%EDGEINDEX(1)%INDEXES(0,iedge1)
632  ierr=2
633  RETURN
634  END IF
635  iedge0=grid1%EDGEINDEX(1)%INDEXES(1,iedge1)
636  IF(iedge0.LT.1.OR.iedge0.GT.grid0%NEDGES) THEN
637  WRITE(*,*) "ERROR: wrong edge index"
638  WRITE(*,*) " IEDGE0, NEDGES0 ",iedge0,grid0%NEDGES
639  ierr=2
640  RETURN
641  END IF
642  icell=grid1%CELLS(1,iedge1)
643  IF(icell.LT.1) THEN
644  icell1_1=icell
645  ELSE
646  ie=indcell(icell)
647  IF(ie.LT.1) THEN
648  WRITE(*,*) "ERROR: wrong cell index"
649  WRITE(*,*) " ICELL not in the index ",icell
650  ierr=3
651  RETURN
652  END IF
653  icell1_1=grid1%CELLINDEX(1)%INDEXES(1,ie)
654  IF(icell1_1.LT.1.OR.icell1_1.GT.grid0%NCELLS) THEN
655  WRITE(*,*) "ERROR: wrong cell index"
656  WRITE(*,*) " ICELL0, NEDGES0 ",icell1_1,grid0%NCELLS
657  ierr=4
658  RETURN
659  END IF
660  END IF
661  icell=grid1%CELLS(2,iedge1)
662  IF(icell.LT.1) THEN
663  icell1_2=icell
664  ELSE
665  ie=indcell(icell)
666  IF(ie.LT.1) THEN
667  WRITE(*,*) "ERROR: wrong cell index"
668  WRITE(*,*) " ICELL not in the index ",icell
669  ierr=3
670  RETURN
671  END IF
672  icell1_2=grid1%CELLINDEX(1)%INDEXES(1,ie)
673  IF(icell1_2.LT.1.OR.icell1_2.GT.grid0%NCELLS) THEN
674  WRITE(*,*) "ERROR: wrong cell index"
675  WRITE(*,*) " ICELL0, NEDGES0 ",icell1_2,grid0%NCELLS
676  ierr=4
677  RETURN
678  END IF
679  END IF
680  icell0_1=grid0%CELLS(1,iedge0)
681  icell0_2=grid0%CELLS(2,iedge0)
682  IF(icell0_1.NE.icell1_1) THEN
683  IF(icell0_1.NE.icell1_2) THEN
684  WRITE(*,*) "ERROR: inconsistent edge and cell indices"
685  WRITE(*,*) " IEDGE1, IEDGE0 ", iedge1, iedge0
686  WRITE(*,*) " ICELL0_1, ICELL0_2 ", icell0_1,icell0_2
687  WRITE(*,*) " ICELL1_1, ICELL1_2 ", icell1_1,icell1_2
688  ierr=5
689  RETURN
690  END IF
691  ELSEIF(icell0_2.NE.icell1_2) THEN
692  WRITE(*,*) "ERROR: inconsistent edge and cell indices"
693  WRITE(*,*) " IEDGE1, IEDGE0 ", iedge1, iedge0
694  WRITE(*,*) " ICELL0_1, ICELL0_2 ", icell0_1,icell0_2
695  WRITE(*,*) " ICELL1_1, ICELL1_2 ", icell1_1,icell1_2
696  ierr=6
697  RETURN
698  END IF
699  END do! DO IEDGE1=1,GRID1%NEDGES
700 
701  CALL gridman_index_compare(grid1%EDGEINDEX(2),grid1%EDGEINDEX(1),
702  c res,ierr)
703  IF(ierr.NE.0) RETURN
704  IF(res.NE.0) THEN
705  WRITE(*,*) "ERROR: edge mappings are different"
706  ierr=10
707  RETURN
708  END IF
709  CALL gridman_index_compare(grid1%CELLINDEX(2),grid1%CELLINDEX(1),
710  c res,ierr)
711  IF(ierr.NE.0) RETURN
712  IF(res.NE.0) THEN
713  WRITE(*,*) "ERROR: cell mappings are different"
714  ierr=11
715  RETURN
716  END IF
717 
718  DEALLOCATE(indcell)
719 
720 C TEST INDEXING OF THE POLYGON SEGMENTS
721  IF(grid1%NEDGEINDEX.GT.2) THEN
722  ALLOCATE(dl(grid1%NEDGES))
723  CALL gridman_grid2d_lengths(grid1,dl,ierr)
724  IF(ierr.NE.0) RETURN
725  DO ip=1,np-1
726  ls=sqrt((xp(ip)-xp(ip+1))**2+(yp(ip)-yp(ip+1))**2)
727  le=0.
728  DO ie=1,grid1%EDGEINDEX(3)%NELEMENTS
729  iseg=grid1%EDGEINDEX(3)%INDEXES(1,ie)
730  IF(iseg.EQ.ip) THEN
731  iedge=grid1%EDGEINDEX(3)%INDEXES(0,ie)
732  le=le+dl(iedge)*grid1%UNIT2SI
733  END IF
734  END DO
735  IF(le.GT.ls+gridman_tol) THEN
736  WRITE(*,*) "ERROR: incorrect indices of segments"
737  WRITE(*,*) " IP, LS, LE ",ip,ls,le
738  ierr=12
739  RETURN
740  END IF
741  END DO
742  DEALLOCATE(dl)
743  END IF
744 
745 
746  END SUBROUTINE test_mappings
gridman::gridman_check
logical, save, public gridman_check
Switch to enforce extra checks of input parameters.
Definition: gridman.f:133
gridman_lib
Explicit interfaces to GRIDMAN subroutines and functions.
Definition: gridman.f:251
gridman_index_compare
subroutine gridman_index_compare(INDEX1, INDEX2, RES, IERR)
Compare two index objects.
Definition: index.f:423
gridman_grid2d_cut
subroutine gridman_grid2d_cut(CUTGRID, GRID, XP, YP, NP, LEX, IERR)
Select part of a 2D grid cut by polygon.
Definition: cut.f:38
gridman::gridman_grid
Data-type which describes a grid as a set of edges, methods in grid.f.
Definition: gridman.f:168
gridman_grid_deallocate
subroutine gridman_grid_deallocate(GRID, IERR)
Deallocate grid object.
Definition: grid1.f:184
gridman_index_elmap
subroutine gridman_index_elmap(ELMAP, NMIN, NMAX, INDEX, IERR)
Create an aray which maps elements into the table of indices.
Definition: index.f:643
gridman_point_in_polygon
logical function gridman_point_in_polygon(XP, YP, NP, X0, Y0)
Find if point lies inside a closed polygon.
Definition: geom.f:259
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_grid_compare
subroutine gridman_grid_compare(GRID1, GRID2, RES, IERR)
Compare two grid objects.
Definition: grid1.f:1068
Generated on Thu Oct 5 2017 09:49:57 for GRIDMAN by   doxygen 1.8.9.1