Télécharger prom.eso

Retour à la liste

Numérotation des lignes :

  1. C PROM SOURCE PV 20/03/24 21:20:24 10554
  2.  
  3. SUBROUTINE PROM(IPMODE,IPCHE,IPCHT,ISUP,IPOUT)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. EXTERNAL SHAPE
  7. *-------------------------------------------------------------------- *
  8. * *
  9. * Sous-programme appelé par l'opérateur PROI *
  10. * Projection d'un chamelem sur le support d un modele par *
  11. * minimisation de l integrale (Ue-Us)**2 ds sur chaque element *
  12. * Ue evalue a partir des fonctions de formes du modele associe *
  13. * au champ d origine si il est fourni ( sinon a partir des *
  14. * fonctions de transformation geometriques de l element) *
  15. * *
  16. * L integrale est exprime à partir des fonctions de forme associees*
  17. * au modèle support recepteur *
  18. *---------------------------------------------------------------------*
  19. * *
  20. * ENTREE : *
  21. * IPMODE modele sur lequel on va projeter *
  22. * IPCHE MCHAML de caracteristique (si le maodele comprends *
  23. * des coques) *
  24. * IPCHT MCHAML a projeter (obligatoirement defini aux noeuds)
  25. * ISUP remis a 1 pour le moment (on travaille aux noeuds)*
  26. * *
  27. * SORTIE : *
  28. * IPOUT MCHAML resultat *
  29. * *
  30. *---------------------------------------------------------------------*
  31. *
  32.  
  33. -INC PPARAM
  34. -INC CCOPTIO
  35. -INC SMCOORD
  36. -INC SMMODEL
  37. -INC SMCHAML
  38. -INC SMELEME
  39. -INC SMCHPOI
  40. -INC SMLCHPO
  41. -INC SMLMOTS
  42. -INC TMTRAV
  43. -INC SMINTE
  44. -INC CCREEL
  45. *
  46. PARAMETER ( O1=1.D0,O2=2.D0)
  47. character*4 nomm
  48. segment snomip
  49. character*4 nomip(0)
  50. endsegment
  51. * SEGMENT MDATA
  52. * INTEGER IPCA(*),IPEXI(*)
  53. * ENDSEGMENT
  54. segment ipcb(nsschp)
  55. segment icpr(2,nbpts)
  56. segment icpr1(nbpts)
  57. C
  58. segment info
  59. integer infell(JG)
  60. endsegment
  61. C
  62. segment sicoq
  63. integer icocoq(nbsous),imfr(nbsous),iiele(nbsous)
  64. integer iminte(nbsous),imint1(nbsous),ibnap(nbsous)
  65. integer inbg(nbsous),ibbas(nbsous)
  66. endsegment
  67. C
  68. segment WRK5
  69. REAL*8 W(200),A(NBN1,NBN1),B(NBN1,NBN1),XPU(3),RHS(NBN1)
  70. REAL*8 SHPXXX(6,NBNN,200),XX(3,200),VAL(NBN1),POIDS(200)
  71. INTEGER NDD
  72. endsegment
  73. C
  74. SEGMENT WRK4
  75. REAL*8 BPSS(3,3),XEL(3,NBN1),XU(3),SHP(6,NBN1),S(6)
  76. REAL*8 TXR(3,3,NBN1),XE(3,NBG)
  77. ENDSEGMENT
  78. SEGMENT VECT
  79. REAL*8 VEC1(IDIM)
  80. REAL*8 VEC2(IDIM)
  81. REAL*8 VECN(IDIM,ibbas(isous))
  82. ENDSEGMENT
  83.  
  84. C= Logique a TRUE si mode de calcul est axisymetrique 1D/2D ou Fourier
  85. LOGICAL LOGAXI
  86. C= Quelques constantes (2.Pi et 4.Pi)
  87. PARAMETER (X2Pi =6.283185307179586476925286766559D0)
  88. PARAMETER (X4Pi=12.566370614359172953850573533118D0)
  89.  
  90. C on remet isup a 1 pour le moment
  91. ISUP=1
  92. MMODEL=IPMODE
  93. SEGACT,MMODEL
  94. NSOUS=KMODEL(/1)
  95. C
  96. ndec=0
  97. call lirent(ndec,0,iret)
  98. if(ndec.eq.0) ndec = 3
  99. *
  100. segact mcoord*mod
  101. NBNOE = nbpts
  102.  
  103. if(ipche.NE.0) then
  104. MCHELM=IPCHE
  105. segact MCHELM
  106. NCHAM=ICHAML(/1)
  107. endif
  108. *
  109. * Création du maillage provisoire support pour interpolation
  110. *
  111. NBSOUS = NSOUS
  112. NBREF = 0
  113. NBELEM = 0
  114. NBNN = 0
  115. C ipt2 contiendra les super couche ou points de gauss
  116. SEGINI IPT2
  117. SEGINI SICOQ
  118. c listmots des phases
  119. ilphmo = -1
  120. jgn = 8
  121. jgm = nsous
  122. segini mlmots
  123. ilphmo = mlmots
  124. jgm = 1
  125. *
  126. * Boucle sur l'ensemble des sous zones du modeles sur lequel on va
  127. * projetter pour initialisations
  128. nponou=0
  129. npcoq =0
  130. DO 30 ISOUS = 1, NSOUS
  131. *
  132. IMODEL = KMODEL(ISOUS)
  133. SEGACT,IMODEL
  134. MELEME = IMAMOD
  135. SEGACT MELEME
  136. MELE = NEFMOD
  137. NBELE1= NUM(/2)
  138. NBN1 = NUM(/1)
  139. ibnap(isous)=1
  140. C
  141. C COQUE INTEGREE OU PAS ?
  142. C
  143. IF(INFMOD(/1).NE.0)THEN
  144. NPINT=INFMOD(1)
  145. if( npint.ne.0)IBNAP(isous) =NPINT
  146. ELSE
  147. NPINT=0
  148. C ibnap(isous)=1
  149. ENDIF
  150.  
  151. if( infmod(/1).lt.2+isup) then
  152. call elquoi(mele,npint,isup,ipinf,imodel)
  153. if (ierr.ne.0) then
  154. segsup mchelm
  155. iret=0
  156. return
  157. endif
  158. info=ipinf
  159. mfr =infell(13)
  160. iele =infell(14)
  161. minte=infell(11)
  162. minte1=infell(12)
  163. nbg=infell(6)
  164. segsup info
  165. else
  166. minte=infmod(isup+2)
  167. minte1=infMOD(12)
  168. nbg=infele(6)
  169. mfr =infele(13)
  170. iele =infele(14)
  171. endif
  172. imfr(isous)=mfr
  173. iiele(isous)=iele
  174. iminte(isous)=minte
  175. imint1(isous)=minte1
  176. inbg(isous)=nbg
  177. * write(6,*) 'mfr iele mele minte minte1 nbg',mfr,iele,
  178. * & mele,minte,minte1,nbg
  179. NBSOUS = 0
  180. NBREF = 0
  181. NBELEM = NBELE1
  182.  
  183. IF(mfr.eq.3.or.mfr.eq.5.or.mfr.eq.9) THEN
  184. C------------------coques
  185. IF(IPCHE.eq.0) THEN
  186. C pour le coques il faut les caracteristiques
  187. call erreur(404)
  188. * segsup info
  189. return
  190. ENDIF
  191. C
  192. icocoq(isous)=1
  193. nbnn = ndec*ndec
  194. if(iele.lt.4) nbnn=ndec
  195. ibbas(isous)= nbnn
  196. if(idim.eq.3.or.(iele.lt.4.and.idim.eq.2)) then
  197. ibnap(isous)=ndec
  198. nbnn= nbnn*ndec
  199. inbg(isous)=nbnn
  200. endif
  201.  
  202. ELSE
  203. C--------------- massifs
  204.  
  205. IF (mfr.EQ.1) THEN
  206. IF (IDIM.EQ.3) THEN
  207. NBNN=ndec*ndec*ndec
  208. ELSE IF (IDIM.EQ.2) THEN
  209. NBNN=ndec*ndec
  210. C* ELSE IF (IDIM.EQ.1) THEN
  211. ELSE
  212. NBNN=ndec
  213. ENDIF
  214. ELSE
  215. NBNN=ndec*ndec
  216. ENDIF
  217.  
  218. C petite rectif utile pour les elements speciaux
  219. if(iele.eq.32) then
  220. C les polygones
  221. nbnn =nbnn*nbn1
  222. elseif(iele.eq.25.or.iele.eq.26) then
  223. C les pyramides
  224. nbnn = nbnn*4
  225. endif
  226. C
  227. inbg(isous)= NBNN
  228. ENDIF
  229. *
  230. SEGINI IPT1
  231. IPT1.ITYPEL = 28
  232. ipt2.lisous(isous) = ipt1
  233. nponou=nponou+nbnn*nbelem
  234.  
  235. if (isous.eq.1) then
  236. mots(1) = conmod(17:24)
  237. else
  238. do ipl = 1,jgm
  239. if (mots(ipl).eq.conmod(17:24)) goto 27
  240. enddo
  241. jgm = jgm + 1
  242. mots(jgm) = conmod(17:24)
  243. 27 continue
  244. endif
  245.  
  246. * segsup info
  247. 30 CONTINUE
  248. segadj mlmots
  249. C-----------------------------------------------------------
  250. C
  251. C
  252. C on va fabriquer des points support provisoires a l interieur
  253. C de chaque element des maillages du modele recepteur
  254. C
  255. C
  256. C----------------------------------------------------------
  257. C write(6,*) ' nombre de points crees ' ,nponou
  258. C (DIP) ajustement de MCOORD pour les noeuds provisoires
  259. C il sera remis a sa taille initiale a la fin du programme
  260. NBPTS = NBNOE + nponou
  261. SEGADJ MCOORD
  262. *-----------------------------------------------------------
  263. * Boucle sur l'ensemble des sous zones du modeles
  264. *------------------------------------------------------------
  265. inupo =NBNOE
  266. DO 100 ISOUS = 1,NSOUS
  267. IMODEL = KMODEL(ISOUS)
  268. segact imodel
  269. MELEME = IMAMOD
  270. segact meleme
  271. nbele1= num(/2)
  272. nbn1 = num(/1)
  273. ipt1=ipt2.lisous(isous)
  274. segact ipt1*mod
  275. NNB = ibbas(isous)
  276. NBG = inbg(isous)
  277. SEGINI WRK4
  278. MELVA1 = 0
  279. MELVA2 = 0
  280. IF(icocoq(isous).ne.0) THEN
  281. C -------------------------------la sz est une coque ------
  282.  
  283. NUCHA = 0
  284. DO 15, NUCH = 1, NCHAM
  285. *
  286. IF ( CONCHE(NUCH).EQ.CONMOD.AND.
  287. & IMACHE(NUCH).EQ.IMAMOD) NUCHA = NUCH
  288. *
  289. 15 CONTINUE
  290. *
  291. IF (NUCHA.NE.0) THEN
  292. MCHAML=ICHAML(NUCHA)
  293. SEGACT,MCHAML
  294. *
  295. MELVA1 = 0
  296. MELVA2 = 0
  297. NCOMP = IELVAL(/1)
  298. DO 20, I = 1, NCOMP
  299. IF (NOMCHE(I)(1:4).EQ.'EPAI') THEN
  300. MELVA1 = IELVAL(I)
  301. SEGACT, MELVA1
  302. ELSEIF (NOMCHE(I)(1:4).EQ.'EXCE') THEN
  303. MELVA2 = IELVAL(I)
  304. SEGACT, MELVA2
  305. ENDIF
  306. 20 CONTINUE
  307. ENDIF
  308. ENDIF
  309. C
  310. MELE = NEFMOD
  311. c coque integree ou pas ?
  312. if(infmod(/1).ne.0)then
  313. npint=infmod(1)
  314. else
  315. npint=0
  316. endif
  317. C
  318. nbnap=ibnap(isous)
  319. mfr =imfr(isous)
  320. iele =iiele(isous)
  321. minte=iminte(isous)
  322. minte1=imint1(isous)
  323. nbg = inbg(isous)
  324. C
  325. segact minte
  326. if(mfr.eq.5) segact minte1
  327.  
  328. nnb = ibbas(isous)
  329. C
  330. ipo=0
  331. *-------------------------------------------------------------------
  332. * Pavage de l element de reference
  333. *-------------------------------------------------------------------
  334. NBNN = inbg(isous)
  335. segini WRK5
  336. iwkr= wrk5
  337. iele=iiele(isous)
  338. CALL PAVPOI(WRK5,iele,ndec,nbn1)
  339. segini vect
  340. *-------------------------------------------------------------------
  341. * Boucle sur les elements de la sous zone du modele
  342. *-------------------------------------------------------------------
  343.  
  344. DO IEL=1,NBELE1
  345. CALL DOXE(XCOOR,IDIM,NBN1,NUM,IEL,XEL)
  346.  
  347. C si c est une coque pouvant etre excentree
  348. C recherche des normales aux noeuds
  349. if(icocoq(isous).ne.0) then
  350. C
  351. IF(IELE.EQ.4) THEN
  352. C --------coq3 dkt dst ----------------
  353. do 98 i=1,idim
  354. VEC1(i)=xel(i,2)-xel(i,1)
  355. 98 VEC2(i)=xel(i,3)-xel(i,1)
  356. vecn(1,1) = vec1(2)*vec2(3)-vec1(3)*vec2(2)
  357. vecn(2,1) = vec1(3)*vec2(1)-vec1(1)*vec2(3)
  358. vecn(3,1) = vec1(1)*vec2(2)-vec1(2)*vec2(1)
  359. vnor =sqrt(vecn(1,1)*vecn(1,1)+vecn(2,1)*vecn(2,1)
  360. & +vecn(3,1)*vecn(3,1))
  361. do 91 i=1,idim
  362. 91 vecn(i,1)=vecn(i,1)/vnor
  363. do 99 i=1,idim
  364. do 99 j=2,ibbas(isous)
  365. 99 vecn(i,j)=vecn(i,1)
  366. ELSEIF (IELE.EQ.8) THEN
  367. C --------------- coq4 ----------------
  368. call cq4loc(xel,xe,BPSS,irrt,1)
  369. C ici c est une estimation du plan moyen
  370. * write(6,fmt='(3E12.5)') ((bpss(i,j),i=1,3),j=1,3)
  371. do i=1,3
  372. do ip=1,ibbas(isous)
  373. vecn(i,ip)=bpss(3,i)
  374. enddo
  375. enddo
  376. ELSEIF(IELE.eq.2) THEN
  377. C ---------------- coq2 --------------
  378. vnor=0.D0
  379. do 92 i=1,idim
  380. VEC1(i)=xel(i,2)-xel(i,1)
  381. 92 vnor=vnor + vec1(i)*vec1(i)
  382. vnor = sqrt(vnor)
  383. vecn(1,1)=-vec1(2)/vnor
  384. vecn(2,1)=vec1(1)/vnor
  385. vecn(1,2)=vecn(1,1)
  386. vecn(2,2)=vecn(2,1)
  387. ELSEIF (iele.eq.5.or.iele.eq.10) THEN
  388. C preparation pour coques integrees
  389. SURFX=XZERO
  390. SURFY=XZERO
  391. SURFZ=XZERO
  392. ENDIF
  393. endif
  394. C ------------------- boucle sur les points d integration
  395. inuo=inupo
  396. C write(6,*) ' ndd et inupo ' ,ndd,inupo
  397. do igau=1,ndd
  398. C
  399. do j=1,idim
  400. xpu(j)=xx(j,igau)
  401. enddo
  402. C write(6,*)(xpu(J),j=1,idim)
  403. C
  404. if(mele.GE.111.and.mele.le.122) then
  405. C les plolygones
  406. call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  407. else
  408. CALL SHAPE(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  409. endif
  410.  
  411. C PROJECTION DANS L ELEMENT REEL
  412. DO 501 IK = 1,IDIM
  413. XU(IK)= XZERO
  414. DO 501 JK = 1,nbn1
  415. XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK))
  416. 501 CONTINUE
  417. C
  418. if ( mfr.eq.5) then
  419. C coques integrees ds l epaisseur
  420. DO I=1,6
  421. S(I)=XZERO
  422. CONTINUE
  423. ENDDO
  424. DO INOE=1,NBN1
  425. S(1)=S(1)+SHP(2,INOE)*XEL(2,INOE)
  426. S(2)=S(2)+SHP(3,INOE)*XEL(3,INOE)
  427. S(3)=S(3)+SHP(3,INOE)*XEL(2,INOE)
  428. S(4)=S(4)+SHP(2,INOE)*XEL(3,INOE)
  429. S(5)=S(5)+SHP(3,INOE)*XEL(1,INOE)
  430. S(6)=S(6)+SHP(2,INOE)*XEL(1,INOE)
  431. ENDDO
  432. SURFX=S(1)*S(2)-S(3)*S(4)
  433. SURFY=S(4)*S(5)-S(2)*S(6)
  434. SURFZ=S(6)*S(3)-S(5)*S(1)
  435. SURF=SQRT(SURFX**2+SURFY**2+SURFZ**2)
  436. C write(6,*) ' normales au pt igau du plan de base '
  437. vecn(1,igau)=surfx/surf
  438. vecn(2,igau)=surfy/surf
  439. vecn(3,igau)=surfz/surf
  440. endif
  441. C
  442. IF(ICOCOQ(ISOUS).EQ.0) THEN
  443. C elements standards
  444. inupo=inupo+1
  445. ipt1.num(igau,iel)=inupo
  446. do i=1,idim
  447. xcoor((inupo-1)*(idim+1)+i)=XU(i)
  448. enddo
  449. C write(6,3335) (xu(i),i=1,idim)
  450. C3335 format(3e12.5)
  451. ELSE
  452. C traitement special des coques
  453. if(melva1.ne.0) then
  454. ibmn =MIN(iel,melva1.velche(/2))
  455. igmn =MIN(igau,melva1.velche(/1))
  456. epai = melva1.velche(igmn,ibmn)
  457. else
  458. epai =xzero
  459. endif
  460. C
  461. if(melva2.ne.0) then
  462. ibmn =MIN(iel,melva2.velche(/2))
  463. igmn =MIN(igau,melva2.velche(/1))
  464. exce = melva2.velche(igmn,ibmn)
  465. else
  466. exce =xzero
  467. endif
  468.  
  469. dp = epai/ibnap(isous)
  470. dnor = exce -(epai + dp)/2.D0
  471. C
  472. do inap=1,ibnap(isous)
  473. dnor=dnor+dp
  474. ipos = (inap-1)*ibbas(isous)+igau
  475. ipt1.num(ipos,iel)=inuo+ipos
  476. do i=1,idim
  477. xcoor((inuo+ipos-1)*(idim+1)+i)= XU(i)+dnor*vecn(i,igau)
  478. enddo
  479. enddo
  480. inupo=inupo+ibnap(isous)
  481. ENDIF
  482.  
  483. C ---------------fin de la boucle sur les points de d integration
  484. enddo
  485.  
  486. C fin de la boucle sur les elements
  487. enddo
  488. segsup wrk5,vect
  489. SEGSUP WRK4
  490. * segsup info
  491. *
  492. 100 CONTINUE
  493. C fin de creation des noeuds supports provisoires
  494. 2003 format(3I4,2X,4e12.5,I4)
  495. C--------------//////////////////////////////////------------------
  496. C maintenant in faut un meleme poi1 de tous les points sur lesquels
  497. C on doit interpoler une valeur
  498. segini icpr1
  499. ipt5=ipt2
  500. nbelem=0
  501. do io=1,max(1,ipt2.lisous(/1))
  502. if (ipt2.lisous(/1).ne.0) ipt5=ipt2.lisous(io)
  503. segact ipt5
  504. do ip=1,ipt5.num(/1)
  505. do il=1,ipt5.num(/2)
  506. ipt=ipt5.num(ip,il)
  507. if (icpr1(ipt).eq.0) then
  508. nbelem=nbelem+1
  509. icpr1(ipt)=nbelem
  510. endif
  511. enddo
  512. enddo
  513. enddo
  514. nbnn=1
  515. nbsous=0
  516. nbref=0
  517. ** write (6,*) ' nbelem ',nbelem
  518. segini ipt4
  519. ipgeom=ipt4
  520. nbelem=0
  521. do i=1,icpr1(/1)
  522. if (icpr1(i).ne.0) then
  523. nbelem=nbelem+1
  524. ** write (6,*) ' i nbelem icpr ',i,nbelem,icpr1(i)
  525. ipt4.num(1,nbelem)=i
  526. endif
  527. enddo
  528. segsup icpr1
  529. C----------------------------------------------------------------
  530. C----------------------------------------------------------------
  531. *-------------on est pret a faire l interpolation sur ipgeom ----
  532. C
  533. *write(6,*) ' maillage apres changer poi1 ' ,ipgeom
  534. C call ecmail(ipgeom)
  535. isort=0
  536. CALL PRO2(IPGEOM,IPCHT,isort,IPOUT,ilphmo)
  537. if (ierr.ne.0) return
  538. C----------------------------------------------------------------
  539. C----------------------------------------------------------------
  540. C
  541. C Definition de "constantes" selon le mode de calcul
  542. C Modes de calcul 1D/2D axisymetrique
  543. IF (IFOMOD.EQ.0.OR.IFOMOD.EQ.4) THEN
  544. LOGAXI=.TRUE.
  545. CteAxi=X2Pi
  546. C Mode de calcul 2D Fourier
  547. ELSE IF (IFOMOD.EQ.1) THEN
  548. LOGAXI=.TRUE.
  549. IF (NIFOUR.EQ.0) THEN
  550. CteAxi=X2Pi
  551. ELSE
  552. CteAxi=XPi
  553. ENDIF
  554. ELSE
  555. LOGAXI=.FALSE.
  556. ENDIF
  557. C----------------------------------------------------------------
  558. C --- il faut maintenant reconstituer les MCHAML
  559. C --- a partir du chpo construit sur ipgeom ----------------
  560. mlchpo = ipout
  561. segact mlchpo
  562. mchel5= ipcht
  563. segact mchel5
  564. N1=NSOUS
  565. L1=12
  566. N3=6
  567. SEGINI MCHEL1
  568. * MCHEL1.TITCHE='SCALAIRE'
  569. MCHEL1.TITCHE = mchel5.titche
  570. MCHEL1.IFOCHE=IFOUR
  571. *//////////////////////////////
  572.  
  573. * boucle phases modele
  574. DO 7000 IPHAS = 1,JGM
  575. * write(6,*) 'iphas ',iphas
  576. segini icpr,snomip
  577. * mdata=isort
  578. * segact mdata
  579. mchpoi=ichpoi(iphas)
  580. segact mchpoi
  581. nsschp=ipchp(/1)
  582. segini ipcb
  583.  
  584. C ip=0
  585. do i=1,nsschp
  586. * mchpoi=ipca(i)
  587. C call ecchpo(mchpoi)
  588. * segact mchpoi
  589. msoupo=ipchp(i)
  590. ipcb(i)=msoupo
  591. segact msoupo
  592. mpoval=ipoval
  593. segact mpoval
  594. nc = nocomp(/2)
  595. C
  596. C tableau des general des composantes
  597. do ic = 1,nc
  598. nomm= nocomp(ic)
  599. if(nomip(/2).eq.0) then
  600. nomip(**)=nomm
  601. else
  602. do k=1,nomip(/2)
  603. if(nomm.eq.nomip(k)) goto 4322
  604. enddo
  605. nomip(**)=nomm
  606. 4322 continue
  607. endif
  608. enddo
  609. * write(6,*) (nomip(l),l=1,nomip(/2))
  610. * reperage de la position ses points dans le chpo
  611. ipt5=igeoc
  612. segact ipt5
  613. * ip=0
  614. do j=1,ipt5.num(/2)
  615. * ip=ip+1
  616. icpr(1,ipt5.num(1,j))=j
  617. icpr(2,ipt5.num(1,j))= i
  618. enddo
  619. enddo
  620. C i=icpr(2,k) point k venant du msoupo du msoupo i dans ipcb
  621.  
  622. inupo=0
  623. segact ipt2
  624. DO 2100 isous = 1,nsous
  625. imodel= KMODEL(isous)
  626. segact imodel
  627. if (conmod(17:24).ne.mots(iphas)) goto 2100
  628. meleme = imamod
  629. segact meleme
  630. idi2=idim
  631. ipt5=ipt2.lisous(isous)
  632. segact ipt5
  633. if(icocoq(isous).ne.0) idi2 = idim-1
  634. mfr=imfr(isous)
  635. * création du nouveau segment MCHAML
  636. N2=nomip(/2)
  637. SEGINI MCHAML
  638. MCHEL1.IMACHE(isous)=MELEME
  639. MCHEL1.ICHAML(isous)=MCHAML
  640. MCHEL1.CONCHE(isous)=CONMOD
  641. MCHEL1.INFCHE(isous,1)=0
  642. MCHEL1.INFCHE(isous,2)=0
  643. MCHEL1.INFCHE(isous,3)=NIFOUR
  644. MCHEL1.INFCHE(isous,4)=0
  645. MCHEL1.INFCHE(isous,5)=0
  646. MCHEL1.INFCHE(isous,6)=1
  647. *
  648. N1EL=NUM(/2)
  649. N2PTEL=0
  650. N2EL=0
  651. C--------------------------------------------------------
  652. C boucle sur les composantes
  653. C--------------------------------------------------------
  654. NBNN=num(/1)
  655. N1PTEL= NBNN
  656. NBN1=NBNN
  657.  
  658. C
  659. SEGINI WRK5,wrk4
  660. iele = iiele(isous)
  661. iwrk= wrk5
  662. CALL PAVPOI(WRK5,iele,ndec,nbnn)
  663. C ----- remplissage de shpxxx ( shptot)
  664. DO IGAU=1,NDD
  665. do j=1,idim
  666. xpu(j)=xx(j,igau)
  667. enddo
  668. if(mele.GE.111.and.mele.le.122) then
  669. call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  670. else
  671. call shape(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  672. endif
  673. do id=1,6
  674. do no=1,nbnn
  675. shpxxx(id,no,igau)=shp(id,no)
  676. enddo
  677. enddo
  678. ENDDO
  679. N1PTEL=NBNN
  680. C ----------- creation des melvals recepteurs -----
  681. DO icomp=1,NOMIP(/2)
  682. SEGINI MELVAL
  683. NOMCHE(icomp)=nomip(icomp)
  684. TYPCHE(icomp)='REAL*8'
  685. IELVAL(icomp)=MELVAL
  686. ENDDO
  687. C--------------------------------------------------
  688. C-------- boucle sur les elements ----------
  689. AXIS=O1
  690. DO 2210 iel=1,num(/2)
  691. CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XEL)
  692. CALL ZERO(A,NBNN,NBNN)
  693. C SSOM1=XZERO
  694. C
  695. C
  696. DO inap=1,ibnap(isous)
  697. DO 40 IGAU=1,ndd
  698. do i=1,idim
  699. xpu(i)=xx(i,igau)
  700. enddo
  701. if(mele.GE.111.and.mele.le.122) then
  702. C les plolygones
  703. call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  704. else
  705. call shape(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  706. endif
  707. CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC)
  708. C Modes de calcul 1D/2D axisymetrique et 2D Fourier
  709. IF (LOGAXI) THEN
  710. CALL DISTRR(XEL,SHP,NBNN,RR)
  711. AXIS=CteAxi*RR
  712. C Mode de calcul 1D spherique
  713. ELSE IF (IFOMOD.EQ.5) THEN
  714. CALL DISTRR(XEL,SHP,NBNN,RR)
  715. AXIS=X4Pi*RR*RR
  716. ENDIF
  717. C
  718. C on veut calculer la fonction au point d integration projete
  719. C PROJECTION
  720. DO IK = 1,IDIM
  721. XU(IK)= XZERO
  722. DO JK = 1,NBNN
  723. XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK))
  724. ENDDO
  725. ENDDO
  726. C
  727. CPZA=poids(igau)*DJAC*AXIS
  728. DO 50 INO1=1,NBNN
  729. DO 60 INO2=1,NBNN
  730. Q = SHPXXX(1,INO1,IGAU)*SHPXXX(1,INO2,IGAU)*CPZA
  731. CCC Q = SHPXXX(1,INO1,IGAU)*SHPXXX(1,INO2,IGAU)
  732. CCC 1 *poids(igau)*DJAC*AXIS
  733. A(INO1,INO2)=A(INO1,INO2)+Q
  734. 60 CONTINUE
  735. 50 CONTINUE
  736. 40 CONTINUE
  737. enddo
  738. C
  739. C ---- fin preparation matrice 'masse ' sur element courant
  740. C do i=1,NBNN
  741. C do j=1,NBNN garder eventuellement pour calculer
  742. C B(i,j)=A(i,j) la difference sur les integrales
  743. C enddo
  744. C enddo
  745. C
  746. call invalm(A,NBNN,NBNN,KERRE,1.D-10)
  747. C ------------------ boucle sur les composantes
  748. DO 2220 icomp=1,nomip(/2)
  749. idebco = inupo
  750. do i=1,nbnn
  751. RHS(i) =XZERO
  752. enddo
  753. C --------- boucle sur les points d integration
  754. C
  755. AXIS = O1
  756. DO 42 inap=1,ibnap(isous)
  757. DO 41 IGAU=1,ndd
  758. do i=1,idim
  759. xpu(i)=xx(i,igau)
  760. enddo
  761. CALL SHAPE(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  762. CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC)
  763. C Modes de calcul 1D/2D axisymetrique et 2D Fourier
  764. IF (LOGAXI) THEN
  765. CALL DISTRR(XEL,SHP,NBNN,RR)
  766. AXIS=CteAxi*RR
  767. C Mode de calcul 1D spherique
  768. ELSE IF (IFOMOD.EQ.5) THEN
  769. CALL DISTRR(XEL,SHP,NBNN,RR)
  770. AXIS=X4Pi*RR*RR
  771. ENDIF
  772.  
  773. ipop=(inap-1)*ndd+igau
  774. inupo=ipt5.num(ipop,iel)
  775.  
  776. jh=icpr (1,inupo)
  777. ipa=icpr(2,inupo)
  778. C il faut verifier si vpocha existe pour ce point
  779. msoupo = ipcb(ipa)
  780. mpoval=ipoval
  781. do l=1,nocomp(/2)
  782. if(nocomp(l).eq.nomip(icomp)) then
  783. vvv=vpocha(jh,icomp)
  784. goto 4555
  785. endif
  786. enddo
  787. vvv=XZERO
  788. 4555 continue
  789. C write(6,7777) inupo,jh,poids(igau),vvv
  790. C
  791. melval=ielval(icomp)
  792. CPZA=vvv*poids(igau)*DJAC*AXIS
  793. DO INO1 = 1,NBNN
  794. RHS(ino1)=RHS(ino1)+CPZA*SHPXXX(1,ino1,igau)
  795. CCC RHS(ino1)=RHS(ino1)+vvv*poids(igau)*
  796. CCC 1 DJAC*SHPXXX(1,INO1,IGAU)*AXIS
  797. ENDDO
  798. 41 continue
  799. 42 continue
  800.  
  801. 7777 format(2I4,2X,2E12.5)
  802. C SSOM1 =SSOM1+SURF*poids(igau)*vvv
  803. C
  804. do i=1,NBNN
  805. val(i)=XZERO
  806. do j=1,NBNN
  807. val(i)=val(i)+A(i,j)*RHS(J)
  808. enddo
  809. enddo
  810. C
  811. do I=1,NBNN
  812. velche(i,iel)=val(i)
  813. enddo
  814. C
  815. * if(icomp.lt.nomip(/2)) inupo =idebco
  816. C ---- fin de la boucle sur les composantes
  817. 2220 continue
  818. C fin de la boucle sur les elements
  819. 2210 continue
  820. C
  821. do i=1,ielval(/1)
  822. ipc1 = ielval(i)
  823. call comred(ipc1)
  824. ielval(i)=ipc1
  825. melval = ipc1
  826. enddo
  827.  
  828. C---------------- fin du traitement des massifs
  829. segsup wrk5,wrk4,ipt5
  830.  
  831. C segsup ipt3
  832. C fin de la boucle sur les sous zones du modele
  833. 2100 continue
  834. C destruction des chpo intermediaires
  835. do i=1,ipcb(/1)
  836. msoupo=ipcb(i)
  837. cgf ipt5= igeoc (correction 7284)
  838. mpoval= ipoval
  839. * mchpoi=ipca(i)
  840. segsup mpoval,msoupo
  841. cgf segsup ipt5 (correction 7284)
  842. enddo
  843. segsup mchpoi
  844. segsup ipcb,snomip
  845. segsup icpr
  846. 7000 CONTINUE
  847. C (fdp) re-ajustement de MCOORD a sa taille initiale
  848. NBPTS = NBNOE
  849. SEGADJ MCOORD
  850. C retrait des maillages temporaires du pre-conditionnement
  851. c (leurs numero de noeuds depasse la taille de MCOOR)
  852. call crech1b
  853. C (fdp) suppression du maillage temporaire IPT1
  854. ipt1=ipgeom
  855. segsup ipt1
  856. segact ipt2
  857. do io=1,ipt2.lisous(/1)
  858. ipt5=ipt2.lisous(io)
  859. * segsup ipt5
  860. enddo
  861. segsup ipt2
  862. segsup ipt4
  863. if(ipche.ne.0) then
  864. mchelm=ipche
  865. endif
  866. IPOUT=MCHEL1
  867. segsup sicoq
  868.  
  869. return
  870. END
  871.  
  872.  
  873.  
  874.  

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