Télécharger la plaquette Cast3M


fr en

Télécharger rot3m.eso

Retour à la liste

Numérotation des lignes : 

rot3m

  1. C ROT3M     SOURCE    OF166741  25/02/21    21:18:22     12166          
  2.  
  3. ************************************************************************
  4. *
  5. *                             R O T 3 M
  6. *                             ---------
  7. *
  8. * FONCTION:
  9. * ---------
  10. *     CALCUL DE LA MATRICE DE MUTUELLES POUR L'ELEMENT ROT3
  11. *
  12. * PARAMETRES:   (E)=ENTREE   (S)=SORTIE   (+ = CONTENU DANS UN COMMUN)
  13. * -----------
  14. *
  15. *     NEF     (E)  NUMERO DE L'ELEMENT-FINI DANS NOMTP (VOIR CCHAMP)
  16. *     IMAIL   (E)  NUMERO DU MAILLAGE ELEMENTAIRE CONSIDERE,DANS
  17. *                  L'OBJET MODELE
  18. *     IPMODE  (E)  POINTEUR SUR UN SEGMENT IMODEL
  19. *     IPCHEM  (E)  POINTEUR SUR LE CHAMELEM DE CARACTERISTIQUE
  20. *     IPSUPJ  (E)  POINTEUR SUR LE MAILLAGE SUPPORT DES COURANTS DE FOUCAULT
  21. *     IPRIGI (E/S) POINTEUR SUR L'OBJET RESULTAT,DE TYPE RIGIDITE
  22. *
  23. * VARIABLES:
  24. * ----------
  25. *
  26. *     NBNN     NOMBRE DE NOEUDS DANS L'ELEMENT CONSIDERE
  27. *     NEF      NUMERO DE L'ELEMENT FINI DANS NOMTP (VOIR CCHAMP)
  28. *     NBELEM   NOMBRE D'ELEMENTS DANS LE MAILLAGE ELEMENTAIRE
  29. *     NBPGAU   NOMBRE DE POINTS DE GAUSS DANS L'ELEMENT-FINI
  30. *     NDIM     NOMBRE DE LIGNES DE LA MATRICE GRADIENT
  31. *     XE(3,NBNN)      COORDONNEES DE L'ELEMENT DANS LE REPERE GLOBAL
  32. *     SHP(6,NBNN)     TABLEAU DE TRAVAIL
  33. *     GRAD(NDIM,NBNN) MATRICE GRADIENT DES FONCTIONS DE FORME BIDIM.
  34. *     VALMAT(NMATR)  TABLEAU DE TRAVAIL
  35. *
  36. * AUTEUR, DATE DE CREATION:
  37. * -------------------------
  38. *
  39. *     YANN STEPHAN , FEVRIER 1997 (COPIE DE ROT3R)
  40. *
  41. ************************************************************************
  42.  
  43.       SUBROUTINE ROT3M(NEF,IMAIL,IPMODE,IPCHEM,IPSUPJ,XMATRI,
  44.      $                 ICPR,ICPR2)
  45.  
  46.       IMPLICIT INTEGER(I-N)
  47.       IMPLICIT REAL*8(A-H,O-Z)
  48.  
  49. -INC PPARAM
  50. -INC CCOPTIO
  51. -INC CCREEL
  52. -INC CCHAMP
  53.  
  54. -INC SMCOORD
  55. -INC SMINTE
  56. -INC SMMODEL
  57. -INC SMRIGID
  58. -INC SMELEME
  59. -INC SMCHAML
  60.  
  61. -INC TMPTVAL
  62.  
  63.       SEGMENT MAXT
  64.         REAL*8 RA(NBNN,NBNN)
  65.       ENDSEGMENT
  66.  
  67.       SEGMENT,MMAT1
  68.         REAL*8 VALMAT(NMATR)
  69.         REAL*8 XE(3,NBNN),XE1(3,NBNN)
  70.         REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN,NBPGAU)
  71.         REAL*8 COSD1(3),COSD2(3)
  72.       ENDSEGMENT
  73.       POINTEUR MMAT2.MMAT1,MMATX.MMAT1
  74.  
  75.       SEGMENT NOTYPE
  76.         CHARACTER*16 TYPE(NBTYPE)
  77.       ENDSEGMENT
  78.  
  79.       SEGMENT SGAUSS
  80.         REAL*8 XGAUSS(3,NBPGAU)
  81.         REAL*8 DX(NBPGAU)
  82.       ENDSEGMENT
  83.       POINTEUR SGX.SGAUSS,SGY.SGAUSS
  84.  
  85.       SEGMENT ICPR(NA)
  86.       SEGMENT ICPR2(NA)
  87.  
  88.       CHARACTER*8 CNM
  89.       CHARACTER*(NCONCH) CONM
  90.       PARAMETER (NINF=3)
  91.       INTEGER INFOS(NINF)
  92.       LOGICAL SELF,NEAR
  93.  
  94. *     PERMEABILITE DU VIDE SUR 4PI
  95.       DATA PM0S4P/1.D-7/
  96.  
  97.       IMODEL = IPMODE
  98.  
  99.       CONM   = imodel.CONMOD
  100.       IPMAIL = imodel.IMAMOD
  101.       CNM    = imodel.CMATEE
  102.  
  103. *     RECUPERATION DES CARACTERISTIQUES GEOMETRIQUES DU MAILLAGE
  104. *     ELEMENTAIRE
  105.  
  106.       MELEME = IPMAIL
  107.       NBNN   = meleme.NUM(/1)
  108.       NBELEM = meleme.NUM(/2)
  109. *
  110. *     REMLIR LE TABLEAU INFOS (INFORMATIONS SUR ELEMENT)
  111.       INFOS(1)=0
  112.       INFOS(2)=0
  113.       INFOS(3)=NIFOUR
  114.  
  115. *     RECUPERATION DES CARACTERISTIQUES D'INTEGRATION
  116. *     POUR LA MATRICE MUTUELLE (RIGIDITE) DE L'ELEMENT
  117. *     FINI LIE A NOTRE MAILLAGE
  118.  
  119.       if (infmod(/1).lt.5) then
  120.         write(ioimp,*) 'ROT3 - INFMOD(/1) < 5'
  121.         call erreur(5)
  122.       endif
  123.       MINTE  = imodel.INFMOD(5)
  124.       NBPGAU = minte.POIGAU(/1)
  125.  
  126. *
  127. *     RECHERCHE LES POINTEURS DES SEGMENTS MELVAL QUI CORRESPONDENT
  128. *     A LA PERMEABILITE ET L'EPAISSEUR DES ELEMENTS
  129. *
  130.       nomid  = 0
  131.       nbrobl = 0
  132.       nbrfac = 0
  133.       IF (CNM.EQ.'ISOTROPE'.OR.CNM.EQ.'ORTHOTRO') THEN
  134.         NBROBL=1
  135.         SEGINI,NOMID
  136.         LESOBL(1)='PERM'
  137.       ELSE
  138.         CALL ERREUR(251)
  139.         RETURN
  140.       ENDIF
  141.       NMATO = nbrobl
  142.       NMATF = nbrfac
  143.       NMATR = NMATO+NMATF
  144.       MOMATR = nomid
  145.  
  146.       NBTYPE=1
  147.       SEGINI,NOTYPE
  148.       TYPE(1) ='REAL*8'
  149.       MOTYPE = NOTYPE
  150.  
  151.       CALL KOMCHA(IPCHEM,IPMAIL,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT)
  152.       IF (IERR.NE.0) RETURN
  153.  
  154.       MPTVAL = IVAMAT
  155.  
  156.       SEGINI SGX,SGY
  157.  
  158.       NDIM  = IDIM
  159.       NDIM2 = IDIM-1
  160.       SEGINI,MAXT,MMAT1,MMAT2
  161.       MMATX = MMAT1
  162. *
  163. *     CALCUL DE LA DISTANCE DE REFERENCE
  164. *
  165.       DREF = 0.D0
  166.       DO IEL = 1, NBELEM
  167.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  168.         CALL MAXLT3(NBNN,XE,DARET)
  169.         DREF = MAX(DREF,DARET)
  170.       ENDDO
  171.  
  172.       IPT1 = IPSUPJ
  173.       NBSOUJ = IPT1.LISOUS(/1)
  174.       NBSOU1 = MAX(1,NBSOUJ)
  175.  
  176. *     BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL
  177.  
  178.       DO 10 IEL=1,NBELEM
  179.  
  180.          MMAT1=MMATX
  181. *
  182. *        ON CHERCHE LES COORDONNEES DES NOEUDS DE L'ELEMENT IEL,
  183. *        DANS LE REPERE GLOBAL
  184. *
  185.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  186. *
  187. *        CALCUL DES COORDONNEES DES NOEUDS DANS LE REPERE LOCAL DE L'
  188. *        ELEMENT COQUE
  189. *
  190.          CALL COQLOC(NBNN,XE,COSD1,COSD2,XE1)
  191. *
  192. *        ON CALCULE LES FONCTIONS DE FORME AUX POINTS DE GAUSS
  193. *
  194.          CALL ELGAUS(MINTE,MMAT1,SGX,IFOIS,IFOI2)
  195. *
  196. *      LE JACOBIEN EST NEGATIF ,MAILLAGE INCORRECT
  197.         IF(IFOIS.NE.0.AND.IFOIS.NE.NBPGAU)THEN
  198.           INTERR(1)=IEL
  199.           CALL ERREUR(195)
  200.           GO TO 999
  201. *
  202. *      CAS OU LE JACOBIEN EST TRES PETIT
  203.         ELSEIF(IFOI2.EQ.NBPGAU)THEN
  204.           INTERR(1)=IEL
  205.           CALL ERREUR (259)
  206.           GO TO 999
  207.         ENDIF
  208. *
  209. *        ON BOUCLE SUR LE MAILLAGE SUPPORT DE COURANTS
  210.         IPT1 = IPSUPJ
  211.         IPT2 = IPT1
  212.         DO 110 ISOUJ = 1, NBSOU1
  213.           IF (NBSOUJ.NE.0) IPT2 = IPT1.LISOUS(ISOUJ)
  214.           NBELJ=IPT2.NUM(/2)
  215.           NBNNJ=IPT2.NUM(/1)
  216.           NBNNT=NBNN+NBNNJ
  217.           NLIGRP=NBNN
  218.           NLIGRD=NBNN
  219.           DO 111 IELJ=1,NBELJ
  220.  
  221.             DO IX=1,NBNN
  222.               DO JX=1,NBNN
  223.                 RA(JX,IX) = 0.D0
  224.               ENDDO
  225.             ENDDO
  226. *
  227.             NEAR=.FALSE.
  228. *
  229.             MMAT1=MMAT2
  230.             SGAUSS = SGY
  231. *
  232. *            ON CHERCHE LES COORDONNEES DES NOEUDS DE L'ELEMENT IEL,
  233. *            DANS LE REPERE GLOBAL
  234. *
  235.              CALL DOXE(XCOOR,IDIM,NBNN,IPT2.NUM,IELJ,XE)
  236. *
  237. *            CALCUL DES COORDONNEES DES NOEUDS DANS LE REPERE LOCAL DE L'
  238. *            ELEMENT COQUE
  239. *
  240.              CALL COQLOC(NBNN,XE,COSD1,COSD2,XE1)
  241. *
  242.              CALL ELGAUS(MINTE,MMAT1,SGAUSS,JFOIS,JFOI2)
  243. *
  244.              IF (JFOIS.NE.0.AND.JFOIS.NE.NBPGAU)THEN
  245. *
  246. *            LE JACOBIEN EST NEGATIF ,MAILLAGE INCORRECT
  247.                INTERR(1)=IEL
  248.                CALL ERREUR(195)
  249.                GO TO 999
  250.              ELSEIF(JFOI2.EQ.NBPGAU)THEN
  251. *
  252. *              CAS OU LE JACOBIEN EST TRES PETIT
  253. *
  254.                INTERR(1)=IEL
  255.                CALL ERREUR (259)
  256.                GO TO 999
  257.              ENDIF
  258. *
  259. *            CALCUL DE LA DISTANCE ENTRE LES DEUX ELEMENTS
  260. *
  261.              CALL S2ELT3(NBNN,NUM,IEL,IPT2.NUM,IELJ,NEAR,SELF)
  262.              CALL D2ELT3(NBNN,XE,MMATX.XE,DT3)
  263.              NEAR=NEAR.OR.(DT3.LE.DREF)
  264. *
  265. *            BOUCLE SUR LES POINTS DE GAUSS MAILLAGE 1
  266. *
  267.              DO 22 IGAU=1,NBPGAU
  268. *
  269.              MMAT1=MMATX
  270.              SGAUSS=SGX
  271. *
  272. *            ON CHERCHE LES VALEURS DE LA PERMEABILITE
  273. *
  274.              MPTVAL=IVAMAT
  275.              DO  30 IM=1,NMATR
  276.              IF(IVAL(IM).NE.0)THEN
  277.                MELVAL=IVAL(IM)
  278.                IBMN=MIN(IEL,VELCHE(/2))
  279.                IGMN=MIN(IGAU,VELCHE(/1))
  280.                VALMAT(IM)=VELCHE(IGMN,IBMN)
  281.              ELSE
  282.               VALMAT(IM)=0
  283.              ENDIF
  284.    30        CONTINUE
  285.              PERM=VALMAT(1)
  286.              XK=PERM*PM0S4P*DX(IGAU)
  287. *
  288. *            BOUCLE SUR LES POINTS DE GAUSS MAILLAGE 2
  289. *
  290.              DO 23 JGAU=1,NBPGAU
  291. *
  292.                MMAT1 = MMAT2
  293.                SGAUSS = SGY
  294. *
  295.                IF (SELF) THEN
  296.                  IF (JGAU.GT.1) GOTO 23
  297.                  CALL SELFT3(SGX.XGAUSS(1,IGAU),NBNN,XE,QQ)
  298.                  YK=XK*QQ
  299.                ELSE IF (NEAR) THEN
  300.                  IF(JGAU.GT.1) GO TO 23
  301.                  CALL NEART3(SGX.XGAUSS(1,IGAU),NBNN,MMATX.XE,XE,QQ)
  302.                  YK=XK*QQ
  303.                ELSE
  304.                  DIST=0.D0
  305.                  DO I=1,IDIM
  306.                    DIST=DIST+(SGX.XGAUSS(I,IGAU)-XGAUSS(I,JGAU))**2
  307.                  ENDDO
  308.                  DIST=SQRT(DIST)
  309.                  YK=XK*DX(JGAU)/DIST
  310.                ENDIF
  311. *
  312. *        ON AJOUTE LE PRODUIT K*DJAC*TRANSPOSEE(GRADX)*GRADY
  313. *        POUR LE POINT DE GAUSS CONSIDERE,A LA MATRICE RE
  314. *
  315.                MMAT1=MMATX
  316.                CALL WXWYSR(GRAD(1,1,IGAU),MMAT2.GRAD(1,1,JGAU),
  317.      &                     YK,NBNN,NBNNJ,IDIM,RA)
  318.    23        CONTINUE
  319.  
  320.    22      CONTINUE
  321.  
  322. *  realisation de l'assemblage
  323.        DO IX = 1, IPT2.NUM(/1)
  324.          IA = IPT2.NUM(IX,IELJ)
  325.          IB = ICPR2(IA)
  326.          DO JX = 1, NUM(/1)
  327.            IC = NUM(JX,IEL)
  328.            ID = ICPR(IC)
  329.            RE(IB,ID,1) = RA(IX,JX) + RE(IB,ID,1)
  330.          ENDDO
  331.        ENDDO
  332. *
  333.   111    CONTINUE
  334.   110  CONTINUE
  335. *
  336.  10   CONTINUE
  337. *     END DO
  338.  
  339. *     on symetrise la matrice
  340.       DO IX = 1, RE(/1)
  341.         DO JX = 1, IX
  342.           XP = 0.5D0 * ( RE(IX,JX,1) + RE(JX,IX,1) )
  343.           RE(IX,JX,1)=XP
  344.           RE(JX,IX,1)=XP
  345.         ENDDO
  346.       ENDDO
  347. *
  348. *     DESACTIVATION DES SEGMENTS
  349. *
  350.  999  CONTINUE
  351.       MMAT1=MMATX
  352.       SEGSUP,MMAT1,MMAT2,MAXT
  353.       MPTVAL=IVAMAT
  354.       SEGSUP,MPTVAL
  355.       NOMID=MOMATR
  356.       SEGSUP,NOMID,NOTYPE
  357.       SEGSUP,SGX,SGY
  358.  
  359. C      RETURN
  360.       END
  361.  
  362.  
  363.  

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