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