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

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