OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
alefvm_sfint3.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!|| alefvm_sfint3 ../engine/source/ale/alefvm/alefvm_sfint3.F
25!||--- called by ------------------------------------------------------
26!|| alefvm_main ../engine/source/ale/alefvm/alefvm_main.F
27!||--- uses -----------------------------------------------------
28!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
29!|| alefvm_mod ../common_source/modules/ale/alefvm_mod.F
30!|| element_mod ../common_source/modules/elements/element_mod.F90
31!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
32!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
33!||====================================================================
34 SUBROUTINE alefvm_sfint3 (
35 1 IXS, NV46 , ALE_CONNECT, IALEFVM_FLG,
36 2 IPM, IPARG, NG ,
37 3 X , IAD22, NEL )
38C-----------------------------------------------
39C D e s c r i p t i o n
40C-----------------------------------------------
41C 'alefvm' is related to a collocated scheme (built from FVM and based on Godunov scheme)
42C which was temporarily introduced for experimental option /INTER/TYPE22 (FSI coupling with cut cell method)
43C This cut cell method is not completed, abandoned, and is not an official option.
44C There is no other use for this scheme which is automatically enabled when /INTER/TYPE22 is defined (INT22>0 => IALEFVM=1).
45C
46C This subroutine is treating an uncut cell.
47C-----------------------------------------------
48C M o d u l e s
49C-----------------------------------------------
50 USE alefvm_mod
51 USE i22bufbric_mod !debug
52 USE i22tri_mod
54 use element_mod , only : nixs
55C-----------------------------------------------
56C I m p l i c i t T y p e s
57C-----------------------------------------------
58#include "implicit_f.inc"
59C-----------------------------------------------
60C G l o b a l P a r a m e t e r s
61C-----------------------------------------------
62#include "mvsiz_p.inc"
63C-----------------------------------------------
64C C o m m o n B l o c k s
65C-----------------------------------------------
66#include "vect01_c.inc"
67#include "inter22.inc"
68#include "param_c.inc"
69C-----------------------------------------------
70C D e s c r i p t i o n
71C-----------------------------------------------
72C This subroutines computes internal forces for
73C finite volume scheme (IALEFVM==1)
74C
75C If option is not detected in input file then
76C subroutine is unplugged
77C-----------------------------------------------
78C D u m m y A r g u m e n t s
79C-----------------------------------------------
80 INTEGER NEL
81 INTEGER :: IXS(NIXS,*),NV46,IALEFVM_FLG, IPM(NPROPMI,*),IPARG(NPARG,*),NG
82 my_real :: X(3,*)
83 my_real :: iad22(*)
84 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
85C-----------------------------------------------
86C L o c a l V a r i a b l e s
87C-----------------------------------------------
88 INTEGER :: IAD2, IAD3
89 INTEGER :: I, II, IV, J, JV, IMAT, ILAW,IFLG_ALE,IFLG_EUL
90 INTEGER :: NIN, IB, IPRES_MOM
91 INTEGER :: NC1,NC2,NC3,NC4,NC5,NC6,NC7,NC8
92 my_real :: f0(3,mvsiz), fface(3,nv46,mvsiz)
93 my_real :: nx(6,mvsiz), ny(6,mvsiz), nz(6,mvsiz), p1,p2,denom,pf
94 my_real :: theta, m1, m2, mf
95 my_real ::
96 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
97 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz),
98 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz), z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz),
99 . z_1, z_2, u1n1, u2n1
100
101C-----------------------------------------------
102C P r e - C o n d i t i o n s
103C-----------------------------------------------
104 IF(alefvm_param%IEnabled==0) RETURN
105 imat = ixs(1,1+nft)
106 ilaw = ipm(2,imat)
107 ialefvm_flg = ipm(251,imat)
108 IF(ialefvm_flg <= 1)RETURN
109 IF(ilaw==11) RETURN
110 iflg_ale = iparg(7,ng)
111 iflg_eul = iparg(11,ng)
112C-----------------------------------------------
113C S o u r c e L i n e s
114C-----------------------------------------------
115
116 nin = 1
117
118 !-------------------------------------------------------------!
119 ! NORMAL VECTOR FOR ALE !
120 !-------------------------------------------------------------!
121! IF(jale==1 .OR. jeul==0)THEN !lagrangian bricks also
122 DO i=1,nel
123 ii = i + nft
124 !---8 local node numbers NC1 TO NC8 for solid element I ---!
125 nc1=ixs(2,ii)
126 nc2=ixs(3,ii)
127 nc3=ixs(4,ii)
128 nc4=ixs(5,ii)
129 nc5=ixs(6,ii)
130 nc6=ixs(7,ii)
131 nc7=ixs(8,ii)
132 nc8=ixs(9,ii)
133 !
134 !---Coordinates of the 8 nodes
135 x1(i)=x(1,nc1)
136 y1(i)=x(2,nc1)
137 z1(i)=x(3,nc1)
138 !
139 x2(i)=x(1,nc2)
140 y2(i)=x(2,nc2)
141 z2(i)=x(3,nc2)
142 !
143 x3(i)=x(1,nc3)
144 y3(i)=x(2,nc3)
145 z3(i)=x(3,nc3)
146 !
147 x4(i)=x(1,nc4)
148 y4(i)=x(2,nc4)
149 z4(i)=x(3,nc4)
150 !
151 x5(i)=x(1,nc5)
152 y5(i)=x(2,nc5)
153 z5(i)=x(3,nc5)
154 !
155 x6(i)=x(1,nc6)
156 y6(i)=x(2,nc6)
157 z6(i)=x(3,nc6)
158 !
159 x7(i)=x(1,nc7)
160 y7(i)=x(2,nc7)
161 z7(i)=x(3,nc7)
162 !
163 x8(i)=x(1,nc8)
164 y8(i)=x(2,nc8)
165 z8(i)=x(3,nc8)
166 ENDDO
167 DO i=1,nel
168 ! Face-1
169 nx(1,i)=(y3(i)-y1(i))*(z2(i)-z4(i)) - (z3(i)-z1(i))*(y2(i)-y4(i))
170 ny(1,i)=(z3(i)-z1(i))*(x2(i)-x4(i)) - (x3(i)-x1(i))*(z2(i)-z4(i))
171 nz(1,i)=(x3(i)-x1(i))*(y2(i)-y4(i)) - (y3(i)-y1(i))*(x2(i)-x4(i))
172 ! Face-2
173 nx(2,i)=(y7(i)-y4(i))*(z3(i)-z8(i)) - (z7(i)-z4(i))*(y3(i)-y8(i))
174 ny(2,i)=(z7(i)-z4(i))*(x3(i)-x8(i)) - (x7(i)-x4(i))*(z3(i)-z8(i))
175 nz(2,i)=(x7(i)-x4(i))*(y3(i)-y8(i)) - (y7(i)-y4(i))*(x3(i)-x8(i))
176 ! Face-3
177 nx(3,i)=(y6(i)-y8(i))*(z7(i)-z5(i)) - (z6(i)-z8(i))*(y7(i)-y5(i))
178 ny(3,i)=(z6(i)-z8(i))*(x7(i)-x5(i)) - (x6(i)-x8(i))*(z7(i)-z5(i))
179 nz(3,i)=(x6(i)-x8(i))*(y7(i)-y5(i)) - (y6(i)-y8(i))*(x7(i)-x5(i))
180 ! Face-4
181 nx(4,i)=(y2(i)-y5(i))*(z6(i)-z1(i)) - (z2(i)-z5(i))*(y6(i)-y1(i))
182 ny(4,i)=(z2(i)-z5(i))*(x6(i)-x1(i)) - (x2(i)-x5(i))*(z6(i)-z1(i))
183 nz(4,i)=(x2(i)-x5(i))*(y6(i)-y1(i)) - (y2(i)-y5(i))*(x6(i)-x1(i))
184 ! Face-5
185 nx(5,i)=(y7(i)-y2(i))*(z6(i)-z3(i)) - (z7(i)-z2(i))*(y6(i)-y3(i))
186 ny(5,i)=(z7(i)-z2(i))*(x6(i)-x3(i)) - (x7(i)-x2(i))*(z6(i)-z3(i))
187 nz(5,i)=(x7(i)-x2(i))*(y6(i)-y3(i)) - (y7(i)-y2(i))*(x6(i)-x3(i))
188 ! Face-6
189 nx(6,i)=(y8(i)-y1(i))*(z4(i)-z5(i)) - (z8(i)-z1(i))*(y4(i)-y5(i))
190 ny(6,i)=(z8(i)-z1(i))*(x4(i)-x5(i)) - (x8(i)-x1(i))*(z4(i)-z5(i))
191 nz(6,i)=(x8(i)-x1(i))*(y4(i)-y5(i)) - (y8(i)-y1(i))*(x4(i)-x5(i))
192 ENDDO
193! ENDIF
194
195 !-------------------------------------------------------------!
196 ! Assembling Forces !
197 !-------------------------------------------------------------!
198
199 !IPRES_MOM = 0,1,2,3,4 : (P1+P2)/2
200 !IPRES_MOM = 5 : (rho1c1P1+rho2c2P2)/(rho1c1+rho2c2) + (rho1c1*rho2c2/(rho1c1+rho2c2)) <uel-uadj,nel> acoustic solver
201
202 ipres_mom = alefvm_param%ISOLVER
203
204! print *, "IPRES_MOM, SFINT3", IPRES_MOM
205
206 SELECT CASE (ipres_mom)
207 CASE(5) ! acoustic solver
208 DO i=1,nel
209 ii = i + nft
210 iad2 = ale_connect%ee_connect%iad_connect(ii)
211 p1 = alefvm_buffer%F_FACE(1,4,ii)
212 m1 = alefvm_buffer%F_FACE(1,5,ii)
213 z_1 = alefvm_buffer%F_FACE(1,3,ii) ! impedance
214 DO j=1,nv46
215 iv = ale_connect%ee_connect%connected(iad2 + j - 1)
216 IF(iv > 0)THEN
217 !--ADJACENT ELEM DOES EXIST
218 iad3 = ale_connect%ee_connect%iad_connect(iv)
219 DO jv=1,nv46
220 IF(ale_connect%ee_connect%connected(iad3 + jv - 1)==ii)EXIT
221 ENDDO
222 u1n1 = + alefvm_buffer%F_FACE(3, j,ii)
223 u2n1 = - alefvm_buffer%F_FACE(3,jv,iv) !UN2 = <U2,N2> = <U2,-N1>
224 z_2 = alefvm_buffer%F_FACE(1, 3,iv)
225 p2 = alefvm_buffer%F_FACE(1, 4,iv)
226 m2 = alefvm_buffer%F_FACE(1, 5,iv)
227 denom = z_1 + z_2
228 mf = min(m1,m2)
229 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
230 pf = (z_1*p2 + z_2*p1)/denom + theta*(z_1*z_2*(u1n1-u2n1)/denom)
231 ELSE
232 !--ADJACENT ELEM DOES NOT EXIST
233 !SLIDING RWALL BC
234 u1n1 = + alefvm_buffer%F_FACE(3, j,ii)
235 !Pf = P1
236 mf = m1
237 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
238 pf = p1 + theta*half*z_1*u1n1
239 !Pf = P1
240 ENDIF
241 fface(1,j,i) = -half*pf*nx(j,i) !N= 2.SURF.n , where abs(n)=1 => abs(N) = 2S , this lines leads to FFACE(1J,I) = Pf.S.n
242 fface(2,j,i) = -half*pf*ny(j,i) ! negative sign because int(div.sigma,DV) => int(-P,dS)
243 fface(3,j,i) = -half*pf*nz(j,i)
244 enddo!next J
245 enddo!next I
246 CASE DEFAULT
247 DO i=1,nel
248 ii = i + nft
249 iad2 = ale_connect%ee_connect%iad_connect(ii)
250 p1 = alefvm_buffer%F_FACE(1,4,ii)
251 DO j=1,nv46
252 iv = ale_connect%ee_connect%connected(iad2 + j - 1)
253 IF(iv > 0)THEN
254 !--ADJACENT ELEM DOES EXIST
255 iad3 = ale_connect%ee_connect%iad_connect(iv)
256 DO jv=1,nv46
257 IF(ale_connect%ee_connect%connected(iad3 + jv - 1)==ii)EXIT
258 ENDDO
259 p2 = alefvm_buffer%F_FACE(1,4,iv)
260 pf = half*(p1+p2)
261 ELSE
262 !--ADJACENT ELEM DOES NOT EXIST
263 !SLIDING RWALL BC
264 pf = p1
265 ENDIF
266 fface(1,j,i) = -half*pf*nx(j,i) ! negative sign because int(div.sigma,DV) => int(-P,dS)
267 fface(2,j,i) = -half*pf*ny(j,i)
268 fface(3,j,i) = -half*pf*nz(j,i)
269 !print *, "P1,P2,Pf", P1,P2,Pf
270 !print *, " normal ", NX(J,I), NY(J,I), NZ(J,I)
271c write ( *, FMT='(A,I10,A,I3,A,3F20.12,A,F20.12)')
272c ." id=", ixs(11,ii)," ncycle=", NCYCLE, " F=", FFACE(1:3,J,I) , "Pf=", Pf
273 enddo!next J
274 enddo!next I
275 END SELECT
276
277 !-------------------------------------------------------------!
278 ! INTEGRAL ON EACH FACE from Integral(DIV(SIGMA),Volume) !
279 !-------------------------------------------------------------!
280 DO i=1,nel
281 ii = i + nft
282 IF(int22/=0)THEN
283 ib = nint(iad22(i))
284 IF(ib>0) cycle
285 ENDIF
286 f0(1,i) = sum(fface(1,1:nv46,i))
287 f0(2,i) = sum(fface(2,1:nv46,i))
288 f0(3,i) = sum(fface(3,1:nv46,i))
289
290 ! write (*,FMT='(A,I6,A,3F30.16)') " brick ID=", ixs(11,II)," Fint=", F0(1:3,I)
291
292
293 alefvm_buffer%FINT_CELL(1,ii) = f0(1,i)
294 alefvm_buffer%FINT_CELL(2,ii) = f0(2,i)
295 alefvm_buffer%FINT_CELL(3,ii) = f0(3,i)
296
297 !Fcell already contains gravity force, so do not erase but add
298 alefvm_buffer%FCELL(1,ii) = alefvm_buffer%FCELL(1,ii) + f0(1,i) + alefvm_buffer%FEXT_CELL(1,ii)
299 alefvm_buffer%FCELL(2,ii) = alefvm_buffer%FCELL(2,ii) + f0(2,i) + alefvm_buffer%FEXT_CELL(2,ii)
300 alefvm_buffer%FCELL(3,ii) = alefvm_buffer%FCELL(3,ii) + f0(3,i) + alefvm_buffer%FEXT_CELL(3,ii)
301 enddo!next I
302
303
304
305 !-------------------------------------------------------------!
306 ! DBUG OUTPUT !
307 !-------------------------------------------------------------!
308! if(IALEFVM_OUTP_FINT /= 0)then
309!
310! if(itask==0)then
311! print *, " |----alefvm_sfint3.F-----|"
312! print *, " | THREAD INFORMATION |"
313! print *, " |------------------------|"
314! print *, " NCYCLE =", NCYCLE
315!
316! do i=1,nb
317! ii = nft + i
318! IF(INT22/=0)THEN
319! IB = NINT(IAD22(I))
320! IF(IB>0) CYCLE
321! ENDIF
322! print *, " brique=", ixs(11,nft+i)
323! write(*,FMT='(A34,6A26)') " ",
324! . "#--- internal force -----#"
325! write (*,FMT='(A,8E26.14)') " force-X =", FCELL(1,II)
326! write (*,FMT='(A,8E26.14)') " force-Y =", FCELL(2,II)
327! write (*,FMT='(A,8E26.14)') " force-Z =", FCELL(3,II)
328! print *, " "
329! enddo
330!
331! endif!if(itask)
332!
333! endif!if(IALEFVM_OUTP_FINT /= 0)
334 !-----------------------------------------!
335
336 RETURN
337 END
subroutine alefvm_sfint3(ixs, nv46, ale_connect, ialefvm_flg, ipm, iparg, ng, x, iad22, nel)
#define min(a, b)
Definition macros.h:20
type(alefvm_buffer_), target alefvm_buffer
Definition alefvm_mod.F:120
type(alefvm_param_), target alefvm_param
Definition alefvm_mod.F:121