1      SUBROUTINE ELIMIN
2C
3      INCLUDE 'com_coor.f'
4      INCLUDE 'com_faces.f'
5      INCLUDE 'com_options.f'
6C
7      INTEGER NODE(27),IPP(4,4,9),NUMF(6),IP8(4,6),NOD(8),IP827(8,8)
8     &       ,IP9(4,4),IP4(2,4),IP4T(3,4),IP6(3,4)
9     &       ,NFPR(6),IBID(28),NONO(8),IFAFA(6)
10      REAL*4  PRES(6),XBID(28)
11      EQUIVALENCE(IBID(2),NODE)
12      EQUIVALENCE(IBID,XBID)
13      EQUIVALENCE(NOD,NODE)
14      LOGICAL*4 ZTRI,EGAL
15      CHARACTER*128 CBID
16      DATA IPP / 25,18, 6,13,  25,17, 1, 9,
17     &           25, 9, 2,18,  25,13, 5,17,
18     &           23,14, 7,15,  23,16, 5,13,
19     &           23,13, 6,14,  23,15, 8,16,
20     &           22,19, 7,14,  22,18, 2,10,
21     &           22,10, 3,19,  22,14, 6,18,
22     &           21,11, 3,10,  21, 9, 1,12,
23     &           21,12, 4,11,  21,10, 2, 9,
24     &           26,15, 7,19,  26,11, 4,20,
25     &           26,20, 8,15,  26,19, 3,11,
26     &           24,16, 8,20,  24,12, 1,17,
27     &           24,17, 5,16,  24,20, 4,12,
28     &           27,23,14,22,  27,24,16,23,
29     &           27,21,12,24,  27,22,10,21,
30     &           27,26,15,23,  27,23,13,25,
31     &           27,25, 9,21,  27,21,11,26,
32     &           27,22,19,26,  27,26,20,24,
33     &           27,24,17,25,  27,25,18,22/
34      DATA NUMF / 3,6,2,5,4,1 /
35      DATA IP4 / 4,1 , 2,3 , 1,2 , 3,4 /
36      DATA IP4T / 1,3,2
37     &           ,2,3,4
38     &           ,1,2,4
39     &           ,1,4,3 /
40      DATA IP6  / 6,1,4
41     &           ,4,2,5
42     &           ,5,3,6
43     &           ,4,5,6 /
44      DATA IP8  / 5,8,4,1
45     &           ,2,3,7,6
46     &           ,5,1,2,6
47     &           ,4,8,7,3
48     &           ,1,4,3,2
49     &           ,8,5,6,7  /
50      DATA IP9  / 9,5,2,6
51     &           ,9,6,3,7
52     &           ,9,7,4,8
53     &           ,9,8,1,5  /
54      DATA IP827 /  1, 9,21,12,17,25,27,24
55     &            , 9, 2,10,21,25,18,22,27
56     &            ,21,10, 3,11,27,22,19,26
57     &            ,12,21,11, 4,24,27,26,20
58     &            ,17,25,27,24, 5,13,23,16
59     &            ,25,18,22,27,13, 6,14,23
60     &            ,27,22,19,26,23,14, 7,15
61     &            ,24,27,26,20,16,23,15, 8  /
62      DATA US3 / .33333333333333333333 /
63C
64      IF (ISTDOUT.EQ.0) THEN
65        IF (ILANG.EQ.0) THEN
66          PRINT*,'Fin de lecture des noeuds. Lecture des �l�ments...'
67        ELSE
68          PRINT*,'End of nodes reading. Beginning elements reading..'
69        ENDIF
70      ENDIF
71      DO I=1,NUMNP
72        ISD(I) = 0
73        INDEX(I) = 0
74      ENDDO
75      IF (NBCORN.NE.0) THEN
76        DO I=1,NUMNP
77          RAYON2(I) = 0.
78        ENDDO
79      ENDIF
80C
81      NDSMAX = 0
82      IFACEXT = 0
83      IVOL   = 0
84      NEL    = 0
85      ICPT   = 0
86      NFC    = 0
87      IGROUP = 0
88      NUMSD  = 0
89      NUMEL  = 0
90      NUMNP2 = NUMNP
91      MODTET = 0
92      ICOIN  = 0
93      IVIEUX = 0
94 1000 CONTINUE
95      IF (ICOURB.EQ.-5) THEN
96        IGROUP = 1
97        ICPT0 = 1
98        NELEM = NELDEL
99        NEL = NELEM
100        IF (ICOURXYZ.EQ.0) THEN
101          NDS = 3
102          DO N=1,NEL
103            NN = (N-1)*NDS+1
104            CALL QUALITY(NODEL(NN),N,NDS)
105          ENDDO
106        ELSEIF(ICOURXYZ.EQ.1) THEN
107          NDS = 2
108        ELSEIF(ICOURXYZ.EQ.2) THEN
109C
110C trait� dans elimin3
111C
112          ICPT = NF
113          GOTO 9999
114        ENDIF
115        NDSEL = NDS
116      ELSE
117        IF (IARCH.EQ.0) THEN
118          READ(IFAC,END=9999,ERR=9999) N,NTYP,NELEM,NDS
119        ELSE
120          CALL litrecbin(IFAC,NODE,LLL,0)
121          IF (LLL.LE.0) GOTO 9999
122          N     = NODE(1)
123          NTYP  = NODE(2)
124          NELEM = NODE(3)
125          NDS   = NODE(4)
126        ENDIF
127        IF (NELEM.LE.0) GOTO 1000
128        NDSEL = NDS
129        IF (NDS*(NEL+NELEM).GT.NEMAX) CALL TROPDEPOINTS(NEL+NELEM,NDS,0)
130CC maintenant on lit le no de sd directement
131CCC        IGROUP = IGROUP+1
132        IGROUP = N
133        ICPT0  = ICPT+1
134        NEL = NEL + NELEM
135      ENDIF
136      NMODP = MAX(1,NELEM/50)
137      IF (NMODP.EQ.1) NMODP = NFMAX
138C
139C element a 27 noeuds
140C
141      IF (NDS.GE.27) THEN
142        IF (ILANG.EQ.0) THEN
143          ELEMENTS = 'Hexa�dres 27 noeuds'
144        ELSE
145          ELEMENTS = 'Hexaedrons 27 nodes'
146        ENDIF
147        MODTET = 40
148        LELEM  = 19
149        IOPREF = 1
150        NBFAC  = 24
151        NFEXT  = 6
152        ICOIN  = 1
153        IOPT   = 0
154        DO N=1,NELEM
155          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
156     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
157          IF (IARCH.EQ.0) THEN
158            READ(IFAC) NE,(NODE(J),J=1,NDS)
159          ELSE
160            CALL litrecbin(IFAC,IBID,LLL,0)
161          ENDIF
162          IMAT = IREF3(NODE(27))
163          NUMEL = NUMEL+1
164          CALL ECRNOD(NODE,NODEL,NUMEL,NDS)
165          DO I=1,NDS
166            IF (ISD(NODE(I)).EQ.0) THEN
167              ISD(NODE(I)) = IGROUP
168            ELSEIF(ISD(NODE(I)).NE.IGROUP) THEN
169              ISD(NODE(I)) = -1
170            ENDIF
171          ENDDO
172          DO I=9,20
173            INDEX(NODE(I)) = INDEX(NODE(I))+1
174          ENDDO
175          DO I=21,26
176            INDEX(NODE(I)) = -1
177          ENDDO
178C
179C Decomposition de l'element a 27 noeuds en 24 faces exterieures
180C
181          DO IFCC = 1,NFEXT
182            IF (IELIMI.EQ.0) IOPT = 1
183            NUMFAC = NUMF(IFCC)
184            DO K=1,4
185              CALL TRIAGE(NODE(IPP(1,K,IFCC))
186     &                   ,NODE(IPP(2,K,IFCC))
187     &                   ,NODE(IPP(3,K,IFCC))
188     &                   ,NODE(IPP(4,K,IFCC))
189     &                   ,ICPT,ICPT0,NFC,IOPT,NUMEL,NUMFAC,NODE(27))
190              ISD2(ICPT) = IGROUP
191              ISD3(ICPT) = IMAT
192            ENDDO
193            IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC)
194          ENDDO
195C
196C Calcul des rayons de courbure et puissances dioptriques pour les yeux
197C
198          DO K=1,6
199            IF (ICOR(NODE(20+K)).NE.0) CALL RAYFAC(NODE,K)
200          ENDDO
201C
202C Decomposition en tetraedres pour isosurfaces
203C
204          DO K=1,8
205            DO KK=1,8
206              NONO(KK) = NODE(IP827(KK,K))
207            ENDDO
208            CALL TETRA(NONO)
209          ENDDO
210        ENDDO
211C
212C 4 noeuds....
213C
214      ELSEIF(NDS.EQ.4) THEN
215C
216C ... 2d --> quadrangles
217C
218        IF (I2D.NE.0.OR.NTYP.EQ.3) THEN
219          NBFAC = 1
220          IPREFC = -1
221          IOPREF = 0
222          IF (I2D.EQ.0) THEN
223            ELEMENTS = 'Quadrangles 3D Q1'
224          ELSE
225            ELEMENTS = 'Quadrangles 2D Q1'
226          ENDIF
227          LELEM = 17
228          DO N=1,NELEM
229          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
230     &           CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
231            ICPT = ICPT+1
232            IF (IARCH.EQ.0) THEN
233              IF (IVIEUX.EQ.0) THEN
234                READ(IFAC,END=3131,ERR=3131) NE,(NOD(J),J=1,NDS)
235     &                                    ,NPR,(NFPR(K),PRES(K),K=1,NPR)
236                GOTO 3132
237 3131           NPR = 0
238                IVIEUX = 1
239 3132           NUMEL = NUMEL+1
240              ELSE
241                READ(IFAC) NE,(NOD(J),J=1,NDS)
242                NUMEL = NUMEL+1
243              ENDIF
244            ELSE
245              CALL litrecbin(IFAC,IBID,LLL,0)
246cc              NE = IBID(1)
247              IF (LLL/4.GT.NDS+1) THEN
248                NPR = IBID(NDS+2)
249                IF (NPR.GT.0) THEN
250                  DO J=1,NPR
251                    NFPR(J) = IBID(NDS+1+2*J)
252                    PRES(J) = XBID(NDS+2*(J+1))
253                  ENDDO
254                ENDIF
255              ELSE
256                NPR = 0
257              ENDIF
258              NUMEL = NUMEL+1
259            ENDIF
260            IF (NPR.GT.0) THEN
261              DO J=1,NPR
262                IF (PRES(J).EQ.9999.) THEN
263                  IR = 5
264                ELSEIF(PRES(J).EQ.-9999.) THEN
265                  IR = 1
266                ELSE
267                  IR = MOD(NINT(PRES(J)),16)
268                  IF (IR.LT.0) IR = IR+16
269                ENDIF
270                IREF(NOD(IP4(1,NFPR(J)))) = IR
271                IREF(NOD(IP4(2,NFPR(J)))) = IR
272              ENDDO
273            ENDIF
274            DO J=1,NDS
275              ICLAS(J,ICPT) = NOD(J)
276            ENDDO
277            CALL ECRNOD(ICLAS(1,ICPT),NODEL,NUMEL,NDS)
278            DO J=1,NDS
279              NN = ICLAS(J,ICPT)
280              INDEX(NN) = INDEX(NN)+1
281              IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT
282              IF (ISD(NN).EQ.0) THEN
283                ISD(NN) = IGROUP
284              ELSEIF(ISD(NN).NE.IGROUP) THEN
285                ISD(NN) = -1
286              ENDIF
287            ENDDO
288            ICLAS(5,ICPT) = NUMEL
289            ICLAS(6,ICPT) = 0
290            NUMNP2 = NUMNP2+1
291            X(NUMNP2) = .25*( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT))
292     &                       +X(ICLAS(3,ICPT))+X(ICLAS(4,ICPT)) )
293            Y(NUMNP2) = .25*( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT))
294     &                       +Y(ICLAS(3,ICPT))+Y(ICLAS(4,ICPT)) )
295            Z(NUMNP2) = .25*( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT))
296     &                       +Z(ICLAS(3,ICPT))+Z(ICLAS(4,ICPT)) )
297            DEPX(NUMNP2) = .25*(DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT))
298     &                         +DEPX(ICLAS(3,ICPT))+DEPX(ICLAS(4,ICPT)))
299            DEPY(NUMNP2) = .25*(DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT))
300     &                         +DEPY(ICLAS(3,ICPT))+DEPY(ICLAS(4,ICPT)))
301            DEPZ(NUMNP2) = .25*(DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT))
302     &                         +DEPZ(ICLAS(3,ICPT))+DEPZ(ICLAS(4,ICPT)))
303            ICLAS(7,ICPT) = NUMNP2
304            ISD2(ICPT)    = IGROUP
305          ENDDO
306C
307C ... 3d --> tetraedres
308C
309        ELSE
310          IOPREF = 1
311          NBFAC = 4
312          MODTET = 1
313          IF (ILANG.EQ.0) THEN
314            ELEMENTS = 'Tetra�dres'
315            LELEM    = 10
316          ELSE
317            ELEMENTS = 'Tetraedrons'
318            LELEM    = 11
319          ENDIF
320          IF (IELIMI.EQ.0) THEN
321            IOPT = 4
322          ELSE
323            IOPT = 0
324          ENDIF
325          IF (NEL.EQ.NELEM) CALL CHECKFACEXT(CBID,LBID,IFACEXT)
326          IBORLU = 0
327 4436     CONTINUE
328          DO N=1,NELEM
329            IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
330     &           CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
331            IF (IFACEXT.GT.0.AND.IBORLU.EQ.0) THEN
332              READ(IFACEXT,*,END=4434,ERR=4434)
333     &             NB,NFA,(IFAFA(I),I=1,NFA)
334              IBORLU = 1
335              GOTO 4435
336 4434         IF (N.EQ.1) THEN
337                IFACEXT = -1
338                IF (ILANG.EQ.0) THEN
339                  PRINT*,'*** Erreur � la lecture de '//CBID(1:LBID)//
340     &                 ' - �l�ment',N
341                ELSE
342                  PRINT*,'*** Error reading '//CBID(1:LBID)//
343     &                 ' - element',N
344                ENDIF
345                NUMEL  = 0
346                NUMNP2 = NUMNP
347                GOTO 4436
348              ENDIF
349 4435         CONTINUE
350            ENDIF
351            IF (IARCH.EQ.0) THEN
352              READ(IFAC) NE,N1,N2,N3,N4
353            ELSE
354              CALL litrecbin(IFAC,IBID,LLL,0)
355cc              NE = IBID(1)
356              N1 = IBID(2)
357              N2 = IBID(3)
358              N3 = IBID(4)
359              N4 = IBID(5)
360            ENDIF
361            IVOL = IVOL+1
362            NTET(1,IVOL) = N1
363            NTET(2,IVOL) = N2
364            NTET(3,IVOL) = N3
365            NTET(4,IVOL) = N4
366            NUMEL = NUMEL+1
367            CALL ECRNOD(NTET(1,IVOL),NODEL,NUMEL,4)
368            NUMNP2 = NUMNP2+1
369            X(NUMNP2) = .25*( X(N1)+X(N2)+X(N3)+X(N4) )
370            Y(NUMNP2) = .25*( Y(N1)+Y(N2)+Y(N3)+Y(N4) )
371            Z(NUMNP2) = .25*( Z(N1)+Z(N2)+Z(N3)+Z(N4) )
372            DEPX(NUMNP2) = .25*( DEPX(N1)+DEPX(N2)+DEPX(N3)+DEPX(N4) )
373            DEPY(NUMNP2) = .25*( DEPY(N1)+DEPY(N2)+DEPY(N3)+DEPY(N4) )
374            DEPZ(NUMNP2) = .25*( DEPZ(N1)+DEPZ(N2)+DEPZ(N3)+DEPZ(N4) )
375            IF (IFACEXT.GT.0) THEN
376              IF (NB.EQ.NUMEL) THEN
377                IBORLU = 0
378                DO K=1,4
379                  IFLAG = 0
380                  DO KK=1,NFA
381                    IF (IFAFA(KK).EQ.K) IFLAG=KK
382                  ENDDO
383                  IF (IFLAG.NE.0) THEN
384                    IOPT = 0
385                    CALL ORDON2(NTET(IP4T(1,K),IVOL)
386     &                   ,NTET(IP4T(2,K),IVOL)
387     &                   ,NTET(IP4T(3,K),IVOL),M1,M2,M3)
388                    CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,K
389     &                   ,NUMNP2)
390                    ISD2(ICPT) = IGROUP
391                  ENDIF
392                ENDDO
393              ENDIF
394            ELSE
395              ICPT00 = ICPT
396              CALL ORDON2(N1,N4,N3,M1,M2,M3)
397              CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,4,NUMNP2)
398              ISD2(ICPT) = IGROUP
399              CALL ORDON2(N1,N2,N4,M1,M2,M3)
400              CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,3,NUMNP2)
401              ISD2(ICPT) = IGROUP
402              CALL ORDON2(N2,N3,N4,M1,M2,M3)
403              CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,2,NUMNP2)
404              ISD2(ICPT) = IGROUP
405              CALL ORDON2(N1,N3,N2,M1,M2,M3)
406              CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,1,NUMNP2)
407              ISD2(ICPT) = IGROUP
408              DO K=ICPT00,ICPT
409                DO J=1,NDS
410                  NN = ICLAS(J,K)
411                  INDEX(NN) = INDEX(NN)+1
412                  IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = K
413                  JJ = ICLAS(J,K)
414                  IF (ISD(JJ).EQ.0) THEN
415                    ISD(JJ) = IGROUP
416                  ELSEIF(ISD(JJ).NE.IGROUP) THEN
417                    ISD(JJ) = -1
418                  ENDIF
419                ENDDO
420              ENDDO
421              IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC)
422            ENDIF
423          ENDDO
424        ENDIF
425C
426C Triangles P2
427C
428      ELSEIF(NDS.EQ.6) THEN
429C
430C c'est comme ca qu'on dit qu'on fait du Q2 (ou P2). idiot
431C
432        MODTET = 40
433        IOPREF = 1
434        NBFAC  = 4
435        ICOIN  = 1
436        IF (I2D.EQ.0) THEN
437          ELEMENTS = 'Triangles 3D Q2'
438        ELSE
439          ELEMENTS = 'Triangles 2D Q2'
440        ENDIF
441        LELEM = 15
442        DO N=1,NELEM
443          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
444     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
445          IF (IARCH.EQ.0) THEN
446            READ(IFAC) NE,(NODE(J),J=1,NDS)
447          ELSE
448            CALL litrecbin(IFAC,IBID,LLL,0)
449cc            NE = IBID(1)
450          ENDIF
451          NUMEL = NUMEL+1
452          CALL ECRNOD(NODE,NODEL,NUMEL,NDS)
453          DO I=1,4
454            ICPT = ICPT+1
455            DO J=1,3
456              NN = NODE(IP6(J,I))
457              INDEX(NN) = INDEX(NN)+1
458              IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT
459              IF (ISD(NN).EQ.0) THEN
460                ISD(NN) = IGROUP
461              ELSEIF(ISD(NN).NE.IGROUP) THEN
462                ISD(NN) = -1
463              ENDIF
464              ICLAS(J,ICPT) = NN
465            ENDDO
466            ICLAS(4,ICPT) = ICLAS(3,ICPT)
467            ICLAS(5,ICPT) = NUMEL
468            ICLAS(6,ICPT) = 0
469            NUMNP2 = NUMNP2+1
470            X(NUMNP2) = X(ICLAS(1,ICPT))
471            Y(NUMNP2) = Y(ICLAS(1,ICPT))
472            Z(NUMNP2) = Z(ICLAS(1,ICPT))
473            DEPX(NUMNP2) = DEPX(ICLAS(1,ICPT))
474            DEPY(NUMNP2) = DEPY(ICLAS(1,ICPT))
475            DEPZ(NUMNP2) = DEPZ(ICLAS(1,ICPT))
476            ICLAS(7,ICPT) = NUMNP2
477            ISD2(ICPT)    = IGROUP
478          ENDDO
479        ENDDO
480C
481C Triangles P1 et segments
482C
483      ELSEIF(NDS.EQ.3.OR.NDS.EQ.2) THEN
484ccc?        IOPREF = 0
485        IOPREF = 1
486        NBFAC = 1
487        IF (NDS.EQ.3) THEN
488          IF (I2D.EQ.0) THEN
489            ELEMENTS = 'Triangles 3D P1'
490          ELSE
491            ELEMENTS = 'Triangles 2D P1'
492          ENDIF
493          LELEM = 15
494        ELSE
495          IF (I2D.EQ.0) THEN
496            ELEMENTS = 'Segments 3D'
497          ELSE
498            ELEMENTS = 'Segments 2D'
499          ENDIF
500          LELEM = 11
501        ENDIF
502        DO N=1,NELEM
503          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
504     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
505          ICPT = ICPT+1
506          IF (ICOURB.NE.-5) THEN
507            IF (IARCH.EQ.0) THEN
508              READ(IFAC) NE,(ICLAS(J,ICPT),J=1,NDS)
509            ELSE
510              CALL litrecbin(IFAC,IBID,LLL,0)
511cc              NE = IBID(1)
512              DO J=1,NDS
513                ICLAS(J,ICPT) = IBID(J+1)
514              ENDDO
515            ENDIF
516          ENDIF
517          NUMEL = NUMEL+1
518          CALL ECRNOD(ICLAS(1,ICPT),NODEL,NUMEL,NDS)
519          ICLAS(4,ICPT) = ICLAS(3,ICPT)
520          ICLAS(5,ICPT) = NUMEL
521          ICLAS(6,ICPT) = 0
522          NUMNP2 = NUMNP2+1
523          IF (NDS.EQ.3) THEN
524            X(NUMNP2) = ( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT))
525     &                  + X(ICLAS(3,ICPT)) )*US3
526            Y(NUMNP2) = ( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT))
527     &                  + Y(ICLAS(3,ICPT)) )*US3
528            Z(NUMNP2) = ( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT))
529     &                  + Z(ICLAS(3,ICPT)) )*US3
530            DEPX(NUMNP2) = ( DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT))
531     &                     + DEPX(ICLAS(3,ICPT)) )*US3
532            DEPY(NUMNP2) = ( DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT))
533     &                     + DEPY(ICLAS(3,ICPT)) )*US3
534            DEPZ(NUMNP2) = ( DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT))
535     &                     + DEPZ(ICLAS(3,ICPT)) )*US3
536          ELSE
537            X(NUMNP2) = ( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT)) )*.5
538            Y(NUMNP2) = ( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT)) )*.5
539            Z(NUMNP2) = ( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT)) )*.5
540            DEPX(NUMNP2) = (DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT)))*.5
541            DEPY(NUMNP2) = (DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT)))*.5
542            DEPZ(NUMNP2) = (DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT)))*.5
543            ICLAS(3,ICPT) = ICLAS(2,ICPT)
544          ENDIF
545          ICLAS(7,ICPT) = NUMNP2
546          ISD2(ICPT) = IGROUP
547          DO J=1,NDS
548            NN = ICLAS(J,ICPT)
549            INDEX(NN) = INDEX(NN)+1
550            IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT
551            JJ = ICLAS(J,ICPT)
552            IF (ISD(JJ).EQ.0) THEN
553              ISD(JJ) = IGROUP
554            ELSEIF(ISD(JJ).NE.IGROUP) THEN
555              ISD(JJ) = -1
556            ENDIF
557          ENDDO
558        ENDDO
559C
560C Briques a 8 noeuds
561C
562      ELSEIF(NDS.EQ.8) THEN
563        IOPREF = 0
564        NBFAC  = 6
565        MODTET = 5
566        IF (ILANG.EQ.0) THEN
567          ELEMENTS = 'Hexa�dres 8 noeuds'
568        ELSE
569          ELEMENTS = 'Hexaedrons 8 nodes'
570        ENDIF
571        LELEM    = 18
572        IF (IELIMI.EQ.0) THEN
573          IOPT = 8
574        ELSE
575          IOPT = 0
576        ENDIF
577        IF (NEL.EQ.NELEM) CALL CHECKFACEXT(CBID,LBID,IFACEXT)
578        IBORLU = 0
579 3436   CONTINUE
580        DO N=1,NELEM
581          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
582     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
583          IF (IFACEXT.GT.0.AND.IBORLU.EQ.0) THEN
584            READ(IFACEXT,*,END=3434,ERR=3434) NB,NFA,(IFAFA(I),I=1,NFA)
585            IBORLU = 1
586            GOTO 3435
587 3434       IF (N.EQ.1) THEN
588              IFACEXT = -1
589              IF (ILANG.EQ.0) THEN
590                PRINT*,'*** Erreur � la lecture de '//CBID(1:LBID)//
591     &               ' - �l�ment',N
592              ELSE
593                PRINT*,'*** Error reading '//CBID(1:LBID)//
594     &               ' - element',N
595              ENDIF
596              NUMEL  = 0
597              NUMNP2 = NUMNP
598              GOTO 3436
599            ENDIF
600 3435       CONTINUE
601          ENDIF
602          IF (IARCH.EQ.0) THEN
603            IF (IVIEUX.EQ.0) THEN
604              READ(IFAC,END=3232,ERR=3232) NE,(NODE(J),J=1,NDS)
605     &                                   ,NPR,(NFPR(J),PRES(J),J=1,NPR)
606              GOTO 3233
607 3232         NPR = 0
608              IVIEUX = 1
609 3233         NUMEL = NUMEL+1
610            ELSE
611              READ(IFAC) NE,(NODE(J),J=1,NDS)
612              NUMEL = NUMEL+1
613            ENDIF
614          ELSE
615            CALL litrecbin(IFAC,IBID,LLL,0)
616cc            NE = IBID(1)
617            IF (LLL/4.GT.NDS+1) THEN
618              NPR = IBID(NDS+2)
619              IF (NPR.GT.0) THEN
620                DO J=1,NPR
621                  NFPR(J) = IBID(NDS+1+2*J)
622                  PRES(J) = XBID(NDS+2*(J+1))
623                ENDDO
624              ENDIF
625            ELSE
626              NPR = 0
627            ENDIF
628            NUMEL = NUMEL+1
629          ENDIF
630          IF (NPR.GT.0) THEN
631            DO J=1,NPR
632              IF (PRES(J).EQ.9999.) THEN
633                IR = 5
634              ELSEIF(PRES(J).EQ.-9999.) THEN
635                IR = 1
636              ELSE
637                IR = MOD(NINT(PRES(J)),16)
638                IF (IR.LT.0) IR = IR+16
639              ENDIF
640              IREF(NODE(IP8(1,NFPR(J)))) = IR
641              IREF(NODE(IP8(2,NFPR(J)))) = IR
642            ENDDO
643          ENDIF
644          CALL ECRNOD(NODE,NODEL,NUMEL,NDS)
645          DO I=1,NDS
646            IF (ISD(NODE(I)).EQ.0) THEN
647              ISD(NODE(I)) = IGROUP
648            ELSEIF(ISD(NODE(I)).NE.IGROUP) THEN
649              ISD(NODE(I)) = -1
650            ENDIF
651          ENDDO
652          NUMNP2 = NUMNP2+1
653          X(NUMNP2) =
654     &      .125*( X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4))
655     &           + X(NODE(5))+X(NODE(6))+X(NODE(7))+X(NODE(8)) )
656          Y(NUMNP2) =
657     &      .125*( Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4))
658     &           + Y(NODE(5))+Y(NODE(6))+Y(NODE(7))+Y(NODE(8)) )
659          Z(NUMNP2) =
660     &      .125*( Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4))
661     &           + Z(NODE(5))+Z(NODE(6))+Z(NODE(7))+Z(NODE(8)) )
662          DEPX(NUMNP2) = .125*( DEPX(NODE(1))+DEPX(NODE(2))
663     &                        + DEPX(NODE(3))+DEPX(NODE(4))
664     &                        + DEPX(NODE(5))+DEPX(NODE(6))
665     &                        + DEPX(NODE(7))+DEPX(NODE(8)) )
666          DEPY(NUMNP2) = .125*( DEPY(NODE(1))+DEPY(NODE(2))
667     &                        + DEPY(NODE(3))+DEPY(NODE(4))
668     &                        + DEPY(NODE(5))+DEPY(NODE(6))
669     &                        + DEPY(NODE(7))+DEPY(NODE(8)) )
670          DEPZ(NUMNP2) = .125*( DEPZ(NODE(1))+DEPZ(NODE(2))
671     &                        + DEPZ(NODE(3))+DEPZ(NODE(4))
672     &                        + DEPZ(NODE(5))+DEPZ(NODE(6))
673     &                        + DEPZ(NODE(7))+DEPZ(NODE(8)) )
674          IF (IFACEXT.GT.0) THEN
675            IF (NB.EQ.NUMEL) THEN
676              IBORLU = 0
677              DO K=1,6
678                IFLAG = 0
679                DO KK=1,NFA
680                  IF (IFAFA(KK).EQ.K) IFLAG=KK
681                ENDDO
682                IF (IFLAG.NE.0) THEN
683                  IOPT = 0
684                  CALL ORDONN(NODE(IP8(1,K)),NODE(IP8(2,K))
685     &                       ,NODE(IP8(3,K)),NODE(IP8(4,K)),M1,M2,M3,M4)
686                  CALL TRIAGE(M1,M2,M3,M4,ICPT,ICPT0,NFC,IOPT,NUMEL,K
687     &                       ,NUMNP2)
688                  ISD2(ICPT) = IGROUP
689                ENDIF
690              ENDDO
691            ENDIF
692          ELSE
693            DO K=1,6
694              CALL ORDONN(NODE(IP8(1,K)),NODE(IP8(2,K))
695     &                   ,NODE(IP8(3,K)),NODE(IP8(4,K)),M1,M2,M3,M4)
696              CALL TRIAGE(M1,M2,M3,M4,ICPT,ICPT0,NFC,IOPT,NUMEL,K
697     &                   ,NUMNP2)
698              ISD2(ICPT) = IGROUP
699            ENDDO
700            IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC)
701          ENDIF
702          DO K=1,8
703            INDEX(NODE(K)) = INDEX(NODE(K))+1
704          ENDDO
705          CALL TETRA(NODE)
706        ENDDO
707C
708C 9 noeuds....
709C
710      ELSEIF(NDS.EQ.9) THEN
711C
712C c'est comme ca qu'on dit qu'on fait du Q2 (ou P2). idiot
713C
714        MODTET = 40
715        IOPREF = 1
716        NBFAC  = 4
717        ICOIN  = 1
718        IF (I2D.EQ.0) THEN
719          ELEMENTS = 'Quadrangles 3D Q2'
720        ELSE
721          ELEMENTS = 'Quadrangles 2D Q2'
722        ENDIF
723        LELEM = 17
724        DO N=1,NELEM
725          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
726     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
727          IF (IARCH.EQ.0) THEN
728            READ(IFAC) NE,(NODE(J),J=1,NDS)
729          ELSE
730            CALL litrecbin(IFAC,IBID,LLL,0)
731cc            NE = IBID(1)
732          ENDIF
733          IMAT = IREF3(NODE(9))
734          NUMEL = NUMEL+1
735          CALL ECRNOD(NODE,NODEL,NUMEL,NDS)
736          DO I=1,4
737            ICPT = ICPT+1
738            DO J=1,4
739              NN = NODE(IP9(J,I))
740              INDEX(NN) = INDEX(NN)+1
741              IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT
742              IF (ISD(NN).EQ.0) THEN
743                ISD(NN) = IGROUP
744              ELSEIF(ISD(NN).NE.IGROUP) THEN
745                ISD(NN) = -1
746              ENDIF
747              ICLAS(J,ICPT) = NN
748            ENDDO
749            ICLAS(5,ICPT) = NUMEL
750            ICLAS(6,ICPT) = 0
751            NUMNP2 = NUMNP2+1
752            X(NUMNP2) = X(ICLAS(1,ICPT))
753            Y(NUMNP2) = Y(ICLAS(1,ICPT))
754            Z(NUMNP2) = Z(ICLAS(1,ICPT))
755            DEPX(NUMNP2) = DEPX(ICLAS(1,ICPT))
756            DEPY(NUMNP2) = DEPY(ICLAS(1,ICPT))
757            DEPZ(NUMNP2) = DEPZ(ICLAS(1,ICPT))
758            ICLAS(7,ICPT) = NUMNP2
759            ISD2(ICPT)    = IGROUP
760            ISD3(ICPT)    = IMAT
761          ENDDO
762        ENDDO
763      ELSE
764        IF (ILANG.EQ.0) THEN
765          PRINT*,'***** El�ments �',NDS,' noeuds de type inconnu *****'
766        ELSE
767          PRINT*,'*****',NDS,'-nodes elements are unsupported *****'
768        ENDIF
769        STOP
770      ENDIF
771      IF (ISTDOUT.EQ.0) CALL ENDFLUSH
772      NDSMAX = MAX(NDSMAX,NDS)
773      IF (ICOURB.NE.-5) GOTO 1000
774C
775C Recopie les coordonnees des points des faces selectionnees
776C
777 9999 IF (ICOURB.NE.-5) CLOSE(IFAC)
778      IF (IFACEXT.GT.0) CLOSE(IFACEXT)
779      ZTRI = NDSMAX.EQ.2.OR.NDSMAX.EQ.3.OR.NDSMAX.EQ.6
780     &  .OR.(NDSMAX.EQ.4.AND.I2D.EQ.0.AND.NTYP.NE.3)
781      IF (ZTRI) THEN
782        NDSB = 3
783      ELSE
784        NDSB = 4
785      ENDIF
786      IF (INOPO.NE.0) THEN
787        IF (ISAUVEMESH.EQ.1) THEN
788          CALL EXEC('/bin/mv -vi '//NOM_FICH(1:LONG)//'* .')
789        ENDIF
790        CALL EXEC('/bin/rm -f '//NOM_FICH(1:LONG)//'*')
791      ENDIF
792      NFACE = 0
793      DO 30 N=1,ICPT
794        IF (ICLAS(1,N).NE.0) THEN
795          NFACE = NFACE + 1
796          ISD2(NFACE) = ISD2(N)
797          ISD3(NFACE) = ISD3(N)
798          NUMSD = MAX(NUMSD,ISD2(NFACE))
799          NNUMFA(NFACE) = ICLAS(5,N)
800          NRFAC(NFACE)  = ICLAS(6,N)
801          DO I=1,NDSB
802            NFAC(I,NFACE) = ICLAS(I,N)
803            XF(I,NFACE) = X(ICLAS(I,N))
804            YF(I,NFACE) = Y(ICLAS(I,N))
805            ZF(I,NFACE) = Z(ICLAS(I,N))
806          ENDDO
807          NFAC(NDSB+1,NFACE) = ICLAS(7,N)
808          XF(NDSB+1,NFACE) = X(ICLAS(7,N))
809          YF(NDSB+1,NFACE) = Y(ICLAS(7,N))
810          ZF(NDSB+1,NFACE) = Z(ICLAS(7,N))
811          IF (ZTRI) THEN
812            XMA = MAX(X(ICLAS(1,N)),X(ICLAS(2,N))
813     &               ,X(ICLAS(3,N)))
814            XMI = MIN(X(ICLAS(1,N)),X(ICLAS(2,N))
815     &               ,X(ICLAS(3,N)))
816            ISMIN = MIN( ISD(ICLAS(1,N)),ISD(ICLAS(2,N))
817     &                  ,ISD(ICLAS(3,N)) )
818            ISMAX = MAX( ISD(ICLAS(1,N)),ISD(ICLAS(2,N))
819     &                  ,ISD(ICLAS(3,N)) )
820          ELSE
821            XMA = MAX(X(ICLAS(1,N)),X(ICLAS(2,N))
822     &               ,X(ICLAS(3,N)),X(ICLAS(4,N)))
823            XMI = MIN(X(ICLAS(1,N)),X(ICLAS(2,N))
824     &               ,X(ICLAS(3,N)),X(ICLAS(4,N)))
825            ISMIN = MIN( ISD(ICLAS(1,N)),ISD(ICLAS(2,N))
826     &                  ,ISD(ICLAS(3,N)),ISD(ICLAS(4,N)) )
827            ISMAX = MAX( ISD(ICLAS(1,N)),ISD(ICLAS(2,N))
828     &                  ,ISD(ICLAS(3,N)),ISD(ICLAS(4,N)) )
829          ENDIF
830          IF (EGAL(XMI,XMA)) THEN
831            IF (EGAL(XMI,-XMED).AND.ISYM.EQ.4) THEN
832              IFPLAN(NFACE) = 3
833            ELSE
834              IFPLAN(NFACE) = 0
835            ENDIF
836          ELSE
837            IF (ZTRI) THEN
838              YMA = MAX(Y(ICLAS(1,N)),Y(ICLAS(2,N))
839     &                 ,Y(ICLAS(3,N)))
840              YMI = MIN(Y(ICLAS(1,N)),Y(ICLAS(2,N))
841     &                 ,Y(ICLAS(3,N)))
842            ELSE
843              YMA = MAX(Y(ICLAS(1,N)),Y(ICLAS(2,N))
844     &                 ,Y(ICLAS(3,N)),Y(ICLAS(4,N)))
845              YMI = MIN(Y(ICLAS(1,N)),Y(ICLAS(2,N))
846     &                 ,Y(ICLAS(3,N)),Y(ICLAS(4,N)))
847            ENDIF
848            IF (EGAL(YMI,YMA)) THEN
849              IF (EGAL(YMI,-YMED)) THEN
850                IFPLAN(NFACE) = 2
851              ELSE
852                IFPLAN(NFACE) = 0
853              ENDIF
854            ENDIF
855          ENDIF
856          IF (ISMIN.NE.ISMAX) THEN
857            ITOTO = ISD2(NFACE)
858            IP = 1000
859            IF (ISD(ICLAS(1,N)).LT.0.AND.ISD(ICLAS(2,N)).LT.0) THEN
860              ITOTO = ITOTO + IP
861              IP = IP*10
862            ENDIF
863            IF (ISD(ICLAS(2,N)).LT.0.AND.ISD(ICLAS(3,N)).LT.0) THEN
864              ITOTO = ITOTO + IP*2
865              IP = IP*10
866            ENDIF
867            IF (ZTRI) THEN
868              IF (ISD(ICLAS(3,N)).LT.0.AND.ISD(ICLAS(1,N)).LT.0) THEN
869                ITOTO = ITOTO + IP*3
870              ENDIF
871            ELSE
872              IF (ISD(ICLAS(3,N)).LT.0.AND.ISD(ICLAS(4,N)).LT.0) THEN
873                ITOTO = ITOTO + IP*3
874                IP = IP*10
875              ENDIF
876              IF (ISD(ICLAS(4,N)).LT.0.AND.ISD(ICLAS(1,N)).LT.0) THEN
877                ITOTO = ITOTO + IP*4
878              ENDIF
879            ENDIF
880            ISD2(NFACE) = ITOTO
881          ENDIF
882        ENDIF
883 30   CONTINUE
884      IGROUP = NUMSD
885C
886C Pourquoi c'etait commente ? 6/98
887C
888      IF (I2D.NE.0) NDS = NDSB
889C
890      IF (NDS.EQ.3.OR.NDS.EQ.2) IPREFC = -1
891      IF (IPREFC.NE.-1) IORIENT = 1
892      IF (NDS.EQ.2) IFC = 1
893      IF (ISTDOUT.EQ.0) THEN
894        IF (ILANG.EQ.0) THEN
895          IF (IVIEUX.NE.0)
896     &       PRINT*,'*** Attention !! Vieux format de maillage.........'
897          PRINT*,'Fin de lecture. Recherche de la fronti�re.........'
898        ELSE
899          IF (IVIEUX.NE.0)
900     &       PRINT*,'*** Warning !! Old-style mesh file................'
901          PRINT*,'End of reading. Computing the boundary............'
902        ENDIF
903      ENDIF
904      CALL MINDIAM(NFACE)
905      CALL FROETC(NTYP)
906      NUMBIS = NUMNP2
907      IF (ZTRI) NDS=3
908C
909C Determination des elements du bord
910C
911      DO I=1,NEL
912        IELBOR(I) = 0
913      ENDDO
914      DO I=1,NFACE
915        IF (IELBOR(NNUMFA(I)).EQ.0) THEN
916          IELBOR(NNUMFA(I)) = I
917        ELSE
918          IELBOR(NNUMFA(I)) = -I
919        ENDIF
920      ENDDO
921      NNN = 0
922      DO I=1,NEL
923        IF (IELBOR(I).NE.0) NNN = NNN+1
924      ENDDO
925      IF (I2D.EQ.0.AND.IPS2D.EQ.0) THEN
926        IF (ISTDOUT.EQ.0) THEN
927          IF (ILANG.EQ.0) THEN
928            PRINT*,NNN,' �l�ments de bord sur'
929     &       ,NEL,' (',REAL(NINT(1000.*REAL(NNN)/REAL(NEL)))*.1,' %)'
930          ELSE
931            PRINT*,NNN,' boundary elements among'
932     &       ,NEL,' (',REAL(NINT(1000.*REAL(NNN)/REAL(NEL)))*.1,' %)'
933          ENDIF
934        ENDIF
935      ENDIF
936C
937C Ecriture du fichier .facext
938C (8 noeuds et 4 noeuds tetra seulement pour l'instant)
939C
940      IF (IFACEXT.EQ.-1.AND.IELIMI.EQ.0) THEN
941        IF (CBID(1:5).NE.'/tmp/') THEN
942          CALL PREMIER_LIBRE(IFACEXT)
943          OPEN(IFACEXT,FILE=CBID(1:LBID),IOSTAT=IERR)
944          IF (IERR.EQ.0) THEN
945            DO I=1,NEL
946              ITAB(I) = 0
947            ENDDO
948            DO I=1,NFACE
949              ITAB(NNUMFA(I)) = ITAB(NNUMFA(I))+1
950            ENDDO
951            NN = 0
952            DO I=1,NEL
953              IF (ITAB(I).GT.0) THEN
954                DO J=1,ITAB(I)
955                  NN = NN+1
956                  IFAFA(J) = NRFAC(NN)
957                ENDDO
958                WRITE(IFACEXT,*) I,ITAB(I),(IFAFA(K),K=1,ITAB(I))
959              ENDIF
960            ENDDO
961            CLOSE(IFACEXT)
962          ENDIF
963        ENDIF
964      ENDIF
965C
966      END
967C=======================================================================
968      SUBROUTINE ECRNOD(NODE,NODEL,NUMEL,NDS)
969C
970      INTEGER NODE(*),NODEL(*)
971C
972      II = (NUMEL-1)*NDS
973      DO I=1,NDS
974        NODEL(II+I) = NODE(I)
975      ENDDO
976C
977C calcul des criteres de qualite des EF
978C
979      CALL QUALITY(NODE,NUMEL,NDS)
980      END
981C=======================================================================
982      SUBROUTINE QUALITY(NODE,N,ND)
983      INCLUDE 'com_coor.f'
984      INCLUDE 'com_faces.f'
985C
986      DIMENSION ANG(4)
987      INTEGER   NODE(*),ISUC3(3),ISUC4(4),IPRE4(4),IP8(4,6)
988      DATA ISUC3 / 2,3,1 /
989      DATA ISUC4 / 2,3,4,1 /
990      DATA IPRE4 / 4,1,2,3 /
991      DATA IP8   / 5,8,4,1,
992     &             2,3,7,6,
993     &             5,1,2,6,
994     &             4,8,7,3,
995     &             1,4,3,2,
996     &             8,5,6,7  /
997C
998      RMAX = 0.
999      IF (ND.EQ.3.OR.ND.EQ.6) THEN
1000C
1001C Triangles
1002C
1003        XC = (X(NODE(1))+X(NODE(2))+X(NODE(3)))/3.
1004        YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3)))/3.
1005        ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3)))/3.
1006        DO J=1,3
1007          K = ISUC3(J)
1008          L = ISUC3(K)
1009          ANG(J) = CALCANG(X,Y,Z,NODE,J,K,L)
1010          RR = SQRT((X(NODE(J))-XC)**2
1011     &             +(Y(NODE(J))-YC)**2
1012     &             +(Z(NODE(J))-ZC)**2)
1013          RMAX = AMAX1(RMAX,RR)
1014        ENDDO
1015        QUALITE(N,1) = AMIN1(ANG(1),ANG(2),ANG(3))
1016     &                /AMAX1(ANG(1),ANG(2),ANG(3))
1017      ELSEIF(ND.EQ.4.OR.ND.EQ.9) THEN
1018C
1019C quadrangles (4,9 noeuds) (a revoir pour les tetra)
1020C
1021        XC = (X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4)))*0.25
1022        YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4)))*0.25
1023        ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4)))*0.25
1024        DO J=1,4
1025          K = ISUC4(J)
1026          L = IPRE4(J)
1027          ANG(J) = CALCANG(X,Y,Z,NODE,J,K,L)
1028          RR = SQRT((X(NODE(J))-XC)**2
1029     &             +(Y(NODE(J))-YC)**2
1030     &             +(Z(NODE(J))-ZC)**2)
1031          RMAX = AMAX1(RMAX,RR)
1032        ENDDO
1033        QUALITE(N,1) = AMIN1(ANG(1),ANG(2),ANG(3),ANG(4))
1034     &                /AMAX1(ANG(1),ANG(2),ANG(3),ANG(4))
1035      ELSE
1036C
1037C Hexa 8, 27 noeuds
1038C
1039        ANGMIN = 5.
1040        ANGMAX = 0.
1041        XC = (X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4))
1042     &       +X(NODE(5))+X(NODE(6))+X(NODE(7))+X(NODE(8)))*0.125
1043        YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4))
1044     &       +Y(NODE(5))+Y(NODE(6))+Y(NODE(7))+Y(NODE(8)))*0.125
1045        ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4))
1046     &       +Z(NODE(5))+Z(NODE(6))+Z(NODE(7))+Z(NODE(8)))*0.125
1047        DO I=1,6
1048          DO II=1,4
1049            J = IP8(II,I)
1050            K = IP8(ISUC4(II),I)
1051            L = IP8(IPRE4(II),I)
1052            AA = CALCANG(X,Y,Z,NODE,J,K,L)
1053            ANGMAX = AMAX1(ANGMAX,AA)
1054            ANGMIN = AMIN1(ANGMIN,AA)
1055          ENDDO
1056        ENDDO
1057        QUALITE(N,1) = ANGMIN/ANGMAX
1058        DO J=1,8
1059          RR = SQRT((X(NODE(J))-XC)**2
1060     &             +(Y(NODE(J))-YC)**2
1061     &             +(Z(NODE(J))-ZC)**2)
1062          RMAX = AMAX1(RMAX,RR)
1063        ENDDO
1064      ENDIF
1065C
1066C Pas fini...
1067C
1068      QUALITE(N,2) = RMAX
1069C
1070      END
1071C=======================================================================
1072      REAL*4 FUNCTION CALCANG(X,Y,Z,NODE,J,K,L)
1073      INTEGER   NODE(*)
1074      DIMENSION X(*),Y(*),Z(*)
1075CC      DATA CQVSPI / 57.29577951 /
1076      DATA PI / 3.14159265358979  /
1077C
1078      JJ = NODE(J)
1079      KK = NODE(K)
1080      LL = NODE(L)
1081      U1 = X(KK)-X(JJ)
1082      V1 = Y(KK)-Y(JJ)
1083      W1 = Z(KK)-Z(JJ)
1084      U2 = X(LL)-X(JJ)
1085      V2 = Y(LL)-Y(JJ)
1086      W2 = Z(LL)-Z(JJ)
1087      XN1 = SQRT(U1**2+V1**2+W1**2)
1088      XN2 = SQRT(U2**2+V2**2+W2**2)
1089      IF (XN1*XN2.NE.0.) THEN
1090        ANG = ASIN(MIN(1.,SQRT((U1*V2-U2*V1)**2
1091     &                 +(V1*W2-V2*W1)**2
1092     &                 +(W1*U2-W2*U1)**2)/(XN1*XN2)))
1093        SCAL = U1*U2 + V1*V2 + W1*W2
1094        IF (SCAL.GE.0.) THEN
1095          CALCANG = ANG
1096        ELSE
1097          CALCANG = PI-ANG
1098        ENDIF
1099      ELSE
1100        CALCANG = -1.
1101      ENDIF
1102      END
1103