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