Télécharger la plaquette Cast3M


fr en

Télécharger raff3d.eso

Retour à la liste

Numérotation des lignes : 

raff3d

  1. C RAFF3D    SOURCE    CB215821  26/08/24    21:18:06     12622          
  2.       SUBROUTINE RAFF3D(IPT2,IPT3,ICPR,KARPOS,KARETE,KMILIE,MELVA2,
  3.      .NACREE,KARAF,IPT4,JPLANS,JPLAN3,JPLCOM,JNOEFA,IPT7,JFARAF,KARET2,
  4.      .XDEN)
  5. c entrees de raff3D:
  6. c           - IPT2:   maillage élémentaire à raffiner 
  7. c           - IPT3:   maillage de SURE précédant 
  8. c           - ICPR:   tableau de passage noeuds local/global
  9. c           - KARPOS: tableau du nb d'arretes par noeuds
  10. c           - KARETE: tableau des d'arretes KARETE(i,j)=k 
  11. c           - KMILIE: tableau des hanging nodes associés au arêtes
  12. c                     si ancienne relation 
  13. c                        KMILIE(i,j)=n : le noeud global n est situe
  14. c                        au milieu de l'arete [i,k]
  15. c                     si nouvelle relation
  16. c                        KMILIE(i,j)=-1 : il faut construire un noeud
  17. c                        au milieu de l'arete [i,k]
  18. c           - MELVA2: champ par elem valait 1 pour les élem à raffiner
  19. c                     et 0 pour les autres 
  20. c           - NACREE: nombre de noeuds à créer dans ipt2
  21. c           - KARAF:  nombre d'élément à raffiner dans ipt2
  22. c           - JPLANS : tableau des faces JPLAN(i,j)=k
  23. c           - JPLAN3 : tableau du nombre de relation de compat par face
  24. c           - JPLCOM : tableau du nb de faces par noeuds
  25. c           - JNOEFA : tableau des noeuds et éléments associé au faces
  26. c           - JFARAF : tableau des hanging nodes associés au faces pour
  27. c                      les anciennes relations.
  28. c           - XDEN :  densité aux noeuds en notation globale
  29. c           - KARET2: tableau du nombre d'éléments touchant les arêtes
  30. c sorties: 
  31. c           - KMILIE: tableau des hanging nodes associés au arêtes 
  32. c                     (en sortie: à partir des relations donnée en 
  33. c                     entrée et de celles crées
  34. c                     si anciene relation à garder
  35. c                       KMILIE(i,j)=n : le noeud global n est situe
  36. c                       au milieu de l'arete [i,k]
  37. c                     si anciene relation à supprimer
  38. c                       KMILIE(i,j)=0 : 
  39. c                     si nouvelle relation
  40. c                       KMILIE(i,j)=n : le noeud n à été construit au
  41. c                       milieu de l'arete [i,k]
  42. c           - JPLAN3 : tableau du nombre de relation de compat par face
  43. c                      mis a jour
  44. c           - JFARAF : tableau des hanging nodes associés au faces pour
  45. c                      mis a jour.
  46. c           - MELEME: Maillage élémentaire raffiné à partir de ipt2
  47. c                     sans relations
  48. c           - IPT7 :  Maillage de relations créé incomplet 
  49. c           - XDEN :  densité aux noeuds en notation globale
  50. c                     interpolé aux nouveaux noeuds.
  51.       IMPLICIT REAL*8(A-H,O-Z)
  52.       IMPLICIT INTEGER(I-N)
  53.  
  54. -INC PPARAM
  55. -INC CCOPTIO
  56. -INC SMCOORD
  57. -INC CCGEOME
  58. -INC SMELEME
  59. -INC SMCHPOI
  60. -INC SMMODEL
  61. -INC SMCHAML
  62. C
  63. C======================================================================
  64. C                           Declarations
  65. C======================================================================
  66.       SEGMENT ICPR(nbpts,2)
  67.       SEGMENT XDEN(nbpts)
  68.       SEGMENT KARETE(NBNDS,NCOL)
  69.       SEGMENT KARET2(NBNDS,NCOL)
  70.       SEGMENT KMILIE(NBNDS,NCOL)
  71.       SEGMENT KARPOS(NBNDS)
  72.       SEGMENT JPLANS(JPLA1,JPLA2)
  73.       SEGMENT JPLAN3(JPLA1,JPLA2)
  74.       SEGMENT JPLCOM(JPLA1)
  75.       SEGMENT JNOEFA(JNBFA,5)
  76.       SEGMENT JFARAF(JNBFA,LLLL)
  77.       SEGMENT NUMNOE(INUMNO)
  78.       SEGMENT IWORK2(JNBSOM)
  79.       SEGMENT IWORK1(IJKLMN)
  80.       SEGMENT IWORK4(JNBSOM)
  81.       SEGMENT IWORK3(IJKLMN)
  82.  
  83.       REAL*8 LON68,LON510,LON47,LON29
  84. C
  85. C======================================================================
  86. C                        Initialisations
  87. C======================================================================
  88.       SEGACT JPLANS,JPLCOM,JNOEFA*MOD,JPLAN3*MOD, IPT3
  89.       SEGACT IPT2,ICPR,KARPOS,KARETE,KMILIE*MOD,MELVA2,JFARAF*MOD
  90.       IJKLMN=4
  91.       ipt7=0
  92.       SEGINI IWORK1
  93.       IWORK1(1)=4
  94.       IWORK1(2)=8
  95.       IWORK1(3)=6
  96.       IWORK1(4)=10
  97. C      JNBFA=JNOEFA(/1)
  98. C      LLLL=2
  99. C      IF (IPT2.ITYPEL.EQ.15) LLLL=10
  100. C      IF ((IPT2.ITYPEL.EQ.17).OR.(IPT2.ITYPEL.EQ.24)) LLLL=6
  101. C      LLLL=10
  102. C      SEGINI JFARAF
  103.       NCOMPL=LNELM(1,(IPT2.ITYPEL-1)*2+2)
  104.       IF (NCOMPL.NE.0) GOTO 999
  105.       NBNN=NBNNE(LNELM(2,(IPT2.ITYPEL-1)*2+1))
  106.       INELM=LNELM(1,(IPT2.ITYPEL-1)*2+1)
  107. C
  108. C Creation du squelette du maillage resultat
  109.       NBELEM=IPT2.NUM(/2)-KARAF+INELM*KARAF
  110.       NBSOUS=0
  111.       NBREF=0
  112.       SEGINI IPT4
  113.       IPT4.ITYPEL=LNELM(2,(IPT2.ITYPEL-1)*2+1)
  114.       LPO2=LPOS2(IPT2.ITYPEL)
  115.       segact mcoord*mod
  116.       NBPT0=nbpts
  117.       NBPT1=nbpts
  118.       NBPTS=NBPT1+NACREE+KARAF*NBINTE(IPT2.ITYPEL)
  119.       SEGADJ MCOORD,XDEN
  120.       INUMNO=NBRAF(IPT2.ITYPEL)
  121.       SEGINI NUMNOE
  122.       LCOMP=1
  123. C
  124. C======================================================================
  125. C                   Phase de raffinement 3D
  126. C======================================================================
  127. C=======================================
  128. C A) Boucle sur les elements a raffiner !
  129. C=======================================
  130.       DO 6 IARAF=1,MELVA2.VELCHE(/2)
  131.         IF (MELVA2.VELCHE(1,IARAF).NE.1) THEN
  132.           DO 1 IJKL=1,IPT2.NUM(/1)
  133.             IPT4.NUM(IJKL,NBELEM)=IPT2.NUM(IJKL,IARAF)
  134.    1      CONTINUE
  135.           NBELEM=NBELEM-1
  136.           GOTO 6
  137.         ENDIF
  138.         JCOMPT=0
  139.         JPOS5=LPOS5(IPT2.ITYPEL)
  140. C
  141. C==========================================================
  142. C B) Boucle sur les noeuds a creer pour raffiner l'element !
  143. C==========================================================
  144.         DO 4 I=1,NBRAF(IPT2.ITYPEL)
  145.           JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
  146.           JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
  147.           JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
  148.           LTYPNO=JTYPNO(JPOS5-1+I)
  149. c          write(*,*) 'LTYPNO', LTYPNO
  150. C
  151. C-------------------------------
  152. C ** B.1 / On est sur une arete !
  153. C-------------------------------
  154.           IF (LTYPNO.EQ.0) THEN
  155.             NPTA=IPT2.NUM(LISNOE(JLISN),IARAF)
  156.             NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF)
  157.             NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  158.             NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  159.             DO 2 K=1,MAX(1,KARPOS(NMIN))
  160.               IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K
  161.    2        CONTINUE
  162.             IF (KMILIE(NMIN,NEXIST).GT.0) THEN
  163. c si le noeud existe deja on rempli numnoe avec Kmilie
  164.               NUMNOE(I)=KMILIE(NMIN,NEXIST)
  165.               JCOMPT=JCOMPT+JLONG
  166.               GOTO 4
  167.             ELSE
  168. c sinon on ajoute un noeud a la fin de numnoe et on le met dans  kmilie
  169.               NBPT1=NBPT1+1
  170.               NUMNOE(I)=NBPT1
  171.               KMILIE(NMIN,NEXIST)=NBPT1
  172.             ENDIF
  173.  
  174. c            IF ((NBPT1.eq.10330).or.(NBPT1.eq.10341)) THEN
  175. c               WRITE(*,*) 'NUMNOE(I)', NUMNOE(I) 
  176. c               WRITE(*,*) 'NPTA', NPTA, 'NPTB', NPTB
  177. c               WRITE(*,*) 'NMIN', NMIN, 'NMAX', NMAX ,'NEXI', NEXIST
  178. c               WRITE(*,*) 'KMILIE(NMIN,NEXIST)', KMILIE(NMIN,NEXIST)
  179. c               WRITE(*,*) 'KARETE(NMIN,NEXIST)', KARETE(NMIN,NEXIST)
  180. c            endif
  181. C
  182. C------------------------------
  183. C ** B.2 / On est sur une face !
  184. C------------------------------
  185. C Il faut identifier cette face (numero, sommets...)
  186.           ELSEIF ((LTYPNO.GT.0).AND.(LTYPNO.LT.7)) THEN
  187. C => a) On initialise toutes les donnees relatives a cette face
  188.             JTYPEL=IPT2.ITYPEL
  189.             JLTEL2=LTEL(1,JTYPEL)-1+LTYPNO
  190.             JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO
  191.             JLDEL1=LDEL(1,JLTEL2)
  192.             JTYFAC=IWORK1(JLDEL1)
  193.             JLDEL2=LDEL(2,JLTEL2)
  194.             JNBSOM=NBSOM(JTYFAC)
  195.             JSPOS=NSPOS(JTYFAC)
  196. C => b) On identifie les sommets de la face (n° global)
  197.             SEGINI IWORK2
  198.             DO 10 IAA=1,JNBSOM
  199.               NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF)
  200.               IWORK2(IAA)=NGLOBA
  201.   10        CONTINUE
  202. C => c) On classe ces sommets par ordre croissant (NPTA < NPTB < NPTC)
  203.             NPTA=nbpts+1
  204.             NPTB=NPTA+1
  205.             NPTC=NPTB+1
  206.             DO 11 ICC=1,JNBSOM
  207.               IF (IWORK2(ICC).LT.NPTA) THEN
  208.                 NPTC=NPTB
  209.                 NPTB=NPTA
  210.                 NPTA=IWORK2(ICC)
  211.               ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  212.                 NPTC=NPTB
  213.                 NPTB=IWORK2(ICC)
  214.               ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  215.                 NPTC=IWORK2(ICC)
  216.               ENDIF
  217.   11        CONTINUE
  218. C => d) On passe ces sommets en n° locale
  219.             NPTA1=NPTA
  220.             NPTB1=NPTB
  221.             NPTC1=NPTC
  222.  
  223.             NPTA2=ICPR(NPTA,1)
  224.             NPTB2=ICPR(NPTB,1)
  225.             NPTC2=ICPR(NPTC,1)
  226.             IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  227.               NPTA=NPTA2
  228.               NPTB=MIN(NPTB2,NPTC2)
  229.               NPTC=MAX(NPTB2,NPTC2)
  230.             ENDIF
  231.             IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  232.               NPTA=NPTB2
  233.               NPTB=MIN(NPTA2,NPTC2)
  234.               NPTC=MAX(NPTA2,NPTC2)
  235.             ENDIF
  236.             IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  237.               NPTA=NPTC2
  238.               NPTB=MIN(NPTA2,NPTB2)
  239.               NPTC=MAX(NPTA2,NPTB2)
  240.             ENDIF
  241. C => e) On cherche le numero de la face
  242.             NEXIS2=0
  243.             DO 12 IEE=1,JPLCOM(NPTA)
  244.               MTMP=JPLANS(NPTA,IEE)
  245.               JJ1=JNOEFA(MTMP,1)
  246.               JJ2=JNOEFA(MTMP,2)
  247.               JJ3=JNOEFA(MTMP,3)
  248.               IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  249.                 NEXIS2=IEE
  250.               ENDIF
  251.   12        CONTINUE
  252.             JNUMFA=JPLANS(NPTA,NEXIS2)
  253. C
  254. C---------------------------------
  255. C ** B.3 / Raffinement de la face !
  256. C---------------------------------
  257.             KTEST1=JPLAN3(NPTA,NEXIS2)
  258.             KTEST2=NBINTE(JTYFAC)
  259. C => a) Si la face est de type QUA4 (un seul noeud a creer, au milieu)
  260.             IF (JTYFAC.EQ.8) THEN
  261.               IF (KTEST1.LT.KTEST2) THEN
  262. c    a.1) si le noeud milieu n'existe pas on l'ajoute a la fin de numnoe
  263. c          et on rempli JFARAF et JPLAN3
  264.                 NBPT1=NBPT1+1
  265.                 NUMNOE(I)=NBPT1
  266.                 JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
  267.                 JFARAF(JNUMFA,1)=NBPT1
  268.               ELSEIF (KTEST1.EQ.KTEST2) THEN
  269. c    a.2) si le noeud milieu existe deja on rempli numnoe grace à JFARAF
  270. c         et on met celui-ci a zéro
  271.                 NUMNOE(I)=JFARAF(JNUMFA,1)
  272.                 JFARAF(JNUMFA,1)=0
  273.                 JCOMPT=JCOMPT+JLONG
  274.                 GOTO 4
  275.               ELSE
  276.                 WRITE(*,*) 'ERREUR'
  277.               ENDIF
  278. c              IF ((NBPT1.eq.10248).or.(NBPT1.eq.10276)) THEN
  279. c              IF ((NPTA.eq.235).or.(NPTA.eq.236).or.(NPTA.eq.237)) THEN
  280. c                WRITE(*,*) 'NUMNOE(I)', NUMNOE(I) 
  281. c                WRITE(*,*) 'NPTA', NPTA1, 'NPTB', NPTB1, 'NPTC', NPTC1
  282. c                WRITE(*,*) 'NPTA', NPTA, 'NPTB', NPTB, 'NPTC', NPTC
  283. c                WRITE(*,*) 'JNUMFA', JNUMFA, 'NEXI', NEXIS2
  284. c                WRITE(*,*) 'JPLAN3(NPTA,NEXIS2)', JPLAN3(NPTA,NEXIS2)
  285. c                WRITE(*,*) 'JFARAF(JNUMFA,1)', JFARAF(JNUMFA,1)
  286. c                WRITE(*,*) 'JNOEFA(JNUMFA,1)', JNOEFA(JNUMFA,1)
  287. c                WRITE(*,*) 'JNOEFA(JNUMFA,2)', JNOEFA(JNUMFA,2)
  288. c                WRITE(*,*) 'JNOEFA(JNUMFA,3)', JNOEFA(JNUMFA,3)
  289. c                WRITE(*,*) 'JNOEFA(JNUMFA,4)', JNOEFA(JNUMFA,4)
  290. c                WRITE(*,*) 'JNOEFA(JNUMFA,5)', JNOEFA(JNUMFA,5)
  291. c              endif
  292.             ENDIF
  293. C => b) Si la face n'est pas de type QUA4 (donc de type TRI6 ou QUA8)
  294.             IF (JTYFAC.NE.8) THEN
  295. C    b.1)   Si il manque des noeuds dans la face
  296. C           On ajoute un noeud à la fin de numnoe
  297.               IF (KTEST1.LT.KTEST2) THEN
  298.                 NBPT1=NBPT1+1
  299.                 NUMNOE(I)=NBPT1
  300.                 JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
  301.                 NEXIS3=0
  302.                 XCO2=0.25
  303.                 DO 13 KBB=1,JLONG
  304.                   XCO1=XCOEFF(JPOS1-1+KBB)
  305.                   IF (XCO1.EQ.XCO2) NEXIS3=KBB
  306.   13            CONTINUE
  307.                 JFARAF(JNUMFA,2*(KTEST1)+1)=NBPT1
  308.                 IF (NEXIS3.NE.0) THEN
  309. C    b.1.1)     si le point est sur une arête du triangle interieur 
  310. C               (pour un TRI6)
  311. C               si le point est sur une arête de la coie interieure 
  312. C               (pour un QUA8)
  313.  
  314. C               on renseigne la partie paire de JFARAF qui doit déja 
  315. C               exister ????
  316.                   NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
  317.                   JFARAF(JNUMFA,2*(KTEST1)+2)=NGLOB
  318.                 ELSE
  319. C    b.1.2)     si le point est le milieu d'un QUA8
  320.  
  321. C               on met la partie paire de JFARAF à zéro
  322.                   JFARAF(JNUMFA,2*(KTEST1)+2)=0
  323.                 ENDIF
  324.               ELSEIF (KTEST1.EQ.KTEST2) THEN
  325. C    b.2)     Si il ne manque pas de noeuds dans la face
  326.                 NEXIS3=0
  327.                 NEXIS4=0
  328.                 XCO2=0.25
  329.                 DO 14 KBB=1,JLONG
  330.                   XCO1=XCOEFF(JPOS1-1+KBB)
  331.                   IF (XCO1.EQ.XCO2) NEXIS3=KBB
  332.   14            CONTINUE
  333.                 IF (NEXIS3.NE.0) THEN
  334. C    b.2.1)     si le point est sur une arête du triangle interieur 
  335. C               (pour un TRI6)
  336. C               si le point est sur une arête de la croie interieure 
  337. C               (pour un QUA8)
  338.                   NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
  339.                   DO 15 KAA=2,JFARAF(/2),2
  340.                      IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA
  341.   15              CONTINUE
  342. c                 on rempli numnoe grace à JFARAF et on met la partie 
  343. c                 associée au noeud en question de celui-ci a 0
  344.                   NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1)
  345.                   JFARAF(JNUMFA,NEXIS4-1)=0
  346.                 ELSE
  347. C    b.2.2)     si le point est le milieu d'un QUA8
  348.                   DO 16 KAA=2,JFARAF(/2),2
  349.                      IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA
  350.   16              CONTINUE
  351. c                 on rempli numnoe grace à JFARAF et on met la partie 
  352. c                 associée au noeud en question de celui-ci a 0
  353.                   NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1)
  354.                   JFARAF(JNUMFA,NEXIS4-1)=0
  355.                 ENDIF
  356.                 JCOMPT=JCOMPT+JLONG
  357.                 GOTO 4
  358.               ENDIF
  359.             ENDIF
  360. C
  361. C------------------------------------------------------
  362. C ** B.4 / On est a l'interieur du volume de l'element !
  363. C------------------------------------------------------
  364.           ELSEIF (LTYPNO.EQ.7) THEN
  365.             NBPT1=NBPT1+1
  366.             NUMNOE(I)=NBPT1
  367.           ENDIF
  368. C allocation de plus de mémoire dans xcoord et xden si nécéssaire
  369.           IF (NBPT1.EQ.NBPTS) THEN
  370.             NBPTS=NBPTS+200
  371.             SEGADJ MCOORD,XDEN
  372.           ENDIF
  373. C
  374. C==============================
  375. C C) Creation du nouveau point !
  376. C==============================
  377. C On continue ici que lorsque l'on doit creer un nouveau point
  378.           XPT=0.
  379.           YPT=0.
  380.           ZPT=0.
  381.           XDEN1=0.D0
  382.           XDENMIN = 1000.
  383.           DO 3 J=1,JLONG
  384.             NGLOB=IPT2.NUM(LISNOE(JLISN-1+J),IARAF)
  385.             XINI=XCOOR((NGLOB-1)*(IDIM+1)+1)
  386.             YINI=XCOOR((NGLOB-1)*(IDIM+1)+2)
  387.             ZINI=XCOOR((NGLOB-1)*(IDIM+1)+3)
  388.             XPT=XPT+XINI*XCOEFF(JPOS1-1+J)
  389.             YPT=YPT+YINI*XCOEFF(JPOS1-1+J)
  390.             ZPT=ZPT+ZINI*XCOEFF(JPOS1-1+J)
  391.  
  392.             XDEN1=XDEN1+XDEN(NGLOB)*XCOEFF(JPOS1-1+J)
  393.             IF (XDEN(NGLOB).LT.XDENMIN) XDENMIN = XDEN(NGLOB) 
  394. c          IF (NBPT1.eq.10219) THEN
  395. c            write (*,*) 'nglob', nglob ,'nbpt1', nbpt1
  396. c            write (*,*) 'xpt ', xpt
  397. c            write (*,*) 'ypt ', ypt
  398. c            write (*,*) 'zpt ', zpt
  399. c            write (*,*) 'xden1 ', xden1
  400. c            write (*,*) 'xcoeff1 ', XCOEFF(JPOS1-1+J)
  401. c          endif
  402.    3      CONTINUE
  403.  
  404. c seuil de densité a xdenmin
  405.           IF (XDEN1.LT.XDENMIN) XDEN1 = XDENMIN
  406.  
  407.           XCOOR((NBPT1-1)*(IDIM+1)+1)=XPT
  408.           XCOOR((NBPT1-1)*(IDIM+1)+2)=YPT
  409.           XCOOR((NBPT1-1)*(IDIM+1)+3)=ZPT
  410.           XCOOR((NBPT1-1)*(IDIM+1)+4)=XDEN1
  411.           XDEN(NBPT1)=XDEN1
  412.           JCOMPT=JCOMPT+JLONG
  413.  
  414. C======================================================================
  415.    4    CONTINUE
  416. C======================================================================
  417.  
  418.         JPOS4=LPOS4(IPT2.ITYPEL)
  419. C
  420. C Cas des tetraedres
  421.         LONDIA=0
  422.         JTYPEL=IPT2.ITYPEL
  423.         IF (JTYPEL.EQ.23) THEN
  424.           LPT6=NUMNOE(6)
  425.           LPT5=NUMNOE(5)
  426.           LPT8=NUMNOE(8)
  427.           LPT10=NUMNOE(10)
  428.           XPT6=XCOOR((LPT6-1)*(IDIM+1)+1)
  429.           XPT5=XCOOR((LPT5-1)*(IDIM+1)+1)
  430.           XPT8=XCOOR((LPT8-1)*(IDIM+1)+1)
  431.           XPT10=XCOOR((LPT10-1)*(IDIM+1)+1)
  432.           YPT6=XCOOR((LPT6-1)*(IDIM+1)+2)
  433.           YPT5=XCOOR((LPT5-1)*(IDIM+1)+2)
  434.           YPT8=XCOOR((LPT8-1)*(IDIM+1)+2)
  435.           YPT10=XCOOR((LPT10-1)*(IDIM+1)+2)
  436.           ZPT6=XCOOR((LPT6-1)*(IDIM+1)+3)
  437.           ZPT5=XCOOR((LPT5-1)*(IDIM+1)+3)
  438.           ZPT8=XCOOR((LPT8-1)*(IDIM+1)+3)
  439.           ZPT10=XCOOR((LPT10-1)*(IDIM+1)+3)
  440.           LON68=SQRT((XPT8-XPT6)**2+(YPT8-YPT6)**2+(ZPT8-ZPT6)**2)
  441.           LON510=SQRT((XPT10-XPT5)**2+(YPT10-YPT5)**2+(ZPT10-ZPT5)**2)
  442.           IF (LON510.LT.LON68) LONDIA=4*4
  443.         ENDIF
  444.         IF (JTYPEL.EQ.24) THEN
  445.           LPT2=NUMNOE(2)
  446.           LPT4=NUMNOE(4)
  447.           LPT7=NUMNOE(7)
  448.           LPT9=NUMNOE(9)
  449.           XPT2=XCOOR((LPT2-1)*(IDIM+1)+1)
  450.           XPT4=XCOOR((LPT4-1)*(IDIM+1)+1)
  451.           XPT7=XCOOR((LPT7-1)*(IDIM+1)+1)
  452.           XPT9=XCOOR((LPT9-1)*(IDIM+1)+1)
  453.           YPT2=XCOOR((LPT2-1)*(IDIM+1)+2)
  454.           YPT4=XCOOR((LPT4-1)*(IDIM+1)+2)
  455.           YPT7=XCOOR((LPT7-1)*(IDIM+1)+2)
  456.           YPT9=XCOOR((LPT9-1)*(IDIM+1)+2)
  457.           ZPT2=XCOOR((LPT2-1)*(IDIM+1)+3)
  458.           ZPT4=XCOOR((LPT4-1)*(IDIM+1)+3)
  459.           ZPT7=XCOOR((LPT7-1)*(IDIM+1)+3)
  460.           ZPT9=XCOOR((LPT9-1)*(IDIM+1)+3)
  461.           LON47=SQRT((XPT7-XPT4)**2+(YPT7-YPT4)**2+(ZPT7-ZPT4)**2)
  462.           LON29=SQRT((XPT9-XPT2)**2+(YPT9-YPT2)**2+(ZPT9-ZPT2)**2)
  463.           IF (LON29.LT.LON47) LONDIA=4*10
  464.         ENDIF
  465. C
  466. C===================================
  467. C D) Creation des nouveaux elements !
  468. C===================================
  469. C On remplit la portion de IPT4 relative aux elements crees a partir
  470. C de la division de l'element IARAF (indice de boucle 1).
  471. C Cette portion de IPT4 contient les colonnes dont la valeur s'etend
  472. C de INELM*(LCOMP-1)+1 a INELM*LCOMP.
  473.         DO 1000 J=1,INELM
  474.         DO 5 I=1,IPT4.NUM(/1)
  475.           NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I)
  476.           IF (((JTYPEL.EQ.23).OR.(JTYPEL.EQ.24)).AND.(J.GT.4)) THEN
  477.             NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I+LONDIA)
  478.           ENDIF
  479.           IF (NTEMP.GT.NBNN) THEN
  480.             IPT4.NUM(I,INELM*(LCOMP-1)+J)=NUMNOE(NTEMP-NBNN)
  481.           ELSE
  482.             IPT4.NUM(I,INELM*(LCOMP-1)+J)=IPT2.NUM(NTEMP,IARAF)
  483.           ENDIF
  484.    5    CONTINUE
  485.  1000   CONTINUE
  486.         LCOMP=LCOMP+1
  487.    6  CONTINUE
  488. C
  489.       NBPT =NBPT1
  490. C      SEGADJ MCOORD,XDEN
  491. C
  492. C=======================================================================
  493. C              Preparation du maillage de relations
  494. C=======================================================================
  495.       SEGSUP KARET2
  496.       NBNDS=KARETE(/1)
  497.       NCOL=KARETE(/2)
  498.       SEGINI KARET2
  499. C==================================
  500. C A) Relations dues aux faces (3D) !
  501. C==================================
  502. C Tous les noeuds qui restent dans JFARAF sont a creer en tant que
  503. C relations de conformite ou des bords du domaine
  504.  
  505. C-----------------------------------------------------------------------
  506. C ** A.1 / Prise en compte des relation 
  507. C          d'arretes pour les arretes de chaque faces 
  508. C-----------------------------------------------------------------------
  509.  
  510. C-----------------------------------------------------------------------
  511. C** A.1.1/ faces raffinée
  512. C-----------------------------------------------------------------------
  513.  
  514. c      write (*,*) 'relations dues au faces'
  515.       JTYPEL=IPT2.ITYPEL
  516.       JLTEL1=LTEL(1,JTYPEL)
  517.       JLTEL2=LTEL(2,JTYPEL)
  518.       IF (JTYPEL.NE.14) GOTO 141
  519. c Boucle sur les elements de l'ancien maillage
  520.       DO 121 KELEM1=1,IPT2.NUM(/2)
  521. c Boucle sur les faces de l'ancien maillage
  522.         DO 122 IBB1=1,JLTEL1
  523.           JLDEL1=LDEL(1,JLTEL2-1+IBB1)
  524.           JTYFAC=IWORK1(JLDEL1)
  525.           JLDEL2=LDEL(2,JLTEL2-1+IBB1)
  526.           JNBSOM=NBSOM(JTYFAC)
  527.           JSPOS=NSPOS(JTYFAC)
  528.           SEGINI IWORK2
  529. C
  530. C 1/ trouver le numéro de la face
  531.           DO 123 IAA=1,JNBSOM
  532.               NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM1)
  533.               IWORK2(IAA)=NGLOBA
  534.  123      CONTINUE
  535. C
  536. C 1.1/ Classement des 3 sommets par ordre croissant de n° globale
  537.           NPTA=NBPT +1
  538.           NPTB=NPTA+1
  539.           NPTC=NPTB+1
  540.           DO 124 ICC=1,JNBSOM
  541.             IF (IWORK2(ICC).LT.NPTA) THEN
  542.               NPTC=NPTB
  543.               NPTB=NPTA
  544.               NPTA=IWORK2(ICC)
  545.             ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  546.               NPTC=NPTB
  547.               NPTB=IWORK2(ICC)
  548.             ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  549.               NPTC=IWORK2(ICC)
  550.             ENDIF
  551.  124      CONTINUE
  552. C
  553. C1.2/ Classement des 3 sommets par ordre croissant de n° locale
  554.           NPTA2=ICPR(NPTA,1)
  555.           NPTB2=ICPR(NPTB,1)
  556.           NPTC2=ICPR(NPTC,1)
  557.           IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  558.             NPTA=NPTA2
  559.             NPTB=MIN(NPTB2,NPTC2)
  560.             NPTC=MAX(NPTB2,NPTC2)
  561.           ENDIF
  562.           IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  563.             NPTA=NPTB2
  564.             NPTB=MIN(NPTA2,NPTC2)
  565.             NPTC=MAX(NPTA2,NPTC2)
  566.           ENDIF
  567.           IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  568.             NPTA=NPTC2
  569.             NPTB=MIN(NPTA2,NPTB2)
  570.             NPTC=MAX(NPTA2,NPTB2)
  571.           ENDIF
  572. C
  573. C 1.3/ Recherche du numero de la face
  574.           NEXIS2=0
  575.           DO 125 IEE=1,JPLCOM(NPTA)
  576.             MTMP=JPLANS(NPTA,IEE)
  577.             JJ1=JNOEFA(MTMP,1)
  578.             JJ2=JNOEFA(MTMP,2)
  579.             JJ3=JNOEFA(MTMP,3)
  580.             IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  581.               NEXIS2=IEE
  582.             ENDIF
  583.  125      CONTINUE
  584.  
  585.           JNUMF1=JPLANS(NPTA,NEXIS2)
  586.  
  587. C
  588. C 2/ Exclusion des faces non voulues 
  589.  
  590. C on exlu les faces où un noeud existait déja
  591.           IF (JFARAF(JNUMF1,1).EQ.0) GOTO 122
  592. C on exlu les faces de bords
  593.           IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 122
  594.  
  595. C 3/ Recherche des arretes de la faces pour y mettre des relations 
  596.  
  597.           IF (JTYPEL.NE.14) GOTO 137
  598.           NBAR1=0
  599.           DO 1001 ISS=1, (JNBSOM-1)
  600.           DO 138 JSS=(ISS+1), JNBSOM
  601.  
  602.             NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  603.             NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  604.  
  605.             DO 139 IRR=1,MAX(1,KARPOS(NMIN))                   
  606.               IF (KARETE(NMIN,IRR).EQ.NMAX) THEN
  607.                 NBAR1=NBAR1+1
  608.                 KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1
  609. c                write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  610. c                IF (JFARAF(JNUMF1,1).EQ.15748) THEN
  611. c                  write (*,*) ISS,JSS 
  612. c                  write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!'
  613. c                  write (*,*) IWORK2(ISS), IWORK2(JSS)
  614. c                  write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  615. c                  write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR) 
  616. c                ENDIF
  617.               ENDIF
  618.  139        CONTINUE
  619.  138      CONTINUE 
  620.  1001     CONTINUE
  621. c          WRITE (*,*) 'Nbarrete', NBAR1
  622. c          IF (NBAR1.NE.4) THEN
  623. c             write (*,*) IWORK2(1), IWORK2(2)
  624. c             write (*,*) IWORK2(3), IWORK2(4) 
  625. c          ENDIF
  626.  137      CONTINUE    
  627.  122    CONTINUE        
  628.  121  CONTINUE
  629.  
  630.  141  CONTINUE
  631.  
  632. C-----------------------------------------------------------------------
  633. C** A.1.2/ relations de faces conservées dans le maillage précédent
  634. C-----------------------------------------------------------------------
  635.  
  636.       SEGINI IPT5
  637.       DO 670 MIE=1,MAX(1,IPT3.LISOUS(/1))
  638.         IF (IPT3.LISOUS(/1).EQ.0) THEN
  639.           IPT5=IPT3
  640.         ELSE
  641.           IPT5=IPT3.LISOUS(MIE)
  642.         ENDIF
  643.         SEGACT IPT5
  644. c        write (*,*) 'ipt3 raff 3D' , ipt3, IPT3.ITYPEL, IPT3.ICOLOR(1)
  645. c        write (*,*) 'ipt5 raff 3D' , ipt5, IPT5.ITYPEL, IPT5.ICOLOR(1)
  646.         IF (IPT5.ITYPEL.NE.48) GOTO 670
  647.         IF ((IPT5.ICOLOR(1).NE.3)) GOTO 670
  648.         LIGNUM=IPT2.NUM(/1)
  649.         JNBSOM=4
  650.         SEGINI IWORK2
  651.  
  652. C
  653. C 1/ trouver le numéro de la face
  654. C 1.1/ Classement des 3 sommets par ordre croissant (num globale)
  655.         DO 671,MAA=1,IPT5.NUM(/2)
  656.           DO 673 MBB=1,JNBSOM
  657. c            write (*,*) 'MAA' , MAA, 'MBB', MBB
  658.             IWORK2(MBB)=IPT5.NUM(MBB+1,MAA)
  659.  
  660.  673      CONTINUE
  661.           NPTA=nbpt +1
  662.           NPTB=NPTA+1
  663.           NPTC=NPTB+1
  664.           DO 672 ICC=1,JNBSOM
  665.             IF (IWORK2(ICC).LT.NPTA) THEN
  666.               NPTC=NPTB
  667.               NPTB=NPTA
  668.               NPTA=IWORK2(ICC)
  669.             ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  670.               NPTC=NPTB
  671.               NPTB=IWORK2(ICC)
  672.             ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  673.               NPTC=IWORK2(ICC)
  674.             ENDIF
  675.  672      CONTINUE
  676. C
  677. C 1.2/ Classement des 3 sommets par ordre croissant (num locale)
  678.           NPTA1=NPTA
  679.           NPTB1=NPTB
  680.           NPTC1=NPTC
  681.           NPTA2=ICPR(NPTA,1)
  682.           NPTB2=ICPR(NPTB,1)
  683.           NPTC2=ICPR(NPTC,1)
  684.           IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  685.             NPTA=NPTA2
  686.             NPTB=MIN(NPTB2,NPTC2)
  687.             NPTC=MAX(NPTB2,NPTC2)
  688.           ENDIF
  689.           IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  690.             NPTA=NPTB2
  691.             NPTB=MIN(NPTA2,NPTC2)
  692.             NPTC=MAX(NPTA2,NPTC2)
  693.           ENDIF
  694.           IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  695.             NPTA=NPTC2
  696.             NPTB=MIN(NPTA2,NPTB2)
  697.             NPTC=MAX(NPTA2,NPTB2)
  698.           ENDIF
  699. C
  700. C 1.3/ Recherche du numero de la face
  701.           NEXIS2=0
  702.           DO 674 IEE=1,JPLCOM(NPTA)
  703.             MTMP=JPLANS(NPTA,IEE)
  704.             JJ1=JNOEFA(MTMP,1)
  705.             JJ2=JNOEFA(MTMP,2)
  706.             JJ3=JNOEFA(MTMP,3)
  707.             IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  708.               NEXIS2=IEE
  709.             ENDIF
  710.  674      CONTINUE
  711.           JNUMF1=JPLANS(NPTA,NEXIS2)
  712. c          IF (MAA.EQ.5) THEN
  713. c                  write (*,*) 'MAA', MAA, 'nexis2', nexis2 
  714. c                  write (*,*) IWORK2(1), IWORK2(2), IWORK2(3), IWORK2(4) 
  715. c                  write (*,*) 'JFARAF(JNUMF1,1)', JFARAF(JNUMF1,1)
  716. c                  write (*,*) 'JNOEFA(JNUMF1,5)',JNOEFA(JNUMF1,5) 
  717. c          ENDIF
  718. C
  719. C 2/ Exclusion des faces non voulues
  720.           IF (NEXIS2.EQ.0) THEN
  721. c            WRITE (*,*)'FACE NON TROUVE SURE:', MAA 
  722.             GOTO 671
  723.           ENDIF
  724. C on exlu les faces où un noeud existait déja
  725.           IF (JFARAF(JNUMF1,1).EQ.0) GOTO 671
  726.  
  727. C on exlu les faces de bords
  728.           IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 671
  729.           NBAR1=0
  730. C
  731. C 3/ Recherche des arretes de la faces pour y mettre des relations
  732.  
  733.           DO 1002 ISS=1, (JNBSOM-1)
  734.           DO 135 JSS=(ISS+1), JNBSOM
  735.  
  736.             NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  737.             NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  738.  
  739.             DO 136 IRR=1,MAX(1,KARPOS(NMIN))                   
  740.               IF (KARETE(NMIN,IRR).EQ.NMAX) THEN
  741.                 NBAR1=NBAR1+1
  742.                 KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1
  743. c                write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  744. c                IF (MAA.EQ.5) THEN
  745. c                  write (*,*) ISS,JSS 
  746. c                  write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!'
  747. c                  write (*,*) IWORK2(ISS), IWORK2(JSS)
  748. c                  write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  749. c                  write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR) 
  750. c                ENDIF
  751.               ENDIF
  752.  136        CONTINUE
  753.  135      CONTINUE 
  754.  1002     CONTINUE
  755.  
  756.  671    CONTINUE
  757.  
  758.  670  CONTINUE
  759.       SEGDES IPT5
  760. C---------------------------------
  761. C ** A.2 / Initialisation de IPT5 !
  762. C---------------------------------
  763.  
  764.  
  765.  
  766. C 1/ Comptage du nombre de noeuds soumis a des relations
  767.       NBRELA=0
  768.       DO 114 IHF=1,JFARAF(/1)
  769. C on exlu les faces de bords
  770.         IF (JNOEFA(IHF,5).EQ.0) GOTO 114
  771. C on exlu les faces où un noeud existait déja
  772.         IF (JFARAF(IHF,1).EQ.0) GOTO 114
  773.         DO 113 JF=1,JFARAF(/2),2
  774.           IF (JFARAF(IHF,JF).NE.0) NBRELA=NBRELA+1
  775.  113    CONTINUE
  776.  114  CONTINUE
  777. C
  778. C 2/ Creation de IPT5
  779.       IF (NBRELA.EQ.0) THEN
  780.         IPT5=0
  781.         GOTO 210
  782.       ENDIF
  783. C      IF (IPT2.ITYPEL.EQ.14) NBNN=4
  784. C      IF (IPT2.ITYPEL.EQ.15) NBNN=8
  785. C      IF (IPT2.ITYPEL.EQ.17) NBNN=5
  786. C      IF (IPT2.ITYPEL.EQ.24) NBNN=5
  787.       NBNN=10
  788.       NBELEM=NBRELA
  789.       NBSOUS=0
  790.       NBREF=0
  791.       SEGINI IPT5
  792.       IPT5.ITYPEL=48
  793. C
  794. C 3/ Renseignement des noeuds support des relations
  795.       NBRELA=0
  796.       DO 116 IPF=1,JFARAF(/1)
  797.         IF (JNOEFA(IPF,5).EQ.0) GOTO 116
  798.         IF (JFARAF(IPF,1).EQ.0) GOTO 116
  799.         DO 115 JF=1,JFARAF(/2),2
  800.           IF (JFARAF(IPF,JF).EQ.0) GOTO 115
  801.           NBRELA=NBRELA+1
  802.           IPT5.NUM(1,NBRELA)=JFARAF(IPF,JF)
  803.  115    CONTINUE
  804.  116  CONTINUE
  805. c      write (*,*) 'NBRELA FACE ', NBRELA
  806. C
  807. C-----------------------------------------------------
  808. C ** A.3 / Recherche des noeuds formant les relations !
  809. C-----------------------------------------------------
  810. C 1/ Boucle sur l'ensemble des noeuds a creer
  811.       DO 200 IARAF=1,MELVA2.VELCHE(/2)
  812.         IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 200
  813.         JCOMPT=0
  814.         JPOS5=LPOS5(IPT2.ITYPEL)
  815.         DO 190 I=1,NBRAF(IPT2.ITYPEL)
  816.           JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
  817.           JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
  818.           JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
  819.           LTYPNO=JTYPNO(JPOS5-1+I)
  820.           IF ((LTYPNO.EQ.0).OR.(LTYPNO.EQ.7)) GOTO 189
  821. C
  822. C 2/ Preparation pour trouver le noeud et la face en question
  823.           JTYPEL=IPT2.ITYPEL
  824.           JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO
  825.           JLDEL1=LDEL(1,JLTEL2)
  826.           JTYFAC=IWORK1(JLDEL1)
  827.           JLDEL2=LDEL(2,JLTEL2)
  828.           JNBSOM=NBSOM(JTYFAC)
  829.           JSPOS=NSPOS(JTYFAC)
  830.           SEGINI IWORK2
  831. C
  832. C 3/ Classement des 3 sommets par ordre croissant de n° globale
  833.           DO 100 IAA=1,JNBSOM
  834.             NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF)
  835.             IWORK2(IAA)=NGLOBA
  836.  100      CONTINUE
  837.           NPTA=nbpt +1
  838.           NPTB=NPTA+1
  839.           NPTC=NPTB+1
  840.           DO 110 ICC=1,JNBSOM
  841.             IF (IWORK2(ICC).LT.NPTA) THEN
  842.               NPTC=NPTB
  843.               NPTB=NPTA
  844.               NPTA=IWORK2(ICC)
  845.             ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  846.               NPTC=NPTB
  847.               NPTB=IWORK2(ICC)
  848.             ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  849.               NPTC=IWORK2(ICC)
  850.             ENDIF
  851.  110      CONTINUE
  852. C
  853. C 4/ Classement des 3 sommets par ordre croissant de n° locale
  854.           NPTA2=ICPR(NPTA,1)
  855.           NPTB2=ICPR(NPTB,1)
  856.           NPTC2=ICPR(NPTC,1)
  857.           IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  858.             NPTA=NPTA2
  859.             NPTB=MIN(NPTB2,NPTC2)
  860.             NPTC=MAX(NPTB2,NPTC2)
  861.           ENDIF
  862.           IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  863.             NPTA=NPTB2
  864.             NPTB=MIN(NPTA2,NPTC2)
  865.             NPTC=MAX(NPTA2,NPTC2)
  866.           ENDIF
  867.           IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  868.             NPTA=NPTC2
  869.             NPTB=MIN(NPTA2,NPTB2)
  870.             NPTC=MAX(NPTA2,NPTB2)
  871.           ENDIF
  872. C
  873. C 5/ Recherche du numero de la face
  874.           NEXIS2=0
  875.           DO 120 IEE=1,JPLCOM(NPTA)
  876.             MTMP=JPLANS(NPTA,IEE)
  877.             JJ1=JNOEFA(MTMP,1)
  878.             JJ2=JNOEFA(MTMP,2)
  879.             JJ3=JNOEFA(MTMP,3)
  880.             IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  881.               NEXIS2=IEE
  882.             ENDIF
  883.  120      CONTINUE
  884.           JNUMFA=JPLANS(NPTA,NEXIS2)
  885.  
  886. C
  887. C 6/ Recherche du numero global du point
  888.           IF (JNOEFA(JNUMFA,5).EQ.0) GOTO 189
  889.           IF (JTYFAC.EQ.8) INOEGL=JFARAF(JNUMFA,1)
  890.           IF (JTYFAC.NE.8) THEN
  891.             NEXIS3=0
  892.             NEXIS4=0
  893.             XCO2=0.25
  894.             DO 140 KBB=1,JLONG
  895.               XCO1=XCOEFF(JPOS1-1+KBB)
  896.               IF (XCO1.EQ.XCO2) NEXIS3=KBB
  897.  140        CONTINUE
  898.             IF (NEXIS3.NE.0) THEN
  899.               NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
  900.               DO 150 KAA=2,JFARAF(/2),2
  901.                 IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA
  902.  150          CONTINUE
  903.               INOEGL=JFARAF(JNUMFA,NEXIS4-1)
  904.             ELSE
  905.               DO 160 KAA=2,JFARAF(/2),2
  906.                 IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA
  907.  160          CONTINUE
  908.               INOEGL=JFARAF(JNUMFA,NEXIS4-1)
  909.             ENDIF
  910.           ENDIF
  911.  
  912. C
  913. C------------------------------
  914. C ** A.4 / Remplissage de IPT5 !
  915. C------------------------------
  916. C 1/ Recherche de la position du point dans IPT5
  917.           NEXIS5=0
  918.           DO 170 IGG=1,IPT5.NUM(/2)
  919.             IF (INOEGL.EQ.IPT5.NUM(1,IGG)) NEXIS5=IGG
  920.  170      CONTINUE
  921.           IF (NEXIS5.EQ.0) GOTO 189
  922.           IF (INOEGL.EQ.139) THEN
  923. c            write (*,*) INOEGL, IWORK2(1), IWORK2(2), IWORK2(3)  
  924.           ENDIF
  925.  
  926. C
  927. C 2/ Renseignement des points formant les relations
  928.           DO 180 IHH=1,JLONG
  929.             IPT5.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF)
  930.  180      CONTINUE
  931.           IF (JLONG.EQ.4) IPT5.NUM(10,NEXIS5)=3
  932.           IF (JLONG.EQ.5) IPT5.NUM(10,NEXIS5)=4
  933.           IF (JLONG.EQ.8) THEN
  934.             IF (JPOS1.EQ.16) IPT5.NUM(10,NEXIS5)=5
  935.             IF (JPOS1.EQ.24) IPT5.NUM(10,NEXIS5)=6
  936.           ENDIF
  937.  189      CONTINUE
  938.         JCOMPT=JCOMPT+JLONG
  939.  190    CONTINUE
  940.  200  CONTINUE
  941.  210  CONTINUE
  942. C
  943. C===================================
  944. C B) Relations dues aux aretes (2D) !
  945. C===================================
  946. C On cree un maillage IPT6 contenant tous les noeuds soumis a des
  947. C relations
  948. C
  949.  
  950. c        ILPL=LPL(IPT4.ITYPEL)
  951. c        ILPT=LPT(IPT4.ITYPEL)
  952. c        DO 317 J=1,IPT4.NUM(/2)
  953. c        DO 317 K=1,ILPL*2-1,2
  954. c          NPTA=IPT4.NUM(KSEGM(ILPT+K-1),J)
  955. c          NPTB=IPT4.NUM(KSEGM(ILPT+K),J)
  956. c
  957. c          if ((NPTB.eq.129).and.(NPTA.eq.9087)) then 
  958. c            write(*,*) '!ici!',NBPT0 
  959. c          endif
  960. c          IF((NPTA.GT.NBPT0).OR.(NPTB.GT.NBPT0)) THEN
  961. c             GOTO 317
  962. c          ENDIF
  963. c          NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  964. c          NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  965. c          if ((NMIN.eq.235).and.(NMAX.eq.236)) then
  966. c              write(*,*) 'nmin', nmin, 'nmax', nmax
  967. c              write(*,*) 'NPTA',NPTA, 'NPTB', NPTB
  968. c          endif
  969. c          NEXIST=0
  970. c          DO 316 I=1,MAX(1,KARPOS(NMIN))
  971. c            IF (KARETE(NMIN,I).EQ.NMAX) THEN
  972. c              KARET2(NMIN,I)=KARET2(NMIN,I)+1
  973. c
  974. c              IF (KMILIE(NMIN,I).gt.0) then
  975. c prise en compte des arretes multiples sur le nouveau maillage
  976. c dans l'exemple ci dessous où les noeuds C et D sont des hanging nodes
  977. c L'arette [C B] n'a pas encore été prise en compte dans karet2 c'est 
  978. c ce qui est fait ici.
  979. c
  980. c |                   |
  981. c |                   |
  982. c |A        C    D    |B
  983. c x---------X----X----X
  984. c |         |    |    |
  985. c |         |    |    |
  986.  
  987. c1/ Test sur la premiere arete : NMIN-KMILIE(NMIN,I)
  988. c                NMILI = ICPR(KMILIE(NMIN,I),1)
  989. c                NEXIS2=0
  990. c                NMIN2=MIN(NMIN,NMILI)
  991. c                NMAX2=MAX(NMIN,NMILI)
  992. c                DO 318 I2=1,KMILIE(/2)
  993. c                  IF (KMILIE(NMIN2,I2).GT.0) THEN
  994. c                    KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1
  995. c                  ENDIF
  996. c 318            CONTINUE
  997. C
  998. C 2/ Test sur la deuxieme arete : NMAX-NPTC
  999. c                NEXIS2=0
  1000. c                NMIN2=MIN(NMAX,NMILI)
  1001. c                NMAX2=MAX(NMAX,NMILI)
  1002. c                DO 319 I2=1,KMILIE(/2)
  1003. c                  IF (KMILIE(NMIN2,I2).GT.0) THEN
  1004. c                    KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1
  1005. c                  ENDIF
  1006. c 319            CONTINUE
  1007. c              ENDIF
  1008. c            ENDIF
  1009. c 316      CONTINUE
  1010. c 317    CONTINUE
  1011.  
  1012. C 1/ Comptage du nombre de noeuds soumis a des relations
  1013.       NBELEM=0
  1014.       DO 1003 J=1,KMILIE(/2)
  1015.       DO 27 I=1,KMILIE(/1)
  1016.  
  1017. c        IF (KMILIE(I,J).eq.10330) write(*,*) '!kmilie', I, J, '=10330!'
  1018. c        IF (KMILIE(I,J).eq.10330) write(*,*) 'KARET2(I,J)', KARET2(I,J)
  1019. c        IF (KMILIE(I,J).eq.10341) write(*,*) '!kmilie', I, J, '=10341!'
  1020. c        IF (KMILIE(I,J).eq.10341) write(*,*) 'KARET2(I,J)', KARET2(I,J)
  1021.         IF (KARET2(I,J).EQ.0) GOTO 27
  1022.         IF (KMILIE(I,J).GT.0) NBELEM=NBELEM+1
  1023.  
  1024.   27  CONTINUE
  1025.  1003 CONTINUE
  1026. c      write (*,*) 'NBELEM', NBELEM
  1027. C
  1028. C 2/ Creation de IPT6
  1029.       IPT6=0
  1030.       IF (NBELEM.EQ.0) GOTO 999
  1031.       NBNN=5
  1032.       NBREF=0
  1033.       NBSOUS=0
  1034.       SEGINI IPT6
  1035.       IPT6.ITYPEL=48
  1036. C
  1037. C 3/ Renseignement des noeuds support des relations
  1038.       DO 1004 J=1,KMILIE(/2)
  1039.       DO 28 I=1,KMILIE(/1)
  1040.         IF (KARET2(I,J).EQ.0) GOTO 28
  1041.         IF (KMILIE(I,J).GT.0) THEN
  1042.           NBREF=NBREF+1
  1043.           IPT6.NUM(1,NBREF)=KMILIE(I,J)
  1044.         ENDIF
  1045.    28  CONTINUE
  1046.  1004 CONTINUE
  1047. C       WRITE (*,*) 'nombre de noeuds supports de rela seg', NBREF
  1048. C
  1049. C 4/ Recherche des noeuds formant les relations
  1050.       DO 24 IARAF=1,MELVA2.VELCHE(/2)
  1051.         IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 24
  1052.         JCOMPT=0
  1053.         JPOS5=LPOS5(IPT2.ITYPEL)
  1054.         DO 23 I=1,NBRAF(IPT2.ITYPEL)
  1055.           JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
  1056.           JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
  1057.           JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
  1058.           LTYPNO=JTYPNO(JPOS5-1+I)
  1059.           IF (LTYPNO.NE.0) GOTO 22
  1060.           NPTA=IPT2.NUM(LISNOE(JLISN),IARAF)
  1061.           NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF)
  1062.           NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  1063.           NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  1064.           DO 29 K=1,MAX(1,KARPOS(NMIN))
  1065.             IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K
  1066.   29      CONTINUE
  1067.           IF (KARET2(NMIN,NEXIST).EQ.0) GOTO 22
  1068.           IF (KMILIE(NMIN,NEXIST).EQ.0) GOTO 22
  1069.           NEXIS5=0
  1070.           DO 20 MM=1,IPT6.NUM(/2)
  1071.             INOEGL=KMILIE(NMIN,NEXIST)
  1072.             INRELA=IPT6.NUM(1,MM)
  1073.             IF (INOEGL.EQ.INRELA) NEXIS5=MM
  1074.   20      CONTINUE
  1075. C
  1076. C 5/ Renseignement des noeuds formant les relations
  1077.           DO 21 IHH=1,JLONG
  1078.             IPT6.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF)
  1079.   21      CONTINUE
  1080.           IF (JLONG.EQ.2) IPT6.NUM(5,NEXIS5)=1
  1081.           IF (JLONG.EQ.3) IPT6.NUM(5,NEXIS5)=2
  1082.   22      CONTINUE
  1083.         JCOMPT=JCOMPT+JLONG
  1084.   23    CONTINUE
  1085.   24  CONTINUE
  1086.  
  1087.  444  CONTINUE
  1088.  
  1089. C
  1090. C============================================
  1091. C C) Creation du maillage de relations final !
  1092. C============================================
  1093.       IF (IPT5.EQ.0) THEN
  1094.         IPT7=IPT6
  1095.         GOTO 999
  1096.       ENDIF
  1097.       NBELEM=IPT5.NUM(/2)+IPT6.NUM(/2)
  1098. C      NBNN=MAX(IPT5.NUM(/1),IPT6.NUM(/1))
  1099.       NBNN=10
  1100.       NBREF=0
  1101.       NBSOUS=0
  1102.       SEGINI IPT7
  1103.       IPT7.ITYPEL=48
  1104.       DO 1005 NEO=1,IPT5.NUM(/2)
  1105.       DO 42 MOR=1,10
  1106.         IPT7.NUM(MOR,NEO)=IPT5.NUM(MOR,NEO)
  1107.         IPT7.ICOLOR(NEO)=IPT5.ICOLOR(NEO)
  1108.   42  CONTINUE
  1109.  1005 CONTINUE
  1110.       NN5 = IPT5.NUM(/2)
  1111.       DO 43 NEO=IPT5.NUM(/2)+1,IPT6.NUM(/2)+IPT5.NUM(/2)
  1112.       DO 43 MOR=1,IPT6.NUM(/1)
  1113.         IF (MOR.LT.IPT6.NUM(/1)) THEN
  1114.           IPT7.NUM(MOR,NEO)=IPT6.NUM(MOR,NEO-NN5)
  1115.         ELSE
  1116.           IPT7.NUM(10,NEO)=IPT6.NUM(MOR,NEO-NN5)
  1117.         ENDIF
  1118.         IPT7.ICOLOR(NEO)=IPT6.ICOLOR(NEO-NN5)
  1119.   43  CONTINUE
  1120.       SEGDES IPT5, IPT6
  1121.  
  1122. C=======================================================================
  1123. C                      Fin du programme
  1124. C=======================================================================
  1125.  999  CONTINUE
  1126. c      write (*,*) 'IPT2', IPT2
  1127. c      write (*,*) 'IPT3', IPT3
  1128. c      write (*,*) 'IPT4', IPT4
  1129. c      write (*,*) 'IPT5', IPT5
  1130. c      write (*,*) 'IPT6', IPT6
  1131. c      write (*,*) 'IPT7', IPT7
  1132.       RETURN
  1133.       END
  1134.  
  1135.  
  1136.  
  1137.  
  1138.  
  1139.  
  1140.  
  1141.  
  1142.  
  1143.  
  1144.  
  1145.  
  1146.  
  1147.  
  1148.  
  1149.  
  1150.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales