OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
initwg_solid.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!|| initwg_solid ../starter/source/spmd/domain_decomposition/initwg_solid.F
25!||--- called by ------------------------------------------------------
26!|| initwg ../starter/source/spmd/domain_decomposition/initwg.F
27!||--- calls -----------------------------------------------------
28!|| bidon ../starter/source/system/machine.F
29!||--- uses -----------------------------------------------------
30!|| ddweights_mod ../starter/share/modules1/ddweights_mod.F
31!|| mid_pid_mod ../starter/share/modules1/mid_pid_mod.F
32!||====================================================================
33 SUBROUTINE initwg_solid(WD,PM,GEO,IXS,IGEO,ISOLNOD,
34 . NUMELS,IPM, SIZE_IRUP,
35 . NUMMAT,NUMGEO,
36 . POIN_PART_SOL,MID_PID_SOL,IPARTS,BUFMAT,
37 . MID_OLD,PID_OLD,MLN_OLD,RECHERCHE,ISOL_OLD,
38 . TELT_PRO,TABMP_L,NPART,MAT_PARAM)
39C-----------------------------------------------
40C M o d u l e s
41C-----------------------------------------------
42 USE ddweights_mod
43 USE mid_pid_mod
44 USE matparam_def_mod
45 use element_mod , only : nixs
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50C-----------------------------------------------
51C C o m m o n B l o c k s
52C-----------------------------------------------
53#include "param_c.inc"
54#include "com01_c.inc"
55#include "tablen_c.inc"
56#include "ddspmd_c.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER
61 . NUMELS,
62 . NUMMAT,NUMGEO,
63 . IXS(NIXS,*),IGEO(NPROPGI,NUMGEO),ISOLNOD(*),
64 . IPM(NPROPMI,*),TABMP_L,NPART
65 INTEGER, INTENT(IN) :: SIZE_IRUP
66C REAL OR REAL*8
67 my_real
68 . PM(NPROPM,*), GEO(NPROPG,*),BUFMAT(*)
69 REAL WD(*)
70 INTEGER MID_OLD,PID_OLD,MLN_OLD,RECHERCHE,ISOL_OLD
71 my_real TELT_PRO
72
73 INTEGER, DIMENSION(2,NPART,*), INTENT(IN) :: POIN_PART_SOL
74 INTEGER, DIMENSION(*), INTENT(IN) :: IPARTS
75 TYPE(mid_pid_type), DIMENSION(NUMMAT,*), INTENT(INOUT) :: MID_PID_SOL
76 TYPE(matparam_struct_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM
77C-----------------------------------------------
78 INTEGER OFF, NPN, MID, PID, JHBE, IGT, MLN,
79 . istrain, ithk, ihbe, ipla, issn, mtn, i, j, k,l,
80 . nfunc,mpt,npts,nptt,nptr,nptot,iflag,jsrot,ivisc,
81 . i_mid,i_pid,i_mid_old,i_pid_old,puid,muid,
82 . elm_typ,elm_typ_old,ilaw,ilaw_old,test_mat,
83 . i_pro,isol2,muid_old,puid_old,
84 . test,nfunc1,nfunc2,nfail,irup2,
85 . isol,indi,iad,indi2,mult
86 INTEGER :: INDI3,ADD_OPTION,INDI_OPT_1,INDI_OPT_2
87 INTEGER :: IRUP_TAB(SIZE_IRUP)
88 my_real :: OPT_1,OPT_2
89
90 REAL
91 . wtype(9),fwihbe,fac8,
92 . tabmat(3),tabx(3),timmat,npt,telt,poids,w,
93 . batozmult,tmat,trup,tmatadd,wd_local
94 INTEGER :: FLAG_NICE_NEWTON,FLAG_GURSON,FLAG_NON_LOCAL
95 INTEGER :: SPECIAL_OPTION,SPE_I_1,SPE_I_2,SPE_I_3
96 my_real :: INVTREF,MULT_SPE
97 INTEGER :: INDI4,POIN_PID,POIN_MID,POIN_PART,COST_CHECK,POIN_ELM_TYP
98 my_real :: INVTELT_PRO
99 my_real :: cc,a,b,a1,a2
100 ! thick shell element cost :
101 INTEGER :: OVERCOST_ELM ,ICPR,NUMBER_LAYER
102 INTEGER :: NLAY,COMPOSITE_MID,COMPOSITE_MLN
103 LOGICAL :: COMPOSITE_OPTION
104
105 LOGICAL :: ISMSTR_COST
106 INTEGER :: ISMSTR,ISMSTR_L,ISM0,ICP0
107 my_real :: add_over_cost,visc_prony
108
109 DATA wtype /1.6 ,1. ,1. ,.9 ,1.1 ,1.4 ,0.65 ,.9 ,2.0/
110C-----------------------------------------------
111 overcost_elm = 0
112 telt = 0
113 jsrot = 0
114 nfunc = 0
115 CALL bidon()
116! DD_OPTIMIZATION = 0 --> default case, DD optimized for Broadwell processor
117! DD_OPTIMIZATION = 1 --> DD optimized for Skylake processor
118! DD_OPTIMIZATION = 2 --> DD optimized for Sandy Bridge processor
119! DD_OPTIMIZATION = 3 --> default case for ARM processor, DD optimized for ThunderX2 processor (ARM)
120 IF(dd_optimization==1) THEN
121! Skylake processor
122#include "weights_p4linux964_spmd_avx512.inc"
123 ELSEIF(dd_optimization==2) THEN
124! Sandy Bridge processor
125#include "weights_p4linux964_spmd_sse3.inc"
126 ELSEIF(dd_optimization==3) THEN
127! ThunderX2 processor (ARMV8.0)
128#include "weights_p4linuxa964_spmd.inc"
129 ELSE
130! DEFAULT CASE
131#if ARCH_CPU
132! ThunderX2 processor (ARMV8.0)
133#include "weights_p4linuxa964_spmd.inc"
134#elif 1
135! Broadwell processor
136#include "weights_p4linux964_spmd.inc"
137#endif
138 ENDIF
139 invtref = one/tpsref
140 DO i = 1, numels
141C -------------------------------
142C Element Property initialization
143C -------------------------------
144 npn = 1
145! JHBE=IHBE_D ! IHBE_D is used for shell
146 jhbe = 1 ! not be important for PID=0
147 mid= ixs(1,i)
148 pid= ixs(10,i)
149
150 mln = nint(pm(19,abs(mid)))
151 isol = isolnod(i)
152 wd_local = wd(i)
153 ! -----------------
154 IF(recherche==1) THEN
155 mid = mid_old
156 pid = pid_old
157 mln = mln_old
158 isol = isol_old
159 wd_local = zero
160 ENDIF
161 ! -----------------
162 IF(isol==8) THEN
163 indi3 = 3
164 ELSEIF(isol==10) THEN
165 indi3 = 4
166 ELSEIF(isol==16) THEN
167 indi3 = 5
168 ELSEIF(isol==20) THEN
169 indi3 = 6
170 ELSEIF(isol==6) THEN
171 indi3 = 7
172 ELSEIF(isol==4) THEN
173 indi3 = 8
174 ELSE
175 indi3 = 9
176 ENDIF
177 ! -----------------
178 IF (pid/=0) THEN
179 jhbe = igeo(10,pid)
180 igt = igeo(11,pid)
181 npn = igeo(4,pid)
182 jsrot = igeo(20,pid)
183 ENDIF
184 nfail = mat_param(abs(mid))%NFAIL
185 irup_tab(1:nfail) = 0
186 IF(nfail/=0) THEN ! up to 6 failure models per material
187 DO j=1,nfail
188 irup_tab(j) = mat_param(abs(mid))%FAIL(j)%IRUPT
189 ENDDO
190 ENDIF
191 tmat = 0.
192 trup = 0.
193 tmatadd = 0.
194 visc_prony = 0.
195 opt_1 = zero
196 opt_2 = zero
197 add_option = 0
198 mult = 0
199 flag_non_local = 0
200 special_option = 0
201 spe_i_1 = 1
202 spe_i_2 = 1
203 ! -------------
204 ! check if composite material is used
205 composite_option = .false.
206 IF (igeo(30,pid)>0 .AND. igeo(11,pid)==22) THEN
207 composite_option = .true.
208 ENDIF
209 ! -------------
210
211 ismstr = igeo(5,pid) ! get the value of ismstr
212
213 ismstr_cost = .false.
214 add_over_cost = zero
215 IF((mln<28).OR.(mln==49).OR.(mln==59)) THEN
216 irup2 = 1
217 ELSE
218 irup2 = 2
219 ENDIF
220
221 ismstr_l = ismstr
222 IF(ismstr<1) THEN
223 ism0 = mat_param(abs(mid))%SMSTR
224 icp0 = mat_param(abs(mid))%STRAIN_FORMULATION
225 IF (icp0 ==2.AND.jhbe/=16) THEN
226 IF (ism0==1) THEN
227 ismstr_l = 11
228 ELSE
229 ismstr_l = 10
230 END IF
231 ELSE
232 IF (ism0==1) THEN
233 ismstr_l = 1
234 ELSE
235 ismstr_l = 2
236 END IF
237 END IF
238 IF (mln == 1.AND.jhbe/=16) ismstr_l = 12
239 ENDIF
240
241
242 IF ( mln==1.OR.mln==38.OR.
243 . mln==90.OR.mln==92.OR.mln==94 ) THEN
244 IF (ismstr_l==10.OR.ismstr_l==12) THEN
245 ismstr_cost = .true.
246 indi = 2
247 ELSE
248 indi = 1
249 ENDIF
250 IF (mat_param(abs(mid))%IVISC > 0) THEN
251 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
252 ENDIF
253 ELSEIF (mln==2) THEN
254 cc = pm(43,mid)
255 IF (cc/=0) THEN
256 indi = 2
257 ELSE
258 indi = 1
259 ENDIF
260 IF (mat_param(abs(mid))%IVISC > 0) THEN
261 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
262 ENDIF
263 ! ----------------
264 ! law 25 : 2 sub-option
265 ELSEIF(mln==25) THEN
266 IF (mat_param(abs(mid))%iparam(1)==1) THEN
267 indi = 2
268 ELSE
269 indi = 1
270 ENDIF
271 IF (mat_param(abs(mid))%IVISC > 0) THEN
272 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
273 ENDIF
274 ! ----------------
275 ELSEIF (mln == 36)THEN
276 nfunc = max(ipm(10,mid) - 3,1)
277 IF (nfunc<=2) THEN
278 indi = 1
279 ELSEIF (nfunc>2.AND.nfunc<=7) THEN
280 indi = 2
281 ELSEIF (nfunc>7) THEN
282 indi = 3
283 ENDIF
284 IF (mat_param(abs(mid))%IVISC > 0) THEN
285 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
286 ENDIF
287 ELSEIF (mln==33) THEN
288 nfunc1 = ipm(11,mid)
289 nfunc2 = ipm(12,mid)
290 IF((nfunc1/=0).OR.(nfunc2/=0)) THEN
291 indi = 2
292 ELSE
293 indi = 1
294 ENDIF
295 IF (mat_param(abs(mid))%IVISC > 0) THEN
296 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
297 ENDIF
298 ELSEIF((mln==42).OR.(mln==62).OR.(mln==69)) THEN
299! check the NPRONY model
300 nfunc = 0
301 IF (mln==42) nfunc = mat_param(abs(mid))%IPARAM(2)
302 IF (mln==62) THEN
303 iad = ipm(7,abs(mid))-1
304 nfunc = nint(bufmat(iad+3))
305 END IF
306 IF (nfunc==0) THEN
307 indi = 1
308 ivisc = mat_param(abs(mid))%IVISC
309 IF (ivisc == 1 .or. ivisc == 2) THEN
310 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
311 ENDIF
312 ELSEIF(nfunc==1) THEN
313 indi = 2
314 ELSEIF(nfunc==2) THEN
315 indi = 3
316 ELSEIF(nfunc>2) THEN
317 indi = 3
318 mult = nfunc - 2
319 indi2 = 2
320 ENDIF
321 ELSEIF((mln==82)) THEN
322 iad=ipm(7,abs(mid))-1
323 nfunc=nint(bufmat(iad+1))
324 IF(nfunc<=1) THEN
325 indi = 1
326 IF (mat_param(abs(mid))%IVISC > 0) THEN
327 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
328 ENDIF
329 ELSEIF(nfunc==2) THEN
330 indi = 2
331 ELSEIF(nfunc==3) THEN
332 indi = 3
333 ELSEIF(nfunc>3) THEN
334 indi = 3
335 mult = nfunc - 3
336 indi2 = 2
337 ENDIF
338 ELSEIF(mln==100) THEN
339 ! SPECIAL TREATMENT :
340 ! LAW 100 : material cost = base cost + viscosity cost + N * network cost
341 ! (optional) (optional)
342 ! INDI INDI_OPT_1 INDI_OPT_2
343 indi=1
344 iad=ipm(7,abs(mid))-1
345
346 add_option = 0
347 opt_1 = zero
348 indi_opt_1 = 2
349 opt_2 = zero
350 indi_opt_2 = 2
351
352 ! viscosity flag
353 IF(nint(bufmat(iad+5))>0) THEN
354 opt_1 = one
355 indi_opt_1 = 2
356 add_option = 1
357 ENDIF
358 ! network flag
359 IF(nint(bufmat(iad+1))>0) THEN
360 opt_2 = nint(bufmat(iad+1))
361 indi_opt_2 = 3
362 add_option = 1
363 ! if network is used, then, viscosity is also used
364 opt_1 = one
365 indi_opt_1 = 2
366 add_option = 1
367 ENDIF
368 ELSEIF(mln==104) THEN
369 iad=ipm(7,abs(mid))-1
370 flag_nice_newton=nint(bufmat(iad+11))
371 IF(flag_nice_newton==2) THEN ! Newtow algo
372 indi = 2
373 ELSE ! Nice algo
374 indi = 1
375 ENDIF
376 flag_gurson=nint(bufmat(iad+30))
377 IF(flag_gurson/=0) THEN
378 special_option=1
379 spe_i_1 = 1
380 spe_i_2 = 1
381 ENDIF
382 IF(flag_gurson==1) THEN
383 spe_i_2 = 1
384 ELSEIF(flag_gurson==2) THEN
385 spe_i_2 = 2
386 ELSEIF(flag_gurson==3) THEN
387 spe_i_2 = 3
388 ENDIF
389 flag_non_local = mat_param(abs(mid))%NLOC
390 IF (mat_param(abs(mid))%IVISC > 0) THEN
391 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
392 ENDIF
393 ELSE
394 indi = 1
395 IF (mat_param(abs(mid))%IVISC > 0) THEN
396 visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
397 ENDIF
398 ENDIF
399 IF(ismstr_cost) add_over_cost = elm_over_cost(1)
400
401 mult_spe = 0.
402 spe_i_3 = 1
403 IF(flag_non_local/=0) THEN
404 spe_i_3 = 1
405 mult_spe = 1.
406 ENDIF
407 cost_check = 0
408!****************************************
409! ---------------------------
410! TETRA 4
411! ---------------------------
412 IF (isol==4.AND. (jsrot /= 1)) THEN
413 ! check if the (mid,pid) cost must be initialized from a previous run
414 IF(recherche==0.AND.test_poids/=0) THEN
415 poin_part = iparts(i)
416 poin_mid = poin_part_sol(1,poin_part,6)
417 poin_pid = poin_part_sol(2,poin_part,6)
418 IF(poin_mid/=0.AND.poin_pid/=0) THEN
419 IF(mid_pid_sol(poin_mid,6)%COST1D(poin_pid)/=zero) THEN
420 cost_check = 1
421 poin_elm_typ = 6
422 telt = mid_pid_sol(poin_mid,poin_elm_typ)%COST1D(poin_pid)
423 ENDIF
424 ENDIF
425 ENDIF
426 ! the (mid,pid) cost must be initialized from .inc file
427 IF(cost_check==0) THEN
428 IF( ddweights(1,1,iabs(mid))/=0)THEN
429 tmat = ddweights(1,1,iabs(mid)) * tpsref
430 ELSE
431 IF(mult/=0) tmatadd = mult * (tet4tnl(mln,indi)-tet4tnl(mln,indi2))
432 IF(add_option/=0) tmatadd = opt_1 * tet4tnl(mln,indi_opt_1) + opt_2 * tet4tnl(mln,indi_opt_2)
433 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
434 tmat = tet4tnl(mln,indi) + tmatadd
435 ENDIF
436! --------------
437! Failure
438 IF(nfail/=0) THEN
439 DO j=1,nfail
440 trup = trup + rupture_tet4(irup_tab(j),irup2)
441 ENDDO
442 ENDIF
443! --------------
444 telt = tmat + tet4telt(1) + trup + mult_spe*nlocal_option(spe_i_3) + add_over_cost + visc_prony
445 ENDIF
446!****************************************
447! ---------------------------
448! TETRA 10 or TETRA4 + JSROT
449! ---------------------------
450 ELSEIF ((isol==10).OR.(isol==4.AND. jsrot==1)) THEN
451 ! check if the (mid,pid) cost must be initialized from a previous run
452 IF(recherche==0.AND.test_poids/=0) THEN
453 IF(isol==10) THEN
454 poin_part = iparts(i)
455 poin_mid = poin_part_sol(1,poin_part,2)
456 poin_pid = poin_part_sol(2,poin_part,2)
457 ! if POIN_MID==0 and POIN_PID == 0, the element cost in the .ddw file is 0 --> must be initialized
458 ! from the .inc file
459 IF(poin_mid/=0.AND.poin_pid/=0) THEN
460 IF(mid_pid_sol(poin_mid,2)%COST1D(poin_pid)/=zero) THEN
461 cost_check = 1
462 poin_elm_typ = 2
463 telt = mid_pid_sol(poin_mid,poin_elm_typ)%COST1D(poin_pid)
464 ENDIF
465 ENDIF
466 ELSEIF(isol==4.AND. jsrot==1) THEN
467 poin_part = iparts(i)
468 poin_mid = poin_part_sol(1,poin_part,6)
469 poin_pid = poin_part_sol(2,poin_part,6)
470 ! if POIN_MID==0 and POIN_PID == 0, the element cost in the .ddw file is 0 --> must be initialized
471 ! from the .inc file
472 IF(poin_mid/=0.AND.poin_pid/=0) THEN
473 IF(mid_pid_sol(poin_mid,6)%COST1D(poin_pid)/=zero) THEN
474 cost_check = 1
475 poin_elm_typ = 6
476 telt = mid_pid_sol(poin_mid,poin_elm_typ)%COST1D(poin_pid)
477 ENDIF
478 ENDIF
479 ENDIF
480 ENDIF
481 ! the (mid,pid) cost must be initialized from .inc file
482 IF(cost_check==0) THEN
483 IF( ddweights(1,1,iabs(mid))/=0)THEN
484 tmat = ddweights(1,1,iabs(mid)) * tpsref
485 ELSE
486 IF(mult/=0) tmatadd = mult * (tet10tnl(mln,indi)-tet10tnl(mln,indi2))
487 IF(add_option/=0) tmatadd = opt_1 * tet10tnl(mln,indi_opt_1) + opt_2 * tet10tnl(mln,indi_opt_2)
488 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
489 tmat = tet10tnl(mln,indi) + tmatadd
490 ENDIF
491! --------------
492! failure
493 IF(nfail/=0) THEN
494 DO j=1,nfail
495 trup = trup + rupture_tet10(irup_tab(j),irup2)
496 ENDDO
497 ENDIF
498! --------------
499 IF(isol==10) telt = tet10telt(1)
500 IF(isol==4.AND. jsrot==1) telt = tet4telt(2)
501 telt = tmat + telt + trup + mult_spe*nlocal_option(spe_i_3) + 4.*(add_over_cost + visc_prony)
502 ENDIF
503 ELSE
504!****************************************
505! ---------------------------
506! SOLID ELEMENT
507! ---------------------------
508 ! check if the (mid,pid) cost must be initialized from a previous run
509 IF(recherche==0.AND.test_poids/=0) THEN
510 IF(isol==6) THEN
511 poin_elm_typ = 5
512 ELSEIF(isol==8) THEN
513 poin_elm_typ = 7
514 ELSEIF(isol==16) THEN
515 poin_elm_typ = 3
516 ELSEIF(isol==20) THEN
517 poin_elm_typ = 4
518 ELSE
519 poin_elm_typ = 1
520 ENDIF
521 poin_part = iparts(i)
522 poin_mid = poin_part_sol(1,poin_part,poin_elm_typ)
523 poin_pid = poin_part_sol(2,poin_part,poin_elm_typ)
524 ! if POIN_MID==0 and POIN_PID == 0, the element cost in the .ddw file is 0 --> must be initialized
525 ! from the .inc file
526 IF(poin_mid/=0.AND.poin_pid/=0) THEN
527 IF(mid_pid_sol(poin_mid,poin_elm_typ)%COST1D(poin_pid)/=zero) THEN
528 cost_check = 1
529 telt = mid_pid_sol(poin_mid,poin_elm_typ)%COST1D(poin_pid)
530 ENDIF
531 ENDIF
532 ENDIF
533 ! the (mid,pid) cost must be initialized from .inc file
534 IF(cost_check==0) THEN
535 IF (jhbe==1) THEN
536! Solides ISOLD1
537 IF( ddweights(1,1,iabs(mid))/=0)THEN
538 tmat = ddweights(1,1,iabs(mid)) * tpsref
539 ELSE
540 IF(mult/=0) tmatadd = mult * (sol1tnl(mln,indi)-sol1tnl(mln,indi2))
541 IF(add_option/=0) tmatadd = opt_1 * sol1tnl(mln,indi_opt_1) + opt_2 * sol1tnl(mln,indi_opt_2)
542 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
543 tmat = sol1tnl(mln,indi) + tmatadd
544 ENDIF
545! --------------
546! Failure
547 IF(nfail/=0) THEN
548 DO j=1,nfail
549 trup = trup + rupture_sol(irup_tab(j),irup2)
550 ENDDO
551 ENDIF
552! --------------
553 telt = tmat + soltelt(1) + trup + mult_spe*nlocal_option(spe_i_3) + add_over_cost + visc_prony
554 ELSEIF (jhbe==2) THEN
555! Solides ISOLD2
556 IF( ddweights(1,1,iabs(mid))/=0)THEN
557 tmat = ddweights(1,1,iabs(mid)) * tpsref
558 ELSE
559 IF(mult/=0) tmatadd = mult * (sol1tnl(mln,indi)-sol1tnl(mln,indi2))
560 IF(add_option/=0) tmatadd = opt_1 * sol1tnl(mln,indi_opt_1) + opt_2 * sol1tnl(mln,indi_opt_2)
561 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
562 tmat = sol1tnl(mln,indi) + tmatadd
563 ENDIF
564! --------------
565! Failure
566 IF(nfail/=0) THEN
567 DO j=1,nfail
568 trup = trup + rupture_sol(irup_tab(j),irup2)
569 ENDDO
570 ENDIF
571! --------------
572 telt = tmat + soltelt(2) + trup + mult_spe*nlocal_option(spe_i_3) + add_over_cost + visc_prony
573 ELSEIF (jhbe==24.OR.jhbe==104) THEN
574! Solides ISOLD24 - HEPH
575 IF( ddweights(1,1,iabs(mid))/=0)THEN
576 tmat = ddweights(1,1,iabs(mid)) * tpsref
577 ELSE
578 IF(mult/=0) tmatadd = mult * (sol1tnl(mln,indi)-sol1tnl(mln,indi2))
579 IF(add_option/=0) tmatadd = opt_1 * sol1tnl(mln,indi_opt_1) + opt_2 * sol1tnl(mln,indi_opt_2)
580 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
581 tmat = sol1tnl(mln,indi) + tmatadd
582 ENDIF
583! --------------
584! Failure
585 IF(nfail/=0) THEN
586 DO j=1,nfail
587 trup = trup + rupture_sol(irup_tab(j),irup2)
588 ENDDO
589 ENDIF
590! --------------
591 telt = tmat + soltelt(3) + trup + mult_spe*nlocal_option(spe_i_3) + add_over_cost + visc_prony
592C
593 ELSEIF (jhbe==12) THEN
594! Solids Insd12 - STD 8 NODE INTEG Point
595 IF( ddweights(1,1,iabs(mid))/=0)THEN
596 tmat = ddweights(1,1,iabs(mid)) * tpsref
597 ELSE
598 IF(mult/=0) tmatadd = mult * (sol8tnl(mln,indi)-sol8tnl(mln,indi2))
599 IF(add_option/=0) tmatadd = opt_1 * sol8tnl(mln,indi_opt_1) + opt_2 * sol8tnl(mln,indi_opt_2)
600 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
601 tmat = sol8tnl(mln,indi) + tmatadd
602 ENDIF
603! --------------
604! Failure
605 IF(nfail/=0) THEN
606 DO j=1,nfail
607 trup = trup + rupture_sol12(irup_tab(j),irup2)
608 ENDDO
609 ENDIF
610! --------------
611 telt = tmat + soltelt(4) + trup + mult_spe*nlocal_option(spe_i_3) + 8.*(add_over_cost + visc_prony)
612 ELSEIF ( (jhbe==14.OR.(jhbe>=222.AND.jhbe<=999)).AND.(igt/=20.AND.igt/=21.AND.igt/=22)) THEN
613! Solides HA8
614 mpt = abs(npn)
615 nptr = max(mpt/100,1)
616 npts = max(mod(mpt/10,10),1)
617 nptt = max(mod(mpt,10),1)
618 nptot = npts*nptt*nptr
619
620 IF( ddweights(1,1,iabs(mid))/=0)THEN
621 tmat = ddweights(1,1,iabs(mid)) * tpsref
622 ELSE
623 IF(mult/=0) tmatadd = mult * (sol1tnl(mln,indi)-sol1tnl(mln,indi2))
624 IF(add_option/=0) tmatadd = opt_1 * sol1tnl(mln,indi_opt_1) + opt_2 * sol1tnl(mln,indi_opt_2)
625 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
626 tmat = sol1tnl(mln,indi) + tmatadd
627 ENDIF
628! --------------
629! Failure
630 IF(nfail/=0) THEN
631 DO j=1,nfail
632 trup = trup + rupture_sol(irup_tab(j),irup2)
633 ENDDO
634 ENDIF
635! --------------
636 ! 8 NPT = 222 = reference weight
637 ! if NPT > 8, element weight = reference weight + (NPT-8) * overweight
638 overcost_elm = 0
639 IF(nptot>8) overcost_elm = nptot-8
640 telt = nptot*(tmat+trup+add_over_cost+visc_prony)+soltelt(5) +overcost_elm *soltelt(6) +
641 . mult_spe*nlocal_option(spe_i_3)
642 ELSEIF(jhbe==14.AND.(igt==20.OR.igt==21.OR.igt==22)) THEN
643! Solides Thick shell
644 mpt = abs(npn)
645 nptr = max(mpt/100,1)
646 npts = max(mod(mpt/10,10),1)
647 nptt = max(mod(mpt,10),1)
648 nptot = npts*nptt*nptr
649
650 IF( ddweights(1,1,iabs(mid))/=0)THEN
651 tmat = ddweights(1,1,iabs(mid)) * tpsref
652 ELSE
653 IF(mult/=0) tmatadd = mult * (sol1tnl(mln,indi)-sol1tnl(mln,indi2))
654 IF(add_option/=0) tmatadd = opt_1 * sol1tnl(mln,indi_opt_1) + opt_2 * sol1tnl(mln,indi_opt_2)
655 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
656
657 number_layer = 0
658 ! ---------------------
659 ! check the number of layer
660 IF(igeo(30,pid)>9) THEN ! number of layer > 9 --> 1 int. point in the "layer direction" + number of layer = NLAY
661 number_layer = igeo(30,pid)
662 icpr = igeo(14,pid)
663 ! ICPR = ijk = rst ( i=r / j=s / k=t)
664 IF(icpr==100) THEN ! r direction
665 overcost_elm = npts*nptt
666 ELSEIF(icpr==10) THEN ! sharp
667 overcost_elm = nptt*nptr
668 ELSE ! Directorate
669 overcost_elm = npts*nptr
670 ENDIF
671 ELSE
672 ! number of layer <= 9 --> number of layer = ICPR direction
673 icpr = igeo(14,pid)
674 number_layer = igeo(30,pid)
675 ! ICPR = ijk = rst ( i=r / j=s / k=t)
676 IF(icpr==100) THEN ! r direction
677 number_layer = nptr
678 overcost_elm = npts*nptt
679 ELSEIF(icpr==10) THEN ! Sharp
680 number_layer = npts
681 overcost_elm = nptt*nptr
682 ELSE ! Directorate
683 number_layer = nptt
684 overcost_elm = npts*nptr
685 ENDIF
686 ENDIF
687 ! ---------------------
688 ! check if composite material is used
689 ! --> if yes, add & sum the different material costs
690 IF(composite_option) THEN
691 DO nlay=1,number_layer
692 composite_mid = igeo(100+nlay,pid)
693 composite_mln = nint(pm(19,abs(composite_mid)))
694 tmatadd = tmatadd + sol1tnl(composite_mln,indi)
695 ENDDO
696 tmatadd = tmatadd - sol1tnl(mln,indi)
697 ENDIF
698 tmat = sol1tnl(mln,indi) + tmatadd
699 ENDIF
700! --------------
701! Failure
702 IF(nfail/=0) THEN
703 DO j=1,nfail
704 trup = trup + rupture_sol(irup_tab(j),irup2)
705 ENDDO
706 ENDIF
707! --------------
708 ! 4 NPT = 202 = 4 * reference weight
709 ! if NPT > 4, element weight = reference weight + overcost
710 telt = overcost_elm*(tmat+visc_prony)+nptot*trup +
711 . overcost_elm*number_layer*soltelt(10) + mult_spe*nlocal_option(spe_i_3) +
712 . overcost_elm * add_over_cost
713 ELSEIF(jhbe==15) THEN
714 ! --------------
715 ! Solid Thick shell
716 ! --------------
717 IF( ddweights(1,1,iabs(mid))/=0)THEN
718 tmat = ddweights(1,1,iabs(mid)) * tpsref
719 ELSE
720 IF(mult/=0) tmatadd = mult * (sol1tnl(mln,indi)-sol1tnl(mln,indi2))
721 IF(add_option/=0) tmatadd = opt_1 * sol1tnl(mln,indi_opt_1) + opt_2 * sol1tnl(mln,indi_opt_2)
722 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
723 tmat = sol1tnl(mln,indi) + tmatadd
724 ENDIF
725! --------------
726! Failure
727 IF(nfail/=0) THEN
728 DO j=1,nfail
729 trup = trup + rupture_sol(irup_tab(j),irup2)
730 ENDDO
731 ENDIF
732
733 nptot = abs(npn)
734 ! --------------
735 ! jhbe=15 : weight = element weight + element extracost + mat + rupture
736 telt = nptot*(tmat+trup+visc_prony) + soltelt(11) + nptot*soltelt(12) +
737 . mult_spe*nlocal_option(spe_i_3) + add_over_cost
738 ! --------------
739 ELSEIF (jhbe==17) THEN
740! ISOLIDE=17 H8C Standard 8-nodes compatible
741! solid full integration formulation 2*2*2 integration points, no hourglass
742 IF( ddweights(1,1,iabs(mid))/=0)THEN
743 tmat = ddweights(1,1,iabs(mid)) * tpsref
744 ELSE
745 IF(mult/=0) tmatadd = mult * (sol1tnl(mln,indi)-sol1tnl(mln,indi2))
746 IF(add_option/=0) tmatadd = opt_1 * sol1tnl(mln,indi_opt_1) + opt_2 * sol1tnl(mln,indi_opt_2)
747 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
748 tmat = sol1tnl(mln,indi) + tmatadd
749 ENDIF
750! --------------
751! Failure
752 IF(nfail/=0) THEN
753 DO j=1,nfail
754 trup = trup + rupture_sol(irup_tab(j),irup2)
755 ENDDO
756 ENDIF
757! --------------
758 telt = (tmat+trup+add_over_cost+visc_prony)*8 + soltelt(7) + mult_spe*nlocal_option(spe_i_3)
759 ELSEIF (jhbe==18) THEN
760! ISOLIDE=18
761 IF( ddweights(1,1,iabs(mid))/=0)THEN
762 tmat = ddweights(1,1,iabs(mid)) * tpsref
763 ELSE
764 IF(mult/=0) tmatadd = mult * (sol1tnl(mln,indi)-sol1tnl(mln,indi2))
765 IF(add_option/=0) tmatadd = opt_1 * sol1tnl(mln,indi_opt_1) + opt_2 * sol1tnl(mln,indi_opt_2)
766 IF(special_option/=0) tmatadd = tmatadd + sol_option(spe_i_1,spe_i_2)
767 tmat = sol1tnl(mln,indi) + tmatadd
768 ENDIF
769! --------------
770! Failure
771 IF(nfail/=0) THEN
772 DO j=1,nfail
773 trup = trup + rupture_sol(irup_tab(j),irup2)
774 ENDDO
775 ENDIF
776! --------------
777 telt = (tmat+trup+add_over_cost+visc_prony)*8 + soltelt(9) + mult_spe*nlocal_option(spe_i_3)
778 ELSE
779! --------------
780! Failure
781 IF(nfail/=0) THEN
782 DO j=1,nfail
783 trup = trup + rupture_sol(irup_tab(j),irup2)
784 ENDDO
785 ENDIF
786! --------------
787 telt = sol1tnl(mln,1) + soltelt(1) + trup
788 . + mult_spe*nlocal_option(spe_i_3) + add_over_cost + visc_prony
789 ENDIF
790 ENDIF
791 ENDIF
792C
793 poids = telt * invtref
794 ! --------------------
795 IF(recherche==0) THEN
796 IF (wd_local==0..AND.mln/=0)THEN
797 wd(i) = poids
798 poin_part = iparts(i)
799 IF (isol==4.AND. (jsrot /= 1)) THEN
800 poin_elm_typ = 6
801 ELSEIF( (isol==10).OR.(isol==4.AND. jsrot==1) ) THEN
802 IF(isol==10) THEN
803 poin_elm_typ = 2
804 ELSE
805 poin_elm_typ = 6
806 ENDIF
807 ELSE
808 IF(isol==6) THEN
809 poin_elm_typ = 5
810 ELSEIF(isol==8) THEN
811 poin_elm_typ = 7
812 ELSEIF(isol==16) THEN
813 poin_elm_typ = 3
814 ELSEIF(isol==20) THEN
815 poin_elm_typ = 4
816 ELSE
817 poin_elm_typ = 1
818 ENDIF
819 ENDIF
820 poin_part = iparts(i)
821 poin_mid = poin_part_sol(1,poin_part,poin_elm_typ)
822 poin_pid = poin_part_sol(2,poin_part,poin_elm_typ)
823 IF(poin_mid/=0.AND.poin_pid/=0) mid_pid_sol(poin_mid,poin_elm_typ)%COST1D(poin_pid) = telt
824 ELSEIF(mln==0)THEN
825 wd(i) = 0.0001
826 END IF
827 ELSE
828 telt_pro = telt
829 ENDIF
830 ! --------------------
831 ENDDO
832 RETURN
833 END
initwg_solid
subroutine initwg_solid(wd, pm, geo, ixs, igeo, isolnod, numels, ipm, size_irup, nummat, numgeo, poin_part_sol, mid_pid_sol, iparts, bufmat, mid_old, pid_old, mln_old, recherche, isol_old, telt_pro, tabmp_l, npart, mat_param)
Definition initwg_solid.F:39
max
#define max(a, b)
Definition macros.h:21
ddweights_mod
Definition ddweights_mod.F:40
mid_pid_mod
Definition mid_pid_mod.F:36
bidon
subroutine bidon
Definition machine.F:41
mid_pid_mod::mid_pid_type
Definition mid_pid_mod.F:44
visc_prony
subroutine visc_prony(visc, nprony, nel, nvarvis, uvarvis, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sv1, sv2, sv3, sv4, sv5, sv6, timestep, rho, viscmax, soundsp, nvar_damp)
Definition visc_prony.F:34