OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
q4forc2.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| q4forc2 ../engine/source/elements/solid_2d/quad4/q4forc2.f
25!||--- called by ------------------------------------------------------
26!|| forint ../engine/source/elements/forint.F
27!||--- calls -----------------------------------------------------
28!|| mmain ../engine/source/materials/mat_share/mmain.F90
29!|| q4coor2 ../engine/source/elements/solid_2d/quad4/q4coor2.F
30!|| q4cumu2 ../engine/source/elements/solid_2d/quad4/q4cumu2.F
31!|| q4cumu2p ../engine/source/elements/solid_2d/quad4/q4cumu2p.f
32!|| q4defo2 ../engine/source/elements/solid_2d/quad4/q4defo2.F
33!|| q4defoc2 ../engine/source/elements/solid_2d/quad4/q4defoc2.f
34!|| q4deri2 ../engine/source/elements/solid_2d/quad4/q4deri2.F
35!|| q4deric2 ../engine/source/elements/solid_2d/quad4/q4deric2.F
36!|| q4fint2 ../engine/source/elements/solid_2d/quad4/q4fint2.F
37!|| q4fintc2 ../engine/source/elements/solid_2d/quad4/q4fintc2.F
38!|| q4rcoor2 ../engine/source/elements/solid_2d/quad4/q4rcoor2.F
39!|| q4rrota2 ../engine/source/elements/solid_2d/quad4/q4rrota2.F
40!|| q4vis2 ../engine/source/elements/solid_2d/quad4/q4vis2.F
41!|| q4zero2 ../engine/source/elements/solid_2d/quad4/q4zero2.F
42!|| qbilan ../engine/source/elements/solid_2d/quad/qbilan.F
43!|| qdlen2 ../engine/source/elements/solid_2d/quad/qdlen2.f
44!|| qlagr2 ../engine/source/elements/solid_2d/quad/qlagr2.F
45!|| qmass2 ../engine/source/elements/solid_2d/quad/qmass2.F
46!|| qmass2p ../engine/source/elements/solid_2d/quad/qmass2p.F
47!|| qmassreal2 ../engine/source/elements/solid_2d/quad/qmassreal2.F
48!|| qmassreal2p ../engine/source/elements/solid_2d/quad/qmassreal2p.F
49!|| qrota2 ../engine/source/elements/solid_2d/quad/qrota2.F
50!|| qvolu2 ../engine/source/elements/solid_2d/quad/qvolu2.F
51!|| s8efmoy3 ../engine/source/elements/solid/solide8e/s8efmoy3.F
52!|| s8zsigp3 ../engine/source/elements/solid/solide8z/s8zsigp3.F
53!|| smallb3 ../engine/source/elements/solid/solide/smallb3.F
54!||--- uses -----------------------------------------------------
55!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
56!|| dt_mod ../engine/source/modules/dt_mod.F
57!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
58!|| glob_therm_mod ../common_source/modules/mat_elem/glob_therm_mod.F90
59!|| mat_elem_mod ../common_source/modules/mat_elem/mat_elem_mod.F90
60!|| mmain_mod ../engine/source/materials/mat_share/mmain.F90
61!|| nlocal_reg_mod ../common_source/modules/nlocal_reg_mod.F
62!|| output_mod ../common_source/modules/output/output_mod.F90
63!|| sensor_mod ../common_source/modules/sensor_mod.F90
64!|| table_mod ../engine/share/modules/table_mod.f
65!|| timer_mod ../engine/source/system/timer_mod.F90
66!||====================================================================
67 SUBROUTINE q4forc2(TIMERS, OUTPUT, PM ,GEO ,IC ,X ,A ,
68 2 V ,MS ,W ,FLUX ,FLU1 ,
69 3 VEUL ,FV ,ALE_CONNECT ,IPARG ,NLOC_DMG,
70 4 ELBUF_TAB,TF ,NPF ,BUFMAT,PARTSAV,
71 5 DT2T ,NELTST ,ITYPTST,STIFN ,OFFSET ,
72 6 EANI ,IPARTQ ,NEL ,IADQ ,FSKY ,
73 7 ICP ,NG ,
74 8 IPM ,BUFVOIS,QMV ,GRESAV,GRTH ,
75 9 IGRTH ,TABLE ,IGEO ,ITASK ,IEXPAN ,
76 A MS_2D ,FSKYM ,IOUTPRT,MAT_ELEM,H3D_STRAIN,SZ_BUFVOIS,SNPC,
77 B STF,SBUFMAT,SVIS,NSVOIS,IDTMINS,IRESP,TT,DT1,
78 . IDEL7NG, IDEL7NOK, IDTMIN, MAXFUNC,IMON_MAT,
79 V USERL_AVAIL, impl_s, idyna,
80 . DT , GLOB_THERM,SENSORS)
81C-----------------------------------------------
82C M o d u l e s
83C-----------------------------------------------
84 USE timer_mod
85 USE output_mod, only : output_
86 USE mmain_mod
87 USE table_mod
88 USE mat_elem_mod
89 USE nlocal_reg_mod
90 USE ale_connectivity_mod
91 USE elbufdef_mod
92 USE dt_mod
93 use glob_therm_mod
94 USE sensor_mod
95C-----------------------------------------------
96C I m p l i c i t T y p e s
97C-----------------------------------------------
98#include "implicit_f.inc"
99C-----------------------------------------------
100C G l o b a l P a r a m e t e r s
101C-----------------------------------------------
102#include "mvsiz_p.inc"
103C-----------------------------------------------
104C C o m m o n B l o c k s
105C-----------------------------------------------
106#include "vect01_c.inc"
107#include "com01_c.inc"
108#include "com04_c.inc"
109#include "parit_c.inc"
110#include "param_c.inc"
111C-----------------------------------------------
112C D u m m y A r g u m e n t s
113C-----------------------------------------------
114 TYPE(timer_), INTENT(INOUT) :: TIMERS
115 TYPE(OUTPUT_), INTENT(INOUT) :: OUTPUT
116 my_real, INTENT(IN) :: DT1
117 my_real, INTENT(IN) :: TT
118 INTEGER, INTENT(IN) :: SNPC
119 INTEGER, INTENT(IN) :: STF
120 INTEGER, INTENT(IN) :: SBUFMAT
121 INTEGER, INTENT(IN) :: NSVOIS
122 INTEGER, INTENT(IN) :: IDTMINS
123 INTEGER ,INTENT(IN) :: IRESP
124 INTEGER ,INTENT(IN) :: IDEL7NG
125 INTEGER ,INTENT(INOUT) :: IDEL7NOK
126 integer,dimension(102) :: IDTMIN
127 INTEGER ,INTENT(IN) :: MAXFUNC
128 INTEGER, INTENT(IN) :: IMPL_S
129 INTEGER, INTENT(IN) :: IDYNA
130 INTEGER, INTENT(IN) :: USERL_AVAIL
131 INTEGER, INTENT(IN) :: IMON_MAT
132 INTEGER IC(*), IPARG(NPARG,*), NPF(*), IPARTQ(*),
133 + ipm(*), grth(*), igrth(*), igeo(*), iadq(4,*),itask,ioutprt
134 INTEGER OFFSET, NEL, NELTST, ITYPTST, ICP, NG, IEXPAN,H3D_STRAIN
135 my_real
136 + dt2t
137 my_real
138 + pm(npropm,*), geo(*), x(*),a(*),v(3,*),ms(*), w(*),partsav(*),
139 + flux(4,*),flu1(*),veul(*),fv(*), tf(*),bufmat(*),
140 + fsky(*),stifn(*),eani(*),bufvois(6,*),qmv(8,*),gresav(*),ms_2d(*),
141 + fskym(*)
142 my_real, DIMENSION(MVSIZ,6), INTENT(INOUT) :: svis
143 TYPE (TTABLE) TABLE(*)
144 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP), TARGET :: ELBUF_TAB
145 TYPE (NLOCAL_STR_) , TARGET :: NLOC_DMG
146 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
147 TYPE (MAT_ELEM_) ,INTENT(INOUT) :: MAT_ELEM
148 TYPE (DT_), INTENT(IN) :: DT
149 type (glob_therm_) ,intent(inout) :: glob_therm
150 type (sensors_),INTENT(INOUT) :: SENSORS
151C-----------------------------------------------
152c FUNCTION: Internal force computation of fully-integrated 2D Quad4 element
153c ARGUMENTS: (I: input, O: output, IO: input & output, W: workspace)
154c TYPE NAME FUNCTION
155c I PM ,GEO Material and geometrical property data
156c I IC(7,NUM_QUAD) connectivity and mid,pid integer data
157c I X(3,NUMNOD) co-ordinate
158c IO A(3,NUMNOD) nodal internal force
159c I V(3,NUMNOD) nodal velocity
160c IO MS(NUMNOD) nodal masse
161c I FLUX(4,NEL) flux at each side used w/ ALE or EULER
162c I FLU1,VEUL,IELVS used w/ ALE or EULER
163c I IPARG(NG) element group data
164c I ELBUF() internal element(material) data used w/ ALE or EULER
165c I TF() ,NPF() Radioss function (x=Time) data
166c I BUFMAT() internal material data
167c IO PARTSAV() output use per part
168c IO DT2T smallest elementary time step
169c O NELTST,ITYPTST element type (property type for spring) which determine DT2T
170c IO STIFN(NUMNOD) nodal stiffness to calcul nodal time step
171c IO EANI() anim output vector
172c I IPARTQ() quad element group data (output)
173c I NEL nb of quad element in this group
174c I IADQ(),FSKY() work arrays for special option of internal force assemlage
175c I ICP flag for constant pressure
176c I IPM(NPROPMI,*) MATERIAL DATA (INTEGER)
177c I BUFVOIS() work table for fluid w/ SPMD
178c I QMV(8,) work table used w/ ALE or EULER
179c I GRESAV,GRTH,IGRTH work table used for TH (time history) output
180c I TABLE new alternative Radioss function(table) data
181c I IGEO geometrical property integer data
182C-----------------------------------------------
183C L o c a l V a r i a b l e s
184C-----------------------------------------------
185 INTEGER I,II,IFLAG,IOFFS,ICPG,ISTRAIN
186C FOR THE USAGE OF EMPTY INPUT
187 INTEGER IBID,IBIDON(1),ipr,enum,SZ_IX
188 INTEGER SZ_BUFVOIS
189C INDEX OF ELEMENT INFORMATION IN GLOBAL TABLE "IC"
190 INTEGER LCO
191C SN OF THE FIRST ELEMENT OF THE GROUP IN GLOBAL STORAGE
192 INTEGER NF1
193C MATERIAL SN, CONNECTIVITY, ID, PROPERTY SN OF ELEMENTS
194 INTEGER MAT(MVSIZ),
195 + NC1(MVSIZ),NC2(MVSIZ),NC3(MVSIZ),NC4(MVSIZ),
196 + NGL(MVSIZ),NGEO(MVSIZ)
197C NUMBER AND INDEXES OF INTEGRATION POINTS
198 INTEGER NPTR,NPTS,IR,IS,ILAY
199C
200 my_real,
201 . DIMENSION(:), POINTER :: eint
202 TYPE(g_bufel_) ,POINTER :: GBUF
203 TYPE(l_bufel_) ,POINTER :: LBUF
204C FOR THE USAGE OF EMPTY INPUT
205 my_real
206 + bid(mvsiz),mbid(mvsiz),ehou(mvsiz)
207C
208C LOCAL DIRECTIONS FOR ORTHOTROPIC MATERIAL USAGE
209C STIFNESS AT INTEGRATION POINT AND OF THE ELEMENT
210C NODAL COORDINATES (t+dt)
211C NODAL VELOCITIES (t+dt/2)
212C DIFFERENCES OF "Y", "Z"
213C SUMMERIES OF "Y"
214C AREA, VOLUME/(THICKNESS .OR. 2*PI), CHARACTERISTIC LENGTH
215C TRANSFORMATION MATRIX [R] FOR CO-ROTATIONAL CASE
216C {X}=[R]{X'} <=> {X'}=T([R]){X}
217 my_real
218 + offg(mvsiz),offs(mvsiz),off(mvsiz),
219 + gama(mvsiz,6),
220 + sti(mvsiz),stim(mvsiz),
221 + y1(mvsiz),y2(mvsiz),y3(mvsiz),y4(mvsiz),
222 + z1(mvsiz),z2(mvsiz),z3(mvsiz),z4(mvsiz),
223 + vy1(mvsiz),vy2(mvsiz),vy3(mvsiz),vy4(mvsiz),
224 + vz1(mvsiz),vz2(mvsiz),vz3(mvsiz),vz4(mvsiz),
225 + y12(mvsiz),y34(mvsiz),y13(mvsiz),y24(mvsiz),
226 + y14(mvsiz),y23(mvsiz),
227 + z12(mvsiz),z34(mvsiz),z13(mvsiz),z24(mvsiz),
228 + z14(mvsiz),z23(mvsiz),
229 + y234(mvsiz),y124(mvsiz),yavg(mvsiz),
230 + aire(mvsiz),volu(mvsiz),deltax(mvsiz),
231 + r11(mvsiz),r12(mvsiz),r13(mvsiz),
232 + r21(mvsiz),r22(mvsiz),r23(mvsiz),
233 + r31(mvsiz),r32(mvsiz),r33(mvsiz)
234C
235C USED WITH ALE OR EULER
236C VISCOUS PRESSURE (SPHERICAL, POSITIVE IF PRESSURE, OUTPUT OF "MMAIN")
237C SHAPE DERIVATIVES (dNi/dY, dNi/dZ) AT CENTER
238C Ni/r AT CENTER
239C SHAPE DERIVATIVES (dNi/dY, dNi/dZ) AT INTEGRATION POINT
240C JACOBIAN MATRIX [J] AT INTEGRATION POINT
241C (W*|J| OR r'* W*|J|) AT INTEGRATION POINT
242C (RHO*VOL-RHO0*VOL0')/RHO == RHO0*(VOL0-VOL0')/RHO
243C PLASTICITY STATE INFORMATION FOR PRESSURE COMPUTATION
244C RATE OF STRAIN FOR STRESS UPDATE
245C COMPONENTS OF ROTATION TENSOR
246C COMPONENTS OF OLD STRESS
247C INTERNAL FORCE IN LOCAL STORAGE
248C ELEMENT AVERAGE PRESSURE
249 my_real
250 + vdy(mvsiz),vdz(mvsiz),vdx(mvsiz),vd2(mvsiz),
251 + muvoid(mvsiz),
252 + vis(mvsiz),
253 + qvis(mvsiz),
254 + ssp(mvsiz),
255 + ssp_eq(mvsiz),
256 + sigy(mvsiz),et(mvsiz),
257c + DEFP(MVSIZ),
258 + r3_free(mvsiz),r4_free(mvsiz),
259 + pyc1(mvsiz),pyc2(mvsiz),pzc1(mvsiz),pzc2(mvsiz),
260 + ay(mvsiz),
261 + py1(mvsiz),py2(mvsiz),py3(mvsiz),py4(mvsiz),
262 + pz1(mvsiz),pz2(mvsiz),pz3(mvsiz),pz4(mvsiz),
263 + rx(mvsiz),ry(mvsiz),rz(mvsiz),
264 + sx(mvsiz),sy(mvsiz),sz(mvsiz),
265 + tx(mvsiz),ty(mvsiz),tz(mvsiz),
266 + airn(mvsiz),voln(mvsiz),
267 + dvol(mvsiz),
268 + nu(mvsiz),e0(mvsiz),c1,fac(mvsiz),
269 + eyy(mvsiz),ezz(mvsiz),exx(mvsiz),
270 + eyz(mvsiz),ezx(mvsiz),exy(mvsiz),
271 + wyy(mvsiz),wzz(mvsiz),wxx(mvsiz),
272 + s1(mvsiz),s2(mvsiz),s3(mvsiz),
273 + s4(mvsiz),s5(mvsiz),s6(mvsiz),
274 + fy1(mvsiz),fz1(mvsiz),fy2(mvsiz),fz2(mvsiz),
275 + fy3(mvsiz),fz3(mvsiz),fy4(mvsiz),fz4(mvsiz),
276 + fay(mvsiz),faz(mvsiz),
277 + pp(mvsiz),dsv(mvsiz)
278C
279 my_real
280 + ay1(mvsiz),ay2(mvsiz),ay3(mvsiz),ay4(mvsiz),yh(mvsiz),
281 + fay1(mvsiz),fay2(mvsiz),fay3(mvsiz),fay4(mvsiz),
282 + det(mvsiz),
283 . byz1(mvsiz),byz2(mvsiz),byz3(mvsiz),byz4(mvsiz),
284 . bzy1(mvsiz),bzy2(mvsiz),bzy3(mvsiz),bzy4(mvsiz),
285 + qn1,qn2,qn3,qn4,nuu(mvsiz)
286C
287 my_real varnl(nel)
288 my_real
289 + wi,ksi,eta
290 my_real
291 + w_gauss(9,9),a_gauss(9,9)
292 DATA w_gauss /
293 1 2. ,0. ,0. ,
294 1 0. ,0. ,0. ,
295 1 0. ,0. ,0. ,
296 2 1. ,1. ,0. ,
297 2 0. ,0. ,0. ,
298 2 0. ,0. ,0. ,
299 3 0.555555555555556,0.888888888888889,0.555555555555556,
300 3 0. ,0. ,0. ,
301 3 0. ,0. ,0. ,
302 4 0.347854845137454,0.652145154862546,0.652145154862546,
303 4 0.347854845137454,0. ,0. ,
304 4 0. ,0. ,0. ,
305 5 0.236926885056189,0.478628670499366,0.568888888888889,
306 5 0.478628670499366,0.236926885056189,0. ,
307 5 0. ,0. ,0. ,
308 6 0.171324492379170,0.360761573048139,0.467913934572691,
309 6 0.467913934572691,0.360761573048139,0.171324492379170,
310 6 0. ,0. ,0. ,
311 7 0.129484966168870,0.279705391489277,0.381830050505119,
312 7 0.417959183673469,0.381830050505119,0.279705391489277,
313 7 0.129484966168870,0. ,0. ,
314 8 0.101228536290376,0.222381034453374,0.313706645877887,
315 8 0.362683783378362,0.362683783378362,0.313706645877887,
316 8 0.222381034453374,0.101228536290376,0. ,
317 9 0.081274388361574,0.180648160694857,0.260610696402935,
318 9 0.312347077040003,0.330239355001260,0.312347077040003,
319 9 0.260610696402935,0.180648160694857,0.081274388361574/
320 DATA a_gauss /
321 1 0. ,0. ,0. ,
322 1 0. ,0. ,0. ,
323 1 0. ,0. ,0. ,
324 2 -.577350269189626,0.577350269189626,0. ,
325 2 0. ,0. ,0. ,
326 2 0. ,0. ,0. ,
327 3 -.774596669241483,0. ,0.774596669241483,
328 3 0. ,0. ,0. ,
329 3 0. ,0. ,0. ,
330 4 -.861136311594053,-.339981043584856,0.339981043584856,
331 4 0.861136311594053,0. ,0. ,
332 4 0. ,0. ,0. ,
333 5 -.906179845938664,-.538469310105683,0. ,
334 5 0.538469310105683,0.906179845938664,0. ,
335 5 0. ,0. ,0. ,
336 6 -.932469514203152,-.661209386466265,-.238619186083197,
337 6 0.238619186083197,0.661209386466265,0.932469514203152,
338 6 0. ,0. ,0. ,
339 7 -.949107912342759,-.741531185599394,-.405845151377397,
340 7 0. ,0.405845151377397,0.741531185599394,
341 7 0.949107912342759,0. ,0. ,
342 8 -.960289856497536,-.796666477413627,-.525532409916329,
343 8 -.183434642495650,0.183434642495650,0.525532409916329,
344 8 0.796666477413627,0.960289856497536,0. ,
345 9 -.968160239507626,-.836031107326636,-.613371432700590,
346 9 -.324253423403809,0. ,0.324253423403809,
347 9 0.613371432700590,0.836031107326636,0.968160239507626/
348C
349C-----------------------------------------------
350C S o u r c e L i n e s
351C-----------------------------------------------
352 sz_ix=numelq+numels+nsvois ! Size of IX array (either IXS+NSVOIS or IXQ)
353 ilay = 1
354 gbuf => elbuf_tab(ng)%GBUF
355 ibid = 0
356 ibidon(1) = 0
357 bid(:) = zero
358
359C
360 DO i=lft,llt
361 ezx(i)=zero
362 exy(i)=zero
363 wyy(i)=zero
364 wzz(i)=zero
365 vdy(i)=zero
366 vdz(i)=zero
367 vdx(i)=zero
368 ENDDO
369C
370 IF (isorth == 0) THEN
371 DO i=lft,llt
372 gama(i,1) = one
373 gama(i,2) = zero
374 gama(i,3) = zero
375 gama(i,4) = zero
376 gama(i,5) = one
377 gama(i,6) = zero
378 ENDDO
379 ELSE
380 DO i=lft,llt
381 gama(i,1) = gbuf%GAMA(i )
382 gama(i,2) = gbuf%GAMA(i + nel)
383 gama(i,3) = gbuf%GAMA(i + 2*nel)
384 gama(i,4) = gbuf%GAMA(i + 3*nel)
385 gama(i,5) = gbuf%GAMA(i + 4*nel)
386 gama(i,6) = gbuf%GAMA(i + 5*nel)
387 ENDDO
388 ENDIF
389 isorthg = 0
390 istrain = iparg(44,ng)
391C
392C GATHER NODAL INFORMATION &
393C COMPUTE INTRINSIC ROTATION FOR CO-ROTATIONAL CASE &
394C PROJECT NODAL COORDINATES AND VELOCITES INTO CO-ROTATIONAL SYSTEM
395 lco = 1 + 7*nft
396 nf1 = 1 + nft
397 IF(jcvt==0) THEN
398 CALL q4coor2(
399 1 x, ic(lco), y1, y2,
400 2 y3, y4, z1, z2,
401 3 z3, z4, nc1, nc2,
402 4 nc3, nc4, ngl, mat,
403 5 ngeo, vd2, vis, v,
404 6 vy1, vy2, vy3, vy4,
405 7 vz1, vz2, vz3, vz4,
406 8 yavg, ay, exx, nel,
407 9 jhbe)
408 ELSE
409C----EXX,YAVG ,AY are calculated in global system anyway, more efficient eo be done here
410C----for the moment constant in element
411 CALL q4rcoor2(
412 1 x, ic(lco), y1, y2,
413 2 y3, y4, z1, z2,
414 3 z3, z4, nc1, nc2,
415 4 nc3, nc4, ngl, mat,
416 5 ngeo, vd2, r11, r12,
417 6 r13, r21, r22, r23,
418 7 r31, r32, r33, gama,
419 8 y234, y124, vis, v,
420 9 vy1, vy2, vy3, vy4,
421 a vz1, vz2, vz3, vz4,
422 b yavg, ay, exx, nel,
423 c isorth)
424 ENDIF
425C COMPUTE FAC(*) FOR PRESSURE ACCORDING TO PLASTICITY STATE
426C---- now assumed strain is used, no effect for Isolid17
427 DO i=lft,llt
428 nu(i)=min(half,pm(21,mat(i)))
429 c1 =pm(32,mat(i))
430 e0(i) =three*(one-two*nu(i))*c1
431 ENDDO
432 IF(icp==2) THEN
433 CALL s8zsigp3(lft ,llt ,gbuf%SIG,e0 ,gbuf%PLA,
434 2 fac ,gbuf%G_PLA,nel )
435 DO i=lft,llt
436 nuu(i)=nu(i)+(half-nu(i))*fac(i)
437 ENDDO
438 ELSEIF(icp==1) THEN
439 DO i=lft,llt
440 nuu(i)=half
441 ENDDO
442 ELSE
443 DO i=lft,llt
444 nuu(i)=zero
445 ENDDO
446 ENDIF
447C COMPUTE AREA & VOLUME
448 CALL qvolu2(
449 1 gbuf%OFF,aire, volu, ngl,
450 2 y1, y2, y3, y4,
451 3 z1, z2, z3, z4,
452 4 y234, y124, nel, jmult,
453 5 jcvt)
454C
455C COMPUTE CHARACTERISTIC LENGTH OF EACH ELEMENT FOR DETERMINING TIME STEP
456 CALL qdlen2(y1,y2,y3,y4,z1,z2,z3,z4,aire,deltax,iparg(63,ng))
457C
458C COMPUTE SHAPE DERIVATIVES AT ELEMENT CENTER
459C
460 CALL q4deric2(
461 1 y1, y2, y3, y4,
462 2 z1, z2, z3, z4,
463 3 y12, y34, y13, y24,
464 4 y14, y23, z12, z34,
465 5 z13, z24, z14, z23,
466 6 pyc1, pyc2, pzc1, pzc2,
467 7 aire, volu, yavg, rx,
468 8 ry, rz, sx, sy,
469 9 sz, nel, jhbe)
470C
471C -----ICPG could be cleaned after---
472 icpg = 0
473C IF(ICPG==2) ICPG = 1
474C
475C COMPUTE SHEAR & VOLUMETRIC STRAIN RATE AT ELEMENT CENTER
476C WITH B(t+dt), V(t+dt/2) & CORRECTION
477 CALL q4defoc2(
478 1 vy1, vy2, vy3, vy4,
479 2 vz1, vz2, vz3, vz4,
480 3 pyc1, pyc2, pzc1, pzc2,
481 4 aire, eyz, exx, dsv,
482 5 icpg, nel, jcvt)
483C
484C
485 ioffs = 0
486 DO i=lft,llt
487 offs(i) = ep20
488 offg(i) = gbuf%OFF(i)
489 ENDDO
490C
491C
492C INITIALIZATION BEFORE INTEGRATION LOOP
493 CALL q4zero2(
494 + fy1, fz1, fy2, fz2, fy3, fz3, fy4, fz4,
495 + fay, faz, fay1, fay2, fay3, fay4,
496 + gbuf%SIG,gbuf%EINT,gbuf%RHO,gbuf%QVIS,gbuf%PLA,
497 + gbuf%EPSD,stim,pp,gbuf%G_PLA,gbuf%G_EPSD,nel)
498C
499C ENTER THE INTEGRATION POINTS LOOP -->
500 nptr = 2
501 npts = 2
502 DO 100 ir=1,nptr
503 DO 200 is=1,npts
504 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
505C
506C INITIALIZE WEIGHTING FACTORS
507 ksi = a_gauss(ir,nptr)
508 eta = a_gauss(is,npts)
509 wi = w_gauss(ir,nptr)*w_gauss(is,npts)
510C
511C INITIALIZE INDEX OF ELEMENT DATA ARRAY
512C
513C COMPUTE JACOBIAN & SHAPE DERIVATIVES AT INTEGRATION POINT
514 CALL q4deri2(
515 1 offg, off, ksi, eta,
516 2 wi, yavg, y12, y34,
517 3 y13, y24, y14, y23,
518 4 z12, z34, z13, z24,
519 5 z14, z23, py1, py2,
520 6 py3, py4, pz1, pz2,
521 7 pz3, pz4, pyc1, pyc2,
522 8 pzc1, pzc2, byz1, byz2,
523 9 byz3, byz4, bzy1, bzy2,
524 a bzy3, bzy4, airn, voln,
525 b nuu, nel, jhbe)
526C
527C COMPUTE RATE OF STRAIN AT INTEGRATION POINT WITH B(t+dt),V(t+dt/2) & CORRECTION
528C MODIFY INITIAL VOLUME & INTERNAL ENERGY DENSITY
529 CALL q4defo2(
530 1 py1, py2, py3, py4,
531 2 pz1, pz2, pz3, pz4,
532 3 byz1, byz2, byz3, byz4,
533 4 bzy1, bzy2, bzy3, bzy4,
534 5 vy1, vy2, vy3, vy4,
535 6 vz1, vz2, vz3, vz4,
536 7 eyz, eyy, ezz, exx,
537 8 wxx, r22, r23, ay,
538 9 off, gbuf%OFF, lbuf%VOL, lbuf%EINT,
539 a dsv, icpg, fac, nel,
540 b jcvt)
541C
542C
543C COMPUTE ADDITIONAL VOLUME CHANGE & MODIFY CURRENT DENSITY
544C DVOL = (RHO*VOL-RHO0*VOL0')/RHO = RHO0*(VOL0-VOL0')/RHO
545C RHO' = RHO0*VOL0'/VOL
546 CALL qlagr2(
547 1 pm, lbuf%VOL, lbuf%RHO, lbuf%EINT,
548 2 voln, dvol, mat, nel)
549C
550C COPY OLD STRESS AND DO NECESSARY TREATMENT
551 CALL qrota2(
552 1 lbuf%SIG,s1, s2, s3,
553 2 s4, s5, s6, wxx,
554 3 nel, jcvt)
555C
556C UPDATE STRESS
557 CALL mmain(timers, output,
558 1 elbuf_tab, ng, pm, geo,
559 2 ale_connect, ic, iparg,
560 3 v, tf, npf, bufmat,
561 4 sti, x, dt2t, neltst,
562 5 ityptst, offset, nel, w,
563 6 off, ngeo, mat, ngl,
564 7 voln, vd2, dvol, deltax,
565 8 vis, qvis, ssp, s1,
566 9 s2, s3, s4, s5,
567 a s6, eyy, ezz, exx,
568 b eyz, ezx, exy, wyy,
569 c wzz, wxx, rx, ry,
570 d rz, sx, sy, sz,
571 e vdy, vdz, vdx, muvoid,
572 f ssp_eq, aire, sigy, et,
573 g bufvois, lbuf%PLA, r3_free, r4_free,
574 h eyy, ezz, exx, eyz,
575 i ezx, exy, wyy, wzz,
576 j wxx, ipm, gama, bid,
577 k bid, bid, bid, bid,
578 l bid, bid, istrain, bid,
579 m bid, ibidon(1), ilay, mbid,
580 n mbid, ir, is, 1,
581 o table, bid, bid, bid,
582 p bid, iparg(1,ng), igeo, bid,
583 q itask, nloc_dmg, varnl, mat_elem,
584 r h3d_strain, jplasol, jsph, sz_bufvois,
585 s snpc, stf, sbufmat, glob_therm,
586 t svis, sz_ix, iresp,
587 * n2d, th_strain, ngroup, tt,
588 . dt1, ntable, numelq, nummat,
589 . numgeo, numnod, numels,
590 . idel7nok, idtmin, maxfunc,
591 . imon_mat, userl_avail, impl_s,
592 . idyna, dt , bid ,sensors)
593C
594C COMPUTE NODAL INTERNAL FORCE
595 CALL q4fint2(
596 1 lbuf%SIG,ay, fay, faz,
597 2 py1, py2, py3, py4,
598 3 pz1, pz2, pz3, pz4,
599 4 byz1, byz2, byz3, byz4,
600 5 bzy1, bzy2, bzy3, bzy4,
601 6 fy1, fz1, fy2, fz2,
602 7 fy3, fz3, fy4, fz4,
603 8 r22, r23, r32, r33,
604 9 airn, voln, qvis, icpg,
605 a nel, jhbe, jcvt, svis)
606C
607C COMPUTE ELEMENT AVERAGE & SUMMERY DATA
608 CALL s8efmoy3(
609 1 lbuf%SIG, voln, qvis, pp,
610 2 lbuf%EINT, lbuf%RHO, lbuf%QVIS, lbuf%PLA,
611 3 lbuf%EPSD, gbuf%EPSD, gbuf%SIG, gbuf%EINT,
612 4 gbuf%RHO, gbuf%QVIS, gbuf%PLA, volu,
613 5 sti, stim, icpg, off,
614 6 lbuf%VOL, gbuf%VOL, gbuf%G_PLA, gbuf%G_EPSD,
615 7 lbuf%EINTTH,gbuf%EINTTH,iexpan, nel,
616 8 bid, bid,svis,glob_therm%NODADT_THERM,
617 9 gbuf%WPLA, lbuf%WPLA, gbuf%G_WPLA)
618C
619 DO i=lft,llt
620 offg(i)=min(offg(i),off(i))
621 IF (lbuf%OFF(i) > one .AND. gbuf%OFF(i) == one) THEN
622 offs(i) = min(lbuf%OFF(i),offs(i))
623 ioffs = 1
624 ENDIF
625 ENDDO
626C
627200 CONTINUE
628100 CONTINUE
629C EXIT THE INTEGRATION POINTS LOOP <--
630C
631 IF(ioffs==1)THEN
632 DO i=lft,llt
633 IF(offs(i)<=two) THEN
634 gbuf%OFF(i) = offs(i)
635 ENDIF
636 ENDDO
637 DO ir=1,nptr
638 DO is=1,npts
639 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
640 DO i=lft,llt
641 IF (gbuf%OFF(i) > one) lbuf%OFF(i) = gbuf%OFF(i)
642 ENDDO
643 ENDDO
644 ENDDO
645 ENDIF
646 CALL smallb3(
647 1 gbuf%OFF,offg, nel, ismstr)
648C
649C ADD NODAL INTERNAL FORCE FOR CONSTANT PRESSURE AND CONSTANT SHEAR STRAIN
650 CALL q4fintc2(
651 1 pyc1, pyc2, pzc1, pzc2,
652 2 ay, fay, fy1, fz1,
653 3 fy2, fz2, fy3, fz3,
654 4 fy4, fz4, aire, volu,
655 5 gbuf%SIG,pp, icpg, nel,
656 6 jhbe)
657C
658 IF(n2d==1.AND.jhbe==17) THEN
659 CALL q4vis2(
660 1 pm, gbuf%OFF,gbuf%RHO,y1,
661 2 y2, y3, y4, z1,
662 3 z2, z3, z4, vy1,
663 4 vy2, vy3, vy4, vz1,
664 5 vz2, vz3, vz4, py1,
665 6 py2, pz1, pz2, fy1,
666 7 fy2, fy3, fy4, fz1,
667 8 fz2, fz3, fz4, aire,
668 9 ssp, nel)
669 ENDIF
670C OUTPUT ELEMENT INFORMATION TO BALANCE TABLE
671 iflag=mod(ncycle,ncpri)
672 IF(ioutprt>0)THEN
673c
674 IF (mtn == 11) THEN
675 eint => elbuf_tab(ng)%GBUF%EINS(1:nel)
676 ELSE
677 eint => elbuf_tab(ng)%GBUF%EINT(1:nel)
678 ENDIF
679 CALL qbilan(
680 1 partsav, gbuf%OFF, eint, gbuf%RHO,
681 2 gbuf%RK, gbuf%VOL, vy1, vy2,
682 3 vy3, vy4, vz1, vz2,
683 4 vz3, vz4, volu, ipartq,
684 5 ehou, r22, r23, r32,
685 6 r33, gresav, grth, igrth,
686 7 ibidon(1), gbuf%EINTTH,itask, nel,
687 8 jtur, jcvt, igre, sensors,
688 9 gbuf%G_WPLA,gbuf%WPLA)
689 ENDIF
690C
691C PROJECT NODAL INTERNAL FORCE INTO GLOBAL SYSTEM FOR CO-ROTATIONAL CASE
692 IF(jcvt/=0) THEN
693 CALL q4rrota2(
694 1 r22, r32, r23, r33,
695 2 fy1, fy2, fy3, fy4,
696 3 fz1, fz2, fz3, fz4,
697 4 nel)
698 ENDIF
699C
700C ADD TERMS OF INTERNAL FORCE FOR AXISYMMETRIC CASE ONLY
701 IF(n2d==1.AND.jhbe==17) THEN
702 DO i=lft,llt
703 fy1(i) = fy1(i) + fay(i)
704 fy2(i) = fy2(i) + fay(i)
705 fy3(i) = fy3(i) + fay(i)
706 fy4(i) = fy4(i) + fay(i)
707 fz1(i) = fz1(i) + faz(i)
708 fz2(i) = fz2(i) + faz(i)
709 fz3(i) = fz3(i) + faz(i)
710 fz4(i) = fz4(i) + faz(i)
711 ENDDO
712 ENDIF
713C
714C COMPUTE MASS STIFFNESS
715 IF(n2d==1.AND.jhbe==17) THEN
716 IF(iparit==0) THEN
717 CALL qmass2(
718 1 gbuf%OFF,gbuf%RHO,ms, aire,
719 2 nc1, nc2, nc3, nc4,
720 3 nel)
721 ELSE
722 CALL qmass2p(
723 1 gbuf%OFF,gbuf%RHO,aire, fsky,
724 2 fsky, iadq, nel, nft)
725 ENDIF
726 ELSE
727 IF(iparit==0) THEN
728 CALL qmass2(
729 1 gbuf%OFF,gbuf%RHO,ms, volu,
730 2 nc1, nc2, nc3, nc4,
731 3 nel)
732 ELSE
733 CALL qmass2p(
734 1 gbuf%OFF,gbuf%RHO,volu, fsky,
735 2 fsky, iadq, nel, nft)
736 ENDIF
737 ENDIF
738C
739C--------------------------
740C UPDATE OF MASSES : ALE physical masses
741C----------------------------
742 IF (jale+jeul > 0 )THEN
743 IF (iparit == 0)THEN
744 CALL qmassreal2(
745 1 gbuf%OFF,gbuf%RHO,ms_2d, volu,
746 2 nc1, nc2, nc3, nc4,
747 3 nel)
748 ELSE
749 CALL qmassreal2p(
750 1 gbuf%OFF,gbuf%RHO,volu, fskym,
751 2 iadq, nel, nft)
752 ENDIF
753 ENDIF
754C
755C ASSEMBLE NODAL INTERNAL FORCE INTO GLOBAL STORAGE
756 IF(iparit==0) THEN
757 CALL q4cumu2(
758 1 a, stifn, nc1, nc2,
759 2 nc3, nc4, fy1, fz1,
760 3 fy2, fz2, fy3, fz3,
761 4 fy4, fz4, stim, nel)
762 ELSE
763 CALL q4cumu2p(
764 1 fsky, fsky, iadq, fy1,
765 2 fz1, fy2, fz2, fy3,
766 3 fz3, fy4, fz4, stim,
767 4 nel, nft)
768 ENDIF
769C-----------
770 RETURN
771 END
q4coor2
subroutine q4coor2(x, ncp, y1, y2, y3, y4, z1, z2, z3, z4, nc1, nc2, nc3, nc4, ngl, mat, ngeo, vd2, vis, v, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, yavg, ay, exx, nel, jhbe)
Definition q4coor2.F:38
min
#define min(a, b)
Definition macros.h:20
mmain
subroutine mmain(pm, elbuf_str, ix, nix, x, geo, iparg, nel, skew, bufmat, ipart, ipartel, nummat, matparam, imat, ipm, ngl, pid, npf, tf, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, rx, ry, rz, sx, sy, sz, gama, voln, dvol, s1, s2, s3, s4, s5, s6, dxx, dyy, dzz, d4, d5, d6, wxx, wyy, wzz)
Definition mmain.F:43
ale_connectivity_mod
Definition ale_connectivity_mod.F:223
dt_mod
Definition dt_mod.F:79
nlocal_reg_mod
Definition nlocal_reg_mod.F:117
table_mod
Definition table_mod.F:112
q4cumu2
subroutine q4cumu2(e, stifn, nc1, nc2, nc3, nc4, fy1, fz1, fy2, fz2, fy3, fz3, fy4, fz4, stim, nel)
Definition q4cumu2.F:33
q4cumu2p
subroutine q4cumu2p(fsky, fskyv, iadq, fy1, fz1, fy2, fz2, fy3, fz3, fy4, fz4, stim, nel, nft)
Definition q4cumu2p.F:33
q4defo2
subroutine q4defo2(py1, py2, py3, py4, pz1, pz2, pz3, pz4, byz1, byz2, byz3, byz4, bzy1, bzy2, bzy3, bzy4, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, eyz, eyy, ezz, exx, wxx, r22, r23, ay, off, offg, volo, eint, dsv, icp, fac, nel, jcvt)
Definition q4defo2.F:40
q4defoc2
subroutine q4defoc2(vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, pyc1, pyc2, pzc1, pzc2, aire, eyz, exx, dsv, icp, nel, jcvt)
Definition q4defoc2.F:34
q4deric2
subroutine q4deric2(y1, y2, y3, y4, z1, z2, z3, z4, y12, y34, y13, y24, y14, y23, z12, z34, z13, z24, z14, z23, pyc1, pyc2, pzc1, pzc2, aire, volu, yavg, rx, ry, rz, sx, sy, sz, nel, jhbe)
Definition q4deric2.F:39
q4fint2
subroutine q4fint2(sig, ay, fay, faz, py1, py2, py3, py4, pz1, pz2, pz3, pz4, byz1, byz2, byz3, byz4, bzy1, bzy2, bzy3, bzy4, fy1, fz1, fy2, fz2, fy3, fz3, fy4, fz4, r22, r23, r32, r33, air, vol, qvis, icp, nel, jhbe, jcvt, svis)
Definition q4fint2.F:39
q4fintc2
subroutine q4fintc2(pyc1, pyc2, pzc1, pzc2, ay, fay, fy1, fz1, fy2, fz2, fy3, fz3, fy4, fz4, aire, volu, sigm, pp, icp, nel, jhbe)
Definition q4fintc2.F:35
q4forc2
subroutine q4forc2(timers, output, pm, geo, ic, x, a, v, ms, w, flux, flu1, veul, fv, ale_connect, iparg, nloc_dmg, elbuf_tab, tf, npf, bufmat, partsav, dt2t, neltst, ityptst, stifn, offset, eani, ipartq, nel, iadq, fsky, icp, ng, ipm, bufvois, qmv, gresav, grth, igrth, table, igeo, itask, iexpan, ms_2d, fskym, ioutprt, mat_elem, h3d_strain, sz_bufvois, snpc, stf, sbufmat, svis, nsvois, idtmins, iresp, tt, dt1, idel7ng, idel7nok, idtmin, maxfunc, imon_mat, userl_avail, impl_s, idyna, dt, glob_therm, sensors)
Definition q4forc2.F:81
q4rrota2
subroutine q4rrota2(r22, r23, r32, r33, y1, y2, y3, y4, z1, z2, z3, z4, nel)
Definition q4rrota2.F:34
q4vis2
subroutine q4vis2(pm, off, rho, y1, y2, y3, y4, z1, z2, z3, z4, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, py1, py2, pz1, pz2, t11, t12, t13, t14, t21, t22, t23, t24, area, cxx, nel)
Definition q4vis2.F:38
q4zero2
subroutine q4zero2(fy1, fz1, fy2, fz2, fy3, fz3, fy4, fz4, fay, faz, fay1, fay2, fay3, fay4, sigm, eintm, rhom, qm, eplasm, epsdm, stin, pp, g_pla, g_epsd, nel)
Definition q4zero2.F:33
qbilan
subroutine qbilan(partsav, off, eint, rho, rk, vol, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, vnew, ipartq, ehou, r22, r23, r32, r33, gresav, grth, igrth, iexpan, eintth, itask, nel, jtur, jcvt, igre, sensors, g_wpla, wpla)
Definition qbilan.F:45
qlagr2
subroutine qlagr2(pm, vol0, rho, eint, voln, dvol, mat, nel)
Definition qlagr2.F:32
qmass2
subroutine qmass2(off, rho, ms, area, nc1, nc2, nc3, nc4, nel)
Definition qmass2.F:37
qmass2p
subroutine qmass2p(off, rho, area, fsky, fskyv, iadq, nel, nft)
Definition qmass2p.F:33
qmassreal2
subroutine qmassreal2(off, rho, ms_2d, vnew, nc1, nc2, nc3, nc4, nel)
Definition qmassreal2.F:41
qmassreal2p
subroutine qmassreal2p(off, rho, vnew, fskym, iadq, nel, nft)
Definition qmassreal2p.F:33
qrota2
subroutine qrota2(sig, s1, s2, s3, s4, s5, s6, wyz, nel, jcvt)
Definition qrota2.F:34
qvolu2
subroutine qvolu2(off, aire, volu, ngl, y1, y2, y3, y4, z1, z2, z3, z4, y234, y124, nel, jmult, jcvt)
Definition qvolu2.F:43
s8efmoy3
subroutine s8efmoy3(sigor, vol, qvis, pp, eint, rho, q, defp, epsd, epsdm, sigm, eintm, rhom, qm, defpm, volg, sti, stin, icp, off, vol0, vol0g, g_pla, g_epsd, eintth, eintthm, iexpan, nel, conde, conden, svis, nodadt_therm, g_wpla, l_wpla, g_wpla_flag)
Definition s8efmoy3.F:40
s8zsigp3
subroutine s8zsigp3(lft, llt, sig, e0, defp, fac, g_pla, nel)
Definition s8zsigp3.F:37
smallb3
subroutine smallb3(offg, off, nel, ismstr)
Definition smallb3.F:44
q4rcoor2
subroutine q4rcoor2(x, ixq, ngl, mxt, pid, ix1, ix2, ix3, ix4, y1, y2, y3, y4, z1, z2, z3, z4, yavg, y234, y124, sy, sz, ty, tz, e1y, e1z, e2y, e2z)
Definition q4coor2.F:33
q4deri2
subroutine q4deri2(vol, ksi, eta, wi, y12, y34, y13, y24, y14, y23, z12, z34, z13, z24, z14, z23, y1, y2, y3, y4, yavg, ihbe, ngl)
Definition q4deri2.F:36
qdlen2
subroutine qdlen2(iparg, aire, deltax, y1, y2, y3, y4, z1, z2, z3, z4)
Definition qdlen2.F:39