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