GRIDMAN
grid managment library
grid2.f
Go to the documentation of this file.
1 C> @file grid2.f
2 C> Extra and advanced methods of the data-type GRIDMAN_GRID, see gridman.f
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> Combine two grid objects into one
25 C*****************************************************************************************
26 C>
27 C> WARNING: if GRID exists it will be overwritten
28 C>
29 C> Units of GRID1 are used in combined grid
30  SUBROUTINE gridman_grid_combine(GRID,GRID1,GRID2,IERR)
31  USE gridman,ONLY:gridman_grid,gridman_sp,gridman_unit,
32  u gridman_dbg,gridman_check
33  USE gridman_lib,ONLY:gridman_grid_check,gridman_grid_allocate
34  IMPLICIT NONE
35 C> Resulting combined grid
36  TYPE(gridman_grid) :: GRID
37 C> First grid to be combined
38  TYPE(gridman_grid) :: GRID1
39 C> Second grid to be combined
40  TYPE(gridman_grid) :: GRID2
41 C> Error code
42  INTEGER,INTENT(OUT) :: IERR
43 
44  INTEGER :: TYPE,NEDGEINDEX,NCELLINDEX,RES,IC,I,II,IERR0
45  INTEGER(GRIDMAN_SP) :: NCELLS,NEDGES,NPOINTS,
46  i iedge,iedge2,nedges1,ncells1
47 
48  IF(gridman_dbg) WRITE(gridman_unit,*)
49  w "Starting GRIDMAN_GRID_COMBINE"
50 
51  ierr=0
52 
53  IF(gridman_check) THEN
54  CALL gridman_grid_check(grid1,res,ierr)
55  IF(res.NE.0.OR.ierr.GT.0) THEN
56  ierr=100
57  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_COMBINE: ",
58  w "incorrect GRID1"
59  RETURN
60  END IF
61  CALL gridman_grid_check(grid2,res,ierr)
62  IF(res.NE.0.OR.ierr.GT.0) THEN
63  ierr=100
64  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_COMBINE: ",
65  w "incorrect GRID2"
66  RETURN
67  END IF
68  END IF
69 
70  IF(grid1%TYPE.NE.grid2%TYPE.OR.
71  f grid1%EDIM.NE.grid2%EDIM.OR.
72  f grid1%PDIM.NE.grid2%PDIM) THEN
73  ierr=100
74  WRITE(gridman_unit,*)
75  w "ERROR in GRIDMAN_GRID_COMBINE: grids have different type"
76  WRITE(gridman_unit,*)
77  w " GRID1: TYPE, EDIM, PDIM ",grid1%TYPE,grid1%EDIM,grid1%PDIM
78  WRITE(gridman_unit,*)
79  w " GRID2: TYPE, EDIM, PDIM ",grid2%TYPE,grid2%EDIM,grid2%PDIM
80  RETURN
81  END IF
82 
83  TYPE=grid1%TYPE
84  ncells=grid1%NCELLS+grid2%NCELLS
85  nedges=grid1%NEDGES+grid2%NEDGES
86  npoints=grid1%NPOINTS+grid2%NPOINTS
87  nedgeindex=grid1%NEDGEINDEX+grid2%NEDGEINDEX
88  ncellindex=grid1%NCELLINDEX+grid2%NCELLINDEX
89 
90  CALL gridman_grid_allocate(grid,TYPE,NEDGES,NPOINTS,NCELLS,IERR,
91  s nedgeindex,ncellindex)
92  IF(ierr.NE.0) GOTO 1000
93 
94  grid%CELLS(:,1:grid1%NEDGES)=grid1%CELLS(:,1:grid1%NEDGES)
95  ncells1=grid1%NCELLS
96  nedges1=grid1%NEDGES
97  DO iedge2=1,grid2%NEDGES
98  iedge=iedge2+nedges1
99  grid%CELLS(:,iedge)=grid2%CELLS(:,iedge2)
100  DO ic=1,2
101  IF(grid%CELLS(ic,iedge).GT.0) THEN
102  grid%CELLS(ic,iedge)=grid%CELLS(ic,iedge)+ncells1
103  ELSEIF(grid%CELLS(ic,iedge).LT.0) THEN
104  grid%CELLS(ic,iedge)=grid%CELLS(ic,iedge)-nedges1
105  END IF
106  END DO
107  END DO
108  !SHIFT FOR CELL INDICES
109  grid%POINTS(:,1:grid1%NEDGES)=grid1%POINTS(:,1:grid1%NEDGES)
110  grid%POINTS(:,grid1%NEDGES+1:nedges)=
111  = grid2%POINTS(:,1:grid2%NEDGES)+grid1%NPOINTS !SHIFT FOR POINT INDICES
112 
113 C UNITS OF THE FIRST GRID ARE TAKEN
114  grid%UNIT2SI=grid1%UNIT2SI
115  grid%UNITS=grid1%UNITS
116 
117  grid%X(:,1:grid1%NPOINTS)=grid1%X(:,1:grid1%NPOINTS)
118  grid%X(:,grid1%NPOINTS+1:npoints)=
119  = grid2%X(:,1:grid2%NPOINTS)*grid%UNIT2SI/grid2%UNIT2SI
120 
121  DO i=1,grid1%NEDGEINDEX
122  CALL gridman_index_copy(grid%EDGEINDEX(i),
123  c grid1%EDGEINDEX(i),ierr)
124  IF(ierr.NE.0) GOTO 1000
125  END DO
126  DO i=1,grid2%NEDGEINDEX
127  ii=i+grid1%NEDGEINDEX
128  CALL gridman_index_copy(grid%EDGEINDEX(ii),
129  c grid2%EDGEINDEX(ii),ierr)
130  IF(ierr.NE.0) GOTO 1000
131 C SHIFT INDICES OF EDGES
132  grid%EDGEINDEX(ii)%INDEXES(0,:)=
133  = grid%EDGEINDEX(ii)%INDEXES(0,:)+grid1%NEDGES
134  END DO
135 
136  DO i=1,grid1%NCELLINDEX
137  CALL gridman_index_copy(grid%CELLINDEX(i),
138  c grid1%CELLINDEX(i),ierr)
139  IF(ierr.NE.0) GOTO 1000
140  END DO
141  DO i=1,grid2%NCELLINDEX
142  ii=i+grid1%NCELLINDEX
143  CALL gridman_index_copy(grid%CELLINDEX(ii),
144  c grid2%CELLINDEX(ii),ierr)
145  IF(ierr.NE.0) GOTO 1000
146 C SHIFT INDICES OF CELLS
147  grid%CELLINDEX(ii)%INDEXES(0,:)=
148  = grid%CELLINDEX(ii)%INDEXES(0,:)+grid1%NCELLS
149  END DO
150 
151  IF(gridman_check) THEN
152  CALL gridman_grid_check(grid,res,ierr)
153  IF(res.NE.0.OR.ierr.GT.0) THEN
154  ierr=100
155  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_COMBINE: ",
156  w "incorrect resulting grid"
157  CALL gridman_grid_deallocate(grid,ierr0)
158  RETURN
159  END IF
160  END IF
161 
162  IF(gridman_dbg) WRITE(gridman_unit,*)
163  w "GRIDMAN_GRID_COMBINE finished"
164  RETURN
165 
166  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID_COMBINE is terminated"
167 
168  END SUBROUTINE gridman_grid_combine
169 
170 C *****************************************************************************************
171 C> Create a list of edges which belong to each cell
172 C *****************************************************************************************
173 C>
174 C> WARNING: array EDGES will be overwritten if already exists
175 C>
176  SUBROUTINE gridman_grid_cells(EDGES,GRID,IERR)
177  USE gridman,ONLY:gridman_grid,gridman_sp,gridman_unit,
178  u gridman_dbg,gridman_check,gridman_indlist
179  USE gridman_lib,ONLY:gridman_grid_check,gridman_grid_allocate,
180  u gridman_indlist_deallocate
181  IMPLICIT NONE
182 
183 C> Resulting list of edges
184  TYPE(gridman_indlist) :: EDGES
185 C> Grid for which this list is generated
186  TYPE(gridman_grid) :: GRID
187 C> Error code
188  INTEGER,INTENT(OUT) :: IERR
189 
190  INTEGER(GRIDMAN_SP) :: ICELL,IEDGE,IND
191  INTEGER :: ST,IC,RES
192 
193  IF(gridman_dbg)
194  w WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_CELLS"
195 
196  ierr=0
197 
198  IF(grid%NCELLS.LT.1) THEN
199  edges%N=0
200  edges%L=0
201  RETURN
202  END IF
203 
204  IF(gridman_check) THEN
205  CALL gridman_grid_check(grid,res,ierr)
206  IF(res.NE.0.OR.ierr.GT.0) THEN
207  ierr=100
208  WRITE(gridman_unit,*)
209  w "ERROR in GRIDMAN_GRID_CELLS: ",
210  w "incorrect grid object"
211  RETURN
212  END IF
213  END IF
214 
215  CALL gridman_indlist_deallocate(edges,ierr)
216  IF(ierr.NE.0) GOTO 1000
217 
218  ALLOCATE(edges%IFIRST(grid%NCELLS),
219  a edges%ILAST(grid%NCELLS),stat=st)
220  IF(st.NE.0) THEN
221  ierr=200
222  WRITE(gridman_unit,*)
223  w "ERROR in GRIDMAN_GRID_CELLS: ",
224  w "can not allocate IFIRST, ILAST"
225  RETURN
226  END IF
227 
228 C NUMBER OF EDGES FOR EACH CELL
229  edges%L=0
230  edges%ILAST=0
231  DO iedge=1,grid%NEDGES
232  DO ic=1,2
233  icell=grid%CELLS(ic,iedge)
234  IF(icell.GT.0) THEN
235  IF(icell.GT.grid%NCELLS) THEN
236  ierr=100
237  WRITE(gridman_unit,*)
238  w "ERROR in GRIDMAN_GRID_CELLS: ",
239  w "cell index is out of bounds"
240  WRITE(gridman_unit,*) " IEDGE, CELLS, NCELLS ",
241  w iedge,grid%CELLS(1:2,iedge),grid%NCELLS
242  GOTO 100
243  ELSE
244  edges%ILAST(icell)=edges%ILAST(icell)+1
245  edges%L=edges%L+1
246  END IF
247  END IF
248  END DO
249  END DO
250  edges%N=grid%NCELLS
251 
252 C FIRST INDEX WHICH CORRESPONDS TO EACH CELL IN ARRAY IND
253  edges%IFIRST(1)=1
254  DO icell=2,edges%N
255  edges%IFIRST(icell)=edges%IFIRST(icell-1)+edges%ILAST(icell-1)
256  END DO
257 
258  ALLOCATE(edges%IND(edges%L),stat=st)
259  IF(st.NE.0) THEN
260  ierr=200
261  WRITE(gridman_unit,*)
262  w "ERROR in GRIDMAN_GRID_CELLS: ",
263  w "can not allocate IND"
264  GOTO 100
265  END IF
266 
267 C LIST OF EDGES FOR EACH CELL
268  edges%ILAST=edges%IFIRST-1
269  DO iedge=1,grid%NEDGES
270  DO ic=1,2
271  icell=grid%CELLS(ic,iedge)
272  IF(icell.GT.0) THEN
273  edges%ILAST(icell)=edges%ILAST(icell)+1
274  ind=edges%ILAST(icell)
275  IF (ind.GT.edges%L) THEN
276  ierr=400
277  WRITE(gridman_unit,*)
278  w "ERROR in GRIDMAN_GRID_CELLS: internal error"
279  WRITE(gridman_unit,*) " IEDGE, ICELL, IND ",iedge,icell,ind
280  GOTO 100
281  END IF
282  edges%IND(ind)=iedge
283  END IF
284  END DO
285  END DO
286 
287  CALL check
288 
289  IF(gridman_dbg)
290  f WRITE(gridman_unit,*) "GRIDMAN_GRID_CELLS finished"
291 
292  RETURN
293 
294  100 CALL gridman_indlist_deallocate(edges,ierr)
295  RETURN
296 
297  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID_CELLS terminated"
298 
299  CONTAINS
300 C
301 C CHECK THE CREATED TABLE
302 C
303  SUBROUTINE check
304 
305  INTEGER(GRIDMAN_SP) :: ICELL1,ICELL2
306 
307  DO icell=1,edges%N-1
308  IF(edges%IFIRST(icell+1).NE.edges%ILAST(icell)+1) THEN
309  ierr=400
310  WRITE(gridman_unit,*)
311  w "ERROR in GRIDMAN_GRID_CELLS: internal error"
312  WRITE(gridman_unit,*) " ICELL, IFIRST, ILAST ",
313  w icell, edges%IFIRST(icell+1),
314  w edges%ILAST(icell)+1
315  RETURN
316  END IF
317  END DO
318  IF(edges%ILAST(edges%N).NE.edges%L) THEN
319  ierr=400
320  WRITE(gridman_unit,*)
321  w "ERROR in GRIDMAN_GRID_CELLS: internal error"
322  WRITE(gridman_unit,*) " L, ILAST ",edges%L,edges%ILAST(edges%N)
323  RETURN
324  END IF
325 
326  DO icell=1,edges%N
327  DO ind=edges%IFIRST(icell),edges%ILAST(icell)
328  iedge=edges%IND(ind)
329  IF(iedge.LT.1.OR.iedge.GT.grid%NEDGES) THEN
330  ierr=400
331  WRITE(gridman_unit,*)
332  w "ERROR in GRIDMAN_GRID_CELLS: internal error"
333  WRITE(gridman_unit,*) " ICELL, IEDGE, NEDGES ",
334  w icell, iedge, grid%NEDGES
335  RETURN
336  END IF
337  icell1=grid%CELLS(1,iedge)
338  icell2=grid%CELLS(2,iedge)
339  IF(icell1.NE.icell.AND.icell2.NE.icell) THEN
340  ierr=400
341  WRITE(gridman_unit,*)
342  w "ERROR in GRIDMAN_GRID_CELLS: internal error"
343  WRITE(gridman_unit,*) " ICELL, IEDGE, ICELL1, ICELL2 ",
344  w icell, iedge, icell1, icell2
345  RETURN
346  END IF
347  END DO
348  END DO
349  END SUBROUTINE check
350 
351  END SUBROUTINE gridman_grid_cells
352 
353 
354 C *****************************************************************************************
355 C> Create table of edges connected to each point
356 C *****************************************************************************************
357 C>
358 C> WARNING: array EDGES will be overwritten if already exists
359 C>
360  SUBROUTINE gridman_grid_points(EDGES,GRID,IERR)
361  USE gridman,ONLY:gridman_grid,gridman_sp,gridman_unit,
362  u gridman_dbg,gridman_check,gridman_indlist
363  USE gridman_lib,ONLY:gridman_grid_check,gridman_grid_allocate,
364  u gridman_indlist_deallocate
365  IMPLICIT NONE
366 C> Resulting list of edges
367  TYPE(gridman_indlist) :: EDGES
368 C> Grid for which this list is generated
369  TYPE(gridman_grid) :: GRID
370 C> Error code
371  INTEGER,INTENT(OUT) :: IERR
372 
373  INTEGER(GRIDMAN_SP) :: NPOINTS,IEDGE,IPOINT,IND
374  INTEGER :: RES,ST,IP
375 
376  IF(gridman_dbg)
377  w WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_POINTS"
378 
379  ierr=0
380 
381  IF(gridman_check) THEN
382  CALL gridman_grid_check(grid,res,ierr)
383  IF(res.NE.0.OR.ierr.GT.0) THEN
384  ierr=100
385  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_POINTS ",
386  w "incorrect grid object"
387  RETURN
388  END IF
389  END IF
390 
391  npoints=grid%NPOINTS
392 
393  IF(npoints.LT.1) THEN
394  ierr=100
395  WRITE(*,*) "ERROR in GRIDMAN_GRID_POINTS: ",
396  w "number of points is out of range"
397  WRITE(*,*) " NPOINTS ",npoints
398  RETURN
399  END IF
400 
401  CALL gridman_indlist_deallocate(edges,ierr)
402  IF(ierr.NE.0) GOTO 1000
403 
404  ALLOCATE(edges%IFIRST(npoints),edges%ILAST(npoints),stat=st)
405  IF(st.NE.0) THEN
406  ierr=200
407  WRITE(gridman_unit,*)
408  w "ERROR in GRIDMAN_GRID_POINTS: ",
409  w "can not allocate IFIRST, ILAST"
410  RETURN
411  END IF
412 
413 C NUMBER OF EDGES CONNECTED TO EACH POINT
414  edges%L=0
415  edges%ILAST=0
416  DO iedge=1,grid%NEDGES
417  DO ip=1,grid%EDIM
418  ipoint=grid%POINTS(ip,iedge)
419  IF(ipoint.GT.0) THEN
420  IF(ipoint.GT.npoints) THEN
421  ierr=100
422  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_POINTS: ",
423  w "point index is out of bounds"
424  WRITE(gridman_unit,*) " IEDGE, IP, IPOINT, NPOINTS ",
425  w iedge, ip, ipoint, npoints
426  GOTO 100
427  END IF
428  edges%ILAST(ipoint)=edges%ILAST(ipoint)+1
429  edges%L=edges%L+1
430  END IF
431  END DO
432  END DO
433  edges%N=npoints
434 
435 C FIRST INDEX WHICH CORRESPONDS TO EACH POINT IN ARRAY IND
436  edges%IFIRST(1)=1
437  DO ipoint=2,edges%N
438  edges%IFIRST(ipoint)=edges%IFIRST(ipoint-1)+edges%ILAST(ipoint-1)
439  END DO
440 
441  ALLOCATE(edges%IND(edges%L),stat=st)
442  IF(st.NE.0) THEN
443  ierr=200
444  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_POINTS: ",
445  w "can not allocate IND"
446  GOTO 100
447  END IF
448 
449 C LIST OF EDGES FOR EACH POINT
450  edges%ILAST=edges%IFIRST-1
451  DO iedge=1,grid%NEDGES
452  DO ip=1,grid%EDIM
453  ipoint=grid%POINTS(ip,iedge)
454  IF(ipoint.GT.0) THEN
455  edges%ILAST(ipoint)=edges%ILAST(ipoint)+1
456  ind=edges%ILAST(ipoint)
457  IF (ind.GT.edges%L) THEN
458  ierr=400
459  WRITE(gridman_unit,*)
460  w "ERROR in GRIDMAN_GRID_POINTS: internal error"
461  WRITE(gridman_unit,*) " IEDGE, IPOINT, IND ",iedge,ipoint,ind
462  GOTO 100
463  END IF
464  edges%IND(ind)=iedge
465  END IF
466  END DO
467  END DO
468 
469  CALL check
470 
471  IF(gridman_dbg)
472  f WRITE(gridman_unit,*) "GRIDMAN_GRID_POINTS finished"
473 
474  RETURN
475 
476 100 CALL gridman_indlist_deallocate(edges,ierr)
477  RETURN
478 
479 1000 WRITE(gridman_unit,*) "GRIDMAN_GRID_POINTS terminated"
480 
481  CONTAINS
482 
483 C
484 C CHECK THE RESULTING TABLE
485 C
486  SUBROUTINE check
487 
488  INTEGER(GRIDMAN_SP) :: IPOINT1
489 
490  DO ipoint=1,edges%N-1
491  IF(edges%IFIRST(ipoint+1).NE.edges%ILAST(ipoint)+1) THEN
492  ierr=400
493  WRITE(gridman_unit,*)
494  w "ERROR in GRIDMAN_GRID_POINTS: internal error"
495  WRITE(gridman_unit,*) " IPOINT, IFIRST, ILAST ",
496  w ipoint,edges%IFIRST(ipoint+1),
497  w edges%ILAST(ipoint)+1
498  RETURN
499  END IF
500  END DO
501  IF(edges%ILAST(edges%N).NE.edges%L) THEN
502  ierr=400
503  WRITE(gridman_unit,*)
504  w "ERROR in GRIDMAN_GRID_POINTS: internal error"
505  WRITE(gridman_unit,*) " L, ILAST ",
506  w edges%L,edges%ILAST(edges%N)
507  RETURN
508  END IF
509 
510  DO ipoint=1,edges%N
511  DO ind=edges%IFIRST(ipoint),edges%ILAST(ipoint)
512  iedge=edges%IND(ind)
513  IF(iedge.LT.1.OR.iedge.GT.grid%NEDGES) THEN
514  ierr=400
515  WRITE(gridman_unit,*)
516  w "ERROR in GRIDMAN_GRID_POINTS: internal error"
517  WRITE(gridman_unit,*) " IPOINT, IEDGE, NEDGES ",
518  w ipoint, iedge, grid%NEDGES
519  RETURN
520  END IF
521  DO ip=1,grid%EDIM
522  ipoint1=grid%POINTS(ip,iedge)
523  IF(ipoint1.EQ.ipoint) GOTO 77
524  END DO
525  ierr=400
526  WRITE(gridman_unit,*)
527  w "ERROR in GRIDMAN_GRID_POINTS: internal error"
528  WRITE(gridman_unit,*) " IPOINT, IEDGE, IPOINTS ",
529  w ipoint, iedge, grid%POINTS(:,iedge)
530  RETURN
531  77 CONTINUE
532  END DO
533  END DO
534  END SUBROUTINE check
535 
536  END SUBROUTINE gridman_grid_points
537 
538 C *****************************************************************************************
539 C> Eliminate cells from GRIDMAN_GRID object
540 C>
541 C> WARNING: object GRID_NEW will be overwritten if already exists
542 C *****************************************************************************************
543  SUBROUTINE gridman_grid_eliminate_cells(GRID_NEW,GRID,LTAKE,IERR)
544  USE gridman,ONLY:gridman_grid,gridman_sp,gridman_unit,
545  u gridman_dbg,gridman_check
546  USE gridman_lib,ONLY:gridman_index_allocate,gridman_grid_check
547  IMPLICIT NONE
548  INTRINSIC trim
549 C> Resulting grid with eliminated cells
550  TYPE(gridman_grid) :: GRID_NEW
551 C> Original grid
552  TYPE(gridman_grid) :: GRID
553 C> List of cells which have to be included/eliminated
554 C>
555 C> If LTAKE(ICELL)=.TRUE. then ICELL is taken into GRID_NEW,
556 C> otherwise ICELL is eliminated
557  LOGICAL,INTENT(IN) :: LTAKE(grid%ncells)
558 C> Error code
559  INTEGER,INTENT(OUT) :: IERR
560 
561 C LOCAL VARIABLES
562  INTEGER(GRIDMAN_SP) :: NCELLS,NEDGES,NPOINTS,
563  i icell,iedge,ipoint,
564  i iedge0,icell0,ipoint0,
565  i ie,ie0,nelements
566  INTEGER(GRIDMAN_SP) :: INEWCELL(grid%ncells),
567  i inewedge(grid%NEDGES),
568  i inewpoint(grid%NPOINTS)
569  LOGICAL :: LTAKE_POINT(grid%npoints)
570  INTEGER :: IC,IP,II,RES
571 
572  IF(gridman_dbg)
573  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_ELIMINATE_CELLS"
574 
575  ierr=0
576 
577  IF(grid%NCELLS.LT.1) THEN
578  WRITE(gridman_unit,*)
579  w "WARNING from GRIDMAN_GRID_ELIMINATE_CELLS: ",
580  w "grid has no cells"
581  ierr=-100
582  RETURN
583  END IF
584 
585  IF(gridman_check) THEN
586  CALL gridman_grid_check(grid,res,ierr)
587  IF(res.NE.0.OR.ierr.GT.0) THEN
588  WRITE(gridman_unit,*)
589  w "ERROR in GRIDMAN_GRID_ELIMINATE_CELLS: ",
590  w "incorrect input grid"
591  ierr=100
592  RETURN
593  END IF
594  END IF
595 
596  !MAPPING OF CELLS
597  ncells=0
598  inewcell=0
599  DO icell=1,grid%NCELLS
600  IF(ltake(icell)) THEN
601  ncells=ncells+1
602  inewcell(icell)=ncells
603  END IF
604  END DO
605 
606  !MAPPING OF EDGES
607  nedges=0
608  inewedge=0
609  DO iedge=1,grid%NEDGES
610  IF(take_edge(iedge)) THEN
611  nedges=nedges+1
612  inewedge(iedge)=nedges
613  END IF
614  END DO
615 
616  !MAPPING OF POINTS
617  ltake_point=.false.
618  DO iedge=1,grid%NEDGES
619  IF(inewedge(iedge).EQ.0) cycle
620  DO ip=1,grid%EDIM
621  ipoint=grid%POINTS(ip,iedge)
622  IF(ipoint.GT.0) THEN
623  ltake_point(ipoint)=.true.
624  ELSE
625  GOTO 100
626  END IF
627  END DO
628  END DO
629  npoints=0
630  inewpoint=0
631  DO ipoint=1,grid%NPOINTS
632  IF(ltake_point(ipoint)) THEN
633  npoints=npoints+1
634  inewpoint(ipoint)=npoints
635  END IF
636  END DO
637 
638  CALL gridman_grid_allocate(grid_new,grid%TYPE,
639  s nedges,npoints,ncells,ierr,
640  s grid%NEDGEINDEX,grid%NCELLINDEX)
641 
642 C EDGES AND POINTS
643  DO iedge0=1,grid%NEDGES
644  iedge=inewedge(iedge0)
645  IF(iedge.LT.1) cycle
646  DO ic=1,2
647  icell0=grid%CELLS(ic,iedge0)
648  IF(icell0.GT.0) THEN
649  icell=inewcell(icell0)
650  ELSEIF(icell0.LT.0) THEN
651  icell=inewedge(-icell0)
652  ELSE
653  icell=0
654  END IF
655  grid_new%CELLS(ic,iedge)=icell
656  END DO
657  DO ip=1,grid%EDIM
658  ipoint0=grid%POINTS(ip,iedge0)
659  IF(ipoint0.GT.0) THEN
660  ipoint=inewpoint(ipoint0)
661  IF(ipoint.GT.0) THEN
662  grid_new%POINTS(ip,iedge)=ipoint
663  ELSE
664  GOTO 200
665  END IF
666  ELSE
667  GOTO 200
668  END IF
669  END DO
670  END DO
671 
672 C COORDINATES
673  DO ipoint0=1,grid%NPOINTS
674  ipoint=inewpoint(ipoint0)
675  IF(ipoint.GT.0) grid_new%X(:,ipoint)= grid%X(:,ipoint0)
676  END DO
677 
678 C DESCRIPTION AND UNITS
679  grid_new%DESCRIPTION=trim(grid%DESCRIPTION)//
680  / " >> after GRIDMAN_GRID_ELIMINATE_CELLS"
681  grid_new%UNIT2SI=grid%UNIT2SI
682  grid_new%UNITS=grid%UNITS
683 
684 C EDGE INDICES
685  DO ii=1,grid%NEDGEINDEX
686  nelements=0
687  DO ie0=1,grid%EDGEINDEX(ii)%NELEMENTS
688  iedge0=grid%EDGEINDEX(ii)%INDEXES(0,ie0)
689  IF(iedge0.LT.0.OR.iedge0.GT.grid%NEDGES) GOTO 300
690  iedge=inewedge(iedge0)
691  IF(iedge.GT.0) nelements=nelements+1
692  END DO
693  CALL gridman_index_allocate(grid_new%EDGEINDEX(ii),
694  c grid%EDGEINDEX(ii)%NINDEX,
695  c nelements,ierr)
696  IF(ierr.NE.0) GOTO 1000
697  ie=0
698  DO ie0=1,grid%EDGEINDEX(ii)%NELEMENTS
699  iedge0=grid%EDGEINDEX(ii)%INDEXES(0,ie0)
700  iedge=inewedge(iedge0)
701  IF(iedge.GT.0) THEN
702  ie=ie+1
703  grid_new%EDGEINDEX(ii)%INDEXES(0,ie)=iedge
704  grid_new%EDGEINDEX(ii)%INDEXES(1:,ie)=
705  = grid%EDGEINDEX(ii)%INDEXES(1:,ie0)
706  END IF
707  END DO
708  grid_new%EDGEINDEX(ii)%DESCRIPTION=grid%EDGEINDEX(ii)%DESCRIPTION
709  grid_new%EDGEINDEX(ii)%COLUMNS=grid%EDGEINDEX(ii)%COLUMNS
710  END DO
711 
712 C CELL INDICES
713  DO ii=1,grid%NCELLINDEX
714  nelements=0
715  DO ie0=1,grid%CELLINDEX(ii)%NELEMENTS
716  icell0=grid%CELLINDEX(ii)%INDEXES(0,ie0)
717  IF(icell0.LT.0.OR.icell0.GT.grid%NCELLS) GOTO 400
718  icell=inewcell(icell0)
719  IF(icell.GT.0) nelements=nelements+1
720  END DO
721  CALL gridman_index_allocate(grid_new%CELLINDEX(ii),
722  c grid%CELLINDEX(ii)%NINDEX,
723  c nelements,ierr)
724  IF(ierr.NE.0) GOTO 1000
725  ie=0
726  DO ie0=1,grid%CELLINDEX(ii)%NELEMENTS
727  icell0=grid%CELLINDEX(ii)%INDEXES(0,ie0)
728  icell=inewcell(icell0)
729  IF(icell.GT.0) THEN
730  ie=ie+1
731  grid_new%CELLINDEX(ii)%INDEXES(0,ie)=icell
732  grid_new%CELLINDEX(ii)%INDEXES(1:,ie)=
733  = grid%CELLINDEX(ii)%INDEXES(1:,ie0)
734  END IF
735  END DO
736  grid_new%CELLINDEX(ii)%DESCRIPTION=grid%CELLINDEX(ii)%DESCRIPTION
737  grid_new%CELLINDEX(ii)%COLUMNS=grid%CELLINDEX(ii)%COLUMNS
738  END DO
739 
740  CALL gridman_grid_check(grid_new,res,ierr)
741  IF(res.NE.0.OR.ierr.GT.0) THEN
742  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_CELLS: ",
743  w "incorrect resulting grid"
744  ierr=400
745  RETURN
746  END IF
747 
748  IF(gridman_dbg)
749  f WRITE(gridman_unit,*) "GRIDMAN_GRID_ELIMINATE_CELLS completed"
750 
751  RETURN
752 
753  100 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_CELLS: ",
754  w "incorrect point index"
755  WRITE(gridman_unit,*) " IEDGE, IP, IPOINT ",
756  w iedge, ip, ipoint
757  ierr=100
758  RETURN
759 
760  200 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_CELLS: ",
761  w "internal error"
762  WRITE(gridman_unit,*) " IEDGE, IEDGE0, IP, IPOINT0, IPOINT ",
763  w iedge, iedge0, ip, ipoint0, ipoint
764  ierr=400
765  RETURN
766 
767  300 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_CELLS: ",
768  w "edge index is out of range"
769  WRITE(gridman_unit,*) " II, IE0, IEDGE0, NEDGES0 ",
770  w ii, ie0, iedge0, grid%NEDGES
771  ierr=100
772  RETURN
773 
774  400 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_CELLS: ",
775  w "cell index is out of range"
776  WRITE(gridman_unit,*) " II, IE0, ICELL0, NCELLS0 ",
777  w ii, ie0, icell0, grid%NCELLS
778  ierr=100
779  RETURN
780 
781  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID_ELIMINATE_CELLS terminated"
782 
783  CONTAINS
784 
785 C EDGE IS TAKEN INTO NEW GRID: TAKE_EDGE=.TRUE., IF
786 C ICELL1>0 .AND. LTAKE(ICELL1) .OR.
787 C ICELL2>0 .AND. LTAKE(ICELL2) .OR.
788 C ICELL1<=0 .AND. ICELL2<=0
789  FUNCTION take_edge(IEDGE)
790  LOGICAL :: TAKE_EDGE
791  INTEGER(GRIDMAN_SP),INTENT(IN) :: IEDGE
792  INTEGER(GRIDMAN_SP) :: ICELL1,ICELL2
793 
794  icell1=grid%CELLS(1,iedge)
795  icell2=grid%CELLS(2,iedge)
796 
797  take_edge=.false.
798  IF(icell1.GT.0) THEN
799  IF(ltake(icell1)) take_edge=.true.
800  END IF
801  IF(icell2.GT.0) THEN
802  IF(ltake(icell2)) take_edge=.true.
803  END IF
804  IF(icell1.LT.1.AND.icell2.LT.1) take_edge=.true.
805 
806  END FUNCTION take_edge
807 
808  END SUBROUTINE gridman_grid_eliminate_cells
809 
810 C *****************************************************************************************
811 C> Eliminate edges the GRIDMAN_GRID object
812 C>
813 C> WARNING: object GRID will be overwritten if already exists
814 C *****************************************************************************************
815  SUBROUTINE gridman_grid_eliminate_edges(GRID_NEW,GRID,LTAKE,IERR)
816  USE gridman,ONLY:gridman_grid,gridman_sp,gridman_unit,
817  u gridman_dbg,gridman_check
818  USE gridman_lib,ONLY:gridman_index_allocate,gridman_grid_check,
819  u gridman_index_copy,gridman_index_deallocate
820  IMPLICIT NONE
821  INTRINSIC trim
822 C> Resulting grid with removed cells
823  TYPE(gridman_grid) :: GRID_NEW
824 C> Original grid
825  TYPE(gridman_grid) :: GRID
826 C> List of edges which have to be included/eliminated, LTAKE(GRID%NEDGES)
827 C>
828 C> If LTAKE(IEDGE)=.TRUE., then IEDGE is taken into GRID_NEW,
829 C> otherwise IEDGE is eliminated
830  LOGICAL,INTENT(IN) :: LTAKE(grid%nedges)
831 C> Error code
832  INTEGER,INTENT(OUT) :: IERR
833 
834 C LOCAL VARIABLES
835  INTEGER(GRIDMAN_SP) :: IEDGE,IPOINT,NEDGES,NPOINTS,
836  i iedge1,ipoint1,ie,ie1,nelements
837  INTEGER(GRIDMAN_SP) :: INEWEDGE(grid%nedges),
838  i inewpoint(grid%NPOINTS)
839  LOGICAL :: LTAKE_POINT(grid%npoints)
840  INTEGER :: IP,RES,II,IERR0
841 
842  IF(gridman_dbg)
843  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_ELIMINATE_EDGES"
844 
845  ierr=0
846 
847  IF(gridman_check) THEN
848  CALL gridman_grid_check(grid,res,ierr)
849  IF(res.NE.0.OR.ierr.GT.0) THEN
850  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_EDGES: ",
851  w "incorrect input grid object"
852  ierr=100
853  RETURN
854  END IF
855  END IF
856 
857  inewedge=0
858  nedges=0
859  DO iedge=1,grid%NEDGES
860  IF(ltake(iedge)) THEN
861  nedges=nedges+1
862  inewedge(iedge)=nedges
863  END IF
864  END DO
865 
866  IF(nedges.LT.1) THEN
867  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_EDGES: ",
868  w "no edges left - the resulting grid would be empty"
869  ierr=100
870  RETURN
871  END IF
872 
873  ltake_point=.false.
874  DO iedge=1,grid%NEDGES
875  IF(inewedge(iedge).GT.0) THEN
876  DO ip=1,grid%EDIM
877  ipoint=grid%POINTS(ip,iedge)
878  ltake_point(ipoint)=.true.
879  END DO
880  END IF
881  END DO
882 
883  inewpoint=0
884  npoints=0
885  DO ipoint=1,grid%NPOINTS
886  IF(ltake_point(ipoint)) THEN
887  npoints=npoints+1
888  inewpoint(ipoint)=npoints
889  END IF
890  END DO
891 
892  CALL gridman_grid_allocate(grid_new,grid%TYPE,
893  c nedges,npoints,grid%NCELLS,ierr,
894  c grid%NEDGEINDEX,grid%NCELLINDEX)
895  IF(ierr.NE.0)
896  c WRITE(gridman_unit,*) "GRIDMAN_GRID_ELIMINATE_EDGES terminated"
897 
898  DO iedge=1,grid%NEDGES
899  iedge1=inewedge(iedge)
900  IF(iedge1.GT.0) THEN
901  grid_new%CELLS(:,iedge1)=grid%CELLS(:,iedge)
902  DO ip=1,grid%EDIM
903  ipoint=grid%POINTS(ip,iedge)
904  ipoint1=inewpoint(ipoint)
905  IF(ipoint1.LT.1) GOTO 400
906  grid_new%POINTS(ip,iedge1)=ipoint1
907  END DO
908  END IF
909  END DO
910  DO ipoint=1,grid%NPOINTS
911  ipoint1=inewpoint(ipoint)
912  IF(ipoint1.GT.0) grid_new%X(:,ipoint1)=grid%X(:,ipoint)
913  END DO
914 
915  grid_new%DESCRIPTION=trim(grid%DESCRIPTION)//
916  / " >> after GRIDMAN_GRID_ELIMINATE_EDGES"
917 
918  grid_new%UNITS=grid%UNITS
919  grid_new%UNIT2SI=grid%UNIT2SI
920 
921  DO ii=1,grid%NCELLINDEX
922  CALL gridman_index_copy(grid_new%CELLINDEX(ii),
923  c grid%CELLINDEX(ii),ierr)
924  IF(ierr.NE.0) GOTO 200
925  END DO
926 
927  DO ii=1,grid%NEDGEINDEX
928  nelements=0
929  DO ie=1,grid%EDGEINDEX(ii)%NELEMENTS
930  iedge=grid%EDGEINDEX(ii)%INDEXES(0,ie)
931  IF(iedge.LT.1.OR.iedge.GT.grid%NEDGES) GOTO 100
932  iedge1=inewedge(iedge)
933  IF(iedge1.GT.0) nelements=nelements+1
934  END DO
935  CALL gridman_index_allocate(grid_new%EDGEINDEX(ii),
936  c grid%EDGEINDEX(ii)%NINDEX,
937  c nelements,ierr)
938  IF(ierr.NE.0) GOTO 210
939  ie1=0
940  DO ie=1,grid%EDGEINDEX(ii)%NELEMENTS
941  iedge=grid%EDGEINDEX(ii)%INDEXES(0,ie)
942  iedge1=inewedge(iedge)
943  IF(iedge1.GT.0) THEN
944  ie1=ie1+1
945  grid_new%EDGEINDEX(ii)%INDEXES(0,ie1)=iedge1
946  grid_new%EDGEINDEX(ii)%INDEXES(1:,ie1)=
947  = grid%EDGEINDEX(ii)%INDEXES(1:,ie)
948  END IF
949  END DO
950  grid_new%EDGEINDEX(ii)%DESCRIPTION=grid%EDGEINDEX(ii)%DESCRIPTION
951  grid_new%EDGEINDEX(ii)%COLUMNS=grid%EDGEINDEX(ii)%COLUMNS
952  END DO
953 
954  CALL gridman_grid_check(grid_new,res,ierr)
955  IF(res.NE.0.OR.ierr.GT.0) THEN
956  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_EDGES: ",
957  w "incorrect resulting grid"
958  CALL gridman_grid_deallocate(grid_new,ierr0)
959  ierr=400
960  RETURN
961  END IF
962 
963  IF(gridman_dbg)
964  f WRITE(gridman_unit,*) "GRIDMAN_GRID_ELIMINATE_EDGES finished"
965 
966  RETURN
967 
968  100 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_EDGES: ",
969  w "incorerct input edge index"
970  WRITE(gridman_unit,*) " II, IE, IEDGE ", ii, ie, iedge
971  CALL gridman_grid_deallocate(grid_new,ierr0)
972  RETURN
973  210 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_EDGES: ",
974  w "creation of celledge index failed, II ",ii
975  CALL gridman_grid_deallocate(grid_new,ierr0)
976  RETURN
977  200 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_EDGES: ",
978  w "creation of cell index failed, II ",ii
979  CALL gridman_grid_deallocate(grid_new,ierr0)
980  RETURN
981  400 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ELIMINATE_EDGES: ",
982  w "internal error "
983  WRITE(gridman_unit,*) " Mismatch between INEWEDGE and INEWPOINT"
984  WRITE(gridman_unit,*) " IEDGE, IEDGE1 ",iedge,iedge1
985  WRITE(gridman_unit,*) " IP, IPOINT, IPOINT1 ",ip,ipoint,ipoint1
986  ierr=400
987  CALL gridman_grid_deallocate(grid_new,ierr0)
988  RETURN
989 
990  END SUBROUTINE gridman_grid_eliminate_edges
991 
992 
993 C *****************************************************************************************
994 C> Remove edges which do not belong to any cell from the GRIDMAN_GRID object
995 C>
996 C> WARNING: object GRID will be overwritten if already exists
997 C *****************************************************************************************
998  SUBROUTINE gridman_grid_remove_free_edges(GRID_NEW,GRID,IERR)
999  USE gridman,ONLY:gridman_grid,gridman_unit,gridman_dbg,gridman_sp
1000  USE gridman_lib,ONLY:gridman_grid_eliminate_edges
1001  IMPLICIT NONE
1002  INTRINSIC trim
1003 C> Resulting grid with removed cells
1004  TYPE(gridman_grid) :: GRID_NEW
1005 C> Original grid
1006  TYPE(gridman_grid) :: GRID
1007 C> Error code
1008  INTEGER,INTENT(OUT) :: IERR
1009 
1010  LOGICAL :: LTAKE(grid%nedges)
1011  INTEGER(GRIDMAN_SP) :: IEDGE
1012 
1013  IF(gridman_dbg)
1014  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_REMOVE_FREE_EDGES"
1015 
1016  DO iedge=1,grid%NEDGES
1017  IF(grid%CELLS(1,iedge).LT.1.AND.
1018  f grid%CELLS(2,iedge).LT.1) THEN
1019  ltake(iedge)=.false.
1020  ELSE
1021  ltake(iedge)=.true.
1022  END IF
1023  END DO
1024 
1025  CALL gridman_grid_eliminate_edges(grid_new,grid,ltake,ierr)
1026  IF(ierr.NE.0)
1027  f WRITE(gridman_unit,*)
1028  w "GRIDMAN_GRID_REMOVE_FREE_EDGES terminated"
1029 
1030  IF(gridman_dbg)
1031  f WRITE(gridman_unit,*) "GRIDMAN_GRID_REMOVE_FREE_EDGES finished"
1032 
1033  END SUBROUTINE gridman_grid_remove_free_edges
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::gridman_unit
integer, save, public gridman_unit
Index of the standard output unit.
Definition: gridman.f:116
gridman_grid_eliminate_edges
subroutine gridman_grid_eliminate_edges(GRID_NEW, GRID, LTAKE, IERR)
Eliminate edges the GRIDMAN_GRID object.
Definition: grid2.f:816
gridman_grid_eliminate_cells
subroutine gridman_grid_eliminate_cells(GRID_NEW, GRID, LTAKE, IERR)
Eliminate cells from GRIDMAN_GRID object.
Definition: grid2.f:544
gridman_lib
Explicit interfaces to GRIDMAN subroutines and functions.
Definition: gridman.f:251
gridman_grid_allocate
subroutine gridman_grid_allocate(GRID, TYPE, NEDGES, NPOINTS, NCELLS, IERR, NEDGEINDEX, NCELLINDEX)
Allocate GRIDMAN_GRID object.
Definition: grid1.f:30
gridman_index_copy
subroutine gridman_index_copy(INDEX2, INDEX1, IERR)
Create a copy of the index object.
Definition: index.f:538
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_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_grid_deallocate
subroutine gridman_grid_deallocate(GRID, IERR)
Deallocate grid object.
Definition: grid1.f:184
gridman_index_allocate
subroutine gridman_index_allocate(INDEX, NINDEX, NELEMENTS, IERR)
Allocate index object.
Definition: index.f:28
gridman::gridman_dbg
logical, save, public gridman_dbg
Switch for debugging mode.
Definition: gridman.f:122
gridman_grid_combine
subroutine gridman_grid_combine(GRID, GRID1, GRID2, IERR)
Combine two grid objects into one.
Definition: grid2.f:31
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_index_deallocate
subroutine gridman_index_deallocate(INDEX, IERR)
Allocate index object.
Definition: index.f:98
gridman_grid_remove_free_edges
subroutine gridman_grid_remove_free_edges(GRID_NEW, GRID, IERR)
Remove edges which do not belong to any cell from the GRIDMAN_GRID object.
Definition: grid2.f:999
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