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

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