File: N:\mfix\model\des\stl_preproc_des_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module name: stl_functions_des                                      !
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           IMPLICIT NONE
13     
14     ! Use this module only to define functions and subroutines.
15           CONTAINS
16     
17     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
18     !                                                                      !
19     !  Subroutine: DES_STL_PREPROCESSING                                   !
20     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
21     !                                                                      !
22     !  Purpose:                                                            !
23     !                                                                      !
24     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
25           SUBROUTINE DES_STL_PREPROCESSING
26     
27     ! Flag to for STL defined geometry
28           use cutcell, only: use_stl
29     ! Number of facets from STL files, (plus DES generated)
30           use stl, only: N_FACETS, N_FACETS_DES
31     ! Start/End position of different STLs
32           use stl, only: STL_START, STL_END
33     ! All STLS
34           use stl, only: ALL_STL
35     ! STLs read from geometry files
36           use stl, only: BASE_STL
37     ! STLs for user specified walls (NSW, PSW, FSW)
38           use stl, only: BCWALLS_STL
39     ! STLs for impermeable surfaces
40           use stl, only: IMPRMBL_STL
41     ! STLs for default walls
42           use stl, only: DEFAULT_STL
43     
44           use stl_dbg_des
45           use error_manager
46     
47           IMPLICIT NONE
48     
49     ! Pre-procssing for the des in order to assign facets to grid cells.
50           WRITE(ERR_MSG,"('Pre-Processing geometry for DES.')")
51           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
52     
53     ! Process the STL files
54           N_FACETS_DES = merge(N_FACETS, 0, USE_STL)
55     ! Store the Start/End of the base STLs from geometry files
56           STL_START(BASE_STL)=1;   STL_END(BASE_STL)=N_FACETS_DES
57     
58     ! Process stair-step geometries
59           CALL CONVERT_BC_WALLS_TO_STL
60     ! Process stair-step geometries
61           CALL CONVERT_IMPERMEABLE_IS_TO_STL
62     ! Process default walls
63           CALL CONVERT_DEFAULT_WALLS_TO_STL
64     
65     ! Bin the STL to the DES grid.
66           CALL BIN_FACETS_TO_DG
67     
68     ! Some functions for debugging.
69     !      CALL STL_DBG_WRITE_FACETS(BASE_STL)
70     !      CALL STL_DBG_WRITE_FACETS(BCWALLS_STL)
71     !      CALL STL_DBG_WRITE_FACETS(IMPRMBL_STL)
72     !      CALL STL_DBG_WRITE_FACETS(DEFAULT_STL)
73     !      CALL STL_DBG_WRITE_FACETS(ALL_STL)
74     !      CALL STL_DBG_WRITE_STL_FROM_DG(STL_TYPE=BASE_STL)
75     
76     ! Pre-procssing for the des in order to assign facets to grid cells.
77           WRITE(ERR_MSG,"('DES geometry pre-processing complete.')")
78           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
79     
80           RETURN
81           END SUBROUTINE DES_STL_PREPROCESSING
82     
83     
84     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
85     !                                                                      !
86     !  Subroutine: BIN_FACETS_TO_DG                                        !
87     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
88     !                                                                      !
89     !  Purpose:                                                            !
90     !                                                                      !
91     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
92           Subroutine BIN_FACETS_TO_DG
93     
94           use desgrid, only: DG_IJKSIZE2
95           use desgrid, only: DG_IEND2, DG_ISTART2
96           use desgrid, only: DG_JEND2, DG_JSTART2
97           use desgrid, only: DG_KEND2, DG_KSTART2
98           use desgrid, only: dg_dxinv, dg_dyinv, dg_dzinv
99     
100           use stl, only: FACETS_AT_DG
101     
102           use geometry, only: XLENGTH, YLENGTH, ZLENGTH, DO_K
103           use stl, only: N_FACETS_DES
104           use stl, only: VERTEX
105     
106           use stl, only: TOL_STL
107           use param1, only: ZERO, ONE
108     
109           use desgrid, only: DG_FUNIJK
110           use desgrid, only: IofPOS, JofPOS, KofPOS
111           use desgrid, only: dg_is_ON_myPE_plus1layers
112     
113           IMPLICIT NONE
114     
115     ! DES Grid cell index.
116           INTEGER :: IJK, IJK2
117     ! Loop counters:
118           INTEGER :: I1, I2, II  ! X-axis
119           INTEGER :: J1, J2, JJ  ! Y-axis
120           INTEGER :: K1, K2, KK  ! Z-axis
121           INTEGER :: NN          ! STLs
122     ! Generic accumulator
123           INTEGER :: COUNT_FAC
124     
125     ! Maximum and minimum extents of the indexed STL
126           DOUBLE PRECISION:: X1,Y1,Z1
127           DOUBLE PRECISION:: X2,Y2,Z2
128     
129     ! Allocate the data storage array.
130           IF(.not.allocated(FACETS_AT_DG)) &
131              allocate(FACETS_AT_DG(DG_IJKSIZE2))
132     
133           FACETS_AT_DG(:)%COUNT = 0
134     
135           DO NN = 1,N_FACETS_DES
136     
137              X1 = minval(VERTEX(1:3,1,NN))
138              X2 = maxval(VERTEX(1:3,1,NN))
139              Y1 = minval(VERTEX(1:3,2,NN))
140              Y2 = maxval(VERTEX(1:3,2,NN))
141              Z1 = minval(VERTEX(1:3,3,NN))
142              Z2 = maxval(VERTEX(1:3,3,NN))
143     
144              I1 = DG_IEND2
145              I2 = DG_ISTART2
146              IF(X2>=-TOL_STL .AND. X1<=XLENGTH+TOL_STL) THEN
147                 I1 = max(iofpos(X1)-1, dg_istart2)
148                 I2 = min(iofpos(X2)+1, dg_iend2)
149              ENDIF
150     
151              J1 = DG_JEND2
152              J2 = DG_JSTART2
153              IF(Y2>=-TOL_STL .AND. Y1<=YLENGTH+TOL_STL) THEN
154                 J1 = max(jofpos(Y1)-1, dg_jstart2)
155                 J2 = min(jofpos(Y2)+1, dg_jend2)
156              ENDIF
157     
158              K1 = DG_KEND2
159              K2 = DG_KSTART2
160              IF(DO_K) THEN
161                 IF(Z2>=-TOL_STL .AND. Z1<=ZLENGTH+TOL_STL) THEN
162                    K1 = max(kofpos(Z1)-1, dg_kstart2)
163                    K2 = min(kofpos(Z2)+1, dg_kend2)
164                 ENDIF
165              ENDIF
166     
167              DO KK=K1,K2
168              DO JJ=J1,J2
169              DO II=I1,I2
170                 IF(dg_is_ON_myPE_plus1layers(II,JJ,KK)) THEN
171                    IJK = DG_FUNIJK(II,JJ,KK)
172                    CALL ADD_FACET_FOR_DES(II,JJ,KK,IJK,NN)
173                 ENDIF
174              ENDDO
175              ENDDO
176              ENDDO
177     
178           ENDDO
179     
180           RETURN
181           END SUBROUTINE BIN_FACETS_TO_DG
182     
183     
184     
185     
186     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
187     !                                                                      !
188     !  Subroutine: ADD_FACET_FOR_DES                                       !
189     !  Author: Rahul Garg                                  Date: 24-Oct-13 !
190     !                                                                      !
191     !  Purpose: Add facets to DES grid cells.                              !
192     !                                                                      !
193     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
194           SUBROUTINE ADD_FACET_FOR_DES(I,J,K,IJK,N)
195     
196           use geometry, only: DO_K
197     
198           use desgrid, only: dg_dxinv, dg_xstart, dg_istart1
199           use desgrid, only: dg_dyinv, dg_ystart, dg_jstart1
200           use desgrid, only: dg_dzinv, dg_zstart, dg_kstart1
201     
202           use discretelement, only: MAX_RADIUS
203     
204           use stl, only: VERTEX
205     
206           use stl_functions_des, only: TRI_BOX_OVERLAP
207     
208           use param1, only: ZERO, HALF, ONE
209           use error_manager
210     
211           IMPLICIT NONE
212     
213     ! DES grid index and facet index
214           INTEGER, INTENT(IN) :: I,J,K,IJK, N
215     
216     ! Center of DES grid cell and half size. Note that a buffer is added to
217     ! the half size to make the cell appear a little larger. This ensures
218     ! that paricles near the edge 'see' STLs that are nearby but do not
219     ! directly intersect the DES grid cell contain the particle center.
220           DOUBLE PRECISION :: CENTER(3), HALFSIZE(3)
221     ! Flag: STL intersects the DES grid cell
222           LOGICAL :: OVERLAP
223     ! DES grid cell dimensions
224           DOUBLE PRECISION :: lDX, lDY, lDZ
225     ! Buffer to ensure all particle-STL collisions are captured.
226           DOUBLE PRECISION :: BUFFER
227     ! Legacy variable - should be removed
228           INTEGER ::  CURRENT_COUNT
229     
230           BUFFER = 1.1d0*MAX_RADIUS
231     
232           lDX = ONE/DG_DXINV
233           lDY = ONE/DG_DYINV
234           lDZ = ONE/DG_DZINV
235     
236           CENTER(1) = dg_xstart + (dble(I-dg_istart1)+HALF)*lDX
237           HALFSIZE(1) = HALF*lDX + BUFFER
238     
239           CENTER(2) = dg_ystart + (dble(J-dg_jstart1)+HALF)*lDY
240           HALFSIZE(2) = HALF*lDY + BUFFER
241     
242           IF(DO_K)THEN
243              CENTER(3) = dg_zstart + (dble(K-dg_kstart1)+HALF)*lDZ
244              HALFSIZE(3) = HALF*lDZ + BUFFER
245           ELSE
246              CENTER(3) = HALF*lDZ
247              HALFSIZE(3) = HALF*lDZ
248           ENDIF
249     
250           CALL TRI_BOX_OVERLAP(CENTER, HALFSIZE, VERTEX(:,:,N), OVERLAP)
251     
252           IF(OVERLAP) CALL ADD_FACET(IJK, N)
253     
254           RETURN
255           END SUBROUTINE ADD_FACET_FOR_DES
256     
257     
258     
259     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
260     !                                                                      !
261     !  Subroutine: ADD_FACET                                               !
262     !                                                                      !
263     !                                                                      !
264     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
265           SUBROUTINE ADD_FACET(IJK, FACET_ID)
266     
267           use stl, only: VERTEX
268           use stl, only: FACETS_AT_DG
269           use param1, only: ZERO
270     
271           implicit none
272     
273           INTEGER, INTENT(IN) :: IJK, facet_id
274     
275           INTEGER, ALLOCATABLE :: int_tmp(:)
276           DOUBLE PRECISION, ALLOCATABLE :: real_tmp(:)
277     
278           INTEGER :: lSIZE, II
279           DOUBLE PRECISION :: smallest_extent, min_temp, max_temp
280     
281     
282           IF(FACETS_AT_DG(IJK)%COUNT > 0) THEN
283     
284              DO II=1, FACETS_AT_DG(IJK)%COUNT
285                 IF(FACET_ID == FACETS_AT_DG(IJK)%ID(II)) RETURN
286              ENDDO
287     
288              FACETS_AT_DG(IJK)%COUNT = FACETS_AT_DG(IJK)%COUNT+1
289     
290              lSIZE = size(FACETS_AT_DG(IJK)%ID)
291              IF(FACETS_AT_DG(IJK)%COUNT +1> lSIZE) THEN
292                 allocate(int_tmp(2*lSIZE)); int_tmp=0
293                 int_tmp(1:lSIZE) = FACETS_AT_DG(IJK)%ID(1:lSIZE)
294                 call move_alloc(int_tmp,FACETS_AT_DG(IJK)%ID)
295     
296                 allocate(int_tmp(2*lSIZE)); int_tmp=0
297                 int_tmp(1:lSIZE) = FACETS_AT_DG(IJK)%DIR(1:lSIZE)
298                 call move_alloc(int_tmp, FACETS_AT_DG(IJK)%DIR)
299     
300                 allocate(real_tmp(2*lSIZE)); real_tmp=ZERO
301                 real_tmp(1:lSIZE) = FACETS_AT_DG(IJK)%MIN(1:lSIZE)
302                 call move_alloc(real_tmp, FACETS_AT_DG(IJK)%MIN)
303     
304                 allocate(real_tmp(2*lSIZE)); real_tmp=ZERO
305                 real_tmp(1:lSIZE) = FACETS_AT_DG(IJK)%MAX(1:lSIZE)
306                 call move_alloc(real_tmp, FACETS_AT_DG(IJK)%MAX)
307              ENDIF
308     
309           ELSE
310              FACETS_AT_DG(IJK)%COUNT = 1
311              IF(.not.allocated(FACETS_AT_DG(IJK)%ID)) &
312                 allocate(FACETS_AT_DG(IJK)%ID(4))
313              IF(.not.allocated(FACETS_AT_DG(IJK)%DIR)) &
314                 allocate(FACETS_AT_DG(IJK)%DIR(4))
315              IF(.not.allocated(FACETS_AT_DG(IJK)%MIN)) &
316                 allocate(FACETS_AT_DG(IJK)%MIN(4))
317              IF(.not.allocated(FACETS_AT_DG(IJK)%MAX)) &
318                 allocate(FACETS_AT_DG(IJK)%MAX(4))
319           ENDIF
320     
321           FACETS_AT_DG(IJK)%ID(FACETS_AT_DG(IJK)%COUNT) = FACET_ID
322     
323           SMALLEST_EXTENT = HUGE(0.0)
324     
325           DO II=1,3
326              MIN_TEMP = MINVAL(VERTEX(:,II,FACET_ID))
327              MAX_TEMP = MAXVAL(VERTEX(:,II,FACET_ID))
328              IF(ABS(MAX_TEMP - MIN_TEMP) < SMALLEST_EXTENT ) THEN
329                 FACETS_AT_DG(IJK)%DIR(FACETS_AT_DG(IJK)%COUNT) = II
330                 FACETS_AT_DG(IJK)%MIN(FACETS_AT_DG(IJK)%COUNT) = MIN_TEMP
331                 FACETS_AT_DG(IJK)%MAX(FACETS_AT_DG(IJK)%COUNT) = MAX_TEMP
332                 SMALLEST_EXTENT = ABS(MAX_TEMP - MIN_TEMP)
333              ENDIF
334           ENDDO
335     
336           RETURN
337           END SUBROUTINE ADD_FACET
338     
339     
340     
341     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
342     !                                                                      !
343     !  Subroutine: CONVERT_BC_WALLS_TO_STL                                 !
344     !  Author: J.Musser                                   Date: 03-Nov-15  !
345     !                                                                      !
346     !  Purpose: Convert user specified walls to STLs for particle-wall     !
347     !  collision detection.                                                !
348     !                                                                      !
349     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
350           Subroutine CONVERT_BC_WALLS_TO_STL
351     
352           use geometry, only: ZLENGTH, DO_K
353     
354           use bc, only: BC_DEFINED, BC_TYPE_ENUM, FREE_SLIP_WALL, NO_SLIP_WALL, PAR_SLIP_WALL
355           use bc, only: BC_I_w, BC_I_e
356           use bc, only: BC_J_s, BC_J_n
357           use bc, only: BC_K_b, BC_K_t
358     
359           use stl, only: N_FACETS_DES
360           use stl, only: STL_START, STL_END, BCWALLS_STL
361     
362           use discretelement, only: XE, YN, ZT
363     
364           use param, only: DIMENSION_BC
365           USE param1, only: ZERO
366     
367           IMPLICIT NONE
368     
369     ! Loop counter.
370           INTEGER :: BCV
371     ! Extents of the BC region with respect to the fluid grid.
372           DOUBLE PRECISION :: lXw, lXe, lYs, lYn, lZb, lZt
373     
374           STL_START(BCWALLS_STL)=N_FACETS_DES+1
375     
376           DO BCV=1, DIMENSION_BC
377              IF(.NOT.BC_DEFINED(BCV)) CYCLE
378     
379              IF(BC_TYPE_ENUM(BCV) == FREE_SLIP_WALL .OR.   &
380                 BC_TYPE_ENUM(BCV) == NO_SLIP_WALL   .OR.   &
381                 BC_TYPE_ENUM(BCV) == PAR_SLIP_WALL) THEN
382     
383                 lXw = XE(BC_I_w(BCV)-1); lXe = XE(BC_I_e(BCV))
384                 lYs = YN(BC_J_s(BCV)-1); lYn = YN(BC_J_n(BCV))
385                 IF(DO_K) THEN
386                    lZb = ZT(BC_K_b(BCV)-1); lZt = ZT(BC_K_t(BCV))
387                 ELSE
388                    lZb = ZERO; lZt = ZLENGTH
389                 ENDIF
390                 CALL GENERATE_STL_BOX(lXw, lXe, lYs, lYn, lZb, lZt)
391              ENDIF
392           ENDDO
393           STL_END(BCWALLS_STL)=N_FACETS_DES
394     
395           RETURN
396           END SUBROUTINE CONVERT_BC_WALLS_TO_STL
397     
398     
399     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
400     !                                                                      !
401     !  Subroutine: CONVERT_IMPERMEABLE_IS_TO_STL                           !
402     !  Author: J.Musser                                   Date: 03-Nov-15  !
403     !                                                                      !
404     !  Purpose: Convert user specified impermeable surfaces to STLs for    !
405     !  particle-wall collision detection.                                  !
406     !                                                                      !
407     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
408           SUBROUTINE CONVERT_IMPERMEABLE_IS_TO_STL
409     
410           use geometry, only: DO_K, ZLENGTH
411     
412           use is, only: IS_DEFINED, IS_TYPE
413           use is, only: IS_I_w, IS_I_e
414           use is, only: IS_J_s, IS_J_n
415           use is, only: IS_K_b, IS_K_t
416     
417           use stl, only: N_FACETS_DES
418           use stl, only: STL_START, STL_END, IMPRMBL_STL
419           use discretelement, only: XE, YN, ZT
420     
421           use param, only: DIMENSION_IS
422           USE param1, only: ZERO
423     
424           use error_manager
425     
426           IMPLICIT NONE
427     
428     ! Loop counter.
429           INTEGER :: ISV
430     ! Extents of the BC region with respect to the fluid grid.
431           DOUBLE PRECISION :: lXw, lXe, lYs, lYn, lZb, lZt
432     
433           STL_START(IMPRMBL_STL)=N_FACETS_DES+1
434     
435           DO ISV=1, DIMENSION_IS
436              IF(.NOT.IS_DEFINED(ISV)) CYCLE
437     
438              IF(trim(IS_TYPE(ISV)) == 'IMPERMEABLE') THEN
439     
440                 lXw = XE(IS_I_w(ISV)-1); lXe = XE(IS_I_e(ISV))
441                 lYs = YN(IS_J_s(ISV)-1); lYn = YN(IS_J_n(ISV))
442                 IF(DO_K) THEN
443                    lZb = ZT(IS_K_b(ISV)-1); lZt = ZT(IS_K_t(ISV))
444                 ELSE
445                    lZb = ZERO; lZt = ZLENGTH
446                 ENDIF
447     
448                 CALL GENERATE_STL_BOX(lXw, lXe, lYs, lYn, lZb, lZt)
449              ELSE
450                 CALL INIT_ERR_MSG('CONVERT_IMPERMEABLE_IS_TO_STL')
451                 WRITE(ERR_MSG,1000) ISV, IS_TYPE(ISV)
452                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
453              ENDIF
454           ENDDO
455     
456           STL_END(IMPRMBL_STL)=N_FACETS_DES
457     
458      1000 FORMAT("Error 1000: DES simulations do not support the ",/       &
459              'specified IS TYPE:',/3x,'IS: ',I3,/3x,'IS_TYPE=',A)
460     
461           RETURN
462           END SUBROUTINE CONVERT_IMPERMEABLE_IS_TO_STL
463     
464     
465     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
466     !                                                                      !
467     !  Subroutine: CONVERT_DEFAULT_WALLS_TO_STL                            !
468     !  Author: J.Musser                                   Date: 03-Nov-15  !
469     !                                                                      !
470     !  Purpose: Convert user specified walls to STLs for particle-wall     !
471     !  collision detection.                                                !
472     !                                                                      !
473     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
474           Subroutine CONVERT_DEFAULT_WALLS_TO_STL
475     
476           USE geometry, only: DO_K
477           USE geometry, only: XLENGTH, YLENGTH, ZLENGTH
478           use stl, only: VERTEX, NORM_FACE
479           use stl, only: N_FACETS_DES
480           use stl, only: STL_START, STL_END, DEFAULT_STL
481     
482           use discretelement, only: DES_PERIODIC_WALLS_X
483           use discretelement, only: DES_PERIODIC_WALLS_Y
484           use discretelement, only: DES_PERIODIC_WALLS_Z
485     
486           USE param1, only: ZERO, ONE
487     
488           IMPLICIT NONE
489     
490           STL_START(DEFAULT_STL)=N_FACETS_DES+1
491     
492     ! West Face
493           IF(.NOT.DES_PERIODIC_WALLS_X)THEN
494              N_FACETS_DES = N_FACETS_DES+1
495              VERTEX(1,:,N_FACETS_DES) = (/ZERO, ZERO, ZERO/)
496              VERTEX(2,:,N_FACETS_DES) = (/ZERO, 2*YLENGTH, ZERO/)
497              VERTEX(3,:,N_FACETS_DES) = (/ZERO, ZERO, 2*ZLENGTH/)
498              NORM_FACE(:,N_FACETS_DES) = (/ONE, ZERO, ZERO/)
499     
500     ! East Face
501              N_FACETS_DES = N_FACETS_DES+1
502              VERTEX(1,:,N_FACETS_DES) = (/XLENGTH, ZERO, ZERO/)
503              VERTEX(2,:,N_FACETS_DES) = (/XLENGTH, 2*YLENGTH, ZERO/)
504              VERTEX(3,:,N_FACETS_DES) = (/XLENGTH, ZERO, 2*ZLENGTH/)
505              NORM_FACE(:,N_FACETS_DES) = (/-ONE, ZERO, ZERO/)
506           ENDIF
507     
508     ! South Face
509           IF(.NOT.DES_PERIODIC_WALLS_Y)THEN
510              N_FACETS_DES = N_FACETS_DES+1
511              VERTEX(1,:,N_FACETS_DES) = (/ZERO, ZERO, ZERO/)
512              VERTEX(2,:,N_FACETS_DES) = (/2*XLENGTH, ZERO, ZERO/)
513              VERTEX(3,:,N_FACETS_DES) = (/ZERO, ZERO, 2*ZLENGTH/)
514              NORM_FACE(:,N_FACETS_DES) = (/ZERO, ONE, ZERO/)
515     
516     ! North Face
517              N_FACETS_DES = N_FACETS_DES+1
518              VERTEX(1,:,N_FACETS_DES) = (/ZERO, YLENGTH, ZERO/)
519              VERTEX(2,:,N_FACETS_DES) = (/2*XLENGTH, YLENGTH, ZERO/)
520              VERTEX(3,:,N_FACETS_DES) = (/ZERO, YLENGTH, 2*ZLENGTH/)
521              NORM_FACE(:,N_FACETS_DES) = (/ZERO, -ONE, ZERO/)
522           ENDIF
523     
524     ! Bottom Face
525           IF(.NOT.DES_PERIODIC_WALLS_Z .AND. DO_K) THEN
526              N_FACETS_DES = N_FACETS_DES+1
527              VERTEX(1,:,N_FACETS_DES) = (/ZERO, ZERO, ZERO/)
528              VERTEX(2,:,N_FACETS_DES) = (/2*XLENGTH, ZERO, ZERO/)
529              VERTEX(3,:,N_FACETS_DES) = (/ZERO, 2*YLENGTH, ZERO/)
530              NORM_FACE(:,N_FACETS_DES) = (/ZERO, ZERO, ONE/)
531     
532     ! Top Face
533              N_FACETS_DES = N_FACETS_DES+1
534              VERTEX(1,:,N_FACETS_DES) = (/ZERO, ZERO, ZLENGTH/)
535              VERTEX(2,:,N_FACETS_DES) = (/2*XLENGTH, ZERO, ZLENGTH/)
536              VERTEX(3,:,N_FACETS_DES) = (/ZERO, 2*YLENGTH, ZLENGTH/)
537              NORM_FACE(:,N_FACETS_DES) = (/ZERO, ZERO, -ONE/)
538           ENDIF
539     
540           STL_END(DEFAULT_STL)=N_FACETS_DES
541     
542           RETURN
543           END SUBROUTINE CONVERT_DEFAULT_WALLS_TO_STL
544     
545           END MODULE STL_PREPROC_DES
546     
547     
548     
549     
550     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
551     !                                                                      !
552     !  Subroutine: GENERATE_STL_BOX                                        !
553     !  Author: J.Musser                                   Date: 03-Nov-15  !
554     !                                                                      !
555     !  Purpose: Given the six corners of a box, create the 12 STLs needed  !
556     !  to define the geometry.                                             !
557     !                                                                      !
558     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
559           SUBROUTINE GENERATE_STL_BOX(pXw, pXe, pYs, pYn, pZb, pZt)
560     
561           use stl, only: VERTEX, NORM_FACE
562           use stl, only: N_FACETS_DES
563     
564           use geometry, only: DO_K
565     
566           use param1, only: ZERO, ONE
567     
568           IMPLICIT NONE
569     
570           DOUBLE PRECISION, INTENT(IN) :: pXw, pXe, pYs, pYn, pZb, pZt
571     
572     ! West Face
573           N_FACETS_DES = N_FACETS_DES+1
574           VERTEX(1,:,N_FACETS_DES) = (/pXw, pYs, pZb/)
575           VERTEX(2,:,N_FACETS_DES) = (/pXw, pYn, pZb/)
576           VERTEX(3,:,N_FACETS_DES) = (/pXw, pYn, pZt/)
577           NORM_FACE(:,N_FACETS_DES) = (/-ONE, ZERO, ZERO/)
578     
579           IF(DO_K)THEN
580              N_FACETS_DES = N_FACETS_DES+1
581              VERTEX(1,:,N_FACETS_DES) = (/pXw, pYs, pZt/)
582              VERTEX(2,:,N_FACETS_DES) = (/pXw, pYs, pZb/)
583              VERTEX(3,:,N_FACETS_DES) = (/pXw, pYn, pZt/)
584              NORM_FACE(:,N_FACETS_DES) = (/-ONE, ZERO, ZERO/)
585           ENDIF
586     
587     ! East Face
588           N_FACETS_DES = N_FACETS_DES+1
589           VERTEX(3,:,N_FACETS_DES) = (/pXe, pYs, pZb/)
590           VERTEX(2,:,N_FACETS_DES) = (/pXe, pYn, pZb/)
591           VERTEX(1,:,N_FACETS_DES) = (/pXe, pYn, pZt/)
592           NORM_FACE(:,N_FACETS_DES) = (/ONE, ZERO, ZERO/)
593     
594           IF(DO_K) THEN
595              N_FACETS_DES = N_FACETS_DES+1
596              VERTEX(3,:,N_FACETS_DES) = (/pXe, pYs, pZt/)
597              VERTEX(2,:,N_FACETS_DES) = (/pXe, pYs, pZb/)
598              VERTEX(1,:,N_FACETS_DES) = (/pXe, pYn, pZt/)
599              NORM_FACE(:,N_FACETS_DES) = (/ONE, ZERO, ZERO/)
600           ENDIF
601     
602     ! South Face
603           N_FACETS_DES = N_FACETS_DES+1
604           VERTEX(1,:,N_FACETS_DES) = (/pXw, pYs, pZb/)
605           VERTEX(2,:,N_FACETS_DES) = (/pXe, pYs, pZb/)
606           VERTEX(3,:,N_FACETS_DES) = (/pXw, pYs, pZt/)
607           NORM_FACE(:,N_FACETS_DES) = (/ZERO, -ONE, ZERO/)
608     
609           IF(DO_K) THEN
610              N_FACETS_DES = N_FACETS_DES+1
611              VERTEX(1,:,N_FACETS_DES) = (/pXe, pYs, pZt/)
612              VERTEX(2,:,N_FACETS_DES) = (/pXe, pYs, pZb/)
613              VERTEX(3,:,N_FACETS_DES) = (/pXw, pYs, pZb/)
614              NORM_FACE(:,N_FACETS_DES) = (/ZERO, -ONE, ZERO/)
615           ENDIF
616     
617     ! North Face
618           N_FACETS_DES = N_FACETS_DES+1
619           VERTEX(3,:,N_FACETS_DES) = (/pXw, pYn, pZb/)
620           VERTEX(2,:,N_FACETS_DES) = (/pXe, pYn, pZb/)
621           VERTEX(1,:,N_FACETS_DES) = (/pXw, pYn, pZt/)
622           NORM_FACE(:,N_FACETS_DES) = (/ZERO, ONE, ZERO/)
623     
624           IF(DO_K) THEN
625              N_FACETS_DES = N_FACETS_DES+1
626              VERTEX(3,:,N_FACETS_DES) = (/pXe, pYn, pZt/)
627              VERTEX(2,:,N_FACETS_DES) = (/pXe, pYn, pZb/)
628              VERTEX(1,:,N_FACETS_DES) = (/pXw, pYn, pZb/)
629              NORM_FACE(:,N_FACETS_DES) = (/ZERO, ONE, ZERO/)
630           ENDIF
631     
632     ! Bottom Face
633           IF(DO_K)THEN
634              N_FACETS_DES = N_FACETS_DES+1
635              VERTEX(1,:,N_FACETS_DES) = (/pXw, pYs, pZb/)
636              VERTEX(2,:,N_FACETS_DES) = (/pXe, pYs, pZb/)
637              VERTEX(3,:,N_FACETS_DES) = (/pXe, pYn, pZb/)
638              NORM_FACE(:,N_FACETS_DES) = (/ZERO, ZERO, -ONE/)
639     
640              N_FACETS_DES = N_FACETS_DES+1
641              VERTEX(1,:,N_FACETS_DES) = (/pXe, pYn, pZb/)
642              VERTEX(2,:,N_FACETS_DES) = (/pXw, pYn, pZb/)
643              VERTEX(3,:,N_FACETS_DES) = (/pXw, pYs, pZb/)
644              NORM_FACE(:,N_FACETS_DES) = (/ZERO, ZERO, -ONE/)
645     
646     ! Top Face
647              N_FACETS_DES = N_FACETS_DES+1
648              VERTEX(3,:,N_FACETS_DES) = (/pXw, pYs, pZb/)
649              VERTEX(2,:,N_FACETS_DES) = (/pXe, pYs, pZb/)
650              VERTEX(1,:,N_FACETS_DES) = (/pXe, pYn, pZb/)
651              NORM_FACE(:,N_FACETS_DES) = (/ZERO, ZERO, ONE/)
652     
653              N_FACETS_DES = N_FACETS_DES+1
654              VERTEX(3,:,N_FACETS_DES) = (/pXe, pYn, pZb/)
655              VERTEX(2,:,N_FACETS_DES) = (/pXw, pYn, pZb/)
656              VERTEX(1,:,N_FACETS_DES) = (/pXw, pYs, pZb/)
657              NORM_FACE(:,N_FACETS_DES) = (/ZERO, ZERO, ONE/)
658           ENDIF
659     
660           RETURN
661           END SUBROUTINE GENERATE_STL_BOX
662