OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ini_fxbody.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "units_c.inc"
#include "scr03_c.inc"
#include "scr05_c.inc"
#include "scr17_c.inc"
#include "fxbcom.inc"
#include "sysunit.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine ini_fxbody (fxbipm, fxbrpm, fxbnod, fxbglm, fxbcpm, fxbcps, fxblm, fxbfls, fxbdls, fxbmod, itab, x, ms, in, fxb_matrix, fxb_matrix_add, fxb_last_adress, icode, nom_opt)

Function/Subroutine Documentation

◆ ini_fxbody()

subroutine ini_fxbody ( integer, dimension(nbipm,*) fxbipm,
fxbrpm,
integer, dimension(*) fxbnod,
fxbglm,
fxbcpm,
fxbcps,
fxblm,
fxbfls,
fxbdls,
fxbmod,
integer, dimension(*) itab,
x,
ms,
in,
fxb_matrix,
integer, dimension(4,*) fxb_matrix_add,
integer, dimension(*) fxb_last_adress,
integer, dimension(*) icode,
integer, dimension(lnopt1,*) nom_opt )

Definition at line 39 of file ini_fxbody.F.

43C-----------------------------------------------
44C M o d u l e s
45C-----------------------------------------------
46 USE message_mod
47 USE inoutfile_mod
48 USE names_and_titles_mod , ONLY : nchartitle
49 USE format_mod , ONLY : fmw_10i
50C-----------------------------------------------
51C I m p l i c i t T y p e s
52C-----------------------------------------------
53#include "implicit_f.inc"
54C-----------------------------------------------
55C C o m m o n B l o c k s
56C-----------------------------------------------
57#include "com04_c.inc"
58#include "units_c.inc"
59#include "scr03_c.inc"
60#include "scr05_c.inc"
61#include "scr17_c.inc"
62#include "fxbcom.inc"
63#include "sysunit.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 INTEGER ITAB(*),FXBIPM(NBIPM,*),FXBNOD(*),FXB_MATRIX_ADD(4,*),FXB_LAST_ADRESS(*),ICODE(*),NOM_OPT(LNOPT1,*)
68 my_real x(3,*),ms(*),in(*),
69 . fxbrpm(*), fxbglm(*), fxbcpm(*), fxbcps(*),
70 . fxblm(*), fxbfls(*), fxbdls(*), fxbmod(*),
71 . fxb_matrix(*)
72C-----------------------------------------------
73C L o c a l V a r i a b l e s
74C-----------------------------------------------
75 INTEGER NFX,ID,IDMAST,NMOD,NMST,NBNO,NME,NTR,ADRGLM,
76 . ADRCP,ADRLM,ADRFLS,ADRDLS,ADRVAR,ADRRPM,ADRMCD,ADRMOD,IMOD,INO,I,LEN,
77 . NLIG,NRES,ILIG,ADRCP2,IR,ADRNOD,NUMNO(10),ISHELL,
78 . IC,J,INFO,IANIM, IMIN, IMAX, REF, ERR, K, L, IDN, I1,I2, J2, REF2,J1,SIZE_MAX,IDUM1,IDUM2,IDUM3,
79 . IDUM4,FLAG,IC1,IC2,BCS(6),SIZE_STIFF,ADR_STIFF,IL1,IL2,IDOF1,IDOF2,NMOD_MAX,NMST_F,NMR,
80 . ADR_MASS,SIZE_MASS,RELOOP,ADRMOD0,ADRLM0,ADRFLS0,LENCP0,LENLM0,LENFLS0,LENDLS0,LENVAR0,
81 . NSAV_MODES,NSAV_MODESN,NN
82 my_real
83 . freq,beta,omega,dtc1,dtc2,bid,xm(3),zz,rho_invmk,unsn,norm,fac,fac2,mstot,tole,
84 . kt,kr,rr,skew(9),rdum1,min_diag_stif,fac2x,fac2y,fac2z,dt_mini,dt_cst,omega_min,
85 . inmin,inmax,omega_max,max_freq,alpha,freq_range,min_cut_off_freq
86 CHARACTER(LEN=NCHARTITLE) :: TITR
87 CHARACTER :: MESS*40,MESS1*40,NWLINE*100,FXBFILE*100
88 INTEGER, DIMENSION(:,:), ALLOCATABLE :: TABSL,ITAG_DOF
89 INTEGER, DIMENSION(:), ALLOCATABLE :: INDEX,LISTN
90 INTEGER IFLAGI1,IFLAGDBL,IRB,CPT,REST,NLIGNE,M,INFO2,LWORK,NBDYN,NDDL,IM,NMOD0,NMST0
91C
92 INTEGER RPM_L,ADRMODE_L,ADRLM_L,ADRFLS_L,ADRGLM_L,ADRCP_L,ADRCP2_L,NSNGLR
93C
94 INTEGER NBNOD,IROT,IDAMP,IBLO,IFILE,REDUX,PRINTOUT
95 my_real, DIMENSION(:),ALLOCATABLE :: vt,vbid
96 my_real, DIMENSION(:,:),ALLOCATABLE :: mass,stif,modes,modest,modes_rb,modes_rbt,stift,
97 . mtemp,mtemp2,mtemp3,tt,t,vectr,modes_sav
98C
99 DOUBLE PRECISION WORK1(1000)
100 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: WORK,WR,WI
101 DOUBLE PRECISION, DIMENSION(:,:),ALLOCATABLE :: EIG_R,EIG_L,INVMK
102
103 INTEGER :: LEN_TMP_NAME
104 CHARACTER(len=2148) :: TMP_NAME
105
106C-----------------------------------------------
107C E x t e r n a l F u n c t i o n s
108C-----------------------------------------------
109 INTEGER USR2SYS
110 my_real
111 . prscal
112 DATA mess/'FLEXIBLE BODY : NODES '/
113 DATA mess1/'FLEXIBLE BODY DEFINITION '/
114C=====================================================================
115C
116 printout = 0
117 adrmod = fxb_last_adress(1)
118 adrglm = fxb_last_adress(2)
119 adrcp = fxb_last_adress(3)
120 adrcp2 = fxb_last_adress(3)
121 adrlm = fxb_last_adress(4)
122 adrfls = fxb_last_adress(5)
123 adrdls = fxb_last_adress(6)
124 adrvar = fxb_last_adress(7)
125 adrrpm = fxb_last_adress(8)
126 adrmcd = fxb_last_adress(9)
127C
128 freq_range = hundred
129C Default min cut off frequency 1000Hz
130 min_cut_off_freq = 1000.0 / fac_time
131C
132 DO nfx=1,nfxbody
133C
134 adrnod = fxbipm(6,nfx)
135C
136 IF (fxbipm(41,nfx)==2) THEN
137C
138 IF (printout == 0) WRITE(iout,1000)
139 printout = 1
140 CALL fretitl2(titr,nom_opt(lnopt1-ltitr+1,nfx),ltitr)
141C
142 id = fxbipm(1,nfx)
143 nbnod = fxbipm(3,nfx)
144 nmod = fxbipm(4,nfx)
145 nmst = fxbipm(5,nfx)
146 adrnod = fxbipm(6,nfx)
147 ishell = fxbipm(16,nfx)
148 nme = fxbipm(17,nfx)
149 iblo = fxbipm(28,nfx)
150 ifile = fxbipm(29,nfx)
151 ianim = fxbipm(36,nfx)
152 size_stiff = fxbipm(42,nfx)
153 size_mass = fxbipm(43,nfx)
154 adr_stiff = fxbipm(44,nfx)
155 adr_mass = fxbipm(45,nfx)
156 nddl = 6*nbnod
157 size_max = max(nmod,nme)
158C
159 ALLOCATE(itag_dof(6,nbnod))
160 ALLOCATE(stif(nddl,nddl),mass(nddl,nddl))
161 ALLOCATE(modes_rb(nddl,nme),modes_rbt(nme,nddl))
162 ALLOCATE(modest(nmod,nddl),modes(nddl,nmod))
163 ALLOCATE(vectr(nddl,6),mtemp(nddl,size_max),mtemp2(size_max,size_max))
164 ALLOCATE(mtemp3(size_max,size_max),vt(nddl))
165C
166CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
167CCCCC INITIAL SKEW CCCCC
168CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
169C
170 fxbipm(14,nfx) = adrrpm
171C
172 skew(1:9) = zero
173 skew(1) = one
174 skew(5) = one
175 skew(9) = one
176 DO i=1,9
177 fxbrpm(adrrpm+i) = skew(i)
178 ENDDO
179C
180 adrrpm=adrrpm+12
181C
182C-- Information of boundary nodes is not used for now
183 DO i=1,nbnod
184 IF (fxbnod(adrnod+i-1) < 0) THEN
185 fxbnod(adrnod+i-1) = abs(fxbnod(adrnod+i-1))
186 ENDIF
187 ENDDO
188 fxbipm(18,nfx) = nbnod
189C
190CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
191CCCCC STIFFNESS MATRIX CCCCC
192CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
193C
194 min_diag_stif = ep30
195 stif(1:nddl,1:nddl) = zero
196 DO i=1,size_stiff
197 i1 = fxb_matrix_add(1,adr_stiff+i-1)
198 i2 = fxb_matrix_add(2,adr_stiff+i-1)
199 idof1 = fxb_matrix_add(3,adr_stiff+i-1)
200 idof2 = fxb_matrix_add(4,adr_stiff+i-1)
201 stif(6*(i1-1)+idof1,6*(i2-1)+idof2) = fxb_matrix(adr_stiff+i-1)
202 stif(6*(i2-1)+idof2,6*(i1-1)+idof1) = fxb_matrix(adr_stiff+i-1)
203 IF ((6*(i2-1)+idof2)==(6*(i1-1)+idof1)) min_diag_stif = min(min_diag_stif,fxb_matrix(adr_stiff+i-1))
204 ENDDO
205C
206CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
207CCCCC MASS MATRIX CCCCC
208CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
209C
210 mass(1:nddl,1:nddl) = zero
211C
212 IF (size_mass == 0) THEN
213C-- Mass matrix automatically generated from nodal masses and inertia
214 DO i=1,nbnod
215 idn = fxbnod(adrnod+i-1)
216 fxbnod(adrnod+i-1) = idn
217 mass(6*(i-1)+1,6*(i-1)+1) = ms(idn)
218 mass(6*(i-1)+2,6*(i-1)+2) = ms(idn)
219 mass(6*(i-1)+3,6*(i-1)+3) = ms(idn)
220 mass(6*(i-1)+4,6*(i-1)+4) = in(idn)
221 mass(6*(i-1)+5,6*(i-1)+5) = in(idn)
222 mass(6*(i-1)+6,6*(i-1)+6) = in(idn)
223 ENDDO
224 ELSE
225 DO i=1,size_mass
226 i1 = fxb_matrix_add(1,adr_mass+i-1)
227 i2 = fxb_matrix_add(2,adr_mass+i-1)
228 idof1 = fxb_matrix_add(3,adr_mass+i-1)
229 idof2 = fxb_matrix_add(4,adr_mass+i-1)
230 mass(6*(i1-1)+idof1,6*(i2-1)+idof2) = fxb_matrix(adr_mass+i-1)
231 mass(6*(i2-1)+idof2,6*(i1-1)+idof1) = fxb_matrix(adr_mass+i-1)
232 ENDDO
233C-- Mass of connections nodes is added to the mass matrix
234 DO i=1,nbnod
235 idn = fxbnod(adrnod+i-1)
236 fxbnod(adrnod+i-1) = idn
237 mass(6*(i-1)+1,6*(i-1)+1) = mass(6*(i-1)+1,6*(i-1)+1)+ms(idn)
238 mass(6*(i-1)+2,6*(i-1)+2) = mass(6*(i-1)+2,6*(i-1)+2)+ms(idn)
239 mass(6*(i-1)+3,6*(i-1)+3) = mass(6*(i-1)+3,6*(i-1)+3)+ms(idn)
240 mass(6*(i-1)+4,6*(i-1)+4) = mass(6*(i-1)+4,6*(i-1)+4)+in(idn)
241 mass(6*(i-1)+5,6*(i-1)+5) = mass(6*(i-1)+5,6*(i-1)+5)+in(idn)
242 mass(6*(i-1)+6,6*(i-1)+6) = mass(6*(i-1)+6,6*(i-1)+6)+in(idn)
243 ENDDO
244 ENDIF
245
246CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
247CCCCC MODES RBODY CCCCC
248CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
249C
250 fxbipm(7,nfx) = adrmod
251C
252 IF (iblo == 0) THEN
253C
254 xm = zero
255 mstot = zero
256 DO i=1,nbnod
257 idn = fxbnod(adrnod+i-1)
258 xm(1) = xm(1) + x(1,idn)*ms(idn)
259 xm(2) = xm(2) + x(2,idn)*ms(idn)
260 xm(3) = xm(3) + x(3,idn)*ms(idn)
261 mstot = mstot + ms(idn)
262 ENDDO
263 DO i=1,3
264 xm(i) = xm(i)/max(em20,mstot)
265 ENDDO
266C
267 modes_rb(1:nddl,1:nme) = zero
268 modes_rbt(1:nme,1:nddl) = zero
269C
270C---------- Modes 1/.../9
271 DO i=1,3
272 DO j=1,3
273 DO k=1,nbnod
274 idn = fxbnod(adrnod+k-1)
275 modes_rb(6*(k-1)+j,3*(i-1)+j) = x(i,idn) - xm(i)
276 ENDDO
277 ENDDO
278 ENDDO
279C---------- Modes 10/11/12
280 DO i=1,3
281 DO k=1,nbnod
282 idn = fxbnod(adrnod+k-1)
283 modes_rb(6*(k-1)+i,9+i) = 1-(x(1,idn)-xm(1))-(x(2,idn)-xm(2))-(x(3,idn)-xm(3))
284 ENDDO
285 ENDDO
286C---------- Modes rotation
287 IF (nme > 12) THEN
288 DO i=1,3
289 DO k=1,nbnod
290 modes_rb(6*(k-1)+3+i,12+i) = one
291 ENDDO
292 ENDDO
293 ENDIF
294C --------- Matrice modes RB transposee
295 DO i=1,nme
296 DO j=1,nddl
297 modes_rbt(i,j) = modes_rb(j,i)
298 ENDDO
299 ENDDO
300C
301 DO i=1,nme
302 DO j=1,nddl
303 fxbmod(adrmod) = modes_rb(j,i)
304 adrmod=adrmod+1
305 ENDDO
306 ENDDO
307C
308C---------- SCALE FACTOR ON LOCAL MODES FOR MASS CONDITIONING
309 fac2x=prscal(modes_rb(1,1), modes_rb(1,1), nddl, mass)
310 fac2y=prscal(modes_rb(1,4), modes_rb(1,4), nddl, mass)
311 fac2z=prscal(modes_rb(1,7), modes_rb(1,7), nddl, mass)
312 fac2 = max(fac2x,fac2y,fac2z)
313 fac=sqrt(fac2)
314C
315 ELSE
316C
317 fac2 = one
318 fac=sqrt(fac2)
319C
320 ENDIF
321C
322CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
323CCCCC STATIC MODES CCCCC
324CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
325C
326C
327 modes(1:nddl,1:nmod) = zero
328 modest(1:nmod,1:nddl) = zero
329 itag_dof(1:6,1:nbnod) = 0
330C
331C--- Tag of connected DOF ---
332 DO i=1,size_stiff
333 il1 = fxb_matrix_add(1,adr_stiff+i-1)
334 il2 = fxb_matrix_add(2,adr_stiff+i-1)
335 idof1 = fxb_matrix_add(3,adr_stiff+i-1)
336 idof2 = fxb_matrix_add(4,adr_stiff+i-1)
337C -- Craig-bampton with only boundary nodes - one static mode per dof in connected
338 itag_dof(idof1,il1) = 1
339 itag_dof(idof2,il2) = 1
340 ENDDO
341C
342 nmst = 0
343 DO i=1,nbnod
344 idn = fxbnod(adrnod+i-1)
345C--- Tag BCS ---
346 ic = icode(idn)
347 ic1=ic/512
348 ic2=(ic-512*ic1)/64
349 bcs(1) = ic1/4
350 bcs(2) = (ic1-4*bcs(1))/2
351 bcs(3) = ic1-4*bcs(1)-2*bcs(2)
352 bcs(4) = ic2/4
353 bcs(5) = (ic2-4*bcs(4))/2
354 bcs(6) = ic2-4*bcs(4)-2*bcs(5)
355
356C--- On static mode generated per connected DOF ---
357 DO k=1,6
358 IF ((bcs(k)==0).AND.(itag_dof(k,i)>0)) THEN
359 nmst = nmst + 1
360 modes(6*(i-1)+k,nmst) = one
361 ENDIF
362 ENDDO
363 ENDDO
364C
365CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
366CCCCC COMPLEMENT MODES - used for othogonalisation with RB modes CCCCC
367CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
368C
369 nmr = 0
370 IF (iblo == 0) THEN
371C
372 cpt = 0
373 DO i=1,nbnod
374 idn = fxbnod(adrnod+i-1)
375 DO j=1,6
376 DO im=1,6
377 vectr(cpt+j,im)=zero
378 ENDDO
379 SELECT CASE (j)
380 CASE (1)
381 vectr(cpt+j,1)=one
382 vectr(cpt+j,5)=x(3,idn)-xm(3)
383 vectr(cpt+j,6)=-(x(2,idn)-xm(2))
384 CASE (2)
385 vectr(cpt+j,2)=one
386 vectr(cpt+j,4)=-(x(3,idn)-xm(3))
387 vectr(cpt+j,6)=x(1,idn)-xm(1)
388 CASE (3)
389 vectr(cpt+j,3)=one
390 vectr(cpt+j,4)=x(2,idn)-xm(2)
391 vectr(cpt+j,5)=-(x(1,idn)-xm(1))
392 CASE (4)
393 vectr(cpt+j,4)=one
394 CASE (5)
395 vectr(cpt+j,5)=one
396 CASE (6)
397 vectr(cpt+j,6)=one
398 END SELECT
399 ENDDO
400 cpt=cpt+6
401 ENDDO
402C
403C--- Orthonormalisation of VECTR with mass matrix ---
404 nmr = 6
405 IF (nmst > 0) CALL orthrg(vectr,mass,nddl,nmr)
406C
407 ENDIF
408c
409 IF ((iblo==0).AND.(nmst > 0)) THEN
410C--- Orthonormalisation of VECTR and vector of static modes with mass matrix ---
411 CALL orthsr(modes,vectr,mass,nddl,nmst,nmr)
412 ENDIF
413C
414CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
415CCCCCCCCCCCCCCCC LOOP -> reloop = 0 - all the modes CCCCCCCCCCCCCCCCC
416CCCCCCCCCCCCCCCC reloop = 1 - reduced modal bases CCCCCCCCCCCCCCCCC
417CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
418C
419 reloop = 0
420 nsav_modes = 0
421 nmst0 = nmst
422 nmod0 = nmod
423 adrmod0 = adrmod
424 adrlm0 = adrlm
425 adrfls0 = adrfls
426 lencp0 =lencp
427 lenlm0 =lenlm
428 lenfls0 =lenfls
429 lendls0 =lendls
430 lenvar0 =lenvar
431C
432 DO reloop = 0,1
433C
434 IF (reloop > 0) THEN
435 nmst = nsav_modes
436 nmod = nsav_modes
437 adrmod = adrmod0
438 adrlm = adrlm0
439 adrfls = adrfls0
440 lencp =lencp0
441 lenlm =lenlm0
442 lenfls =lenfls0
443 lendls =lendls0
444 lenvar =lenvar0
445 modes(1:nddl,1:nmod) = modes_sav(1:nddl,1:nmod)
446 ENDIF
447C
448CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
449CCCCC STATIC MODES - ORTHONORMALIZATION CCCCC
450CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
451C
452 IF (nmst > 0) THEN
453 nmst_f = 0
454C TOLE = MIN_DIAG_STIF*MIN_DIAG_STIF*EM5*EM5
455 tole = em04
456 CALL orthst(modes,mass,nddl,nmst,nmst_f,tole)
457 nmst = nmst_f
458 ENDIF
459
460C-- Number of static modes is updated
461 nmod = nmst
462 fxbipm(4,nfx) = nmod
463 fxbipm(5,nfx) = nmst
464C
465 ntr = 9
466 lencp =lencp -ntr*nmod0*nme+ntr*nmod*nme
467 lenlm =lenlm -nmod0+nmod
468 lenfls=lenfls-nmst0*(2*nmod0-nmst0+1)/2+nmst*(2*nmod-nmst+1)/2
469 lendls=lendls-nmod0+nmst0+nmod+nmst
470 lenvar=lenvar-nmod0+nmod
471 fxbipm(38,nfx)=min(nmod,fxbipm(38,nfx))
472C
473 DO i=1,nddl
474 DO j=1,nmod
475 modest(j,i)=modes(i,j)
476 ENDDO
477 ENDDO
478C
479 DO i=1,nmst
480 DO j=1,nddl
481 fxbmod(adrmod) = fac*modes(j,i)
482 adrmod=adrmod+1
483 ENDDO
484 ENDDO
485C
486CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
487CCCCC REDUCED DIAG MASS MATRIX CCCCC
488CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
489C
490 fxbipm(10,nfx) = adrlm
491C
492 mtemp(1:nddl,1:nmod) = matmul(mass(1:nddl,1:nddl),modes(1:nddl,1:nmod))
493 mtemp2(1:nmod,1:nmod) = matmul(modest(1:nmod,1:nddl),mtemp(1:nddl,1:nmod))
494C
495 DO i=1,nmod
496 fxblm(adrlm) = fac2*mtemp2(i,i)
497 adrlm=adrlm+1
498 ENDDO
499C
500CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
501CCCCC REDUCED FULL STIFF MATRIX CCCCC
502CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
503C
504 fxbipm(11,nfx) = adrfls
505C
506 mtemp(1:nddl,1:nmst) = matmul(stif(1:nddl,1:nddl),modes(1:nddl,1:nmst))
507 mtemp3(1:nmst,1:nmst) = matmul(modest(1:nmst,1:nddl),mtemp(1:nddl,1:nmst))
508C
509 DO i=1,nmst
510 DO j=1,nmst
511CC-- Mat sym trian sup
512 IF (j >= i) THEN
513 fxbfls(adrfls) = fac2*mtemp3(i,j)
514 adrfls=adrfls+1
515 ENDIF
516 ENDDO
517 ENDDO
518C
519CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
520CCCCC COMPUTATION OF MAX FREQUENCY CCCCC
521CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
522C
523 ALLOCATE(wi(nmst),wr(nmst),eig_r(nmst,nmst),eig_l(nmst,nmst),invmk(nmst,nmst))
524C
525 DO i=1,nmst
526 DO j=1,nmst
527 IF (i==j) THEN
528 mtemp2(i,j) = one/max(em20,mtemp2(i,j))
529 ELSE
530 mtemp2(i,j) = zero
531 ENDIF
532 ENDDO
533 ENDDO
534C
535 invmk(1:nmst,1:nmst) = matmul(mtemp2(1:nmst,1:nmst),mtemp3(1:nmst,1:nmst))
536C
537#ifndef WITHOUT_LINALG
538 lwork = -1
539 CALL dgeev( 'N','V',nmst,invmk,nmst,wr,wi,eig_l,nmst,eig_r,nmst,work1,lwork,info2)
540 lwork = max( 1000, int( work1( 1 ) ) )
541 ALLOCATE(work(lwork))
542C
543 CALL dgeev( 'N','V',nmst,invmk,nmst,wr,wi,eig_l,nmst,eig_r,nmst,work,lwork,info2)
544#else
545 info2 = 1
546#endif
547C
548 IF( info2>0 ) THEN
549 WRITE(*,*)'The algorithm failed to compute eigenvalues.'
550 stop
551 END IF
552C
553 rho_invmk = zero
554 omega_min = ep30
555 DO i=1,nmst
556 IF (wi(i)==zero) THEN
557 omega_min = min(omega_min,sqrt(wr(i)))
558 ENDIF
559 ENDDO
560C
561 omega_max = max(freq_range*omega_min,two*pi*min_cut_off_freq)
562 dt_mini = two/omega_max
563 nsav_modesn = 0
564 DO i=1,nmst
565 IF (wi(i)==zero) THEN
566 rho_invmk = max(rho_invmk,abs(wr(i)))
567 IF ((two/sqrt(abs(wr(i))) > dt_mini)) THEN
568 nsav_modesn = nsav_modesn+1
569 ENDIF
570 ENDIF
571 ENDDO
572C
573 IF ((nmst - nsav_modesn > 0).AND.(reloop==0)) THEN
574 nsav_modes = nsav_modesn
575 ALLOCATE(modes_sav(nddl,nsav_modes))
576 modes_sav(1:nddl,1:nsav_modes) = zero
577 nn = 0
578 DO i=1,nmst
579 IF ((two/sqrt(abs(wr(i))) > dt_mini).AND.((wi(i)==zero))) THEN
580 nn = nn + 1
581 DO j=1,nddl
582 DO k=1,nmst
583 modes_sav(j,nn)=modes_sav(j,nn)+modes(j,k)*eig_r(k,i)
584 ENDDO
585 ENDDO
586 ENDIF
587 ENDDO
588 ENDIF
589C
590 DEALLOCATE(wi,wr,eig_r,eig_l,work,invmk)
591C
592 fxbrpm(fxbipm(14,nfx))=zep9*two/sqrt(rho_invmk)
593C
594 IF (nmst - nsav_modesn == 0) EXIT
595C
596CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
597CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
598CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
599C
600 ENDDO
601C
602 IF (ALLOCATED(modes_sav)) DEALLOCATE (modes_sav)
603C
604CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
605CCCCC DAMPING CCCCC
606CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
607C
608 idamp = 0
609C
610 IF (idamp == 0) THEN
611 fxbrpm(adrrpm)=zero
612 fxbrpm(adrrpm+1)=zero
613 adrrpm=adrrpm+2
614 ENDIF
615 fxbrpm(adrrpm)=zero
616 fxbrpm(adrrpm+1)=zero
617 adrrpm=adrrpm+2
618C
619CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
620CCCCC MASS MATRIX PROJ ON RB MODES CCCCC
621CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
622C
623 fxbipm(8,nfx) = adrglm
624C
625 mtemp(1:nddl,1:nme) = matmul(mass(1:nddl,1:nddl),modes_rb(1:nddl,1:nme))
626 mtemp2(1:nme,1:nme) = matmul(modes_rbt(1:nme,1:nddl),mtemp(1:nddl,1:nme))
627C
628 DO i=1,nme
629 DO j=1,nme
630CC-- Mat sym trian sup
631 IF (j >= i) THEN
632 fxbglm(adrglm) = mtemp2(i,j)
633 adrglm=adrglm+1
634 ENDIF
635 ENDDO
636 ENDDO
637C
638CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
639CCCCC MASS MATRIX COUPLED PROJECTION CCCCC
640CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
641C
642 fxbipm(9,nfx) = adrcp
643C
644 DO i = 1,3
645 DO j = 1,3
646C
647 DO k=1,nmod
648 vt(1:nddl) = zero
649 cpt = 0
650 DO m = 1,nbnod
651 vt(cpt+i) = modes(cpt+j,k)
652 vt(cpt+3+i) = modes(cpt+3+j,k)
653 cpt = cpt + 6
654 ENDDO
655 DO l=1,nme
656 mtemp(l,k)=fac*prscal(vt,modes_rb(1,l),nddl,mass)
657 ENDDO
658 ENDDO
659C
660 DO l=1,nme
661 DO k=1,nmod
662 fxbcpm(adrcp) = mtemp(l,k)
663 adrcp=adrcp+1
664 ENDDO
665 ENDDO
666C
667 ENDDO
668 ENDDO
669C
670CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
671CCCCC STIFF MATRIX COUPLED PROJECTION CCCCC
672CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
673C
674 fxbipm(9,nfx) = adrcp2
675C
676 DO i = 1,3
677 DO j = 1,3
678C
679 DO k=1,nme
680 vt(1:nddl) = zero
681 cpt = 0
682 DO m = 1,nbnod
683 vt(cpt+i) = modes_rb(cpt+j,k)
684 vt(cpt+3+i) = modes_rb(cpt+3+j,k)
685 cpt = cpt + 6
686 ENDDO
687 DO l=1,nmod
688 mtemp(k,l)=fac*prscal(modes(1,l),vt,nddl,stif)
689 ENDDO
690 ENDDO
691C
692 DO l=1,nme
693 DO k=1,nmod
694 fxbcps(adrcp2) = mtemp(l,k)
695 adrcp2=adrcp2+1
696 ENDDO
697 ENDDO
698C
699 ENDDO
700 ENDDO
701C
702CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
703CCCCC ADDITIONAL ADDRESSES CCCCC
704CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
705C
706 fxbipm(13,nfx) = adrvar
707 adrvar=adrvar+nmod+nme
708 fxbipm(15,nfx) = adrmcd
709 adrmcd=adrmcd+nme*nme
710C
711CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
712CCC PRINTOUT CONTROL CCCCC
713CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
714C
715 WRITE(iout,1100) id,trim(titr),nmst,nme,nbnod,fxbrpm(fxbipm(14,nfx))
716 WRITE(iout,1200)
717 WRITE(iout,fmt=fmw_10i)(itab(fxbnod(adrnod+i-1)),i=1,nbnod)
718C
719 IF (ipri >= 5) THEN
720C
721C-- FXB file generation for full output
722C
723 rpm_l = fxbipm(14,nfx)
724 adrmode_l = fxbipm(7,nfx)
725 adrlm_l = fxbipm(10,nfx)
726 adrfls_l = fxbipm(11,nfx)
727 adrglm_l = fxbipm(8,nfx)
728 adrcp_l = fxbipm(9,nfx)
729 adrcp2_l = fxbipm(9,nfx)
730C
731 ref= 11000
732 IF (nfx < 10) THEN
733 WRITE(tmp_name,'(A,I1,A)') outfile_name(1:outfile_name_len)//"FXB_CONTROL",nfx,"_0000.fxb"
734 len_tmp_name = outfile_name_len + 21
735 ELSE
736 WRITE(tmp_name,'(A,I2,A)') outfile_name(1:outfile_name_len)//"FXB_CONTROL",nfx,"_0000.fxb"
737 len_tmp_name = outfile_name_len + 22
738 ENDIF
739C
740 OPEN(unit=ref,file=tmp_name(1:len_tmp_name),
741 . access='SEQUENTIAL',form='FORMATTED',status='UNKNOWN')
742C
743 WRITE(ref,'(A)') '# FLEXIBLE BODY DATA'
744 WRITE(ref,'(A)') '#--1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|--10---|'
745 WRITE(ref,'(A)') '# Param'
746 WRITE(ref,'(A)') "# Nbmod Nbstat Nbnod Irot Idamp Iblo Ifile Intype"
747C
748 WRITE(ref,'(7I8)') nmod,nmst,nbnod,ishell,idamp,iblo,ifile
749 WRITE(ref,'(A)') '# Nodes'
750C
751 cpt = 0
752 zz = nbnod / 10
753 nligne = ceiling(zz)
754 DO i=1,nligne
755 WRITE(ref,'(10I8)') (itab(fxbnod(adrnod+cpt+k-1)),k=1,10)
756 cpt = cpt + 10
757 ENDDO
758 rest = nbnod-10*nligne
759 IF (rest > 0) WRITE(ref,'(10I8)') (itab(fxbnod(adrnod+cpt+k-1)),k=1,rest)
760C
761 max_freq = half*omega_max/pi
762 WRITE(ref,'(A)') '# Mrot11 Mrot12 Mrot13 Mrot21 Mrot22'
763 WRITE(ref,'(5(1PE16.9))') fxbrpm(rpm_l+1),fxbrpm(rpm_l+2),fxbrpm(rpm_l+3),fxbrpm(rpm_l+4),fxbrpm(rpm_l+5)
764 WRITE(ref,'(A)') '# Mrot23 Mrot31 Mrot32 Mrot33 Freq'
765 WRITE(ref,'(5(1PE16.9))') fxbrpm(rpm_l+6),fxbrpm(rpm_l+7),fxbrpm(rpm_l+8),fxbrpm(rpm_l+9),max_freq
766C
767C --- block for damping
768C
769 IF (idamp == 1) THEN
770 WRITE(ref,'(A)') '# RAYLEIGH DAMPING'
771 WRITE(ref,'(A)') ' '
772 WRITE(ref,'(A)') ' '
773 WRITE(ref,'(2(1PE16.9))') fxbrpm(adrrpm-4),fxbrpm(adrrpm-3)
774 ENDIF
775C
776 WRITE(ref,'(A)') '# BLOCK5 - RB MODES'
777C
778C --- block5 - Modes RBODY
779C
780 IF (iblo == 0) THEN
781C
782 DO i=1,nme
783 WRITE(ref,'(A)') '# X Y Z XX YY'
784 WRITE(ref,'(A)') '# ZZ'
785 DO j=1,nbnod
786 WRITE(ref,'(5(1PE16.9))') (fxbmod(adrmode_l+k-1),k=1,5)
787 WRITE(ref,'((1PE16.9))') fxbmod(adrmode_l+6-1)
788 adrmode_l = adrmode_l + 6
789 ENDDO
790 ENDDO
791C
792 ENDIF
793C
794 WRITE(ref,'(A)') '# BLOCK6 - REDUCED MODES ROTATION '
795 WRITE(ref,'(A)') '# X Y Z XX YY'
796 WRITE(ref,'(A)') '# ZZ'
797C
798C --- block7 - Reduced Modes
799C
800 WRITE(ref,'(A)') '# BLOCK7 - REDUCED MODES'
801C
802 DO i=1,nmst
803 WRITE(ref,'(A)') '# X Y Z XX YY'
804 WRITE(ref,'(A)') '# ZZ'
805 DO j=1,nbnod
806 WRITE(ref,'(5(1PE16.9))') (fxbmod(adrmode_l+k-1),k=1,5)
807 WRITE(ref,'((1PE16.9))') fxbmod(adrmode_l+6-1)
808 adrmode_l = adrmode_l + 6
809 ENDDO
810 ENDDO
811C
812C --- block8 - Diag Mass matrix
813C
814 WRITE(ref,'(A)') '# BLOCK8 - REDUCED DIAG MASS MATRIX'
815 WRITE(ref,'(A)') '# X Y Z XX YY'
816 WRITE(ref,'(A)') '# ZZ'
817C
818 zz = nmod / 5
819 nligne = ceiling(zz)
820 DO i=1,nligne
821 WRITE(ref,'(5(1PE16.9))') (fxblm(adrlm_l+k-1),k=1,5)
822 adrlm_l = adrlm_l + 5
823 ENDDO
824 rest = nmod-5*nligne
825 IF (rest > 0) THEN
826 WRITE(ref,'(5(1PE16.9))') (fxblm(adrlm_l+k-1),k=1,rest)
827 adrlm_l = adrlm_l + 1
828 ENDIF
829C
830C --- block9 - Full Stiff matrix
831C
832 WRITE(ref,'(A)') '# BLOCK9 - REDUCED FULL STIFF MATRIX'
833 WRITE(ref,'(a)') '# X Y Z XX YY'
834 WRITE(ref,'(A)') '# ZZ'
835C
836 size_max = ((nmst+1)*nmst)/2
837 zz = size_max / 5
838 nligne = ceiling(zz)
839 DO i=1,nligne
840 WRITE(ref,'(5(1PE16.9))') (fxbfls(adrfls_l+k-1),k=1,5)
841 adrfls_l = adrfls_l + 5
842 ENDDO
843 rest = size_max-5*nligne
844 IF (rest > 0) THEN
845 WRITE(ref,'(5(1PE16.9))') (fxbfls(adrfls_l+k-1),k=1,rest)
846 adrfls_l = adrfls_l + 1
847 ENDIF
848C
849C --- block11 - Reduced Mass matrix on rbody modes - diag
850C
851 WRITE(ref,'(A)') '# BLOCK11 - RB MODES RB MATRIX'
852 WRITE(ref,'(A)') '# X Y Z XX YY'
853 WRITE(ref,'(A)') '# ZZ'
854C
855 size_max = ((nme+1)*nme)/2
856 zz = size_max / 5
857 nligne = ceiling(zz)
858 DO i=1,nligne
859 WRITE(ref,'(5(1PE16.9))') (fxbglm(adrglm_l+k-1),k=1,5)
860 adrglm_l = adrglm_l + 5
861 ENDDO
862 rest = size_max-5*nligne
863 IF (rest > 0) THEN
864 WRITE(ref,'(5(1PE16.9))') (fxbglm(adrglm_l+k-1),k=1,rest)
865 adrglm_l = adrglm_l + 1
866 ENDIF
867C
868C --- block12 - Coupled mass projection
869C
870 WRITE(ref,'(A)') '# BLOCK12 - COUPLE MASS PROJECTION'
871 WRITE(ref,'(A)') '# X Y Z XX YY'
872 WRITE(ref,'(A)') '# ZZ'
873C
874 DO l=1,9
875 WRITE(ref,'(A,7I8)') '#',l
876 size_max = nme*nmod
877 zz = size_max / 5
878 nligne = ceiling(zz)
879 DO i=1,nligne
880 WRITE(ref,'(5(1PE16.9))') (fxbcpm(adrcp_l+k-1),k=1,5)
881 adrcp_l = adrcp_l + 5
882 ENDDO
883 rest = size_max-5*nligne
884 IF (rest > 0) THEN
885 WRITE(ref,'(5(1PE16.9))') (fxbcpm(adrcp_l+k-1),k=1,rest)
886 adrcp_l = adrcp_l + 1
887 ENDIF
888 ENDDO
889C
890C --- block13 - Coupled stiff projection
891C
892 WRITE(ref,'(A)') '# BLOCK13 - COUPLE STIFF PROJECTION'
893 WRITE(ref,'(A)') '# X Y Z XX YY'
894 WRITE(ref,'(A)') '# ZZ'
895C
896 DO l=1,9
897 WRITE(ref,'(A,7I8)') '#',l
898 size_max = nme*nmod
899 zz = size_max / 5
900 nligne = ceiling(zz)
901 DO i=1,nligne
902 WRITE(ref,'(5(1PE16.9))') (fxbcps(adrcp2_l+k-1),k=1,5)
903 adrcp2_l = adrcp2_l + 5
904 ENDDO
905 rest = size_max-5*nligne
906 IF (rest > 0) THEN
907 WRITE(ref,'(5(1PE16.9))') (fxbcps(adrcp2_l+k-1),k=1,rest)
908 adrcp2_l = adrcp2_l + 1
909 ENDIF
910 ENDDO
911C
912 WRITE(ref,'(A)') '#--1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|--10---|'
913C
914 CLOSE(unit=ref,status='KEEP')
915C
916CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
917CCC FIN PRINTOUT CONTROL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
918CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
919C
920 ENDIF
921C
922 DEALLOCATE(itag_dof,stif,mass,modes_rb,modes_rbt)
923 DEALLOCATE(modest,modes,vectr,mtemp,mtemp2,mtemp3,vt)
924C
925 ENDIF
926C
927 ENDDO
928C
929 RETURN
930 CALL ancmsg(msgid=566,
931 . msgtype=msgerror,
932 . anmode=aninfo,
933 . i1=id,
934 . c1=titr,
935 . c2=fxbfile,
936 . c3=nwline)
937C
938 RETURN
939C
9401000 FORMAT(/
941 . ' FLEXIBLE BODY INITIALIZATION FROM DMIG'/
942 . ' ---------------------- ')
943C
9441100 FORMAT( /6x,'FLEXIBLE BODY ID ',i10,1x,a
945 . /10x,'NUMBER OF MODES ',i10
946 . /10x,'NUMBER OF RIGID BODY MODES ',i10
947 . /10x,'NUMBER OF NODES ',i10
948 . /10x,'STABILITY TIME STEP ',1pe10.3)
949C
9501200 FORMAT(/
951 . ' LIST OF NODES')
952C
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha
#define alpha
Definition eval.h:35
dgeev
subroutine dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition dgeev.f:192
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
inoutfile_mod
Definition inoutfile_mod.F:96
inoutfile_mod::outfile_name
character(len=outfile_char_len) outfile_name
Definition inoutfile_mod.F:111
inoutfile_mod::outfile_name_len
integer outfile_name_len
Definition inoutfile_mod.F:109
message_mod
Definition message_mod.F:1249
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003
orthsr
subroutine orthsr(vects, vectr, mas, nddl, nms, nmr)
Definition ortho_normalization.F:80
orthrg
subroutine orthrg(vect, mas, nddl, nb_modes)
Definition ortho_normalization.F:32
orthst
subroutine orthst(vects, mas, nddl, nms, nmsf, tole)
Definition ortho_normalization.F:125
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
fretitl2
subroutine fretitl2(titr, iasc, l)
Definition freform.F:804