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
75 use element_mod , only : nixs,nixc,nixq,nixtg
76C-----------------------------------------------
77C I m p l i c i t T y p e s
78C-----------------------------------------------
79#include "implicit_f.inc"
80C-----------------------------------------------
81C G l o b a l P a r a m e t e r s
82C-----------------------------------------------
83#include "mvsiz_p.inc"
84C-----------------------------------------------
85C C o m m o n B l o c k s
86C-----------------------------------------------
87#include "vect01_c.inc"
88#include "param_c.inc"
89#include "com01_c.inc"
90#include "com04_c.inc"
91#include "scr17_c.inc"
92#include "tabsiz_c.inc"
93C-----------------------------------------------
94C D u m m y A r g u m e n t s
95C-----------------------------------------------
96 INTEGER, INTENT(IN),TARGET :: IPART(LIPART1,*)
97 INTEGER, INTENT(IN) :: IPARTS(NUMELS), IPARTQ(NUMELQ), IPARTTG(NUMELTG)
98 INTEGER, INTENT(IN) :: IPM(NPROPMI, NUMMAT),ITAB(NUMNOD)
99 INTEGER, INTENT(IN) :: IPARG(NPARG,NGROUP),NV46
100 INTEGER, INTENT(IN), TARGET :: IXS(NIXS, NUMELS), IXQ(NIXQ, NUMELQ), IXTG(NIXTG, NUMELTG)
101 TYPE(elbuf_struct_),DIMENSION(NGROUP),TARGET :: ELBUF_TAB
102 my_real, INTENT(IN) :: x(3,numnod), pm(npropm,nummat), bufmat(*)
103 TYPE(multi_fvm_struct), INTENT(IN) :: MULTI_FVM
104 TYPE (SURF_), DIMENSION(NSURF),INTENT(IN) :: IGRSURF
105 TYPE(t_ebcs_tab), INTENT(INOUT) :: EBCS_TAB
106C-----------------------------------------------
107 TYPE (GROUP_) , DIMENSION(NGRPART) :: IGRPART
108 TYPE(t_ale_connectivity), INTENT(INOUT) :: ALE_CONNECTIVITY
109 INTEGER,INTENT(IN)::NPF(SNPC)
110 my_real,INTENT(IN)::TF(STF)
111 TYPE(matparam_struct_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM
112C-----------------------------------------------
113C L o c a l V a r i a b l e s
114C-----------------------------------------------
115 my_real :: grav0,bx,by,bz,nx,ny,nz,dotprod,dotp,rho0,ngx,ngy,ngz,alpha,interp(3),vala,valb,norm
116 my_real :: z(3), depth(mvsiz), n(3,8),vec(3), delta_p(mvsiz),psurf,pgrav(mvsiz),vfrac,p(3,1:4)
117 my_real :: lambda,diag1(3),diag2(3),tol,dist,b(3),volfrac,dept
118 my_real,ALLOCATABLE,DIMENSION(:,:) :: min_coor, max_coor, n_surf, zf
119 INTEGER :: K,KK,IGRP,ISURF,IGRAV,NPART_IN_GROUP,MLW,NG,IAD0,II,NEL,J,IE,I,IERROR,IMAT,M(MVSIZ),Q,NSEG,ISEG
120 INTEGER :: USERID
121 INTEGER :: NOD(8), MATID, IADBUF,IFORM,NEL2,LIST(MVSIZ),NBMAT,MATLAW,NIX,UID,SURF_TYPE,IELTYP
122 INTEGER, ALLOCATABLE, DIMENSION(:) :: IP
123 TYPE(g_bufel_) ,POINTER :: GBUF
124 TYPE(l_bufel_), POINTER :: LBUF
125 LOGICAL :: lCYCLE,lERROR,lTRIA,lFOUND_PROJ,lTEST,l_PLANAR_SURF,l_USER_SURF
126 INTEGER, DIMENSION(:, :), POINTER :: IX
127 character*64 :: ERRMSG
128C-----------------------------------------------
129C S o u r c e L i n e s
130C-----------------------------------------------
131
132 CALL inigrav_part_list( ipart ,igrpart, ebcs_tab ) !list related parts to apply gravity loading with /ebcs/nrf
133
134 DO k=1,ninigrav
135 igrp = inigrv(1,k)
136 isurf = inigrv(2,k)
137 igrav = inigrv(3,k)
138 uid = inigrv(4,k)
139 bx = linigrav(01,k)
140 by = linigrav(02,k)
141 bz = linigrav(03,k)
142 nx = linigrav(04,k)
143 ny = linigrav(05,k)
144 nz = linigrav(06,k)
145 grav0 = linigrav(07,k)
146 ngx = linigrav(08,k)
147 ngy = linigrav(09,k)
148 ngz = linigrav(10,k)
149 psurf = linigrav(11,k)
150 IF (grav0==zero)cycle !nothing to do if g=0 at time=0 => Delta_P=0 => P=P(t=0)=P0
151 !READING GRPART DATA
152 IF(igrp/=0)THEN
153 npart_in_group = igrpart(igrp)%NENTITY
154 ALLOCATE(ip(npart_in_group),stat=ierror)
155 IF(ierror/=0)CALL ancmsg(msgid = 268,
156 . anmode = aninfo,
157 . msgtype = msgerror,
158 . c1 = 'INIGRAV')
159 ip(1:npart_in_group) = ipart(4,igrpart(igrp)%ENTITY(1:npart_in_group))
160 ELSE
161 !by default treat all parts
162 iad0 = 0
163 npart_in_group = 0
164 ENDIF
165 l_planar_surf = .false.
166 l_user_surf = .false.
167 surf_type=-1 !-1:unaffected
168 IF(isurf>0)surf_type=igrsurf(isurf)%TYPE
169 IF(surf_type==200 .OR. surf_type == -1)THEN
170 l_planar_surf=.true.
171 ELSEIF(surf_type == 0) THEN
172 l_user_surf=.true.
173 nseg = igrsurf(isurf)%NSEG
174 !STORING DATA FROM USER SURFACE (BOX,NORMAL,CENTROID)
175 ALLOCATE(max_coor(3,nseg),min_coor(3,nseg),n_surf(3,nseg),zf(3,nseg))
176 DO iseg=1,nseg
177 ltria=.false.
178 IF(igrsurf(isurf)%NODES(iseg,4)==0)ltria=.true.
179 IF(ltria)THEN
180 p(1:3,4)=p(1:3,3)
181 ENDIF
182 !face coordinates
183 p(1,1:4)=x(1,igrsurf(isurf)%NODES(iseg,1:4))
184 p(2,1:4)=x(2,igrsurf(isurf)%NODES(iseg,1:4))
185 p(3,1:4)=x(3,igrsurf(isurf)%NODES(iseg,1:4))
186 !box (to skip cell centroid if too far from it)
187 min_coor(1,iseg)=min(p(1,1),p(1,2),p(1,3),p(1,4))
188 min_coor(2,iseg)=min(p(2,1),p(2,2),p(2,3),p(2,4))
189 min_coor(3,iseg)=min(p(3,1),p(3,2),p(3,3),p(3,4))
190 max_coor(1,iseg)=max(p(1,1),p(1,2),p(1,3),p(1,4))
191 max_coor(2,iseg)=max(p(2,1),p(2,2),p(2,3),p(2,4))
192 max_coor(3,iseg)=max(p(3,1),p(3,2),p(3,3),p(3,4))
193 !face normal
194 diag1(1)=(p(1,3)-p(1,1))
195 diag1(2)=(p(2,3)-p(2,1))
196 diag1(3)=(p(3,3)-p(3,1))
197 diag2(1)=(p(1,4)-p(1,2))
198 diag2(2)=(p(2,4)-p(2,2))
199 diag2(3)=(p(3,4)-p(3,2))
200 n_surf(1,iseg)= diag1(2)*diag2(3)-diag1(3)*diag2(2)
201 n_surf(2,iseg)=-diag1(1)*diag2(3)+diag1(3)*diag2(1)
202 n_surf(3,iseg)= diag1(1)*diag2(2)-diag1(2)*diag2(1)
203 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))
204 n_surf(1,iseg)=n_surf(1,iseg)/norm
205 n_surf(2,iseg)=n_surf(2,iseg)/norm
206 n_surf(3,iseg)=n_surf(3,iseg)/norm
207 !face centroid
208 IF(ltria)THEN
209 zf(1,iseg)=third*sum(p(1,1:3))
210 zf(2,iseg)=third*sum(p(2,1:3))
211 zf(3,iseg)=third*sum(p(3,1:3))
212 ELSE
213 zf(1,iseg)=fourth*sum(p(1,1:4))
214 zf(2,iseg)=fourth*sum(p(2,1:4))
215 zf(3,iseg)=fourth*sum(p(3,1:4))
216 ENDIF
217 ENDDO
218 !define a tolerance to check if cell centroid too far from face (box)
219 iseg=1
220 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))
221 ENDIF
222
223
224 DO ng=1,ngroup
225 mlw = iparg(1,ng)
226 ity = iparg(5,ng)
227 nel = iparg(2,ng)
228 nft = iparg(3,ng)
229 iad = iparg(4,ng) - 1
230 gbuf => elbuf_tab(ng)%GBUF
231 pgrav(1:nel)=zero
232 IF(mlw==0)cycle
233 IF( (ity == 1. .AND.n2d==0) .OR. ((ity == 2. .OR. ity == 7).AND.n2d/=0) )THEN
234 IF(ity==1)THEN
235 imat = ixs(1,1+nft)
236 m(1:nel) = iparts(1+nft:nel+nft)
237 ELSEIF(ity==2)THEN
238 imat = ixq(1,1+nft)
239 m(1:nel) = ipartq(1+nft:nel+nft)
240 ELSEIF(ity==7)THEN
241 imat = ixtg(1,1+nft)
242 m(1:nel) = iparttg(1+nft:nel+nft)
243 ENDIF
244 !check if PART must be treated or not
245 list(1:mvsiz) = 0
246 nel2 = 0
247 !There may be several shares per group.All the elements have even mat and even prop.
248 ! Creation of a list list (1: nel 2) c'1: nel with the elements of the shares to be treated
249 !===================================================================!
250 ! ELEMENT SELECTION
251 !===================================================================!
252 IF(igrp/=0)THEN
253 lcycle = .true.
254 DO i=1,nel
255 DO q=1,npart_in_group
256 IF(ipart(4,m(i))==ip(q))THEN
257 nel2 = nel2+1
258 list(nel2) = i
259 lcycle = .false.
260 EXIT
261 ENDIF
262 ENDDO
263 enddo!next I
264 IF(lcycle)cycle
265 ELSE
266 nel2=nel
267 DO i=1,nel
268 list(i) = i !(/1:NEL/) !identite
269 ENDDO
270 ENDIF
271 ! LIST is the sublist of 1:NEL of elements whose PART is to be processed with INIGRAV
272
273 ! This group is going to be treated since PART_ID in GRPART_ID
274 !===================================================================!
275 ! CENTROID CALCULATION
276 !===================================================================!
277 depth(:) = ep20
278 DO kk=1,nel2
279 i = list(kk)
280 ie = i + nft
281 userid = 0
282 IF(n2d==0)THEN
283 n(1:3,1) = x(1:3,ixs(2,ie))
284 n(1:3,2) = x(1:3,ixs(3,ie))
285 n(1:3,3) = x(1:3,ixs(4,ie))
286 n(1:3,4) = x(1:3,ixs(5,ie))
287 n(1:3,5) = x(1:3,ixs(6,ie))
288 n(1:3,6) = x(1:3,ixs(7,ie))
289 n(1:3,7) = x(1:3,ixs(8,ie))
290 n(1:3,8) = x(1:3,ixs(9,ie))
291 z(1) = one_over_8*(sum(n(1,1:8)))
292 z(2) = one_over_8*(sum(n(2,1:8)))
293 z(3) = one_over_8*(sum(n(3,1:8)))
294 userid = ixs(nixs,ie)
295 ELSEIF(ity==2)THEN
296 n(1:3,1) = x(1:3,ixq(2,ie))
297 n(1:3,2) = x(1:3,ixq(3,ie))
298 n(1:3,3) = x(1:3,ixq(4,ie))
299 n(1:3,4) = x(1:3,ixq(5,ie))
300 z(1) = fourth*(sum(n(1,1:4)))
301 z(2) = fourth*(sum(n(2,1:4)))
302 z(3) = fourth*(sum(n(3,1:4)))
303 userid = ixq(nixq,ie)
304 ELSEIF(ity==7)THEN
305 n(1:3,1) = x(1:3,ixtg(2,ie))
306 n(1:3,2) = x(1:3,ixtg(3,ie))
307 n(1:3,3) = x(1:3,ixtg(4,ie))
308 n(1:3,4) = zero
309 z(1) = third*(sum(n(1,1:3)))
310 z(2) = third*(sum(n(2,1:3)))
311 z(3) = third*(sum(n(3,1:3)))
312 userid = ixtg(nixtg,ie)
313 ENDIF
314 !===================================================================!
315 IF(l_planar_surf)THEN
316 !===================================================================!
317 ! --- DISTANCE CALCULATION ---
318 !-PLANAR SURFACE : IGRSURF(ISURF)%TYPE = 200 OR ISURF=0 !
319 !===================================================================!
320 !---NORMAL : NX,NY,NZ
321 !---BASIS POINT : BX,BY,BZ
322 !intersection point such as INTERP = Z + ALPHA . NG et <B-INTERP,N> = 0, donc :
323 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
324 interp(1) = z(1) + alpha*ngx
325 interp(2) = z(2) + alpha*ngy
326 interp(3) = z(3) + alpha*ngz
327 ! Depth = <perf_P, n>
328 !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))
329 !DEPTH(I) = SQRT(DEPTH(I))
330 depth(kk)= (interp(1)-z(1))*ngx+(interp(2)-z(2))*ngy+(interp(3)-z(3))*ngz
331 depth(kk) = -depth(kk)
332 lerror=.false.
333 lfound_proj=.true.
334 ELSEIF(l_user_surf) THEN
335 !===================================================================!
336 ! --- DISTANCE CALCULATION ---
337 !-USER SURFACE : IGRSURF(ISURF)%TYPE=0 , ISURF/=0 !
338 !===================================================================!
339 !LOOP ON CELL, THEN LOOP ON FACE : complexity is O(NCELL*NFACE), can be optimized with space partitioning
340 lerror=.false.
341 lfound_proj=.false.
342 DO iseg=1,nseg
343
344 IF(ngx==zero)THEN
345 IF(z(1)+tol < min_coor(1,iseg) .OR. z(1)-tol > max_coor(1,iseg))cycle
346 ENDIF
347 IF(ngy==zero)THEN
348 IF(z(2)+tol < min_coor(2,iseg) .OR. z(2)-tol > max_coor(2,iseg))cycle
349 ENDIF
350 IF(ngz==zero)THEN
351 IF(z(3)+tol < min_coor(3,iseg) .OR. z(3)-tol > max_coor(3,iseg))cycle
352 ENDIF
353 dotp = n_surf(1,iseg)*ngx + n_surf(2,iseg)*ngy + n_surf(3,iseg)*ngz
354 IF(abs(dotp)<=em04)THEN
355 CALL ancmsg(msgid = 90,
356 . anmode = aninfo,
357 . msgtype = msgerror,
358 . i1=uid,
359 . c1="INPUT SURFACE HAS INCOMPATIBLE SLOPE WITH GRAVITY DIRECTION"
360 . )
361 lerror=.true.
362 EXIT
363 ENDIF
364 p(1,1:4)=x(1,igrsurf(isurf)%NODES(iseg,1:4))
365 p(2,1:4)=x(2,igrsurf(isurf)%NODES(iseg,1:4))
366 p(3,1:4)=x(3,igrsurf(isurf)%NODES(iseg,1:4))
367 ! Let INTERP the normal projection point on plane generated by current face ISEG
368 ! INTERP = Z + LAMBDA*NG (colinear to gravity direction , NG gravity vector, Z cell centroid, LAMBDA unknown scalar)
369 ! <INTERP-B,N_SURF> = 0 (must satisfying plane equation, B basis point, Z cell centroid)
370 ! => INTERP = ...
371 b(1:3)=zf(1:3,iseg)
372 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))
373 IF(dist<=tol*tol)THEN
374 !use another basis point if face centroid superimposed with cell centroid
375 b(1:3)=p(1:3,1)
376 ENDIF
377 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)
378 lambda = lambda / dotp
379 interp(1) = z(1)+lambda*ngx
380 interp(2) = z(2)+lambda*ngy
381 interp(3) = z(3)+lambda*ngz
382 !CHECK INTERP IS INSIDE FACE, OTHERWISE SKIP IT
383 CALL check_is_on_segment(p(1:3,1:4),b,interp,ltria,ltest)
384 IF(.NOT.ltest)cycle
385 dept = -((interp(1)-z(1))*ngx+(interp(2)-z(2))*ngy+(interp(3)-z(3))*ngz)
386 IF(abs(dept) < abs(depth(kk)))THEN
387 depth(kk) = dept
388 ENDIF
389 lfound_proj=.true.
390 ENDDO
391 IF(lerror)EXIT
392 IF(.NOT.lfound_proj)THEN
393 errmsg = "UNABLE TO PROJECT ON SURFACE CENTROID FROM CELL ID= "
394 WRITE(errmsg(52:62),fmt='(I10)')userid
395 CALL ancmsg(msgid = 90,
396 . anmode = aninfo,
397 . msgtype = msgerror,
398 . i1=uid,
399 . c1=errmsg(1:62)
400 . )
401 EXIT
402 ENDIF
403 ENDIF !l_PLANAR_SURF elif l_USER_SURF
404 ENDDO !next I
405 IF(lerror)EXIT
406 IF(.NOT.lfound_proj)EXIT
407 !===================================================================!
408 ! INITIALISATION DE L'ETAT INITIAL !
409 !===================================================================!
410 IF (n2d == 0) THEN
411 matid = ixs(1, 1 + nft)
412 ELSEIF(ity==2)THEN
413 matid = ixq(1, 1 + nft)
414 ELSEIF(ity==7)THEN
415 matid = ixtg(1, 1 + nft)
416 ENDIF
417 iadbuf = max(1,ipm(7, matid)) ! Adress of data for the material in the material buffer
418 !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
419 SELECT CASE (mlw)
420 !=======================================!
421 ! 3,4,6,49 (ANY EOS) !
422 !=======================================!
423 CASE (3, 4, 6, 10, 49) !Material law defined with Equations Of States
424 rho0 = pm(1,matid)
425 CALL inigrav_eos(nel,nel2, ng, matid, ipm, grav0, rho0, depth, pm, bufmat(iadbuf),
426 . elbuf_tab, psurf, list , pgrav, 0 ,mlw, npf, tf ,nummat,mat_param)
427 !=======================================!
428 ! MM-ALE 37 !
429 !=======================================!
430 CASE (37)
431 CALL inigrav_m37(nel,nel2, ng, matid, ipm, grav0, depth, pm, bufmat(iadbuf), elbuf_tab, psurf, list)
432 !=======================================!
433 ! MM-ALE 51 !
434 !=======================================!
435 CASE (51)
436 IF (n2d == 0) THEN ! HEXA
437 ix => ixs(1:nixs, 1:numels)
438 nix = nixs
439 ELSEIF (ity == 2) then! QUADS
440 ix => ixq(1:nixq, 1:numelq)
441 nix = nixq
442 ELSEIF (ity == 7) then! TRIAS
443 ix => ixtg(1:nixtg, 1:numeltg)
444 nix = nixtg
445 ENDIF
446 iform = bufmat(iadbuf + 31 - 1)
447 IF (iform /= 3) THEN
448 CALL inigrav_m51(nel , nel2, ng, matid, ipm, grav0, depth, pm, bufmat(iadbuf), elbuf_tab,
449 . psurf, list, ale_connectivity, nv46, ix , nix , nft , bufmat, iparg)
450 ELSE
451 CALL ancmsg(msgid = 84,
452 . anmode = aninfo,
453 . msgtype = msgwarning
454 . )
455 ENDIF
456 !=======================================!
457 ! MULTI-FLUID 151 !
458 !=======================================!
459 CASE(151)
460
461 IF (n2d == 0) THEN ! HEXA
462 ix => ixs(1:nixs, lft + nft:llt + nft)
463 ELSEIF (ity == 2) then! QUADS
464 ix => ixq(1:nixq, lft + nft:llt + nft)
465 ELSEIF (ity == 7) THEN ! trias
466 ix => ixtg(1:nixtg, lft + nft:llt + nft)
467 ENDIF
468 nbmat = mat_param(ix(1, 1))%MULTIMAT%NB
469 !mixture density : in case of inigrav with (water+vapor) vapor eos must use mixture density
470 rho0=zero
471 DO imat=1,nbmat
472 matid = mat_param(ix(1, 1))%MULTIMAT%MID(imat)
473 volfrac= mat_param(ix(1, 1))%MULTIMAT%VFRAC(imat)
474 matlaw = ipm(2, matid)
475 rho0 = rho0+pm(1,matid)*volfrac
476 ENDDO
477
478 DO imat=1,nbmat
479 matid = mat_param(ix(1, 1))%MULTIMAT%MID(imat)
480 matlaw = ipm(2, matid)
481 CALL inigrav_eos(nel, nel2 , ng , matid, ipm, grav0, rho0, depth, pm, bufmat(iadbuf), elbuf_tab,
482 . psurf, list , pgrav, imat , mlw, npf,tf,nummat,mat_param)
483 ENDDO
484 ! GLOBAL STATE
485 DO kk=1,nel2
486 i = list(kk)
487 IF(nbmat==1)THEN
488 ! warning in this case GBUF => LBUF so dont reinit GBUF to 0.
489 !LBUF => ELBUF_TAB(NG)%BUFLY(IMAT)%LBUF(1, 1, 1)
490 !GBUF%RHO(I) = LBUF%RHO(I)
491 !GBUF%EINT(I) = LBUF%EINT(I)
492 gbuf%SIG(i + 0 * nel) = -pgrav(kk)
493 gbuf%SIG(i + 1 * nel) = -pgrav(kk)
494 gbuf%SIG(i + 2 * nel) = -pgrav(kk)
495 ELSE
496 gbuf%RHO(i) = zero
497 gbuf%EINT(i) = zero
498 DO imat=1,nbmat
499 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
500 vfrac = lbuf%VOL(i) / gbuf%VOL(i)
501 gbuf%EINT(i) = gbuf%EINT(i) + vfrac*lbuf%EINT(i)
502 gbuf%RHO(i) = gbuf%RHO(i) + vfrac*lbuf%RHO(i)
503 ENDDO
504 gbuf%SIG(i + 0 * nel) = -pgrav(kk)
505 gbuf%SIG(i + 1 * nel) = -pgrav(kk)
506 gbuf%SIG(i + 2 * nel) = -pgrav(kk)
507 ENDIF
508 ENDDO
509 !=======================================!
510 ! DEFAULT !
511 !=======================================!
512 CASE DEFAULT
513 IF (n2d == 0) THEN ! HEXA
514 ix => ixs(1:nixs, 1:numels)
515 nix = nixs
516 ELSEIF (ity == 2) then! QUADS
517 ix => ixq(1:nixq, 1:numelq)
518 nix = nixq
519 ELSEIF (ity == 7) then! TRIAS
520 ix => ixtg(1:nixtg, 1:numeltg)
521 nix = nixtg
522 ENDIF
523 CALL ancmsg(msgid = 89,
524 . anmode = aninfo,
525 . msgtype = msgwarning,
526 . i1 = uid,
527 . i2 = ipart(4,m(i)),
528 . i3 = mlw)
529 EXIT
530 END SELECT
531
532 ENDIF !(ITY == 1. .OR. ITY == 2) .AND.MLW /= 0)
533 ENDDO !next NG
534
535 IF(ALLOCATED(max_coor))DEALLOCATE(max_coor)
536 IF(ALLOCATED(min_coor))DEALLOCATE(min_coor)
537 IF(ALLOCATED(n_surf))DEALLOCATE(n_surf)
538 IF(ALLOCATED(zf))DEALLOCATE(zf)
539 IF(igrp/=0)THEN
540 IF(ALLOCATED(ip))DEALLOCATE(ip)
541 ENDIF
542
543 ENDDO !next K
544
545 RETURN
546 END
547
548
549
550
551!||====================================================================
552!|| check_is_on_segment ../starter/source/initial_conditions/inigrav/inigrav_load.F
553!||--- called by ------------------------------------------------------
554!|| inigrav_load ../starter/source/initial_conditions/inigrav/inigrav_load.F
555!||====================================================================
556 SUBROUTINE check_is_on_segment(P,Z,H,lTRIA,lTEST)
557C-----------------------------------------------
558C Description
559C-----------------------------------------------
560C Check if given point H is on quadrangle P(:,1),P(:,2),P(:,3),P(:,4)
561C using triangulation with midpoint Z(:)
562!
563! P(:,1)
564! +
565! /|\
566! / | \
567! / | \
568! / | \
569! / 1 | 4 \
570! / | \
571! P(:,2) +------Z+-----+B = P(:,4)
572! \ | /
573! \ 2 | 3 /
574! \ | /
575! \ | /
576! \ | /
577! \|/
578! +
579! A = P(:,3)
580
581!This function determines if it is inside
582! or outside.
583! ANCIEN CRITERE
584! 1 :interieur
585! 0 :exterieur
586!
587! Z+-----+B CRITERIA :
588! | H / ----------
589! | + / (za^zh).(zh^zb) & (az^ah).(ah^ab) > 0
590! | / <=> / (ZA^U ).( U^ZB) & (AZ^V ).( V^AB) > 0
591! | / | U=ZH
592! |/ \ V=AH
593! A+
594!
595! New criterion (much better CPU performance):
596! Coordinates of H in the (ZA,ZB) frame must be
597! -TOL < x < 1+ TOL
598! -TOL < y < 1+ TOL
599
600C-----------------------------------------------
601C M o d u l e s
602C-----------------------------------------------
603C
604C-----------------------------------------------
605C I m p l i c i t T y p e s
606C-----------------------------------------------
607#include "implicit_f.inc"
608C-----------------------------------------------
609C G l o b a l P a r a m e t e r s
610C-----------------------------------------------
611C
612C-----------------------------------------------
613C C o m m o n B l o c k s
614C-----------------------------------------------
615C
616C-----------------------------------------------
617C D u m m y A r g u m e n t s
618C-----------------------------------------------
619 my_real,INTENT(IN) :: p(3,4),h(3),z(3)
620 LOGICAL,INTENT(IN) :: lTRIA
621 LOGICAL, INTENT(OUT) :: lTEST
622C-----------------------------------------------
623C L o c a l V a r i a b l e s
624C-----------------------------------------------
625 INTEGER ISONTRIA(4),IA,IB,NTRIA,I
626 my_real
627 . PS1,PS2,PS3,ZH(3),ZB(3),ZA(3),
628 . TOL,
629 . norm_za_2,norm_zb_2,
630 . coef,
631 . t, s, tolcrit
632C-----------------------------------------------
633C S o u r c e L i n e s
634C-----------------------------------------------
635 ntria=4
636 IF(ltria)ntria=1
637 isontria(1:4)=0
638 tolcrit = em06
639 tol = tolcrit
640 ia=0
641 ib=1
642 zh(1) = h(1) - z(1)
643 zh(2) = h(2) - z(2)
644 zh(3) = h(3) - z(3)
645
646 DO i=1,ntria
647 ia=ia+1
648 ib=ib+1
649 IF(ia>4)ia=1
650 IF(ib>4)ib=1
651 za(1)=p(1,ia)-z(1)
652 za(2)=p(2,ia)-z(2)
653 za(3)=p(3,ia)-z(3)
654 zb(1)=p(1,ib)-z(1)
655 zb(2)=p(2,ib)-z(2)
656 zb(3)=p(3,ib)-z(3)
657
658 norm_za_2 = (za(1)*za(1)+za(2)*za(2)+za(3)*za(3))
659 norm_zb_2 = (zb(1)*zb(1)+zb(2)*zb(2)+zb(3)*zb(3))
660
661 ps1 = za(1)*zh(1)+za(2)*zh(2)+za(3)*zh(3)
662 ps2 = zb(1)*zh(1)+zb(2)*zh(2)+zb(3)*zh(3)
663 ps3 = zb(1)*za(1)+zb(2)*za(2)+zb(3)*za(3)
664
665 coef = one-ps3*ps3/norm_za_2/norm_zb_2
666 coef = one / coef
667
668 t = coef * ( ps2/norm_zb_2 - ps3*ps1/norm_za_2/norm_zb_2)
669 s = coef * ( -ps3*ps2/norm_za_2/norm_zb_2 + ps1/norm_za_2)
670
671 !IF(IFLG_DB == 1)THEN
672 ! print *, "coor ZA,ZB =", t,s
673 ! write(*,fmt='(A12,3F30.16)')"*createnode ",Z(1:3)
674 ! write(*,fmt='(A12,3F30.16)')"*createnode ",H(1:3)
675 !ENDIF
676
677 IF(t>=-tol)THEN
678 IF(s>=-tol )THEN
679 IF(s+t<=one+tol)THEN
680 isontria(i) = 1
681 ENDIF
682 ENDIF
683 ENDIF
684
685 ENDDO
686
687 ltest = .false.
688 IF(sum(isontria(1:ntria)) >= 1) ltest = .true.
689
690 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:557
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
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:225
groupdef_mod
Definition groupdef_mod.F:668
inigrav::inigrv
integer, dimension(:,:), allocatable inigrv
Definition inigrav_mod.F:38
message_mod
Definition message_mod.F:1257
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:895
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:256