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.  
  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.  

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