OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
s10for_edgec.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!|| s10for_edgec ../engine/source/elements/solid/solide10/s10for_edgec.F
25!||--- called by ------------------------------------------------------
26!|| s10for_distor ../engine/source/elements/solid/solide10/s10for_distor.F
27!||====================================================================
28 SUBROUTINE s10for_edgec(
29 . STI, LL , STI_C,
30 . XX , YY , ZZ ,
31 . XX0, YY0, ZZ0,
32 . VX , VY , VZ ,
33 . FOR_T1, FOR_T2, FOR_T3,
34 . FOR_T4, FOR_T5, FOR_T6,
35 . FOR_T7, FOR_T8, FOR_T9,
36 . FOR_T10, TOL_E, IFCE ,
37 . IFCTL , NEL ,E_DISTOR,
38 . DT1 )
39C-----------------------------------------------
40C M o d u l e s
41C-----------------------------------------------
42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46C-----------------------------------------------
47C G l o b a l P a r a m e t e r s
48C-----------------------------------------------
49#include "mvsiz_p.inc"
50C-----------------------------------------------
51C C o m m o n B l o c k s
52C-----------------------------------------------
53C D u m m y A r g u m e n t s
54C-----------------------------------------------
55 INTEGER, INTENT(IN) :: NEL
56 INTEGER, INTENT (OUT) :: IFCTL
57 INTEGER, DIMENSION(MVSIZ), INTENT(OUT) :: IFCE
58 my_real, INTENT(IN) :: TOL_E,DT1
59 my_real, DIMENSION(MVSIZ), INTENT(IN) :: LL
60 DOUBLE PRECISION, DIMENSION(MVSIZ,10), INTENT(IN) ::
61 . XX, YY, ZZ,
62 . XX0,YY0,ZZ0
63 my_real, DIMENSION(MVSIZ), INTENT(INOUT) ::STI,STI_C
64 my_real, DIMENSION(MVSIZ,3), INTENT (INOUT) :: FOR_T10,
65 . for_t1, for_t2, for_t3,
66 . for_t4, for_t5, for_t6,
67 . for_t7, for_t8, for_t9
68 my_real, DIMENSION(MVSIZ,10), INTENT(IN) ::
69 . vx, vy, vz
70 my_real, DIMENSION(NEL), INTENT(INOUT) :: e_distor
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 my_real
75 . stif(mvsiz),dmin(mvsiz),dmin2(mvsiz),fn(mvsiz,6),
76 . xm(mvsiz,6),ym(mvsiz,6),zm(mvsiz,6),ks(mvsiz),
77 . njx(mvsiz,6),njy(mvsiz,6),njz(mvsiz,6),fact,
78 . vxm(mvsiz,6),vym(mvsiz,6),vzm(mvsiz,6),
79 . dx,dy,dz,d_max,lam,s2,norm,lam0,f_q,fac,leng,
80 . nx,ny,nz,fx,fy,fz,f_t,ds2(mvsiz,6),tol2,kts,fnj,
81 . ddx,ddy,ddz,ddn
82 INTEGER I,J,JJ,IFC2(MVSIZ,6),IFF
83C-----------------------------------------------
84C P r e - C o n d i t i o n s
85C-----------------------------------------------
86C-----------------------------------------------
87C S o u r c e C o d e
88C-----------------------------------------------
89C---- 5(1,2),6(2,3),7(1,3),8(1,4),9(2,4),10(3,4)
90 ifctl = 0
91 DO i=1,nel
92 IF (sti_c(i)==zero) cycle
93 xm(i,1) = half*(xx(i,1)+xx(i,2))
94 ym(i,1) = half*(yy(i,1)+yy(i,2))
95 zm(i,1) = half*(zz(i,1)+zz(i,2))
96 xm(i,2) = half*(xx(i,3)+xx(i,2))
97 ym(i,2) = half*(yy(i,3)+yy(i,2))
98 zm(i,2) = half*(zz(i,3)+zz(i,2))
99 xm(i,3) = half*(xx(i,3)+xx(i,1))
100 ym(i,3) = half*(yy(i,3)+yy(i,1))
101 zm(i,3) = half*(zz(i,3)+zz(i,1))
102 xm(i,4) = half*(xx(i,4)+xx(i,1))
103 ym(i,4) = half*(yy(i,4)+yy(i,1))
104 zm(i,4) = half*(zz(i,4)+zz(i,1))
105 xm(i,5) = half*(xx(i,4)+xx(i,2))
106 ym(i,5) = half*(yy(i,4)+yy(i,2))
107 zm(i,5) = half*(zz(i,4)+zz(i,2))
108 xm(i,6) = half*(xx(i,4)+xx(i,3))
109 ym(i,6) = half*(yy(i,4)+yy(i,3))
110 zm(i,6) = half*(zz(i,4)+zz(i,3))
111 stif(i) = sti_c(i)
112 dmin(i) = tol_e*ll(i)
113 dmin2(i) = dmin(i)*dmin(i)
114 ifce(i)=0
115 ENDDO
116C---- first sorting
117 iff = 0
118 ifc2 = 0
119 DO j=1,6
120 jj = j + 4
121 DO i=1,nel
122 IF (sti_c(i)==zero) cycle
123 dx = xx(i,jj) - xm(i,j)
124 dy = yy(i,jj) - ym(i,j)
125 dz = zz(i,jj) - zm(i,j)
126 d_max = two*max(abs(dx),abs(dy),abs(dz))
127 IF (d_max>dmin(i)) THEN
128 ifc2(i,j) = 1
129 iff = 1
130 END IF
131 ENDDO
132 ENDDO
133 IF (iff >0) THEN
134 DO i=1,nel
135 IF (sti_c(i)==zero) cycle
136 vxm(i,1) = half*(vx(i,1)+vx(i,2))
137 vym(i,1) = half*(vy(i,1)+vy(i,2))
138 vzm(i,1) = half*(vz(i,1)+vz(i,2))
139 vxm(i,2) = half*(vx(i,3)+vx(i,2))
140 vym(i,2) = half*(vy(i,3)+vy(i,2))
141 vzm(i,2) = half*(vz(i,3)+vz(i,2))
142 vxm(i,3) = half*(vx(i,3)+vx(i,1))
143 vym(i,3) = half*(vy(i,3)+vy(i,1))
144 vzm(i,3) = half*(vz(i,3)+vz(i,1))
145 vxm(i,4) = half*(vx(i,4)+vx(i,1))
146 vym(i,4) = half*(vy(i,4)+vy(i,1))
147 vzm(i,4) = half*(vz(i,4)+vz(i,1))
148 vxm(i,5) = half*(vx(i,4)+vx(i,2))
149 vym(i,5) = half*(vy(i,4)+vy(i,2))
150 vzm(i,5) = half*(vz(i,4)+vz(i,2))
151 vxm(i,6) = half*(vx(i,4)+vx(i,3))
152 vym(i,6) = half*(vy(i,4)+vy(i,3))
153 vzm(i,6) = half*(vz(i,4)+vz(i,3))
154 ENDDO
155 tol2 = tol_e*tol_e
156C
157 DO i=1,nel
158 IF (ifc2(i,1)==0) cycle
159 dx = xx0(i,1)-xx0(i,2)
160 dy = yy0(i,1)-yy0(i,2)
161 dz = zz0(i,1)-zz0(i,2)
162 ds2(i,1) = dx*dx+dy*dy+dz*dz
163 ENDDO
164 DO i=1,nel
165 IF (ifc2(i,2)==0) cycle
166 dx = xx0(i,3)-xx0(i,2)
167 dy = yy0(i,3)-yy0(i,2)
168 dz = zz0(i,3)-zz0(i,2)
169 ds2(i,2) = dx*dx+dy*dy+dz*dz
170 ENDDO
171 DO i=1,nel
172 IF (ifc2(i,3)==0) cycle
173 dx = xx0(i,1)-xx0(i,3)
174 dy = yy0(i,1)-yy0(i,3)
175 dz = zz0(i,1)-zz0(i,3)
176 ds2(i,3) = dx*dx+dy*dy+dz*dz
177 ENDDO
178 DO i=1,nel
179 IF (ifc2(i,4)==0) cycle
180 dx = xx0(i,1)-xx0(i,4)
181 dy = yy0(i,1)-yy0(i,4)
182 dz = zz0(i,1)-zz0(i,4)
183 ds2(i,4) = dx*dx+dy*dy+dz*dz
184 ENDDO
185 DO i=1,nel
186 IF (ifc2(i,5)==0) cycle
187 dx = xx0(i,4)-xx0(i,2)
188 dy = yy0(i,4)-yy0(i,2)
189 dz = zz0(i,4)-zz0(i,2)
190 ds2(i,5) = dx*dx+dy*dy+dz*dz
191 ENDDO
192 DO i=1,nel
193 IF (ifc2(i,6)==0) cycle
194 dx = xx0(i,3)-xx0(i,4)
195 dy = yy0(i,3)-yy0(i,4)
196 dz = zz0(i,3)-zz0(i,4)
197 ds2(i,6) = dx*dx+dy*dy+dz*dz
198 END DO
199C----
200 fn = zero
201 f_t = three
202 f_q = 0.67
203 lam0 = tol_e
204 DO j=1,6
205 jj = j + 4
206 DO i=1,nel
207 IF (ifc2(i,j)==0) cycle
208 dx = xx(i,jj) - xm(i,j)
209 dy = yy(i,jj) - ym(i,j)
210 dz = zz(i,jj) - zm(i,j)
211 s2 = dx*dx+dy*dy+dz*dz
212 IF (s2>tol2*ds2(i,j)) THEN
213 ifctl = 1
214 lam = sqrt(s2/ds2(i,j))
215 norm=one/sqrt(s2)
216 njx(i,j) = dx*norm
217 njy(i,j) = dy*norm
218 njz(i,j) = dz*norm
219 leng = sqrt(ds2(i,j))
220 fac = f_q*(lam+one)*(lam+one)
221 kts = f_t*(fac+one)*sti_c(i)
222 fn(i,j) = kts*(lam-lam0)*leng
223 fact = f_t*(three*fac+one)
224 stif(i) = max(stif(i),fact*sti_c(i)) ! will be added to STI not max()
225!
226 ddx = vx(i,jj) - vxm(i,j)
227 ddy = vy(i,jj) - vym(i,j)
228 ddz = vz(i,jj) - vzm(i,j)
229 ddn = dt1*(ddx*njx(i,j)+ddy*njy(i,j)+ddz*njz(i,j))
230 e_distor(i) = e_distor(i) - fn(i,j)*ddn
231 END IF
232 ENDDO
233 ENDDO
234C ---- 5 (1,2)
235 j = 1
236 DO i=1,nel
237 IF (fn(i,j)==zero) cycle
238 fnj = half*fn(i,j)
239 nx = njx(i,j)
240 ny = njy(i,j)
241 nz = njz(i,j)
242 fx = fnj*nx
243 fy = fnj*ny
244 fz = fnj*nz
245 for_t1(i,1) = for_t1(i,1) + fx
246 for_t1(i,2) = for_t1(i,2) + fy
247 for_t1(i,3) = for_t1(i,3) + fz
248 for_t2(i,1) = for_t2(i,1) + fx
249 for_t2(i,2) = for_t2(i,2) + fy
250 for_t2(i,3) = for_t2(i,3) + fz
251!
252 for_t5(i,1) = for_t5(i,1) - two*fx
253 for_t5(i,2) = for_t5(i,2) - two*fy
254 for_t5(i,3) = for_t5(i,3) - two*fz
255 ifce(i)=j
256 ENDDO
257C---- 6(2,3)
258 j = 2
259 DO i=1,nel
260 IF (fn(i,j)==zero) cycle
261 fnj = half*fn(i,j)
262 nx = njx(i,j)
263 ny = njy(i,j)
264 nz = njz(i,j)
265 fx = fnj*nx
266 fy = fnj*ny
267 fz = fnj*nz
268 for_t3(i,1) = for_t3(i,1) + fx
269 for_t3(i,2) = for_t3(i,2) + fy
270 for_t3(i,3) = for_t3(i,3) + fz
271 for_t2(i,1) = for_t2(i,1) + fx
272 for_t2(i,2) = for_t2(i,2) + fy
273 for_t2(i,3) = for_t2(i,3) + fz
274!
275 for_t6(i,1) = for_t6(i,1) - two*fx
276 for_t6(i,2) = for_t6(i,2) - two*fy
277 for_t6(i,3) = for_t6(i,3) - two*fz
278 ifce(i)=j
279 ENDDO
280C---- 7(1,3)
281 j = 3
282 DO i=1,nel
283 IF (fn(i,j)==zero) cycle
284 fnj = half*fn(i,j)
285 nx = njx(i,j)
286 ny = njy(i,j)
287 nz = njz(i,j)
288 fx = fnj*nx
289 fy = fnj*ny
290 fz = fnj*nz
291 for_t3(i,1) = for_t3(i,1) + fx
292 for_t3(i,2) = for_t3(i,2) + fy
293 for_t3(i,3) = for_t3(i,3) + fz
294 for_t1(i,1) = for_t1(i,1) + fx
295 for_t1(i,2) = for_t1(i,2) + fy
296 for_t1(i,3) = for_t1(i,3) + fz
297 for_t7(i,1) = for_t7(i,1) - two*fx
298 for_t7(i,2) = for_t7(i,2) - two*fy
299 for_t7(i,3) = for_t7(i,3) - two*fz
300 ifce(i)=j
301 ENDDO
302C---- 8(1,4)
303 j = 4
304 DO i=1,nel
305 IF (fn(i,j)==zero) cycle
306 fnj = half*fn(i,j)
307 nx = njx(i,j)
308 ny = njy(i,j)
309 nz = njz(i,j)
310 fx = fnj*nx
311 fy = fnj*ny
312 fz = fnj*nz
313 for_t1(i,1) = for_t1(i,1) + fx
314 for_t1(i,2) = for_t1(i,2) + fy
315 for_t1(i,3) = for_t1(i,3) + fz
316 for_t4(i,1) = for_t4(i,1) + fx
317 for_t4(i,2) = for_t4(i,2) + fy
318 for_t4(i,3) = for_t4(i,3) + fz
319 for_t8(i,1) = for_t8(i,1) - two*fx
320 for_t8(i,2) = for_t8(i,2) - two*fy
321 for_t8(i,3) = for_t8(i,3) - two*fz
322 ifce(i)=j
323 ENDDO
324C---- 9(2,4)
325 j = 5
326 DO i=1,nel
327 IF (fn(i,j)==zero) cycle
328 fnj = half*fn(i,j)
329 nx = njx(i,j)
330 ny = njy(i,j)
331 nz = njz(i,j)
332 fx = fnj*nx
333 fy = fnj*ny
334 fz = fnj*nz
335 for_t4(i,1) = for_t4(i,1) + fx
336 for_t4(i,2) = for_t4(i,2) + fy
337 for_t4(i,3) = for_t4(i,3) + fz
338 for_t2(i,1) = for_t2(i,1) + fx
339 for_t2(i,2) = for_t2(i,2) + fy
340 for_t2(i,3) = for_t2(i,3) + fz
341 for_t9(i,1) = for_t9(i,1) - two*fx
342 for_t9(i,2) = for_t9(i,2) - two*fy
343 for_t9(i,3) = for_t9(i,3) - two*fz
344 ifce(i)=j
345 ENDDO
346C----10(3,4)
347 j = 6
348 DO i=1,nel
349 IF (fn(i,j)==zero) cycle
350 fnj = half*fn(i,j)
351 nx = njx(i,j)
352 ny = njy(i,j)
353 nz = njz(i,j)
354 fx = fnj*nx
355 fy = fnj*ny
356 fz = fnj*nz
357 for_t4(i,1) = for_t4(i,1) + fx
358 for_t4(i,2) = for_t4(i,2) + fy
359 for_t4(i,3) = for_t4(i,3) + fz
360 for_t3(i,1) = for_t3(i,1) + fx
361 for_t3(i,2) = for_t3(i,2) + fy
362 for_t3(i,3) = for_t3(i,3) + fz
363 for_t10(i,1) = for_t10(i,1) - two*fx
364 for_t10(i,2) = for_t10(i,2) - two*fy
365 for_t10(i,3) = for_t10(i,3) - two*fz
366 ifce(i)=j
367 ENDDO
368 DO i=1,nel
369 IF (sti_c(i)==zero) cycle
370 sti(i) = max(sti(i),stif(i))
371 ENDDO
372 END IF !(IFF >0) THEN
373
374
375 RETURN
376 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
max
#define max(a, b)
Definition macros.h:21
s10for_edgec
subroutine s10for_edgec(sti, ll, sti_c, xx, yy, zz, xx0, yy0, zz0, vx, vy, vz, for_t1, for_t2, for_t3, for_t4, for_t5, for_t6, for_t7, for_t8, for_t9, for_t10, tol_e, ifce, ifctl, nel, e_distor, dt1)
Definition s10for_edgec.F:39