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