Télécharger mhybr2.eso

Retour à la liste

Numérotation des lignes :

  1. C MHYBR2 SOURCE BP208322 15/06/22 21:20:52 8543
  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. -INC CCOPTIO
  67. -INC CCREEL
  68. -INC SMCOORD
  69. -INC SMINTE
  70. -INC SMMODEL
  71. -INC SMRIGID
  72. -INC SMELEME
  73. -INC SMCHAML
  74. *
  75. SEGMENT MAXE
  76. REAL*8 TXR(IDIM,IDIM),XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM)
  77. ENDSEGMENT
  78. *
  79. SEGMENT MMAT1
  80. REAL*8 CEL(NBDDL,NBDDL),CEL1(NBDDL,NBDDL),XE(3,NBNN)
  81. REAL*8 SHP(6,NBNN),SHY(IDIM,NBDDL),ZJAC(IDIM,IDIM)
  82. 4d)v style="font: normal normal 1em/1.2em monospace; margin:0; padding:0; background:none; vertical-align:top;"> REAL*8 CMAT(IDIM,IDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM)
  • INTEGER ICSTO(NBDDL)
  • ENDSEGMENT
  • *
  • SEGMENT NOTYPE
  • CHARACTER*16 TYPE(NBTYPE)
  • ENDSEGMENT
  • *
  • SEGMENT MVELCH
  • REAL*8 VALMAT(NV1)
  • ENDSEGMENT
  • *
  • SEGMENT MPTVAL
  • INTEGER IPOS(NS) , NSOF(NS)
  • INTEGER IVAL(NCOSOU)
  • CHARACTER*16 TYVAL(NCOSOU)
  • ENDSEGMENT
  • *
  • SEGMENT HYBSTO
  • REAL*8 HYBASE(NDIM,NBDDL,NBPP)
  • ENDSEGMENT
  • *
  • CHARACTER*8 CNM
  • CHARACTER*(NCONCH) CONM
  • PARAMETER(NINF=3)
  • INTEGER INFOS(NINF)
  • LOGICAL lsupma
  • *
  • * Recup. des caracteristiques geometriques du maillage elementaire
  • * et du maillage hybride dual
  • *
  • IMODEL = IPMODE
  • SEGACT IMODEL
  • CONM = CONMOD
  • IPMAIL = IMAMOD
  • MELEME = IPGEOS
  • SEGACT MELEME
  • NBNN = NUM(/1)
  • NBELEM = NUM(/2)
  • NEFHYB = NEFMOD
  • NEF = NUMGEO(NEFHYB)
  • MFR = NUMMFR(NEFHYB)
  • *
  • MRIGID = IPRIGI
  • SEGACT MRIGID
  • IPT1 = IRIGEL(1,IMAIL)
  • SEGACT IPT1
  • NBDDL = IPT1.NUM(/1)
  • SEGDES IPT1
  • *
  • * Recup. des caracteristiques d'integration de l'EF support geometrique
  • * de l'EF hybride
  • *
  • CALL TSHAPE(NEF,'GAUSS',IPINTE)
  • IF (IERR.NE.0) THEN
  • SEGDES IMODEL , MELEME
  • RETURN
  • ENDIF
  • *
  • * Recup. des fonctions de bases hybrides
  • *
  • CALL HYSHPT(NEFHYB,NBDDL,IPINTE,IPTHYB)
  • IF (IERR.NE.0) THEN
  • SEGDES IMODEL , MELEME
  • RETURN
  • ENDIF
  • *
  • * Activation des segments "d'integration"
  • *
  • MINTE = IPINTE
  • SEGACT MINTE
  • NBPGAU = POIGAU(/1)
  • HYBSTO = IPTHYB
  • SEGACT HYBSTO
  • *
  • * Recup. des caracteristiques d'integration au centre de l'EF
  • *
  • CALL RESHPT(1,NBNN,NEF,NEF,0,IPT1,IRT1)
  • MINTE1 = IPT1
  • SEGACT MINTE1
  • *
  • * Initialisation des chapeaux de l'objet rigidité
  • *
  • xMATRI = IRIGEL(4,IMAIL)
  • SEGACT xMATRI*MOD
  • NLIGRP = NBDDL
  • NLIGRD = NBDDL
  • *
  • * Recherche du nom du materiau ...trope, de son numero, de sa nature
  • *
  • NFOR = FORMOD(/2)
  • NMAT = MATMOD(/2)
  • CALL NOMATE(FORMOD,NFOR,MATMOD,NMAT,CNM,INM,INT)
  • IF (IERR.NE.0) THEN
  • SEGDES IMODEL , MELEME
  • SEGDES MINTE , MINTE1
  • SEGSUP xMATRI , MRIGID , HYBSTO
  • RETURN
  • ENDIF
  • *
  • * Remplissage du tableau INFOS (informations sur element).
  • *
  • INFOS(1) = 0
  • INFOS(2) = 0
  • INFOS(3) = NIFOUR
  • *
  • * Remplissage de MOMATR : Noms des composantes obligatoires et
  • * facultatives contenus dans IPCHEM.
  • *
  • if(lnomid(6).ne.0) then
  • nomid=lnomid(6)
  • segact nomid
  • momatr=nomid
  • nmatr=lesobl(/2)
  • nmatf=lesfac(/2)
  • lsupma=.false.
  • else
  • lsupma=.true.
  • CALL IDMATR(MFR,IMODEL,MOMATR,NMATR,NMATF)
  • endif
  • NMATT = NMATR
  • NV1 = NMATT
  • *
  • * Initialisation de NOTYPE : Type des infos contenus dans IPCHEM.
  • *
  • NBTYPE = 1
  • SEGINI NOTYPE
  • MOTYPE = NOTYPE
  • TYPE(1) ='REAL*8'
  • *
  • * Verification des informations transmises et recuperation de IVAMAT,
  • * pointeur vers le segment MPTVAL contenant les pointeurs vers les
  • * composantes obligatoires du MCHAML de caracteristiques.
  • *
  • CALL KOMCHA(IPCHEM,IPGEOS,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT)
  • SEGSUP NOTYPE
  • SEGACT MELEME
  • IF (IERR.NE.0) THEN
  • SEGDES MELEME
  • SEGDES IMODEL
  • SEGDES MINTE , MINTE1
  • SEGSUP xMATRI , MRIGID , HYBSTO
  • RETURN
  • ENDIF
  • *
  • * Initialisation des tableaux de travail
  • *
  • NDIM = IDIM * (IDIM+1)
  • SEGINI MMAT1 , MVELCH , MAXE
  • *
  • *-------------------------------------------------------
  • * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL
  • *-------------------------------------------------------
  • *
  • DO 30 IEL=1,NBELEM
  • *
  • * Initialisations
  • *
  • IFOIS = 0
  • IFOI2 = 0
  • CALL ZERO(CEL,NBDDL,NBDDL)
  • *
  • * Recuperation des coordonnees globales des noeuds de l'element IEL
  • *
  • CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  • *
  • * Calcul des axes locaux dans le cas orthotrope ou anisotrope
  • *
  • IF (INM.EQ.2.OR.INM.EQ.3) THEN
  • NBSH = MINTE1.SHPTOT(/2)
  • CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR)
  • if (nbsh.eq.-1) then
  •  
  • call erreur(525)
  • return
  • endif
  • ENDIF
  •  
  • *--------------------------------
  • * BOUCLE SUR LES POINTS DE GAUSS
  • *--------------------------------
  • DO 20 IGAU=1,NBPGAU
  • *
  • * Calcul de la matrice jacobienne de la fonction de passage du
  • * repere local au repere global, de son determinant ET recup.
  • * des fonctions de base hybride au point de gauss IGAU.
  • *
  • CALL MHYBR3(IGAU,NBNN,NBDDL,NDIM,IDIM,IDIM,XE,HYBASE,
  • S SHPTOT,SHY,SHP,ZJAC,DJAC)
  • *
  • * Controle du maillage
  • *
  • IF (DJAC.LT.0.D0) IFOIS = IFOIS + 1
  • IF (ABS(DJAC).LT.XPETIT) THEN
  • IFOI2 = IFOI2 + 1
  • DJAC = XPETIT
  • ENDIF
  • *
  • * Calcul du poids d'integration global affecte dans DJAC.
  • *
  • DJAC = POIGAU(IGAU) / ABS(DJAC)
  • *
  • * Recherche des valeurs des composantes de la matrice de permeabilite :
  • * VALMAT(im) contient la im° composante au point de gauss IGAU.
  • *
  • MPTVAL = IVAMAT
  • DO 10 IM=1,NMATT
  • IF (IVAL(IM).NE.0) THEN
  • MELVAL = IVAL(IM)
  • IBMN = MIN(IEL,VELCHE(/2))
  • IGMN = MIN(IGAU,VELCHE(/1))
  • VALMAT(IM) = VELCHE(IGMN,IBMN)
  • ELSE
  • VALMAT(IM) = 0.D0
  • ENDIF
  • 10 CONTINUE
  • *
  • *= Passage de la matrice de permeabilite du repere local au global
  • *
  • CALL MATGLO(CMAT,CMAT1,CMAT2,TXR,XLOC,XGLOB,VALMAT,
  • S IDIM,IDIM,INM,IFOMOD)
  • *
  • *- Calcul de la contribution du point de gauss IGAU a la matrice
  • *- elementaire CEL de l'element IEL :
  • *- POIGAU/DJAC * transpose( ZJAC*SHY ) *inv(CMAT)* ( ZJAC*SHY )
  • *- On calcule CMAT2=inv(CMAT) avec INVRS puis
  • *- on calcule CMAT1=transpose(ZJAC)*CMAT2*ZJAC avec PRODT puis
  • *- on somme POIGAU/DJAC * transp.(SHY)*CMAT1*SHY avec BDBST.
  • *
  • CALL INVRS(CMAT,IDIM,CMAT2,CJAC)
  • IF (CJAC.EQ.0.D0) THEN
  • CALL ERREUR(406)
  • INTERR(1) = IEL
  • CALL ERREUR(259)
  • ENDIF
  • IF (IERR.NE.0) THEN
  • SEGSUP xMATRI , MRIGID
  • GOTO 40
  • ENDIF
  • CALL PRODT(CMAT1,CMAT2,ZJAC,IDIM,IDIM)
  • CALL BDBST(SHY,DJAC,CMAT1,NBDDL,IDIM,CEL)
  • 20 CONTINUE
  • *
  • * Le jacobien est negatif --> MAILLAGE INCORRECT
  • *
  • IF (IFOIS.NE.0.AND.IFOIS.NE.NBPGAU) THEN
  • INTERR(1) = IEL
  • CALL ERREUR(195)
  • SEGSUP xMATRI , MRIGID
  • GOTO 40
  • ENDIF
  • *
  • * Le jacobien est tres petit --> MAILLAGE INCORRECT
  • *
  • IF (IFOI2.EQ.NBPGAU) THEN
  • INTERR(1) = IEL
  • CALL ERREUR(259)
  • SEGSUP xMATRI , MRIGID
  • GOTO 40
  • ENDIF
  • *
  • * Lump pour les carrés
  • IF(ILUMP.GT.0)THEN
  • *
  • DO 60 II= 1,NBDDL
  • SCLI=0.D0
  • DO 65 JJ=1,NBDDL
  • SCLI= ABS(CEL(II,JJ))+SCLI
  • CEL(II,JJ)=0.D0
  • 65 CONTINUE
  • CEL(II,II)=SCLI
  • 60 CONTINUE
  • ENDIF
  • *
  • * Inversion de la matrice masse hybride
  • *
  • CALL INVER(CEL,NBDDL,ICRIT,CEL1,ICSTO,XPETIT)
  • IF (ICRIT.EQ.1) THEN
  • MOTERR(1:8) = 'DARCY'
  • CALL ERREUR(700)
  • SEGSUP xMATRI , MRIGID
  • GOTO 40
  • ENDIF
  • *
  • * Remplissage de XMATRI
  • *
  • * SEGINI XMATRI
  • * IMATTT(IEL) = XMATRI
  • CALL REMPMT(CEL,NBDDL,RE(1,1,iel))
  • * SEGDES XMATRI
  • 30 CONTINUE
  • *
  • * Desactivation des segments
  • *
  • SEGDES xMATRI , MRIGID
  • 40 CONTINUE
  • SEGSUP MMAT1 , MVELCH , MAXE , HYBSTO
  • SEGDES MELEME
  • SEGDES IMODEL
  • SEGDES MINTE , MINTE1
  • *
  • MPTVAL = IVAMAT
  • DO 50 I=1,NMATT
  • MELVAL = IVAL(I)
  • SEGDES MELVAL
  • 50 CONTINUE
  • SEGSUP MPTVAL
  • NOMID = MOMATR
  • if(lsupma)SEGSUP NOMID
  • *
  • RETURN
  • END
  •  
  •  
  •  
  •  
  •  
  •  
  • © Cast3M 2003 - Tous droits réservés.
    Mentions légales