Télécharger la plaquette Cast3M


fr en

Télécharger mhybr2.eso

Retour à la liste

Numérotation des lignes : 

mhybr2

  1. C MHYBR2    SOURCE    CB215821  24/04/12    21:16:42     11897          
  2.       SUBROUTINE MHYBR2(IMAIL,IPMODE,IPCHEM,IPRIGI,IPGEOS,ILUMP)
  3. C-----------------------------------------------------------------------
  4. C Calcul de l'inverse de la matrice masse hybride.
  5. C Traitement du cas des elements finis massifs a integration numerique
  6. C pour un maillage elementaire et une formulation hybride.
  7. C-----------------------------------------------------------------------
  8. C
  9. C---------------------------
  10. C Parametres Entree/Sortie :
  11. C---------------------------
  12. C
  13. C E/  IMAIL   : Numero du maillage elementaire considere,
  14. C               dans l'objet modele.
  15. C E/  IPMODE  : Pointeur sur un segment IMODEL.
  16. C E/  IPCHEM  : Pointeur sur un chamelem de caracteristiques.
  17. C E/  IPGEOS  : Pointeur sur le maillage sommet
  18. C E/S IPRIGI  : Pointeur sur l'objet rigidite resultat.
  19. C
  20. C----------------------
  21. C Variables en COMMON :
  22. C----------------------
  23. C
  24. C E/  XCOOR   : VOIR SMCOORD
  25. C E/  IERR    : VOIR CCOPTIO
  26. C E/  IDIM    : VOIR CCOPTIO
  27. C E/  INTERR  : VOIR CCOPTIO
  28. C E/  IFOMOD  : VOIR CCOPTIO
  29. C E/  XPETIT  : VOIR CCREEL
  30. C
  31. C----------------------
  32. C Tableaux de travail :
  33. C----------------------
  34. C
  35. C     NBNN    : Nombre de noeuds dans l'element considere
  36. C     NEFHYB  : Numero de l'element fini dans NOMTP.
  37. C     NEF     : Numero de l'element fini support geometrique
  38. C               dans NOMTP (voir CCHAMP)
  39. C     NBELEM  : Nombre d'element dans le maillage elementaire
  40. C     NBPGAU  : Nombre de points de gauss pour l'element fini NEF
  41. C     CEL     : Matrice de conductivite elementaire
  42. C     XE      : Coordonnees des noeuds dans le repere global
  43. C     CMAT    : Matrice de permeabilite dans le repere global
  44. C     SHP     : Tableau de travail contenant les fonctions de forme au
  45. C               point de gauss utilise + derivees
  46. C     SHY     : Contient les fonctions de base hybride en un point
  47. C               mais pas les derivees de la fonction de base.
  48. C     VALMAT  : Valeurs des coeff. de la matrice CMAT et des
  49. C               cosinus directeurs des axes d'ortho. / repere local
  50. C     XGLOB   : Cosinus directeurs des axes d'ortho. / repere global
  51. C     XLOC    : Cosinus directeurs des axes d'ortho. / repere local
  52. C     TXR     : Cosinus directeurs des axes locaux / repere global
  53. C
  54. C
  55. C-----------------------------------------------------------------------
  56. C
  57. C Langage : ESOPE + FORTRAN77
  58. C
  59. C Auteurs : F.DABBENE 08/93
  60. C
  61. C-----------------------------------------------------------------------
  62.       IMPLICIT INTEGER(I-N)
  63.       IMPLICIT REAL*8(A-H,O-Z)
  64. *
  65. -INC CCHAMP
  66.  
  67. -INC PPARAM
  68. -INC CCOPTIO
  69. -INC CCREEL
  70. -INC SMCOORD
  71. -INC SMINTE
  72. -INC SMMODEL
  73. -INC SMRIGID
  74. -INC SMELEME
  75. -INC SMCHAML
  76. *
  77.       SEGMENT MAXE
  78.          REAL*8 TXR(IDIM,IDIM),XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM)
  79.       ENDSEGMENT
  80. *
  81.       SEGMENT MMAT1
  82.          REAL*8  CEL(NBDDL,NBDDL),CEL1(NBDDL,NBDDL),XE(3,NBNN)
  83.          REAL*8  SHP(6,NBNN),SHY(IDIM,NBDDL),ZJAC(IDIM,IDIM)
  84.          REAL*8  CMAT(IDIM,IDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM)
  85.          INTEGER ICSTO(NBDDL)
  86.       ENDSEGMENT
  87. *
  88.       SEGMENT NOTYPE
  89.          CHARACTER*16 TYPE(NBTYPE)
  90.       ENDSEGMENT
  91. *
  92.       SEGMENT MVELCH
  93.          REAL*8 VALMAT(NV1)
  94.       ENDSEGMENT
  95. *
  96.       SEGMENT MPTVAL
  97.          INTEGER      IPOS(NS) , NSOF(NS)
  98.          INTEGER      IVAL(NCOSOU)
  99.          CHARACTER*16 TYVAL(NCOSOU)
  100.       ENDSEGMENT
  101. *
  102.       SEGMENT HYBSTO
  103.          REAL*8 HYBASE(NDIM,NBDDL,NBPP)
  104.       ENDSEGMENT
  105. *
  106.       CHARACTER*8  CNM
  107.       CHARACTER*(NCONCH) CONM
  108.       PARAMETER(NINF=3)
  109.       INTEGER INFOS(NINF)
  110.       LOGICAL lsupma
  111. *
  112. * Recup. des caracteristiques geometriques du maillage elementaire
  113. * et du maillage hybride dual
  114. *
  115.       IMODEL = IPMODE
  116.       SEGACT IMODEL
  117.       CONM   = CONMOD
  118.       IPMAIL = IMAMOD
  119.       MELEME = IPGEOS
  120.       SEGACT MELEME
  121.       NBNN   = NUM(/1)
  122.       NBELEM = NUM(/2)
  123.       NEFHYB = NEFMOD
  124.       NEF    = NUMGEO(NEFHYB)
  125.       MFR    = NUMMFR(NEFHYB)
  126. *
  127.       MRIGID = IPRIGI
  128.       SEGACT MRIGID
  129.       IPT1   = IRIGEL(1,IMAIL)
  130.       SEGACT IPT1
  131.       NBDDL  = IPT1.NUM(/1)
  132.       SEGDES IPT1
  133. *
  134. * Recup. des caracteristiques d'integration de l'EF support geometrique
  135. * de l'EF hybride
  136. *
  137.       CALL TSHAPE(NEF,'GAUSS',IPINTE)
  138.       IF (IERR.NE.0) THEN
  139.          SEGDES IMODEL , MELEME
  140.          RETURN
  141.       ENDIF
  142. *
  143. * Recup. des fonctions de bases hybrides
  144. *
  145.       CALL HYSHPT(NEFHYB,NBDDL,IPINTE,IPTHYB)
  146.       IF (IERR.NE.0) THEN
  147.          SEGDES IMODEL , MELEME
  148.          RETURN
  149.       ENDIF
  150. *
  151. * Activation des segments "d'integration"
  152. *
  153.       MINTE  = IPINTE
  154.       SEGACT MINTE
  155.       NBPGAU = POIGAU(/1)
  156.       HYBSTO = IPTHYB
  157.       SEGACT HYBSTO
  158. *
  159. * Recup. des caracteristiques d'integration au centre de l'EF
  160. *
  161.       CALL RESHPT(1,NBNN,NEF,NEF,0,IPT1,IRT1)
  162.       MINTE1 = IPT1
  163.       SEGACT MINTE1
  164. *
  165. * Initialisation des chapeaux de l'objet rigidité
  166. *
  167.       xMATRI = IRIGEL(4,IMAIL)
  168.       SEGACT xMATRI*MOD
  169.       NLIGRP = NBDDL
  170.       NLIGRD = NBDDL
  171. *
  172. * Recherche du nom du materiau ...trope, de son numero, de sa nature
  173. *
  174.       NFOR = FORMOD(/2)
  175.       NMAT = MATMOD(/2)
  176.       CALL NOMATE(FORMOD,NFOR,MATMOD,NMAT,CNM,INM,INT)
  177.       IF (IERR.NE.0) THEN
  178.          SEGDES IMODEL , MELEME
  179.          SEGDES MINTE  , MINTE1
  180.          SEGSUP xMATRI , MRIGID , HYBSTO
  181.          RETURN
  182.       ENDIF
  183. *
  184. * Remplissage du tableau INFOS (informations sur element).
  185. *
  186.       INFOS(1) = 0
  187.       INFOS(2) = 0
  188.       INFOS(3) = NIFOUR
  189. *
  190. * Remplissage de MOMATR : Noms des composantes obligatoires et
  191. * facultatives contenus dans IPCHEM.
  192. *
  193.       if(lnomid(6).ne.0) then
  194.          nomid=lnomid(6)
  195.          segact nomid
  196.          momatr=nomid
  197.          nmatr=lesobl(/2)
  198.          nmatf=lesfac(/2)
  199.          lsupma=.false.
  200.       else
  201.          lsupma=.true.
  202.          CALL IDMATR(MFR,IMODEL,MOMATR,NMATR,NMATF)
  203.       endif
  204.       NMATT = NMATR
  205.       NV1   = NMATT
  206. *
  207. * Initialisation de NOTYPE : Type des infos contenus dans IPCHEM.
  208. *
  209.       NBTYPE = 1
  210.       SEGINI NOTYPE
  211.       MOTYPE  = NOTYPE
  212.       TYPE(1) ='REAL*8'
  213. *
  214. * Verification des informations transmises et recuperation de IVAMAT,
  215. * pointeur vers le segment MPTVAL contenant les pointeurs vers les
  216. * composantes obligatoires du MCHAML de caracteristiques.
  217. *
  218.       CALL KOMCHA(IPCHEM,IPGEOS,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT)
  219.       SEGSUP NOTYPE
  220.       SEGACT MELEME
  221.       IF (IERR.NE.0) THEN
  222.          SEGDES MELEME
  223.          SEGDES IMODEL
  224.          SEGDES MINTE , MINTE1
  225.          SEGSUP xMATRI , MRIGID , HYBSTO
  226.          RETURN
  227.       ENDIF
  228. *
  229. * Initialisation des tableaux de travail
  230. *
  231.       NDIM = IDIM * (IDIM+1)
  232.       SEGINI MMAT1 , MVELCH , MAXE
  233. *
  234. *-------------------------------------------------------
  235. * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL
  236. *-------------------------------------------------------
  237. *
  238.       DO 30 IEL=1,NBELEM
  239. *
  240. * Initialisations
  241. *
  242.          IFOIS = 0
  243.          IFOI2 = 0
  244.          CALL ZERO(CEL,NBDDL,NBDDL)
  245. *
  246. * Recuperation des coordonnees globales des noeuds de l'element IEL
  247. *
  248.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  249. *
  250. * Calcul des axes locaux dans le cas orthotrope ou anisotrope
  251. *
  252.          IF (INM.EQ.2.OR.INM.EQ.3) THEN
  253.             NBSH = MINTE1.SHPTOT(/2)
  254.             CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR)
  255.             if (nbsh.eq.-1) then
  256.  
  257.              call erreur(525)
  258.              return
  259.             endif
  260.          ENDIF
  261.  
  262. *--------------------------------
  263. * BOUCLE SUR LES POINTS DE GAUSS
  264. *--------------------------------
  265.          DO 20 IGAU=1,NBPGAU
  266. *
  267. * Calcul de la matrice jacobienne de la fonction de passage du
  268. * repere local au repere global, de son determinant ET recup.
  269. * des fonctions de base hybride au point de gauss IGAU.
  270. *
  271.             CALL MHYBR3(IGAU,NBNN,NBDDL,NDIM,IDIM,IDIM,XE,HYBASE,
  272.      S                  SHPTOT,SHY,SHP,ZJAC,DJAC)
  273. *
  274. * Controle du maillage
  275. *
  276.             IF (DJAC.LT.0.D0) IFOIS = IFOIS + 1
  277.             IF (ABS(DJAC).LT.XPETIT) THEN
  278.                IFOI2 = IFOI2 + 1
  279.                DJAC  = XPETIT
  280.             ENDIF
  281. *
  282. * Calcul du poids d'integration global affecte dans DJAC.
  283. *
  284.             DJAC = POIGAU(IGAU) / ABS(DJAC)
  285. *
  286. * Recherche des valeurs des composantes de la matrice de permeabilite :
  287. * VALMAT(im) contient la im° composante au point de gauss IGAU.
  288. *
  289.             MPTVAL = IVAMAT
  290.             DO 10 IM=1,NMATT
  291.                IF (IVAL(IM).NE.0) THEN
  292.                   MELVAL     = IVAL(IM)
  293.                   IBMN       = MIN(IEL,VELCHE(/2))
  294.                   IGMN       = MIN(IGAU,VELCHE(/1))
  295.                   VALMAT(IM) = VELCHE(IGMN,IBMN)
  296.                ELSE
  297.                   VALMAT(IM) = 0.D0
  298.                ENDIF
  299.   10        CONTINUE
  300. *
  301. *= Passage de la matrice de permeabilite du repere local au global
  302. *
  303.             CALL MATGLO(CMAT,CMAT1,CMAT2,TXR,XLOC,XGLOB,VALMAT,
  304.      S                  IDIM,IDIM,INM,IFOMOD)
  305. *
  306. *- Calcul de la contribution du point de gauss IGAU a la matrice
  307. *- elementaire CEL de l'element IEL :
  308. *-    POIGAU/DJAC * transpose( ZJAC*SHY ) *inv(CMAT)* ( ZJAC*SHY )
  309. *- On calcule CMAT2=inv(CMAT)                    avec INVRS puis
  310. *- on calcule CMAT1=transpose(ZJAC)*CMAT2*ZJAC   avec PRODT puis
  311. *- on somme POIGAU/DJAC * transp.(SHY)*CMAT1*SHY avec BDBST.
  312. *
  313.             CALL INVRS(CMAT,IDIM,CMAT2,CJAC)
  314.             IF (CJAC.EQ.0.D0) THEN
  315.                CALL ERREUR(406)
  316.                INTERR(1) = IEL
  317.                CALL ERREUR(259)
  318.             ENDIF
  319.             IF (IERR.NE.0) THEN
  320.                SEGSUP xMATRI , MRIGID
  321.                GOTO 40
  322.             ENDIF
  323.             CALL PRODT(CMAT1,CMAT2,ZJAC,IDIM,IDIM)
  324.             CALL BDBST(SHY,DJAC,CMAT1,NBDDL,IDIM,CEL)
  325.   20     CONTINUE
  326. *
  327. * Le jacobien est negatif --> MAILLAGE INCORRECT
  328. *
  329.          IF (IFOIS.NE.0.AND.IFOIS.NE.NBPGAU) THEN
  330.             INTERR(1) = IEL
  331.             CALL ERREUR(195)
  332.             SEGSUP xMATRI , MRIGID
  333.             GOTO 40
  334.          ENDIF
  335. *
  336. * Le jacobien est tres petit --> MAILLAGE INCORRECT
  337. *
  338.          IF (IFOI2.EQ.NBPGAU) THEN
  339.             INTERR(1) = IEL
  340.             CALL ERREUR(259)
  341.             SEGSUP xMATRI , MRIGID
  342.             GOTO 40
  343.          ENDIF
  344. *
  345. * Lump pour les carrés
  346.       IF(ILUMP.GT.0)THEN
  347. *
  348.          DO 60 II= 1,NBDDL
  349.          SCLI=0.D0
  350.          DO 65 JJ=1,NBDDL
  351.          SCLI= ABS(CEL(II,JJ))+SCLI
  352.          CEL(II,JJ)=0.D0
  353.    65    CONTINUE
  354.          CEL(II,II)=SCLI
  355.    60    CONTINUE
  356.       ENDIF
  357. *
  358. * Inversion de la matrice masse hybride
  359. *
  360.          CALL INVER(CEL,NBDDL,ICRIT,CEL1,ICSTO,XPETIT)
  361.          IF (ICRIT.EQ.1) THEN
  362.             MOTERR(1:8) = 'DARCY'
  363.             CALL ERREUR(700)
  364.             SEGSUP xMATRI , MRIGID
  365.             GOTO 40
  366.          ENDIF
  367. *
  368. * Remplissage de XMATRI
  369. *
  370. *         SEGINI XMATRI
  371. *         IMATTT(IEL) = XMATRI
  372.          CALL REMPMT(CEL,NBDDL,RE(1,1,iel))
  373. *         SEGDES XMATRI
  374.   30  CONTINUE
  375. *
  376. * Desactivation des segments
  377. *
  378.       SEGDES xMATRI , MRIGID
  379.   40  CONTINUE
  380.       SEGSUP MMAT1 , MVELCH , MAXE , HYBSTO
  381.       SEGDES MELEME
  382.       SEGDES IMODEL
  383.       SEGDES MINTE , MINTE1
  384. *
  385.       MPTVAL = IVAMAT
  386.       DO 50 I=1,NMATT
  387.          MELVAL = IVAL(I)
  388.          SEGDES MELVAL
  389.   50  CONTINUE
  390.       SEGSUP MPTVAL
  391.       NOMID  = MOMATR
  392.       if(lsupma)SEGSUP NOMID
  393. *
  394.       RETURN
  395.       END
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  

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