1 : C MENISM SOURCE CHAT 11/03/16 21:27:52 6902
2 : SUBROUTINE MENISM
3 : *********************************************************************
4 : *
5 : * T_CHPO1 T_CHML1=MESM MOD1 BLOQ1 (('TOUT') or 'ROTA');
6 : *
7 : * CALCUL DES MECHANISMES ELEMENTAIRES D'UNE STRUCTURE
8 : * COMPOSEE DE POUTRE.
9 : *
10 : * P.PEGON (ISPRA) AOUT 1996
11 : *********************************************************************
12 : IMPLICIT INTEGER(I-N)
13 : IMPLICIT REAL*8 (A-H,O-Z)
14 : -INC CCOPTIO
15 : -INC SMRIGID
16 : -INC SMMODEL
17 : -INC SMCOORD
18 : -INC SMELEME
19 : -INC SMCHAML
20 : -INC SMCHPOI
21 : -INC SMTABLE
22 : SEGMENT ICPR(XCOOR(/1)/(IDIM+1))
23 : SEGMENT WORK0
24 : REAL*8 C(MRIGI,2*NCOMP)
25 : REAL*8 Q(2*NCOMP,2*NCOMP)
26 : REAL*8 Q1(2*NCOMP,2*NCOMP+MECAN)
27 : REAL*8 C1(MRIGI,2*NCOMP)
28 : REAL*8 C2(MRIGI,2*NCOMP+MECAN)
29 : REAL*8 XE(3,2)
30 : REAL*8 V1(2*NCOMP)
31 : REAL*8 TH(MRIGI)
32 : ENDSEGMENT
33 : SEGMENT WORK1
34 : REAL*8 A(NEQUA,NINCO)
35 : ENDSEGMENT
36 : SEGMENT WORK2
37 : REAL*8 B(NINCO,NINCO)
38 : ENDSEGMENT
39 : SEGMENT WORK3
40 : REAL*8 B1(NNNCO,NMECA)
41 : ENDSEGMENT
42 : *
43 : PARAMETER(NCOMP=3,NECAN=3,MRIGI=3)
44 : CHARACTER*4 COMP(NCOMP),CCOMP,CANI(NECAN)
45 : PARAMETER (NMOCLE=2)
46 : CHARACTER*4 MOTCLE(NMOCLE)
47 : *
48 : DATA COMP/'UX ','UY ','RZ '/
49 : DATA CANI/'RZP1','RZP2','UP2 '/
50 : DATA MOTCLE/'TOUT','ROTA'/
51 : *
52 : * Lecture des donnees
53 : *
53 : CALL LIROBJ('MMODEL ',MMODEL,1,IRET)
54 : IF(IRET.EQ.0)RETURN
54 : CALL LIROBJ('RIGIDITE',MRIGID,1,IRET)
55 : IF(IRET.EQ.0)RETURN
55 : CALL LIRMOT(MOTCLE,NMOCLE,IVAL,0)
56 : IF(IVAL.EQ.0.OR.IVAL.EQ.1) THEN
57 : MECAN=3
58 : ELSE IF(IVAL.EQ.2)THEN
59 : MECAN=2
60 : ENDIF
61 : *
62 : * On verifie que l'on est en 2D et que les elements sont de
63 : * poutre ou de timo
64 : *
65 : IF (IDIM.NE.2)THEN
66 : WRITE(IOIMP,*)'MENISM: ne fonctionne que en 2D'
66 : CALL ERREUR(19)
67 : RETURN
68 : ENDIF
69 : *
70 : IRET=0
71 : SEGACT,MMODEL
72 : NZONE=KMODEL(/1)
73 : DO IZONE=1,NZONE
74 : IMODEL=KMODEL(IZONE)
75 : SEGACT,IMODEL
76 : IF(NEFMOD.NE.29.AND.NEFMOD.NE.84)THEN
77 : IRET=IRET+1
78 : ENDIF
79 : ENDDO
80 : IF(IRET.GT.0)THEN
81 : WRITE(IOIMP,*)
82 : > 'MENISM: ne fonctionne que pour des POUT ou des TIMO'
82 : CALL ERREUR(16)
83 : GOTO 9999
84 : ENDIF
85 : *
86 : * On numerote les points support actifs
87 : *
88 : SEGINI,ICPR
89 : DO I=1,XCOOR(/1)/(IDIM+1)
90 : ICPR(I)=0
91 : ENDDO
92 : NPOI=0
93 : NBEL=0
94 : DO IZONE=1,NZONE
95 : IMODEL=KMODEL(IZONE)
96 : MELEME=IMAMOD
97 : SEGACT,MELEME
98 : NBELEM=NUM(/2)
99 : NBEL=NBEL+NBELEM
100 : DO I=1,NUM(/1)
101 : DO J=1,NBELEM
102 : IKI=NUM(I,J)
103 : IF (ICPR(IKI).EQ.0)THEN
104 : NPOI=NPOI+1
105 : ICPR(IKI)=NPOI
106 : ENDIF
107 : ENDDO
108 : ENDDO
109 : ENDDO
110 : *
111 : * Verification que les c.l. s'appliquent a des points
112 : * et a des d.o.f.s du probleme
113 : *
114 : NLIM=0
115 : SEGACT,MRIGID
116 : NRI=IRIGEL(/1)
117 : NR=IRIGEL(/2)
118 : DO I=1,NR
119 : MELEME=IRIGEL(1,I)
120 : DESCR =IRIGEL(3,I)
121 : xMATRI=IRIGEL(4,I)
122 : NEGALI=IRIGEL(6,I)
123 : IF (NEGALI.NE.0)THEN
124 : WRITE(IOIMP,*)
125 : > 'MENISM: on ne veut que des rigidites d"egalites'
125 : CALL ERREUR(21)
126 : GOTO 9997
127 : ENDIF
128 : SEGACT,MELEME
129 : IF (ITYPEL.NE.22)THEN
130 : WRITE(IOIMP,*)
131 : > 'MENISM: on ne veut que des rigidites de blocage'
131 : CALL ERREUR(21)
132 : GOTO 9997
133 : ENDIF
134 : NBELEM=NUM(/2)
135 : NBNN=NUM(/1)
136 : NLIM=NLIM+NBELEM
137 : DO J=1,NBELEM
138 : DO K=2,NBNN
139 : IF(ICPR(NUM(K,J)).EQ.0)THEN
140 : WRITE(IOIMP,*)'MENISM: ',
141 : > 'un point soumis a blocage n"est pas dans le modele'
141 : CALL ERREUR(21)
142 : GOTO 9997
143 : ENDIF
144 : ENDDO
145 : ENDDO
146 : SEGACT,xMATRI
147 : NELRIG=re(/3)
148 : IF(NBELEM.NE.NELRIG)THEN
149 : WRITE(IOIMP,*)
150 : > 'MENISM: il y a un probleme sioux dans une des rigidites'
150 : CALL ERREUR(21)
151 : GOTO 9997
152 : ENDIF
153 : SEGACT,DESCR
154 : NLIGRP=LISINC(/2)
155 : IF(NBNN.NE.NLIGRP)THEN
156 : WRITE(IOIMP,*)
157 : > 'MENISM: il y a un probleme sioux dans une des rigidites'
157 : CALL ERREUR(21)
158 : GOTO 9997
159 : ENDIF
160 : DO K=2,NLIGRP
161 : CCOMP=LISINC(K)
162 : DO 1 J=1,NCOMP
163 : IF(COMP(J).EQ.CCOMP)GOTO 2
164 : 1 CONTINUE
165 : WRITE(IOIMP,*)
166 : > 'MENISM: un des d.o.f. de blocage n"est pas dans le modele'
166 : CALL ERREUR(21)
167 : GOTO 9997
168 : 2 CONTINUE
169 : ENDDO
170 : ENDDO
171 : *
172 : * On est en mesure de commencer le travail en se souvenant que
173 : *
174 : * NCOMP=nb de d.o.f par point
175 : * MECAN=nb de mecanisme par element
176 : * MRIGI=nb de contrainte rigide
177 : * NPOI=nb de point support de NCOMP dof
178 : * NBEL=nb d'elements
179 : * NLIM=nb de condition aux limites
180 : *
181 : * En consequence on peut allouer le pb avec
182 : *
183 : * NINCO=nb d'inconnues
184 : * NEQUA=nb d'equations
185 : *
186 : NINCO=NCOMP*NPOI+MECAN*NBEL
187 : NEQUA=MRIGI*NBEL+NLIM
188 : NNNCO=NCOMP*NPOI
189 : NMECA=NINCO-NEQUA
190 : *
191 : * On commence par la matrice A
192 : *
193 : SEGINI,WORK1
193 : CALL ZERO(A,NEQUA,NINCO)
194 : *
195 : * Boucle sur les elements
196 : *
197 : SEGINI,WORK0
198 : IEL=0
199 : DO IZONE=1,NZONE
200 : IMODEL=KMODEL(IZONE)
201 : MELEME=IMAMOD
202 : DO I=1,NUM(/2)
203 : IEL=IEL+1
204 : *
205 : * Calcul de C2 elementaire
206 : *
206 : CALL DOXE(XCOOR,IDIM,2,NUM,I,XE)
206 : CALL TIPOCQ(XE, C,Q,MRIGI,2*NCOMP)
206 : CALL TIPOQ1(Q,2*NCOMP, Q1,2*NCOMP+MECAN)
206 : CALL MATMAT(C,Q1,MRIGI,2*NCOMP,2*NCOMP+MECAN,C2)
207 : *
208 : * Assemblage de C2 dans A
209 : *
210 : * 1) noeud 1 de l'element
211 : *
212 : ILIGN=MRIGI*(IEL-1)
213 : ICOLO=NCOMP*(ICPR(NUM(1,I))-1)
214 : DO J=1,MRIGI
215 : DO K=1,NCOMP
216 : A(ILIGN+J,ICOLO+K)=C2(J,K)
217 : ENDDO
218 : ENDDO
219 : *
220 : * 2) noeud 2 de l'element
221 : *
222 : ICOLO=NCOMP*(ICPR(NUM(2,I))-1)
223 : DO J=1,MRIGI
224 : DO K=1,NCOMP
225 : A(ILIGN+J,ICOLO+K)=C2(J,K+NCOMP)
226 : ENDDO
227 : ENDDO
228 : *
229 : * 3) degres suplementaires de l'element
230 : *
231 : ICOLO=NCOMP*NPOI+MECAN*(IEL-1)
232 : DO J=1,MRIGI
233 : DO K=1,MECAN
234 : A(ILIGN+J,ICOLO+K)=C2(J,K+2*NCOMP)
235 : ENDDO
236 : ENDDO
237 : ENDDO
238 : ENDDO
239 : *
240 : * Assemblage des conditions aux limites
241 : *
242 : ILIGN=MRIGI*NBEL
243 : NRI=IRIGEL(/1)
244 : NR=IRIGEL(/2)
245 : DO I=1,NR
246 : MELEME=IRIGEL(1,I)
247 : DESCR=IRIGEL(3,I)
248 : xMATRI=IRIGEL(4,I)
249 : DO J=1,NUM(/2)
250 : ILIGN=ILIGN+1
251 : * XMATRI=IMATTT(J)
252 : * SEGACT,XMATRI
253 : DO K=2,NUM(/1)
254 : CCOMP=LISINC(K)
255 : DO 3 L=1,NCOMP
256 : IF(COMP(L).EQ.CCOMP)GOTO 4
257 : 3 CONTINUE
258 : 4 CONTINUE
259 : ICOLO=NCOMP*(ICPR(NUM(K,J))-1)
260 : A(ILIGN,ICOLO+L)=RE(1,K,j)
261 : ENDDO
262 : * SEGDES,XMATRI
263 : ENDDO
264 : ENDDO
265 : *
266 : * On alloue B et on appelle le solveur
267 : *
268 : SEGINI,WORK2
268 : CALL GAJOME(A,NEQUA,NINCO,B,IRET,IOIMP)
269 : SEGSUP,WORK1
270 : IF (IRET.NE.0)THEN
271 : SEGSUP,WORK2
272 : WRITE(IOIMP,*)
273 : > 'MENISM: il y a un mecanisme global'
273 : CALL ERREUR(26)
274 : GOTO 9997
275 : ENDIF
276 : *
277 : * On ne garde que la partie active de B dans B1
278 : *
279 : SEGINI,WORK3
280 : DO J=1,NMECA
281 : DO I=1,NNNCO
282 : B1(I,J)=B(I,J+NEQUA)
283 : ENDDO
284 : ENDDO
285 : SEGSUP,WORK2
286 : *
287 : * On prepare les tables de sortie
288 : *
289 : M=NMECA+1
290 : SEGINI,MTABLE
291 : MLOTAB=M
292 : MTABTI(1)='MOT'
292 : CALL POSCHA('SOUSTYPE',IRET)
293 : MTABII(1)=IRET
294 : MTABTV(1)='MOT'
294 : CALL POSCHA('MECANISMES_PAR_NOEUDS',IRET)
295 : MTABIV(1)=IRET
296 : DO I=1,NMECA
297 : MTABTI(I+1)='ENTIER'
298 : MTABII(I+1)=I
299 : MTABTV(I+1)='CHPOINT '
300 : MTABIV(I+1)=0
301 : ENDDO
302 : MTAB1=MTABLE
303 : *
304 : SEGINI,MTABLE=MTAB1
304 : CALL POSCHA('MECANISMES_PAR_ELEMENTS',IRET)
305 : MTABIV(1)=IRET
306 : DO I=1,NMECA
307 : MTABTV(I+1)='MCHAML '
308 : ENDDO
309 : MTAB2=MTABLE
310 : *
311 : * Maillage support des chanpoint solution
312 : *
313 : NBSOUS=0
314 : NBREF=0
315 : NBELEM=NPOI
316 : NBNN=1
317 : SEGINI IPT1
318 : IPT1.ITYPEL=1
319 : DO I=1,XCOOR(/1)/(IDIM+1)
320 : IF (ICPR(I).NE.0)THEN
321 : IPT1.NUM(1,ICPR(I))=I
322 : ENDIF
323 : ENDDO
324 : SEGDES,IPT1
325 : *
326 : * Generation des chanpoint solution
327 : *
328 : DO I=1,NMECA
329 : IF(I.EQ.1)THEN
330 : * prototype MCHPOI
331 : NAT=2
332 : NSOUPO=1
333 : SEGINI,MCHPOI
334 : MOCHDE=' CHPOINT CREE PAR MENISM'
335 : MTYPOI='MECANISM'
336 : JATTRI(1)=1
337 : IFOPOI=IFOUR
338 : MCHPO1=MCHPOI
339 : * prototype MSOUPO
340 : NC=NCOMP
341 : SEGINI,MSOUPO
342 : DO J=1,NCOMP
343 : NOCOMP(J)=COMP(J)
344 : ENDDO
345 : IGEOC=IPT1
346 : MSOUP1=MSOUPO
347 : ELSE
348 : SEGINI,MCHPOI=MCHPO1
349 : SEGINI,MSOUPO=MSOUP1
350 : ENDIF
351 : IPCHP(1)=MSOUPO
352 : SEGDES,MCHPOI
353 : N=NPOI
354 : SEGINI,MPOVAL
355 : IPOVAL=MPOVAL
356 : SEGDES,MSOUPO
357 : DO J=1,NPOI
358 : DO K=1,NCOMP
359 : VPOCHA(J,K)=B1(NCOMP*(J-1)+K,I)
360 : ENDDO
361 : ENDDO
362 : SEGDES,MPOVAL
363 : MTAB1.MTABIV(I+1)=MCHPOI
364 : ENDDO
365 : SEGDES,MTAB1
366 : *
367 : * Generation des mchaml solution
368 : *
369 : DO I=1,NMECA
370 : *
371 : * Chapeau
372 : *
373 : IF(I.EQ.1)THEN
374 : * prototype MCHAML
375 : N2=MECAN
376 : SEGINI,MCHAML
377 : DO J=1,MECAN
378 : NOMCHE(J)=CANI(J)//' '
379 : TYPCHE(J)='REAL*8 '
380 : ENDDO
381 : MCHAM1=MCHAML
382 : * prototype MCHELM
383 : N1=NZONE
384 : L1=9
385 : N3=6
386 : SEGINI,MCHELM
387 : TITCHE='MECANISME'
388 : IFOCHE=IFOUR
389 : DO IZONE=1,NZONE
390 : IMODEL=KMODEL(IZONE)
391 : IMACHE(IZONE)=IMAMOD
392 : INFCHE(IZONE,1)=1
393 : INFCHE(IZONE,2)=0
394 : INFCHE(IZONE,3)=NIFOUR
395 : INFCHE(IZONE,4)=0
396 : INFCHE(IZONE,5)=0
397 : INFCHE(IZONE,6)=0
398 : ENDDO
399 : MCHEL1=MCHELM
400 : *
401 : N1PTEL=1
402 : N2PTEL=0
403 : N2EL=0
404 : ELSE
405 : SEGINI,MCHELM=MCHEL1
406 : ENDIF
407 : *
408 : * Boucle sur les elements
409 : *
410 : DO IZONE=1,NZONE
411 : SEGINI,MCHAML=MCHAM1
412 : ICHAML(IZONE)=MCHAML
413 : IMODEL=KMODEL(IZONE)
414 : MELEME=IMAMOD
415 : N1EL=NUM(/2)
416 : DO J=1,MECAN
417 : SEGINI,MELVAL
418 : IELVAL(J)=MELVAL
419 : ENDDO
420 : DO J=1,N1EL
421 : *
422 : * Calcul de C1 elementaire
423 : *
423 : CALL DOXE(XCOOR,IDIM,2,NUM,J,XE)
423 : CALL TIPOCQ(XE, C,Q,MRIGI,2*NCOMP)
423 : CALL MATMAT(C,Q,MRIGI,2*NCOMP,2*NCOMP,C1)
424 : *
425 : * On remplit V1
426 : *
427 : * 1) noeud 1 de l'element
428 : *
429 : ICOLO=NCOMP*(ICPR(NUM(1,J))-1)
430 : DO K=1,NCOMP
431 : V1(K)=B1(ICOLO+K,I)
432 : ENDDO
433 : *
434 : * 2) noeud 2 de l'element
435 : *
436 : ICOLO=NCOMP*(ICPR(NUM(2,J))-1)
437 : DO K=1,NCOMP
438 : V1(K+NCOMP)=B1(ICOLO+K,I)
439 : ENDDO
440 : *
441 : * Calcul de THeta elementaire
442 : *
442 : CALL MATVE1(C1,V1,MRIGI,2*NCOMP,TH,2)
443 : *
444 : * Remplissage des MELVAL
445 : *
446 : DO K=1,MECAN
447 : MELVAL=IELVAL(K)
448 : VELCHE(1,J)=TH(K)
449 : ENDDO
450 : ENDDO
451 : DO J=1,MECAN
452 : MELVAL=IELVAL(J)
453 : SEGDES,MELVAL
454 : ENDDO
455 : SEGDES,MCHAML
456 : ENDDO
457 : SEGDES,MCHELM
458 : MTAB2.MTABIV(I+1)=MCHELM
459 : ENDDO
460 : SEGDES,MTAB2
461 : SEGSUP,MCHAM1
462 : *
463 : * A faire ...
464 : *
465 : * WARNING!!!!!!!!! SIGNE DES ROTATIONS!!!!!!! fait en modifiant C
466 : SEGSUP,WORK3
467 : SEGSUP,WORK0
468 : *
469 : * On rend la main a GIBIANE
470 : *
470 : CALL ECROBJ('TABLE',MTAB2)
470 : CALL ECROBJ('TABLE',MTAB1)
471 : *
472 : * Desactivation de la rigidite (imatri, descripteur et maillage)
473 : *
474 : 9997 DO I=1,NR
475 : MELEME=IRIGEL(1,I)
476 : SEGDES,MELEME
477 : DESCR=IRIGEL(3,I)
478 : SEGDES,DESCR
479 : xMATRI=IRIGEL(4,I)
480 : SEGDES,xMATRI
481 : ENDDO
482 : SEGDES,MRIGID
483 : *
484 : * Supression de ICPR
485 : *
486 : SEGSUP,ICPR
487 : *
488 : * Desactivation des maillages du modele
489 : *
490 : DO IZONE=1,NZONE
491 : IMODEL=KMODEL(IZONE)
492 : MELEME=IMAMOD
493 : SEGDES,MELEME
494 : ENDDO
495 : *
496 : * Desactivation du modele
497 : *
498 : 9999 DO IZONE=1,NZONE
499 : IMODEL=KMODEL(IZONE)
500 : SEGDES,IMODEL
501 : ENDDO
502 : SEGDES,MMODEL
503 : *
504 : RETURN
505 : END
506 :
507 :
508 :
509 :
© Cast3M 2003 - All rights reserved.
Disclaimer