OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
gradient_reconstruction.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!|| gradient_reconstruction ../engine/source/ale/alemuscl/gradient_reconstruction.F
25!||--- called by ------------------------------------------------------
26!|| ale51_gradient_reconstruction ../engine/source/ale/alemuscl/ale51_gradient_reconstruction.F
27!||--- calls -----------------------------------------------------
28!|| cg ../engine/source/ale/alemuscl/conjugate_gradient.F
29!||--- uses -----------------------------------------------------
30!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
31!|| alemuscl_mod ../common_source/modules/ale/alemuscl_mod.F
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!|| segvar_mod ../engine/share/modules/segvar_mod.F
34!||====================================================================
35 SUBROUTINE gradient_reconstruction(IXS, X, ALE_CONNECT, ITRIMAT,SEGVAR)
36C-----------------------------------------------
37C D e s c r i p t i o n
38C This subroutine computes a gradient of the scalar field value in each
39C element:
40C mean square approximation of the gradient in the face related
41C neighborhood of the element
42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE alemuscl_mod
46 USE segvar_mod
47 USE ale_connectivity_mod
48 use element_mod , only :nixs
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "vect01_c.inc"
57#include "com04_c.inc"
58C-----------------------------------------------
59C D u m m y A r g u m e n t s
60C-----------------------------------------------
61 INTEGER, INTENT(IN) :: IXS(NIXS,NUMELS)
62 my_real, INTENT(IN) :: x(3,numnod)
63 INTEGER, INTENT(IN) :: ITRIMAT
64 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
65C-----------------------------------------------
66C L o c a l V a r i a b l e s
67C-----------------------------------------------
68 INTEGER :: I, II, KK, IAD2, LGTH
69 my_real :: xk, yk, zk, xl, yl, zl, xf, yf, zf
70 my_real :: valk, vall
71 my_real :: mat(3, 3), rhs(3), sol(3)
72 INTEGER :: VOIS_ID
73 INTEGER :: FACE_TO_NODE_LOCAL_ID(6, 4)
74 my_real :: norm(3), a(3), b(3), c(3), surf, surf1, surf2
75 TYPE(t_segvar) :: SEGVAR
76C-----------------------------------------------
77C S o u r c e L i n e s
78C-----------------------------------------------
79!!! Once for all, associate node local id to a face number
80!!! Face 1
81 face_to_node_local_id(1, 1) = 1 ; face_to_node_local_id(1, 2) = 4
82 face_to_node_local_id(1, 3) = 3 ; face_to_node_local_id(1, 4) = 2
83!!! Face 2
84 face_to_node_local_id(2, 1) = 3 ; face_to_node_local_id(2, 2) = 4
85 face_to_node_local_id(2, 3) = 8 ; face_to_node_local_id(2, 4) = 7
86!!! Face 3
87 face_to_node_local_id(3, 1) = 5 ; face_to_node_local_id(3, 2) = 6
88 face_to_node_local_id(3, 3) = 7 ; face_to_node_local_id(3, 4) = 8
89!!! Face 4
90 face_to_node_local_id(4, 1) = 1 ; face_to_node_local_id(4, 2) = 2
91 face_to_node_local_id(4, 3) = 6 ; face_to_node_local_id(4, 4) = 5
92!!! Face 5
93 face_to_node_local_id(5, 1) = 2 ; face_to_node_local_id(5, 2) = 3
94 face_to_node_local_id(5, 3) = 7 ; face_to_node_local_id(5, 4) = 6
95!!! Face 6
96 face_to_node_local_id(6, 1) = 1 ; face_to_node_local_id(6, 2) = 5
97 face_to_node_local_id(6, 3) = 8 ; face_to_node_local_id(6, 4) = 4
98
99 DO i = lft, llt
100 ii = i + nft
101 !!! Reset mat, rhs
102 mat(1:3, 1:3) = zero ; rhs(1:3) = zero
103 !!! Value of the target function in the element
104 valk = alemuscl_buffer%VOLUME_FRACTION(ii,itrimat)
105 xk = alemuscl_buffer%ELCENTER(ii,1) ;
106 yk = alemuscl_buffer%ELCENTER(ii,2) ;
107 zk = alemuscl_buffer%ELCENTER(ii,3)
108 !!! IXS(2:9, II) : Node global ID
109 iad2 = ale_connect%ee_connect%iad_connect(ii)
110 lgth = ale_connect%ee_connect%iad_connect(ii+1)-iad2
111 DO kk = 1, lgth
112 vois_id = ale_connect%ee_connect%connected(iad2 + kk - 1)
113 !IF (VOIS_ID /= 0 .AND. VOIS_ID <= NUMELS) THEN
114 IF (vois_id > 0) THEN
115 !!! Value of the target function in the current neighbor
116 vall = alemuscl_buffer%VOLUME_FRACTION(vois_id,itrimat)
117 xl = alemuscl_buffer%ELCENTER(vois_id,1) ;
118 yl = alemuscl_buffer%ELCENTER(vois_id,2) ;
119 zl = alemuscl_buffer%ELCENTER(vois_id,3) ;
120 ELSE
121 IF(vois_id == 0) THEN
122 vall = valk
123 ELSE
124 !vois_id<0 : means EBCS), -vois_id is seg_id
125 vall = segvar%PHASE_ALPHA(itrimat,-vois_id)
126 ENDIF
127 xf = zero
128 yf = zero
129 zf = zero
130
131 a(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 1) + 1, ii))
132 b(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 2) + 1, ii))
133 c(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 3) + 1, ii))
134
135 norm(1) = (b(2) - a(2)) * (c(3) - a(3)) - (b(3) - a(3)) * (c(2) - a(2))
136 norm(2) = (b(3) - a(3)) * (c(1) - a(1)) - (b(1) - a(1)) * (c(3) - a(3))
137 norm(3) = (b(1) - a(1)) * (c(2) - a(2)) - (b(2) - a(2)) * (c(1) - a(1))
138
139 surf1 = half * abs(sqrt(norm(1) * norm(1) + norm(2) * norm(2) + norm(3) * norm(3)))
140 xf = surf1 * third * (a(1) + b(1) + c(1))
141 yf = surf1 * third * (a(2) + b(2) + c(2))
142 zf = surf1 * third * (a(3) + b(3) + c(3))
143
144 a(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 1) + 1, ii))
145 b(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 3) + 1, ii))
146 c(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 4) + 1, ii))
147
148 norm(1) = (b(2) - a(2)) * (c(3) - a(3)) - (b(3) - a(3)) * (c(2) - a(2))
149 norm(2) = (b(3) - a(3)) * (c(1) - a(1)) - (b(1) - a(1)) * (c(3) - a(3))
150 norm(3) = (b(1) - a(1)) * (c(2) - a(2)) - (b(2) - a(2)) * (c(1) - a(1))
151
152 surf2 = half * abs(sqrt(norm(1) * norm(1) + norm(2) * norm(2) + norm(3) * norm(3)))
153 xf = xf + surf2 * third * (a(1) + b(1) + c(1))
154 yf = yf + surf2 * third * (a(2) + b(2) + c(2))
155 zf = zf + surf2 * third * (a(3) + b(3) + c(3))
156
157 surf = surf1 + surf2
158 xf = xf / surf
159 yf = yf / surf
160 zf = zf / surf
161 !!! Build face centroid
162c$$$ DO JJ = 1, 4
163c$$$ NODE_ID = IXS(FACE_TO_NODE_LOCAL_ID(KK, JJ) + 1, II)
164c$$$ XF = XF + FOURTH * X(1, NODE_ID)
165c$$$ YF = YF + FOURTH * X(2, NODE_ID)
166c$$$ ZF = ZF + FOURTH * X(3, NODE_ID)
167c$$$ ENDDO
168 xl = two * xf - alemuscl_buffer%ELCENTER(ii,1)
169 yl = two * yf - alemuscl_buffer%ELCENTER(ii,2)
170 zl = two * zf - alemuscl_buffer%ELCENTER(ii,3)
171 ENDIF
172
173 !!! Incrementing mat and rhs
174 rhs(1) = rhs(1) + (valk - vall) * (xl - xk)
175 rhs(2) = rhs(2) + (valk - vall) * (yl - yk)
176 rhs(3) = rhs(3) + (valk - vall) * (zl - zk)
177 mat(1, 1) = mat(1, 1) + (xl - xk) * (xl - xk)
178 mat(1, 2) = mat(1, 2) + (xl - xk) * (yl - yk)
179 mat(1, 3) = mat(1, 3) + (xl - xk) * (zl - zk)
180 mat(2, 1) = mat(2, 1) + (yl - yk) * (xl - xk)
181 mat(2, 2) = mat(2, 2) + (yl - yk) * (yl - yk)
182 mat(2, 3) = mat(2, 3) + (yl - yk) * (zl - zk)
183 mat(3, 1) = mat(3, 1) + (zl - zk) * (xl - xk)
184 mat(3, 2) = mat(3, 2) + (zl - zk) * (yl - yk)
185 mat(3, 3) = mat(3, 3) + (zl - zk) * (zl - zk)
186 ENDDO
187
188 CALL cg(3, mat, rhs, sol, 3, em10)
189 !CALL DIRECT_SOLVE(mat, rhs ,sol)
190 !!! Solution goes to the gradient
191 alemuscl_buffer%GRAD(ii,1,itrimat) = -sol(1)
192 alemuscl_buffer%GRAD(ii,2,itrimat) = -sol(2)
193 alemuscl_buffer%GRAD(ii,3,itrimat) = -sol(3)
194 ENDDO ! I = LFT, LLT
195C-----------------------------------------------
196 END SUBROUTINE gradient_reconstruction
cg
subroutine cg(dim, mat, rhs, sol, max_iter, tol)
Definition conjugate_gradient.F:29
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
gradient_reconstruction
subroutine gradient_reconstruction(ixs, x, ale_connect, itrimat, segvar)
Definition gradient_reconstruction.F:36
ale_connectivity_mod
Definition ale_connectivity_mod.F:225
alemuscl_mod
Definition alemuscl_mod.F:55
alemuscl_mod::alemuscl_buffer
type(alemuscl_buffer_) alemuscl_buffer
Definition alemuscl_mod.F:75
segvar_mod
Definition segvar_mod.F:63
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:256
segvar_mod::t_segvar
Definition segvar_mod.F:71