OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
inigrav_load.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!|| inigrav_load ../starter/source/initial_conditions/inigrav/inigrav_load.F
25!||--- called by ------------------------------------------------------
26!|| initia ../starter/source/elements/initia/initia.f
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| check_is_on_segment ../starter/source/initial_conditions/inigrav/inigrav_load.F
30!|| inigrav_eos ../starter/source/initial_conditions/inigrav/inigrav_eos.F
31!|| inigrav_m37 ../starter/source/initial_conditions/inigrav/inigrav_m37.F
32!|| inigrav_m51 ../starter/source/initial_conditions/inigrav/inigrav_m51.F
33!|| inigrav_part_list ../starter/source/initial_conditions/inigrav/inigrav_part_list.F
34!||--- uses -----------------------------------------------------
35!|| inigrav ../starter/share/modules1/inigrav_mod.F
36!|| message_mod ../starter/share/message_module/message_mod.F
37!||====================================================================
38 SUBROUTINE inigrav_load(ELBUF_TAB ,IPART ,IGRPART ,IPARG ,IPARTTG ,
39 1 IPARTS ,IPARTQ ,X ,IXS ,IXQ ,
40 2 IXTG ,PM ,IPM ,BUFMAT ,MULTI_FVM,
41 3 ALE_CONNECTIVITY ,NV46 ,IGRSURF ,ITAB ,EBCS_TAB,
42 4 NPF ,TF ,MAT_PARAM)
43C-----------------------------------------------
44C Description
45C-----------------------------------------------
46C This subroutine is computing distance from a given surface along vector provided by /GRAV option.
47C Surface can be planar or provided by user.
48C
49C REMINDER :
50!
51! IGRSURF(IGS)%ELEM(J) :: element attached to the segment(J) of the surface
52! IGRSURF(IGS)%ELTYP(J) :: type of element attached to the segment of the
53! ITYP = 0 - surf of segments
54! ITYP = 1 - surf of solids
55! ITYP = 2 - surf of quads
56! ITYP = 3 - surf of SH4N
57! ITYP = 4 - line of trusses
58! ityp = 5 - line of beams
59! ITYP = 6 - line of springs
60! ITYP = 7 - surf of SH3N
61! ITYP = 8 - line of XELEM (nstrand element)
62! ITYP = 101 - ISOGEOMETRIQUE
63!
64C-----------------------------------------------
65C M o d u l e s
66C-----------------------------------------------
67 USE elbufdef_mod
68 USE inigrav
69 USE message_mod
70 USE multi_fvm_mod
71 USE groupdef_mod
72 USE ebcs_mod
73 USE ale_connectivity_mod
74 USE matparam_def_mod
75C-----------------------------------------------
76C I m p l i c i t T y p e s
77C-----------------------------------------------
78#include "implicit_f.inc"
79C-----------------------------------------------
80C G l o b a l P a r a m e t e r s
81C-----------------------------------------------
82#include "mvsiz_p.inc"
83C-----------------------------------------------
84C C o m m o n B l o c k s
85C-----------------------------------------------
86#include "vect01_c.inc"
87#include "param_c.inc"
88#include "com01_c.inc"
89#include "com04_c.inc"
90#include "scr17_c.inc"
91#include "tabsiz_c.inc"
92C-----------------------------------------------
93C D u m m y A r g u m e n t s
94C-----------------------------------------------
95 INTEGER, INTENT(IN),TARGET :: IPART(LIPART1,*)
96 INTEGER, INTENT(IN) :: IPARTS(NUMELS), IPARTQ(NUMELQ), IPARTTG(NUMELTG)
97 INTEGER, INTENT(IN) :: IPM(NPROPMI, NUMMAT),ITAB(NUMNOD)
98 INTEGER, INTENT(IN) :: IPARG(NPARG,NGROUP),NV46
99 INTEGER, INTENT(IN), TARGET :: IXS(NIXS, NUMELS), IXQ(NIXQ, NUMELQ), IXTG(NIXTG, NUMELTG)
100 TYPE(elbuf_struct_),DIMENSION(NGROUP),TARGET :: ELBUF_TAB
101 my_real, INTENT(IN) :: x(3,numnod), pm(npropm,nummat), bufmat(*)
102 TYPE(multi_fvm_struct), INTENT(IN) :: MULTI_FVM
103 TYPE (SURF_), DIMENSION(NSURF),INTENT(IN) :: IGRSURF
104 TYPE(t_ebcs_tab), INTENT(INOUT) :: EBCS_TAB
105C-----------------------------------------------
106 TYPE (GROUP_) , DIMENSION(NGRPART) :: IGRPART
107 TYPE(t_ale_connectivity), INTENT(INOUT) :: ALE_CONNECTIVITY
108 INTEGER,INTENT(IN)::NPF(SNPC)
109 my_real,INTENT(IN)::TF(STF)
110 TYPE(matparam_struct_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM
111C-----------------------------------------------
112C L o c a l V a r i a b l e s
113C-----------------------------------------------
114 my_real :: grav0,bx,by,bz,nx,ny,nz,dotprod,dotp,rho0,ngx,ngy,ngz,alpha,interp(3),vala,valb,norm
115 my_real :: z(3), depth(mvsiz), n(3,8),vec(3), delta_p(mvsiz),psurf,pgrav(mvsiz),vfrac,p(3,1:4)
116 my_real :: lambda,diag1(3),diag2(3),tol,dist,b(3),volfrac
117 my_real,ALLOCATABLE,DIMENSION(:,:) :: min_coor, max_coor, n_surf, zf
118 INTEGER :: K,KK,IGRP,ISURF,IGRAV,NPART_IN_GROUP,MLW,NG,IAD0,II,NEL,J,IE,I,IERROR,IMAT,M(MVSIZ),Q,NSEG,ISEG
119 INTEGER :: USERID
120 INTEGER :: NOD(8), MATID, IADBUF,IFORM,NEL2,LIST(MVSIZ),NBMAT,MATLAW,NIX,UID,SURF_TYPE,IELTYP
121 INTEGER, ALLOCATABLE, DIMENSION(:) :: IP
122 TYPE(g_bufel_) ,POINTER :: GBUF
123 TYPE(l_bufel_), POINTER :: LBUF
124 LOGICAL :: lCYCLE,lERROR,lTRIA,lFOUND_PROJ,lTEST,l_PLANAR_SURF,l_USER_SURF
125 INTEGER, DIMENSION(:, :), POINTER :: IX
126 character*64 :: ERRMSG
127C-----------------------------------------------
128C S o u r c e L i n e s
129C-----------------------------------------------
130
131 CALL inigrav_part_list( ipart ,igrpart, ebcs_tab ) !list related parts to apply gravity loading with /ebcs/nrf
132
133 DO k=1,ninigrav
134 igrp = inigrv(1,k)
135 isurf = inigrv(2,k)
136 igrav = inigrv(3,k)
137 uid = inigrv(4,k)
138 bx = linigrav(01,k)
139 by = linigrav(02,k)
140 bz = linigrav(03,k)
141 nx = linigrav(04,k)
142 ny = linigrav(05,k)
143 nz = linigrav(06,k)
144 grav0 = linigrav(07,k)
145 ngx = linigrav(08,k)
146 ngy = linigrav(09,k)
147 ngz = linigrav(10,k)
148 psurf = linigrav(11,k)
149 IF (grav0==zero)cycle !nothing to do if g=0 at time=0 => Delta_P=0 => P=P(t=0)=P0
150 !READING GRPART DATA
151 IF(igrp/=0)THEN
152 npart_in_group = igrpart(igrp)%NENTITY
153 ALLOCATE(ip(npart_in_group),stat=ierror)
154 IF(ierror/=0)CALL ancmsg(msgid = 268,
155 . anmode = aninfo,
156 . msgtype = msgerror,
157 . c1 = 'INIGRAV')
158 ip(1:npart_in_group) = ipart(4,igrpart(igrp)%ENTITY(1:npart_in_group))
159 ELSE
160 !by default treat all parts
161 iad0 = 0
162 npart_in_group = 0
163 ENDIF
164 l_planar_surf = .false.
165 l_user_surf = .false.
166 surf_type=-1 !-1:unaffected
167 IF(isurf>0)surf_type=igrsurf(isurf)%TYPE
168 IF(surf_type==200 .OR. surf_type == -1)THEN
169 l_planar_surf=.true.
170 ELSEIF(surf_type == 0) THEN
171 l_user_surf=.true.
172 nseg = igrsurf(isurf)%NSEG
173 !STORING DATA FROM USER SURFACE (BOX,NORMAL,CENTROID)
174 ALLOCATE(max_coor(3,nseg),min_coor(3,nseg),n_surf(3,nseg),zf(3,nseg))
175 DO iseg=1,nseg
176 ltria=.false.
177 IF(igrsurf(isurf)%NODES(iseg,4)==0)ltria=.true.
178 IF(ltria)THEN
179 p(1:3,4)=p(1:3,3)
180 ENDIF
181 !face coordinates
182 p(1,1:4)=x(1,igrsurf(isurf)%NODES(iseg,1:4))
183 p(2,1:4)=x(2,igrsurf(isurf)%NODES(iseg,1:4))
184 p(3,1:4)=x(3,igrsurf(isurf)%NODES(iseg,1:4))
185 !box (to skip cell centroid if too far from it)
186 min_coor(1,iseg)=min(p(1,1),p(1,2),p(1,3),p(1,4))
187 min_coor(2,iseg)=min(p(2,1),p(2,2),p(2,3),p(2,4))
188 min_coor(3,iseg)=min(p(3,1),p(3,2),p(3,3),p(3,4))
189 max_coor(1,iseg)=max(p(1,1),p(1,2),p(1,3),p(1,4))
190 max_coor(2,iseg)=max(p(2,1),p(2,2),p(2,3),p(2,4))
191 max_coor(3,iseg)=max(p(3,1),p(3,2),p(3,3),p(3,4))
192 !face normal
193 diag1(1)=(p(1,3)-p(1,1))
194 diag1(2)=(p(2,3)-p(2,1))
195 diag1(3)=(p(3,3)-p(3,1))
196 diag2(1)=(p(1,4)-p(1,2))
197 diag2(2)=(p(2,4)-p(2,2))
198 diag2(3)=(p(3,4)-p(3,2))
199 n_surf(1,iseg)= diag1(2)*diag2(3)-diag1(3)*diag2(2)
200 n_surf(2,iseg)=-diag1(1)*diag2(3)+diag1(3)*diag2(1)
201 n_surf(3,iseg)= diag1(1)*diag2(2)-diag1(2)*diag2(1)
202 norm = sqrt(n_surf(1,iseg)*n_surf(1,iseg)+n_surf(2,iseg)*n_surf(2,iseg)+n_surf(3,iseg)*n_surf(3,iseg))
203 n_surf(1,iseg)=n_surf(1,iseg)/norm
204 n_surf(2,iseg)=n_surf(2,iseg)/norm
205 n_surf(3,iseg)=n_surf(3,iseg)/norm
206 !face centroid
207 IF(ltria)THEN
208 zf(1,iseg)=third*sum(p(1,1:3))
209 zf(2,iseg)=third*sum(p(2,1:3))
210 zf(3,iseg)=third*sum(p(3,1:3))
211 ELSE
212 zf(1,iseg)=fourth*sum(p(1,1:4))
213 zf(2,iseg)=fourth*sum(p(2,1:4))
214 zf(3,iseg)=fourth*sum(p(3,1:4))
215 ENDIF
216 ENDDO
217 !define a tolerance to check if cell centroid too far from face (box)
218 iseg=1
219 tol=em02*max(max_coor(1,iseg)-min_coor(1,iseg), max_coor(2,iseg)-min_coor(2,iseg), max_coor(3,iseg)-min_coor(3,iseg))
220 ENDIF
221
222
223 DO ng=1,ngroup
224 mlw = iparg(1,ng)
225 ity = iparg(5,ng)
226 nel = iparg(2,ng)
227 nft = iparg(3,ng)
228 iad = iparg(4,ng) - 1
229 gbuf => elbuf_tab(ng)%GBUF
230 pgrav(1:nel)=zero
231 IF(mlw==0)cycle
232 IF( (ity == 1. .AND.n2d==0) .OR. ((ity == 2. .OR. ity == 7).AND.n2d/=0) )THEN
233 IF(ity==1)THEN
234 imat = ixs(1,1+nft)
235 m(1:nel) = iparts(1+nft:nel+nft)
236 ELSEIF(ity==2)THEN
237 imat = ixq(1,1+nft)
238 m(1:nel) = ipartq(1+nft:nel+nft)
239 ELSEIF(ity==7)THEN
240 imat = ixtg(1,1+nft)
241 m(1:nel) = iparttg(1+nft:nel+nft)
242 ENDIF
243 !check if PART must be treated or not
244 list(1:mvsiz) = 0
245 nel2 = 0
246 !il peut y avoir plusieurs PARTS par groupe. Tous les elements ont meme mat et meme PROP.
247 ! creation dune list LIST(1:NEL2) c 1:NEL avec les elements des parts a traiter
248 !===================================================================!
249 ! SELECTION DES ELEMENTS !
250 !===================================================================!
251 IF(igrp/=0)THEN
252 lcycle = .true.
253 DO i=1,nel
254 DO q=1,npart_in_group
255 IF(ipart(4,m(i))==ip(q))THEN
256 nel2 = nel2+1
257 list(nel2) = i
258 lcycle = .false.
259 EXIT
260 ENDIF
261 ENDDO
262 enddo!next I
263 IF(lcycle)cycle
264 ELSE
265 nel2=nel
266 DO i=1,nel
267 list(i) = i !(/1:NEL/) !identite
268 ENDDO
269 ENDIF
270 ! LIST est la sous liste de 1:NEL des elements dont la PART est a traiter avec INIGRAV
271
272 ! This group is going to be treated since PART_ID in GRPART_ID
273 !===================================================================!
274 ! CALCUL DES CENTROIDES !
275 !===================================================================!
276 DO kk=1,nel2
277 i = list(kk)
278 ie = i + nft
279 userid = 0
280 IF(n2d==0)THEN
281 n(1:3,1) = x(1:3,ixs(2,ie))
282 n(1:3,2) = x(1:3,ixs(3,ie))
283 n(1:3,3) = x(1:3,ixs(4,ie))
284 n(1:3,4) = x(1:3,ixs(5,ie))
285 n(1:3,5) = x(1:3,ixs(6,ie))
286 n(1:3,6) = x(1:3,ixs(7,ie))
287 n(1:3,7) = x(1:3,ixs(8,ie))
288 n(1:3,8) = x(1:3,ixs(9,ie))
289 z(1) = one_over_8*(sum(n(1,1:8)))
290 z(2) = one_over_8*(sum(n(2,1:8)))
291 z(3) = one_over_8*(sum(n(3,1:8)))
292 userid = ixs(nixs,ie)
293 ELSEIF(ity==2)THEN
294 n(1:3,1) = x(1:3,ixq(2,ie))
295 n(1:3,2) = x(1:3,ixq(3,ie))
296 n(1:3,3) = x(1:3,ixq(4,ie))
297 n(1:3,4) = x(1:3,ixq(5,ie))
298 z(1) = fourth*(sum(n(1,1:4)))
299 z(2) = fourth*(sum(n(2,1:4)))
300 z(3) = fourth*(sum(n(3,1:4)))
301 userid = ixq(nixq,ie)
302 ELSEIF(ity==7)THEN
303 n(1:3,1) = x(1:3,ixtg(2,ie))
304 n(1:3,2) = x(1:3,ixtg(3,ie))
305 n(1:3,3) = x(1:3,ixtg(4,ie))
306 n(1:3,4) = zero
307 z(1) = third*(sum(n(1,1:3)))
308 z(2) = third*(sum(n(2,1:3)))
309 z(3) = third*(sum(n(3,1:3)))
310 userid = ixtg(nixtg,ie)
311 ENDIF
312 !===================================================================!
313 IF(l_planar_surf)THEN
314 !===================================================================!
315 ! --- CALCUL DES DISTANCES --- !
316 !-PLANAR SURFACE : IGRSURF(ISURF)%TYPE = 200 OR ISURF=0 !
317 !===================================================================!
318 !---NORMAL : NX,NY,NZ
319 !---BASIS POINT : BX,BY,BZ
320 !intersection point such as INTERP = Z + ALPHA . NG et <B-INTERP,N> = 0, donc :
321 alpha = (bx-z(1))*nx+(by-z(2))*ny+(bz-z(3))*nz / (ngx*nx+ngy*ny+ngz*nz) !non nul car verification a la lecture
322 interp(1) = z(1) + alpha*ngx
323 interp(2) = z(2) + alpha*ngy
324 interp(3) = z(3) + alpha*ngz
325 ! profondeur = <INTERP_P , N>
326 !DEPTH(I) = (INTERP(1)-Z(1))*(INTERP(1)-Z(1))+(INTERP(2)-Z(2))*(INTERP(2)-Z(2))+(INTERP(3)-Z(3))*(INTERP(3)-Z(3))
327 !DEPTH(I) = SQRT(DEPTH(I))
328 depth(kk)= (interp(1)-z(1))*ngx+(interp(2)-z(2))*ngy+(interp(3)-z(3))*ngz
329 depth(kk) = -depth(kk)
330 lerror=.false.
331 lfound_proj=.true.
332 ELSEIF(l_user_surf) THEN
333 !===================================================================!
334 ! --- calcul des distances --- !
335 !-USER SURFACE : IGRSURF(ISURF)%TYPE=0 , ISURF/=0 !
336 !===================================================================!
337 !LOOP ON CELL, THEN LOOP ON FACE : complexity is O(NCELL*NFACE), can be optimized with space partitioning
338 lerror=.false.
339 lfound_proj=.false.
340 DO iseg=1,nseg
341
342 IF(ngx==zero)THEN
343 IF(z(1)+tol < min_coor(1,iseg) .OR. z(1)-tol > max_coor(1,iseg))cycle
344 ENDIF
345 IF(ngy==zero)THEN
346 IF(z(2)+tol < min_coor(2,iseg) .OR. z(2)-tol > max_coor(2,iseg))cycle
347 ENDIF
348 IF(ngz==zero)THEN
349 IF(z(3)+tol < min_coor(3,iseg) .OR. z(3)-tol > max_coor(3,iseg))cycle
350 ENDIF
351 dotp = n_surf(1,iseg)*ngx + n_surf(2,iseg)*ngy + n_surf(3,iseg)*ngz
352 IF(abs(dotp)<=em04)THEN
353 CALL ancmsg(msgid = 90,
354 . anmode = aninfo,
355 . msgtype = msgerror,
356 . i1=uid,
357 . c1="INPUT SURFACE HAS INCOMPATIBLE SLOPE WITH GRAVITY DIRECTION"
358 . )
359 lerror=.true.
360 EXIT
361 ENDIF
362 p(1,1:4)=x(1,igrsurf(isurf)%NODES(iseg,1:4))
363 p(2,1:4)=x(2,igrsurf(isurf)%NODES(iseg,1:4))
364 p(3,1:4)=x(3,igrsurf(isurf)%NODES(iseg,1:4))
365 ! Let INTERP the normal projection point on plane generated by current face ISEG
366 ! interp = z + lambda*ng(colinear to gravity direction , ng gravity vector, z cell centroid, lambda unknown scalar)
367 ! <INTERP-B,N_SURF> = 0 (must satisfying plane equation, B basis point, Z cell centroid)
368 ! => INTERP = ...
369 b(1:3)=zf(1:3,iseg)
370 dist=(b(1)-z(1))*(b(1)-z(1)) + (b(2)-z(2))*(b(2)-z(2)) + (b(3)-z(3))*(b(3)-z(3))
371 IF(dist<=tol*tol)THEN
372 !use another basis point if face centroid superimposed with cell centroid
373 b(1:3)=p(1:3,1)
374 ENDIF
375 lambda = (b(1)-z(1))*n_surf(1,iseg) + (b(2)-z(2))*n_surf(2,iseg) + (b(3)-z(3))*n_surf(3,iseg)
376 lambda = lambda / dotp
377 interp(1) = z(1)+lambda*ngx
378 interp(2) = z(2)+lambda*ngy
379 interp(3) = z(3)+lambda*ngz
380 !CHECK INTERP IS INSIDE FACE, OTHERWISE SKIP IT
381 CALL check_is_on_segment(p(1:3,1:4),b,interp,ltria,ltest)
382 IF(.NOT.ltest)cycle
383 depth(kk)= (interp(1)-z(1))*ngx+(interp(2)-z(2))*ngy+(interp(3)-z(3))*ngz
384 depth(kk) = -depth(kk)
385 lfound_proj=.true.
386 ENDDO
387 IF(lerror)EXIT
388 IF(.NOT.lfound_proj)THEN
389 errmsg = "UNABLE TO PROJECT ON SURFACE CENTROID FROM CELL ID= "
390 WRITE(errmsg(52:62),fmt='(I10)')userid
391 CALL ancmsg(msgid = 90,
392 . anmode = aninfo,
393 . msgtype = msgerror,
394 . i1=uid,
395 . c1=errmsg(1:62)
396 . )
397 EXIT
398 ENDIF
399 ENDIF !l_PLANAR_SURF elif l_USER_SURF
400 ENDDO !next I
401 IF(lerror)EXIT
402 IF(.NOT.lfound_proj)EXIT
403 !===================================================================!
404 ! INITIALISATION DE L'ETAT INITIAL !
405 !===================================================================!
406 IF (n2d == 0) THEN
407 matid = ixs(1, 1 + nft)
408 ELSEIF(ity==2)THEN
409 matid = ixq(1, 1 + nft)
410 ELSEIF(ity==7)THEN
411 matid = ixtg(1, 1 + nft)
412 ENDIF
413 iadbuf = max(1,ipm(7, matid)) ! Adress of data for the material in the material buffer
414 !NEL2 is size of subgroup in LIST(1:NEL2), NEL is size of element group needed to shift indexes in MBUF%VAR and GBUF%SIG
415 SELECT CASE (mlw)
416 !=======================================!
417 ! 3,4,6,49 (ANY EOS) !
418 !=======================================!
419 CASE (3, 4, 6, 10, 49) !Material law defined with Equations Of States
420 rho0 = pm(1,matid)
421 CALL inigrav_eos(nel,nel2, ng, matid, ipm, grav0, rho0, depth, pm, bufmat(iadbuf),
422 . elbuf_tab, psurf, list , pgrav, 0 ,mlw, npf, tf ,nummat,mat_param)
423 !=======================================!
424 ! MM-ALE 37 !
425 !=======================================!
426 CASE (37)
427 CALL inigrav_m37(nel,nel2, ng, matid, ipm, grav0, depth, pm, bufmat(iadbuf), elbuf_tab, psurf, list)
428 !=======================================!
429 ! MM-ALE 51 !
430 !=======================================!
431 CASE (51)
432 IF (n2d == 0) THEN ! HEXA
433 ix => ixs(1:nixs, 1:numels)
434 nix = nixs
435 ELSEIF (ity == 2) then! QUADS
436 ix => ixq(1:nixq, 1:numelq)
437 nix = nixq
438 ELSEIF (ity == 7) then! TRIAS
439 ix => ixtg(1:nixtg, 1:numeltg)
440 nix = nixtg
441 ENDIF
442 iform = bufmat(iadbuf + 31 - 1)
443 IF (iform /= 3) THEN
444 CALL inigrav_m51(nel , nel2, ng, matid, ipm, grav0, depth, pm, bufmat(iadbuf), elbuf_tab,
445 . psurf, list, ale_connectivity, nv46, ix , nix , nft , bufmat, iparg)
446 ELSE
447 CALL ancmsg(msgid = 84,
448 . anmode = aninfo,
449 . msgtype = msgwarning
450 . )
451 ENDIF
452 !=======================================!
453 ! MULTI-FLUID 151 !
454 !=======================================!
455 CASE(151)
456
457 IF (n2d == 0) THEN ! HEXA
458 ix => ixs(1:nixs, lft + nft:llt + nft)
459 ELSEIF (ity == 2) then! QUADS
460 ix => ixq(1:nixq, lft + nft:llt + nft)
461 ELSEIF (ity == 7) THEN ! TRIAS
462 ix => ixtg(1:nixtg, lft + nft:llt + nft)
463 ENDIF
464 nbmat = mat_param(ix(1, 1))%MULTIMAT%NB
465 !mixture density : in case of inigrav with (water+vapor) vapor eos must use mixture density
466 rho0=zero
467 DO imat=1,nbmat
468 matid = mat_param(ix(1, 1))%MULTIMAT%MID(imat)
469 volfrac= mat_param(ix(1, 1))%MULTIMAT%VFRAC(imat)
470 matlaw = ipm(2, matid)
471 rho0 = rho0+pm(1,matid)*volfrac
472 ENDDO
473
474 DO imat=1,nbmat
475 matid = mat_param(ix(1, 1))%MULTIMAT%MID(imat)
476 matlaw = ipm(2, matid)
477 CALL inigrav_eos(nel, nel2 , ng , matid, ipm, grav0, rho0, depth, pm, bufmat(iadbuf), elbuf_tab,
478 . psurf, list , pgrav, imat , mlw, npf,tf,nummat,mat_param)
479 ENDDO
480 ! GLOBAL STATE
481 DO kk=1,nel2
482 i = list(kk)
483 IF(nbmat==1)THEN
484 ! warning in this case GBUF => LBUF so dont reinit GBUF to 0.
485 !LBUF => ELBUF_TAB(NG)%BUFLY(IMAT)%LBUF(1, 1, 1)
486 !GBUF%RHO(I) = LBUF%RHO(I)
487 !GBUF%EINT(I) = LBUF%EINT(I)
488 gbuf%SIG(i + 0 * nel) = -pgrav(kk)
489 gbuf%SIG(i + 1 * nel) = -pgrav(kk)
490 gbuf%SIG(i + 2 * nel) = -pgrav(kk)
491 ELSE
492 gbuf%RHO(i) = zero
493 gbuf%EINT(i) = zero
494 DO imat=1,nbmat
495 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
496 vfrac = lbuf%VOL(i) / gbuf%VOL(i)
497 gbuf%EINT(i) = gbuf%EINT(i) + vfrac*lbuf%EINT(i)
498 gbuf%RHO(i) = gbuf%RHO(i) + vfrac*lbuf%RHO(i)
499 ENDDO
500 gbuf%SIG(i + 0 * nel) = -pgrav(kk)
501 gbuf%SIG(i + 1 * nel) = -pgrav(kk)
502 gbuf%SIG(i + 2 * nel) = -pgrav(kk)
503 ENDIF
504 ENDDO
505 !=======================================!
506 ! DEFAULT !
507 !=======================================!
508 CASE DEFAULT
509 IF (n2d == 0) THEN ! HEXA
510 ix => ixs(1:nixs, 1:numels)
511 nix = nixs
512 ELSEIF (ity == 2) then! QUADS
513 ix => ixq(1:nixq, 1:numelq)
514 nix = nixq
515 ELSEIF (ity == 7) then! TRIAS
516 ix => ixtg(1:nixtg, 1:numeltg)
517 nix = nixtg
518 ENDIF
519 CALL ancmsg(msgid = 89,
520 . anmode = aninfo,
521 . msgtype = msgwarning,
522 . i1 = uid,
523 . i2 = ipart(4,m(i)),
524 . i3 = mlw)
525 EXIT
526 END SELECT
527
528 ENDIF !(ITY == 1. .OR. ITY == 2) .AND.MLW /= 0)
529 ENDDO !next NG
530
531 IF(ALLOCATED(max_coor))DEALLOCATE(max_coor)
532 IF(ALLOCATED(min_coor))DEALLOCATE(min_coor)
533 IF(ALLOCATED(n_surf))DEALLOCATE(n_surf)
534 IF(ALLOCATED(zf))DEALLOCATE(zf)
535 IF(igrp/=0)THEN
536 IF(ALLOCATED(ip))DEALLOCATE(ip)
537 ENDIF
538
539 ENDDO !next K
540
541 RETURN
542 END
543
544
545
546
547!||====================================================================
548!|| check_is_on_segment ../starter/source/initial_conditions/inigrav/inigrav_load.F
549!||--- called by ------------------------------------------------------
550!|| inigrav_load ../starter/source/initial_conditions/inigrav/inigrav_load.F
551!||====================================================================
552 SUBROUTINE check_is_on_segment(P,Z,H,lTRIA,lTEST)
553C-----------------------------------------------
554C Description
555C-----------------------------------------------
556C Check if given point H is on quadrangle P(:,1),P(:,2),P(:,3),P(:,4)
557C using triangulation with midpoint Z(:)
558!
559! P(:,1)
560! +
561! /|\
562! / | \
563! / | \
564! / | \
565! / 1 | 4 \
566! / | \
567! P(:,2) +------Z+-----+B = P(:,4)
568! \ | /
569! \ 2 | 3 /
570! \ | /
571! \ | /
572! \ | /
573! \|/
574! +
575! A = P(:,3)
576
577!Il s'agit de savoir avec cette fonction s'il est a l'interieur
578! ou a l'exterieur.
579! ANCIEN CRITERE
580! 1 :interieur
581! 0 :exterieur
582!
583! Z+-----+B CRITERIA :
584! | H / ----------
585! | + / (ZA^ZH).(ZH^ZB) & (AZ^AH).(AH^AB) > 0
586! | / <=> / (za^u ).( u^zb) & (az^v ).( v^ab) > 0
587! | / | U=ZH
588! |/ \ V=AH
589! A+
590!
591! NOUVEAU CRITERE (bien meilleure performance CPU) :
592! Coordonnes de H dans le repere (ZA,ZB) doit etre
593! -TOL < x < 1+ TOL
594! -TOL < y < 1+ TOL
595
596C-----------------------------------------------
597C M o d u l e s
598C-----------------------------------------------
599C
600C-----------------------------------------------
601C I m p l i c i t T y p e s
602C-----------------------------------------------
603#include "implicit_f.inc"
604C-----------------------------------------------
605C G l o b a l P a r a m e t e r s
606C-----------------------------------------------
607C
608C-----------------------------------------------
609C C o m m o n B l o c k s
610C-----------------------------------------------
611C
612C-----------------------------------------------
613C D u m m y A r g u m e n t s
614C-----------------------------------------------
615 my_real,INTENT(IN) :: p(3,4),h(3),z(3)
616 LOGICAL,INTENT(IN) :: lTRIA
617 LOGICAL, INTENT(OUT) :: lTEST
618C-----------------------------------------------
619C L o c a l V a r i a b l e s
620C-----------------------------------------------
621 INTEGER ISONTRIA(4),IA,IB,NTRIA,I
622 my_real
623 . PS1,PS2,PS3,ZH(3),ZB(3),ZA(3),
624 . TOL,
625 . norm_za_2,norm_zb_2,
626 . coef,
627 . t, s, tolcrit
628C-----------------------------------------------
629C S o u r c e L i n e s
630C-----------------------------------------------
631 ntria=4
632 IF(ltria)ntria=1
633 isontria(1:4)=0
634 tolcrit = em06
635 tol = tolcrit
636 ia=0
637 ib=1
638
639 DO i=1,ntria
640 ia=ia+1
641 ib=ib+1
642 IF(ia>4)ia=1
643 IF(ib>4)ib=1
644 za(1)=p(1,ia)-z(1)
645 za(2)=p(2,ia)-z(2)
646 za(3)=p(3,ia)-z(3)
647 zb(1)=p(1,ib)-z(1)
648 zb(2)=p(2,ib)-z(2)
649 zb(3)=p(3,ib)-z(3)
650
651 norm_za_2 = (za(1)*za(1)+za(2)*za(2)+za(3)*za(3))
652 norm_zb_2 = (zb(1)*zb(1)+zb(2)*zb(2)+zb(3)*zb(3))
653
654 zh(1) = h(1) - z(1)
655 zh(2) = h(2) - z(2)
656 zh(3) = h(3) - z(3)
657
658 ps1 = za(1)*zh(1)+za(2)*zh(2)+za(3)*zh(3)
659
660 ps2 = zb(1)*zh(1)+zb(2)*zh(2)+zb(3)*zh(3)
661
662 ps3 = zb(1)*zh(1)+zb(2)*za(2)+zb(3)*za(3)
663
664 coef = one-ps3*ps3/norm_za_2/norm_zb_2
665 coef = one / coef
666
667 t = coef * ( ps2/norm_zb_2 - ps3*ps1/norm_za_2/norm_zb_2)
668 s = coef * ( -ps3*ps2/norm_za_2/norm_zb_2 + ps1/norm_za_2)
669
670 !IF(iflg_db == 1)THEN
671 ! print *, "coor ZA,ZB =", t,s
672 ! write(*,fmt='(A12,3F30.16)')"*createnode ",Z(1:3)
673 ! write(*,fmt='(A12,3F30.16)')"*createnode ",H(1:3)
674 !ENDIF
675
676 IF(t>=-tol)THEN
677 IF(s>=-tol )THEN
678 IF(s+t<=one+tol)THEN
679 isontria(i) = 1
680 ENDIF
681 ENDIF
682 ENDIF
683
684 ENDDO
685
686 ltest=.false.
687 IF(sum(isontria(1:ntria))>=1)ltest=.true.
688
689 END SUBROUTINE check_is_on_segment
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
inigrav_eos
subroutine inigrav_eos(nelg, nel, ng, matid, ipm, grav0, rho0, depth, pm, bufmat, elbuf_tab, psurf, list, pres, imat, mlw, npf, tf, nummat, mat_param)
Definition inigrav_eos.F:32
inigrav_load
subroutine inigrav_load(elbuf_tab, ipart, igrpart, iparg, iparttg, iparts, ipartq, x, ixs, ixq, ixtg, pm, ipm, bufmat, multi_fvm, ale_connectivity, nv46, igrsurf, itab, ebcs_tab, npf, tf, mat_param)
Definition inigrav_load.F:43
check_is_on_segment
subroutine check_is_on_segment(p, z, h, ltria, ltest)
Definition inigrav_load.F:553
inigrav_m37
subroutine inigrav_m37(nelg, nel, ng, matid, ipm, grav0, depth, pm, bufmat, elbuf_tab, psurf, list)
Definition inigrav_m37.F:30
inigrav_m51
subroutine inigrav_m51(nelg, nel, ng, matid, ipm, grav0, depth, pm, bufmat, elbuf_tab, psurf, list, ale_connectivity, nv46, ix, nix, nft, bufmatg, iparg)
Definition inigrav_m51.F:33
inigrav_part_list
subroutine inigrav_part_list(ipart, igrpart, ebcs_tab)
Definition inigrav_part_list.F:32
initia
subroutine initia(iparg, elbuf, ms, in, v, x, ixs, ixq, ixc, ixt, ixp, ixr, detonators, geo, pm, rby, npby, lpby, npc, npts, pld, veul, ale_connectivity, skew, fill, ipart, itab, sensors, skvol, ixtg, thk, nloc_dmg, group_param_tab, glob_therm, igrnod, igrsurf, bufsf, vr, bufmat, xlas, las, dtelem, mss, msq, msc, mst, msp, msr, mstg, ptg, inc, nod2eltg, knod2eltg, inp, inr, intg, index, itri, kxx, ixx, xelemwa, iwa, nod2elq, knod2elq, nod2els, knod2els, kxsp, ixsp, nod2sp, ispcond, icode, iskew, iskn, ispsym, xframe, isptag, spbuf, mssx, nsigi, npbyl, lpbyl, rbyl, msnf, mssf, nsigsh, igeo, ipm, nsigs, nsigsph, vns, vnsx, stc, stt, stp, str, sttg, stur, bns, bnsx, volnod, bvolnod, etnod, nshnod, stifint, fxbdep, fxbvit, fxbacc, fxbipm, fxbrpm, fxbelm, fxbsig, fxbmod, ins, ptshel, ptsh3n, ptsol, ptquad, wma, ptsph, fxbnod, mbufel, mdepl, fxani, numel, nsigrs, sh4tree, sh3tree, mcp, temp, imerge2, iadmerge2, slnrbm, nslnrbm, rmstifn, rmstifr, ms_layer, zi_layer, itag, itagel, mcpc, mcptg, xrefc, xreftg, xrefs, mssa, msrt, irbe2, lrbe2, inivol, kvol, nbsubmat, ixs10, ixs16, ixs20, totaddmas, ipmas, stifn, msz2, itagn, sitage, itage, ixr_kj, elbuf_tab, nom_opt, ptr_nopt_rbe2, ptr_nopt_adm, ptr_nopt_fun, sol2sph, irst, sh3trim, xfem_tab, kxig3d, ixig3d, msig3d, knot, nctrlmax, wige, stack, rnoise, drape, sh4ang, sh3ang, geo_stack, igeo_stack, stifintr, strc, strp, strr, strtg, perturb, itagnd, nativ_sms, iloadp, facload, ptspri, nsigbeam, ptbeam, nsigtruss, pttruss, multi_fvm, sigi, sigsh, sigsp, sigsph, sigrs, sigbeam, sigtruss, strsglob, straglob, orthoglob, isigsh, iyldini, ksigsh3, fail_ini, iusolyld, iuser, iddlevel, inimap1d, inimap2d, func2d, fvm_inivel, tagprt_sms, igrbric, igrquad, igrsh4n, igrsh3n, igrpart, totmas, knotlocpc, knotlocel, vnige, bnige, fxbglm, fxbcpm, fxbcps, fxblm, fxbfls, fxbdls, fxb_matrix, fxb_matrix_add, fxb_last_adress, ptr_nopt_fxb, r_skew, knod2el1d, nod2el1d, ebcs_tab, rby_iniaxis, alea, knod2elc, nod2elc, dr, slrbody, drapeg, ipari, intbuf_tab, interfaces, mat_param, npreload_a, preload_a, fail_fractal, fail_brokmann, defaults, ndamp_freq_range, dampr, ibeam_vector, rbeam_vector, ikine)
Definition initia.F:188
interp
subroutine interp(tf, tt, npoint, f, tg)
Definition interp.F:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
ale_connectivity_mod
Definition ale_connectivity_mod.F:223
groupdef_mod
Definition groupdef_mod.F:662
inigrav::inigrv
integer, dimension(:,:), allocatable inigrv
Definition inigrav_mod.F:38
message_mod
Definition message_mod.F:1249
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
starter
program starter
Definition starter.F:39
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:254