OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
press_seg3.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!|| press_seg3 ../engine/source/loads/general/load_pcyl/press_seg3.F
25!||--- called by ------------------------------------------------------
26!|| h3d_pre_skin_scalar ../engine/source/output/h3d/h3d_results/h3d_skin_scalar.F
27!|| h3d_skin_vector ../engine/source/output/h3d/h3d_results/h3d_skin_vector.F
28!|| pressure_cyl ../engine/source/loads/general/load_pcyl/pressure_cyl.F
29!||--- calls -----------------------------------------------------
30!|| table_interp_dydx ../engine/source/tools/curve/table_tools.F
31!||--- uses -----------------------------------------------------
32!|| pload_cyl_mod ../common_source/modules/loads/pload_cyl_mod.F
33!|| table_mod ../engine/share/modules/table_mod.F
34!||====================================================================
35 SUBROUTINE press_seg3(A ,B ,C ,N1 ,DIR ,
36 . IFUNC ,TABLE ,XFACR ,XFACT ,PRESS )
37C-----------------------------------------------
38C M o d u l e s
39C-----------------------------------------------
40 USE table_mod
41 USE pload_cyl_mod
42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46C-----------------------------------------------
47C C o m m o n B l o c k s
48C-----------------------------------------------
49#include "com04_c.inc"
50#include "com08_c.inc"
51C-----------------------------------------------
52C D u m m y A r g u m e n t s
53C-----------------------------------------------
54 INTEGER ,INTENT(IN) :: IFUNC
55 my_real ,INTENT(IN) :: xfacr !< radius scale factor
56 my_real ,INTENT(IN) :: xfact !< time scale factor
57 my_real ,DIMENSION(3) ,INTENT(IN) :: a,b,c,n1,dir
58 TYPE (TTABLE), DIMENSION(NTABLE) ,INTENT(IN) :: TABLE
59 my_real ,INTENT(OUT) :: press
60C-----------------------------------------------
61C L o c a l V a r i a b l e s
62C-----------------------------------------------
63 INTEGER :: NPT
64 my_real :: a11,a12,a21,a22,b1,b2,det,alpha,beta,gamma,f1,f2,r,s,
65 . r1,r2,r3,rmax,p1,p2,p3,dydx,nx,ny,nz,norm,cosp,sinp
66 my_real, DIMENSION(3) :: ab,ac,an,p,ap,bp,cp
67 my_real, DIMENSION(2) :: xx
68
69
70c-----------------------------------------------
71c 1) find intersection point P between plane surface A-B-C and line N1-N2
72c 2) calculate radial coordinates of each segment node in cylindrical
73c coord system defined by direction axis
74c (distance between point P and each segment node projected
75c on plane perpendicular to axis)
76c=======================================================================
77 ab(1) = b(1) - a(1)
78 ab(2) = b(2) - a(2)
79 ab(3) = b(3) - a(3)
80 ac(1) = c(1) - a(1)
81 ac(2) = c(2) - a(2)
82 ac(3) = c(3) - a(3)
83 an(1) = n1(1) - a(1)
84 an(2) = n1(2) - a(2)
85 an(3) = n1(3) - a(3)
86 alpha = dir(1)
87 beta = dir(2)
88 gamma = dir(3)
89 npt = SIZE(table(ifunc)%X(1)%VALUES)
90 rmax = table(ifunc)%X(1)%VALUES(npt) * xfacr
91 p1 = zero
92 p2 = zero
93 p3 = zero
94c
95 IF (alpha /= zero) THEN
96 f1 = beta / alpha
97 f2 = gamma / alpha
98 b1 = an(2) - f1 * an(1)
99 b2 = an(3) - f2 * an(1)
100 a11 = ab(2) - f1 * ab(1)
101 a12 = ac(2) - f1 * ac(1)
102 a21 = ab(3) - f2 * ab(1)
103 a22 = ac(3) - f2 * ac(1)
104 ELSE IF (beta /= zero) THEN
105 f1 = alpha / beta
106 f2 = gamma / beta
107 b1 = an(1) - f1 * an(2)
108 b2 = an(3) - f2 * an(2)
109 a11 = ab(1) - f1 * ab(2)
110 a12 = ac(1) - f1 * ac(2)
111 a21 = ab(3) - f2 * ab(2)
112 a22 = ac(3) - f2 * ac(2)
113 ELSE IF (gamma /= zero) THEN
114 f1 = alpha / gamma
115 f2 = beta / gamma
116 b1 = an(1) - f1 * an(3)
117 b2 = an(2) - f2 * an(3)
118 a11 = ab(1) - f1 * ab(3)
119 a12 = ac(1) - f1 * ac(3)
120 a21 = ab(2) - f2 * ab(3)
121 a22 = ac(2) - f2 * ac(3)
122 ELSE
123 f1 = 0
124 f2 = 0
125 b1 = 0
126 b2 = 0
127 a11 = 0
128 a12 = 0
129 a21 = 0
130 a22 = 0
131 END IF
132c
133 det = a11*a22 - a12*a21
134 IF(det .NE. zero) THEN
135 r = (a22 * b1 - a12 * b2) / det
136 s = (a11 * b2 - a21 * b1) / det
137 ELSE
138 r = 0
139 s = 0
140 END IF
141 ! P = projection of axe origin N1 on ABC plane following direction DIR
142 p(1)= a(1) + r*ab(1) + s*ac(1)
143 p(2)= a(2) + r*ab(2) + s*ac(2)
144 p(3)= a(3) + r*ab(3) + s*ac(3)
145 ap(1) = p(1) - a(1)
146 ap(2) = p(2) - a(2)
147 ap(3) = p(3) - a(3)
148 bp(1) = p(1) - b(1)
149 bp(2) = p(2) - b(2)
150 bp(3) = p(3) - b(3)
151 cp(1) = p(1) - c(1)
152 cp(2) = p(2) - c(2)
153 cp(3) = p(3) - c(3)
154 r1 = sqrt(ap(1)**2 + ap(2)**2 + ap(3)**2)
155 r2 = sqrt(bp(1)**2 + bp(2)**2 + bp(3)**2)
156 r3 = sqrt(cp(1)**2 + cp(2)**2 + cp(3)**2)
157 cosp = (dir(1) * ap(1) + dir(2) * ap(2) + dir(3) * ap(3)) / max(r1,em20)
158 sinp = max(sqrt(one - cosp**2), em20)
159 p1 = zero
160 p2 = zero
161 p3 = zero
162 xx(2) = tt / xfact ! Time
163c
164c pressure in 1st point of 3N sub-segment
165 cosp = (dir(1) * ap(1) + dir(2) * ap(2) + dir(3) * ap(3)) / max(r1,em20)
166 sinp = max(sqrt(one - cosp**2), em20)
167 r1 = r1 * sinp
168 IF (r1 <= rmax) THEN
169 xx(1) = r1/xfacr
170 CALL table_interp_dydx(table(ifunc),xx,2,p1,dydx)
171 END IF
172c pressure in 2nd point of 3N sub-segment
173 cosp = (dir(1) * bp(1) + dir(2) * bp(2) + dir(3) * bp(3)) / max(r2,em20)
174 sinp = max(sqrt(one - cosp**2), em20)
175 r2 = r2 * sinp
176 IF (r2 <= rmax) THEN
177 xx(1) = r2/xfacr
178 CALL table_interp_dydx(table(ifunc),xx,2,p2,dydx)
179 END IF
180c pressure in 3rd point of 3N sub-segment
181 cosp = (dir(1) * cp(1) + dir(2) * cp(2) + dir(3) * cp(3)) / max(r3,em20)
182 sinp = max(sqrt(one - cosp**2), em20)
183 r3 = r3 * sinp
184 IF (r3 <= rmax) THEN
185 xx(1) = r3/xfacr
186 CALL table_interp_dydx(table(ifunc),xx,2,p3,dydx)
187 END IF
188 press = (p1 + p2 + p3) * third
189 IF (press > zero) THEN
190c calculate cos(angle) between cylinder axis and segment normal to scale pressure
191 nx = ab(2) * ac(3) - ab(3) * ac(2)
192 ny = ab(3) * ac(1) - ab(1) * ac(3)
193 nz = ab(1) * ac(2) - ab(2) * ac(1)
194 norm = sqrt(nx**2 + ny**2 + nz**2)
195 cosp = abs(nx * alpha + ny * beta + nz * gamma) / norm
196 press = press * cosp
197 END IF
198c-----------
199 RETURN
200 END
201
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
table_interp_dydx
subroutine table_interp_dydx(table, xx, xxdim, yy, dydx)
Definition table_tools.F:419
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21
pload_cyl_mod
Definition pload_cyl_mod.F:35
table_mod
Definition table_mod.F:113
press_seg3
subroutine press_seg3(a, b, c, n1, dir, ifunc, table, xfacr, xfact, press)
Definition press_seg3.F:37