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

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