File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/get_stl_data.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 SUBROUTINE GET_MSH_DATA
28
29
30
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
53 INTEGER,PARAMETER :: MAX_ZONES = 1000
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
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
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
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
152 ENDIF
153
154 ENDDO
155
156 CALL SKIP(333,3)
157
158
159
160 READ(333,*) (BUFF_CHAR(I),I=1,4)
161
162 CALL TEXT_HEX2INT(BUFF_CHAR(3),N_FACES)
163
164
165 = 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
224 (:,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
261
262
263
264 = 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
286 (:,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
312
313 = 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
335 (:,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
365
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
377 = 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
393
394 = 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
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
556
557
558 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 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
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653 SUBROUTINE GET_STL_DATA
654
655
656
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
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
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
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
763 ENDIF
764
765 ELSE
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
776
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
796
797 = 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
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
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
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
862 = 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.
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.
898 ENDIF
899
900
901 IF(IGNORE_CURRENT_FACET) THEN
902 IGNORED_FACETS = IGNORED_FACETS + 1
903 ELSE
904
905 = 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
929
930
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
999
1000
1001
1002
1003
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
1037
1038
1039 CONTINUE
1040 WRITE (*, 1500)
1041 CALL MFIX_EXIT(myPE)
1042 920 CONTINUE
1043 WRITE (*, 1600)
1044 CALL MFIX_EXIT(myPE)
1045
1046
1047
1048
1049 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
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088 SUBROUTINE EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,NODE,f_stl,CLIP_FLAG,BCID)
1089
1090
1091
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
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
1129
1130
1131
1132
1133
1134
1135
1136
1137 RETURN
1138
1139 ENDIF
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150 = 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
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
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
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217 = .FALSE.
1218
1219
1220
1221
1222 = NORM_FACE(1,FACET)
1223 NFy = NORM_FACE(2,FACET)
1224 NFz = NORM_FACE(3,FACET)
1225
1226
1227
1228
1229 = xb - xa
1230 naby = yb - ya
1231 nabz = zb - za
1232
1233 dot_denom = NFx*nabx + NFy*naby + NFz*nabz
1234
1235
1236
1237
1238 IF(dabs(dot_denom)<TOL_STL) THEN
1239
1240 INTERSECT_FLAG = .FALSE.
1241
1242 RETURN
1243
1244 ELSE
1245
1246
1247
1248
1249
1250 = 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
1260
1261
1262
1263
1264
1265 IF((t>=ZERO).AND.(t<=ONE)) THEN
1266
1267
1268 = 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.
1276 = Px
1277 = Py
1278 zc = Pz
1279 ELSE
1280 INTERSECT_FLAG = .FALSE.
1281 RETURN
1282 ENDIF
1283
1284 ELSE
1285
1286 INTERSECT_FLAG = .FALSE.
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.
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
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
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
1353
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
1371
1372 = NORM_FACE(1,FACET)
1373 NFy = NORM_FACE(2,FACET)
1374 NFz = NORM_FACE(3,FACET)
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384 = 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
1399
1400
1401
1402
1403 = NFx * V2x + NFy * V2y + NFz * V2z
1404
1405 IF(dabs(DH)>TOL_STL_DP) THEN
1406 = .FALSE.
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