File: /nfs/home/0/users/jenkins/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)
108     
109           IF(MyPE == PE_IO) WRITE(*,2000)'MSH file opened. Starting reading data...'
110     
111           OPEN(UNIT=444, FILE='checkgeometry.stl')
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')
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           INTEGER ::FILE_UNIT,I
589           DO I = 1,N_SKIP
590              READ(FILE_UNIT,*)
591           END DO
592           RETURN
593           END SUBROUTINE SKIP
594     
595           SUBROUTINE TEXT_HEX2INT(STRING,INT)
596     
597           CHARACTER(LEN=10)  :: INTERNAL
598           CHARACTER(LEN=*)   :: STRING
599           CHARACTER(LEN=10)  :: CLEANED_STRING
600           INTEGER :: INT
601     
602           CLEANED_STRING = STRING
603           IF(STRING(1:1)=='(') CLEANED_STRING = STRING(2:)
604     
605           WRITE(INTERNAL,fmt='(A10)') TRIM(CLEANED_STRING)
606           READ(INTERNAL,FMT='(Z10)') INT
607     
608           RETURN
609     
610           END SUBROUTINE TEXT_HEX2INT
611     
612           SUBROUTINE TEXT_DEC2INT(STRING,INT)
613     
614           CHARACTER(LEN=10)  :: INTERNAL
615           CHARACTER(LEN=*)   :: STRING
616           CHARACTER(LEN=10)  :: CLEANED_STRING
617           INTEGER :: INT
618     
619           CLEANED_STRING = STRING
620           IF(STRING(1:1)=='(') CLEANED_STRING = STRING(2:)
621     
622           WRITE(INTERNAL,fmt='(A10)') TRIM(CLEANED_STRING)
623           READ(INTERNAL,FMT='(I10)') INT
624     
625           RETURN
626     
627           END SUBROUTINE TEXT_DEC2INT
628     
629     
630     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
631     !                                                                      C
632     !  Module name: GET_STL_DATA                                           C
633     !  Purpose: reads face verticesa and normal vectors from an STL file   C
634     !                                                                      C
635     !  Author: Jeff Dietiker                              Date: 30-JAN-09  C
636     !  Reviewer:                                          Date: **-***-**  C
637     !                                                                      C
638     !  Revision Number:                                                    C
639     !  Purpose:                                                            C
640     !  Author:                                            Date: dd-mmm-yy  C
641     !  Reviewer:                                          Date: dd-mmm-yy  C
642     !                                                                      C
643     !  Literature/Document References:                                     C
644     !                                                                      C
645     !  Variables referenced:                                               C
646     !                                                                      C
647     !  Variables modified:                                                 C
648     !                                                                      C
649     !  Local variables:                                                    C
650     !                                                                      C
651     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
652     !
653           SUBROUTINE GET_STL_DATA
654     
655     !-----------------------------------------------
656     !   M o d u l e s
657     !-----------------------------------------------
658     
659           USE param
660           USE param1
661           USE physprop
662           USE fldvar
663           USE run
664           USE scalars
665           USE funits
666           USE rxns
667           USE compar
668           USE mpi_utility
669           USE progress_bar
670           USE stl
671           USE vtk
672           USE quadric
673           USE constant
674           USE bc
675           IMPLICIT NONE
676     
677           INTEGER :: N,IGNORED_FACETS
678           LOGICAL :: PRESENT,KEEP_READING,IGNORE_CURRENT_FACET
679           DOUBLE PRECISION ::v1x,v1y,v1z
680           DOUBLE PRECISION ::v2x,v2y,v2z
681           DOUBLE PRECISION ::v3x,v3y,v3z
682           DOUBLE PRECISION ::x12,y12,z12,x13,y13,z13,dp,d12,d13
683           DOUBLE PRECISION ::cos_angle,cos_small_angle
684           DOUBLE PRECISION ::n1,n2,n3,NORM
685           DOUBLE PRECISION ::ABSTRANS
686           CHARACTER(LEN=32) ::TEST_CHAR,BUFF_CHAR
687           CHARACTER(LEN=32) ::geometryfile(0:DIMENSION_BC)
688     
689           INTEGER :: BCV,NUMBER_OF_GEOMETRY_FILES,NUMBER_OF_BC_PATCHES
690           INTEGER :: BC_PATCH(0:DIMENSION_BC),L2,NF1,NF2
691           LOGICAL :: BC_PATCH_FOUND_IN_STL(DIMENSION_BC)
692     
693           CHARACTER(LEN=100) :: FNAME
694           integer :: stl_unit, nf
695     
696           BC_PATCH_FOUND_IN_STL = .FALSE.
697     
698           geometryfile(0) = 'geometry.stl'
699           NUMBER_OF_GEOMETRY_FILES = 0
700           NUMBER_OF_BC_PATCHES     = 0
701     
702     !     DETERMINE WHICH BOUNDARY CONDITIONS NEED STL FILE
703     
704           IF(MyPE == PE_IO) THEN
705              WRITE(*,100)'From Get_Stl_Data: Analysing BCs in mfix.dat...'
706              WRITE(*,120)'BC_ID','BC_TYPE'
707           ENDIF
708           DO BCV = 1, DIMENSION_BC
709     
710              IF(BC_TYPE(BCV)(1:2) == 'CG') THEN
711     
712                 NUMBER_OF_GEOMETRY_FILES = NUMBER_OF_GEOMETRY_FILES + 1
713                 NUMBER_OF_BC_PATCHES     = NUMBER_OF_BC_PATCHES + 1
714     
715     
716                 BC_PATCH(NUMBER_OF_GEOMETRY_FILES) = BCV
717                 WRITE(geometryfile(NUMBER_OF_GEOMETRY_FILES),200) 'geometry_',BCV
718     
719                 IF(MyPE == PE_IO) WRITE(*,130)BCV,BC_TYPE(BCV)
720     
721              ENDIF
722           ENDDO
723     
724           IF(MyPE == PE_IO) WRITE(*,110)'Number of CG BCs in mfix.dat = ',NUMBER_OF_BC_PATCHES
725     
726     
727     
728     
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     
738     ! VERIFY THAT THERE IS AT LEAST ONE STL FILE TO READ
739     
740           IF(NUMBER_OF_BC_PATCHES==0) THEN
741              IF(MyPE == PE_IO) THEN
742                 WRITE(*,100) 'ERROR: NO CARTESIAN GRID BOUNDARY CONDITION SPECIFIED.'
743                 WRITE(*,100) 'AT LEAST ONE BC_TYPE MUST START WITH CG (FOR EXAMPLE CG_NSW)'
744                 WRITE(*,100) 'RUN ABORTED'
745                 CALL MFIX_EXIT(MYPE)
746              ENDIF
747     
748           ELSEIF(NUMBER_OF_BC_PATCHES==1) THEN
749     
750              INQUIRE(FILE='geometry.stl',EXIST=PRESENT)
751              IF(PRESENT) THEN
752     
753     !            IF(STL_BC_ID==BC_PATCH(1)) THEN
754                    IF(MyPE == PE_IO) THEN
755                       WRITE(*,110) 'The file geometry.stl exists and the following  BC patch was found: ', BC_PATCH(1)
756                    ENDIF
757     
758                    NF1 = 1
759                    NF2 = 1
760                    geometryfile(NF1)='geometry.stl'
761     
762     !            ENDIF
763              ENDIF
764     
765           ELSE  ! More than one CG BC type
766              INQUIRE(FILE='geometry.stl',EXIST=PRESENT)
767              IF(PRESENT) THEN
768                 IF(MyPE == PE_IO) THEN
769                    WRITE(*,100) 'The file geometry.stl exists and several CG BC types are defined.'
770                    WRITE(*,100) 'All BC patches will be read from geometry.stl.'
771                 ENDIF
772     
773                 geometryfile(0)='geometry.stl'
774                 NF1 = 0
775                 NF2 = 0                       ! This invokes a special treatment
776                                               ! to reading multiple solids
777     
778              ELSE
779     
780                 NF1 = 1
781                 NF2 = NUMBER_OF_GEOMETRY_FILES
782                 IF(MyPE == PE_IO) THEN
783                    WRITE(*,100) 'The file geometry.stl does not exist and several CG_BC types are defined.'
784                    WRITE(*,100) 'Each BC patch will be read from geometry_BCID.stl.'
785                    WRITE(*,100) 'where BCID is 4-paded integer'
786                 ENDIF
787     
788              ENDIF
789     
790     
791     
792           ENDIF
793     
794     
795     ! START READING EACH STL FILE, ONE FOR EACH BC_TYPE
796     
797           N_FACETS = 0
798           IGNORED_FACETS = 0
799     
800     
801           DO N = NF1, NF2
802     
803              IF(MyPE == PE_IO) WRITE(*,2000) 'Reading geometry from '//TRIM(geometryfile(N))//' ...'
804     
805              INQUIRE(FILE=TRIM(geometryfile(N)),EXIST=PRESENT)
806              IF(.NOT.PRESENT) THEN
807                 IF(MyPE == PE_IO) THEN
808                    WRITE(*,"('(PE ',I3,'): input data file, ',A11,' is missing: run aborted')") &
809                    myPE,TRIM(geometryfile(N))
810                 ENDIF
811                 CALL MFIX_EXIT(MYPE)
812              ENDIF
813     
814     
815     !
816     !     OPEN geometry.stl ASCII FILE
817     !
818              OPEN(UNIT=333, FILE=TRIM(geometryfile(N)), STATUS='OLD', ERR=910)
819     
820              IF(MyPE == PE_IO) WRITE(*,2000)'STL file opened. Starting reading data...'
821     
822              KEEP_READING = .TRUE.
823     
824              DO WHILE(KEEP_READING)
825     
826                 READ(333,*,ERR=920,END=930) TEST_CHAR
827     !            print *,'TEST_CHAR=',TEST_CHAR
828                 IF(TRIM(TEST_CHAR) == 'solid'.AND.N==0) THEN
829     
830                    BACKSPACE(333)
831     
832                    READ(333,*,ERR=920,END=930) BUFF_CHAR,BUFF_CHAR
833     
834                    L2=LEN(TRIM(BUFF_CHAR))-3
835     
836                    READ(BUFF_CHAR(L2:L2+4),fmt='(I4.4)') BC_PATCH(N)
837     
838                    IF(MyPE == PE_IO)  WRITE(*,140) 'Found BC patch #', BC_PATCH(N), ' in STL file.'
839     
840                    IF(BC_TYPE(BC_PATCH(N))(1:2)/='CG') THEN
841                       IF(MyPE == PE_IO)  THEN
842                          WRITE(*,110) 'This BC patch does not mach a CG BC in mfix.dat:',BC_PATCH(N)
843                          WRITE(*,100)'Please Correct mfix.dat and/or stl file amd try again'
844                          CALL MFIX_EXIT(myPE)
845                       ENDIF
846                    ENDIF
847     
848                    BC_PATCH_FOUND_IN_STL(BC_PATCH(N)) = .TRUE.
849     
850                 ELSEIF(TRIM(TEST_CHAR) == 'facet') THEN
851     
852                    BACKSPACE(333)
853                    IGNORE_CURRENT_FACET = .FALSE.
854     
855                    READ(333,*,ERR=920,END=930) BUFF_CHAR,BUFF_CHAR,N1,N2,N3  ! Read unit normal vector
856                    READ(333,*,ERR=920,END=930)
857                    READ(333,*,ERR=920,END=930) BUFF_CHAR, V1x,V1y,V1z
858                    READ(333,*,ERR=920,END=930) BUFF_CHAR, V2x,V2y,V2z
859                    READ(333,*,ERR=920,END=930) BUFF_CHAR, V3x,V3y,V3z
860     
861                    N1 = N1 * OUT_STL_VALUE  ! Reverse unit vector if needed (this will switch fluid and blocked cells)
862                    N2 = N2 * OUT_STL_VALUE
863                    N3 = N3 * OUT_STL_VALUE
864     
865                    x12 = V2x - V1x
866                    y12 = V2y - V1y
867                    z12 = V2z - V1z
868     
869                    x13 = V3x - V1x
870                    y13 = V3y - V1y
871                    z13 = V3z - V1z
872     
873     
874                    dp  = x12*x13 + y12*y13 + z12*z13
875                    d12 = sqrt(x12**2+y12**2+z12**2)
876                    d13 = sqrt(x13**2+y13**2+z13**2)
877     
878                    IF((d12*d13)>TOL_STL) THEN
879                       cos_angle = dp/(d12*d13)
880                    ELSE
881                       cos_angle = ONE
882                    ENDIF
883     
884                    cos_small_angle = dcos(STL_SMALL_ANGLE / 180.0 * PI)
885     
886                    IF(DABS(cos_angle)>cos_small_angle) THEN
887                       IGNORE_CURRENT_FACET = .TRUE.    ! Ignore small facets
888                    ENDIF
889     
890                    NORM = sqrt(N1**2+N2**2+N3**2)
891     
892                    IF(NORM>TOL_STL) THEN
893                       N1 = N1 /NORM
894                       N2 = N2 /NORM
895                       N3 = N3 /NORM
896                    ELSE
897                       IGNORE_CURRENT_FACET = .TRUE.  ! Ignore facets with zero normal vector
898                    ENDIF
899     
900     
901                    IF(IGNORE_CURRENT_FACET) THEN
902                       IGNORED_FACETS = IGNORED_FACETS + 1
903                    ELSE                                                      ! Save vertex coordinates for valid facets
904                                                                              ! and performs translation
905                       N_FACETS = N_FACETS + 1
906     
907                       NORM_FACE(1,N_FACETS) = N1
908                       NORM_FACE(2,N_FACETS) = N2
909                       NORM_FACE(3,N_FACETS) = N3
910     
911                       VERTEX(1,1,N_FACETS) = SCALE_STL*V1x + TX_STL
912                       VERTEX(1,2,N_FACETS) = SCALE_STL*V1y + TY_STL
913                       VERTEX(1,3,N_FACETS) = SCALE_STL*V1z + TZ_STL
914     
915                       VERTEX(2,1,N_FACETS) = SCALE_STL*V2x + TX_STL
916                       VERTEX(2,2,N_FACETS) = SCALE_STL*V2y + TY_STL
917                       VERTEX(2,3,N_FACETS) = SCALE_STL*V2z + TZ_STL
918     
919                       VERTEX(3,1,N_FACETS) = SCALE_STL*V3x + TX_STL
920                       VERTEX(3,2,N_FACETS) = SCALE_STL*V3y + TY_STL
921                       VERTEX(3,3,N_FACETS) = SCALE_STL*V3z + TZ_STL
922     
923                       BC_ID_STL_FACE(N_FACETS) = BC_PATCH(N)
924     
925     
926                    ENDIF
927     
928     !            ELSEIF(TRIM(TEST_CHAR) == 'endsolid') THEN
929     
930     !               KEEP_READING = .FALSE.
931     
932                 ENDIF
933     
934              ENDDO
935     
936     
937     930      IF(MyPE == PE_IO)  WRITE(*,100) 'Done reading file.'
938     
939              CLOSE(333)
940     
941     
942     
943     
944           ENDDO
945     
946     
947           IF(myPE==0.AND.NF2==0) THEN
948              DO N = 1,NUMBER_OF_BC_PATCHES
949     
950                 IF(.NOT.BC_PATCH_FOUND_IN_STL(BC_PATCH(N))) THEN
951                    WRITE (*, 140)'Did not find BC patch: ',BC_PATCH(N) , ' in stl file'
952                    WRITE(*,100)'Please correct mfix.dat and/or stl file amd try again'
953                    CALL MFIX_EXIT(myPE)
954                 ENDIF
955     
956              ENDDO
957           ENDIF
958     
959     
960     
961     
962     
963           XMIN_STL = MINVAL(VERTEX(:,1,1:N_FACETS))
964           XMAX_STL = MAXVAL(VERTEX(:,1,1:N_FACETS))
965           YMIN_STL = MINVAL(VERTEX(:,2,1:N_FACETS))
966           YMAX_STL = MAXVAL(VERTEX(:,2,1:N_FACETS))
967           ZMIN_STL = MINVAL(VERTEX(:,3,1:N_FACETS))
968           ZMAX_STL = MAXVAL(VERTEX(:,3,1:N_FACETS))
969     
970           IF(MyPE == PE_IO) THEN
971              WRITE(*,2000)'STL file(s) successfully read.'
972              WRITE(*,110)'Total number of facets read =',N_FACETS + IGNORED_FACETS
973              WRITE(*,110)'Number of valid facets      =',N_FACETS
974              WRITE(*,110)'Number of ignored facets    =',IGNORED_FACETS
975              WRITE(*,100)'Geometry range from stl file:'
976              IF(SCALE_STL/=ONE) THEN
977                 WRITE(*,5000)'Scaling factor:',SCALE_STL
978              ENDIF
979              ABSTRANS = dabs(TX_STL)+dabs(TY_STL)+dabs(TZ_STL)
980              IF(ABSTRANS>TOL_STL) THEN
981                 WRITE(*,3000)'Translation vector (X,Y,Z)=',TX_STL,TY_STL,TZ_STL
982              ENDIF
983              WRITE(*,4000)'x-range = ', XMIN_STL,XMAX_STL
984              WRITE(*,4000)'y-range = ', YMIN_STL,YMAX_STL
985              WRITE(*,4000)'z-range = ', ZMIN_STL,ZMAX_STL
986              WRITE(*,4000)''
987           ENDIF
988     
989           XMIN_STL = XMIN_STL - 10.0*TOL_STL
990           XMAX_STL = XMAX_STL + 10.0*TOL_STL
991           YMIN_STL = YMIN_STL - 10.0*TOL_STL
992           YMAX_STL = YMAX_STL + 10.0*TOL_STL
993           ZMIN_STL = ZMIN_STL - 10.0*TOL_STL
994           ZMAX_STL = ZMAX_STL + 10.0*TOL_STL
995     
996     
997     
998     !      IF(XMIN_STL<ZERO) XMIN_STL=ZERO
999     !      IF(XMAX_STL>XLENGTH) XMAX_STL=XLENGTH
1000     !      IF(YMIN_STL<ZERO) YMIN_STL=ZERO
1001     !      IF(YMAX_STL>YLENGTH) YMAX_STL=YLENGTH
1002     !      IF(ZMIN_STL<ZERO) ZMIN_STL=ZERO
1003     !      IF(ZMAX_STL>ZLENGTH) ZMAX_STL=ZLENGTH
1004     
1005           IF(mype.eq.pe_io.and..false.) then
1006              WRITE(fname,'(A,"_FACETS_READ", ".stl")') &
1007              TRIM(RUN_NAME)
1008              open(stl_unit, file = fname, form='formatted')
1009              write(stl_unit,*)'solid vcg'
1010     
1011     
1012              do nf = 1, n_facets
1013     
1014                 write(stl_unit,*) '   facet normal ', NORM_FACE(:,NF)
1015                 write(stl_unit,*) '      outer loop'
1016                 write(stl_unit,*) '         vertex ', VERTEX(1,1:3,NF)
1017                 write(stl_unit,*) '         vertex ', VERTEX(2,1:3,NF)
1018                 write(stl_unit,*) '         vertex ', VERTEX(3,1:3,NF)
1019                 write(stl_unit,*) '      endloop'
1020                 write(stl_unit,*)'   endfacet'
1021              enddo
1022     
1023     
1024              write(stl_unit,*)'endsolid vcg'
1025              close(stl_unit, status = 'keep')
1026              IF(MyPE == PE_IO) THEN
1027                 WRITE(*,100) 'The file FACETS_READ.stl was sucessfully written.'
1028                 WRITE(*,100) 'and is provided for convenience (it is not used).'
1029              ENDIF
1030     
1031           endif
1032     
1033           RETURN
1034     
1035     !======================================================================
1036     !     HERE IF AN ERROR OCCURED OPENNING/READING THE FILE
1037     !======================================================================
1038     !
1039      910  CONTINUE
1040           WRITE (*, 1500)
1041           CALL MFIX_EXIT(myPE)
1042      920  CONTINUE
1043           WRITE (*, 1600)
1044           CALL MFIX_EXIT(myPE)
1045     ! 930  CONTINUE
1046     !     WRITE (*, 1700)
1047     !     CALL MFIX_EXIT(myPE)
1048     !
1049      1500 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
1050           'Unable to open stl file',/1X,70('*')/)
1051      1600 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
1052           'Error while reading stl file',/1X,70('*')/)
1053      1700 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
1054           'End of file reached while reading stl file',/1X,70('*')/)
1055      2000 FORMAT(1X,A)
1056      2010 FORMAT(1X,A,I4,A)
1057      3000 FORMAT(1X,A,'(',F10.4,';',F10.4,';',F10.4,')')
1058      4000 FORMAT(1X,A,F10.4,' to ',F10.4)
1059      5000 FORMAT(1X,A,F10.4)
1060           END SUBROUTINE GET_STL_DATA
1061     
1062     
1063     
1064     
1065     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1066     !                                                                      C
1067     !  Module name: EVAL_STL_FCT_AT                                        C
1068     !  Purpose: Assigns a value to f_stl at any NODE of cell IJK           C
1069     !                                                                      C
1070     !  Author: Jeff Dietiker                              Date: 30-JAN-09  C
1071     !  Reviewer:                                          Date: **-***-**  C
1072     !                                                                      C
1073     !  Revision Number:                                                    C
1074     !  Purpose:                                                            C
1075     !  Author:                                            Date: dd-mmm-yy  C
1076     !  Reviewer:                                          Date: dd-mmm-yy  C
1077     !                                                                      C
1078     !  Literature/Document References:                                     C
1079     !                                                                      C
1080     !  Variables referenced:                                               C
1081     !                                                                      C
1082     !  Variables modified:                                                 C
1083     !                                                                      C
1084     !  Local variables:                                                    C
1085     !                                                                      C
1086     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1087     !
1088           SUBROUTINE EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,NODE,f_stl,CLIP_FLAG,BCID)
1089     
1090     !-----------------------------------------------
1091     !   M o d u l e s
1092     !-----------------------------------------------
1093     
1094           USE param
1095           USE param1
1096           USE physprop
1097           USE fldvar
1098           USE run
1099           USE scalars
1100           USE funits
1101           USE rxns
1102           USE compar
1103           USE mpi_utility
1104           USE progress_bar
1105     !      USE quadric
1106           USE cutcell
1107           USE stl
1108           USE functions
1109           IMPLICIT NONE
1110     
1111           INTEGER :: BCID
1112           LOGICAL :: CLIP_FLAG
1113           DOUBLE PRECISION :: f_stl
1114     
1115           CHARACTER (LEN=*) :: TYPE_OF_CELL
1116           INTEGER :: IJK,IJKC,NODE
1117     
1118           IF(NODE==15) print*,'Warning: eval stl at node 15'
1119     
1120           IF(N_FACETS < 1) RETURN
1121     
1122           f_stl = UNDEFINED
1123     
1124           IJKC = IJK_OF_NODE(NODE)
1125     
1126           IF(IJKC<1.OR.IJKC>DIMENSION_3) THEN
1127     
1128     !         print*,'myPE        =',myPE
1129     !         print*,'DIMENSION_3 =',DIMENSION_3
1130     !         print*,'IJKC        =',IJKC
1131     !         print*,'IJK         =',IJK
1132     !         print*,'I,J,K       =',I_OF(IJK),J_OF(IJK),K_OF(IJK)
1133     !         print*,'NODE        =',NODE
1134     
1135     !     Leave f_stl as undefined if we are out of bounds
1136     !     This can happen at I=1 and NODE=4 for example
1137              RETURN
1138     
1139           ENDIF
1140     
1141     !!      IF(NODE/=15.AND.NODE/=0) THEN
1142     !!         IF(SNAP(IJKC)) THEN
1143     !!            f_stl = ZERO
1144     !!!            print*,'stlfc snapped=',IJKC,NODE,F_AT(IJKC)
1145     !!            RETURN
1146     !!         ENDIF
1147     !!      ENDIF
1148     
1149     
1150           f_stl = F_AT(IJKC)
1151     
1152           RETURN
1153     
1154     
1155           CLIP_FLAG = .TRUE.
1156     
1157     
1158           RETURN
1159     
1160     
1161           END SUBROUTINE EVAL_STL_FCT_AT
1162     
1163     
1164     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1165     !                                                                      C
1166     !  Module name: intersect_line_with_facet                              C
1167     !  Purpose: Finds the intersection between a facet                     C
1168     !           and the line (xa,ya,za) and (xb,yb,zb).                    C
1169     !                                                                      C
1170     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1171     !  Reviewer:                                          Date:            C
1172     !                                                                      C
1173     !  Revision Number #                                  Date: ##-###-##  C
1174     !  Author: #                                                           C
1175     !  Purpose: #                                                          C
1176     !                                                                      C
1177     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1178       SUBROUTINE INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,FACET,INTERSECT_FLAG,xc,yc,zc)
1179     
1180           USE param
1181           USE param1
1182           USE parallel
1183           USE constant
1184           USE run
1185           USE toleranc
1186           USE geometry
1187           USE indices
1188           USE compar
1189           USE sendrecv
1190           USE quadric
1191           USE STL
1192     
1193           IMPLICIT NONE
1194           DOUBLE PRECISION:: xa,ya,za,xb,yb,zb,xc,yc,zc
1195           INTEGER :: FACET
1196           DOUBLE PRECISION:: NFx,NFy,NFz,nabx,naby,nabz
1197           DOUBLE PRECISION :: dot_denom,dot_num
1198           DOUBLE PRECISION :: VP1Ax,VP1Ay,VP1Az
1199           DOUBLE PRECISION :: Px,Py,Pz
1200           DOUBLE PRECISION :: t
1201           DOUBLE PRECISION :: d_ac,d_bc
1202           LOGICAL :: INSIDE_FACET,INTERSECT_FLAG
1203     
1204     
1205     !======================================================================
1206     ! This subroutine tries to find the intersection of a line AB with one
1207     ! of the facets defined in the stl file
1208     !
1209     !
1210     ! 1) Verify that intersection is possible. Only facets and line AB that are
1211     !    not parallel are acceptable candidates
1212     ! 2) Find intersection point P of line AB with plane containing the facet.
1213     !    Valid intersection point must be between A and B to continue.
1214     ! 3) Verify that point P is inside triangular facet
1215     !======================================================================
1216     
1217           INTERSECT_FLAG = .FALSE.
1218     
1219     !======================================================================
1220     !  Facet normal vector (normalized)
1221     !======================================================================
1222           NFx = NORM_FACE(1,FACET)
1223           NFy = NORM_FACE(2,FACET)
1224           NFz = NORM_FACE(3,FACET)
1225     
1226     !======================================================================
1227     !  AB vector (NOT normalized)
1228     !======================================================================
1229           nabx = xb - xa
1230           naby = yb - ya
1231           nabz = zb - za
1232     
1233           dot_denom = NFx*nabx + NFy*naby + NFz*nabz
1234     !======================================================================
1235     ! 1) Verify that intersection is possible. Only facets and line AB that are
1236     !    not parallel are acceptable candidates
1237     !======================================================================
1238           IF(dabs(dot_denom)<TOL_STL) THEN
1239     
1240              INTERSECT_FLAG = .FALSE.              ! No intersection (facet nornal
1241                                                    ! and AB are perpendicular
1242              RETURN
1243     
1244           ELSE
1245     
1246     !======================================================================
1247     ! 2) Find intersection point P of line AB with plane containing the facet
1248     !    Line AB is parametrized with parameter t (from t=0 at A to t=1 at B)
1249     !======================================================================
1250              VP1Ax = VERTEX(1,1,FACET) - xa
1251              VP1Ay = VERTEX(1,2,FACET) - ya
1252              VP1Az = VERTEX(1,3,FACET) - za
1253     
1254              dot_num = NFx*VP1Ax + NFy*VP1Ay + NFz*VP1Az
1255     
1256              t = dot_num / dot_denom
1257     
1258     !======================================================================
1259     ! 3) Verify that point P is inside triangular facet
1260     !    Facet is parametrized with (u,v), a point is inside triange if:
1261     !    u   >= 0  , and
1262     !    v   >= 0  , and
1263     !    u+v <= 1
1264     !======================================================================
1265              IF((t>=ZERO).AND.(t<=ONE)) THEN       ! Intersection between A and B
1266                                                    ! Now test if intersection point is inside triangle
1267     
1268                 Px = xa + t*nabx
1269                 Py = ya + t*naby
1270                 Pz = za + t*nabz
1271     
1272                 CALL IS_POINT_INSIDE_FACET(Px,Py,Pz,FACET,INSIDE_FACET)
1273     
1274                 IF(INSIDE_FACET) THEN
1275                    INTERSECT_FLAG = .TRUE.              ! Valid intersection
1276                    xc = Px                              ! point P is inside triangle)
1277                    yc = Py
1278                    zc = Pz
1279                 ELSE
1280                    INTERSECT_FLAG = .FALSE.             ! Invalid intersection
1281                    RETURN                               ! point P is outside triangle)
1282                 ENDIF
1283     
1284              ELSE
1285     
1286                 INTERSECT_FLAG = .FALSE.           ! No intersection between A and B
1287                 RETURN
1288     
1289              ENDIF
1290     
1291           ENDIF
1292     
1293     
1294           d_ac = sqrt((xc-xa)**2 + (yc-ya)**2 + (zc-za)**2)
1295           d_bc = sqrt((xc-xb)**2 + (yc-yb)**2 + (zc-zb)**2)
1296     
1297     
1298           IF(d_ac<TOL_STL.OR.d_bc<TOL_STL) THEN
1299              INTERSECT_FLAG = .FALSE.             ! Exclude corner intersection
1300           ENDIF
1301     
1302     
1303      1000 FORMAT(A,3(2X,G12.5))
1304     
1305     
1306           RETURN
1307     
1308           END SUBROUTINE INTERSECT_LINE_WITH_FACET
1309     
1310     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1311     !                                                                      C
1312     !  Module name: IS_POINT_INSIDE_FACET                                  C
1313     !  Purpose: Verifies that point P is inside facet                      C
1314     !                                                                      C
1315     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1316     !  Reviewer:                                          Date:            C
1317     !                                                                      C
1318     !  Revision Number #                                  Date: ##-###-##  C
1319     !  Author: #                                                           C
1320     !  Purpose: #                                                          C
1321     !                                                                      C
1322     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1323       SUBROUTINE IS_POINT_INSIDE_FACET(Px,Py,Pz,FACET,INSIDE_FACET)
1324     
1325           USE param
1326           USE param1
1327           USE parallel
1328           USE constant
1329           USE run
1330           USE toleranc
1331           USE geometry
1332           USE indices
1333           USE compar
1334           USE sendrecv
1335           USE quadric
1336           USE STL
1337     
1338           IMPLICIT NONE
1339           INTEGER :: FACET,VV
1340           DOUBLE PRECISION :: NFx,NFy,NFz
1341           DOUBLE PRECISION :: Px,Py,Pz,V0x,V0y,V0z,V1x,V1y,V1z,V2x,V2y,V2z
1342           DOUBLE PRECISION :: dot00,dot01,dot02,dot11,dot12,dot_check
1343           DOUBLE PRECISION :: Inv_denom
1344           DOUBLE PRECISION :: u,v
1345           LOGICAL :: U_POSITIVE,V_POSITIVE,UPVL1,INSIDE_FACET
1346     
1347     
1348           DOUBLE PRECISION :: Vx,Vy,Vz
1349           DOUBLE PRECISION :: D(3),MINVAL_D,DH
1350     
1351     !======================================================================
1352     !  If point P is very close to one of the Facet vertices,
1353     !  consider that P belongs to facet, and return.
1354     !======================================================================
1355              DO VV = 1,3
1356                 Vx = VERTEX(VV,1,FACET)
1357                 Vy = VERTEX(VV,2,FACET)
1358                 Vz = VERTEX(VV,3,FACET)
1359                 D(VV) = sqrt((Px - Vx)**2 + (Py - Vy)**2 + (Pz - Vz)**2 )
1360              ENDDO
1361     
1362              MINVAL_D = MINVAL(D)
1363     
1364              IF(MINVAL_D < TOL_STL) THEN
1365                 INSIDE_FACET = .TRUE.
1366                 RETURN
1367              ENDIF
1368     
1369     !======================================================================
1370     !  Facet normal vector (normalized)
1371     !======================================================================
1372           NFx = NORM_FACE(1,FACET)
1373           NFy = NORM_FACE(2,FACET)
1374           NFz = NORM_FACE(3,FACET)
1375     
1376     !======================================================================
1377     ! This subroutine verifies that point P is inside triangular facet
1378     ! Facet is parametrized with (u,v), a point is inside triange if:
1379     ! u   >= 0  , and
1380     ! v   >= 0  , and
1381     ! u+v <= 1
1382     !======================================================================
1383     
1384           V0x = VERTEX(2,1,FACET) - VERTEX(1,1,FACET)
1385           V0y = VERTEX(2,2,FACET) - VERTEX(1,2,FACET)
1386           V0z = VERTEX(2,3,FACET) - VERTEX(1,3,FACET)
1387     
1388           V1x = VERTEX(3,1,FACET) - VERTEX(1,1,FACET)
1389           V1y = VERTEX(3,2,FACET) - VERTEX(1,2,FACET)
1390           V1z = VERTEX(3,3,FACET) - VERTEX(1,3,FACET)
1391     
1392           V2x = Px - VERTEX(1,1,FACET)
1393           V2y = Py - VERTEX(1,2,FACET)
1394           V2z = Pz - VERTEX(1,3,FACET)
1395     
1396           dot_check = NFx*V2x + NFy*V2y + NFz*V2z
1397     
1398     !      IF(dabs(dot_check)>TOL_STL_DP) THEN          ! reject points that do not
1399     !         INSIDE_FACET = .FALSE.                 ! belong to plane containing facet
1400     !         RETURN
1401     !      ENDIF
1402     
1403           DH = NFx * V2x + NFy * V2y + NFz * V2z
1404     
1405           IF(dabs(DH)>TOL_STL_DP) THEN          ! reject points that do not
1406              INSIDE_FACET = .FALSE.                 ! belong to plane containing facet
1407              RETURN
1408           ENDIF
1409     
1410           dot00 = V0x*V0x + V0y*V0y + V0z*V0z
1411           dot01 = V0x*V1x + V0y*V1y + V0z*V1z
1412           dot02 = V0x*V2x + V0y*V2y + V0z*V2z
1413           dot11 = V1x*V1x + V1y*V1y + V1z*V1z
1414           dot12 = V1x*V2x + V1y*V2y + V1z*V2z
1415     
1416           Inv_denom = ONE / (dot00*dot11 - dot01*dot01)
1417     
1418           u = (dot11*dot02 - dot01*dot12) * Inv_denom
1419           v = (dot00*dot12 - dot01*dot02) * Inv_denom
1420     
1421           U_POSITIVE = (u>=-TOL_STL)
1422           V_POSITIVE = (v>=-TOL_STL)
1423           UPVL1 = ((u+v)<=ONE+TOL_STL)
1424     
1425           INSIDE_FACET = (U_POSITIVE.AND.V_POSITIVE.AND.UPVL1)
1426     
1427     
1428           RETURN
1429     
1430           END SUBROUTINE IS_POINT_INSIDE_FACET
1431     
1432     
1433     
1434     
1435