The TELEMAC-MASCARET system  trunk
divise.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE divise
3 ! *****************
4 !
5  &(x,y,ikle,ncolor,npoin,nelem,nelmax,nsom2,som2,indicp,indice,
6  & shp,elt,npmax,corr,level)
7 !
8 !***********************************************************************
9 ! PROGICIEL : STBTEL V7.2 J-M JANIN (LNH) 30 87 72 84
10 ! ORIGINE : TELEMAC
11 !***********************************************************************
12 !
13 ! FONCTION : DIVISION PAR 4 DE TOUTES LES MAILLES
14 !
15 !
16 !
17 !
18 !
19 !history J-M HERVOUET (JUBILADO)
20 !+ 24/10/2016
21 !+ V7P2
22 !+ Optimisation in the case of splitting of all elements.
23 !
24 !history A. LEROY (EDF LAB, LNHE)
25 !+ 27/10/2016
26 !+ V7P2
27 !+ Adding optional variables CORR and LEVEL for the automatic
28 !+ refinement procedure.
29 !
30 !-----------------------------------------------------------------------
31 ! ARGUMENTS
32 ! .________________.____.______________________________________________
33 ! | NOM |MODE| ROLE
34 ! |________________|____|______________________________________________
35 ! | X,Y |<-->| COORDONNEES DU MAILLAGE .
36 ! | IKLE |<-->| NUMEROS GLOBAUX DES NOEUDS DE CHAQUE ELEMENT
37 ! | NCOLOR |<-->| TABLEAU DES COULEURS DES POINTS DU MAILLAGE
38 ! | NPOIN |<-->| NOMBRE TOTAL DE NOEUDS DU MAILLAGE
39 ! | NELEM |<-->| NOMBRE TOTAL D'ELEMENTS DU MAILLAGE
40 ! | NELMAX | -->| DIMENSION EFFECTIVE DES TABLEAUX CONCERNANT
41 ! |________________|____|______________________________________________
42 ! MODE : -->(DONNEE NON MODIFIEE), <--(RESULTAT), <-->(DONNEE MODIFIEE)
43 !----------------------------------------------------------------------
44 ! APPELE PAR : STBTEL
45 !***********************************************************************
46 !
48 !
49  IMPLICIT NONE
50 !
51 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
52 !
53  INTEGER, INTENT(IN) :: NSOM2,NELMAX,NPMAX
54  INTEGER, INTENT(INOUT) :: ELT(npmax)
55  INTEGER, INTENT(INOUT) :: NPOIN,NELEM
56  INTEGER, INTENT(INOUT) :: IKLE(nelmax,*),INDICP(*),INDICE(*)
57  INTEGER, INTENT(INOUT) :: NCOLOR(*)
58  DOUBLE PRECISION, INTENT(IN) :: SOM2(10,2)
59  DOUBLE PRECISION, INTENT(INOUT) :: X(*),Y(*),SHP(npmax,3)
60 ! OPTIONAL ARGUMENTS
61  INTEGER, INTENT(INOUT), OPTIONAL :: CORR(nelmax,*)
62  INTEGER, INTENT(IN), OPTIONAL :: LEVEL
63 !
64 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
65 !
66  INTEGER IELEM,IPOIN,ISOM,NO1,NO2,NO3,NP1,NE1,NE2,NE3,I1,I2
67  INTEGER ISEG1, ISEG2, ISEG3
68  INTEGER :: NEW_ELEM, NEW_POIN
69  INTEGER :: ISEG,NO(3),NP(3),NEXT(3),I
70  DOUBLE PRECISION DX,DY
71  LOGICAL FOUND
72  parameter( next = (/ 2,3,1 /) )
73 !
74  INTEGER, ALLOCATABLE :: TAB(:,:,:)
75 !
76 !-----------------------------------------------------------------------
77 !
78  IF(nsom2.GE.3) THEN
79 !
80 !=======================================================================
81 ! LOOKING FOR ELEMENTS TO BE SPLIT INTO 2 OR 4, SPLITTING THEM
82 ! SPLITTING ELEMENTS INTO 4 OR 2 WITHIN A POLYGON
83 ! THIS WILL ADD DUPLICATED POINTS TO BE ELIMINATED THEN IN REMAIL
84 !=======================================================================
85 !
86  shp = -1
87  elt = -1
88  DO ipoin = 1,npoin
89  indicp(ipoin) = 1
90  ENDDO
91  DO isom = 1,nsom2
92  dx = som2(isom+1,1) - som2(isom,1)
93  dy = som2(isom+1,2) - som2(isom,2)
94  DO ipoin = 1,npoin
95  IF (dx*(y(ipoin)-som2(isom,2)).LT.
96  & dy*(x(ipoin)-som2(isom,1))) indicp(ipoin) = 0
97  ENDDO
98  ENDDO
99 !
100  DO ielem = 1,nelem
101  indice(ielem) = indicp(ikle(ielem,1))
102  & + 2*indicp(ikle(ielem,2))
103  & + 4*indicp(ikle(ielem,3))
104  ENDDO
105 !
106  new_poin = 0
107  new_elem = 0
108 !
109  ALLOCATE(tab(npoin,11,2))
110 !
111  DO i = 1,npoin
112  indicp(i) = 0
113  ENDDO
114 !
115  DO ielem = 1,nelem
116 !
117  IF(indice(ielem).EQ.7) THEN
118 !
119 ! ALREADY EXISTING 3 POINTS
120  no(1) = ikle(ielem,1)
121  no(2) = ikle(ielem,2)
122  no(3) = ikle(ielem,3)
123 ! 3 NEW POINTS
124  DO iseg=1,3
125  i1=min(no(iseg),no(next(iseg)))
126  i2=max(no(iseg),no(next(iseg)))
127  found=.false.
128  IF(indicp(i1).GT.0) THEN
129 ! LOOKING FOR AN ALREADY EXISTING SEGMENT STARTING WITH POINT I1
130  DO i=1,indicp(i1)
131  IF(i2.EQ.tab(i1,i,1)) THEN
132 ! FOUND!
133  np(iseg)=tab(i1,i,2)
134  found=.true.
135  EXIT
136  ENDIF
137  ENDDO
138  ENDIF
139 ! CASE OF A NEW SEGMENT
140  IF(.NOT.found) THEN
141  indicp(i1)=indicp(i1)+1
142 ! SECOND POINT OF NEW SEGMENT
143  tab(i1,indicp(i1),1)=i2
144 ! RANK OF MIDDLE POINT OF NEW SEGMENT
145  new_poin = new_poin + 1
146  tab(i1,indicp(i1),2)=npoin + new_poin
147  np(iseg)=npoin + new_poin
148  ENDIF
149  ENDDO
150 ! 3 NEW NUMBERS OF ELEMENTS
151  ne1 = nelem + new_elem + 1
152  ne2 = nelem + new_elem + 2
153  ne3 = nelem + new_elem + 3
154  new_elem = new_elem + 3
155 ! NEW ARRAYS (COORDINATES, COLOUR,...)
156  DO iseg=1,3
157  x(np(iseg)) = 0.5d0 * ( x(no(iseg)) + x(no(next(iseg))) )
158  y(np(iseg)) = 0.5d0 * ( y(no(iseg)) + y(no(next(iseg))) )
159  ncolor(np(iseg)) = ncolor(no(iseg))
160  ENDDO
161 ! OLD ELEMENT NUMBER TAKEN FOR FIRST NEW ONE
162  ikle(ielem,2) = np(1)
163  ikle(ielem,3) = np(3)
164 ! 3 OTHER NEW ELEMENTS
165  ikle(ne1,1) = np(1)
166  ikle(ne1,2) = no(2)
167  ikle(ne1,3) = np(2)
168  ikle(ne2,1) = np(3)
169  ikle(ne2,2) = np(2)
170  ikle(ne2,3) = no(3)
171  ikle(ne3,1) = np(2)
172  ikle(ne3,2) = np(3)
173  ikle(ne3,3) = np(1)
174 ! FILLING SHP AND ELT FOR THE 6 POINTS IN THE TRIANGLE
175 ! (THE GENERAL ALGORITHM IN INTERP IS AWFULLY LONG)
176 ! HERE POINTS WILL BE REACHED SEVERAL TIMES BUT WELL, NEVERMIND
177 ! IT IS MUCH FASTER
178  elt(no(1))=ielem
179  shp(no(1),1)=1.d0
180  shp(no(1),2)=0.d0
181  shp(no(1),3)=0.d0
182  elt(no(2))=ielem
183  shp(no(2),1)=0.d0
184  shp(no(2),2)=1.d0
185  shp(no(2),3)=0.d0
186  elt(no(3))=ielem
187  shp(no(3),1)=0.d0
188  shp(no(3),2)=0.d0
189  shp(no(3),3)=1.d0
190  elt(np(1))=ielem
191  shp(np(1),1)=0.5d0
192  shp(np(1),2)=0.5d0
193  shp(np(1),3)=0.d0
194  elt(np(2))=ielem
195  shp(np(2),1)=0.d0
196  shp(np(2),2)=0.5d0
197  shp(np(2),3)=0.5d0
198  elt(np(3))=ielem
199  shp(np(3),1)=0.5d0
200  shp(np(3),2)=0.d0
201  shp(np(3),3)=0.5d0
202 !
203  ELSEIF (indice(ielem).EQ.3.OR.
204  & indice(ielem).EQ.5.OR.
205  & indice(ielem).EQ.6) THEN
206 !
207  IF(indice(ielem).EQ.3) THEN
208  no1 = ikle(ielem,1)
209  no2 = ikle(ielem,2)
210  no3 = ikle(ielem,3)
211  iseg1 = 1
212  iseg2 = 2
213  iseg3 = 3
214  ELSEIF (indice(ielem).EQ.5) THEN
215  no1 = ikle(ielem,3)
216  no2 = ikle(ielem,1)
217  no3 = ikle(ielem,2)
218  iseg1 = 3
219  iseg2 = 1
220  iseg3 = 2
221  ELSE
222  no1 = ikle(ielem,2)
223  no2 = ikle(ielem,3)
224  no3 = ikle(ielem,1)
225  iseg1 = 2
226  iseg2 = 1
227  iseg3 = 3
228  ENDIF
229  i1=min(no1,no2)
230  i2=max(no1,no2)
231  found=.false.
232  IF(indicp(i1).GT.0) THEN
233 ! LOOKING FOR AN ALREADY EXISTING SEGMENT STARTING WITH POINT I1
234  DO i=1,indicp(i1)
235  IF(i2.EQ.tab(i1,i,1)) THEN
236 ! FOUND!
237  np1=tab(i1,i,2)
238  found=.true.
239  EXIT
240  ENDIF
241  ENDDO
242  ENDIF
243 ! CASE OF A NEW SEGMENT
244  IF(.NOT.found) THEN
245  indicp(i1)=indicp(i1)+1
246 ! SECOND POINT OF NEW SEGMENT
247  tab(i1,indicp(i1),1)=i2
248 ! RANK OF MIDDLE POINT OF NEW SEGMENT
249  new_poin = new_poin + 1
250  tab(i1,indicp(i1),2)=npoin + new_poin
251  np1=npoin+new_poin
252  ENDIF
253 !
254  ne1 = nelem + new_elem + 1
255  x(np1) = 0.5d0 * ( x(no1) + x(no2) )
256  y(np1) = 0.5d0 * ( y(no1) + y(no2) )
257  ncolor(np1) = ncolor(no1)
258  ikle(ielem,1) = no1
259  ikle(ielem,2) = np1
260  ikle(ielem,3) = no3
261  ikle( ne1,1) = no2
262  ikle( ne1,2) = no3
263  ikle( ne1,3) = np1
264 ! FILLING SHP AND ELT FOR THE 4 POINTS IN THE TRIANGLE
265 ! (THE GENERAL ALGORITHM IN INTERP IS AWFULLY LONG)
266 ! HERE POINTS WILL BE REACHED SEVERAL TIMES BUT WELL, NEVERMIND
267 ! IT IS MUCH FASTER
268  elt(no1)=ielem
269  shp(no1,iseg1)=1.d0
270  shp(no1,iseg2)=0.d0
271  shp(no1,iseg3)=0.d0
272  elt(no2)=ielem
273  shp(no2,iseg1)=0.d0
274  shp(no2,iseg2)=1.d0
275  shp(no2,iseg3)=0.d0
276  elt(no3)=ielem
277  shp(no3,iseg1)=0.d0
278  shp(no3,iseg2)=0.d0
279  shp(no3,iseg3)=1.d0
280  elt(np1)=ielem
281  shp(np1,iseg1)=0.5d0
282  shp(np1,iseg2)=0.5d0
283  shp(np1,iseg3)=0.d0
284 
285  new_elem = new_elem + 1
286  ELSE
287 
288 ! FILLING SHP AND ELT FOR THE 3 POINTS IN THE TRIANGLE
289 ! (THE GENERAL ALGORITHM IN INTERP IS AWFULLY LONG)
290 ! HERE POINTS WILL BE REACHED SEVERAL TIMES BUT WELL, NEVERMIND
291 ! IT IS MUCH FASTER
292  no1 = ikle(ielem,1)
293  no2 = ikle(ielem,2)
294  no3 = ikle(ielem,3)
295  elt(no1)=ielem
296  shp(no1,1)=1.d0
297  shp(no1,2)=0.d0
298  shp(no1,3)=0.d0
299  elt(no2)=ielem
300  shp(no2,1)=0.d0
301  shp(no2,2)=1.d0
302  shp(no2,3)=0.d0
303  elt(no3)=ielem
304  shp(no3,1)=0.d0
305  shp(no3,2)=0.d0
306  shp(no3,3)=1.d0
307 !
308  ENDIF
309 !
310  ENDDO !IELEM
311 !
312  npoin = npoin + new_poin
313  nelem = nelem + new_elem
314 
315  DEALLOCATE(tab)
316 !
317  ELSE
318 !
319 !=======================================================================
320 ! SPLITTING ELEMENTS INTO 4 IN ALL THE DOMAIN
321 ! DONE WITHOUT UNDUE DUPLICATIONS, SO NO NEED TO CALL REMAIL AFTER
322 !=======================================================================
323 !
324 ! INDICP : WILL BE THE NUMBER OF SEGMENTS TO WHICH A POINT BELONGS
325 ! AND IS ITS POINT OF SMALLEST RANK
326 ! SUPPOSED HERE NOT TO BE LARGER THAN 11
327 !
328 ! TAB(IPOIN,I,1) : SECOND POINT OF Ith SEGMENT STARTING WITH POINT
329 ! IPOIN
330 !
331 ! TAB(IPOIN,I,2) : MIDDLE POINT OF Ith SEGMENT STARTING WITH POINT
332 ! IPOIN
333 !
334  ALLOCATE(tab(npoin,11,2))
335 !
336  DO ipoin = 1,npoin
337  indicp(ipoin) = 0
338  ENDDO
339 !
340  DO ielem=1,nelem
341 ! ALREADY EXISTING 3 POINTS
342  no(1) = ikle(ielem,1)
343  no(2) = ikle(ielem,2)
344  no(3) = ikle(ielem,3)
345 ! 3 NEW POINTS
346  DO iseg=1,3
347  i1=min(no(iseg),no(next(iseg)))
348  i2=max(no(iseg),no(next(iseg)))
349  found=.false.
350  IF(indicp(i1).GT.0) THEN
351 ! LOOKING FOR AN ALREADY EXISTING SEGMENT STARTING WITH POINT I1
352  DO i=1,indicp(i1)
353  IF(i2.EQ.tab(i1,i,1)) THEN
354 ! FOUND!
355  np(iseg)=tab(i1,i,2)
356  found=.true.
357  EXIT
358  ENDIF
359  ENDDO
360  ENDIF
361 ! CASE OF A NEW SEGMENT
362  IF(.NOT.found) THEN
363  indicp(i1)=indicp(i1)+1
364 ! SECOND POINT OF NEW SEGMENT
365  tab(i1,indicp(i1),1)=i2
366 ! RANK OF MIDDLE POINT OF NEW SEGMENT
367  npoin=npoin+1
368  tab(i1,indicp(i1),2)=npoin
369  np(iseg)=npoin
370  ENDIF
371  ENDDO
372 ! 3 NEW NUMBERS OF ELEMENTS
373  ne1 = nelem + ielem
374  ne2 = 2*nelem + ielem
375  ne3 = 3*nelem + ielem
376 ! NEW ARRAYS (COORDINATES, COLOUR,...)
377  DO iseg=1,3
378  x(np(iseg)) = 0.5d0 * ( x(no(iseg)) + x(no(next(iseg))) )
379  y(np(iseg)) = 0.5d0 * ( y(no(iseg)) + y(no(next(iseg))) )
380  ncolor(np(iseg)) = ncolor(no(iseg))
381  ENDDO
382 ! OLD ELEMENT NUMBER TAKEN FOR FIRST NEW ONE
383  ikle(ielem,2) = np(1)
384  ikle(ielem,3) = np(3)
385 ! 3 OTHER NEW ELEMENTS
386  ikle(ne1,1) = np(1)
387  ikle(ne1,2) = no(2)
388  ikle(ne1,3) = np(2)
389  ikle(ne2,1) = np(3)
390  ikle(ne2,2) = np(2)
391  ikle(ne2,3) = no(3)
392  ikle(ne3,1) = np(2)
393  ikle(ne3,2) = np(3)
394  ikle(ne3,3) = np(1)
395 ! FILLING SHP AND ELT FOR THE 6 POINTS IN THE TRIANGLE
396 ! (THE GENERAL ALGORITHM IN INTERP IS AWFULLY LONG)
397 ! HERE POINTS WILL BE REACHED SEVERAL TIMES BUT WELL, NEVERMIND
398 ! IT IS MUCH FASTER
399  elt(no(1))=ielem
400  shp(no(1),1)=1.d0
401  shp(no(1),2)=0.d0
402  shp(no(1),3)=0.d0
403  elt(no(2))=ielem
404  shp(no(2),1)=0.d0
405  shp(no(2),2)=1.d0
406  shp(no(2),3)=0.d0
407  elt(no(3))=ielem
408  shp(no(3),1)=0.d0
409  shp(no(3),2)=0.d0
410  shp(no(3),3)=1.d0
411  elt(np(1))=ielem
412  shp(np(1),1)=0.5d0
413  shp(np(1),2)=0.5d0
414  shp(np(1),3)=0.d0
415  elt(np(2))=ielem
416  shp(np(2),1)=0.d0
417  shp(np(2),2)=0.5d0
418  shp(np(2),3)=0.5d0
419  elt(np(3))=ielem
420  shp(np(3),1)=0.5d0
421  shp(np(3),2)=0.d0
422  shp(np(3),3)=0.5d0
423  ENDDO
424 !
425 ! FOR AUTOMATIC REFINEMENT PROCEDURE
426 !
427  IF(PRESENT(corr).AND.PRESENT(level)) THEN
428  DO ielem=1,nelem
429  corr(ielem, level) = ielem
430  corr(ielem+ nelem,level) = ielem
431  corr(ielem+2*nelem,level) = ielem
432  corr(ielem+3*nelem,level) = ielem
433  ENDDO
434  ENDIF
435 !
436  nelem=4*nelem
437  DEALLOCATE(tab)
438 !
439  ENDIF
440 !
441 !=======================================================================
442 ! SORTIE LISTING
443 !=======================================================================
444 !
445  WRITE(lu,50) npoin,nelem
446 50 FORMAT(//,1x,'CUTTING ELEMENTS BY 4',
447  & /,1x,'---------------------',/,
448  & /,1x,'NEW NUMBER OF POINTS : ',i9,
449  & /,1x,'NEW NUMBER OF ELEMENTS : ',i9)
450 !
451  RETURN
452  END SUBROUTINE
subroutine divise(X, Y, IKLE, NCOLOR, NPOIN, NELEM, NELMAX, NSOM2, SOM2, INDICP, INDICE, SHP, ELT, NPMAX, CORR, LEVEL)
Definition: divise.f:8