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