OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ratio_fill.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!|| ratio_fill ../starter/source/initial_conditions/inivol/ratio_fill.F
25!||--- called by ------------------------------------------------------
26!|| inisoldist ../starter/source/initial_conditions/inivol/inisoldist.F
27!||--- calls -----------------------------------------------------
28!|| in_out_side ../starter/source/initial_conditions/inivol/in_out_side.F
29!||--- uses -----------------------------------------------------
30!|| inivol_def_mod ../starter/share/modules1/inivol_mod.F
31!||====================================================================
32 SUBROUTINE ratio_fill(
33 . X1 ,X2 ,X3 ,X4 ,X5 ,X6 ,X7 ,X8 ,
34 . Y1 ,Y2 ,Y3 ,Y4 ,Y5 ,Y6 ,Y7 ,Y8 ,
35 . Z1 ,Z2 ,Z3 ,Z4 ,Z5 ,Z6 ,Z7 ,Z8 ,
36 . IDP ,X ,IXS ,IPART_ ,IFILL ,NTRACE ,NTRACE0 ,DIS ,
37 . NSOLTOSF ,NNOD2SURF ,INOD2SURF ,KNOD2SURF ,JMID ,IPHASE ,INPHASE ,KVOL ,
38 . SURF_TYPE ,IAD_BUFR ,BUFSF ,NOD_NORMAL ,ISOLNOD ,NBSUBMAT,FILL_RATIO ,ICUMU ,
39 . NSEG ,SURF_ELTYP,SURF_NODES,NBCONTY ,IDC ,NBIP ,IDSURF ,SWIFTSURF,
40 . SEGTOSURF ,IGRSURF ,IVOLSURF ,NSURF_INVOL ,IXQ ,IXTG ,ITYP ,NEL ,
41 . NUMEL_TOT ,NUM_INIVOL,INIVOL ,I_INIVOL)
42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE groupdef_mod
46 USE inivol_def_mod , ONLY : inivol_struct_
47 use element_mod , only :nixs,nixq,nixtg
48C-----------------------------------------------
49C I m p l i c i t T y p e s
50C-----------------------------------------------
51#include "implicit_f.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56C-----------------------------------------------
57C C o m m o n B l o c k s
58C-----------------------------------------------
59#include "com01_c.inc"
60#include "com04_c.inc"
61C-----------------------------------------------
62C D u m m y A r g u m e n t s
63C-----------------------------------------------
64 INTEGER,INTENT(IN) :: I_INIVOL !< inivol identifier
65 INTEGER,INTENT(IN) :: NUM_INIVOL !< number of inivol options
66 TYPE (INIVOL_STRUCT_), DIMENSION(NUM_INIVOL), INTENT(INOUT) :: INIVOL !< inivol data structure
67 INTEGER, INTENT(IN) :: NEL !< number of element in current group NG
68 INTEGER, INTENT(IN) :: NBSUBMAT !< total number of submaterial
69 my_real, INTENT(INOUT) :: KVOL(NBSUBMAT,NEL) !< working array for volume fractions
70 INTEGER NTRACE,NTRACE0,IFILL,JMID,IDP,NSEG,NBCONTY,IDC,NNOD2SURF,
71 . ISOLNOD,ICUMU,SURF_TYPE,IAD_BUFR,IDSURF,IVOLSURF(NSURF),NUMEL_TOT
72 INTEGER IXS(NIXS,NUMELS),IXQ(NIXQ,NUMELQ),IXTG(NIXTG,NUMELTG),IPART_(*),NSOLTOSF(NBCONTY,*),
73 . inod2surf(nnod2surf,*),knod2surf(*),iphase(nbsubmat+1,numel_tot),
74 . inphase(ntrace,nel),surf_eltyp(nseg),ityp,
75 . surf_nodes(nseg,4),nbip(nbsubmat,*),swiftsurf(nsurf),segtosurf(*),nsurf_invol
76 my_real x1(*),x2(*),x3(*),x4(*),x5(*),x6(*),x7(*),x8(*),
77 . y1(*),y2(*),y3(*),y4(*),y5(*),y6(*),y7(*),y8(*),
78 . z1(*),z2(*),z3(*),z4(*),z5(*),z6(*),z7(*),z8(*),
79 . x(3,*),dis(nsurf_invol,*),bufsf(*),nod_normal(3,*),fill_ratio
80 TYPE (SURF_), DIMENSION(NSURF) :: IGRSURF
81C-----------------------------------------------
82C L o c a l V a r i a b l e s
83C-----------------------------------------------
84 INTEGER I,II,J,JJ,K,N,N1,N2,N3,K1,K2,K3,OK,OK1,OK2,OK3,INOD,
85 . ie,nsh,ipl,ip,ixpl(4),getel,nphase,iph,nip,iad0,
86 . ix(8),npoint,ntrace_tot,jmid_old,isurf
87 INTEGER FULL(MVSIZ),JCT(MVSIZ),TRACEP(MVSIZ),TRACEN(MVSIZ)
88 INTEGER BUFFILL1(NBSUBMAT),BUFFILL2(NBSUBMAT,MVSIZ),ELEM_NUMNOD,ISUBMAT_TO_SUBSTRACT
89
90 my_real :: kvol_bak(nbsubmat)
91 my_real
92 . xmin,ymin,zmin,xmax,ymax,zmax,dx,dy,dz,xx(3,8),
93 . xk(ntrace0),yk(ntrace0),zk(ntrace0),xfas(3,4),
94 . x0,y0,z0,dist,dist_old,tmpsum,xn(3),
95 . xk0(ntrace),yk0(ntrace),zk0(ntrace),
96 . l12(3,3),l23(3,3),l31(3,3),ll(3,3),
97 . coef,aaa(3),bbb(3),ccc(3),cg(3)
98 my_real
99 . xs(ntrace,mvsiz),ys(ntrace,mvsiz),zs(ntrace,mvsiz),
100 . disp(ntrace,mvsiz),xp1,yp1,zp1,xp2,yp2,zp2,aa,bb,cc,
101 . xg,yg,zg,skw(9),dgr,tmp(3),x_prime,y_prime,z_prime,
102 . vf_to_substract
103 my_real :: sumvf
104
105 INTEGER :: d1(4),d2(4),d3(4),d4(4)
106 DATA d1/1,2,3,4/,d2/2,3,4,1/,d3/3,4,1,2/,d4/4,1,2,3/
107C-----------------------------------------------
108C S o u r c e L i n e s
109C-----------------------------------------------
110 elem_numnod = -huge(elem_numnod)
111 ! check part to fill with current phase (JMID)
112 k = 0
113 DO i=1,nel
114 IF (ipart_(i) /= idp) cycle
115 k = k + 1
116 ENDDO
117 IF (k == 0) RETURN
118
119 DO i=1,nel
120 full(i) = 0
121 tracep(i) = 0
122 tracen(i) = 0
123 ENDDO
124
125 disp(1:ntrace,1:mvsiz) = zero
126 jmid_old = 0
127
128 ! search for cut solid elems by containers
129 DO i=1,nel
130 IF (ipart_(i) /= idp) cycle
131 k = 0
132 ok1 = 0
133 ok2 = 0
134 ok3 = 0
135 IF(n2d == 0)THEN
136 IF (isolnod == 4) THEN
137 ix(1) =ixs(2,i)
138 ix(2) =ixs(4,i)
139 ix(3) =ixs(7,i)
140 ix(4) =ixs(6,i)
141 elem_numnod = 4
142 ELSEIF (isolnod == 8) THEN
143 ix(1) =ixs(2,i)
144 ix(2) =ixs(3,i)
145 ix(3) =ixs(4,i)
146 ix(4) =ixs(5,i)
147 ix(5) =ixs(6,i)
148 ix(6) =ixs(7,i)
149 ix(7) =ixs(8,i)
150 ix(8) =ixs(9,i)
151 elem_numnod = 8
152 ELSE
153 cycle !not a solid elem
154 ENDIF
155 ELSE
156 IF(ityp == 7)THEN
157 ix(1) =ixtg(2,i)
158 ix(2) =ixtg(3,i)
159 ix(3) =ixtg(4,i)
160 ix(4) =0
161 elem_numnod = 3
162 ELSEIF(ityp == 2)THEN
163 ix(1) =ixq(2,i)
164 ix(2) =ixq(3,i)
165 ix(3) =ixq(4,i)
166 ix(4) =ixq(5,i)
167 elem_numnod = 4
168 ELSE
169 cycle ! not a solid elem
170 ENDIF
171 ENDIF
172 DO j=1,elem_numnod
173 n = ix(j)
174 IF (dis(ivolsurf(idsurf),n) /= zero) THEN
175 k = k + 1
176 IF (dis(ivolsurf(idsurf),n) > zero) THEN
177 ok1 = ok1 + 1
178 ELSEIF (dis(ivolsurf(idsurf),n) < zero) THEN
179 ok2 = ok2 + 1
180 ENDIF
181 ENDIF
182 IF (dis(ivolsurf(idsurf),n) == zero) ok3 = ok3 + 1
183 ENDDO
184
185 IF (k > 0) THEN
186 IF (ok1 == elem_numnod .OR. (ok1+ok3) == elem_numnod) THEN
187 full(i) = 1
188 ELSEIF (ok2 == elem_numnod .OR. (ok2+ok3) == elem_numnod) THEN
189 full(i) = -1
190 ELSEIF (ok1 > 0 .AND. ok2 > 0) THEN
191 full(i) = 2
192 ENDIF
193 ENDIF ! IF (K > 0)
194 ENDDO ! DO I=1,NEL
195
196 ie = 0
197 DO i=1,nel
198 jct(i) = 0
199 IF(full(i) == 2)THEN
200 ie = ie + 1
201 jct(ie) = i
202 END IF
203 END DO
204 getel = ie
205
206 ! Get trace samples coordinates:
207
208 IF (isolnod == 4 .OR. n2d/=0) THEN
209 DO i=1,nel
210 npoint = 0
211 IF (ipart_(i) /= idp) cycle
212 xx(1,1)=x1(i)
213 xx(2,1)=y1(i)
214 xx(3,1)=z1(i)
215 xx(1,2)=x2(i)
216 xx(2,2)=y2(i)
217 xx(3,2)=z2(i)
218 xx(1,3)=x3(i)
219 xx(2,3)=y3(i)
220 xx(3,3)=z3(i)
221 xx(1,4)=x4(i)
222 xx(2,4)=y4(i)
223 xx(3,4)=z4(i)
224 DO k=1,4
225 ! LEVEL - 1 -> 1 SAMPLE POINT
226 npoint = npoint + 1
227 cg(1) = third*(xx(1,d1(k))+xx(1,d2(k))+xx(1,d3(k)))
228 cg(2) = third*(xx(2,d1(k))+xx(2,d2(k))+xx(2,d3(k)))
229 cg(3) = third*(xx(3,d1(k))+xx(3,d2(k))+xx(3,d3(k)))
230 ! COEF = 1.0/4.0
231 coef = fourth
232 aaa(1) = xx(1,d4(k))+coef*(xx(1,d1(k))-xx(1,d4(k)))
233 aaa(2) = xx(2,d4(k))+coef*(xx(2,d1(k))-xx(2,d4(k)))
234 aaa(3) = xx(3,d4(k))+coef*(xx(3,d1(k))-xx(3,d4(k)))
235 bbb(1) = xx(1,d4(k))+coef*(xx(1,d2(k))-xx(1,d4(k)))
236 bbb(2) = xx(2,d4(k))+coef*(xx(2,d2(k))-xx(2,d4(k)))
237 bbb(3) = xx(3,d4(k))+coef*(xx(3,d2(k))-xx(3,d4(k)))
238 ccc(1) = xx(1,d4(k))+coef*(xx(1,d3(k))-xx(1,d4(k)))
239 ccc(2) = xx(2,d4(k))+coef*(xx(2,d3(k))-xx(2,d4(k)))
240 ccc(3) = xx(3,d4(k))+coef*(xx(3,d3(k))-xx(3,d4(k)))
241 xk0(npoint) = fourth*(aaa(1)+bbb(1)+ccc(1)+xx(1,d4(k)))
242 yk0(npoint) = fourth*(aaa(2)+bbb(2)+ccc(2)+xx(2,d4(k)))
243 zk0(npoint) = fourth*(aaa(3)+bbb(3)+ccc(3)+xx(3,d4(k)))
244 ! LEVEL - 2 -> 4 SAMPLE POINTS ( 4 triangles on level 2 )
245 ! COEF = 3.0/8.0
246 coef = three*one_over_8
247 aaa(1) = xx(1,d4(k))+coef*(xx(1,d1(k))-xx(1,d4(k)))
248 aaa(2) = xx(2,d4(k))+coef*(xx(2,d1(k))-xx(2,d4(k)))
249 aaa(3) = xx(3,d4(k))+coef*(xx(3,d1(k))-xx(3,d4(k)))
250 bbb(1) = xx(1,d4(k))+coef*(xx(1,d2(k))-xx(1,d4(k)))
251 bbb(2) = xx(2,d4(k))+coef*(xx(2,d2(k))-xx(2,d4(k)))
252 bbb(3) = xx(3,d4(k))+coef*(xx(3,d2(k))-xx(3,d4(k)))
253 ccc(1) = xx(1,d4(k))+coef*(xx(1,d3(k))-xx(1,d4(k)))
254 ccc(2) = xx(2,d4(k))+coef*(xx(2,d3(k))-xx(2,d4(k)))
255 ccc(3) = xx(3,d4(k))+coef*(xx(3,d3(k))-xx(3,d4(k)))
256 npoint = npoint + 1
257 xk0(npoint) = third*(aaa(1)+bbb(1)+ccc(1))
258 yk0(npoint) = third*(aaa(2)+bbb(2)+ccc(2))
259 zk0(npoint) = third*(aaa(3)+bbb(3)+ccc(3))
260 l12(1,1) = half*(aaa(1)+bbb(1))
261 l12(2,1) = half*(aaa(2)+bbb(2))
262 l12(3,1) = half*(aaa(3)+bbb(3))
263 l23(1,1) = half*(bbb(1)+ccc(1))
264 l23(2,1) = half*(bbb(2)+ccc(2))
265 l23(3,1) = half*(bbb(3)+ccc(3))
266 l31(1,1) = half*(ccc(1)+aaa(1))
267 l31(2,1) = half*(ccc(2)+aaa(2))
268 l31(3,1) = half*(ccc(3)+aaa(3))
269 npoint = npoint + 1
270 xk0(npoint) = third*(aaa(1)+l12(1,1)+l31(1,1))
271 yk0(npoint) = third*(aaa(2)+l12(2,1)+l31(2,1))
272 zk0(npoint) = third*(aaa(3)+l12(3,1)+l31(3,1))
273 npoint = npoint + 1
274 xk0(npoint) = third*(bbb(1)+l12(1,1)+l23(1,1))
275 yk0(npoint) = third*(bbb(2)+l12(2,1)+l23(2,1))
276 zk0(npoint) = third*(bbb(3)+l12(3,1)+l23(3,1))
277 npoint = npoint + 1
278 xk0(npoint) = third*(ccc(1)+l23(1,1)+l31(1,1))
279 yk0(npoint) = third*(ccc(2)+l23(2,1)+l31(2,1))
280 zk0(npoint) = third*(ccc(3)+l23(3,1)+l31(3,1))
281 ! LEVEL - 3 -> 9 SAMPLE POINTS ( 9 triangles on level 3 )
282 ! COEF = 5.0/8.0
283 coef = five*one_over_8
284 aaa(1) = xx(1,d4(k))+coef*(xx(1,d1(k))-xx(1,d4(k)))
285 aaa(2) = xx(2,d4(k))+coef*(xx(2,d1(k))-xx(2,d4(k)))
286 aaa(3) = xx(3,d4(k))+coef*(xx(3,d1(k))-xx(3,d4(k)))
287 bbb(1) = xx(1,d4(k))+coef*(xx(1,d2(k))-xx(1,d4(k)))
288 bbb(2) = xx(2,d4(k))+coef*(xx(2,d2(k))-xx(2,d4(k)))
289 bbb(3) = xx(3,d4(k))+coef*(xx(3,d2(k))-xx(3,d4(k)))
290 ccc(1) = xx(1,d4(k))+coef*(xx(1,d3(k))-xx(1,d4(k)))
291 ccc(2) = xx(2,d4(k))+coef*(xx(2,d3(k))-xx(2,d4(k)))
292 ccc(3) = xx(3,d4(k))+coef*(xx(3,d3(k))-xx(3,d4(k)))
293 cg(1) = third*(aaa(1)+bbb(1)+ccc(1))
294 cg(2) = third*(aaa(2)+bbb(2)+ccc(2))
295 cg(3) = third*(aaa(3)+bbb(3)+ccc(3))
296 l12(1,1) = third*(two*aaa(1)+bbb(1))
297 l12(2,1) = third*(two*aaa(2)+bbb(2))
298 l12(3,1) = third*(two*aaa(3)+bbb(3))
299 l12(1,2) = third*(aaa(1)+two*bbb(1))
300 l12(2,2) = third*(aaa(2)+two*bbb(2))
301 l12(3,2) = third*(aaa(3)+two*bbb(3))
302 l23(1,1) = third*(two*bbb(1)+ccc(1))
303 l23(2,1) = third*(two*bbb(2)+ccc(2))
304 l23(3,1) = third*(two*bbb(3)+ccc(3))
305 l23(1,2) = third*(bbb(1)+two*ccc(1))
306 l23(2,2) = third*(bbb(2)+two*ccc(2))
307 l23(3,2) = third*(bbb(3)+two*ccc(3))
308 l31(1,1) = third*(two*ccc(1)+aaa(1))
309 l31(2,1) = third*(two*ccc(2)+aaa(2))
310 l31(3,1) = third*(two*ccc(3)+aaa(3))
311 l31(1,2) = third*(ccc(1)+two*aaa(1))
312 l31(2,2) = third*(ccc(2)+two*aaa(2))
313 l31(3,2) = third*(ccc(3)+two*aaa(3))
314 npoint = npoint + 1
315 xk0(npoint) = third*(aaa(1)+l12(1,1)+l31(1,2))
316 yk0(npoint) = third*(aaa(2)+l12(2,1)+l31(2,2))
317 zk0(npoint) = third*(aaa(3)+l12(3,1)+l31(3,2))
318 npoint = npoint + 1
319 xk0(npoint) = third*(bbb(1)+l12(1,2)+l23(1,1))
320 yk0(npoint) = third*(bbb(2)+l12(2,2)+l23(2,1))
321 zk0(npoint) = third*(bbb(3)+l12(3,2)+l23(3,1))
322 npoint = npoint + 1
323 xk0(npoint) = third*(ccc(1)+l23(1,2)+l31(1,1))
324 yk0(npoint) = third*(ccc(2)+l23(2,2)+l31(2,1))
325 zk0(npoint) = third*(ccc(3)+l23(3,2)+l31(3,1))
326 npoint = npoint + 1
327 xk0(npoint) = third*(cg(1)+l12(1,1)+l31(1,2))
328 yk0(npoint) = third*(cg(2)+l12(2,1)+l31(2,2))
329 zk0(npoint) = third*(cg(3)+l12(3,1)+l31(3,2))
330 npoint = npoint + 1
331 xk0(npoint) = third*(cg(1)+l12(1,2)+l23(1,1))
332 yk0(npoint) = third*(cg(2)+l12(2,2)+l23(2,1))
333 zk0(npoint) = third*(cg(3)+l12(3,2)+l23(3,1))
334 npoint = npoint + 1
335 xk0(npoint) = third*(cg(1)+l23(1,2)+l31(1,1))
336 yk0(npoint) = third*(cg(2)+l23(2,2)+l31(2,1))
337 zk0(npoint) = third*(cg(3)+l23(3,2)+l31(3,1))
338 npoint = npoint + 1
339 xk0(npoint) = third*(cg(1)+l12(1,1)+l12(1,2))
340 yk0(npoint) = third*(cg(2)+l12(2,1)+l12(2,2))
341 zk0(npoint) = third*(cg(3)+l12(3,1)+l12(3,2))
342 npoint = npoint + 1
343 xk0(npoint) = third*(cg(1)+l23(1,1)+l23(1,2))
344 yk0(npoint) = third*(cg(2)+l23(2,1)+l23(2,2))
345 zk0(npoint) = third*(cg(3)+l23(3,1)+l23(3,2))
346 npoint = npoint + 1
347 xk0(npoint) = third*(cg(1)+l31(1,1)+l31(1,2))
348 yk0(npoint) = third*(cg(2)+l31(2,1)+l31(2,2))
349 zk0(npoint) = third*(cg(3)+l31(3,1)+l31(3,2))
350 ! LEVEL - 4 -> 16 SAMPLE POINTS ( 16 triangles on level 4 )
351 ! COEF = 7.0/8.0
352 coef = seven*one_over_8
353 aaa(1) = xx(1,d4(k))+coef*(xx(1,d1(k))-xx(1,d4(k)))
354 aaa(2) = xx(2,d4(k))+coef*(xx(2,d1(k))-xx(2,d4(k)))
355 aaa(3) = xx(3,d4(k))+coef*(xx(3,d1(k))-xx(3,d4(k)))
356 bbb(1) = xx(1,d4(k))+coef*(xx(1,d2(k))-xx(1,d4(k)))
357 bbb(2) = xx(2,d4(k))+coef*(xx(2,d2(k))-xx(2,d4(k)))
358 bbb(3) = xx(3,d4(k))+coef*(xx(3,d2(k))-xx(3,d4(k)))
359 ccc(1) = xx(1,d4(k))+coef*(xx(1,d3(k))-xx(1,d4(k)))
360 ccc(2) = xx(2,d4(k))+coef*(xx(2,d3(k))-xx(2,d4(k)))
361 ccc(3) = xx(3,d4(k))+coef*(xx(3,d3(k))-xx(3,d4(k)))
362 cg(1) = third*(aaa(1)+bbb(1)+ccc(1))
363 cg(2) = third*(aaa(2)+bbb(2)+ccc(2))
364 cg(3) = third*(aaa(3)+bbb(3)+ccc(3))
365 npoint = npoint + 1
366 xk0(npoint) = cg(1)
367 yk0(npoint) = cg(2)
368 zk0(npoint) = cg(3)
369
370 l12(1,1) = fourth*(three*aaa(1)+bbb(1))
371 l12(2,1) = fourth*(three*aaa(2)+bbb(2))
372 l12(3,1) = fourth*(three*aaa(3)+bbb(3))
373 l12(1,2) = half*(aaa(1)+bbb(1))
374 l12(2,2) = half*(aaa(2)+bbb(2))
375 l12(3,2) = half*(aaa(3)+bbb(3))
376 l12(1,3) = fourth*(aaa(1)+three*bbb(1))
377 l12(2,3) = fourth*(aaa(2)+three*bbb(2))
378 l12(3,3) = fourth*(aaa(3)+three*bbb(3))
379 l23(1,1) = fourth*(three*bbb(1)+ccc(1))
380 l23(2,1) = fourth*(three*bbb(2)+ccc(2))
381 l23(3,1) = fourth*(three*bbb(3)+ccc(3))
382 l23(1,2) = half*(bbb(1)+ccc(1))
383 l23(2,2) = half*(bbb(2)+ccc(2))
384 l23(3,2) = half*(bbb(3)+ccc(3))
385 l23(1,3) = fourth*(bbb(1)+three*ccc(1))
386 l23(2,3) = fourth*(bbb(2)+three*ccc(2))
387 l23(3,3) = fourth*(bbb(3)+three*ccc(3))
388 l31(1,1) = fourth*(three*ccc(1)+aaa(1))
389 l31(2,1) = fourth*(three*ccc(2)+aaa(2))
390 l31(3,1) = fourth*(three*ccc(3)+aaa(3))
391 l31(1,2) = half*(ccc(1)+aaa(1))
392 l31(2,2) = half*(ccc(2)+aaa(2))
393 l31(3,2) = half*(ccc(3)+aaa(3))
394 l31(1,3) = fourth*(ccc(1)+three*aaa(1))
395 l31(2,3) = fourth*(ccc(2)+three*aaa(2))
396 l31(3,3) = fourth*(ccc(3)+three*aaa(3))
397 ll(1,1) = half*(l12(1,2)+l31(1,2))
398 ll(2,1) = half*(l12(2,2)+l31(2,2))
399 ll(3,1) = half*(l12(3,2)+l31(3,2))
400 ll(1,2) = half*(l12(1,2)+l23(1,2))
401 ll(2,2) = half*(l12(2,2)+l23(2,2))
402 ll(3,2) = half*(l12(3,2)+l23(3,2))
403 ll(1,3) = half*(l23(1,2)+l12(1,2))
404 ll(2,3) = half*(l23(2,2)+l12(2,2))
405 ll(3,3) = half*(l23(3,2)+l12(3,2))
406 npoint = npoint + 1
407 xk0(npoint) = third*(aaa(1)+l12(1,1)+l31(1,3))
408 yk0(npoint) = third*(aaa(2)+l12(2,1)+l31(2,3))
409 zk0(npoint) = third*(aaa(3)+l12(3,1)+l31(3,3))
410 npoint = npoint + 1
411 xk0(npoint) = third*(bbb(1)+l12(1,3)+l23(1,1))
412 yk0(npoint) = third*(bbb(2)+l12(2,3)+l23(2,1))
413 zk0(npoint) = third*(bbb(3)+l12(3,3)+l23(3,1))
414 npoint = npoint + 1
415 xk0(npoint) = third*(ccc(1)+l23(1,3)+l31(1,1))
416 yk0(npoint) = third*(ccc(2)+l23(2,3)+l31(2,1))
417 zk0(npoint) = third*(ccc(3)+l23(3,3)+l31(3,1))
418 npoint = npoint + 1
419 xk0(npoint) = third*(ll(1,1)+l12(1,1)+l31(1,3))
420 yk0(npoint) = third*(ll(2,1)+l12(2,1)+l31(2,3))
421 zk0(npoint) = third*(ll(3,1)+l12(3,1)+l31(3,3))
422 npoint = npoint + 1
423 xk0(npoint) = third*(ll(1,2)+l12(1,3)+l23(1,1))
424 yk0(npoint) = third*(ll(2,2)+l12(2,3)+l23(2,1))
425 zk0(npoint) = third*(ll(3,2)+l12(3,3)+l23(3,1))
426 npoint = npoint + 1
427 xk0(npoint) = third*(ll(1,3)+l23(1,3)+l31(1,1))
428 yk0(npoint) = third*(ll(2,3)+l23(2,3)+l31(2,1))
429 zk0(npoint) = third*(ll(3,3)+l23(3,3)+l31(3,1))
430 npoint = npoint + 1
431 xk0(npoint) = third*(ll(1,1)+l12(1,1)+l12(1,2))
432 yk0(npoint) = third*(ll(2,1)+l12(2,1)+l12(2,2))
433 zk0(npoint) = third*(ll(3,1)+l12(3,1)+l12(3,2))
434 npoint = npoint + 1
435 xk0(npoint) = third*(ll(1,2)+l12(1,2)+l12(1,3))
436 yk0(npoint) = third*(ll(2,2)+l12(2,2)+l12(2,3))
437 zk0(npoint) = third*(ll(3,2)+l12(3,2)+l12(3,3))
438 npoint = npoint + 1
439 xk0(npoint) = third*(ll(1,2)+l23(1,1)+l23(1,2))
440 yk0(npoint) = third*(ll(2,2)+l23(2,1)+l23(2,2))
441 zk0(npoint) = third*(ll(3,2)+l23(3,1)+l23(3,2))
442 npoint = npoint + 1
443 xk0(npoint) = third*(ll(1,3)+l23(1,2)+l23(1,3))
444 yk0(npoint) = third*(ll(2,3)+l23(2,2)+l23(2,3))
445 zk0(npoint) = third*(ll(3,3)+l23(3,2)+l23(3,3))
446 npoint = npoint + 1
447 xk0(npoint) = third*(ll(1,3)+l31(1,1)+l31(1,2))
448 yk0(npoint) = third*(ll(2,3)+l31(2,1)+l31(2,2))
449 zk0(npoint) = third*(ll(3,3)+l31(3,1)+l31(3,2))
450 npoint = npoint + 1
451 xk0(npoint) = third*(ll(1,1)+l31(1,2)+l31(1,3))
452 yk0(npoint) = third*(ll(2,1)+l31(2,2)+l31(2,3))
453 zk0(npoint) = third*(ll(3,1)+l31(3,2)+l31(3,3))
454 npoint = npoint + 1
455 xk0(npoint) = third*(ll(1,1)+ll(1,2)+l12(1,2))
456 yk0(npoint) = third*(ll(2,1)+ll(2,2)+l12(2,2))
457 zk0(npoint) = third*(ll(3,1)+ll(3,2)+l12(3,2))
458 npoint = npoint + 1
459 xk0(npoint) = third*(ll(1,2)+ll(1,3)+l23(1,2))
460 yk0(npoint) = third*(ll(2,2)+ll(2,3)+l23(2,2))
461 zk0(npoint) = third*(ll(3,2)+ll(3,3)+l23(3,2))
462 npoint = npoint + 1
463 xk0(npoint) = third*(ll(1,3)+ll(1,1)+l31(1,2))
464 yk0(npoint) = third*(ll(2,3)+ll(2,1)+l31(2,2))
465 zk0(npoint) = third*(ll(3,3)+ll(3,1)+l31(3,2))
466 ENDDO ! DO K=1,4
467
468 DO j=1,npoint
469 xs(j,i)=xk0(j)
470 ys(j,i)=yk0(j)
471 zs(j,i)=zk0(j)
472 ENDDO
473 ENDDO ! DO I=1,NEL
474
475 ! NTRACE_TOT = NPOINT
476 ELSEIF (isolnod == 8) THEN
477
478 DO i=1,nel
479 xx(1,1)=x1(i)
480 xx(2,1)=y1(i)
481 xx(3,1)=z1(i)
482 xx(1,2)=x2(i)
483 xx(2,2)=y2(i)
484 xx(3,2)=z2(i)
485 xx(1,3)=x3(i)
486 xx(2,3)=y3(i)
487 xx(3,3)=z3(i)
488 xx(1,4)=x4(i)
489 xx(2,4)=y4(i)
490 xx(3,4)=z4(i)
491 xx(1,5)=x5(i)
492 xx(2,5)=y5(i)
493 xx(3,5)=z5(i)
494 xx(1,6)=x6(i)
495 xx(2,6)=y6(i)
496 xx(3,6)=z6(i)
497 xx(1,7)=x7(i)
498 xx(2,7)=y7(i)
499 xx(3,7)=z7(i)
500 xx(1,8)=x8(i)
501 xx(2,8)=y8(i)
502 xx(3,8)=z8(i)
503
504 IF (ipart_(i) /= idp) cycle
505
506 xmin = ep20
507 ymin = ep20
508 zmin = ep20
509 xmax =-ep20
510 ymax =-ep20
511 zmax =-ep20
512
513 DO j=1,8
514 xmin=min(xmin,xx(1,j))
515 ymin=min(ymin,xx(2,j))
516 zmin=min(zmin,xx(3,j))
517 xmax=max(xmax,xx(1,j))
518 ymax=max(ymax,xx(2,j))
519 zmax=max(zmax,xx(3,j))
520 END DO
521
522 dx = (xmax-xmin)/float(ntrace0)
523 dy = (ymax-ymin)/float(ntrace0)
524 dz = (zmax-zmin)/float(ntrace0)
525
526 n1 = ntrace0
527 n2 = ntrace0
528 n3 = ntrace0
529
530 xk(1) = xmin + dx*half
531 yk(1) = ymin + dy*half
532 zk(1) = zmin + dz*half
533
534 DO k1=2,n1
535 xk(k1) = xk(k1-1) + dx
536 yk(k1) = yk(k1-1) + dy
537 zk(k1) = zk(k1-1) + dz
538 ENDDO
539
540 j=0
541 DO k3=1,n3
542 DO k2=1,n2
543 DO k1=1,n1
544 j=j+1
545 xs(j,i)=xk(k1)
546 ys(j,i)=yk(k2)
547 zs(j,i)=zk(k3)
548 ENDDO ! DO K1=1,N1
549 ENDDO ! DO K2=1,N2
550 ENDDO ! DO K3=1,N3
551 ENDDO ! DO I=1,NEL
552 ENDIF ! IF (ISOLNOD == 4)
553
554 IF (isolnod == 4 .OR. n2d/=0) THEN
555 ntrace_tot = npoint
556 ELSEIF (isolnod == 8) THEN
557 ntrace_tot = ntrace
558 ENDIF
559
560 DO ii=1,getel
561 i=jct(ii)
562 IF (ipart_(i) /= idp) cycle
563
564 DO ip=1,ntrace_tot
565 inod = 0
566 dist = zero
567 dist_old = ep20
568
569 IF(surf_type == 101)THEN
570 !--ellipsoid
571 iad0 = iad_bufr
572 dist = zero
573 !get this out from the loop ! does not depends on IP
574 aa = bufsf(iad0+1)
575 bb = bufsf(iad0+2)
576 cc = bufsf(iad0+3)
577 xg = bufsf(iad0+4)
578 yg = bufsf(iad0+5)
579 zg = bufsf(iad0+6)
580 skw(1)=bufsf(iad0+7)
581 skw(2)=bufsf(iad0+8)
582 skw(3)=bufsf(iad0+9)
583 skw(4)=bufsf(iad0+10)
584 skw(5)=bufsf(iad0+11)
585 skw(6)=bufsf(iad0+12)
586 skw(7)=bufsf(iad0+13)
587 skw(8)=bufsf(iad0+14)
588 skw(9)=bufsf(iad0+15)
589 dgr=bufsf(iad0+36)
590 ! transition matrix
591 x_prime = skw(1)*(xs(ip,i)-xg) + skw(4)*(ys(ip,i)-yg) + skw(7)*(zs(ip,i)-zg)
592 y_prime = skw(2)*(xs(ip,i)-xg) + skw(5)*(ys(ip,i)-yg) + skw(8)*(zs(ip,i)-zg)
593 z_prime = skw(3)*(xs(ip,i)-xg) + skw(6)*(ys(ip,i)-yg) + skw(9)*(zs(ip,i)-zg)
594 tmp(1)= abs(x_prime)/aa
595 tmp(2)= abs(y_prime)/bb
596 tmp(3)= abs(z_prime)/cc
597 IF(tmp(1)/=zero)tmp(1)= exp(dgr*log(tmp(1)))
598 IF(tmp(2)/=zero)tmp(2)= exp(dgr*log(tmp(2)))
599 IF(tmp(3)/=zero)tmp(3)= exp(dgr*log(tmp(3)))
600 dist = (tmp(1)+tmp(2)+tmp(3))
601 disp(ip,i) = one-dist
602
603 ELSEIF (surf_type == 200) THEN
604 !-Planar Surface
605 !get this out from the loop ! does not depends on IP
606 iad0 = iad_bufr
607 dist = zero
608
609 xp1 = bufsf(iad0+1)
610 yp1 = bufsf(iad0+2)
611 zp1 = bufsf(iad0+3)
612 xp2 = bufsf(iad0+4)
613 yp2 = bufsf(iad0+5)
614 zp2 = bufsf(iad0+6)
615
616 aa = xp2 - xp1
617 bb = yp2 - yp1
618 cc = zp2 - zp1
619
620 dist = aa*(xs(ip,i)-xp1)+bb*(ys(ip,i)-yp1)+cc*(zs(ip,i)-zp1)
621 tmpsum = sqrt(aa*aa+bb*bb+cc*cc)
622 tmpsum = one/max(em30,tmpsum)
623 dist = dist*tmpsum
624 disp(ip,i) = dist
625
626 ELSE
627 !--surface of eleme,ts
628 IF(n2d == 0)THEN
629 IF (isolnod == 4) THEN
630 ix(1) =ixs(2,i)
631 ix(2) =ixs(4,i)
632 ix(3) =ixs(7,i)
633 ix(4) =ixs(6,i)
634 elem_numnod = 4
635 ELSEIF (isolnod == 8) THEN
636 ix(1) =ixs(2,i)
637 ix(2) =ixs(3,i)
638 ix(3) =ixs(4,i)
639 ix(4) =ixs(5,i)
640 ix(5) =ixs(6,i)
641 ix(6) =ixs(7,i)
642 ix(7) =ixs(8,i)
643 ix(8) =ixs(9,i)
644 elem_numnod = 8
645 ENDIF
646 ELSE
647 IF(ityp == 7)THEN
648 ix(1) =ixtg(2,i)
649 ix(2) =ixtg(3,i)
650 ix(3) =ixtg(4,i)
651 ix(4) =0
652 elem_numnod = 3
653 ELSEIF(ityp == 2)THEN
654 ix(1) =ixq(2,i)
655 ix(2) =ixq(3,i)
656 ix(3) =ixq(4,i)
657 ix(4) =ixq(5,i)
658 elem_numnod = 4
659 ENDIF
660 ENDIF
661
662 DO j=1,elem_numnod
663 n = ix(j)
664! NSOLTOSF(IDC,N) - the shell node of the container IDC, to which
665! the distance of the eulerian node N to the shell
666! container is minimum
667 nsh = nsoltosf(idc,n)
668 IF (nsh <= 0) cycle
669! KNOD2SURF(NSH) - the nb of surfaces of the current container IDC,
670! connected to node NSH
671 DO jj=1,knod2surf(nsh)
672! INOD2SURF(JJ,NSH) - identify the segment of the current surface
673! container
674 ipl = inod2surf(jj,nsh)
675! SEGTOSURF(IPL) - identify the right surface to whom the segment
676! IPL belongs. One node NSH can connect segments coming from
677! different containers surface.
678! Do not forget that containers overwrite phases (superpose)
679 isurf = segtosurf(ipl)
680 ipl = ipl - swiftsurf(isurf)
681 IF (ipl <= 0 .OR. ipl > nseg) cycle
682!
683 ixpl(1) = igrsurf(isurf)%NODES(ipl,1)
684 ixpl(2) = igrsurf(isurf)%NODES(ipl,2)
685 ixpl(3) = igrsurf(isurf)%NODES(ipl,3)
686 ixpl(4) = igrsurf(isurf)%NODES(ipl,4)
687!
688 xfas(1,1) = x(1,ixpl(1))
689 xfas(2,1) = x(2,ixpl(1))
690 xfas(3,1) = x(3,ixpl(1))
691 xfas(1,2) = x(1,ixpl(2))
692 xfas(2,2) = x(2,ixpl(2))
693 xfas(3,2) = x(3,ixpl(2))
694 xfas(1,3) = x(1,ixpl(3))
695 xfas(2,3) = x(2,ixpl(3))
696 xfas(3,3) = x(3,ixpl(3))
697 xfas(1,4) = x(1,ixpl(4))
698 xfas(2,4) = x(2,ixpl(4))
699 xfas(3,4) = x(3,ixpl(4))
700
701 DO k=1,4
702 ! compute min dist from trace sample points to cutting container:
703 x0 = xs(ip,i) - xfas(1,k)
704 y0 = ys(ip,i) - xfas(2,k)
705 z0 = zs(ip,i) - xfas(3,k)
706 dist = x0*x0 + y0*y0 + z0*z0
707 dist = sqrt(dist)
708 IF (dist < dist_old .and. dist > em10) THEN
709 dist_old = dist
710 inod = ixpl(k)
711 ENDIF
712 ENDDO ! DO K=1,4
713 ENDDO ! DO II=1,KNOD2SURF(NSH)
714 ENDDO ! DO J=1,8
715
716 IF (inod == 0) GOTO 122
717
718 IF (dist_old == ep20) dist_old = zero
719 disp(ip,i) = dist_old
720
721 ! get sig n of dist of trace sample points
722 xn(1) = xs(ip,i)
723 xn(2) = ys(ip,i)
724 xn(3) = zs(ip,i)
725 dist = zero
726 CALL in_out_side(
727 . inod ,inod2surf ,knod2surf ,nnod2surf ,x ,
728 . xn ,dist ,nseg ,surf_eltyp,nod_normal,
729 . surf_nodes,swiftsurf ,idsurf ,segtosurf )
730 disp(ip,i) = dist
731
732 ENDIF
733
734 ! START counting trace samples and filling process
735 jmid_old = inphase(ip,i)
736 IF (disp(ip,i) > zero) THEN
737 tracep(i) = tracep(i) + 1
738 IF (ifill == 0) THEN
739 IF (jmid_old /= jmid) THEN
740 inphase(ip,i) = jmid ! get new phase
741 nbip(jmid_old,i) = nbip(jmid_old,i) - 1
742 nbip(jmid,i) = nbip(jmid,i) + 1
743 ENDIF
744 END IF ! IF (IFILL == 0)
745 ELSEIF (disp(ip,i) < zero) THEN
746 tracen(i) = tracen(i) + 1
747 IF (ifill == 1) THEN
748 IF (jmid_old /= jmid) THEN
749 inphase(ip,i) = jmid ! get new phase
750 nbip(jmid_old,i) = nbip(jmid_old,i) - 1
751 nbip(jmid,i) = nbip(jmid,i) + 1
752 ENDIF
753 ENDIF ! IF (IFILL == 1)
754 ELSEIF (disp(ip,i) == zero .and. jmid /= jmid_old) THEN
755 ! eliminate phase within trace sample points on container surface
756 nbip(jmid_old,i) = nbip(jmid_old,i) - 1
757 inphase(ip,i) = 0
758 ENDIF ! IF (DISP(IP,I) > ZERO)
759
760 122 CONTINUE
761
762 ENDDO ! IP=1,NTRACE_TOT
763 ! check / remove old non-existing phase within element :
764 IF (jmid_old > 0)THEN
765 IF(nbip(jmid_old,i) < 0) nbip(jmid_old,i) = 0
766 ENDIF
767
768 ok = 0
769 nphase = iphase(nbsubmat+1,i)
770 k = nphase
771 DO j=1,nphase
772 IF (jmid /= iphase(j,i)) ok = ok + 1
773 ENDDO
774 IF (ok == k) THEN
775 k = k + 1
776 iphase(k,i) = jmid
777 iphase(nbsubmat+1,i) = k ! increasing the nb of phases within cut element
778 ENDIF
779
780 IF (tracep(i) <= 0 .and. tracen(i) > 0) THEN
781 full(i) = -1
782 ELSEIF (tracep(i) > 0 .and. tracen(i) <= 0) THEN
783 full(i) = 1
784 ENDIF
785
786 ENDDO ! DO II=1,GETEL
787
788 ! recompute / reset / overwrite phases inside solid element
789 DO i=1,nel
790 IF (ipart_(i) /= idp) cycle
791 IF (full(i) == 1 .and. ifill == 0) THEN
792 kvol_bak(1:nbsubmat) = kvol(1:nbsubmat,i)
793 DO j=1,nbsubmat
794 iphase(j,i) = 0
795 IF(icumu == 0)kvol(j,i) = zero
796 nbip(j,i) = 0
797 ENDDO
798 iphase(1,i) = jmid
799 iphase(nbsubmat+1,i) = 1
800 nbip(jmid,i) = ntrace_tot
801 DO ip=1,ntrace_tot
802 inphase(ip,i) = jmid
803 ENDDO
804 IF(icumu == 0)kvol(jmid,i)=zero
805 kvol(jmid,i)= kvol(jmid,i)+fill_ratio
806 ! if added volume ratio makes that sum is > 1, then substract from previous filling
807 IF(icumu == -1)THEN
808 ! if added volume ratio makes that sum is > 1, then substract from previous filling
809 sumvf = sum(kvol(1:nbsubmat, i))
810 IF (sumvf > one)THEN
811 IF(idc > 1)THEN
812 ! SUBSTRACT FROM PREVIOUS STEP
813 isubmat_to_substract = inivol(i_inivol)%CONTAINER(idc-1)%SUBMAT_ID
814 vf_to_substract = sumvf-one
815 vf_to_substract = min(vf_to_substract, kvol(isubmat_to_substract,i))
816 kvol(isubmat_to_substract,i) = kvol(isubmat_to_substract,i) - vf_to_substract
817 ELSE
818 ! SUBSTRACT FROM DEFAULT SUBMATERIAL
819 isubmat_to_substract = maxloc(kvol_bak,1)
820 vf_to_substract = sumvf-one
821 vf_to_substract = min(vf_to_substract, kvol(isubmat_to_substract,i))
822 kvol(isubmat_to_substract,i) = kvol(isubmat_to_substract,i) - vf_to_substract
823 ENDIF
824 END IF
825 ENDIF
826 ELSEIF (full(i) == -1 .and. ifill == 1) THEN
827 kvol_bak(1:nbsubmat) = kvol(1:nbsubmat,i)
828 DO j=1,nbsubmat
829 iphase(j,i) = 0
830 IF(icumu == 0)kvol(j,i) = zero
831 nbip(j,i) = 0
832 ENDDO
833 iphase(1,i) = jmid
834 iphase(nbsubmat+1,i) = 1
835 nbip(jmid,i) = ntrace_tot
836 DO ip=1,ntrace_tot
837 inphase(ip,i) = jmid
838 ENDDO
839 IF(icumu == 0)kvol(jmid,i)=zero
840 kvol(jmid,i)= kvol(jmid,i)+fill_ratio
841 IF(icumu == -1)THEN
842 ! if added volume ratio makes that sum is > 1, then substract from previous filling
843 sumvf = sum(kvol(1:nbsubmat, i))
844 IF (sumvf > one)THEN
845 IF(idc > 1)THEN
846 ! SUBSTRACT FROM PREVIOUS STEP
847 isubmat_to_substract = inivol(i_inivol)%CONTAINER(idc-1)%SUBMAT_ID
848 vf_to_substract = sumvf-one
849 vf_to_substract = min(vf_to_substract, kvol(isubmat_to_substract,i))
850 kvol(isubmat_to_substract,i) = kvol(isubmat_to_substract,i) - vf_to_substract
851 ELSE
852 ! SUBSTRACT FROM DEFAULT SUBMATERIAL
853 isubmat_to_substract = maxloc(kvol_bak,1)
854 vf_to_substract = sumvf-one
855 vf_to_substract = min(vf_to_substract, kvol(isubmat_to_substract,i))
856 kvol(isubmat_to_substract,i) = kvol(isubmat_to_substract,i) - vf_to_substract
857 ENDIF
858 END IF
859 ENDIF
860 ELSEIF (full(i) == 2) THEN
861 DO j=1,nbsubmat
862 buffill1(j) = 0
863 buffill2(j,i)= 0
864 ENDDO
865
866 DO j=1,iphase(nbsubmat+1,i)
867 iph = iphase(j,i)
868 IF (iph /= 0) THEN
869 IF (nbip(iph,i) == 0) iphase(j,i) = 0
870 ENDIF
871 ENDDO
872
873 k = 0
874 ok = 0
875 DO j=1,iphase(nbsubmat+1,i)
876 IF (iphase(j,i) /= 0) THEN
877 iph = iphase(j,i)
878 k = k + 1
879 buffill1(k) = iphase(j,i)
880 buffill2(k,i)= nbip(iph,i)
881 ENDIF
882 ENDDO
883
884 IF (iphase(nbsubmat+1,i) > 1) THEN
885 DO j=1,nbsubmat
886 iphase(j,i) = 0
887 nbip(j,i) = 0
888 ENDDO
889
890 DO j=1,k
891 iphase(j,i) = buffill1(j)
892 iph = iphase(j,i)
893 nbip(iph,i) = buffill2(j,i)
894 ENDDO
895 iphase(nbsubmat+1,i) = k
896 ENDIF
897 ENDIF ! IF (FULL(I) == 1 .and. IFILL == 0)
898
899
900 ENDDO ! DO I=1,NEL
901
902 ! COMPUTE VOLUME FRACTION
903 DO i=1,nel
904 IF (ipart_(i) /= idp) cycle
905 nip = 0
906 IF (iphase(nbsubmat+1,i) > 1) THEN
907 DO j=1,iphase(nbsubmat+1,i)
908 iph = iphase(j,i)
909 nip = nip + nbip(iph,i)
910 ENDDO
911 kvol_bak(1:nbsubmat) = kvol(1:nbsubmat,i)
912 DO j=1,iphase(nbsubmat+1,i)
913 iph = iphase(j,i)
914 IF(icumu == 0)kvol(iph,i)=zero
915 kvol(iph,i)= kvol(iph,i)+fill_ratio*float(nbip(iph,i))/float(nip)
916 IF(icumu == -1)THEN
917 ! if added volume ratio makes that sum is > 1, then substract from previous filling
918 sumvf = sum(kvol(1:nbsubmat, i))
919 IF (sumvf > one )THEN
920 IF(idc > 1)THEN
921 ! SUBSTRACT FROM PREVIOUS STEP
922 isubmat_to_substract = inivol(i_inivol)%CONTAINER(idc-1)%SUBMAT_ID
923 vf_to_substract = sumvf-one
924 vf_to_substract = min(vf_to_substract, kvol(isubmat_to_substract,i))
925 kvol(isubmat_to_substract,i) = kvol(isubmat_to_substract,i) - vf_to_substract
926 ELSE
927 ! SUBSTRACT FROM PREVIOUS STEP
928 isubmat_to_substract = maxloc(kvol_bak,1)
929 vf_to_substract = sumvf-one
930 vf_to_substract = min(vf_to_substract, kvol(isubmat_to_substract,i))
931 kvol(isubmat_to_substract,i) = kvol(isubmat_to_substract,i) - vf_to_substract
932 ENDIF
933 END IF
934 ENDIF
935 ENDDO
936 ELSE
937 ENDIF
938 ENDDO ! DO I=1,NEL
939!------------------
940 RETURN
941 END
cg
subroutine cg(dim, mat, rhs, sol, max_iter, tol)
Definition conjugate_gradient.F:29
in_out_side
subroutine in_out_side(inod, inod2surf, knod2surf, nnod2surf, x, xn, dist, nseg, surf_eltyp, nod_normal, surf_nodes, swiftsurf, idsurf, segtosurf)
Definition in_out_side.F:35
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:274
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
groupdef_mod
Definition groupdef_mod.F:668
inivol_def_mod
Definition inivol_mod.F:45
ratio_fill
subroutine ratio_fill(x1, x2, x3, x4, x5, x6, x7, x8, y1, y2, y3, y4, y5, y6, y7, y8, z1, z2, z3, z4, z5, z6, z7, z8, idp, x, ixs, ipart_, ifill, ntrace, ntrace0, dis, nsoltosf, nnod2surf, inod2surf, knod2surf, jmid, iphase, inphase, kvol, surf_type, iad_bufr, bufsf, nod_normal, isolnod, nbsubmat, fill_ratio, icumu, nseg, surf_eltyp, surf_nodes, nbconty, idc, nbip, idsurf, swiftsurf, segtosurf, igrsurf, ivolsurf, nsurf_invol, ixq, ixtg, ityp, nel, numel_tot, num_inivol, inivol, i_inivol)
Definition ratio_fill.F:42
inivol_def_mod::inivol_struct_
Definition inivol_mod.F:70