Télécharger la plaquette Cast3M


fr en

Télécharger dedu.eso

Retour à la liste

Numérotation des lignes : 

dedu

  1. C DEDU      SOURCE    CB215821  24/04/12    21:15:34     11897          
  2.  
  3. C-----------------------------------------------------------------------
  4. C 1ere possibilite
  5. C            OBJ1  =  DEDU OBJ2 MAIT_ANC MAIT_NOU (REGU)
  6. C Deduction d'un maillage a partir d'un autre et du deplacement de
  7. C noeuds maitres
  8. C  AUTEUR : PV
  9. C-----------------------------------------------------------------------
  10. C 2eme possibilite
  11. C            OBJ2 = OBJ1 DEDU 'TRAN' GEO1 GEO2
  12. C 3eme possibilite
  13. C            OBJ2 = OBJ1 DEDU FLOT1 POIN1 (POIN2) 'ROTA' GEO1 GEO2
  14. C OBJ1 : type POINT, MAILLAGE, CHPOINT, MCHAML, MMODEL
  15. C Passage du support de OBJ1 a celui de OBJ2 comme celui de GEO1 a GEO2
  16. C Les points du support de OBJ1 doivent appartenir a GEO1
  17. C On realise eventuellement la rotation des composantes (envoi a PROPER)
  18. C  AUTEUR : KICH
  19. C-----------------------------------------------------------------------
  20. C  10/2003 : Prise en compte du cas IDIM=1 et des modes associes.
  21. C-----------------------------------------------------------------------
  22. C  Remarques :
  23. C  -----------
  24. C    DEDU est appele par l'operateur DEPL (deplac.eso) avec ICAS=1.
  25. C    Ce cas correspond a une transformation affine de l'espace (noeuds
  26. C    maitres = un maillage a (IDIM+1) points).
  27. C-----------------------------------------------------------------------
  28.  
  29.       SUBROUTINE DEDU(ICAS)
  30.  
  31.       IMPLICIT INTEGER(I-N)
  32.       IMPLICIT REAL*8 (A-H,O-Z)
  33.  
  34.  
  35. -INC PPARAM
  36. -INC CCOPTIO
  37. -INC CCREEL
  38. -INC SMCOORD
  39. -INC SMELEME
  40. -INC CCGEOME
  41. -INC SMRIGID
  42. -INC SMCHPOI
  43. -INC SMMODEL
  44. -INC SMCHAML
  45. -INC SMTEXTE
  46.       DIMENSION iuniq(10)
  47.       CHARACTER*4 mcle(5)
  48.       CHARACTER*(LOCOMP) strcomp
  49.  
  50.       SEGMENT kon(nkon,numnp,2)
  51.       SEGMENT xkon(nkon+1,numnp)
  52.       SEGMENT ipxkon(nkon+1,numnp)
  53.       SEGMENT ikon(numnp)
  54.       SEGMENT icpr(nbpts)
  55.       SEGMENT icp2(nbpts)
  56.       SEGMENT idcp(numnp)
  57.       SEGMENT icorr(numnp)
  58.       SEGMENT xxtra
  59.         REAL*8 xdep(3,4),xarr(3,4),ved1(3),ved2(3),ved3(3),vea1(3),
  60.      .         vea2(3),vea3(3),pa(3,3),pa2(2,2),pa3(3,3),vb(3)
  61.       END SEGMENT
  62.  
  63.       DATA mcle / 'REGU','ELIM','TRAN','ROTA','ADAP' /
  64.  
  65.       segact mcoord*mod
  66.       idimp1=IDIM+1
  67.  
  68. C  Mot cle optionnel REGU : coordonnees barycentriques identiques
  69.       CALL LIRMOT(mcle,5,iregu,0)
  70.  
  71. C  2eme possibilite : option TRAN
  72. C --------------------
  73.       IF (iregu.EQ.3) THEN
  74.         CALL PROPER(4)
  75.         RETURN
  76.       ENDIF
  77. C  3eme possibilite : option ROTA
  78. C --------------------
  79. C  Option non disponible en dimension 1 (pas de signification)
  80.       IF (iregu.EQ.4) THEN
  81.         IF (IDIM.EQ.1) THEN
  82.           MOTERR(1:4)=mcle(iregu)
  83.           INTERR(1)=IDIM
  84.           CALL ERREUR(971)
  85.         ELSE
  86.           CALL PROPER(5)
  87.         ENDIF
  88.         RETURN
  89.       ENDIF
  90. C  5eme possibilite : option ADAP
  91. C  On passe la main à la procédure DEDUADAP
  92. C --------------------
  93.       IF (iregu.EQ.5) THEN
  94.          SEGINI MTEXTE
  95.          LTT=8
  96.          MTEXT(1:LTT) ='DEDUADAP'
  97.          NCART=8
  98.          SEGDES MTEXTE
  99.          CALL ECROBJ('TEXTE',MTEXTE)
  100.          RETURN
  101.       ENDIF
  102.  
  103. C  Par defaut : option ELIM (ignore les parties concaves)
  104. C  Correspond au cas ICAS=1 (operateur DEPL 'DEDU')
  105.       IF (iregu.EQ.0) iregu=2
  106.  
  107. C  1ere possibilite : options REGU et ELIM / appel par DEPLAC
  108. C --------------------
  109. C  Lecture du maillage initial
  110.       CALL LIROBJ('MAILLAGE',ipt1,1,iretou)
  111.       IF (ierr.NE.0) RETURN
  112. C  Lecture des points maitres initiaux
  113.       CALL LIROBJ('MAILLAGE',ipt2,0,iretou)
  114.       IF (iretou.NE.0) THEN
  115.         ichpo=0
  116. C  lecture ddes points maitres finaux
  117.         CALL LIROBJ('MAILLAGE',ipt3,1,iretou)
  118.          IF (ierr.NE.0) RETURN
  119.       ELSE
  120. C  Lecture d'un chpoint de deplacement
  121.         ichpo=1
  122.         CALL LIROBJ('CHPOINT ',mchpoi,1,iretou)
  123.         IF (ierr.NE.0) RETURN
  124. C  Création de ipt2 et ipt3 pour se ramener au cas précédent
  125.         nbpts0=nbpts
  126. C  Comptage des points du chpoint concernes par UX UY UZ UR
  127.         SEGACT mchpoi
  128.         jatan=jattri(1)
  129.         nbpcon=0
  130.         DO nso=1,ipchp(/1)
  131.           msoupo=ipchp(nso)
  132.           SEGACT msoupo
  133.           ipt6=igeoc
  134.           SEGACT ipt6
  135.           ipc=0
  136.           DO i=1,nocomp(/2)
  137.             strcomp=nocomp(i)
  138.             IF ((strcomp.EQ.'UX  ').OR.(strcomp.EQ.'UR  ').OR.
  139.      .          ((strcomp.EQ.'UY  ').AND.(IDIM.NE.1)).OR.
  140.      .          ((strcomp.EQ.'UZ  ').AND.(IDIM.NE.1))) ipc=ipc+i
  141.           ENDDO
  142.           IF (ipc.NE.0) nbpcon=nbpcon+ipt6.num(/2)
  143.         ENDDO
  144. C  Test sur nbpcon=0 a ecrire
  145.         nbref=0
  146.         nbsous=0
  147.         nbelem=nbpcon
  148.         nbnn=1
  149.         SEGINI ipt2
  150.         SEGINI ipt3
  151.         ipt2.itypel=1
  152.         ipt3.itypel=1
  153. C  Remplissage de ipt2 et ipt3 - Creation des points de ipt3
  154.         idej=0
  155.         DO 9 nso=1,ipchp(/1)
  156.           msoupo=ipchp(nso)
  157.           ipt6=igeoc
  158.           ipox=0
  159.           ipoy=0
  160.           ipoz=0
  161.           DO i=1,nocomp(/2)
  162.             strcomp=nocomp(i)
  163.             IF ((strcomp.EQ.'UX  ').OR.(strcomp.EQ.'UR  ')) ipox=i
  164.             IF ((strcomp.EQ.'UY  ').AND.(IDIM.NE.1)) ipoy=i
  165.             IF (strcomp.EQ.'UZ  ') THEN
  166.               IF (IDIM.EQ.2) THEN
  167.                 ipoy=i
  168.               ELSE IF (IDIM.GE.3) THEN
  169.                 ipoz=i
  170.               ENDIF
  171.             ENDIF
  172.           ENDDO
  173.           IF ((ipox+ipoy+ipoz).EQ.0) GOTO 9
  174.           mpoval=ipoval
  175.           SEGACT mpoval
  176.           nbpan=nbpts
  177.           nbpts=nbpts+ipt6.num(/2)
  178.           SEGADJ mcoord
  179.           DO nel=1,ipt6.num(/2)
  180.             iep=ipt6.num(1,nel)
  181.             ipt2.num(1,idej+nel)=iep
  182.             ipt3.num(1,idej+nel)=nbpan+nel
  183.             irefep=(iep-1)*idimp1
  184.             iref=(nbpan+nel-1)*idimp1
  185.             IF (ipox.NE.0) THEN
  186.               xcoor(iref+1)=xcoor(irefep+1)+vpocha(nel,ipox)
  187.             ELSE
  188.               xcoor(iref+1)=xcoor(irefep+1)
  189.             ENDIF
  190.             IF (ipoy.NE.0) THEN
  191.               xcoor(iref+2)=xcoor(irefep+2)+vpocha(nel,ipoy)
  192.             ELSE
  193.               IF (IDIM.GE.2) xcoor(iref+2)=xcoor(irefep+2)
  194.             ENDIF
  195.             IF (ipoz.NE.0) THEN
  196.               xcoor(iref+3)=xcoor(irefep+3)+vpocha(nel,ipoz)
  197.             ELSE
  198.               IF (IDIM.GE.3) xcoor(iref+3)=xcoor(irefep+3)
  199.             ENDIF
  200.             xcoor(iref+idimp1)=xcoor(irefep+idimp1)
  201.           ENDDO
  202.           idej=idej+ipt6.num(/2)
  203.           SEGDES ipt6,mpoval,msoupo
  204.  9      CONTINUE
  205.         SEGDES mchpoi
  206.       ENDIF
  207.  
  208. C  Verification compatibilite points maitres initiaux et finaux
  209. C  Cas ichpo = 0 : cette verification est a faire lorsque la
  210. C  transformation est definie par deux maillages.
  211. C  Cas ichpo = 1 : les points maitres initiaux et finaux sont crees
  212. C  en parallele de la meme maniere a partir du CHPO donne. La veri-
  213. C  fication est donc inutile.
  214.       IF (ichpo.EQ.1) GOTO 100
  215.       SEGACT ipt2,ipt3
  216.       IF (ipt2.lisous(/1).NE.ipt3.lisous(/1)) THEN
  217.         CALL ERREUR(16)
  218.         SEGDES ipt2,ipt3
  219.         RETURN
  220.       ENDIF
  221.       ipt7=ipt2
  222.       ipt8=ipt3
  223.       DO i=1,MAX(1,ipt2.lisous(/1))
  224.         IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(i)
  225.         IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(i)
  226.         SEGACT ipt7,ipt8
  227.         moterr=' '
  228.         IF (ipt7.num(/1).NE.ipt8.num(/1)) moterr(1:8)='noeuds'
  229.         IF (ipt7.num(/2).NE.ipt8.num(/2)) moterr(1:8)='elements'
  230.         SEGDES ipt7,ipt8
  231.         IF (moterr.NE.' ') THEN
  232.           CALL ERREUR(550)
  233.           RETURN
  234.         ENDIF
  235.       ENDDO
  236.  
  237.  100  SEGACT MCOORD*MOD
  238. C  Creation du tableau de correspondance icpr donnant pour chaque
  239. C  noeud de ipt2 son noeud image de ipt3
  240.       SEGACT ipt2,ipt3
  241.       SEGINI icpr
  242.       numnp=0
  243.       ipt7=ipt2
  244.       ipt8=ipt3
  245.       DO i=1,MAX(1,ipt2.lisous(/1))
  246.         IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(i)
  247.         IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(i)
  248.         SEGACT ipt7,ipt8
  249.         DO k=1,ipt7.num(/2)
  250.           DO l=1,ipt7.num(/1)
  251.             ip=ipt7.num(l,k)
  252.             IF (icpr(ip).EQ.0) THEN
  253.               numnp=numnp+1
  254.               icpr(ip)=ipt8.num(l,k)
  255.             ENDIF
  256.           ENDDO
  257.         ENDDO
  258.         SEGDES ipt7,ipt8
  259.       ENDDO
  260.       SEGDES ipt2,ipt3
  261.  
  262. C  Est-ce une transformation affine ?
  263.       IF (numnp.LE.idimp1) THEN
  264.         SEGINI xxtra
  265.         numnp=0
  266.         DO i=1,icpr(/1)
  267.           IF (icpr(i).NE.0) THEN
  268.             numnp=numnp+1
  269.             iref=(i-1)*idimp1
  270.             irefep=(icpr(i)-1)*idimp1
  271.             DO j=1,idim
  272.               xdep(j,numnp)=xcoor(iref+j)
  273.               xarr(j,numnp)=xcoor(irefep+j)
  274.             ENDDO
  275.           ENDIF
  276.         ENDDO
  277.         SEGSUP icpr
  278. C  Faut-il creer de nouveaux points ?
  279.         IF (numnp.EQ.1) THEN
  280. C  On en cree un par translation en x
  281.           numnp=numnp+1
  282.           DO j=1,IDIM
  283.             xdep(j,numnp)=xdep(j,numnp-1)
  284.             xarr(j,numnp)=xarr(j,numnp-1)
  285.           ENDDO
  286.           xdep(1,numnp)=xdep(1,numnp)+1.
  287.           xarr(1,numnp)=xarr(1,numnp)+1.
  288.         ENDIF
  289.         DO j=1,IDIM
  290.           ved1(j)=xdep(j,2)-xdep(j,1)
  291.           vea1(j)=xarr(j,2)-xarr(j,1)
  292.         ENDDO
  293.         IF (numnp.EQ.2) THEN
  294. C  En 1D, on ne fait rien.
  295. C  En 2D, on cherche un point sur la perpendiculaire a la droite
  296. C  definie par les deux premiers points, passant par le premier
  297. C  point et a une distance unite (formant un repere orthonorme direct).
  298. C  En 3D, on prend un point dans la direction y.
  299.           numnp=numnp+1
  300.           IF (IDIM.EQ.2) THEN
  301.             xdep(1,3)=ved1(2)
  302.             xdep(2,3)=-ved1(1)
  303.             aa=1./SQRT(xdep(1,3)*xdep(1,3)+xdep(2,3)*xdep(2,3))
  304.             xdep(1,3)=xdep(1,1)+xdep(1,3)*aa
  305.             xdep(2,3)=xdep(2,1)+xdep(2,3)*aa
  306.             xarr(1,3)=vea1(2)
  307.             xarr(2,3)=-vea1(2)
  308.             aa=1./SQRT(xarr(1,3)*xarr(1,3)+xarr(2,3)*xarr(2,3))
  309.             xarr(1,3)=xarr(1,1)+xarr(1,3)*aa
  310.             xarr(2,3)=xarr(2,1)+xarr(2,3)*aa
  311.           ELSE IF (IDIM.EQ.3) THEN
  312.             xdep(1,3)=xdep(1,1)
  313.             xdep(2,3)=xdep(2,1)+1.
  314.             xdep(3,3)=xdep(3,1)
  315.             xarr(1,3)=xarr(1,1)
  316.             xarr(2,3)=xarr(2,1)+1.
  317.             xarr(3,3)=xarr(3,1)
  318.           ENDIF
  319.         ENDIF
  320.         IF (IDIM.NE.1) THEN
  321.           DO j=1,IDIM
  322.             ved2(j)=xdep(j,3)-xdep(j,1)
  323.             vea2(j)=xarr(j,3)-xarr(j,1)
  324.           ENDDO
  325.         ENDIF
  326. C  En 3D, il faut encore un point hors du plan forme par les 3 points
  327.         IF (numnp.LT.idimp1) THEN
  328.           ved3(1)=ved1(2)*ved2(3)-ved1(3)*ved2(2)
  329.           ved3(2)=ved1(3)*ved2(1)-ved1(1)*ved2(3)
  330.           ved3(3)=ved1(1)*ved2(2)-ved1(2)*ved2(1)
  331.           aa=1./SQRT(ved3(1)*ved3(1)+ved3(2)*ved3(2)+ved3(3)*ved3(3))
  332.           vea3(1)=vea1(2)*vea2(3)-vea1(3)*vea2(2)
  333.           vea3(2)=vea1(3)*vea2(1)-vea1(1)*vea2(3)
  334.           vea3(3)=vea1(1)*vea2(2)-vea1(2)*vea2(1)
  335.           bb=1./SQRT(vea3(1)*vea3(1)+vea3(2)*vea3(2)+vea3(3)*vea3(3))
  336.           DO j=1,3
  337.             ved3(j)=ved3(j)*aa
  338.             vea3(j)=vea3(j)*bb
  339.             xdep(j,4)=xdep(j,1)+ved3(j)
  340.             xarr(j,4)=xarr(j,1)+vea3(j)
  341.           ENDDO
  342.         ENDIF
  343.         IF (IDIM.EQ.3) THEN
  344.           DO j=1,3
  345.             ved3(j)=xdep(j,4)-xdep(j,1)
  346.             vea3(j)=xarr(j,4)-xarr(j,1)
  347.           ENDDO
  348.         ENDIF
  349. C  Calcul de la matrice de passage pa du repere ved au repere vea
  350. C  On traite le cas IDIM=1 a part (sinon reso33 serait a modifier).
  351.         IF (IDIM.EQ.1) THEN
  352.           IF ((ved1(1).EQ.0.).OR.(vea1(1).EQ.0.)) THEN
  353.             CALL ERREUR(21)
  354.             SEGSUP xxtra
  355.             IF (ichpo.EQ.1) THEN
  356.               nbpts=nbpts0
  357.               SEGADJ mcoord
  358.               SEGSUP ipt2,ipt3
  359.             ENDIF
  360.             RETURN
  361.           ENDIF
  362.           pa(1,1)=vea1(1)/ved1(1)
  363.         ELSE
  364.           DO j=1,IDIM
  365.             vb(1)=vea1(j)
  366.             vb(2)=vea2(j)
  367.             IF (IDIM.EQ.3) vb(3)=vea3(j)
  368.             pa3(1,1)=ved1(1)
  369.             pa3(1,2)=ved1(2)
  370.             pa3(2,1)=ved2(1)
  371.             pa3(2,2)=ved2(2)
  372.             IF (IDIM.EQ.3) THEN
  373.               pa3(1,3)=ved1(3)
  374.               pa3(2,3)=ved2(3)
  375.               pa3(3,1)=ved3(1)
  376.               pa3(3,2)=ved3(2)
  377.               pa3(3,3)=ved3(3)
  378.             ENDIF
  379.             CALL reso33(pa3,vb,3,IDIM,kerr)
  380. C  Erreur : en 2D -> 3 points alignes, en 3D -> 4 points coplanaires
  381.             IF (kerr.NE.0) THEN
  382.               CALL ERREUR(21)
  383.               SEGSUP xxtra
  384.               IF (ichpo.EQ.1) THEN
  385.                 nbpts=nbpts0
  386.                 SEGADJ mcoord
  387.                 SEGSUP ipt2,ipt3
  388.               ENDIF
  389.               RETURN
  390.             ENDIF
  391.             pa(j,1)=vb(1)
  392.             pa(j,2)=vb(2)
  393.             IF (IDIM.EQ.3) pa(j,3)=vb(3)
  394.           ENDDO
  395.         ENDIF
  396. C  Fabrication de la liste des points a transformer et du maillage
  397. C  resultat (segment MELEME)
  398.         SEGINI icpr
  399.         nbpan=nbpts
  400.         SEGACT ipt1
  401.         SEGINI,ipt4=ipt1
  402.         nbnn=ipt4.num(/1)
  403.         nbelem=ipt4.num(/2)
  404.         nbsous=ipt4.lisous(/1)
  405.         nbref=ipt4.lisref(/1)
  406.         IF (ICAS.NE.1) THEN
  407.           nbref=0
  408.           SEGADJ ipt4
  409.         ENDIF
  410.         ipt8=ipt1
  411.         ipt5=ipt4
  412.         numnp=0
  413.         DO i=1,MAX(1,ipt1.lisous(/1))
  414.           IF (ipt1.lisous(/1).NE.0) THEN
  415.             ipt8=ipt1.lisous(i)
  416.             SEGACT ipt8
  417.             SEGINI,ipt5=ipt8
  418.             nbnn=ipt5.num(/1)
  419.             nbelem=ipt5.num(/2)
  420.             nbsous=ipt5.lisous(/1)
  421.             nbref=ipt5.lisref(/1)
  422.             IF (ICAS.NE.1) THEN
  423.               nbref=0
  424.               SEGADJ ipt5
  425.             ENDIF
  426.             ipt4.lisous(i)=ipt5
  427.           ENDIF
  428.           DO iel=1,ipt8.num(/2)
  429.             DO inn=1,ipt8.num(/1)
  430.               ipt=ipt8.num(inn,iel)
  431.               IF (icpr(ipt).EQ.0) THEN
  432.                 IF (ICAS.EQ.1) THEN
  433.                   icpr(ipt)=ipt
  434.                 ELSE
  435.                   numnp=numnp+1
  436.                   icpr(ipt)=numnp+nbpan
  437.                 ENDIF
  438.               ENDIF
  439.               ipt5.num(inn,iel)=icpr(ipt)
  440.             ENDDO
  441.           ENDDO
  442.           SEGDES ipt8,ipt5
  443.         ENDDO
  444.         SEGDES ipt4,ipt1
  445.         IF (ICAS.EQ.0) THEN
  446.           nbpts=numnp+nbpan
  447.           SEGADJ mcoord
  448.         ENDIF
  449.         DO ipt=1,nbpan
  450.           IF (icpr(ipt).NE.0) THEN
  451.             jno=(ipt-1)*idimp1
  452.             ino=(icpr(ipt)-1)*idimp1
  453.             DO i=1,IDIM
  454.               vb(i)=xcoor(jno+i)-xdep(i,1)
  455.             ENDDO
  456.             DO i=1,IDIM
  457.               aa=xarr(i,1)
  458.               DO j=1,IDIM
  459.                 aa=aa+pa(i,j)*vb(j)
  460.               ENDDO
  461.               xcoor(ino+i)=aa
  462.             ENDDO
  463.             xcoor(ino+idimp1)=xcoor(jno+idimp1)
  464.           ENDIF
  465.         ENDDO
  466.         SEGSUP icpr,xxtra
  467.         CALL ECROBJ('MAILLAGE',ipt4)
  468.         IF (ichpo.EQ.1) THEN
  469.           SEGSUP ipt2,ipt3
  470.         ENDIF
  471.         RETURN
  472.       ENDIF
  473.  
  474. C  Cas d'une transformation non affine
  475. C  Interdit dans le cas ICAS=1 (appel par deplac, operateur DEPL 'DEDU')
  476. C  On recupere le maillage initial si ichpo=1.
  477.       IF (ICAS.EQ.1) THEN
  478.         CALL ERREUR(21)
  479.         IF (ichpo.EQ.1) THEN
  480.           nbpts=nbpts0
  481.           SEGADJ mcoord
  482.           SEGSUP ipt2,ipt3
  483.         ENDIF
  484.         RETURN
  485.       ENDIF
  486. C  CHOIX : Pour l'instant, cas non implemente en DIMEnsion 1.
  487.       IF (IDIM.EQ.1) THEN
  488.         MOTERR(1:4)=mcle(iregu)
  489.         INTERR(1)=IDIM
  490.         CALL ERREUR(971)
  491.         IF (ichpo.EQ.1) THEN
  492.           nbpts=nbpts0
  493.           SEGADJ mcoord
  494.           SEGSUP ipt2,ipt3
  495.         ENDIF
  496.         RETURN
  497.       ENDIF
  498. C  Creation d'un tableau de stockage contenant pour chaque point les
  499. C  aretes des elements le touchant sauf celles qui le contiennent
  500.       SEGACT ipt1
  501.       SEGINI icpr
  502.       numnp=0
  503.       ipt8=ipt1
  504.       DO i=1,MAX(1,ipt1.lisous(/1))
  505.         IF (ipt1.lisous(/1).NE.0) ipt8=ipt1.lisous(i)
  506.         SEGACT ipt8
  507.         DO iel=1,ipt8.num(/2)
  508.           DO inn=1,ipt8.num(/1)
  509.             ipt=ipt8.num(inn,iel)
  510.             IF (icpr(ipt).EQ.0) THEN
  511.               numnp=numnp+1
  512.               icpr(ipt)=numnp
  513.             ENDIF
  514.           ENDDO
  515.         ENDDO
  516.       ENDDO
  517.       SEGINI idcp
  518.       DO i=1,icpr(/1)
  519.         IF (icpr(i).NE.0) idcp(icpr(i))=i
  520.       ENDDO
  521.       nkon=5
  522.       SEGINI kon,ikon
  523.       DO i=1,MAX(1,ipt1.lisous(/1))
  524.         IF (ipt1.lisous(/1).NE.0) ipt8=ipt1.lisous(i)
  525.         DO iel=1,ipt8.num(/2)
  526.           DO inn=1,ipt8.num(/1)
  527.             ip=ipt8.num(inn,iel)
  528.             ipc=icpr(ip)
  529.             DO 41 inp=1,ipt8.num(/1)
  530.               kpt1=ipt8.num(inp,iel)
  531.               IF (kpt1.EQ.ip) GOTO 41
  532.               inps=inp+1
  533.               IF (inps.GT.ipt8.num(/1)) inps=1
  534.               kpt2=ipt8.num(inps,iel)
  535.               IF (kpt2.EQ.ip) GOTO 41
  536.               DO idj=1,ikon(ipc)
  537.                 IF ((kon(idj,ipc,1).EQ.kpt1).AND.
  538.      .              (kon(idj,ipc,2).EQ.kpt2)) GOTO 41
  539.                 IF ((kon(idj,ipc,1).EQ.kpt2).AND.
  540.      .              (kon(idj,ipc,2).EQ.kpt1)) GOTO 41
  541.               ENDDO
  542.               ikon(ipc)=ikon(ipc)+1
  543.               IF (ikon(ipc).GT.nkon) THEN
  544.                 nkon=ikon(ipc)
  545.                 SEGADJ kon
  546.               ENDIF
  547.               kon(ikon(ipc),ipc,1)=kpt1
  548.               kon(ikon(ipc),ipc,2)=kpt2
  549.   41        CONTINUE
  550.           ENDDO
  551.         ENDDO
  552.       ENDDO
  553. C  Fermeture du contour si necessaire (avec 1 segment supplementaire)
  554.       DO 80 ipc=1,kon(/2)
  555.         iuni=0
  556.         DO iko=1,ikon(ipc)
  557.           kpt1=kon(iko,ipc,1)
  558.           kpt2=kon(iko,ipc,2)
  559.           DO ikq=1,ikon(ipc)
  560.             IF (ikq.NE.iko) THEN
  561.               IF (kpt1.EQ.kon(ikq,ipc,1)) GOTO 86
  562.               IF (kpt1.EQ.kon(ikq,ipc,2)) GOTO 86
  563.             ENDIF
  564.           ENDDO
  565. C  kpt1 unique
  566.           iuni=iuni+1
  567.           iuniq(iuni)=kpt1
  568.   86      DO 88 ikq=1,ikon(ipc)
  569.             IF (ikq.EQ.iko) GOTO 88
  570.             IF (kpt2.EQ.kon(ikq,ipc,1)) GOTO 90
  571.             IF (kpt2.EQ.kon(ikq,ipc,2)) GOTO 90
  572.   88      CONTINUE
  573. C  kpt1 unique
  574.           iuni=iuni+1
  575.           iuniq(iuni)=kpt2
  576.   90      CONTINUE
  577.         ENDDO
  578.         IF (iuni.EQ.0) GOTO 80
  579.         ikon(ipc)=ikon(ipc)+1
  580.         IF (ikon(ipc).GT.nkon) THEN
  581.           nkon=ikon(ipc)
  582.           SEGADJ kon
  583.         ENDIF
  584.         kon(ikon(ipc),ipc,1)=iuniq(1)
  585.         kon(ikon(ipc),ipc,2)=iuniq(2)
  586.   80  CONTINUE
  587. C  Calcul des coordonnees barycentriques
  588.       SEGINI xkon,ipxkon
  589.       DO i=1,kon(/2)
  590.         ipcour=0
  591.         jp=idcp(i)
  592.         ipcini=ipcour
  593.   58    CONTINUE
  594.         DO ipc=ipcini+1,ipcour
  595.           ipxkon(ipc,i)=0
  596.           xkon(ipc,i)=0.
  597.         ENDDO
  598.         ipcour=ipcini
  599.         DO 59 j=1,ikon(i)
  600. C  Creation d'un angle
  601.           IF (kon(j,i,1).EQ.0) GOTO 59
  602.           ip1=kon(j,i,1)
  603.           ip2=kon(j,i,2)
  604.           DO 60 k=j+1,ikon(i)
  605.             IF (kon(k,i,1).EQ.0) GOTO 60
  606.             jp1=0
  607.             IF (ip1.EQ.kon(k,i,1)) THEN
  608.               jp1=ip1
  609.               jp2=ip2
  610.               jp3=kon(k,i,2)
  611.             ELSE IF (ip1.EQ.kon(k,i,2)) THEN
  612.               jp1=ip1
  613.               jp2=ip2
  614.               jp3=kon(k,i,1)
  615.             ELSE IF (ip2.EQ.kon(k,i,1)) THEN
  616.               jp1=ip2
  617.               jp2=ip1
  618.               jp3=kon(k,i,2)
  619.             ELSE IF (ip2.EQ.kon(k,i,2)) THEN
  620.               jp1=ip2
  621.               jp2=ip1
  622.               jp3=kon(k,i,1)
  623.             ELSE
  624.               GOTO 60
  625.             ENDIF
  626. *  jp le point en cours de travail
  627. *  jp1 le sommet jp2 et jp3 les deux autres points
  628.             xp=xcoor((jp-1)*idimp1+1)
  629.             yp=xcoor((jp-1)*idimp1+2)
  630.             zp=0.
  631.             IF (IDIM.EQ.3) zp=xcoor((jp-1)*idimp1+3)
  632.             xp1=xcoor((jp1-1)*idimp1+1)
  633.             yp1=xcoor((jp1-1)*idimp1+2)
  634.             zp1=0.
  635.             IF (IDIM.EQ.3) zp1=xcoor((jp1-1)*idimp1+3)
  636.             xp2=xcoor((jp2-1)*idimp1+1)
  637.             yp2=xcoor((jp2-1)*idimp1+2)
  638.             zp2=0.
  639.             IF (IDIM.EQ.3) zp2=xcoor((jp2-1)*idimp1+3)
  640.             xp3=xcoor((jp3-1)*idimp1+1)
  641.             yp3=xcoor((jp3-1)*idimp1+2)
  642.             zp3=0.
  643.             IF (IDIM.EQ.3) zp3=xcoor((jp3-1)*idimp1+3)
  644. *  le  sinus (trafique)
  645.             xv1=xp2-xp1
  646.             yv1=yp2-yp1
  647.             zv1=zp2-zp1
  648.  
  649.             aa=SQRT(xv1*xv1+yv1*yv1+zv1*zv1)
  650.             IF (aa .LT. XPETIT) THEN
  651.               xv1=0.D0
  652.               yv1=0.D0
  653.               zv1=0.D0
  654.             ELSE
  655.               xv1=xv1/aa
  656.               yv1=yv1/aa
  657.               zv1=zv1/aa
  658.             ENDIF
  659.             xv2=xp3-xp1
  660.             yv2=yp3-yp1
  661.             zv2=zp3-zp1
  662.  
  663.             aa=SQRT(xv2*xv2+yv2*yv2+zv2*zv2)
  664.             IF (aa .LT. XPETIT) THEN
  665.               xv2=0.D0
  666.               yv2=0.D0
  667.               zv2=0.D0
  668.             ELSE
  669.               xv2=xv2/aa
  670.               yv2=yv2/aa
  671.               zv2=zv2/aa
  672.             ENDIF
  673.             xv=yv1*zv2-zv1*yv2
  674.             yv=zv1*xv2-xv1*zv2
  675.             zv=xv1*yv2-yv1*xv2
  676.  
  677. *  provisoirement on prend le signe positif
  678.             sinp=MAX(1.D-15,SQRT(xv*xv+yv*yv+zv*zv))
  679. *  recherche du signe comparaison avec jp jp1 vect jp2 jp3
  680.             xw1=xp-xp1
  681.             yw1=yp-yp1
  682.             zw1=zp-zp1
  683.             aa=SQRT(xw1*xw1+yw1*yw1+zw1*zw1)
  684.             IF (aa .LT. XPETIT) THEN
  685.               xw1=0.D0
  686.               yw1=0.D0
  687.               zw1=0.D0
  688.             ELSE
  689.               xw1=xw1/aa
  690.               yw1=yw1/aa
  691.               zw1=zw1/aa
  692.             ENDIF
  693.             xw2=xp3-xp2
  694.             yw2=yp3-yp2
  695.             zw2=zp3-zp2
  696.  
  697.             aa=SQRT(xw2*xw2+yw2*yw2+zw2*zw2)
  698.             IF (aa .LT. XPETIT) THEN
  699.               xw2=0.D0
  700.               yw2=0.D0
  701.               zw2=0.D0
  702.             ELSE
  703.               xw2=xw2/aa
  704.               yw2=yw2/aa
  705.               zw2=zw2/aa
  706.             ENDIF
  707.             xw=yw1*zw2-zw1*yw2
  708.             yw=zw1*xw2-xw1*zw2
  709.             zw=xw1*yw2-yw1*xw2
  710.  
  711.             IF (iregu.NE.2) THEN
  712.               IF ((xv*xw+yv*yw+zv*zw).LT.0.) sinp=-sinp
  713.             ELSE
  714.               IF ((xv*xw+yv*yw+zv*zw).LT.-1.D-15) then
  715. * Essai si le signe est negatif on supprime le coin pour ne travailler
  716. * que dans un convexe
  717.                 kon(k,i,1)=0
  718.                 kon(k,i,2)=0
  719.                 kon(j,i,1)=jp2
  720.                 kon(j,i,2)=jp3
  721.                 GOTO 58
  722.               ENDIF
  723.             ENDIF
  724.             ipcour=ipcour+1
  725.             ipxkon(ipcour,i)=jp1
  726. *  les distances aux cotes
  727.             xv=xp-xp1
  728.             yv=yp-yp1
  729.             zv=zp-zp1
  730.             scal=xv*xv1+yv*yv1+zv*zv1
  731.             xv=xv-scal*xv1
  732.             yv=yv-scal*yv1
  733.             zv=zv-scal*zv1
  734.             dis1=MAX(1.D-20,SQRT(xv*xv+yv*yv+zv*zv))
  735.             xv=xp-xp1
  736.             yv=yp-yp1
  737.             zv=zp-zp1
  738.             scal=xv*xv2+yv*yv2+zv*zv2
  739.             xv=xv-scal*xv2
  740.             yv=yv-scal*yv2
  741.             zv=zv-scal*zv2
  742.             dis2=MAX(1.D-20,SQRT(xv*xv+yv*yv+zv*zv))
  743. C  la coordonnee barycentrique :
  744.             xkon(ipcour,i)=sinp/(dis1*dis2)
  745.             IF (iregu.EQ.1) xkon(ipcour,i)=1.D0
  746.   60      CONTINUE
  747.   59    CONTINUE
  748. *  on normalise
  749.         skon=0.
  750.         DO j=1,ikon(i)+1
  751.           skon=skon+xkon(j,i)
  752.         ENDDO
  753.         DO j=1,ikon(i)+1
  754.         xkon(j,i)=xkon(j,i)/skon
  755.         ENDDO
  756.       ENDDO
  757.       SEGSUP kon,ikon
  758. *  OK on a les coefficients barycentriques.
  759. *  On duplique le maillage puis on cree les nouveaux noeuds
  760. *  Mais d'abord voir les points imposes
  761.       SEGACT ipt2,ipt3
  762.       SEGINI icp2
  763.       ipt7=ipt2
  764.       ipt8=ipt3
  765.       DO isous=1,MAX(1,ipt2.lisous(/1))
  766.         IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(isous)
  767.         IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(isous)
  768.         SEGACT ipt7,ipt8
  769.         DO i=1,ipt7.num(/1)
  770.           DO j=1,ipt7.num(/2)
  771.             ip=ipt7.num(i,j)
  772.             IF (icp2(ip).EQ.0) icp2(ip)=ipt8.num(i,j)
  773.           ENDDO
  774.         ENDDO
  775.         SEGDES ipt7,ipt8
  776.       ENDDO
  777.       SEGDES ipt2,ipt3
  778.       SEGINI icorr
  779.       DO 210 i=1,icpr(/1)
  780.         IF (icpr(i).EQ.0) GOTO 210
  781.         IF (icp2(i).EQ.0) THEN
  782.           nbpts=nbpts+1
  783.           SEGADJ mcoord
  784.           icorr(icpr(i))=nbpts
  785.           DO j=1,idimp1
  786.             xcoor((nbpts-1)*idimp1+j)=xcoor((i-1)*idimp1+j)
  787.           ENDDO
  788.         ELSE
  789.           icorr(icpr(i))=icp2(i)
  790.         ENDIF
  791.  210  CONTINUE
  792. *  On duplique les maillages
  793.       SEGINI,ipt4=ipt1
  794.       nbref=0
  795.       nbsous=ipt4.lisous(/1)
  796.       nbnn=ipt4.num(/1)
  797.       nbelem=ipt4.num(/2)
  798.       SEGADJ ipt4
  799.       ipt7=ipt1
  800.       ipt8=ipt4
  801.       DO isous=1,MAX(1,ipt1.lisous(/1))
  802.         IF (ipt1.lisous(/1).NE.0) THEN
  803.           ipt7=ipt1.lisous(isous)
  804.           SEGINI,ipt8=ipt7
  805.           ipt4.lisous(isous)=ipt8
  806.         ENDIF
  807.         DO i=1,ipt8.num(/1)
  808.           DO j=1,ipt8.num(/2)
  809.             ipt8.num(i,j)=icorr(icpr(ipt8.num(i,j)))
  810.           ENDDO
  811.         ENDDO
  812.         SEGDES ipt7,ipt8
  813.       ENDDO
  814.       SEGDES ipt1,ipt4
  815. C  Maintenant on peut deplacer les points
  816.       nfois=50
  817.       CALL LIRENT(nfois,0,iretou)
  818.       DO ifois=1,nfois
  819.         DO i=1,numnp
  820.           IF (icp2(idcp(i)).EQ.0) THEN
  821.             DO j=1,IDIM
  822.               xaux=0.
  823.               DO k=1,nkon+1
  824.                 IF (ipxkon(k,i).NE.0) xaux=xaux+xkon(k,i)
  825.      .                  *xcoor((icorr(icpr(ipxkon(k,i)))-1)*idimp1+j)
  826.               ENDDO
  827.               xcoor((icorr(i)-1)*idimp1+j)=xaux
  828.             ENDDO
  829.           ENDIF
  830.         ENDDO
  831.       ENDDO
  832.  
  833.       IF (ichpo.EQ.0) THEN
  834.          CALL ECROBJ('MAILLAGE',ipt4)
  835.       ELSE
  836. C  Sortie d'un CHPOINT, les points concernes sont reperes dans icpr
  837.         SEGSUP ipt2,ipt3
  838.         nbelem=icorr(/1)
  839.         nbnn=1
  840.         nbsous=0
  841.         nbref=0
  842.         SEGINI ipt2
  843.         ipt2.itypel=1
  844.         nat=1
  845.         nsoupo=1
  846.         SEGINI mchpo1
  847.         mchpo1.jattri(1)=jatan
  848.         IF (IDIM.EQ.1) nc=1
  849.         IF (IDIM.EQ.2) nc=2
  850.         IF (IDIM.EQ.3) nc=3
  851.         SEGINI msoupo
  852.         mchpo1.ipchp(1)=msoupo
  853.         igeoc=ipt2
  854.         n=nbelem
  855.         SEGINI mpoval
  856.         ipoval=mpoval
  857.         IF (IDIM.EQ.1) THEN
  858.           IF (NIFOUR.LE.9) THEN
  859.             nocomp(1)='UX'
  860.           ELSE
  861.             nocomp(1)='UR'
  862.           ENDIF
  863.         ELSE IF (IDIM.EQ.2) THEN
  864.           IF (IFOUR.GE.0) THEN
  865.             nocomp(1)='UR'
  866.             nocomp(2)='UZ'
  867.           ELSE
  868.             nocomp(1)='UX'
  869.             nocomp(2)='UY'
  870.           ENDIF
  871.         ELSE IF (IDIM.EQ.3) THEN
  872.           nocomp(1)='UX'
  873.           nocomp(2)='UY'
  874.           nocomp(3)='UZ'
  875.         ENDIF
  876.         SEGDES msoupo
  877.         DO i=1,idcp(/1)
  878.           ipt2.num(1,i)=idcp(i)
  879.           ifi=(icorr(i)-1)*idimp1
  880.           ide=(idcp(i)-1)*idimp1
  881.           DO j=1,IDIM
  882.             vpocha(i,j)=xcoor(ifi+j)-xcoor(ide+j)
  883.           ENDDO
  884.         ENDDO
  885.         SEGDES mpoval,ipt2,mchpo1
  886.         CALL ECROBJ('CHPOINT ',mchpo1)
  887.       ENDIF
  888.       SEGSUP icpr,idcp,xkon,ipxkon,icp2,icorr
  889.  
  890.       RETURN
  891.       END
  892.  
  893.  
  894.  
  895.  
  896.  
  897.  
  898.  

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