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  24/04/12    21:17:01     11897          
  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 5 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.         LCOMP=LCOMP+1
  486.    6  CONTINUE
  487. C
  488.       NBPT =NBPT1
  489. C      SEGADJ MCOORD,XDEN
  490. C
  491. C=======================================================================
  492. C              Preparation du maillage de relations
  493. C=======================================================================
  494.       SEGSUP KARET2
  495.       NBNDS=KARETE(/1)
  496.       NCOL=KARETE(/2)
  497.       SEGINI KARET2
  498. C==================================
  499. C A) Relations dues aux faces (3D) !
  500. C==================================
  501. C Tous les noeuds qui restent dans JFARAF sont a creer en tant que
  502. C relations de conformite ou des bords du domaine
  503.  
  504. C-----------------------------------------------------------------------
  505. C ** A.1 / Prise en compte des relation 
  506. C          d'arretes pour les arretes de chaque faces 
  507. C-----------------------------------------------------------------------
  508.  
  509. C-----------------------------------------------------------------------
  510. C** A.1.1/ faces raffinée
  511. C-----------------------------------------------------------------------
  512.  
  513. c      write (*,*) 'relations dues au faces'
  514.       JTYPEL=IPT2.ITYPEL
  515.       JLTEL1=LTEL(1,JTYPEL)
  516.       JLTEL2=LTEL(2,JTYPEL)
  517.       IF (JTYPEL.NE.14) GOTO 141
  518. c Boucle sur les elements de l'ancien maillage
  519.       DO 121 KELEM1=1,IPT2.NUM(/2)
  520. c Boucle sur les faces de l'ancien maillage
  521.         DO 122 IBB1=1,JLTEL1
  522.           JLDEL1=LDEL(1,JLTEL2-1+IBB1)
  523.           JTYFAC=IWORK1(JLDEL1)
  524.           JLDEL2=LDEL(2,JLTEL2-1+IBB1)
  525.           JNBSOM=NBSOM(JTYFAC)
  526.           JSPOS=NSPOS(JTYFAC)
  527.           SEGINI IWORK2
  528. C
  529. C 1/ trouver le numéro de la face
  530.           DO 123 IAA=1,JNBSOM
  531.               NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM1)
  532.               IWORK2(IAA)=NGLOBA
  533.  123      CONTINUE
  534. C
  535. C 1.1/ Classement des 3 sommets par ordre croissant de n° globale
  536.           NPTA=NBPT +1
  537.           NPTB=NPTA+1
  538.           NPTC=NPTB+1
  539.           DO 124 ICC=1,JNBSOM
  540.             IF (IWORK2(ICC).LT.NPTA) THEN
  541.               NPTC=NPTB
  542.               NPTB=NPTA
  543.               NPTA=IWORK2(ICC)
  544.             ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  545.               NPTC=NPTB
  546.               NPTB=IWORK2(ICC)
  547.             ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  548.               NPTC=IWORK2(ICC)
  549.             ENDIF
  550.  124      CONTINUE
  551. C
  552. C1.2/ Classement des 3 sommets par ordre croissant de n° locale
  553.           NPTA2=ICPR(NPTA,1)
  554.           NPTB2=ICPR(NPTB,1)
  555.           NPTC2=ICPR(NPTC,1)
  556.           IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  557.             NPTA=NPTA2
  558.             NPTB=MIN(NPTB2,NPTC2)
  559.             NPTC=MAX(NPTB2,NPTC2)
  560.           ENDIF
  561.           IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  562.             NPTA=NPTB2
  563.             NPTB=MIN(NPTA2,NPTC2)
  564.             NPTC=MAX(NPTA2,NPTC2)
  565.           ENDIF
  566.           IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  567.             NPTA=NPTC2
  568.             NPTB=MIN(NPTA2,NPTB2)
  569.             NPTC=MAX(NPTA2,NPTB2)
  570.           ENDIF
  571. C
  572. C 1.3/ Recherche du numero de la face
  573.           NEXIS2=0
  574.           DO 125 IEE=1,JPLCOM(NPTA)
  575.             MTMP=JPLANS(NPTA,IEE)
  576.             JJ1=JNOEFA(MTMP,1)
  577.             JJ2=JNOEFA(MTMP,2)
  578.             JJ3=JNOEFA(MTMP,3)
  579.             IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  580.               NEXIS2=IEE
  581.             ENDIF
  582.  125      CONTINUE
  583.  
  584.           JNUMF1=JPLANS(NPTA,NEXIS2)
  585.  
  586. C
  587. C 2/ Exclusion des faces non voulues 
  588.  
  589. C on exlu les faces où un noeud existait déja
  590.           IF (JFARAF(JNUMF1,1).EQ.0) GOTO 122
  591. C on exlu les faces de bords
  592.           IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 122
  593.  
  594. C 3/ Recherche des arretes de la faces pour y mettre des relations 
  595.  
  596.           IF (JTYPEL.NE.14) GOTO 137
  597.           NBAR1=0
  598.           DO 138 ISS=1, (JNBSOM-1)
  599.           DO 138 JSS=(ISS+1), JNBSOM
  600.  
  601.             NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  602.             NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  603.  
  604.             DO 139 IRR=1,MAX(1,KARPOS(NMIN))                   
  605.               IF (KARETE(NMIN,IRR).EQ.NMAX) THEN
  606.                 NBAR1=NBAR1+1
  607.                 KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1
  608. c                write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  609. c                IF (JFARAF(JNUMF1,1).EQ.15748) THEN
  610. c                  write (*,*) ISS,JSS 
  611. c                  write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!'
  612. c                  write (*,*) IWORK2(ISS), IWORK2(JSS)
  613. c                  write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  614. c                  write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR) 
  615. c                ENDIF
  616.               ENDIF
  617.  139        CONTINUE
  618.  138      CONTINUE 
  619. c          WRITE (*,*) 'Nbarrete', NBAR1
  620. c          IF (NBAR1.NE.4) THEN
  621. c             write (*,*) IWORK2(1), IWORK2(2)
  622. c             write (*,*) IWORK2(3), IWORK2(4) 
  623. c          ENDIF
  624.  137      CONTINUE    
  625.  122    CONTINUE        
  626.  121  CONTINUE
  627.  
  628.  141  CONTINUE
  629.  
  630. C-----------------------------------------------------------------------
  631. C** A.1.2/ relations de faces conservées dans le maillage précédent
  632. C-----------------------------------------------------------------------
  633.  
  634.       SEGINI IPT5
  635.       DO 670 MIE=1,MAX(1,IPT3.LISOUS(/1))
  636.         IF (IPT3.LISOUS(/1).EQ.0) THEN
  637.           IPT5=IPT3
  638.         ELSE
  639.           IPT5=IPT3.LISOUS(MIE)
  640.         ENDIF
  641.         SEGACT IPT5
  642. c        write (*,*) 'ipt3 raff 3D' , ipt3, IPT3.ITYPEL, IPT3.ICOLOR(1)
  643. c        write (*,*) 'ipt5 raff 3D' , ipt5, IPT5.ITYPEL, IPT5.ICOLOR(1)
  644.         IF (IPT5.ITYPEL.NE.48) GOTO 670
  645.         IF ((IPT5.ICOLOR(1).NE.3)) GOTO 670
  646.         LIGNUM=IPT2.NUM(/1)
  647.         JNBSOM=4
  648.         SEGINI IWORK2
  649.  
  650. C
  651. C 1/ trouver le numéro de la face
  652. C 1.1/ Classement des 3 sommets par ordre croissant (num globale)
  653.         DO 671,MAA=1,IPT5.NUM(/2)
  654.           DO 673 MBB=1,JNBSOM
  655. c            write (*,*) 'MAA' , MAA, 'MBB', MBB
  656.             IWORK2(MBB)=IPT5.NUM(MBB+1,MAA)
  657.  
  658.  673      CONTINUE
  659.           NPTA=nbpt +1
  660.           NPTB=NPTA+1
  661.           NPTC=NPTB+1
  662.           DO 672 ICC=1,JNBSOM
  663.             IF (IWORK2(ICC).LT.NPTA) THEN
  664.               NPTC=NPTB
  665.               NPTB=NPTA
  666.               NPTA=IWORK2(ICC)
  667.             ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  668.               NPTC=NPTB
  669.               NPTB=IWORK2(ICC)
  670.             ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  671.               NPTC=IWORK2(ICC)
  672.             ENDIF
  673.  672      CONTINUE
  674. C
  675. C 1.2/ Classement des 3 sommets par ordre croissant (num locale)
  676.           NPTA1=NPTA
  677.           NPTB1=NPTB
  678.           NPTC1=NPTC
  679.           NPTA2=ICPR(NPTA,1)
  680.           NPTB2=ICPR(NPTB,1)
  681.           NPTC2=ICPR(NPTC,1)
  682.           IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  683.             NPTA=NPTA2
  684.             NPTB=MIN(NPTB2,NPTC2)
  685.             NPTC=MAX(NPTB2,NPTC2)
  686.           ENDIF
  687.           IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  688.             NPTA=NPTB2
  689.             NPTB=MIN(NPTA2,NPTC2)
  690.             NPTC=MAX(NPTA2,NPTC2)
  691.           ENDIF
  692.           IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  693.             NPTA=NPTC2
  694.             NPTB=MIN(NPTA2,NPTB2)
  695.             NPTC=MAX(NPTA2,NPTB2)
  696.           ENDIF
  697. C
  698. C 1.3/ Recherche du numero de la face
  699.           NEXIS2=0
  700.           DO 674 IEE=1,JPLCOM(NPTA)
  701.             MTMP=JPLANS(NPTA,IEE)
  702.             JJ1=JNOEFA(MTMP,1)
  703.             JJ2=JNOEFA(MTMP,2)
  704.             JJ3=JNOEFA(MTMP,3)
  705.             IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  706.               NEXIS2=IEE
  707.             ENDIF
  708.  674      CONTINUE
  709.           JNUMF1=JPLANS(NPTA,NEXIS2)
  710. c          IF (MAA.EQ.5) THEN
  711. c                  write (*,*) 'MAA', MAA, 'nexis2', nexis2 
  712. c                  write (*,*) IWORK2(1), IWORK2(2), IWORK2(3), IWORK2(4) 
  713. c                  write (*,*) 'JFARAF(JNUMF1,1)', JFARAF(JNUMF1,1)
  714. c                  write (*,*) 'JNOEFA(JNUMF1,5)',JNOEFA(JNUMF1,5) 
  715. c          ENDIF
  716. C
  717. C 2/ Exclusion des faces non voulues
  718.           IF (NEXIS2.EQ.0) THEN
  719. c            WRITE (*,*)'FACE NON TROUVE SURE:', MAA 
  720.             GOTO 671
  721.           ENDIF
  722. C on exlu les faces où un noeud existait déja
  723.           IF (JFARAF(JNUMF1,1).EQ.0) GOTO 671
  724.  
  725. C on exlu les faces de bords
  726.           IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 671
  727.           NBAR1=0
  728. C
  729. C 3/ Recherche des arretes de la faces pour y mettre des relations
  730.  
  731.           DO 135 ISS=1, (JNBSOM-1)
  732.           DO 135 JSS=(ISS+1), JNBSOM
  733.  
  734.             NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  735.             NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  736.  
  737.             DO 136 IRR=1,MAX(1,KARPOS(NMIN))                   
  738.               IF (KARETE(NMIN,IRR).EQ.NMAX) THEN
  739.                 NBAR1=NBAR1+1
  740.                 KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1
  741. c                write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  742. c                IF (MAA.EQ.5) THEN
  743. c                  write (*,*) ISS,JSS 
  744. c                  write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!'
  745. c                  write (*,*) IWORK2(ISS), IWORK2(JSS)
  746. c                  write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  747. c                  write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR) 
  748. c                ENDIF
  749.               ENDIF
  750.  136        CONTINUE
  751.  135      CONTINUE 
  752.  
  753.  671    CONTINUE
  754.  
  755.  670  CONTINUE
  756.       SEGDES IPT5
  757. C---------------------------------
  758. C ** A.2 / Initialisation de IPT5 !
  759. C---------------------------------
  760.  
  761.  
  762.  
  763. C 1/ Comptage du nombre de noeuds soumis a des relations
  764.       NBRELA=0
  765.       DO 114 IHF=1,JFARAF(/1)
  766. C on exlu les faces de bords
  767.         IF (JNOEFA(IHF,5).EQ.0) GOTO 114
  768. C on exlu les faces où un noeud existait déja
  769.         IF (JFARAF(IHF,1).EQ.0) GOTO 114
  770.         DO 113 JF=1,JFARAF(/2),2
  771.           IF (JFARAF(IHF,JF).NE.0) NBRELA=NBRELA+1
  772.  113    CONTINUE
  773.  114  CONTINUE
  774. C
  775. C 2/ Creation de IPT5
  776.       IF (NBRELA.EQ.0) THEN
  777.         IPT5=0
  778.         GOTO 210
  779.       ENDIF
  780. C      IF (IPT2.ITYPEL.EQ.14) NBNN=4
  781. C      IF (IPT2.ITYPEL.EQ.15) NBNN=8
  782. C      IF (IPT2.ITYPEL.EQ.17) NBNN=5
  783. C      IF (IPT2.ITYPEL.EQ.24) NBNN=5
  784.       NBNN=10
  785.       NBELEM=NBRELA
  786.       NBSOUS=0
  787.       NBREF=0
  788.       SEGINI IPT5
  789.       IPT5.ITYPEL=48
  790. C
  791. C 3/ Renseignement des noeuds support des relations
  792.       NBRELA=0
  793.       DO 116 IPF=1,JFARAF(/1)
  794.         IF (JNOEFA(IPF,5).EQ.0) GOTO 116
  795.         IF (JFARAF(IPF,1).EQ.0) GOTO 116
  796.         DO 115 JF=1,JFARAF(/2),2
  797.           IF (JFARAF(IPF,JF).EQ.0) GOTO 115
  798.           NBRELA=NBRELA+1
  799.           IPT5.NUM(1,NBRELA)=JFARAF(IPF,JF)
  800.  115    CONTINUE
  801.  116  CONTINUE
  802. c      write (*,*) 'NBRELA FACE ', NBRELA
  803. C
  804. C-----------------------------------------------------
  805. C ** A.3 / Recherche des noeuds formant les relations !
  806. C-----------------------------------------------------
  807. C 1/ Boucle sur l'ensemble des noeuds a creer
  808.       DO 200 IARAF=1,MELVA2.VELCHE(/2)
  809.         IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 200
  810.         JCOMPT=0
  811.         JPOS5=LPOS5(IPT2.ITYPEL)
  812.         DO 190 I=1,NBRAF(IPT2.ITYPEL)
  813.           JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
  814.           JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
  815.           JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
  816.           LTYPNO=JTYPNO(JPOS5-1+I)
  817.           IF ((LTYPNO.EQ.0).OR.(LTYPNO.EQ.7)) GOTO 189
  818. C
  819. C 2/ Preparation pour trouver le noeud et la face en question
  820.           JTYPEL=IPT2.ITYPEL
  821.           JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO
  822.           JLDEL1=LDEL(1,JLTEL2)
  823.           JTYFAC=IWORK1(JLDEL1)
  824.           JLDEL2=LDEL(2,JLTEL2)
  825.           JNBSOM=NBSOM(JTYFAC)
  826.           JSPOS=NSPOS(JTYFAC)
  827.           SEGINI IWORK2
  828. C
  829. C 3/ Classement des 3 sommets par ordre croissant de n° globale
  830.           DO 100 IAA=1,JNBSOM
  831.             NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF)
  832.             IWORK2(IAA)=NGLOBA
  833.  100      CONTINUE
  834.           NPTA=nbpt +1
  835.           NPTB=NPTA+1
  836.           NPTC=NPTB+1
  837.           DO 110 ICC=1,JNBSOM
  838.             IF (IWORK2(ICC).LT.NPTA) THEN
  839.               NPTC=NPTB
  840.               NPTB=NPTA
  841.               NPTA=IWORK2(ICC)
  842.             ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  843.               NPTC=NPTB
  844.               NPTB=IWORK2(ICC)
  845.             ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  846.               NPTC=IWORK2(ICC)
  847.             ENDIF
  848.  110      CONTINUE
  849. C
  850. C 4/ Classement des 3 sommets par ordre croissant de n° locale
  851.           NPTA2=ICPR(NPTA,1)
  852.           NPTB2=ICPR(NPTB,1)
  853.           NPTC2=ICPR(NPTC,1)
  854.           IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  855.             NPTA=NPTA2
  856.             NPTB=MIN(NPTB2,NPTC2)
  857.             NPTC=MAX(NPTB2,NPTC2)
  858.           ENDIF
  859.           IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  860.             NPTA=NPTB2
  861.             NPTB=MIN(NPTA2,NPTC2)
  862.             NPTC=MAX(NPTA2,NPTC2)
  863.           ENDIF
  864.           IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  865.             NPTA=NPTC2
  866.             NPTB=MIN(NPTA2,NPTB2)
  867.             NPTC=MAX(NPTA2,NPTB2)
  868.           ENDIF
  869. C
  870. C 5/ Recherche du numero de la face
  871.           NEXIS2=0
  872.           DO 120 IEE=1,JPLCOM(NPTA)
  873.             MTMP=JPLANS(NPTA,IEE)
  874.             JJ1=JNOEFA(MTMP,1)
  875.             JJ2=JNOEFA(MTMP,2)
  876.             JJ3=JNOEFA(MTMP,3)
  877.             IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  878.               NEXIS2=IEE
  879.             ENDIF
  880.  120      CONTINUE
  881.           JNUMFA=JPLANS(NPTA,NEXIS2)
  882.  
  883. C
  884. C 6/ Recherche du numero global du point
  885.           IF (JNOEFA(JNUMFA,5).EQ.0) GOTO 189
  886.           IF (JTYFAC.EQ.8) INOEGL=JFARAF(JNUMFA,1)
  887.           IF (JTYFAC.NE.8) THEN
  888.             NEXIS3=0
  889.             NEXIS4=0
  890.             XCO2=0.25
  891.             DO 140 KBB=1,JLONG
  892.               XCO1=XCOEFF(JPOS1-1+KBB)
  893.               IF (XCO1.EQ.XCO2) NEXIS3=KBB
  894.  140        CONTINUE
  895.             IF (NEXIS3.NE.0) THEN
  896.               NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
  897.               DO 150 KAA=2,JFARAF(/2),2
  898.                 IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA
  899.  150          CONTINUE
  900.               INOEGL=JFARAF(JNUMFA,NEXIS4-1)
  901.             ELSE
  902.               DO 160 KAA=2,JFARAF(/2),2
  903.                 IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA
  904.  160          CONTINUE
  905.               INOEGL=JFARAF(JNUMFA,NEXIS4-1)
  906.             ENDIF
  907.           ENDIF
  908.  
  909. C
  910. C------------------------------
  911. C ** A.4 / Remplissage de IPT5 !
  912. C------------------------------
  913. C 1/ Recherche de la position du point dans IPT5
  914.           NEXIS5=0
  915.           DO 170 IGG=1,IPT5.NUM(/2)
  916.             IF (INOEGL.EQ.IPT5.NUM(1,IGG)) NEXIS5=IGG
  917.  170      CONTINUE
  918.           IF (NEXIS5.EQ.0) GOTO 189
  919.           IF (INOEGL.EQ.139) THEN
  920. c            write (*,*) INOEGL, IWORK2(1), IWORK2(2), IWORK2(3)  
  921.           ENDIF
  922.  
  923. C
  924. C 2/ Renseignement des points formant les relations
  925.           DO 180 IHH=1,JLONG
  926.             IPT5.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF)
  927.  180      CONTINUE
  928.           IF (JLONG.EQ.4) IPT5.NUM(10,NEXIS5)=3
  929.           IF (JLONG.EQ.5) IPT5.NUM(10,NEXIS5)=4
  930.           IF (JLONG.EQ.8) THEN
  931.             IF (JPOS1.EQ.16) IPT5.NUM(10,NEXIS5)=5
  932.             IF (JPOS1.EQ.24) IPT5.NUM(10,NEXIS5)=6
  933.           ENDIF
  934.  189      CONTINUE
  935.         JCOMPT=JCOMPT+JLONG
  936.  190    CONTINUE
  937.  200  CONTINUE
  938.  210  CONTINUE
  939. C
  940. C===================================
  941. C B) Relations dues aux aretes (2D) !
  942. C===================================
  943. C On cree un maillage IPT6 contenant tous les noeuds soumis a des
  944. C relations
  945. C
  946.  
  947. c        ILPL=LPL(IPT4.ITYPEL)
  948. c        ILPT=LPT(IPT4.ITYPEL)
  949. c        DO 317 J=1,IPT4.NUM(/2)
  950. c        DO 317 K=1,ILPL*2-1,2
  951. c          NPTA=IPT4.NUM(KSEGM(ILPT+K-1),J)
  952. c          NPTB=IPT4.NUM(KSEGM(ILPT+K),J)
  953. c
  954. c          if ((NPTB.eq.129).and.(NPTA.eq.9087)) then 
  955. c            write(*,*) '!ici!',NBPT0 
  956. c          endif
  957. c          IF((NPTA.GT.NBPT0).OR.(NPTB.GT.NBPT0)) THEN
  958. c             GOTO 317
  959. c          ENDIF
  960. c          NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  961. c          NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  962. c          if ((NMIN.eq.235).and.(NMAX.eq.236)) then
  963. c              write(*,*) 'nmin', nmin, 'nmax', nmax
  964. c              write(*,*) 'NPTA',NPTA, 'NPTB', NPTB
  965. c          endif
  966. c          NEXIST=0
  967. c          DO 316 I=1,MAX(1,KARPOS(NMIN))
  968. c            IF (KARETE(NMIN,I).EQ.NMAX) THEN
  969. c              KARET2(NMIN,I)=KARET2(NMIN,I)+1
  970. c
  971. c              IF (KMILIE(NMIN,I).gt.0) then
  972. c prise en compte des arretes multiples sur le nouveau maillage
  973. c dans l'exemple ci dessous où les noeuds C et D sont des hanging nodes
  974. c L'arette [C B] n'a pas encore été prise en compte dans karet2 c'est 
  975. c ce qui est fait ici.
  976. c
  977. c |                   |
  978. c |                   |
  979. c |A        C    D    |B
  980. c x---------X----X----X
  981. c |         |    |    |
  982. c |         |    |    |
  983.  
  984. c1/ Test sur la premiere arete : NMIN-KMILIE(NMIN,I)
  985. c                NMILI = ICPR(KMILIE(NMIN,I),1)
  986. c                NEXIS2=0
  987. c                NMIN2=MIN(NMIN,NMILI)
  988. c                NMAX2=MAX(NMIN,NMILI)
  989. c                DO 318 I2=1,KMILIE(/2)
  990. c                  IF (KMILIE(NMIN2,I2).GT.0) THEN
  991. c                    KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1
  992. c                  ENDIF
  993. c 318            CONTINUE
  994. C
  995. C 2/ Test sur la deuxieme arete : NMAX-NPTC
  996. c                NEXIS2=0
  997. c                NMIN2=MIN(NMAX,NMILI)
  998. c                NMAX2=MAX(NMAX,NMILI)
  999. c                DO 319 I2=1,KMILIE(/2)
  1000. c                  IF (KMILIE(NMIN2,I2).GT.0) THEN
  1001. c                    KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1
  1002. c                  ENDIF
  1003. c 319            CONTINUE
  1004. c              ENDIF
  1005. c            ENDIF
  1006. c 316      CONTINUE
  1007. c 317    CONTINUE
  1008.  
  1009. C 1/ Comptage du nombre de noeuds soumis a des relations
  1010.       NBELEM=0
  1011.       DO 27 J=1,KMILIE(/2)
  1012.       DO 27 I=1,KMILIE(/1)
  1013.  
  1014. c        IF (KMILIE(I,J).eq.10330) write(*,*) '!kmilie', I, J, '=10330!'
  1015. c        IF (KMILIE(I,J).eq.10330) write(*,*) 'KARET2(I,J)', KARET2(I,J)
  1016. c        IF (KMILIE(I,J).eq.10341) write(*,*) '!kmilie', I, J, '=10341!'
  1017. c        IF (KMILIE(I,J).eq.10341) write(*,*) 'KARET2(I,J)', KARET2(I,J)
  1018.         IF (KARET2(I,J).EQ.0) GOTO 27
  1019.         IF (KMILIE(I,J).GT.0) NBELEM=NBELEM+1
  1020.  
  1021.   27  CONTINUE
  1022. c      write (*,*) 'NBELEM', NBELEM
  1023. C
  1024. C 2/ Creation de IPT6
  1025.       IPT6=0
  1026.       IF (NBELEM.EQ.0) GOTO 999
  1027.       NBNN=5
  1028.       NBREF=0
  1029.       NBSOUS=0
  1030.       SEGINI IPT6
  1031.       IPT6.ITYPEL=48
  1032. C
  1033. C 3/ Renseignement des noeuds support des relations
  1034.       DO 28 J=1,KMILIE(/2)
  1035.       DO 28 I=1,KMILIE(/1)
  1036.         IF (KARET2(I,J).EQ.0) GOTO 28
  1037.         IF (KMILIE(I,J).GT.0) THEN
  1038.           NBREF=NBREF+1
  1039.           IPT6.NUM(1,NBREF)=KMILIE(I,J)
  1040.         ENDIF
  1041.    28  CONTINUE
  1042. C       WRITE (*,*) 'nombre de noeuds supports de rela seg', NBREF
  1043. C
  1044. C 4/ Recherche des noeuds formant les relations
  1045.       DO 24 IARAF=1,MELVA2.VELCHE(/2)
  1046.         IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 24
  1047.         JCOMPT=0
  1048.         JPOS5=LPOS5(IPT2.ITYPEL)
  1049.         DO 23 I=1,NBRAF(IPT2.ITYPEL)
  1050.           JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
  1051.           JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
  1052.           JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
  1053.           LTYPNO=JTYPNO(JPOS5-1+I)
  1054.           IF (LTYPNO.NE.0) GOTO 22
  1055.           NPTA=IPT2.NUM(LISNOE(JLISN),IARAF)
  1056.           NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF)
  1057.           NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  1058.           NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  1059.           DO 29 K=1,MAX(1,KARPOS(NMIN))
  1060.             IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K
  1061.   29      CONTINUE
  1062.           IF (KARET2(NMIN,NEXIST).EQ.0) GOTO 22
  1063.           IF (KMILIE(NMIN,NEXIST).EQ.0) GOTO 22
  1064.           NEXIS5=0
  1065.           DO 20 MM=1,IPT6.NUM(/2)
  1066.             INOEGL=KMILIE(NMIN,NEXIST)
  1067.             INRELA=IPT6.NUM(1,MM)
  1068.             IF (INOEGL.EQ.INRELA) NEXIS5=MM
  1069.   20      CONTINUE
  1070. C
  1071. C 5/ Renseignement des noeuds formant les relations
  1072.           DO 21 IHH=1,JLONG
  1073.             IPT6.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF)
  1074.   21      CONTINUE
  1075.           IF (JLONG.EQ.2) IPT6.NUM(5,NEXIS5)=1
  1076.           IF (JLONG.EQ.3) IPT6.NUM(5,NEXIS5)=2
  1077.   22      CONTINUE
  1078.         JCOMPT=JCOMPT+JLONG
  1079.   23    CONTINUE
  1080.   24  CONTINUE
  1081.  
  1082.  444  CONTINUE
  1083.  
  1084. C
  1085. C============================================
  1086. C C) Creation du maillage de relations final !
  1087. C============================================
  1088.       IF (IPT5.EQ.0) THEN
  1089.         IPT7=IPT6
  1090.         GOTO 999
  1091.       ENDIF
  1092.       NBELEM=IPT5.NUM(/2)+IPT6.NUM(/2)
  1093. C      NBNN=MAX(IPT5.NUM(/1),IPT6.NUM(/1))
  1094.       NBNN=10
  1095.       NBREF=0
  1096.       NBSOUS=0
  1097.       SEGINI IPT7
  1098.       IPT7.ITYPEL=48
  1099.       DO 42 NEO=1,IPT5.NUM(/2)
  1100.       DO 42 MOR=1,10
  1101.         IPT7.NUM(MOR,NEO)=IPT5.NUM(MOR,NEO)
  1102.         IPT7.ICOLOR(NEO)=IPT5.ICOLOR(NEO)
  1103.   42  CONTINUE
  1104.       NN5 = IPT5.NUM(/2)
  1105.       DO 43 NEO=IPT5.NUM(/2)+1,IPT6.NUM(/2)+IPT5.NUM(/2)
  1106.       DO 43 MOR=1,IPT6.NUM(/1)
  1107.         IF (MOR.LT.IPT6.NUM(/1)) THEN
  1108.           IPT7.NUM(MOR,NEO)=IPT6.NUM(MOR,NEO-NN5)
  1109.         ELSE
  1110.           IPT7.NUM(10,NEO)=IPT6.NUM(MOR,NEO-NN5)
  1111.         ENDIF
  1112.         IPT7.ICOLOR(NEO)=IPT6.ICOLOR(NEO-NN5)
  1113.   43  CONTINUE
  1114.       SEGDES IPT5, IPT6
  1115.  
  1116. C=======================================================================
  1117. C                      Fin du programme
  1118. C=======================================================================
  1119.  999  CONTINUE
  1120. c      write (*,*) 'IPT2', IPT2
  1121. c      write (*,*) 'IPT3', IPT3
  1122. c      write (*,*) 'IPT4', IPT4
  1123. c      write (*,*) 'IPT5', IPT5
  1124. c      write (*,*) 'IPT6', IPT6
  1125. c      write (*,*) 'IPT7', IPT7
  1126.       RETURN
  1127.       END
  1128.  
  1129.  
  1130.  
  1131.  
  1132.  
  1133.  
  1134.  
  1135.  
  1136.  
  1137.  
  1138.  
  1139.  
  1140.  
  1141.  
  1142.  
  1143.  

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