The TELEMAC-MASCARET system  trunk
mod_handle_sections.f
Go to the documentation of this file.
1 ! **************************
2  MODULE mod_handle_sections
3 ! **************************
4 !
5 !***********************************************************************
6 ! PARTEL
7 !***********************************************************************
8 !
9 !BRIEF Treatment of sections
10 !
11 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
12 
13  IMPLICIT NONE
14  CONTAINS
15 ! **************************
16  SUBROUTINE handle_sections
17 ! **************************
18 !
19  & (namesec, nparts, nelem, ndp, ikle, npoin, f, knogl)
20 !
21 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
22 !| NAMESEC |<--| Name of sections file
23 !| NPARTS |<--| Number of partitions
24 !| NELEM |<--| Number of elements
25 !| NDP |<--| Number of points per element
26 !| IKLE |<->| Connectiviy array
27 !| NPOIN |<--| Number of points
28 !| F |<->| Coordinates
29 !| KNOGL |<--| Global to local numbering
30 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
31 !
32  USE declarations_special
33  USE declarations_partel
34  USE mod_hash_table
35  USE bief, ONLY : chain_type,nbmaxnshare
36 !
37 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
38 !
39  CHARACTER(LEN=PATH_LEN), INTENT(IN) :: NAMESEC
40  INTEGER, INTENT(IN) :: NELEM, NDP
41  INTEGER, INTENT(IN) :: IKLE(nelem,ndp)
42  INTEGER, INTENT(IN) :: NPOIN
43  INTEGER, INTENT(IN) :: NPARTS
44  DOUBLE PRECISION, INTENT(IN) :: F(npoin,2)
45  TYPE(hash_table), INTENT(IN) :: KNOGL
46 !
47 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
48 !
49  INTEGER NSCT
50  type(chain_type), ALLOCATABLE :: chain(:)
51  INTEGER, PARAMETER :: NSEMAX=500 ! MAX NUMBER OF SEGMENTS IN A SECTION
52  INTEGER, ALLOCATABLE :: LISTE(:,:), ANPBEG(:),ANPEND(:)
53  INTEGER :: NSEC, IHOWSEC, ISEC, IELEM, IM(1), IN(1), NPBEG, NPEND
54  INTEGER :: PT, I1,I2,I3, ARR,DEP, ILPREC,ILBEST,ELBEST,IGBEST
55  DOUBLE PRECISION :: XA, YA, DISTB, DISTE, DMINB, DMINE
56  DOUBLE PRECISION :: DIST1, DIST2, DIST3, DIST
57  LOGICAL FOUND
58  INTEGER :: I, IERR, ISEG, J, K, M, N
59  INTEGER :: TEMP
60  CHARACTER(LEN=11) :: EXTENS
61  EXTERNAL extens
62 
63  CHARACTER(LEN=PATH_LEN) :: NAMEOUT
64 !
65  CALL get_free_id(nsct)
66  OPEN (nsct,file=trim(namesec),form='FORMATTED',status='OLD')
67  READ (nsct,*) ! COMMENT LINE
68  READ (nsct,*) nsec, ihowsec
69  IF (.NOT.ALLOCATED(chain)) ALLOCATE (chain(nsec))
70  IF (ihowsec<0) THEN
71  DO isec=1,nsec
72  READ (nsct,*) chain(isec)%DESCR
73  READ (nsct,*) chain(isec)%NPAIR(:)
74  chain(isec)%XYBEG(1)=f(chain(isec)%NPAIR(1),1)
75  chain(isec)%XYBEG(2)=f(chain(isec)%NPAIR(1),2)
76  chain(isec)%XYEND(1)=f(chain(isec)%NPAIR(2),1)
77  chain(isec)%XYEND(2)=f(chain(isec)%NPAIR(2),2)
78  WRITE(lu,*) 'SECTION ',chain(isec)%DESCR
79  WRITE(lu,*) 'BEGINS AT X=',chain(isec)%XYBEG(1),
80  & ' Y=',chain(isec)%XYBEG(2)
81  WRITE(lu,*) 'ENDS AT X=',chain(isec)%XYEND(1),
82  & ' Y=',chain(isec)%XYEND(2)
83  ENDDO
84  ELSE
85  DO isec=1,nsec
86  READ (nsct,*) chain(isec)%DESCR
87  READ (nsct,*)
88  & ( chain(isec)%XYBEG(i), i = 1, SIZE(chain(isec)%XYBEG(:),1) ),
89  & ( chain(isec)%XYEND(i), i = 1, SIZE(chain(isec)%XYEND(:),1) )
90  chain(isec)%NPAIR(:)=0
91  ENDDO
92  ENDIF
93  CLOSE(nsct)
94 !
95 ! IF TERMINAL POINTS GIVEN BY COORDINATES, FIND NEAREST NODES FIRST
96 !
97  WRITE(lu,*) 'NPOIN:',npoin
98  IF(ihowsec.GE.0) THEN
99  DO isec=1,nsec
100  xa=f(1,1)
101  ya=f(1,2)
102  dminb = (chain(isec)%XYBEG(1)-xa)**2
103  & + (chain(isec)%XYBEG(2)-ya)**2
104  dmine = (chain(isec)%XYEND(1)-xa)**2
105  & + (chain(isec)%XYEND(2)-ya)**2
106  chain(isec)%NPAIR(1)=1
107  chain(isec)%NPAIR(2)=1
108  DO i=2,npoin ! COMPUTATIONALLY INTENSIVE
109  xa=f(i,1)
110  ya=f(i,2)
111  distb = (chain(isec)%XYBEG(1)-xa)**2
112  & + (chain(isec)%XYBEG(2)-ya)**2
113  diste = (chain(isec)%XYEND(1)-xa)**2
114  & + (chain(isec)%XYEND(2)-ya)**2
115  IF ( distb < dminb ) THEN
116  chain(isec)%NPAIR(1)=i
117  dminb=distb
118  ENDIF
119  IF ( diste < dmine ) THEN
120  chain(isec)%NPAIR(2)=i
121  dmine=diste
122  ENDIF
123  ENDDO
124  WRITE(lu,'(A,3(1X,I9))')
125  & ' -> SECTION, TERMINAL NODES: ',
126  & isec, chain(isec)%NPAIR(:)
127  ENDDO
128  ELSE
129  DO isec=1,nsec
130  WRITE(lu,'(A,1X,I9,4(1X,1PG13.6))')
131  & ' -> SECTION, TERMINAL COORDINATES: ', isec,
132  & chain(isec)%XYBEG, chain(isec)%XYEND
133  ENDDO
134  ENDIF
135 !
136  WRITE(lu,*) 'NSEC,IHOWSEC: ',nsec,ihowsec
137  WRITE(lu,*) 'ANTICIPATED SECTIONS SUMMARY:'
138  DO isec=1,nsec
139  WRITE(lu,*) chain(isec)%DESCR
140  WRITE(lu,*) chain(isec)%XYBEG(:), chain(isec)%XYEND(:)
141  WRITE(lu,*) chain(isec)%NPAIR(:)
142  ENDDO
143 !
144 ! NOW FOLLOW THE FLUSEC SUBROUTINE IN BIEF TO FIND SECTIONS
145 ! IN THE GLOBAL MESH -> FILL THE FIELD LISTE
146 !
147  ALLOCATE(liste(nsemax,2),stat=ierr) ! WORKHORSE
148  CALL check_allocate(ierr, 'LISTE')
149 !
150  DO isec =1,nsec
151 !
152  dep = chain(isec)%NPAIR(1)
153  arr = chain(isec)%NPAIR(2)
154 !
155  pt = dep
156  iseg = 0
157  dist=(f(dep,1)-f(arr,1))**2+(f(dep,2)-f(arr,2))**2
158 !
159  1010 CONTINUE ! A JUMP POINT
160 !
161  DO ielem =1,nelem
162  i1 = ikle(ielem,1)
163  i2 = ikle(ielem,2)
164  i3 = ikle(ielem,3)
165  IF (pt.EQ.i1.OR.pt.EQ.i2.OR.pt.EQ.i3) THEN
166  dist1 = (f(i1,1)-f(arr,1))**2 + (f(i1,2)-f(arr,2))**2
167  dist2 = (f(i2,1)-f(arr,1))**2 + (f(i2,2)-f(arr,2))**2
168  dist3 = (f(i3,1)-f(arr,1))**2 + (f(i3,2)-f(arr,2))**2
169  IF (dist1.LT.dist) THEN
170  dist = dist1
171  elbest = ielem
172  igbest = i1
173  ilbest = 1
174  IF(i1.EQ.pt) ilprec = 1
175  IF(i2.EQ.pt) ilprec = 2
176  IF(i3.EQ.pt) ilprec = 3
177  ENDIF
178  IF (dist2.LT.dist) THEN
179  dist = dist2
180  elbest = ielem
181  igbest = i2
182  ilbest = 2
183  IF(i1.EQ.pt) ilprec = 1
184  IF(i2.EQ.pt) ilprec = 2
185  IF(i3.EQ.pt) ilprec = 3
186  ENDIF
187  IF(dist3.LT.dist) THEN
188  dist = dist3
189  elbest = ielem
190  igbest = i3
191  ilbest = 3
192  IF(i1.EQ.pt) ilprec = 1
193  IF(i2.EQ.pt) ilprec = 2
194  IF(i3.EQ.pt) ilprec = 3
195  ENDIF
196  ENDIF
197 !
198  END DO ! OVER ELEMENTS
199 !
200  IF (igbest.EQ.pt) THEN
201  WRITE(lu,*)'FLUSEC : ALGORITHM FAILED'
202  CALL plante(1)
203  stop
204  ELSE
205  pt = igbest
206  iseg = iseg + 1
207  IF (iseg.GT.nsemax) THEN
208  WRITE(lu,*) 'TOO MANY SEGMENTS IN A '
209  WRITE(lu,*) 'SECTION. INCREASE NSEMAX'
210  CALL plante(1)
211  stop
212  ENDIF
213  liste(iseg,1) = ikle(elbest,ilprec)
214  liste(iseg,2) = ikle(elbest,ilbest)
215  IF (igbest.NE.arr) GOTO 1010
216  ENDIF
217  chain(isec)%NSEG = iseg
218  ALLOCATE (chain(isec)%LISTE(chain(isec)%NSEG,3), stat=ierr)
219  CALL check_allocate(ierr, 'CHAIN(ISEC)%LISTE')
220  DO iseg=1,chain(isec)%NSEG
221  chain(isec)%LISTE(iseg,1)=liste(iseg,1)
222  chain(isec)%LISTE(iseg,2)=liste(iseg,2)
223  chain(isec)%LISTE(iseg,3)=-1 ! INITIALISE... FOR DEVEL
224  END DO
225  ENDDO ! OVER SECTIONS
226  DEALLOCATE (liste)
227 !
228  ALLOCATE (anpbeg(nbmaxnshare), stat=ierr)
229  CALL check_allocate(ierr, 'ANPBEG')
230  ALLOCATE (anpend(nbmaxnshare), stat=ierr)
231  CALL check_allocate(ierr, 'ANPEND')
232 !
233  DO isec=1,nsec
234  DO iseg=1,chain(isec)%NSEG
235 !
236  npbeg=0
237  DO m=1,nparts
238  IF(hash_table_get(knogl,chain(isec)%LISTE(iseg,1),m)>0)THEN
239  npbeg=npbeg+1
240  END IF
241  END DO
242 
243  npend=0
244  DO m=1,nparts
245  IF(hash_table_get(knogl,chain(isec)%LISTE(iseg,2),m)>0)THEN
246  npend=npend+1
247  END IF
248  END DO
249 !
250  IF (npbeg>nbmaxnshare .OR. npend>nbmaxnshare) THEN
251  WRITE(lu,*) 'NPBEG OR NPEND: ',npbeg,npend
252  WRITE(lu,*) 'ARE LARGER THAN NBMAXNSHARE: ',nbmaxnshare
253  CALL plante(1)
254  stop
255  ENDIF
256 !
257  ! THE NICE AND USUAL CASE WHEN BOTH SEGMENT ENDS
258  ! BELONG TO ONE SUBDOMAIN - ONLY ONE POSITION IN KNOGL
259  IF ( npbeg==1 .AND. npend==1) THEN
260  m=0
261  DO k=1,nparts
262  temp=hash_table_get(knogl,chain(isec)%LISTE(iseg,1),k)
263  IF(m<temp)THEN
264  im(:)=k
265  m=temp
266  END IF
267  END DO
268 
269  !IM(:) = MAXLOC (
270  & ! HASH_TABLE_GET(KNOGL,CHAIN(ISEC)%LISTE(ISEG,1),:) )
271  !IN(:) = MAXLOC (
272  & ! HASH_TABLE_GET(KNOGL,CHAIN(ISEC)%LISTE(ISEG,2),:) )
273  m=0
274  DO k=1,nparts
275  temp=hash_table_get(knogl,chain(isec)%LISTE(iseg,2),k)
276  IF(m<temp)THEN
277  in(:)=k
278  m=temp
279  END IF
280  END DO
281 
282  IF (im(1)==in(1)) THEN
283  chain(isec)%LISTE(iseg,3)=im(1)
284  ELSE ! THEY BELONG TO DIFFERENT SUBDOMAINS? HOW COME?
285  WRITE(lu,*) 'IMPOSSIBLE CASE (1) BY SECTIONS'
286  CALL plante(1)
287  stop
288  ENDIF
289  ! AT LEAST ONE OF THE TERMINAL NODES IS ON THE INTERFACE
290  ! TAKE THE LARGEST COMMON PARTITION NUMBER THEY BOTH BELONG TO
291  ELSE
292  IF (npbeg==1 .AND. npend>1) THEN ! THE SEGMENT'S END TOUCHES THE INTERFACE
293  m=0
294  DO k=1,nparts
295  temp=hash_table_get(knogl,chain(isec)%LISTE(iseg,1),k)
296  IF(m<temp)THEN
297  im(:)=k
298  m=temp
299  END IF
300  END DO
301  !IM(:) = MAXLOC (
302  & ! HASH_TABLE_GET(KNOGL,CHAIN(ISEC)%LISTE(ISEG,1),:) )
303  IF (hash_table_get(knogl,
304  & chain(isec)%LISTE(iseg,2),im(1))>0 ) THEN
305  chain(isec)%LISTE(iseg,3) = im(1)
306  ELSE
307  WRITE(lu,*) 'IMPOSSIBLE CASE (2) BY SECTIONS'
308  CALL plante(1)
309  stop
310  ENDIF
311  ELSE IF (npbeg>1 .AND. npend==1) THEN ! THE SEGMENT'S BEG. TOUCHES THE INTERFACE
312  m=0
313  DO k=1,nparts
314  temp=hash_table_get(knogl,chain(isec)%LISTE(iseg,2),k)
315  IF(m<temp)THEN
316  in(:)=k
317  m=temp
318  END IF
319  END DO
320  !IN(:) = MAXLOC (
321  & ! HASH_TABLE_GET(KNOGL,CHAIN(ISEC)%LISTE(ISEG,2),:) )
322  IF ( hash_table_get(knogl,
323  & chain(isec)%LISTE(iseg,1),in(1))>0 ) THEN
324  chain(isec)%LISTE(iseg,3) = in(1)
325  ELSE
326  WRITE(lu,*) 'IMPOSSIBLE CASE (3) BY SECTIONS'
327  CALL plante(1)
328  stop
329  ENDIF
330  ELSE ! I.E. (NPBEG>1 .AND. NPEND>1) - LIES JUST ON THE INTERFACE OR "A SHORTCUT"
331  anpbeg=0
332  anpend=0
333  i=0
334  DO n=1,nparts
335  IF ( hash_table_get(knogl,
336  & chain(isec)%LISTE(iseg,1),n)>0 ) THEN
337  i=i+1
338  anpbeg(i)=n
339  ENDIF
340  END DO
341  IF (i/=npbeg) WRITE(lu,*) 'OH! I/=NPBEG'
342  i=0
343  DO n=1,nparts
344  IF ( hash_table_get(knogl,
345  & chain(isec)%LISTE(iseg,2),n)>0 ) THEN
346  i=i+1
347  anpend(i)=n
348  ENDIF
349  END DO
350  IF (i/=npend) WRITE(lu,*) 'OH! I/=NPEND'
351 !
352  WRITE(lu,*) 'ANPBEG: ',anpbeg
353  WRITE(lu,*) 'ANPEND: ',anpend
354 !
355  found=.false.
356  DO i=npbeg,1,-1
357  DO j=npend,1,-1
358  IF (anpbeg(i)==anpend(j)) THEN
359  chain(isec)%LISTE(iseg,3) = anpbeg(i)
360  found=.true.
361  EXIT
362  ENDIF
363  END DO
364  IF (found) EXIT
365  END DO
366  IF (.NOT.found) THEN
367  WRITE(lu,*) 'BY SECTION WITH NODES: ',
368  & chain(isec)%LISTE(iseg,1),chain(isec)%LISTE(iseg,2)
369  WRITE(lu,*) 'IMPOSSIBLE CASE (4) BY SECTIONS'
370  CALL plante(1)
371  stop
372  ENDIF
373 
374  ENDIF
375  ENDIF
376  ENDDO
377  ENDDO
378 !
379  DEALLOCATE (anpbeg,anpend)
380 !
381 ! WRITE FILES
382 !
383  DO n=1,nparts
384  nameout=trim(namesec)//extens(nparts-1,n-1)
385 
386  WRITE(lu,*) 'WRITING: ', trim(nameout)
387 
388  OPEN(nsct,file=trim(nameout),form='FORMATTED',status='UNKNOWN')
389  rewind(nsct)
390  WRITE(nsct,*) '# SECTIONS PARTITIONED FOR ',
391  & extens(nparts-1,n-1)
392  WRITE(nsct,*) nsec, 1
393  DO isec=1,nsec
394  WRITE(nsct,*) trim(chain(isec)%DESCR)
395  i=count(chain(isec)%LISTE(:,3)==n)
396  WRITE(nsct,*) i
397  DO iseg=1,chain(isec)%NSEG
398  IF (chain(isec)%LISTE(iseg,3)==n) THEN
399  WRITE(nsct,*)
400  & hash_table_get(knogl,chain(isec)%LISTE(iseg,1),n),
401  & hash_table_get(knogl,chain(isec)%LISTE(iseg,2),n)
402  ENDIF
403  END DO
404  END DO
405  CLOSE(nsct)
406  END DO
407 !
408  END SUBROUTINE
409  END MODULE
mod_handle_sections::handle_sections
subroutine handle_sections(NAMESEC, NPARTS, NELEM, NDP, IKLE, NPOIN, F, KNOGL)
Definition: mod_handle_sections.f:21
declarations_special
Definition: declarations_special.F:3
mod_handle_sections
Definition: mod_handle_sections.f:3
declarations_partel
Definition: declarations_partel.f:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mod_hash_table
Definition: mod_hash_table.f:3
mod_hash_table::hash_table
Definition: mod_hash_table.f:37
mod_hash_table::hash_table_get
integer function hash_table_get(HT, X, Y)
Definition: mod_hash_table.f:341
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:38:16 for The TELEMAC-MASCARET system by doxygen 1.8.13