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

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