File: N:\mfix\model\cartesian_grid\get_stl_data.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: GET_MSH_DATA                                           C
4     !  Purpose: reads face vertices and normal vectors from an MSH file    C
5     !           (generated by Gambit). After reading the geometry, the     C
6     !           data is converted into same format as STL and the          C
7     !           pre-processing switches to STL procedure                   C
8     !                                                                      C
9     !  Author: Jeff Dietiker                              Date: 30-JAN-09  C
10     !  Reviewer:                                          Date: **-***-**  C
11     !                                                                      C
12     !  Revision Number:                                                    C
13     !  Purpose:                                                            C
14     !  Author:                                            Date: dd-mmm-yy  C
15     !  Reviewer:                                          Date: dd-mmm-yy  C
16     !                                                                      C
17     !  Literature/Document References:                                     C
18     !                                                                      C
19     !  Variables referenced:                                               C
20     !                                                                      C
21     !  Variables modified:                                                 C
22     !                                                                      C
23     !  Local variables:                                                    C
24     !                                                                      C
25     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
26     !
27           SUBROUTINE GET_MSH_DATA
28     
29     !-----------------------------------------------
30     !   M o d u l e s
31     !-----------------------------------------------
32     
33           USE bc
34           USE compar
35           USE constant
36           USE exit, only: mfix_exit
37           USE fldvar
38           USE funits
39           USE mpi_utility
40           USE param
41           USE param1
42           USE physprop
43           USE progress_bar
44           USE quadric
45           USE run
46           USE rxns
47           USE scalars
48           USE stl
49           USE vtk
50     
51           IMPLICIT NONE
52     
53           INTEGER,PARAMETER :: MAX_POINTS = 10000000    ! Currently limited to 10 Million, increase if needed
54           INTEGER,PARAMETER :: MAX_ZONES  = 1000        ! Currently limited to 1,000     , increase if needed
55     
56           INTEGER ::I,D,NN,NF
57     
58           INTEGER   :: GRID_DIMENSION, N_POINTS,N_FACES,N_FACE_ZONES
59           INTEGER   :: ZONE,ZONE_ID,N1,N2
60           INTEGER   :: PPFACE
61     
62           LOGICAL :: ALL_POINTS_READ,ALL_FACES_READ
63     
64           INTEGER, ALLOCATABLE, DIMENSION(:,:) ::POINT_ZONE_INFO,FACE_ZONE_INFO
65           DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: POINT_COORDS
66     
67           CHARACTER(LEN=18), ALLOCATABLE, DIMENSION(:,:) :: BC_LABEL_TEXT
68     
69           LOGICAL, ALLOCATABLE, DIMENSION(:) :: BC_ASSIGNED
70     
71           CHARACTER(LEN=18),DIMENSION(10) :: BUFF_CHAR
72     
73           INTEGER :: P1,P2,P3,P4
74           DOUBLE PRECISION, DIMENSION(3) :: NORMAL, VECTMP, VECTMP2
75           DOUBLE PRECISION ::NORM
76     
77           DOUBLE PRECISION ::v1x,v1y,v1z
78           DOUBLE PRECISION ::v2x,v2y,v2z
79           DOUBLE PRECISION ::v3x,v3y,v3z
80           DOUBLE PRECISION ::ABSTRANS
81     
82           LOGICAL :: PRESENT
83     
84           INTEGER :: BC_LABELS_TO_READ, BC_LABELS_READ,NFZ,ZID,L2
85           INTEGER :: N_BC_IGNORED
86           INTEGER :: N_QUAD2TRI
87     
88           ALLOCATE(POINT_COORDS(3,MAX_POINTS))
89           ALLOCATE(POINT_ZONE_INFO(MAX_ZONES,3))
90           ALLOCATE(FACE_ZONE_INFO(MAX_ZONES,3))
91           ALLOCATE(BC_LABEL_TEXT(MAX_ZONES,2))
92           ALLOCATE(BC_ASSIGNED(MAX_ZONES))
93     
94           WRITE(*,2000) 'READING geometry from geometry.msh...'
95     
96           INQUIRE(FILE='geometry.msh',EXIST=PRESENT)
97           IF(.NOT.PRESENT) THEN
98              IF(MyPE == PE_IO) THEN
99                 WRITE(*,"('(PE ',I3,'): input data file, ',A11,' is missing: run aborted')") &
100                 myPE,'geometry.msh'
101              ENDIF
102              CALL MFIX_EXIT(MYPE)
103           ENDIF
104     !
105     !
106     !     OPEN geometry.msh ASCII FILE
107     !
108           OPEN(UNIT=333, FILE='geometry.msh', STATUS='OLD', ERR=910,CONVERT='BIG_ENDIAN')
109     
110           IF(MyPE == PE_IO) WRITE(*,2000)'MSH file opened. Starting reading data...'
111     
112           OPEN(UNIT=444, FILE='checkgeometry.stl',CONVERT='BIG_ENDIAN')
113           write(444,*)'solid vcg'
114     
115           CALL SKIP(333,3)
116     
117     ! Reading grid dimension
118           READ(333,*) (BUFF_CHAR(I),I=1,2)
119     
120           CALL TEXT_HEX2INT(BUFF_CHAR(2)(1:1),GRID_DIMENSION)
121     
122           CALL SKIP(333,1)
123     
124     ! Reading Points
125     
126           READ(333,*) (BUFF_CHAR(I),I=1,4)
127     
128           CALL TEXT_HEX2INT(BUFF_CHAR(4),N_POINTS)
129     
130           ZONE = 1
131           ALL_POINTS_READ = .FALSE.
132     
133           DO WHILE(.NOT.ALL_POINTS_READ)
134     
135              READ(333,*) (BUFF_CHAR(I),I=1,4)
136              CALL TEXT_HEX2INT(BUFF_CHAR(2),ZONE_ID)
137              CALL TEXT_HEX2INT(BUFF_CHAR(3),N1)
138              CALL TEXT_HEX2INT(BUFF_CHAR(4),N2)
139     
140              POINT_ZONE_INFO(ZONE,1) = ZONE_ID
141              POINT_ZONE_INFO(ZONE,2) = N1
142              POINT_ZONE_INFO(ZONE,3) = N2
143     
144              DO NN = N1,N2
145                 READ(333,*) (POINT_COORDS(D,NN),D=1,GRID_DIMENSION)
146              ENDDO
147     
148              IF(N2/=N_POINTS) THEN
149                 ZONE = ZONE + 1
150              ELSE
151                 ALL_POINTS_READ = .TRUE.
152     !            WRITE(*,*)' ALL POINTS READ SUCCESSFULLY.'
153               ENDIF
154     
155           ENDDO
156     
157           CALL SKIP(333,3)
158     
159     ! Reading FACES
160     
161           READ(333,*) (BUFF_CHAR(I),I=1,4)
162     
163           CALL TEXT_HEX2INT(BUFF_CHAR(3),N_FACES)
164     !      WRITE(*,*)' READING FACES ...'
165     
166           ZONE = 1
167           ALL_FACES_READ = .FALSE.
168     
169           BC_LABELS_TO_READ = 0
170           N_BC_IGNORED = 0
171           NF=0
172           N_QUAD2TRI=0
173     
174           DO WHILE(.NOT.ALL_FACES_READ)
175     
176              READ(333,*) (BUFF_CHAR(I),I=1,4)
177     
178              CALL TEXT_HEX2INT(BUFF_CHAR(1)(5:),ZONE_ID)
179              CALL TEXT_HEX2INT(BUFF_CHAR(2),N1)
180              CALL TEXT_HEX2INT(BUFF_CHAR(3),N2)
181     
182              FACE_ZONE_INFO(ZONE,1) = ZONE_ID
183              FACE_ZONE_INFO(ZONE,2) = N1
184              FACE_ZONE_INFO(ZONE,3) = N2
185     
186              BC_LABELS_TO_READ = BC_LABELS_TO_READ + 1
187     
188              IF(IS_CG(BC_TYPE_ENUM(ZONE_ID))) THEN
189     
190                 DO NN = N1,N2
191                    READ(333,*)PPFACE
192                    IF(PPFACE<3.OR.PPFACE>4) THEN
193                       IF(MyPE == PE_IO) WRITE(*,*)PPFACE, 'POINTS PER FACE. EACH FACE MUST HAVE 3 OR 4 POINTS.'
194                       CALL MFIX_EXIT(MYPE)
195     
196                    ELSEIF(PPFACE==3) THEN
197                       BACKSPACE(333)
198                       READ(333,*) (BUFF_CHAR(I),I=1,4)
199                       CALL TEXT_HEX2INT(BUFF_CHAR(2),P1)
200                       CALL TEXT_HEX2INT(BUFF_CHAR(3),P2)
201                       CALL TEXT_HEX2INT(BUFF_CHAR(4),P3)
202     
203                       V1x = POINT_COORDS(1,P1)
204                       V1y = POINT_COORDS(2,P1)
205                       V1z = POINT_COORDS(3,P1)
206     
207                       V2x = POINT_COORDS(1,P2)
208                       V2y = POINT_COORDS(2,P2)
209                       V2z = POINT_COORDS(3,P2)
210     
211                       V3x = POINT_COORDS(1,P3)
212                       V3y = POINT_COORDS(2,P3)
213                       V3z = POINT_COORDS(3,P3)
214     
215                       VECTMP  = POINT_COORDS(:,P2)-POINT_COORDS(:,P1)
216                       VECTMP2 = POINT_COORDS(:,P3)-POINT_COORDS(:,P1)
217                       NORMAL = CROSS_PRODUCT(VECTMP,VECTMP2)
218     
219                       NORM = sqrt(dot_product(normal(:),normal(:)))
220                       NORMAL = NORMAL / NORM
221     
222                       NF = NF + 1
223     
224                       ! Save and Reverse unit vector if needed (this will switch fluid and blocked cells)
225                       NORM_FACE(:,NF) = NORMAL*OUT_MSH_VALUE
226     
227     
228                       VERTEX(1,1,NF) = SCALE_MSH*V1x + TX_MSH
229                       VERTEX(1,2,NF) = SCALE_MSH*V1y + TY_MSH
230                       VERTEX(1,3,NF) = SCALE_MSH*V1z + TZ_MSH
231     
232                       VERTEX(2,1,NF) = SCALE_MSH*V2x + TX_MSH
233                       VERTEX(2,2,NF) = SCALE_MSH*V2y + TY_MSH
234                       VERTEX(2,3,NF) = SCALE_MSH*V2z + TZ_MSH
235     
236                       VERTEX(3,1,NF) = SCALE_MSH*V3x + TX_MSH
237                       VERTEX(3,2,NF) = SCALE_MSH*V3y + TY_MSH
238                       VERTEX(3,3,NF) = SCALE_MSH*V3z + TZ_MSH
239     
240                       BC_ID_STL_FACE(NF) = ZONE_ID
241     
242     
243                       write(444,*) '   facet normal ', NORM_FACE(:,NF)
244                       write(444,*) '      outer loop'
245                       write(444,*) '         vertex ', VERTEX(1,1:3,NF)
246                       write(444,*) '         vertex ', VERTEX(2,1:3,NF)
247                       write(444,*) '         vertex ', VERTEX(3,1:3,NF)
248                       write(444,*) '      endloop'
249                       write(444,*) '   endfacet'
250     
251                    ELSEIF(PPFACE==4) THEN
252                       BACKSPACE(333)
253                       READ(333,*) (BUFF_CHAR(I),I=1,5)
254                       CALL TEXT_HEX2INT(BUFF_CHAR(2),P1)
255                       CALL TEXT_HEX2INT(BUFF_CHAR(3),P2)
256                       CALL TEXT_HEX2INT(BUFF_CHAR(4),P3)
257                       CALL TEXT_HEX2INT(BUFF_CHAR(5),P4)
258     
259                       N_QUAD2TRI = N_QUAD2TRI + 1
260     
261     ! Splitting Quad face 1-2-3-4 into two triangles 1-2-3 and 1-3-4
262     
263     ! First triangle 1-2-3
264     
265                       V1x = POINT_COORDS(1,P1)
266                       V1y = POINT_COORDS(2,P1)
267                       V1z = POINT_COORDS(3,P1)
268     
269                       V2x = POINT_COORDS(1,P2)
270                       V2y = POINT_COORDS(2,P2)
271                       V2z = POINT_COORDS(3,P2)
272     
273                       V3x = POINT_COORDS(1,P3)
274                       V3y = POINT_COORDS(2,P3)
275                       V3z = POINT_COORDS(3,P3)
276     
277                       VECTMP  = POINT_COORDS(:,P2)-POINT_COORDS(:,P1)
278                       VECTMP2 = POINT_COORDS(:,P3)-POINT_COORDS(:,P1)
279                       NORMAL = CROSS_PRODUCT(VECTMP,VECTMP2)
280     
281                       NORM = sqrt(dot_product(normal(:),normal(:)))
282                       NORMAL = NORMAL / NORM
283     
284                       NF = NF + 1
285     
286                       ! Save and Reverse unit vector if needed (this will switch fluid and blocked cells)
287                       NORM_FACE(:,NF) = NORMAL*OUT_MSH_VALUE
288     
289                       VERTEX(1,1,NF) = SCALE_MSH*V1x + TX_MSH
290                       VERTEX(1,2,NF) = SCALE_MSH*V1y + TY_MSH
291                       VERTEX(1,3,NF) = SCALE_MSH*V1z + TZ_MSH
292     
293                       VERTEX(2,1,NF) = SCALE_MSH*V2x + TX_MSH
294                       VERTEX(2,2,NF) = SCALE_MSH*V2y + TY_MSH
295                       VERTEX(2,3,NF) = SCALE_MSH*V2z + TZ_MSH
296     
297                       VERTEX(3,1,NF) = SCALE_MSH*V3x + TX_MSH
298                       VERTEX(3,2,NF) = SCALE_MSH*V3y + TY_MSH
299                       VERTEX(3,3,NF) = SCALE_MSH*V3z + TZ_MSH
300     
301                       BC_ID_STL_FACE(NF) = ZONE_ID
302     
303     
304                       write(444,*) '   facet normal ', NORM_FACE(:,NF)
305                       write(444,*) '      outer loop'
306                       write(444,*) '         vertex ', VERTEX(1,1:3,NF)
307                       write(444,*) '         vertex ', VERTEX(2,1:3,NF)
308                       write(444,*) '         vertex ', VERTEX(3,1:3,NF)
309                       write(444,*) '      endloop'
310                       write(444,*)'   endfacet'
311     
312     ! Second triangle 1-3-4
313     
314                       V1x = POINT_COORDS(1,P1)
315                       V1y = POINT_COORDS(2,P1)
316                       V1z = POINT_COORDS(3,P1)
317     
318                       V2x = POINT_COORDS(1,P3)
319                       V2y = POINT_COORDS(2,P3)
320                       V2z = POINT_COORDS(3,P3)
321     
322                       V3x = POINT_COORDS(1,P4)
323                       V3y = POINT_COORDS(2,P4)
324                       V3z = POINT_COORDS(3,P4)
325     
326                       VECTMP  = POINT_COORDS(:,P3)-POINT_COORDS(:,P1)
327                       VECTMP2 = POINT_COORDS(:,P4)-POINT_COORDS(:,P1)
328                       NORMAL = CROSS_PRODUCT(VECTMP,VECTMP2)
329     
330                       NORM = sqrt(dot_product(normal(:),normal(:)))
331                       NORMAL = NORMAL / NORM
332     
333                       NF = NF + 1
334     
335                       ! Save and Reverse unit vector if needed (this will switch fluid and blocked cells)
336                       NORM_FACE(:,NF) = NORMAL*OUT_MSH_VALUE
337     
338                       VERTEX(1,1,NF) = SCALE_MSH*V1x + TX_MSH
339                       VERTEX(1,2,NF) = SCALE_MSH*V1y + TY_MSH
340                       VERTEX(1,3,NF) = SCALE_MSH*V1z + TZ_MSH
341     
342                       VERTEX(2,1,NF) = SCALE_MSH*V2x + TX_MSH
343                       VERTEX(2,2,NF) = SCALE_MSH*V2y + TY_MSH
344                       VERTEX(2,3,NF) = SCALE_MSH*V2z + TZ_MSH
345     
346                       VERTEX(3,1,NF) = SCALE_MSH*V3x + TX_MSH
347                       VERTEX(3,2,NF) = SCALE_MSH*V3y + TY_MSH
348                       VERTEX(3,3,NF) = SCALE_MSH*V3z + TZ_MSH
349     
350                       BC_ID_STL_FACE(NF) = ZONE_ID
351     
352                       write(444,*) '   facet normal ', NORM_FACE(:,NF)
353                       write(444,*) '      outer loop'
354                       write(444,*) '         vertex ', VERTEX(1,1:3,NF)
355                       write(444,*) '         vertex ', VERTEX(2,1:3,NF)
356                       write(444,*) '         vertex ', VERTEX(3,1:3,NF)
357                       write(444,*) '      endloop'
358                       write(444,*)'   endfacet'
359     
360                    ENDIF
361                 ENDDO
362     
363              ELSE
364                 N_BC_IGNORED = N_BC_IGNORED + 1
365     !            WRITE(*,*)'BOUNDARY CONDITION FOR ZONE',ZONE_ID,' IS NOT DEFINED.'
366     !            WRITE(*,*)'THIS ZONE IS IGNORED'
367                 CALL SKIP(333,N2-N1+1)
368              ENDIF
369     
370     
371              IF(N2/=N_FACES) THEN
372                 ZONE = ZONE + 1
373                 ALL_FACES_READ = .FALSE.
374                 CALL SKIP(333,1)
375              ELSE
376                 ALL_FACES_READ = .TRUE.
377     !            WRITE(*,*)' ALL FACES READ SUCCESSFULLY.'
378                 N_FACE_ZONES = ZONE
379               ENDIF
380     
381           ENDDO
382     
383           write(444,*)'endsolid vcg'
384     
385           close(444)
386     
387           IF(MyPE == PE_IO) THEN
388              WRITE(*,*) ' The file check_geometry.stl was sucessfully written.'
389              WRITE(*,*) ' This is the equivalent of geometry.msh in STL format,'
390              WRITE(*,*) ' and is provided for convenience (it is not used).'
391           ENDIF
392     
393     ! Reading Boundary condition labels
394     
395           BC_LABELS_READ = 0
396           BC_ASSIGNED = .FALSE.
397     
398     
399           DO WHILE(BC_LABELS_READ < BC_LABELS_TO_READ)
400     
401              READ(333,*) BUFF_CHAR(1)
402     
403              IF(BUFF_CHAR(1)=='(45') THEN
404                 BACKSPACE(333)
405                 READ(333,*) (BUFF_CHAR(I),I=1,4)
406     
407                 CALL TEXT_DEC2INT(BUFF_CHAR(2),ZONE_ID)
408     
409                 DO ZONE = 1,N_FACE_ZONES
410                    IF(FACE_ZONE_INFO(ZONE,1)==ZONE_ID) THEN
411                       BC_LABELS_READ = BC_LABELS_READ + 1
412                       BC_LABEL_TEXT(ZONE,1) = TRIM(BUFF_CHAR(3))
413                       BC_LABEL_TEXT(ZONE,2) = TRIM(BUFF_CHAR(4))
414                       IF(IS_CG(BC_TYPE_ENUM(ZONE_ID))) THEN
415                          BC_ASSIGNED(ZONE) = .TRUE.
416                       ENDIF
417                    ENDIF
418                 ENDDO
419     
420              ENDIF
421           ENDDO
422     
423           IF(MyPE == PE_IO) THEN
424     
425              WRITE(*,*)' Summary of data read from geometry.msh file:'
426              WRITE(*,*)'======================================================================'
427              WRITE(*,*)' DIMENSION   POINTS       FACES       ZONES'
428              WRITE(*,*)'======================================================================'
429     
430              WRITE(*,1020) GRID_DIMENSION,N_POINTS,N_FACES,N_FACE_ZONES
431              WRITE(*,*)''
432              WRITE(*,*)' BOUNDARY CONDITION DETECTED (INFO EXTRACTED FROM .MSH FILE):'
433              WRITE(*,*)'======================================================================'
434              WRITE(*,*)'  ZONE  BC_TYPE(MFIX)  FACES   BC_TYPE(GAMBIT)      BC LABEL'
435              WRITE(*,*)'======================================================================'
436     
437              N_FACES = 0
438              DO ZONE = 1,N_FACE_ZONES
439                 IF(BC_ASSIGNED(ZONE)) THEN
440                    ZID = FACE_ZONE_INFO(ZONE,1)
441                    NFZ = FACE_ZONE_INFO(ZONE,3)-FACE_ZONE_INFO(ZONE,2) + 1
442                    N_FACES = N_FACES + NFZ
443                    L2=LEN(TRIM(BC_LABEL_TEXT(ZONE,2)))-4
444                    WRITE(*,1000) ZID,BC_TYPE_ENUM(ZID),NFZ,BC_LABEL_TEXT(ZONE,1),BC_LABEL_TEXT(ZONE,2)(1:L2)
445                 ENDIF
446              ENDDO
447     
448              WRITE(*,*)''
449     
450              DO ZONE = 1,N_FACE_ZONES
451                 IF(.NOT.BC_ASSIGNED(ZONE)) THEN
452                    ZID = FACE_ZONE_INFO(ZONE,1)
453                    NFZ = FACE_ZONE_INFO(ZONE,3)-FACE_ZONE_INFO(ZONE,2) + 1
454                    L2=LEN(TRIM(BC_LABEL_TEXT(ZONE,2)))-4
455                    WRITE(*,1000) ZID,'NOT USED',NFZ,BC_LABEL_TEXT(ZONE,1),BC_LABEL_TEXT(ZONE,2)(1:L2)
456                 ENDIF
457              ENDDO
458     
459              WRITE(*,*)'======================================================================'
460              WRITE(*,*)' PLEASE VERIFY THAT BOUNDARY CONDITIONS ARE CORRECTLY ASSIGNED.'
461              WRITE(*,*)' MODIFY BC_TYPE IN mfix.dat IF NECESSARY.'
462              WRITE(*,*)''
463     
464           ENDIF
465     
466     
467           XMIN_MSH = MINVAL(VERTEX(:,1,1:N_FACES))
468           XMAX_MSH = MAXVAL(VERTEX(:,1,1:N_FACES))
469           YMIN_MSH = MINVAL(VERTEX(:,2,1:N_FACES))
470           YMAX_MSH = MAXVAL(VERTEX(:,2,1:N_FACES))
471           ZMIN_MSH = MINVAL(VERTEX(:,3,1:N_FACES))
472           ZMAX_MSH = MAXVAL(VERTEX(:,3,1:N_FACES))
473     
474           IF(MyPE == PE_IO) THEN
475              WRITE(*,2000)'MSH file successfully read.'
476              WRITE(*,*)' Total number of faces used as boundary faces =',N_FACES
477              IF(N_QUAD2TRI>0) WRITE(*,*)' Number of quad faces split into triangles    =',N_QUAD2TRI
478     
479     
480              WRITE(*,*)' RANGE OF MSH FILE:'
481              IF(SCALE_MSH/=ONE) THEN
482                 WRITE(*,5000)' AFTER SCALING BY A FACTOR OF ',SCALE_MSH
483              ENDIF
484              ABSTRANS = dabs(TX_MSH)+dabs(TY_MSH)+dabs(TZ_MSH)
485              IF(ABSTRANS>TOL_MSH) THEN
486                 WRITE(*,3000)' AFTER TRANSLATION OF (X,Y,Z)=',TX_MSH,TY_MSH,TZ_MSH
487              ENDIF
488              WRITE(*,4000)'X-RANGE = ', XMIN_MSH,XMAX_MSH
489              WRITE(*,4000)'Y-RANGE = ', YMIN_MSH,YMAX_MSH
490              WRITE(*,4000)'Z-RANGE = ', ZMIN_MSH,ZMAX_MSH
491              WRITE(*,4000)''
492           ENDIF
493     
494           XMIN_MSH = XMIN_MSH - 10.0*TOL_MSH
495           XMAX_MSH = XMAX_MSH + 10.0*TOL_MSH
496           YMIN_MSH = YMIN_MSH - 10.0*TOL_MSH
497           YMAX_MSH = YMAX_MSH + 10.0*TOL_MSH
498           ZMIN_MSH = ZMIN_MSH - 10.0*TOL_MSH
499           ZMAX_MSH = ZMAX_MSH + 10.0*TOL_MSH
500     
501           N_FACETS = N_FACES
502     
503           N_FACETS = NF
504     
505           XMIN_STL = XMIN_MSH
506           XMAX_STL = XMAX_MSH
507           YMIN_STL = YMIN_MSH
508           YMAX_STL = YMAX_MSH
509           ZMIN_STL = ZMIN_MSH
510           ZMAX_STL = ZMAX_MSH
511     
512           OUT_STL_VALUE = OUT_MSH_VALUE
513     
514           CLOSE(333)
515     
516           OPEN(UNIT=444, FILE='msgeometry.stl',CONVERT='BIG_ENDIAN')
517     
518              DO ZONE = 1,N_FACE_ZONES
519                 IF(BC_ASSIGNED(ZONE)) THEN
520                    ZID = FACE_ZONE_INFO(ZONE,1)
521                    L2=LEN(TRIM(BC_LABEL_TEXT(ZONE,2)))-4
522     !               WRITE(*,1000) ZID,BC_TYPE(ZID),NFZ,BC_LABEL_TEXT(ZONE,1),BC_LABEL_TEXT(ZONE,2)(1:L2)
523                    print*,'=======>  Zone ID= ', ZID,BC_LABEL_TEXT(ZONE,2)(1:L2)
524                    write(444,6000) 'solid', BC_LABEL_TEXT(ZONE,2)(1:L2),ZID
525                    DO NF=1,N_FACETS
526                       IF(BC_ID_STL_FACE(NF) == ZID) THEN
527     
528                       write(444,*) '   facet normal ', NORM_FACE(:,NF)
529                       write(444,*) '      outer loop'
530                       write(444,*) '         vertex ', VERTEX(1,1:3,NF)
531                       write(444,*) '         vertex ', VERTEX(2,1:3,NF)
532                       write(444,*) '         vertex ', VERTEX(3,1:3,NF)
533                       write(444,*) '      endloop'
534                       write(444,*)'   endfacet'
535     
536     
537                       ENDIF
538                    ENDDO
539     
540                    write(444,6000) 'endsolid', BC_LABEL_TEXT(ZONE,2)(1:L2),ZID
541     
542                 ENDIF
543              ENDDO
544     
545           close(444)
546     
547           DEALLOCATE(POINT_COORDS)
548           DEALLOCATE(POINT_ZONE_INFO)
549           DEALLOCATE(FACE_ZONE_INFO)
550           DEALLOCATE(BC_LABEL_TEXT)
551           DEALLOCATE(BC_ASSIGNED)
552     
553           RETURN
554     
555     !======================================================================
556     !     HERE IF AN ERROR OCCURED OPENNING/READING THE FILE
557     !======================================================================
558     !
559      910  CONTINUE
560           WRITE (*, 1500)
561           CALL MFIX_EXIT(myPE)
562      920  CONTINUE
563           WRITE (*, 1600)
564           CALL MFIX_EXIT(myPE)
565      930  CONTINUE
566           WRITE (*, 1700)
567           CALL MFIX_EXIT(myPE)
568     
569      1000 FORMAT(1X,I4,7X,A8,I8,1X,A20,1X,A20)
570      1010 FORMAT(1X,I4,5X,A20,1X,A20)
571      1020 FORMAT(5X,I3,2X,3(I10,2X))
572     !
573      1500 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
574           'Unable to open stl file',/1X,70('*')/)
575      1600 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
576           'Error while reading stl file',/1X,70('*')/)
577      1700 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
578           'End of file reached while reading stl file',/1X,70('*')/)
579      2000 FORMAT(1X,A)
580      2010 FORMAT(1X,A,I4,A)
581      3000 FORMAT(1X,A,'(',F10.4,';',F10.4,';',F10.4,')')
582      4000 FORMAT(1X,A,F10.4,' to ',F10.4)
583      5000 FORMAT(1X,A,F10.4)
584      6000 FORMAT(A,1X,A,'_',I4.4)
585           END SUBROUTINE GET_MSH_DATA
586     
587     
588           SUBROUTINE SKIP(FILE_UNIT,N_SKIP)
589             IMPLICIT NONE
590             INTEGER, INTENT(IN) ::FILE_UNIT,N_SKIP
591             INTEGER :: I
592             DO I = 1,N_SKIP
593                READ(FILE_UNIT,*)
594             END DO
595             RETURN
596           END SUBROUTINE SKIP
597     
598           SUBROUTINE TEXT_HEX2INT(STRING,INT)
599     
600           CHARACTER(LEN=10)  :: INTERNAL
601           CHARACTER(LEN=*)   :: STRING
602           CHARACTER(LEN=10)  :: CLEANED_STRING
603           INTEGER :: INT
604     
605           CLEANED_STRING = STRING
606           IF(STRING(1:1)=='(') CLEANED_STRING = STRING(2:)
607     
608           WRITE(INTERNAL,fmt='(A10)') TRIM(CLEANED_STRING)
609           READ(INTERNAL,FMT='(Z10)') INT
610     
611           RETURN
612     
613           END SUBROUTINE TEXT_HEX2INT
614     
615           SUBROUTINE TEXT_DEC2INT(STRING,INT)
616     
617           CHARACTER(LEN=10)  :: INTERNAL
618           CHARACTER(LEN=*)   :: STRING
619           CHARACTER(LEN=10)  :: CLEANED_STRING
620           INTEGER :: INT
621     
622           CLEANED_STRING = STRING
623           IF(STRING(1:1)=='(') CLEANED_STRING = STRING(2:)
624     
625           WRITE(INTERNAL,fmt='(A10)') TRIM(CLEANED_STRING)
626           READ(INTERNAL,FMT='(I10)') INT
627     
628           RETURN
629     
630           END SUBROUTINE TEXT_DEC2INT
631     
632     
633     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
634     !                                                                      C
635     !  Module name: GET_STL_DATA                                           C
636     !  Purpose: reads face verticesa and normal vectors from an STL file   C
637     !                                                                      C
638     !  Author: Jeff Dietiker                              Date: 30-JAN-09  C
639     !  Reviewer:                                          Date: **-***-**  C
640     !                                                                      C
641     !  Revision Number:                                                    C
642     !  Purpose:                                                            C
643     !  Author:                                            Date: dd-mmm-yy  C
644     !  Reviewer:                                          Date: dd-mmm-yy  C
645     !                                                                      C
646     !  Literature/Document References:                                     C
647     !                                                                      C
648     !  Variables referenced:                                               C
649     !                                                                      C
650     !  Variables modified:                                                 C
651     !                                                                      C
652     !  Local variables:                                                    C
653     !                                                                      C
654     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
655     !
656           SUBROUTINE GET_STL_DATA
657     
658     !-----------------------------------------------
659     !   M o d u l e s
660     !-----------------------------------------------
661     
662           USE bc
663           USE compar
664           USE constant
665           USE exit, only: mfix_exit
666           USE fldvar
667           USE funits
668           USE mpi_utility
669           USE param
670           USE param1
671           USE physprop
672           USE progress_bar
673           USE quadric
674           USE run
675           USE rxns
676           USE scalars
677           USE stl
678           USE vtk
679           IMPLICIT NONE
680     
681           INTEGER :: NN,IGNORED_FACETS
682           LOGICAL :: PRESENT,KEEP_READING,IGNORE_CURRENT_FACET
683           DOUBLE PRECISION ::v1x,v1y,v1z
684           DOUBLE PRECISION ::v2x,v2y,v2z
685           DOUBLE PRECISION ::v3x,v3y,v3z
686           DOUBLE PRECISION ::x12,y12,z12,x13,y13,z13,dp,d12,d13
687           DOUBLE PRECISION ::cos_angle,cos_small_angle
688           DOUBLE PRECISION ::n1,n2,n3,NORM
689           DOUBLE PRECISION ::ABSTRANS
690           CHARACTER(LEN=32) ::TEST_CHAR,BUFF_CHAR
691           CHARACTER(LEN=255) ::geometryfile(0:DIMENSION_BC)
692     
693           INTEGER :: BCV,NUMBER_OF_GEOMETRY_FILES,NUMBER_OF_BC_PATCHES
694           INTEGER :: BC_PATCH(0:DIMENSION_BC),L2,NF1,NF2
695           LOGICAL :: BC_PATCH_FOUND_IN_STL(DIMENSION_BC)
696     
697           CHARACTER(LEN=100) :: FNAME
698           integer :: stl_unit, nf
699     
700           BC_PATCH_FOUND_IN_STL = .FALSE.
701     
702           geometryfile(0) = 'geometry.stl'
703           NUMBER_OF_GEOMETRY_FILES = 0
704           NUMBER_OF_BC_PATCHES     = 0
705     
706     !     DETERMINE WHICH BOUNDARY CONDITIONS NEED STL FILE
707     
708           IF(MyPE == PE_IO) THEN
709              WRITE(*,100)'From Get_Stl_Data: Analysing BCs in mfix.dat...'
710              WRITE(*,120)'BC_ID','BC_TYPE'
711           ENDIF
712           DO BCV = 1, DIMENSION_BC
713     
714              IF(IS_CG(BC_TYPE_ENUM(BCV))) THEN
715     
716                 NUMBER_OF_GEOMETRY_FILES = NUMBER_OF_GEOMETRY_FILES + 1
717                 NUMBER_OF_BC_PATCHES     = NUMBER_OF_BC_PATCHES + 1
718     
719     
720                 BC_PATCH(NUMBER_OF_GEOMETRY_FILES) = BCV
721                 WRITE(geometryfile(NUMBER_OF_GEOMETRY_FILES),200) 'geometry_',BCV
722     
723                 IF(MyPE == PE_IO) WRITE(*,130)BCV,BC_TYPE(BCV)
724     
725              ENDIF
726           ENDDO
727     
728           IF(MyPE == PE_IO) WRITE(*,110)'Number of CG BCs in mfix.dat = ',NUMBER_OF_BC_PATCHES
729     
730     100  FORMAT(1X,A)
731     110  FORMAT(1X,A,I6)
732     120  FORMAT(1X,A6,4X,A7)
733     130  FORMAT(1X,I6,4X,A6)
734     140  FORMAT(1X,A,I6,A)
735     200  FORMAT(A,I4.4,".stl")
736     
737     ! VERIFY THAT THERE IS AT LEAST ONE STL FILE TO READ
738     
739           IF(NUMBER_OF_BC_PATCHES==0) THEN
740              IF(MyPE == PE_IO) THEN
741                 WRITE(*,100) 'ERROR: NO CARTESIAN GRID BOUNDARY CONDITION SPECIFIED.'
742                 WRITE(*,100) 'AT LEAST ONE BC_TYPE MUST START WITH CG (FOR EXAMPLE CG_NSW)'
743                 WRITE(*,100) 'RUN ABORTED'
744                 CALL MFIX_EXIT(MYPE)
745              ENDIF
746     
747           ELSEIF(NUMBER_OF_BC_PATCHES==1) THEN
748     
749              INQUIRE(FILE='geometry.stl',EXIST=PRESENT)
750              IF(PRESENT) THEN
751     
752     !            IF(STL_BC_ID==BC_PATCH(1)) THEN
753                    IF(MyPE == PE_IO) THEN
754                       WRITE(*,110) 'The file geometry.stl exists and the following  BC patch was found: ', BC_PATCH(1)
755                    ENDIF
756     
757                    NF1 = 1
758                    NF2 = 1
759                    geometryfile(NF1)='geometry.stl'
760     
761     !            ENDIF
762              ENDIF
763     
764           ELSE  ! More than one CG BC type
765              INQUIRE(FILE='geometry.stl',EXIST=PRESENT)
766              IF(PRESENT) THEN
767                 IF(MyPE == PE_IO) THEN
768                    WRITE(*,100) 'The file geometry.stl exists and several CG BC types are defined.'
769                    WRITE(*,100) 'All BC patches will be read from geometry.stl.'
770                 ENDIF
771     
772                 geometryfile(0)='geometry.stl'
773                 NF1 = 0
774                 NF2 = 0                       ! This invokes a special treatment
775                                               ! to reading multiple solids
776     
777              ELSE
778     
779                 NF1 = 1
780                 NF2 = NUMBER_OF_GEOMETRY_FILES
781                 IF(MyPE == PE_IO) THEN
782                    WRITE(*,100) 'The file geometry.stl does not exist and several CG_BC types are defined.'
783                    WRITE(*,100) 'Each BC patch will be read from geometry_BCID.stl.'
784                    WRITE(*,100) 'where BCID is 4-paded integer'
785                 ENDIF
786     
787              ENDIF
788     
789     
790     
791           ENDIF
792     
793     
794     ! START READING EACH STL FILE, ONE FOR EACH BC_TYPE
795     
796           N_FACETS = 0
797           IGNORED_FACETS = 0
798     
799     
800           DO NN = NF1, NF2
801     
802              IF(MyPE == PE_IO) WRITE(*,2000) 'Reading geometry from '//TRIM(geometryfile(NN))//' ...'
803     
804              INQUIRE(FILE=TRIM(geometryfile(NN)),EXIST=PRESENT)
805              IF(.NOT.PRESENT) THEN
806                 IF(MyPE == PE_IO) THEN
807                    WRITE(*,"('(PE ',I3,'): input data file, ',A11,' is missing: run aborted')") &
808                    myPE,TRIM(geometryfile(NN))
809                 ENDIF
810                 CALL MFIX_EXIT(MYPE)
811              ENDIF
812     
813     
814     !
815     !     OPEN geometry.stl ASCII FILE
816     !
817              OPEN(UNIT=333, FILE=TRIM(geometryfile(NN)), STATUS='OLD', ERR=910,CONVERT='BIG_ENDIAN')
818     
819              IF(MyPE == PE_IO) WRITE(*,2000)'STL file opened. Starting reading data...'
820     
821              KEEP_READING = .TRUE.
822     
823              DO WHILE(KEEP_READING)
824     
825                 READ(333,*,ERR=920,END=930) TEST_CHAR
826     !            print *,'TEST_CHAR=',TEST_CHAR
827                 IF(TRIM(TEST_CHAR) == 'solid'.AND.NN==0) THEN
828     
829                    BACKSPACE(333)
830     
831                    READ(333,*,ERR=920,END=930) BUFF_CHAR,BUFF_CHAR
832     
833                    L2=LEN(TRIM(BUFF_CHAR))-3
834     
835                    READ(BUFF_CHAR(L2:L2+4),fmt='(I4.4)') BC_PATCH(NN)
836     
837                    IF(MyPE == PE_IO)  WRITE(*,140) 'Found BC patch #', BC_PATCH(NN), ' in STL file.'
838     
839                    IF(.NOT.IS_CG(BC_TYPE_ENUM(BC_PATCH(NN)))) THEN
840                       IF(MyPE == PE_IO)  THEN
841                          WRITE(*,110) 'This BC patch does not mach a CG BC in mfix.dat:',BC_PATCH(NN)
842                          WRITE(*,100)'Please Correct mfix.dat and/or stl file amd try again'
843                          CALL MFIX_EXIT(myPE)
844                       ENDIF
845                    ENDIF
846     
847                    BC_PATCH_FOUND_IN_STL(BC_PATCH(NN)) = .TRUE.
848     
849                 ELSEIF(TRIM(TEST_CHAR) == 'facet') THEN
850     
851                    BACKSPACE(333)
852                    IGNORE_CURRENT_FACET = .FALSE.
853     
854                    READ(333,*,ERR=920,END=930) BUFF_CHAR,BUFF_CHAR,N1,N2,N3  ! Read unit normal vector
855                    READ(333,*,ERR=920,END=930)
856                    READ(333,*,ERR=920,END=930) BUFF_CHAR, V1x,V1y,V1z
857                    READ(333,*,ERR=920,END=930) BUFF_CHAR, V2x,V2y,V2z
858                    READ(333,*,ERR=920,END=930) BUFF_CHAR, V3x,V3y,V3z
859     
860                    N1 = N1 * OUT_STL_VALUE  ! Reverse unit vector if needed (this will switch fluid and blocked cells)
861                    N2 = N2 * OUT_STL_VALUE
862                    N3 = N3 * OUT_STL_VALUE
863     
864                    x12 = V2x - V1x
865                    y12 = V2y - V1y
866                    z12 = V2z - V1z
867     
868                    x13 = V3x - V1x
869                    y13 = V3y - V1y
870                    z13 = V3z - V1z
871     
872     
873                    dp  = x12*x13 + y12*y13 + z12*z13
874                    d12 = sqrt(x12**2+y12**2+z12**2)
875                    d13 = sqrt(x13**2+y13**2+z13**2)
876     
877                    IF((d12*d13)>TOL_STL) THEN
878                       cos_angle = dp/(d12*d13)
879                    ELSE
880                       cos_angle = ONE
881                    ENDIF
882     
883                    cos_small_angle = dcos(STL_SMALL_ANGLE / 180.0 * PI)
884     
885                    IF(DABS(cos_angle)>cos_small_angle) THEN
886                       IGNORE_CURRENT_FACET = .TRUE.    ! Ignore small facets
887                    ENDIF
888     
889                    NORM = sqrt(N1**2+N2**2+N3**2)
890     
891                    IF(NORM>TOL_STL) THEN
892                       N1 = N1 /NORM
893                       N2 = N2 /NORM
894                       N3 = N3 /NORM
895                    ELSE
896                       IGNORE_CURRENT_FACET = .TRUE.  ! Ignore facets with zero normal vector
897                    ENDIF
898     
899     
900                    IF(IGNORE_CURRENT_FACET) THEN
901                       IGNORED_FACETS = IGNORED_FACETS + 1
902                    ELSE                                                      ! Save vertex coordinates for valid facets
903                                                                              ! and performs translation
904                       N_FACETS = N_FACETS + 1
905     
906                       NORM_FACE(1,N_FACETS) = N1
907                       NORM_FACE(2,N_FACETS) = N2
908                       NORM_FACE(3,N_FACETS) = N3
909     
910                       VERTEX(1,1,N_FACETS) = SCALE_STL*V1x + TX_STL
911                       VERTEX(1,2,N_FACETS) = SCALE_STL*V1y + TY_STL
912                       VERTEX(1,3,N_FACETS) = SCALE_STL*V1z + TZ_STL
913     
914                       VERTEX(2,1,N_FACETS) = SCALE_STL*V2x + TX_STL
915                       VERTEX(2,2,N_FACETS) = SCALE_STL*V2y + TY_STL
916                       VERTEX(2,3,N_FACETS) = SCALE_STL*V2z + TZ_STL
917     
918                       VERTEX(3,1,N_FACETS) = SCALE_STL*V3x + TX_STL
919                       VERTEX(3,2,N_FACETS) = SCALE_STL*V3y + TY_STL
920                       VERTEX(3,3,N_FACETS) = SCALE_STL*V3z + TZ_STL
921     
922                       BC_ID_STL_FACE(N_FACETS) = BC_PATCH(NN)
923     
924     
925                    ENDIF
926     
927     !            ELSEIF(TRIM(TEST_CHAR) == 'endsolid') THEN
928     
929     !               KEEP_READING = .FALSE.
930     
931                 ENDIF
932     
933              ENDDO
934     
935     
936     930      IF(MyPE == PE_IO)  WRITE(*,100) 'Done reading file.'
937     
938              CLOSE(333)
939     
940     
941     
942     
943           ENDDO
944     
945     
946           IF(myPE==0.AND.NF2==0) THEN
947              DO NN = 1,NUMBER_OF_BC_PATCHES
948     
949                 IF(.NOT.BC_PATCH_FOUND_IN_STL(BC_PATCH(NN))) THEN
950                    WRITE (*, 140)'Did not find BC patch: ',BC_PATCH(NN) , ' in stl file'
951                    WRITE(*,100)'Please correct mfix.dat and/or stl file amd try again'
952                    CALL MFIX_EXIT(myPE)
953                 ENDIF
954     
955              ENDDO
956           ENDIF
957     
958     
959     
960     
961     
962           XMIN_STL = MINVAL(VERTEX(:,1,1:N_FACETS))
963           XMAX_STL = MAXVAL(VERTEX(:,1,1:N_FACETS))
964           YMIN_STL = MINVAL(VERTEX(:,2,1:N_FACETS))
965           YMAX_STL = MAXVAL(VERTEX(:,2,1:N_FACETS))
966           ZMIN_STL = MINVAL(VERTEX(:,3,1:N_FACETS))
967           ZMAX_STL = MAXVAL(VERTEX(:,3,1:N_FACETS))
968     
969           IF(MyPE == PE_IO) THEN
970              WRITE(*,2000)'STL file(s) successfully read.'
971              WRITE(*,110)'Total number of facets read =',N_FACETS + IGNORED_FACETS
972              WRITE(*,110)'Number of valid facets      =',N_FACETS
973              WRITE(*,110)'Number of ignored facets    =',IGNORED_FACETS
974              WRITE(*,100)'Geometry range from stl file:'
975              IF(SCALE_STL/=ONE) THEN
976                 WRITE(*,5000)'Scaling factor:',SCALE_STL
977              ENDIF
978              ABSTRANS = dabs(TX_STL)+dabs(TY_STL)+dabs(TZ_STL)
979              IF(ABSTRANS>TOL_STL) THEN
980                 WRITE(*,3000)'Translation vector (X,Y,Z)=',TX_STL,TY_STL,TZ_STL
981              ENDIF
982              WRITE(*,4000)'x-range = ', XMIN_STL,XMAX_STL
983              WRITE(*,4000)'y-range = ', YMIN_STL,YMAX_STL
984              WRITE(*,4000)'z-range = ', ZMIN_STL,ZMAX_STL
985              WRITE(*,4000)''
986           ENDIF
987     
988           XMIN_STL = XMIN_STL - 10.0*TOL_STL
989           XMAX_STL = XMAX_STL + 10.0*TOL_STL
990           YMIN_STL = YMIN_STL - 10.0*TOL_STL
991           YMAX_STL = YMAX_STL + 10.0*TOL_STL
992           ZMIN_STL = ZMIN_STL - 10.0*TOL_STL
993           ZMAX_STL = ZMAX_STL + 10.0*TOL_STL
994     
995     
996     
997     !      IF(XMIN_STL<ZERO) XMIN_STL=ZERO
998     !      IF(XMAX_STL>XLENGTH) XMAX_STL=XLENGTH
999     !      IF(YMIN_STL<ZERO) YMIN_STL=ZERO
1000     !      IF(YMAX_STL>YLENGTH) YMAX_STL=YLENGTH
1001     !      IF(ZMIN_STL<ZERO) ZMIN_STL=ZERO
1002     !      IF(ZMAX_STL>ZLENGTH) ZMAX_STL=ZLENGTH
1003     
1004           IF(mype.eq.pe_io.and..false.) then
1005              WRITE(fname,'(A,"_FACETS_READ", ".stl")') &
1006              TRIM(RUN_NAME)
1007              open(stl_unit, file = fname, form='formatted',convert='big_endian')
1008              write(stl_unit,*)'solid vcg'
1009     
1010     
1011              do nf = 1, n_facets
1012     
1013                 write(stl_unit,*) '   facet normal ', NORM_FACE(:,NF)
1014                 write(stl_unit,*) '      outer loop'
1015                 write(stl_unit,*) '         vertex ', VERTEX(1,1:3,NF)
1016                 write(stl_unit,*) '         vertex ', VERTEX(2,1:3,NF)
1017                 write(stl_unit,*) '         vertex ', VERTEX(3,1:3,NF)
1018                 write(stl_unit,*) '      endloop'
1019                 write(stl_unit,*)'   endfacet'
1020              enddo
1021     
1022     
1023              write(stl_unit,*)'endsolid vcg'
1024              close(stl_unit, status = 'keep')
1025              IF(MyPE == PE_IO) THEN
1026                 WRITE(*,100) 'The file FACETS_READ.stl was sucessfully written.'
1027                 WRITE(*,100) 'and is provided for convenience (it is not used).'
1028              ENDIF
1029     
1030           endif
1031     
1032           RETURN
1033     
1034     !======================================================================
1035     !     HERE IF AN ERROR OCCURED OPENNING/READING THE FILE
1036     !======================================================================
1037     !
1038      910  CONTINUE
1039           WRITE (*, 1500)
1040           CALL MFIX_EXIT(myPE)
1041      920  CONTINUE
1042           WRITE (*, 1600)
1043           CALL MFIX_EXIT(myPE)
1044     ! 930  CONTINUE
1045     !     WRITE (*, 1700)
1046     !     CALL MFIX_EXIT(myPE)
1047     !
1048      1500 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
1049           'Unable to open stl file',/1X,70('*')/)
1050      1600 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
1051           'Error while reading stl file',/1X,70('*')/)
1052      1700 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
1053           'End of file reached while reading stl file',/1X,70('*')/)
1054      2000 FORMAT(1X,A)
1055      2010 FORMAT(1X,A,I4,A)
1056      3000 FORMAT(1X,A,'(',F10.4,';',F10.4,';',F10.4,')')
1057      4000 FORMAT(1X,A,F10.4,' to ',F10.4)
1058      5000 FORMAT(1X,A,F10.4)
1059           END SUBROUTINE GET_STL_DATA
1060     
1061     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1062     !                                                                      C
1063     !  Module name: EVAL_STL_FCT_AT                                        C
1064     !  Purpose: Assigns a value to f_stl at any NODE of cell IJK           C
1065     !                                                                      C
1066     !  Author: Jeff Dietiker                              Date: 30-JAN-09  C
1067     !  Reviewer:                                          Date: **-***-**  C
1068     !                                                                      C
1069     !  Revision Number:                                                    C
1070     !  Purpose:                                                            C
1071     !  Author:                                            Date: dd-mmm-yy  C
1072     !  Reviewer:                                          Date: dd-mmm-yy  C
1073     !                                                                      C
1074     !  Literature/Document References:                                     C
1075     !                                                                      C
1076     !  Variables referenced:                                               C
1077     !                                                                      C
1078     !  Variables modified:                                                 C
1079     !                                                                      C
1080     !  Local variables:                                                    C
1081     !                                                                      C
1082     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1083     !
1084           SUBROUTINE EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,NODE,f_stl,CLIP_FLAG,BCID)
1085     
1086     !-----------------------------------------------
1087     !   M o d u l e s
1088     !-----------------------------------------------
1089     
1090           USE param
1091           USE param1
1092           USE physprop
1093           USE fldvar
1094           USE run
1095           USE scalars
1096           USE funits
1097           USE rxns
1098           USE compar
1099           USE mpi_utility
1100           USE progress_bar
1101     !      USE quadric
1102           USE cutcell
1103           USE stl
1104           USE functions
1105           IMPLICIT NONE
1106     
1107           INTEGER :: BCID
1108           LOGICAL :: CLIP_FLAG
1109           DOUBLE PRECISION :: f_stl
1110     
1111           CHARACTER (LEN=*) :: TYPE_OF_CELL
1112           INTEGER :: IJK,IJKC,NODE
1113     
1114           IF(NODE==15) print*,'Warning: eval stl at node 15'
1115     
1116           IF(N_FACETS < 1) RETURN
1117     
1118           f_stl = UNDEFINED
1119     
1120           IJKC = IJK_OF_NODE(NODE)
1121     
1122           IF(IJKC<1.OR.IJKC>DIMENSION_3) THEN
1123     
1124     !         print*,'myPE        =',myPE
1125     !         print*,'DIMENSION_3 =',DIMENSION_3
1126     !         print*,'IJKC        =',IJKC
1127     !         print*,'IJK         =',IJK
1128     !         print*,'I,J,K       =',I_OF(IJK),J_OF(IJK),K_OF(IJK)
1129     !         print*,'NODE        =',NODE
1130     
1131     !     Leave f_stl as undefined if we are out of bounds
1132     !     This can happen at I=1 and NODE=4 for example
1133              RETURN
1134     
1135           ENDIF
1136     
1137     !!      IF(NODE/=15.AND.NODE/=0) THEN
1138     !!         IF(SNAP(IJKC)) THEN
1139     !!            f_stl = ZERO
1140     !!!            print*,'stlfc snapped=',IJKC,NODE,F_AT(IJKC)
1141     !!            RETURN
1142     !!         ENDIF
1143     !!      ENDIF
1144     
1145     
1146           f_stl = F_AT(IJKC)
1147     
1148           RETURN
1149     
1150     
1151           CLIP_FLAG = .TRUE.
1152     
1153     
1154           RETURN
1155     
1156     
1157           END SUBROUTINE EVAL_STL_FCT_AT
1158     
1159     
1160     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1161     !                                                                      C
1162     !  Module name: intersect_line_with_facet                              C
1163     !  Purpose: Finds the intersection between a facet                     C
1164     !           and the line (xa,ya,za) and (xb,yb,zb).                    C
1165     !                                                                      C
1166     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1167     !  Reviewer:                                          Date:            C
1168     !                                                                      C
1169     !  Revision Number #                                  Date: ##-###-##  C
1170     !  Author: #                                                           C
1171     !  Purpose: #                                                          C
1172     !                                                                      C
1173     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1174       SUBROUTINE INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,FACET,INTERSECT_FLAG,xc,yc,zc)
1175     
1176           USE param
1177           USE param1
1178           USE parallel
1179           USE constant
1180           USE run
1181           USE toleranc
1182           USE geometry
1183           USE indices
1184           USE compar
1185           USE sendrecv
1186           USE quadric
1187           USE STL
1188     
1189           IMPLICIT NONE
1190           DOUBLE PRECISION:: xa,ya,za,xb,yb,zb,xc,yc,zc
1191           INTEGER :: FACET
1192           DOUBLE PRECISION:: NFx,NFy,NFz,nabx,naby,nabz
1193           DOUBLE PRECISION :: dot_denom,dot_num
1194           DOUBLE PRECISION :: VP1Ax,VP1Ay,VP1Az
1195           DOUBLE PRECISION :: Px,Py,Pz
1196           DOUBLE PRECISION :: tt
1197           DOUBLE PRECISION :: d_ac,d_bc
1198           LOGICAL :: INSIDE_FACET,INTERSECT_FLAG
1199     
1200     
1201     !======================================================================
1202     ! This subroutine tries to find the intersection of a line AB with one
1203     ! of the facets defined in the stl file
1204     !
1205     !
1206     ! 1) Verify that intersection is possible. Only facets and line AB that are
1207     !    not parallel are acceptable candidates
1208     ! 2) Find intersection point P of line AB with plane containing the facet.
1209     !    Valid intersection point must be between A and B to continue.
1210     ! 3) Verify that point P is inside triangular facet
1211     !======================================================================
1212     
1213           INTERSECT_FLAG = .FALSE.
1214     
1215     !======================================================================
1216     !  Facet normal vector (normalized)
1217     !======================================================================
1218           NFx = NORM_FACE(1,FACET)
1219           NFy = NORM_FACE(2,FACET)
1220           NFz = NORM_FACE(3,FACET)
1221     
1222     !======================================================================
1223     !  AB vector (NOT normalized)
1224     !======================================================================
1225           nabx = xb - xa
1226           naby = yb - ya
1227           nabz = zb - za
1228     
1229           dot_denom = NFx*nabx + NFy*naby + NFz*nabz
1230     !======================================================================
1231     ! 1) Verify that intersection is possible. Only facets and line AB that are
1232     !    not parallel are acceptable candidates
1233     !======================================================================
1234           IF(dabs(dot_denom)<TOL_STL) THEN
1235     
1236              INTERSECT_FLAG = .FALSE.              ! No intersection (facet nornal
1237                                                    ! and AB are perpendicular
1238              RETURN
1239     
1240           ELSE
1241     
1242     !======================================================================
1243     ! 2) Find intersection point P of line AB with plane containing the facet
1244     !    Line AB is parametrized with parameter t (from t=0 at A to t=1 at B)
1245     !======================================================================
1246              VP1Ax = VERTEX(1,1,FACET) - xa
1247              VP1Ay = VERTEX(1,2,FACET) - ya
1248              VP1Az = VERTEX(1,3,FACET) - za
1249     
1250              dot_num = NFx*VP1Ax + NFy*VP1Ay + NFz*VP1Az
1251     
1252              tt = dot_num / dot_denom
1253     
1254     !======================================================================
1255     ! 3) Verify that point P is inside triangular facet
1256     !    Facet is parametrized with (u,v), a point is inside triange if:
1257     !    u   >= 0  , and
1258     !    v   >= 0  , and
1259     !    u+v <= 1
1260     !======================================================================
1261              IF((tt>=ZERO).AND.(tt<=ONE)) THEN       ! Intersection between A and B
1262                                                    ! Now test if intersection point is inside triangle
1263     
1264                 Px = xa + tt*nabx
1265                 Py = ya + tt*naby
1266                 Pz = za + tt*nabz
1267     
1268                 CALL IS_POINT_INSIDE_FACET(Px,Py,Pz,FACET,INSIDE_FACET)
1269     
1270                 IF(INSIDE_FACET) THEN
1271                    INTERSECT_FLAG = .TRUE.              ! Valid intersection
1272                    xc = Px                              ! point P is inside triangle)
1273                    yc = Py
1274                    zc = Pz
1275                 ELSE
1276                    INTERSECT_FLAG = .FALSE.             ! Invalid intersection
1277                    RETURN                               ! point P is outside triangle)
1278                 ENDIF
1279     
1280              ELSE
1281     
1282                 INTERSECT_FLAG = .FALSE.           ! No intersection between A and B
1283                 RETURN
1284     
1285              ENDIF
1286     
1287           ENDIF
1288     
1289     
1290           d_ac = sqrt((xc-xa)**2 + (yc-ya)**2 + (zc-za)**2)
1291           d_bc = sqrt((xc-xb)**2 + (yc-yb)**2 + (zc-zb)**2)
1292     
1293     
1294           IF(d_ac<TOL_STL.OR.d_bc<TOL_STL) THEN
1295              INTERSECT_FLAG = .FALSE.             ! Exclude corner intersection
1296           ENDIF
1297     
1298     
1299      1000 FORMAT(A,3(2X,G12.5))
1300     
1301     
1302           RETURN
1303     
1304           END SUBROUTINE INTERSECT_LINE_WITH_FACET
1305     
1306     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1307     !                                                                      C
1308     !  Module name: IS_POINT_INSIDE_FACET                                  C
1309     !  Purpose: Verifies that point P is inside facet                      C
1310     !                                                                      C
1311     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1312     !  Reviewer:                                          Date:            C
1313     !                                                                      C
1314     !  Revision Number #                                  Date: ##-###-##  C
1315     !  Author: #                                                           C
1316     !  Purpose: #                                                          C
1317     !                                                                      C
1318     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1319       SUBROUTINE IS_POINT_INSIDE_FACET(Px,Py,Pz,FACET,INSIDE_FACET)
1320     
1321           USE param
1322           USE param1
1323           USE parallel
1324           USE constant
1325           USE run
1326           USE toleranc
1327           USE geometry
1328           USE indices
1329           USE compar
1330           USE sendrecv
1331           USE quadric
1332           USE STL
1333     
1334           IMPLICIT NONE
1335           INTEGER :: FACET,VV
1336           DOUBLE PRECISION :: NFx,NFy,NFz
1337           DOUBLE PRECISION :: Px,Py,Pz,V0x,V0y,V0z,V1x,V1y,V1z,V2x,V2y,V2z
1338           DOUBLE PRECISION :: dot00,dot01,dot02,dot11,dot12,dot_check
1339           DOUBLE PRECISION :: Inv_denom
1340           DOUBLE PRECISION :: u,v
1341           LOGICAL :: U_POSITIVE,V_POSITIVE,UPVL1,INSIDE_FACET
1342     
1343           DOUBLE PRECISION :: Vx,Vy,Vz
1344           DOUBLE PRECISION :: D(3),MINVAL_D,DH
1345     
1346     !======================================================================
1347     !  If point P is very close to one of the Facet vertices,
1348     !  consider that P belongs to facet, and return.
1349     !======================================================================
1350              DO VV = 1,3
1351                 Vx = VERTEX(VV,1,FACET)
1352                 Vy = VERTEX(VV,2,FACET)
1353                 Vz = VERTEX(VV,3,FACET)
1354                 D(VV) = sqrt((Px - Vx)**2 + (Py - Vy)**2 + (Pz - Vz)**2 )
1355              ENDDO
1356     
1357              MINVAL_D = MINVAL(D)
1358     
1359              IF(MINVAL_D < TOL_STL) THEN
1360                 INSIDE_FACET = .TRUE.
1361                 RETURN
1362              ENDIF
1363     
1364     !======================================================================
1365     !  Facet normal vector (normalized)
1366     !======================================================================
1367           NFx = NORM_FACE(1,FACET)
1368           NFy = NORM_FACE(2,FACET)
1369           NFz = NORM_FACE(3,FACET)
1370     
1371     !======================================================================
1372     ! This subroutine verifies that point P is inside triangular facet
1373     ! Facet is parametrized with (u,v), a point is inside triange if:
1374     ! u   >= 0  , and
1375     ! v   >= 0  , and
1376     ! u+v <= 1
1377     !======================================================================
1378     
1379           V0x = VERTEX(2,1,FACET) - VERTEX(1,1,FACET)
1380           V0y = VERTEX(2,2,FACET) - VERTEX(1,2,FACET)
1381           V0z = VERTEX(2,3,FACET) - VERTEX(1,3,FACET)
1382     
1383           V1x = VERTEX(3,1,FACET) - VERTEX(1,1,FACET)
1384           V1y = VERTEX(3,2,FACET) - VERTEX(1,2,FACET)
1385           V1z = VERTEX(3,3,FACET) - VERTEX(1,3,FACET)
1386     
1387           V2x = Px - VERTEX(1,1,FACET)
1388           V2y = Py - VERTEX(1,2,FACET)
1389           V2z = Pz - VERTEX(1,3,FACET)
1390     
1391           dot_check = NFx*V2x + NFy*V2y + NFz*V2z
1392     
1393     !      IF(dabs(dot_check)>TOL_STL_DP) THEN          ! reject points that do not
1394     !         INSIDE_FACET = .FALSE.                 ! belong to plane containing facet
1395     !         RETURN
1396     !      ENDIF
1397     
1398           DH = NFx * V2x + NFy * V2y + NFz * V2z
1399     
1400           IF(dabs(DH)>TOL_STL_DP) THEN          ! reject points that do not
1401              INSIDE_FACET = .FALSE.                 ! belong to plane containing facet
1402              RETURN
1403           ENDIF
1404     
1405           dot00 = V0x*V0x + V0y*V0y + V0z*V0z
1406           dot01 = V0x*V1x + V0y*V1y + V0z*V1z
1407           dot02 = V0x*V2x + V0y*V2y + V0z*V2z
1408           dot11 = V1x*V1x + V1y*V1y + V1z*V1z
1409           dot12 = V1x*V2x + V1y*V2y + V1z*V2z
1410     
1411           Inv_denom = ONE / (dot00*dot11 - dot01*dot01)
1412     
1413           u = (dot11*dot02 - dot01*dot12) * Inv_denom
1414           v = (dot00*dot12 - dot01*dot02) * Inv_denom
1415     
1416           U_POSITIVE = (u>=-TOL_STL)
1417           V_POSITIVE = (v>=-TOL_STL)
1418           UPVL1 = ((u+v)<=ONE+TOL_STL)
1419     
1420           INSIDE_FACET = (U_POSITIVE.AND.V_POSITIVE.AND.UPVL1)
1421     
1422           RETURN
1423     
1424           END SUBROUTINE IS_POINT_INSIDE_FACET
1425