File: RELATIVE:/../../../mfix.git/model/des/stl_preproc_des_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module name: des_stl_functions                                      !
4     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
5     !                                                                      !
6     !  Purpose: This module containd routines for geometric interaction    !
7     !  required for STL files.                                             !
8     !                                                                      !
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
10           MODULE STL_PREPROC_DES
11     
12           use des_stl_functions, only: TestTriangleAABB
13     
14           use desgrid
15           use error_manager
16     
17           IMPLICIT NONE
18     
19     ! Use this module only to define functions and subroutines.
20           CONTAINS
21     
22     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
23     !                                                                      !
24     !  Subroutine: des_stl_preprocessing                                   !
25     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
26     !                                                                      !
27     !  Purpose:                                                            !
28     !                                                                      !
29     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
30           SUBROUTINE DES_STL_PREPROCESSING
31     
32           USE cutcell, only: use_stl
33     
34           USE discretelement, only: STL_FACET_TYPE
35           USE discretelement, only: FACET_TYPE_NORMAL
36           USE discretelement, only: FACET_TYPE_MI
37           USE discretelement, only: FACET_TYPE_PO
38           USE discretelement, only: COUNT_FACET_TYPE_NORMAL
39           USE discretelement, only: COUNT_FACET_TYPE_MI
40           USE discretelement, only: COUNT_FACET_TYPE_PO
41     
42           USE desgrid
43     
44           USE param
45           USE stl
46     
47           implicit none
48     
49           INTEGER :: IJK, COUNT, NF
50     
51     ! Pre-procssing for the des in order to assign facets to grid cells.
52           WRITE(ERR_MSG,"('Pre-Processing geometry for DES.')")
53           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
54     
55           CALL ALLOCATE_DES_STL_ARRAYS
56     
57           N_FACETS_DES = 0
58     
59     ! Set N_facets_des to add any more facets needed by dem and not to.
60     ! contaminate the Eulerian-Eulerian CG stuff
61           IF(USE_STL) N_FACETS_DES = N_FACETS
62     
63     ! Set all the facets generated so far as normal facet types
64           STL_FACET_TYPE(1:N_FACETS) = FACET_TYPE_NORMAL
65     
66     
67           CALL BIN_FACETS_TO_GRID_DES
68     
69           DO IJK = 1, DG_IJKSIZE2
70              COUNT = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
71              IF(COUNT.eq.0) DEALLOCATE(LIST_FACET_AT_DES(IJK)%FACET_LIST)
72           ENDDO
73     
74           COUNT_FACET_TYPE_NORMAL = 0
75           COUNT_FACET_TYPE_PO     = 0
76           COUNT_FACET_TYPE_MI     = 0
77           DO NF = 1, N_FACETS_DES
78              IF(STL_FACET_TYPE(NF) == FACET_TYPE_NORMAL) &
79                 COUNT_FACET_TYPE_NORMAL = COUNT_FACET_TYPE_NORMAL + 1
80     
81              IF(STL_FACET_TYPE(NF).EQ.FACET_TYPE_MI) &
82                 COUNT_FACET_TYPE_MI = COUNT_FACET_TYPE_MI + 1
83     
84              IF(STL_FACET_TYPE(NF).EQ.FACET_TYPE_PO) &
85                 COUNT_FACET_TYPE_PO = COUNT_FACET_TYPE_PO + 1
86           ENDDO
87     
88     !      CALL DEBUG_WRITE_GRID_FACEINFO
89     !      CALL DEBUG_WRITE_STL_FROM_GRID_FACET(WRITE_FACETS_EACH_CELL=.FALSE.)
90     
91     !      CALL DEBUG_WRITE_ALL_READIN_FACETS
92     
93           RETURN
94           END SUBROUTINE DES_STL_PREPROCESSING
95     
96     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
97     !                                                                      !
98     !  Subroutine: ALLOCATE_DES_STL_ARRAYS                                 !
99     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
100     !                                                                      !
101     !  Purpose:                                                            !
102     !                                                                      !
103     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
104           SUBROUTINE ALLOCATE_DES_STL_ARRAYS
105     
106           USE desgrid
107           USE discretelement, only: STL_Facet_type
108           USE discretelement, only: CELLNEIGHBOR_FACET
109           USE discretelement, only: CELLNEIGHBOR_FACET_MAX
110           USE discretelement, only: CELLNEIGHBOR_FACET_NUM
111           USE discretelement, only: CELLNEIGHBOR_FACET_NUM
112           USE stl, only: LIST_FACET_AT_DES
113           USE stl, only: NO_NEIGHBORING_FACET_DES
114           USE param
115           USE stl, only: DIM_STL, MAX_FACETS_PER_CELL_DES
116     
117           IMPLICIT NONE
118     
119           integer :: ii,ijk
120     
121           ALLOCATE(LIST_FACET_AT_DES(DG_IJKSIZE2))
122     
123           Allocate(  CELLNEIGHBOR_FACET (DG_IJKSIZE2) )
124           Allocate(  CELLNEIGHBOR_FACET_MAX (DG_IJKSIZE2) )
125           Allocate(  CELLNEIGHBOR_FACET_NUM (DG_IJKSIZE2) )
126           DO II = 1, DG_IJKSIZE2
127              CELLNEIGHBOR_FACET_MAX(II) = 4
128              allocate(CELLNEIGHBOR_FACET(II)%P(CELLNEIGHBOR_FACET_MAX(II)))
129              allocate(CELLNEIGHBOR_FACET(II)%EXTENTDIR(CELLNEIGHBOR_FACET_MAX(II)))
130              allocate(CELLNEIGHBOR_FACET(II)%EXTENTMIN(CELLNEIGHBOR_FACET_MAX(II)))
131              allocate(CELLNEIGHBOR_FACET(II)%EXTENTMAX(CELLNEIGHBOR_FACET_MAX(II)))
132              CELLNEIGHBOR_FACET_NUM(II) = 0
133           ENDDO
134     
135           DO IJK = 1, DG_IJKSIZE2
136              LIST_FACET_AT_DES(IJK)%COUNT_FACETS = 0
137              ALLOCATE(LIST_FACET_AT_DES(IJK)%FACET_LIST(MAX_FACETS_PER_CELL_DES))
138           ENDDO
139     
140           ALLOCATE(NO_NEIGHBORING_FACET_DES(DG_IJKSIZE2))
141           NO_NEIGHBORING_FACET_DES = .false.
142     
143           Allocate(stl_facet_type(DIM_STL))
144     
145           END SUBROUTINE ALLOCATE_DES_STL_ARRAYS
146     
147     
148     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
149     !                                                                      !
150     !  Subroutine: bin_facets_to_grid_des                                  !
151     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
152     !                                                                      !
153     !  Purpose:                                                            !
154     !                                                                      !
155     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
156           Subroutine BIN_FACETS_TO_GRID_DES
157     
158           USE desgrid
159           USE discretelement, only: dimn
160           USE functions
161           USE geometry
162           USE indices
163           USE mpi_utility
164           USE param1, only: zero
165           Use stl, only: tol_stl, n_facets_des
166           Use stl, only: vertex, NO_NEIGHBORING_FACET_DES, LIST_FACET_AT_DES
167     
168           IMPLICIT NONE
169     
170     ! DES Grid cell index.
171           INTEGER :: IJK, IJK2
172     ! Loop counters:
173           INTEGER :: I, I1, I2, II  ! X-axis
174           INTEGER :: J, J1, J2, JJ  ! Y-axis
175           INTEGER :: K, K1, K2, KK  ! Z-axis
176           INTEGER :: N              ! STLs
177     ! Generic accumulator
178           INTEGER :: COUNT_FAC
179     
180     ! Maximum and minimum extents of the indexed STL
181           DOUBLE PRECISION:: X1,Y1,Z1
182           DOUBLE PRECISION:: X2,Y2,Z2
183     
184     
185     
186           DO N = 1,N_FACETS_DES
187     
188              X1 = minval(VERTEX(1:3,1,N))
189              X2 = maxval(VERTEX(1:3,1,N))
190              Y1 = minval(VERTEX(1:3,2,N))
191              Y2 = maxval(VERTEX(1:3,2,N))
192              Z1 = minval(VERTEX(1:3,3,N))
193              Z2 = maxval(VERTEX(1:3,3,N))
194     
195              I1 = DG_IEND2
196              I2 = DG_ISTART2
197              IF(X2>=ZERO .AND. X1<=XLENGTH+TOL_STL) THEN
198                 I1 = iofpos(X1)
199                 I2 = iofpos(1/dg_dxinv+X2)
200              ENDIF
201     
202              J1 = DG_JEND2
203              J2 = DG_JSTART2
204              IF(Y2>=ZERO .AND. Y1<=YLENGTH+TOL_STL) THEN
205                 J1 = jofpos(Y1)
206                 J2 = jofpos(1/dg_dyinv+Y2)
207              ENDIF
208     
209              K1 = DG_KEND2
210              K2 = DG_KSTART2
211              IF(DO_K) THEN
212                 IF(Z2>=ZERO .AND. Z1<=ZLENGTH+TOL_STL) THEN
213                    K1 = kofpos(Z1)
214                    K2 = kofpos(1/dg_dzinv+Z2)
215                 ENDIF
216              ENDIF
217     
218     
219              DO K=K1,K2
220              DO J=J1,J2
221              DO I=I1,I2
222                 IF(dg_is_ON_myPE_plus1layers(I,J,K)) THEN
223                    IJK = DG_FUNIJK(I,J,K)
224                    CALL ADD_FACET_FOR_DES(I,J,K,IJK,N)
225                 ENDIF
226              ENDDO
227              ENDDO
228              ENDDO
229     
230           ENDDO
231     
232     
233     ! JM: I don't think this is needed.
234           DO K = DG_KSTART1, DG_KEND1
235           DO J = DG_JSTART1, DG_JEND1
236           DO I = DG_ISTART1, DG_IEND1
237     
238              IJK = DG_FUNIJK(I, J, K)
239     
240              I1 =  MAX( I - 1, DG_ISTART2)
241              I2 =  MIN( I + 1, DG_IEND2)
242     
243              J1 =  MAX( J - 1, DG_JSTART2)
244              J2 =  MIN (J + 1, DG_JEND2)
245     
246              K1 =  MAX( K - 1, DG_KSTART2)
247              K2 =  MIN (K + 1, DG_KEND2)
248     
249              NO_NEIGHBORING_FACET_DES(IJK)  = .FALSE.
250     
251              COUNT_FAC = 0
252              DO KK = K1, k2
253              DO JJ = J1, j2
254              DO II = I1, i2
255                 IJK2 = DG_FUNIJK(II, JJ, KK)
256                 COUNT_FAC = COUNT_FAC + LIST_FACET_AT_DES(IJK2)%COUNT_FACETS
257              ENDDO
258              ENDDO
259              ENDDO
260     
261              IF(COUNT_FAC == 0) NO_NEIGHBORING_FACET_DES(IJK) = .TRUE.
262     
263           ENDDO
264           ENDDO
265           ENDDO
266     
267           RETURN
268           END SUBROUTINE BIN_FACETS_TO_GRID_DES
269     
270     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
271     !                                                                      !
272     !  Subroutine: ADD_FACET_FOR_DES                                       !
273     !  Author: Rahul Garg                                  Date: 24-Oct-13 !
274     !                                                                      !
275     !  Purpose: Add facets to DES grid cells.                              !
276     !                                                                      !
277     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
278           SUBROUTINE ADD_FACET_FOR_DES(I,J,K,IJK,N)
279     
280           USE param1, only: zero, one
281           use geometry, only: DO_K
282     
283           use desgrid, only: dg_dxinv
284           use desgrid, only: dg_dyinv
285           use desgrid, only: dg_dzinv
286     
287           Use discretelement, only: MAX_RADIUS
288     
289           Use stl
290           use error_manager
291     
292           IMPLICIT NONE
293     
294     ! DES grid index and facet index
295           INTEGER, INTENT(IN) :: I,J,K,IJK, N
296     
297     ! Center of DES grid cell and half size. Note that a buffer is added to
298     ! the half size to make the cell appear a little larger. This ensures
299     ! that paricles near the edge 'see' STLs that are nearby but do not
300     ! directly intersect the DES grid cell contain the particle center.
301           DOUBLE PRECISION :: CENTER(3), HALFSIZE(3)
302     ! Flag: STL intersects the DES grid cell
303           LOGICAL :: OVERLAP
304     ! DES grid cell dimensions
305           DOUBLE PRECISION :: lDX, lDY, lDZ
306     ! Buffer to ensure all particle-STL collisions are captured.
307           DOUBLE PRECISION :: BUFFER
308     ! Legacy variable - should be removed
309           INTEGER ::  CURRENT_COUNT
310     
311     
312           BUFFER = 1.1d0*MAX_RADIUS
313     
314           lDX = ONE/DG_DXINV
315           lDY = ONE/DG_DYINV
316           lDZ = ONE/DG_DZINV
317     
318           CENTER(1) = (dble(I-2)+HALF)*lDX
319           HALFSIZE(1) = HALF*lDX + BUFFER
320     
321           CENTER(2) = (dble(J-2)+HALF)*lDY
322           HALFSIZE(2) = HALF*lDY + BUFFER
323     
324           IF(DO_K)THEN
325              CENTER(3) = (dble(K-2)+HALF)*lDZ
326              HALFSIZE(3) = HALF*lDZ + BUFFER
327           ELSE
328              CENTER(3) = HALF*lDZ
329              HALFSIZE(3) = HALF*lDZ
330           ENDIF
331     
332           CALL TRI_BOX_OVERLAP(CENTER, HALFSIZE, VERTEX(:,:,N), OVERLAP)
333     
334           IF(OVERLAP) THEN
335     
336              CALL ADD_FACET(IJK, N)
337              CURRENT_COUNT = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
338     
339     ! These variables are not really needed. They are currently used
340     ! by the PIC model and should be moved to the new forms used by
341     ! the DEM wall collision routine.
342              IF(CURRENT_COUNT .LT. MAX_FACETS_PER_CELL_DES) THEN
343                 LIST_FACET_AT_DES(IJK)%COUNT_FACETS = CURRENT_COUNT+1
344                 LIST_FACET_AT_DES(IJK)%FACET_LIST(CURRENT_COUNT+1) = N
345              ELSE
346     
347              ENDIF
348           ENDIF
349     
350           RETURN
351           END SUBROUTINE ADD_FACET_FOR_DES
352     
353     
354     !......................................................................!
355     ! Subroutine TRI_BOX_OVERLAP                                           !
356     ! Author: J.Musser                                   Date: 10-22-2015  !
357     !                                                                      !
358     ! Purpose: Determine if a box (DES grid cell) intersects the triangle  !
359     !    (SLT). Note that the DES grid size is slightly increased to       !
360     !    capture STLs near the boarder of the cell. Otherwise, some        !
361     !    collisions could ve over looked.                                  !
362     !                                                                      !
363     ! Author: Tomas Akenine-Moller                   Accessed: 10-22-2015  !
364     ! REF: http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/    !
365     !         code/tribox2.txt                                             !
366     !......................................................................!
367           SUBROUTINE TRI_BOX_OVERLAP(pCENTER, pHALFSIZE, pVERTS, pOVERLAP)
368     
369           IMPLICIT NONE
370     
371           DOUBLE PRECISION, INTENT(IN) :: pCENTER(3), pHALFSIZE(3)
372           DOUBLE PRECISION, INTENT(IN) :: pVERTS(3,3)
373           LOGICAL, INTENT(OUT) :: pOVERLAP
374     
375           DOUBLE PRECISION :: v0(3), v1(3), v2(3)
376           DOUBLE PRECISION :: fex, fey, fez
377           DOUBLE PRECISION :: normal(3), e0(3), e1(3), e2(3)
378     
379           pOVERLAP = .FALSE.
380     
381           v0 = pVERTS(1,:) - pCENTER
382           v1 = pVERTS(2,:) - pCENTER
383           v2 = pVERTS(3,:) - pCENTER
384     
385           e0 = v1-v0
386           e1 = v2-v1
387           e2 = v0-v2
388     
389           fex = abs(e0(1))
390           fey = abs(e0(2))
391           fez = abs(e0(3))
392     
393           if(ATEST_X01(e0(3),e0(2),fez,fey)) return
394           if(ATEST_Y02(e0(3),e0(1),fez,fex)) return
395           if(ATEST_Z12(e0(2),e0(1),fey,fex)) return
396     
397           fex = abs(e1(1))
398           fey = abs(e1(2))
399           fez = abs(e1(3))
400     
401           if(ATEST_X01(e1(3),e1(2),fez,fey)) return
402           if(ATEST_Y02(e1(3),e1(1),fez,fex)) return
403           if(ATEST_Z0 (e1(2),e1(1),fey,fex)) return
404     
405           fex = abs(e2(1))
406           fey = abs(e2(2))
407           fez = abs(e2(3))
408     
409           if(ATEST_X2 (e2(3),e2(2),fez,fey)) return
410           if(ATEST_Y1 (e2(3),e2(1),fez,fex)) return
411           if(ATEST_Z12(e2(2),e2(1),fey,fex)) return
412     
413           if(findMin(v0(1),v1(1),v2(1)) > phalfsize(1) .OR. &
414              findMax(v0(1),v1(1),v2(1)) <-phalfsize(1)) return
415     
416           if(findMin(v0(2),v1(2),v2(2)) > phalfsize(2) .OR. &
417              findMax(v0(2),v1(2),v2(2)) <-phalfsize(2)) return
418     
419           if(findMin(v0(3),v1(3),v2(3)) > phalfsize(3) .OR. &
420              findMax(v0(3),v1(3),v2(3)) <-phalfsize(3)) return
421     
422     
423           normal(1) = e0(2)*e1(3)-e0(3)*e1(2)
424           normal(2) = e0(3)*e1(1)-e0(1)*e1(3)
425           normal(3) = e0(1)*e1(2)-e0(2)*e1(1)
426     
427           if(.NOT.planeBoxOverlap(normal,v0,phalfsize)) return
428     
429           pOVERLAP = .TRUE.
430     
431           RETURN
432     
433           CONTAINS
434     
435     !``````````````````````````````````````````````````````````````````````!
436     !                                                                      !
437     !                                                                      !
438     !                                                                      !
439     !``````````````````````````````````````````````````````````````````````!
440           LOGICAL FUNCTION planeBoxOverlap(norm, vert, maxbox)
441     
442           double precision :: norm(3), vert(3), maxbox(3)
443     
444           integer :: lc
445           double precision :: vmin(3), vmax(3), v
446     
447           do lc=1,3
448              v=vert(lc)
449              if(norm(lc) > 0.0d0) then
450                 vmin(lc) = -maxbox(lc) - v
451                 vmax(lc) =  maxbox(lc) - v
452              else
453                 vmin(lc) = maxbox(lc) - v
454                 vmax(lc) =-maxbox(lc) - v
455              endif
456           enddo
457     
458           if(dot_product(norm,vmin) > 0.0d0) RETURN
459           planeBoxOverlap=(dot_product(norm,vmax) >= 0.0d0)
460     
461           RETURN
462           END FUNCTION planeBoxOverlap
463     
464     !``````````````````````````````````````````````````````````````````````!
465     !                                                                      !
466     !                                                                      !
467     !                                                                      !
468     !``````````````````````````````````````````````````````````````````````!
469           DOUBLE PRECISION FUNCTION findMin(x0,x1,x2)
470     
471           double precision :: x0,x1,x2
472     
473           findMin = x0
474     
475           if(x1<findMin) findMin=x1
476           if(x2<findMin) findMin=x2
477     
478           RETURN
479           END FUNCTION findMin
480     
481     !``````````````````````````````````````````````````````````````````````!
482     !                                                                      !
483     !                                                                      !
484     !                                                                      !
485     !``````````````````````````````````````````````````````````````````````!
486           DOUBLE PRECISION FUNCTION findMax(x0,x1,x2)
487     
488           double precision :: x0,x1,x2
489     
490           findMax = x0
491     
492           if(x1>findMax) findMax=x1
493           if(x2>findMax) findMax=x2
494     
495           RETURN
496           END FUNCTION findMax
497     
498     !``````````````````````````````````````````````````````````````````````!
499     !                                                                      !
500     !                                                                      !
501     !                                                                      !
502     !``````````````````````````````````````````````````````````````````````!
503           LOGICAL FUNCTION ATEST_X01(a,b,fa,fb)
504     
505           double precision :: a, b, fa, fb
506           double precision :: lMin, lMax, p0, p2, rad
507     
508           p0 = a*v0(2) - b*v0(3)
509           p2 = a*v2(2) - b*v2(3)
510     
511           if(p0<p2) then; lMIN=p0; lMAX=p2
512           else; lMIN=p2; lMAX=p0; endif
513     
514           rad=fa*phalfsize(2) + fb*phalfsize(3)
515           ATEST_X01=(lmin>rad .OR. lmax<-rad)
516     
517           END FUNCTION ATEST_X01
518     
519     !``````````````````````````````````````````````````````````````````````!
520     !                                                                      !
521     !                                                                      !
522     !                                                                      !
523     !``````````````````````````````````````````````````````````````````````!
524           LOGICAL FUNCTION ATEST_X2(a,b,fa,fb)
525     
526           double precision :: a, b, fa, fb
527           double precision :: lMin, lMax, p0, p1, rad
528     
529           p0 = a*v0(2) - b*v0(3)
530           p1 = a*v1(2) - b*v1(3)
531     
532           if(p0<p1) then; lMIN=p0; lMAX=p1
533           else; lMIN=p1; lMAX=p0; endif
534     
535           rad=fa*phalfsize(2) + fb*phalfsize(3)
536           ATEST_X2=(lmin>rad .OR. lmax<-rad)
537     
538           END FUNCTION ATEST_X2
539     
540     !``````````````````````````````````````````````````````````````````````!
541     !                                                                      !
542     !                                                                      !
543     !                                                                      !
544     !``````````````````````````````````````````````````````````````````````!
545           LOGICAL FUNCTION ATEST_Y02(a,b,fa,fb)
546     
547           double precision :: a, b, fa, fb
548           double precision :: lMin, lMax, p0, p2, rad
549     
550           p0 = -a*v0(1) + b*v0(3)
551           p2 = -a*v2(1) + b*v2(3)
552     
553           if(p0<p2) then; lMIN=p0; lMAX=p2
554           else; lMIN=p2; lMAX=p0; endif
555     
556           rad=fa*phalfsize(1) + fb*phalfsize(3)
557           ATEST_Y02=(lmin>rad .OR. lmax<-rad)
558     
559           END FUNCTION ATEST_Y02
560     
561     !``````````````````````````````````````````````````````````````````````!
562     !                                                                      !
563     !                                                                      !
564     !                                                                      !
565     !``````````````````````````````````````````````````````````````````````!
566           LOGICAL FUNCTION ATEST_Y1(a,b,fa,fb)
567     
568           double precision :: a, b, fa, fb
569           double precision :: lMin, lMax, p0, p1, rad
570     
571           p0 = -a*v0(1) + b*v0(3)
572           p1 = -a*v1(1) + b*v1(3)
573     
574           if(p0<p1) then; lMIN=p0; lMAX=p1
575           else; lMIN=p1; lMAX=p0; endif
576     
577           rad=fa*phalfsize(1) + fb*phalfsize(3)
578           ATEST_Y1=(lmin>rad .OR. lmax<-rad)
579     
580           END FUNCTION ATEST_Y1
581     
582     !``````````````````````````````````````````````````````````````````````!
583     !                                                                      !
584     !                                                                      !
585     !                                                                      !
586     !``````````````````````````````````````````````````````````````````````!
587           LOGICAL FUNCTION ATEST_Z12(a,b,fa,fb)
588     
589           double precision :: a, b, fa, fb
590           double precision :: lMin, lMax, p1, p2, rad
591     
592           p1 = a*v1(1) - b*v1(2)
593           p2 = a*v2(1) - b*v2(2)
594     
595           if(p2<p1) then; lMIN=p2; lMAX=p1
596           else; lMIN=p1; lMAX=p2; endif
597     
598           rad=fa*phalfsize(1) + fb*phalfsize(2)
599           ATEST_Z12=(lmin>rad .OR. lmax<-rad)
600     
601           END FUNCTION ATEST_Z12
602     
603     !``````````````````````````````````````````````````````````````````````!
604     !                                                                      !
605     !                                                                      !
606     !                                                                      !
607     !``````````````````````````````````````````````````````````````````````!
608           LOGICAL FUNCTION ATEST_Z0(a,b,fa,fb)
609     
610           double precision :: a, b, fa, fb
611           double precision :: lMin, lMax, p0, p1, rad
612     
613           p0 = a*v0(1) - b*v0(2)
614           p1 = a*v1(1) - b*v1(2)
615     
616           if(p0<p1) then; lMIN=p0; lMAX=p1
617           else; lMIN=p1; lMAX=p0; endif
618     
619           rad=fa*phalfsize(1) + fb*phalfsize(2)
620           ATEST_Z0=(lmin>rad .OR. lmax<-rad)
621     
622           END FUNCTION ATEST_Z0
623     
624           END SUBROUTINE TRI_BOX_OVERLAP
625     
626     
627     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
628     !                                                                      !
629     !  Subroutine: DEBUG_WRITE_GRID_FACEINFO                               !
630     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
631     !                                                                      !
632     !  Purpose:                                                            !
633     !                                                                      !
634     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
635     
636           SUBROUTINE DEBUG_WRITE_GRID_FACEINFO
637           USE param
638           USE param1
639           USE run
640           USE geometry
641           USE indices
642           USE compar
643           USE stl
644           USE functions
645     
646           IMPLICIT NONE
647     
648           INTEGER :: CELL_ID, I, J, K, COUNT, COUNT_FACETS
649     
650           CHARACTER(LEN=100) :: FN
651     
652           IF(NODESI*NODESJ*NODESK == 1) THEN
653              WRITE(FN,'("FACETS_DG_GRID.DAT")')
654           ELSE
655              WRITE(FN,'("FACETS_DG_GRID_",I5.5,".DAT")') myPE
656           ENDIF
657     
658           OPEN(1001, file = TRIM(FN), form ='formatted',CONVERT='BIG_ENDIAN')
659     
660           DO K=DG_KSTART2, DG_KEND2
661           DO J=DG_JSTART2, DG_JEND2
662           DO I=DG_ISTART2, DG_IEND2
663     
664              CELL_ID = DG_FUNIJK(I,J,K)
665              COUNT_FACETS =  LIST_FACET_AT_DES(CELL_ID)%COUNT_FACETS
666              IF(COUNT_FACETS.eq.0) cycle
667     
668              WRITE(1001,2000) CELL_ID, I, J, K,  COUNT_FACETS
669     
670      2000 FORMAT(50('*'),/2X,'CELL IJK, I, J, K =        = ', i20, 2x,     &
671              4(2x,i10),/2X,'TOTAL FACETS',18(' '),'= ',3(2x, i10))
672     
673              DO COUNT = 1, COUNT_FACETS
674                 WRITE(1001, '(2x, i20)')                                   &
675                    LIST_FACET_AT_DES(CELL_ID)%FACET_LIST(COUNT)
676              ENDDO
677     
678           ENDDO
679           ENDDO
680           ENDDO
681     
682           CLOSE(1001, STATUS = "keep")
683     
684           RETURN
685           END SUBROUTINE DEBUG_WRITE_GRID_FACEINFO
686     
687     
688     
689     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
690     !                                                                      !
691     !  Subroutine: DEBUG_write_all_readin_facets                           !
692     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
693     !                                                                      !
694     !  Purpose:                                                            !
695     !                                                                      !
696     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
697           Subroutine  DEBUG_write_all_readin_facets
698     
699           USE stl
700           USE compar
701     
702           IMPLICIT NONE
703     
704           INTEGER :: N
705     !      CHARACTER(LEN=100) :: FN
706     
707     
708           WRITE(ERR_MSG, 2000) N_FACETS, N_FACETS_DES
709           CALL FLUSH_ERR_MSG(HEADER=.FALSE.,FOOTER=.FALSE.)
710     
711      2000 FORMAT('Saving STL geometry on DES grid: FACETS_TO_DG.stl',      &
712              11x,'Facet Count:',I10,/2x,'DES Grid Facet Count:',I10)
713     
714           OPEN(UNIT=444, FILE='FACETS_TO_DG.stl',CONVERT='BIG_ENDIAN')
715           write(444,*)'solid vcg'
716     
717           DO N = 1, N_FACETS_DES
718              write(444,*) '   facet normal ', NORM_FACE(1:3,N)
719              write(444,*) '      outer loop'
720              write(444,*) '         vertex ', VERTEX(1,1:3,N)
721              write(444,*) '         vertex ', VERTEX(2,1:3,N)
722              write(444,*) '         vertex ', VERTEX(3,1:3,N)
723              write(444,*) '      endloop'
724              write(444,*) '   endfacet'
725           ENDDO
726           write(444,*)'endsolid vcg'
727           close(444)
728     
729     
730     
731           END Subroutine DEBUG_write_all_readin_facets
732     
733     
734     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
735     !                                                                      !
736     !  Subroutine: DEBUG_write_stl_from_grid_facet                         !
737     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
738     !                                                                      !
739     !  Purpose:                                                            !
740     !                                                                      !
741     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
742           SUBROUTINE DEBUG_WRITE_STL_FROM_GRID_FACET(WRITE_FACETS_EACH_CELL)
743     
744           use run
745           USE stl
746     
747           USE functions
748           USE geometry
749           USE indices
750           USE compar
751           USE discretelement, only: STL_FACET_TYPE, FACET_TYPE_NORMAL, &
752                FACET_TYPE_MI, FACET_TYPE_PO, &
753                COUNT_FACET_TYPE_PO
754     
755           IMPLICIT NONE
756     
757           LOGICAL, INTENT(IN),optional  :: WRITE_FACETS_EACH_CELL
758           INTEGER ::  CELL_ID, N, I, J, K, COUNT, COUNT_FACETS, W_UNIT
759           CHARACTER(LEN=200) :: FN, FN_PO
760           LOGICAL :: write_each_cell
761           LOGICAL, DIMENSION(:), allocatable :: FACET_WRITTEN
762     
763           ALLOCATE (FACET_WRITTEN(DIM_STL))
764     
765           write_each_cell = .false.
766           if(present(WRITE_FACETS_EACH_CELL)) then
767              write_each_cell = WRITE_FACETS_EACH_CELL
768           endif
769     
770           FACET_WRITTEN = .false.
771     
772           IF(NODESI*NODESJ*NODESK == 1) THEN
773              WRITE(FN,'("GEOM_DG.stl")')
774              WRITE(FN_PO,'("GEOM_DG_PO.stl")')
775           ELSE
776              WRITE(FN,'("GEOM_DG_",I5.5,".stl")') MYPE
777              WRITE(FN_PO,'("GEOM_DG_PO_",I5.5,".stl")')  MYPE
778           ENDIF
779     
780           IF(COUNT_FACET_TYPE_PO.GE.1) THEN
781              OPEN(UNIT=443, FILE=TRIM(FN_PO),CONVERT='BIG_ENDIAN')
782              write(443,*)'solid vcg'
783           endif
784     
785           OPEN(UNIT=444, FILE=trim(FN),CONVERT='BIG_ENDIAN')
786           write(444,*)'solid vcg'
787     
788           DO K=DG_KSTART2, DG_KEND2
789           DO J=DG_JSTART2, DG_JEND2
790           DO I=DG_ISTART2, DG_IEND2
791     
792              CELL_ID = DG_FUNIJK(I,J,K)
793              COUNT_FACETS =  LIST_FACET_AT_DES(CELL_ID)%COUNT_FACETS
794     
795              IF(COUNT_FACETS.EQ.0) CYCLE
796     
797              IF(WRITE_EACH_CELL) THEN
798     
799                 WRITE(FN,'("GEOM_DG",3("_",I3.3),"_",I8.8,".stl")')        &
800                    I, J, K, CELL_ID
801     
802     
803                 WRITE(FN_PO,'("GEOM_DG_PO",3("_",I3.3),"_",I8.8,".stl")') &
804                    I,J,K,CELL_ID
805     
806                 OPEN(UNIT=446, FILE=FN,CONVERT='BIG_ENDIAN')
807                 WRITE(446,*)'solid vcg'
808     
809                 OPEN(UNIT=445, FILE=FN_PO,CONVERT='BIG_ENDIAN')
810                 WRITE(445,*)'solid vcg'
811              ENDIF
812     
813              DO COUNT = 1, COUNT_FACETS
814                 N = LIST_FACET_AT_DES(CELL_ID)%FACET_LIST(COUNT)
815     
816                 IF(WRITE_EACH_CELL) THEN
817                    IF(STL_FACET_TYPE(N).EQ.FACET_TYPE_NORMAL) THEN
818                       W_UNIT  = 446
819                    ELSE
820                       W_UNIT  = 445
821                    ENDIF
822     
823                    write(w_unit,*) '   facet normal ', NORM_FACE(:,N)
824                    write(w_unit,*) '      outer loop'
825                    write(w_unit,*) '         vertex ', VERTEX(1,1:3,N)
826                    write(w_unit,*) '         vertex ', VERTEX(2,1:3,N)
827                    write(w_unit,*) '         vertex ', VERTEX(3,1:3,N)
828                    write(w_unit,*) '      endloop'
829                    write(w_unit,*) '   endfacet'
830     
831                 ENDIF
832     
833                 IF (FACET_WRITTEN(N)) CYCLE
834     
835                 IF(STL_FACET_TYPE(N).EQ.FACET_TYPE_NORMAL) THEN
836                    W_UNIT  = 444
837                 ELSE
838                    W_UNIT  = 443
839                 ENDIF
840     
841                 write(w_unit,*) '   facet normal ', NORM_FACE(:,N)
842                 write(w_unit,*) '      outer loop'
843                 write(w_unit,*) '         vertex ', VERTEX(1,1:3,N)
844                 write(w_unit,*) '         vertex ', VERTEX(2,1:3,N)
845                 write(w_unit,*) '         vertex ', VERTEX(3,1:3,N)
846                 write(w_unit,*) '      endloop'
847                 write(w_unit,*) '   endfacet'
848                 facet_written(N) = .true.
849              ENDDO
850     
851              if(write_each_cell) then
852                 write(445,*)'endsolid vcg'
853                 close(445)
854     
855                 write(446,*)'endsolid vcg'
856                 close(446)
857              ENDIF
858     
859              ENDDO
860           ENDDO
861           ENDDO
862           write(444,*)'endsolid vcg'
863           close(444)
864     
865           if(count_facet_type_po.gt.0) then
866              write(443,*)'endsolid vcg'
867              close(443)
868           ENDIF
869     
870           DEALLOCATE (FACET_WRITTEN)
871     
872     !      IF(MyPE == PE_IO) THEN
873     !         WRITE(*,*) ' The file geometry_from_grid_facets.stl was sucessfully written.'
874     !         WRITE(*,*) ' This is the based on facets stored grid wise for DES modules'
875     !         WRITE(*,*) ' and is provided for convenience and debugging (it is not used).'
876     !      ENDIF
877           END Subroutine  DEBUG_write_stl_from_grid_facet
878     
879     
880     
881     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
882     !                                                                      !
883     !  Subroutine: ADD_FACET                                               !
884     !                                                                      !
885     !                                                                      !
886     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
887           SUBROUTINE ADD_FACET(CELL_ID, FACET_ID)
888     
889           Use discretelement
890           USE stl
891     
892           implicit none
893     
894           INTEGER, INTENT(IN) :: cell_id, facet_id
895     
896           INTEGER, DIMENSION(:), ALLOCATABLE :: int_tmp
897           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: real_tmp
898     
899           INTEGER :: lSIZE2, ii
900           DOUBLE PRECISION :: smallest_extent, min_temp, max_temp
901     
902           IF(STL_FACET_TYPE(facet_id) /= FACET_TYPE_NORMAL) RETURN
903     
904           DO II = 1, CELLNEIGHBOR_FACET_NUM(CELL_ID)
905              IF(FACET_ID .EQ. CELLNEIGHBOR_FACET(CELL_ID)%P(II)) RETURN
906           ENDDO
907     
908           CELLNEIGHBOR_FACET_NUM(CELL_ID) = &
909              CELLNEIGHBOR_FACET_NUM(CELL_ID) + 1
910     
911           NO_NEIGHBORING_FACET_DES(CELL_ID)  = .FALSE.
912     
913           IF(cellneighbor_facet_num(cell_id) > &
914              cellneighbor_facet_max(cell_id)) THEN
915     
916              cellneighbor_facet_max(cell_id) = &
917              2*cellneighbor_facet_max(cell_id)
918     
919              lSIZE2 = size(cellneighbor_facet(cell_id)%p)
920              allocate(int_tmp(cellneighbor_facet_max(cell_id)))
921              int_tmp(1:lSIZE2) = cellneighbor_facet(cell_id)%p(1:lSIZE2)
922              call move_alloc(int_tmp,cellneighbor_facet(cell_id)%p)
923     
924              lSIZE2 = size(cellneighbor_facet(cell_id)%extentdir)
925              allocate(int_tmp(cellneighbor_facet_max(cell_id)))
926              int_tmp(1:lSIZE2) = &
927                 cellneighbor_facet(cell_id)%extentdir(1:lSIZE2)
928              call move_alloc(int_tmp,cellneighbor_facet(cell_id)%extentdir)
929     
930              lSIZE2 = size(cellneighbor_facet(cell_id)%extentmin)
931              allocate(real_tmp(cellneighbor_facet_max(cell_id)))
932              real_tmp(1:lSIZE2) = &
933                 cellneighbor_facet(cell_id)%extentmin(1:lSIZE2)
934              call move_alloc(real_tmp,cellneighbor_facet(cell_id)%extentmin)
935     
936              lSIZE2 = size(cellneighbor_facet(cell_id)%extentmax)
937              allocate(real_tmp(cellneighbor_facet_max(cell_id)))
938              real_tmp(1:lSIZE2) = &
939                 cellneighbor_facet(cell_id)%extentmax(1:lSIZE2)
940              call move_alloc(real_tmp,cellneighbor_facet(cell_id)%extentmax)
941     
942           ENDIF
943     
944           CELLNEIGHBOR_FACET(CELL_ID)%&
945              P(CELLNEIGHBOR_FACET_NUM(CELL_ID)) = FACET_ID
946           SMALLEST_EXTENT = HUGE(0.0)
947     
948           DO II=1,3
949              MIN_TEMP = MINVAL(VERTEX(:,II,FACET_ID))
950              MAX_TEMP = MAXVAL(VERTEX(:,II,FACET_ID))
951     
952              IF(ABS(MAX_TEMP - MIN_TEMP) < SMALLEST_EXTENT ) THEN
953                  CELLNEIGHBOR_FACET(CELL_ID)%&
954                     EXTENTDIR(CELLNEIGHBOR_FACET_NUM(CELL_ID)) = II
955                  CELLNEIGHBOR_FACET(CELL_ID)%&
956                     EXTENTMIN(CELLNEIGHBOR_FACET_NUM(CELL_ID)) = MIN_TEMP
957                  CELLNEIGHBOR_FACET(CELL_ID)%&
958                     EXTENTMAX(CELLNEIGHBOR_FACET_NUM(CELL_ID)) = MAX_TEMP
959                  SMALLEST_EXTENT = ABS(MAX_TEMP - MIN_TEMP)
960              ENDIF
961           ENDDO
962     
963           END SUBROUTINE ADD_FACET
964     
965     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
966     !                                                                      !
967     !  SUBROUTINE: CHECK_IF_PARTICLE_OVERLAPS_STL                          !
968     !  Authors: Rahul Garg                               Date: 21-Mar-2014 !
969     !                                                                      !
970     !  Purpose: This subroutine is special written to check if a particle  !
971     !          overlaps any of the STL faces. The routine exits on         !
972     !          detecting an overlap. It is called after initial            !
973     !          generation of lattice configuration to remove out of domain !
974     !          particles                                                   !
975     !                                                                      !
976     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
977           SUBROUTINE CHECK_IF_PARTICLE_OVERLAPS_STL(POS, fI, fJ, fK, REMOVE)
978     
979           USE run
980           USE param1
981           USE discretelement, only: dimn, xe, yn, zt
982           USE geometry
983           USE constant
984           USE cutcell
985           USE indices
986           USE stl
987           USE compar
988           USE des_stl_functions
989           use desgrid
990           USE functions
991     
992           Implicit none
993     
994           DOUBLE PRECISION, INTENT(IN) :: POS(DIMN)
995           INTEGER, INTENT(IN) :: fI, fJ, fK
996           LOGICAL, INTENT(OUT) :: REMOVE
997     
998     ! Integers mapping the fluid cell corners to DES Grid indices.
999           INTEGER :: I1, I2, J1, J2, K1, K2
1000     
1001           INTEGER I, J, K, IJK, NF, LC
1002     
1003           DOUBLE PRECISION :: LINE_T
1004           DOUBLE PRECISION :: RADSQ, DIST(3)
1005     
1006           REMOVE = .TRUE.
1007     
1008           I1 = IofPOS(XE(fI-1))
1009           I2 = IofPOS(XE( fI ))
1010     
1011           J1 = JofPOS(YN(fJ-1))
1012           J2 = JofPOS(YN( fJ ))
1013     
1014           K1 = KofPOS(ZT(fK-1))
1015           K2 = KofPOS(ZT( fK ))
1016     
1017           RADSQ = (1.05d0*MAX_RADIUS)**2
1018     
1019           DO K = K1, K2
1020           DO J = J1, J2
1021           DO I = I1, I2
1022     
1023              IF(.NOT.DG_is_ON_myPE_plus1layers(I,J,K)) CYCLE
1024     
1025              IJK = DG_FUNIJK(I,J,K)
1026     
1027     ! The point is on the non-fluid side of the plane if t>0
1028              DO LC = 1, LIST_FACET_AT_DES(IJK)%COUNT_FACETS
1029                 NF = LIST_FACET_AT_DES(IJK)%FACET_LIST(LC)
1030     
1031                 CALL INTERSECTLNPLANE(POS, NORM_FACE(:,NF), &
1032                    VERTEX(1,:,NF), NORM_FACE(:,NF), LINE_T)
1033     
1034     ! Orthogonal projection puts the point outside of the domain or less than
1035     ! one particle radius to the facet.
1036                 DIST = LINE_T*NORM_FACE(:,NF)
1037                 IF(LINE_T > ZERO .OR. dot_product(DIST,DIST)<=RADSQ) RETURN
1038     
1039              ENDDO
1040           ENDDO
1041           ENDDO
1042           ENDDO
1043     
1044           REMOVE = .FALSE.
1045     
1046           RETURN
1047           END SUBROUTINE CHECK_IF_PARTICLE_OVERLAPS_STL
1048     
1049     
1050     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1051     !                                                                      !
1052     !                                                                      !
1053     !                                                                      !
1054     !                                                                      !
1055     !                                                                      !
1056     !                                                                      !
1057     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1058           SUBROUTINE write_this_facet_and_part(FID, PID)
1059           USE run
1060           USE param1
1061           USE discretelement
1062           USE geometry
1063           USE compar
1064           USE constant
1065           USE cutcell
1066           USE funits
1067           USE indices
1068           USE physprop
1069           USE parallel
1070           USE stl
1071           USE des_stl_functions
1072           Implicit none
1073           !facet id and particle id
1074           Integer, intent(in) :: fid, pid
1075           Integer :: stl_unit, vtp_unit , k
1076           CHARACTER(LEN=100) :: stl_fname, vtp_fname
1077           real :: temp_array(3)
1078     
1079           stl_unit = 1001
1080           vtp_unit = 1002
1081     
1082           WRITE(vtp_fname,'(A,"_OFFENDING_PARTICLE",".vtp")') TRIM(RUN_NAME)
1083           WRITE(stl_fname,'(A,"_STL_FACE",".stl")') TRIM(RUN_NAME)
1084     
1085           open(vtp_unit, file = vtp_fname, form='formatted',convert='big_endian')
1086           open(stl_unit, file = stl_fname, form='formatted',convert='big_endian')
1087     
1088           write(vtp_unit,"(a)") '<?xml version="1.0"?>'
1089           write(vtp_unit,"(a,es24.16,a)") '<!-- time =',s_time,'s -->'
1090           write(vtp_unit,"(a,a)") '<VTKFile type="PolyData"',&
1091                ' version="0.1" byte_order="LittleEndian">'
1092           write(vtp_unit,"(3x,a)") '<PolyData>'
1093           write(vtp_unit,"(6x,a,i10.10,a,a)")&
1094                '<Piece NumberOfPoints="',1,'" NumberOfVerts="0" ',&
1095                'NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">'
1096           write(vtp_unit,"(9x,a)")&
1097                '<PointData Scalars="Diameter" Vectors="Velocity">'
1098           write(vtp_unit,"(12x,a)")&
1099                '<DataArray type="Float32" Name="Diameter" format="ascii">'
1100           write (vtp_unit,"(15x,es13.6)") (2*des_radius(pid))
1101           write(vtp_unit,"(12x,a)") '</DataArray>'
1102     
1103           temp_array = zero
1104           temp_array(:) = des_vel_new(:,pid)
1105           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
1106                'Name="Velocity" NumberOfComponents="3" format="ascii">'
1107           write (vtp_unit,"(15x,3(es13.6,3x))")&
1108                ((temp_array(k)),k=1,3)
1109           write(vtp_unit,"(12x,a,/9x,a)") '</DataArray>','</PointData>'
1110           ! skip cell data
1111           write(vtp_unit,"(9x,a)") '<CellData></CellData>'
1112     
1113           temp_array = zero
1114           temp_array(1:dimn) = des_pos_new(1:dimn, pid)
1115           write(vtp_unit,"(9x,a)") '<Points>'
1116           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
1117                'Name="Position" NumberOfComponents="3" format="ascii">'
1118           write (vtp_unit,"(15x,3(es13.6,3x))")&
1119                ((temp_array(k)),k=1,3)
1120           write(vtp_unit,"(12x,a,/9x,a)")'</DataArray>','</Points>'
1121           ! Write tags for data not included (vtp format style)
1122           write(vtp_unit,"(9x,a,/9x,a,/9x,a,/9x,a)")'<Verts></Verts>',&
1123                '<Lines></Lines>','<Strips></Strips>','<Polys></Polys>'
1124           write(vtp_unit,"(6x,a,/3x,a,/a)")&
1125                '</Piece>','</PolyData>','</VTKFile>'
1126     
1127           !Now write the facet info
1128     
1129           write(stl_unit,*)'solid vcg'
1130     
1131           write(stl_unit,*) '   facet normal ', NORM_FACE(1:3,FID)
1132           write(stl_unit,*) '      outer loop'
1133           write(stl_unit,*) '         vertex ', VERTEX(1,1:3,FID)
1134           write(stl_unit,*) '         vertex ', VERTEX(2,1:3,FID)
1135           write(stl_unit,*) '         vertex ', VERTEX(3,1:3,FID)
1136           write(stl_unit,*) '      endloop'
1137           write(stl_unit,*) '   endfacet'
1138     
1139           write(stl_unit,*)'endsolid vcg'
1140     
1141           close(vtp_unit, status = 'keep')
1142           close(stl_unit, status = 'keep')
1143     
1144           end SUBROUTINE write_this_facet_and_part
1145     
1146     
1147     
1148     
1149     
1150     
1151     
1152     
1153     !----------------------------------------------------------------------!
1154     !                                                                      !
1155     !                                                                      !
1156     !                                                                      !
1157     !----------------------------------------------------------------------!
1158           SUBROUTINE write_stls_this_dg(dg)
1159     
1160           Use usr
1161     
1162     ! STL Vertices
1163           use stl, only: VERTEX
1164     ! STL Facet normals
1165           use stl, only: NORM_FACE
1166     
1167           use discretelement, only: CELLNEIGHBOR_FACET
1168           use discretelement, only: CELLNEIGHBOR_FACET_NUM
1169     
1170           IMPLICIT NONE
1171     !-----------------------------------------------
1172           integer, intent(in) :: dg
1173           integer :: lc, nf
1174     
1175           logical :: EXISTS
1176           character(len=6) :: IDX
1177     
1178           write(idx,"(I6.6)") dg
1179           open(unit=555, file='dg_'//idx//'.stl', status='UNKNOWN')
1180     
1181           write(555,*) 'solid vcg'
1182           DO lc=1, cellneighbor_facet_num(dg)
1183     
1184              NF = cellneighbor_facet(dg)%p(lc)
1185     
1186              write(555,*) '   facet normal ', NORM_FACE(:,NF)
1187              write(555,*) '      outer loop'
1188              write(555,*) '         vertex ', VERTEX(1,1:3,NF)
1189              write(555,*) '         vertex ', VERTEX(2,1:3,NF)
1190              write(555,*) '         vertex ', VERTEX(3,1:3,NF)
1191              write(555,*) '      endloop'
1192              write(555,*) '   endfacet'
1193           enddo
1194           close(555)
1195     
1196           RETURN
1197           END SUBROUTINE write_stls_this_dg
1198     
1199     
1200     !----------------------------------------------------------------------!
1201     !                                                                      !
1202     !                                                                      !
1203     !                                                                      !
1204     !----------------------------------------------------------------------!
1205           SUBROUTINE write_this_stl(lc)
1206     
1207     
1208     ! STL Vertices
1209           use stl, only: VERTEX
1210     ! STL Facet normals
1211           use stl, only: NORM_FACE
1212           use compar, only: myPE
1213     
1214           IMPLICIT NONE
1215     !-----------------------------------------------
1216           integer, intent(in) :: lc
1217     
1218           logical :: EXISTS
1219           character(len=4) :: IDX
1220           character(len=4) :: IPE
1221     
1222     
1223           write(idx,"(I4.4)") LC
1224           write(ipe,"(I4.4)") myPE
1225           open(unit=555, file='geo_'//idx//'_'//IPE//'.stl',&
1226              status='UNKNOWN')
1227           write(555,*) 'solid vcg'
1228           write(555,*) '   facet normal ', NORM_FACE(:,LC)
1229           write(555,*) '      outer loop'
1230           write(555,*) '         vertex ', VERTEX(1,1:3,LC)
1231           write(555,*) '         vertex ', VERTEX(2,1:3,LC)
1232           write(555,*) '         vertex ', VERTEX(3,1:3,LC)
1233           write(555,*) '      endloop'
1234           write(555,*) '   endfacet'
1235           close(555)
1236     
1237     
1238           RETURN
1239           END SUBROUTINE write_this_stl
1240     
1241           END MODULE STL_PREPROC_DES
1242     
1243     
1244