File: /nfs/home/0/users/jenkins/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: DES_CONVERT_BOX_TO_FACETS
35           USE discretelement, only: STL_FACET_TYPE
36           USE discretelement, only: FACET_TYPE_NORMAL
37           USE discretelement, only: FACET_TYPE_MI
38           USE discretelement, only: FACET_TYPE_PO
39           USE discretelement, only: COUNT_FACET_TYPE_NORMAL
40           USE discretelement, only: COUNT_FACET_TYPE_MI
41           USE discretelement, only: COUNT_FACET_TYPE_PO
42     
43           USE desgrid
44     
45           USE param
46           USE stl
47     
48           implicit none
49     
50           INTEGER :: IJK, COUNT, NF
51     
52     ! Pre-procssing for the des in order to assign facets to grid cells.
53           WRITE(ERR_MSG,"('Pre-Processing geometry for DES.')")
54           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
55     
56           CALL ALLOCATE_DES_STL_ARRAYS
57     
58           N_FACETS_DES = 0
59     
60     ! Set N_facets_des to add any more facets needed by dem and not to.
61     ! contaminate the Eulerian-Eulerian CG stuff
62           IF(USE_STL) N_FACETS_DES = N_FACETS
63     
64     ! Set all the facets generated so far as normal facet types
65           STL_FACET_TYPE(1:N_FACETS) = FACET_TYPE_NORMAL
66     
67     
68           CALL BIN_FACETS_TO_GRID_DES
69     
70           DO IJK = 1, DG_IJKSIZE2
71              COUNT = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
72              IF(COUNT.eq.0) DEALLOCATE(LIST_FACET_AT_DES(IJK)%FACET_LIST)
73           ENDDO
74     
75           COUNT_FACET_TYPE_NORMAL = 0
76           COUNT_FACET_TYPE_PO     = 0
77           COUNT_FACET_TYPE_MI     = 0
78           DO NF = 1, N_FACETS_DES
79              IF(STL_FACET_TYPE(NF) == FACET_TYPE_NORMAL) &
80                 COUNT_FACET_TYPE_NORMAL = COUNT_FACET_TYPE_NORMAL + 1
81     
82              IF(STL_FACET_TYPE(NF).EQ.FACET_TYPE_MI) &
83                 COUNT_FACET_TYPE_MI = COUNT_FACET_TYPE_MI + 1
84     
85              IF(STL_FACET_TYPE(NF).EQ.FACET_TYPE_PO) &
86                 COUNT_FACET_TYPE_PO = COUNT_FACET_TYPE_PO + 1
87           ENDDO
88     
89     !      CALL DEBUG_WRITE_GRID_FACEINFO
90     !      CALL DEBUG_WRITE_STL_FROM_GRID_FACET(WRITE_FACETS_EACH_CELL=.FALSE.)
91     
92     !      CALL DEBUG_WRITE_ALL_READIN_FACETS
93     
94           RETURN
95           END SUBROUTINE DES_STL_PREPROCESSING
96     
97     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
98     !                                                                      !
99     !  Subroutine: ALLOCATE_DES_STL_ARRAYS                                 !
100     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
101     !                                                                      !
102     !  Purpose:                                                            !
103     !                                                                      !
104     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
105           SUBROUTINE ALLOCATE_DES_STL_ARRAYS
106     
107           USE desgrid
108           USE discretelement, only: STL_Facet_type
109           USE discretelement, only: CELLNEIGHBOR_FACET
110           USE discretelement, only: CELLNEIGHBOR_FACET_MAX
111           USE discretelement, only: CELLNEIGHBOR_FACET_NUM
112           USE discretelement, only: CELLNEIGHBOR_FACET_NUM
113           USE stl, only: LIST_FACET_AT_DES
114           USE stl, only: NO_NEIGHBORING_FACET_DES
115           USE param
116           USE stl, only: DIM_STL, MAX_FACETS_PER_CELL_DES
117     
118           IMPLICIT NONE
119     
120           integer :: ii,ijk
121     
122           ALLOCATE(LIST_FACET_AT_DES(DG_IJKSIZE2))
123     
124           Allocate(  CELLNEIGHBOR_FACET (DG_IJKSIZE2) )
125           Allocate(  CELLNEIGHBOR_FACET_MAX (DG_IJKSIZE2) )
126           Allocate(  CELLNEIGHBOR_FACET_NUM (DG_IJKSIZE2) )
127           DO II = 1, DG_IJKSIZE2
128              CELLNEIGHBOR_FACET_MAX(II) = 4
129              allocate(CELLNEIGHBOR_FACET(II)%P(CELLNEIGHBOR_FACET_MAX(II)))
130              allocate(CELLNEIGHBOR_FACET(II)%EXTENTDIR(CELLNEIGHBOR_FACET_MAX(II)))
131              allocate(CELLNEIGHBOR_FACET(II)%EXTENTMIN(CELLNEIGHBOR_FACET_MAX(II)))
132              allocate(CELLNEIGHBOR_FACET(II)%EXTENTMAX(CELLNEIGHBOR_FACET_MAX(II)))
133              CELLNEIGHBOR_FACET_NUM(II) = 0
134           ENDDO
135     
136           DO IJK = 1, DG_IJKSIZE2
137              LIST_FACET_AT_DES(IJK)%COUNT_FACETS = 0
138              ALLOCATE(LIST_FACET_AT_DES(IJK)%FACET_LIST(MAX_FACETS_PER_CELL_DES))
139           ENDDO
140     
141           ALLOCATE(NO_NEIGHBORING_FACET_DES(DG_IJKSIZE2))
142           NO_NEIGHBORING_FACET_DES = .false.
143     
144           Allocate(stl_facet_type(DIM_STL))
145     
146           END SUBROUTINE ALLOCATE_DES_STL_ARRAYS
147     
148     
149     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
150     !                                                                      !
151     !  Subroutine: bin_facets_to_grid_des                                  !
152     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
153     !                                                                      !
154     !  Purpose:                                                            !
155     !                                                                      !
156     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
157           Subroutine BIN_FACETS_TO_GRID_DES
158     
159           USE desgrid
160           USE discretelement, only: dimn, xe, yn, zt
161           USE functions
162           USE geometry
163           USE indices
164           USE mpi_utility
165           USE param1, only: zero
166           Use stl, only: tol_stl, n_facets_des
167           Use stl, only: vertex, NO_NEIGHBORING_FACET_DES, LIST_FACET_AT_DES
168     
169           IMPLICIT NONE
170     
171     ! DES Grid cell index.
172           INTEGER :: IJK, IJK2
173     ! Loop counters:
174           INTEGER :: I, I1, I2, II  ! X-axis
175           INTEGER :: J, J1, J2, JJ  ! Y-axis
176           INTEGER :: K, K1, K2, KK  ! Z-axis
177           INTEGER :: N              ! STLs
178     ! Generic accumulator
179           INTEGER :: COUNT_FAC
180     
181     ! Maximum and minimum extents of the indexed STL
182           DOUBLE PRECISION:: X1,Y1,Z1
183           DOUBLE PRECISION:: X2,Y2,Z2
184     
185     
186     
187           DO N = 1,N_FACETS_DES
188     
189              X1 = minval(VERTEX(1:3,1,N))
190              X2 = maxval(VERTEX(1:3,1,N))
191              Y1 = minval(VERTEX(1:3,2,N))
192              Y2 = maxval(VERTEX(1:3,2,N))
193              Z1 = minval(VERTEX(1:3,3,N))
194              Z2 = maxval(VERTEX(1:3,3,N))
195     
196              I1 = DG_IEND2
197              I2 = DG_ISTART2
198              IF(X2>=ZERO .AND. X1<=XLENGTH+TOL_STL) THEN
199                 I1 = iofpos(X1)
200                 I2 = iofpos(1/dg_dxinv+X2)
201              ENDIF
202     
203              J1 = DG_JEND2
204              J2 = DG_JSTART2
205              IF(Y2>=ZERO .AND. Y1<=YLENGTH+TOL_STL) THEN
206                 J1 = jofpos(Y1)
207                 J2 = jofpos(1/dg_dyinv+Y2)
208              ENDIF
209     
210              K1 = DG_KEND2
211              K2 = DG_KSTART2
212              IF(DO_K) THEN
213                 IF(Z2>=ZERO .AND. Z1<=ZLENGTH+TOL_STL) THEN
214                    K1 = kofpos(Z1)
215                    K2 = kofpos(1/dg_dzinv+Z2)
216                 ENDIF
217              ENDIF
218     
219     
220              DO K=K1,K2
221              DO J=J1,J2
222              DO I=I1,I2
223                 IF(dg_is_ON_myPE_plus1layers(I,J,K)) THEN
224                    IJK = DG_FUNIJK(I,J,K)
225                    CALL ADD_FACET_FOR_DES(I,J,K,IJK,N)
226                 ENDIF
227              ENDDO
228              ENDDO
229              ENDDO
230     
231           ENDDO
232     
233           DO K = DG_KSTART1, DG_KEND1
234           DO J = DG_JSTART1, DG_JEND1
235           DO I = DG_ISTART1, DG_IEND1
236     
237              IJK = DG_FUNIJK(I, J, K)
238     
239              I1 =  MAX( I - 1, DG_ISTART2)
240              I2 =  MIN( I + 1, DG_IEND2)
241     
242              J1 =  MAX( J - 1, DG_JSTART2)
243              J2 =  MIN (J + 1, DG_JEND2)
244     
245              K1 =  MAX( K - 1, DG_KSTART2)
246              K2 =  MIN (K + 1, DG_KEND2)
247     
248              NO_NEIGHBORING_FACET_DES(IJK)  = .FALSE.
249     
250              COUNT_FAC = 0
251              DO KK = K1, k2
252              DO JJ = J1, j2
253              DO II = I1, i2
254                 IJK2 = DG_FUNIJK(II, JJ, KK)
255                 COUNT_FAC = COUNT_FAC + LIST_FACET_AT_DES(IJK2)%COUNT_FACETS
256              ENDDO
257              ENDDO
258              ENDDO
259     
260              IF(COUNT_FAC == 0) NO_NEIGHBORING_FACET_DES(IJK) = .TRUE.
261     
262           ENDDO
263           ENDDO
264           ENDDO
265     
266           RETURN
267           END SUBROUTINE BIN_FACETS_TO_GRID_DES
268     
269     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
270     !                                                                      C
271     !  Module name: ADD_FACET_for_des                                      C
272     !  Purpose: Add facet to list in IJK scalar cell for the               C
273     !           discrete modules.                                          C
274     !                                                                      C
275     !  Author: Rahul Garg                              Date: 24-Oct-13     C
276     !  Reviewer:                                          Date:            C
277     !                                                                      C
278     !  Revision Number #                                  Date: ##-###-##  C
279     !  Author: #                                                           C
280     !  Purpose: #                                                          C
281     !                                                                      C
282     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
283           SUBROUTINE ADD_FACET_FOR_DES(I,J,K,IJK,N)
284     
285           USE param1, only: zero, one
286           use geometry, only: DO_K, ZLENGTH
287     
288           use run, only: RUN_NAME
289     
290           use compar, only: myPE
291           use compar, only: NODESi, NODESj, NODESk
292     
293           use desgrid, only: dg_xstart, dg_dxinv
294           use desgrid, only: dg_ystart, dg_dyinv
295           use desgrid, only: dg_zstart, dg_dzinv
296     
297           Use stl
298           use error_manager
299     
300           IMPLICIT NONE
301     
302           INTEGER, INTENT(IN) :: I,J,K,IJK, N
303           !LOCAL VARIABLES
304     
305           INTEGER ::  CURRENT_COUNT, COUNT
306     
307           double precision ::   box_origin(3), box_extents(3)
308           Logical :: sa_exist
309           integer :: sep_axis
310     
311           CHARACTER(LEN=100) :: FNAME
312           integer :: stl_unit, fid
313     
314           stl_unit = 1001
315     
316           BOX_ORIGIN(1) = DG_XSTART + (I-DG_ISTART1)/DG_DXINV
317           BOX_EXTENTS(1) = 1.0/DG_DXINV
318     
319           BOX_ORIGIN(2) = DG_YSTART + (J-DG_JSTART1)/DG_DYINV
320           BOX_EXTENTS(2) = 1.0/DG_DYINV
321     
322           IF(DO_K)THEN
323              BOX_ORIGIN(3) = DG_ZSTART + (K-DG_KSTART1)/DG_DZINV
324              BOX_EXTENTS(3) = 1.0/DG_DZINV
325           ELSE
326              BOX_ORIGIN(3) = 0.0D0
327              BOX_EXTENTS(3) = ZLENGTH
328           ENDIF
329     
330     ! Use the separating axis test to determine if the cell and facet cannot
331     ! intersect one another.
332           CALL TESTTRIANGLEAABB(VERTEX(:,:,N), NORM_FACE(:,N),             &
333              BOX_ORIGIN(:), BOX_EXTENTS(:), SA_EXIST, SEP_AXIS,I,J,K )
334     
335     ! A sparating axis exists so the cell and the triangle do not intersect.
336           IF(SA_EXIST) RETURN
337     
338           CURRENT_COUNT = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
339     
340           IF(CURRENT_COUNT .LT. MAX_FACETS_PER_CELL_DES) THEN
341     
342              LIST_FACET_AT_DES(IJK)%COUNT_FACETS = CURRENT_COUNT+1
343              LIST_FACET_AT_DES(IJK)%FACET_LIST(CURRENT_COUNT+1) = N
344     
345           ELSE
346              CALL INIT_ERR_MSG("des_stl_functions_mod::add_facets_for_des")
347     
348              WRITE(err_msg, 200) MAX_FACETS_PER_CELL_DES, IJK, &
349                 I, J, K, mype, .FALSE.
350              CALL flush_err_msg(footer = .false.)
351     
352     
353              write(err_msg, *) "current_list for this cell is"
354              CALL flush_err_msg(header = .false., footer = .false.)
355     
356     
357              IF(nodesI*nodesJ*nodesK.gt.1) then
358     
359                 WRITE(fname,'(A,"_TROUBLE_CELL",A, I5.5, 3(A,I5.5), ".stl")') &
360                 TRIM(RUN_NAME), '_pid', mype, '_I', I, '_J', J, '_K', K
361              else
362     
363                 WRITE(fname,'(A,"_TROUBLE_CELL", 3(A,I5.5), ".stl")') &
364                 TRIM(RUN_NAME), '_I', I, '_J', J, '_K', K
365              endif
366     
367              open(stl_unit, file = fname, form='formatted')
368              write(stl_unit,*)'solid vcg'
369     
370              DO COUNT  = 1, CURRENT_COUNT
371     
372                  !write(*, '(/,I10)') LIST_FACET_AT(IJK)%FACET_LIST(COUNT)
373                 FID = LIST_FACET_AT_DES(IJK)%FACET_LIST(COUNT)
374                 write(err_msg, '(I10)') FID
375     
376                 CALL flush_err_msg(header = .false., footer = .false.)
377     
378     
379                  write(stl_unit,*) '   facet normal ', NORM_FACE(1:3,FID)
380                  write(stl_unit,*) '      outer loop'
381                  write(stl_unit,*) '         vertex ', VERTEX(1,1:3,FID)
382                  write(stl_unit,*) '         vertex ', VERTEX(2,1:3,FID)
383                  write(stl_unit,*) '         vertex ', VERTEX(3,1:3,FID)
384                  write(stl_unit,*) '      endloop'
385                  write(stl_unit,*) '   endfacet'
386     
387               ENDDO
388     
389               write(stl_unit,*)'endsolid vcg'
390               close(stl_unit, status = 'keep')
391     
392     
393               write(err_msg, *) 'Stopping'
394               CALL flush_err_msg(header = .false., abort = .true.)
395           ENDIF
396     
397     
398     
399      200  FORMAT('ERROR MESSAGE FROM CUT_CELL_PREPROCESSING', /10x,        &
400              'INCREASE MAX_FACETS_PER_CELL_DES from the current value of', &
401              I3, /10x,'Happening for cell IJK, I, J, K = ', 4(2x, i5),/10X,&
402              'mype, Is on myPe? ', I6, L2, /10X,'see the file ',&
403              'TROUBLE_CELL for all the current facets in this cell')
404     
405           END SUBROUTINE ADD_FACET_FOR_DES
406     
407     
408     
409     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
410     !                                                                      !
411     !  Subroutine: DEBUG_WRITE_GRID_FACEINFO                               !
412     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
413     !                                                                      !
414     !  Purpose:                                                            !
415     !                                                                      !
416     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
417     
418           SUBROUTINE DEBUG_WRITE_GRID_FACEINFO
419           USE param
420           USE param1
421           USE run
422           USE geometry
423           USE indices
424           USE compar
425           USE stl
426           USE functions
427     
428           IMPLICIT NONE
429     
430           INTEGER :: CELL_ID, I, J, K, COUNT, COUNT_FACETS
431     
432           CHARACTER(LEN=100) :: FN
433     
434           IF(NODESI*NODESJ*NODESK == 1) THEN
435              WRITE(FN,'("FACETS_DG_GRID.DAT")')
436           ELSE
437              WRITE(FN,'("FACETS_DG_GRID_",I5.5,".DAT")') myPE
438           ENDIF
439     
440           OPEN(1001, file = TRIM(FN), form ='formatted')
441     
442           DO K=DG_KSTART2, DG_KEND2
443           DO J=DG_JSTART2, DG_JEND2
444           DO I=DG_ISTART2, DG_IEND2
445     
446              CELL_ID = DG_FUNIJK(I,J,K)
447              COUNT_FACETS =  LIST_FACET_AT_DES(CELL_ID)%COUNT_FACETS
448              IF(COUNT_FACETS.eq.0) cycle
449     
450              WRITE(1001,2000) CELL_ID, I, J, K,  COUNT_FACETS
451     
452      2000 FORMAT(50('*'),/2X,'CELL IJK, I, J, K =        = ', i20, 2x,     &
453              4(2x,i10),/2X,'TOTAL FACETS',18(' '),'= ',3(2x, i10))
454     
455              DO COUNT = 1, COUNT_FACETS
456                 WRITE(1001, '(2x, i20)')                                   &
457                    LIST_FACET_AT_DES(CELL_ID)%FACET_LIST(COUNT)
458              ENDDO
459     
460           ENDDO
461           ENDDO
462           ENDDO
463     
464           CLOSE(1001, STATUS = "keep")
465     
466           RETURN
467           END SUBROUTINE DEBUG_WRITE_GRID_FACEINFO
468     
469     
470     
471     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
472     !                                                                      !
473     !  Subroutine: DEBUG_write_all_readin_facets                           !
474     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
475     !                                                                      !
476     !  Purpose:                                                            !
477     !                                                                      !
478     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
479           Subroutine  DEBUG_write_all_readin_facets
480     
481           USE stl
482           USE compar
483     
484           IMPLICIT NONE
485     
486           INTEGER :: N
487     !      CHARACTER(LEN=100) :: FN
488     
489     
490           WRITE(ERR_MSG, 2000) N_FACETS, N_FACETS_DES
491           CALL FLUSH_ERR_MSG(HEADER=.FALSE.,FOOTER=.FALSE.)
492     
493      2000 FORMAT('Saving STL geometry on DES grid: FACETS_TO_DG.stl',      &
494              11x,'Facet Count:',I10,/2x,'DES Grid Facet Count:',I10)
495     
496           OPEN(UNIT=444, FILE='FACETS_TO_DG.stl')
497           write(444,*)'solid vcg'
498     
499           DO N = 1, N_FACETS_DES
500              write(444,*) '   facet normal ', NORM_FACE(1:3,N)
501              write(444,*) '      outer loop'
502              write(444,*) '         vertex ', VERTEX(1,1:3,N)
503              write(444,*) '         vertex ', VERTEX(2,1:3,N)
504              write(444,*) '         vertex ', VERTEX(3,1:3,N)
505              write(444,*) '      endloop'
506              write(444,*) '   endfacet'
507           ENDDO
508           write(444,*)'endsolid vcg'
509           close(444)
510     
511     
512     
513           END Subroutine DEBUG_write_all_readin_facets
514     
515     
516     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
517     !                                                                      !
518     !  Subroutine: DEBUG_write_stl_from_grid_facet                         !
519     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
520     !                                                                      !
521     !  Purpose:                                                            !
522     !                                                                      !
523     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
524           SUBROUTINE DEBUG_WRITE_STL_FROM_GRID_FACET(WRITE_FACETS_EACH_CELL)
525     
526           use run
527           USE stl
528     
529           USE functions
530           USE geometry
531           USE indices
532           USE compar
533           USE discretelement, only: STL_FACET_TYPE, FACET_TYPE_NORMAL, &
534                FACET_TYPE_MI, FACET_TYPE_PO, COUNT_FACET_TYPE_NORMAL, &
535                COUNT_FACET_TYPE_MI, COUNT_FACET_TYPE_PO
536     
537           IMPLICIT NONE
538     
539           LOGICAL, INTENT(IN),optional  :: WRITE_FACETS_EACH_CELL
540           INTEGER ::  CELL_ID, N, I, J, K, COUNT, COUNT_FACETS, W_UNIT
541           CHARACTER(LEN=200) :: FN, FN_PO
542           LOGICAL :: write_each_cell
543           LOGICAL, DIMENSION(:), allocatable :: FACET_WRITTEN
544     
545           ALLOCATE (FACET_WRITTEN(DIM_STL))
546     
547           write_each_cell = .false.
548           if(present(WRITE_FACETS_EACH_CELL)) then
549              write_each_cell = WRITE_FACETS_EACH_CELL
550           endif
551     
552           FACET_WRITTEN = .false.
553     
554           IF(NODESI*NODESJ*NODESK == 1) THEN
555              WRITE(FN,'("GEOM_DG.stl")')
556              WRITE(FN_PO,'("GEOM_DG_PO.stl")')
557           ELSE
558              WRITE(FN,'("GEOM_DG_",I5.5,".stl")') MYPE
559              WRITE(FN_PO,'("GEOM_DG_PO_",I5.5,".stl")')  MYPE
560           ENDIF
561     
562           IF(COUNT_FACET_TYPE_PO.GE.1) THEN
563              OPEN(UNIT=443, FILE=TRIM(FN_PO))
564              write(443,*)'solid vcg'
565           endif
566     
567           OPEN(UNIT=444, FILE=trim(FN))
568           write(444,*)'solid vcg'
569     
570           DO K=DG_KSTART2, DG_KEND2
571           DO J=DG_JSTART2, DG_JEND2
572           DO I=DG_ISTART2, DG_IEND2
573     
574              CELL_ID = DG_FUNIJK(I,J,K)
575              COUNT_FACETS =  LIST_FACET_AT_DES(CELL_ID)%COUNT_FACETS
576     
577              IF(COUNT_FACETS.EQ.0) CYCLE
578     
579              IF(WRITE_EACH_CELL) THEN
580     
581                 WRITE(FN,'("GEOM_DG",3("_",I3.3),"_",I8.8,".stl")')        &
582                    I, J, K, CELL_ID
583     
584     
585                 WRITE(FN_PO,'("GEOM_DG_PO",3("_",I3.3),"_",I8.8,".stl")') &
586                    I,J,K,CELL_ID
587     
588                 OPEN(UNIT=446, FILE=FN)
589                 WRITE(446,*)'solid vcg'
590     
591                 OPEN(UNIT=445, FILE=FN_PO)
592                 WRITE(445,*)'solid vcg'
593              ENDIF
594     
595              DO COUNT = 1, COUNT_FACETS
596                 N = LIST_FACET_AT_DES(CELL_ID)%FACET_LIST(COUNT)
597     
598                 IF(WRITE_EACH_CELL) THEN
599                    IF(STL_FACET_TYPE(N).EQ.FACET_TYPE_NORMAL) THEN
600                       W_UNIT  = 446
601                    ELSE
602                       W_UNIT  = 445
603                    ENDIF
604     
605                    write(w_unit,*) '   facet normal ', NORM_FACE(:,N)
606                    write(w_unit,*) '      outer loop'
607                    write(w_unit,*) '         vertex ', VERTEX(1,1:3,N)
608                    write(w_unit,*) '         vertex ', VERTEX(2,1:3,N)
609                    write(w_unit,*) '         vertex ', VERTEX(3,1:3,N)
610                    write(w_unit,*) '      endloop'
611                    write(w_unit,*) '   endfacet'
612     
613                 ENDIF
614     
615                 IF (FACET_WRITTEN(N)) CYCLE
616     
617                 IF(STL_FACET_TYPE(N).EQ.FACET_TYPE_NORMAL) THEN
618                    W_UNIT  = 444
619                 ELSE
620                    W_UNIT  = 443
621                 ENDIF
622     
623                 write(w_unit,*) '   facet normal ', NORM_FACE(:,N)
624                 write(w_unit,*) '      outer loop'
625                 write(w_unit,*) '         vertex ', VERTEX(1,1:3,N)
626                 write(w_unit,*) '         vertex ', VERTEX(2,1:3,N)
627                 write(w_unit,*) '         vertex ', VERTEX(3,1:3,N)
628                 write(w_unit,*) '      endloop'
629                 write(w_unit,*) '   endfacet'
630                 facet_written(N) = .true.
631              ENDDO
632     
633              if(write_each_cell) then
634                 write(445,*)'endsolid vcg'
635                 close(445)
636     
637                 write(446,*)'endsolid vcg'
638                 close(446)
639              ENDIF
640     
641              ENDDO
642           ENDDO
643           ENDDO
644           write(444,*)'endsolid vcg'
645           close(444)
646     
647           if(count_facet_type_po.gt.0) then
648              write(443,*)'endsolid vcg'
649              close(443)
650           ENDIF
651     
652           DEALLOCATE (FACET_WRITTEN)
653     
654     !      IF(MyPE == PE_IO) THEN
655     !         WRITE(*,*) ' The file geometry_from_grid_facets.stl was sucessfully written.'
656     !         WRITE(*,*) ' This is the based on facets stored grid wise for DES modules'
657     !         WRITE(*,*) ' and is provided for convenience and debugging (it is not used).'
658     !      ENDIF
659           END Subroutine  DEBUG_write_stl_from_grid_facet
660     
661     
662           END MODULE STL_PREPROC_DES
663     
664     
665