File: N:\mfix\model\des\stl_functions_des_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module name: stl_functions_des                                      !
4     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
5     !                                                                      !
6     !  Purpose: This module containd routines for geometric interaction    !
7     !  required for STL files.                                             !
8     !                                                                      !
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
10           MODULE stl_functions_des
11     
12           IMPLICIT NONE
13     
14     ! Use this module only to define functions and subroutines.
15           CONTAINS
16     
17     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
18     !                                                                      !
19     !  Subroutine: ClosestPtPointTriangle                                  !
20     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
21     !                                                                      !
22     !  Purpose:                                                            !
23     !                                                                      !
24     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
25           Subroutine ClosestPtPointTriangle(pointp, points, closest_point)
26           USE param1, only: zero, one
27           USE discretelement, only: dimn
28     
29           IMPLICIT NONE
30     
31     ! points are the three nodes of the triangle
32     ! point p is the sphere center
33           double precision, intent(in), dimension(3,3) :: points
34           double precision, intent(in), dimension(dimn) :: pointp
35           double precision, intent(out), dimension(dimn) ::  closest_point
36     ! Local variables
37           double precision, dimension(3) :: pointa, pointb, pointc
38           double precision, dimension(dimn) :: ab, ac, ap, bp,cp
39           double precision :: d1, d2, d3, d4, vc, v, d5, d6, vb, w, va, denom
40     
41           pointa = points(1,:)
42           pointb = points(2,:)
43           pointc = points(3,:)
44     
45           ab = pointb - pointa
46           ac = pointc - pointa
47           ap = pointp - pointa
48           d1 = DOT_PRODUCT(ab, ap)
49           d2 = DOT_PRODUCT(ac, ap)
50     
51           IF(d1 <= Zero .AND. d2 <= zero) then
52              closest_point = pointa
53              return
54           end if
55     
56     ! Check if P in vertex region outside B
57           bp = pointp - pointb
58           d3 = DOT_PRODUCT(ab, bp);
59           d4 = DOT_PRODUCT(ac, bp);
60           if (d3 >= zero .and. d4 <= d3) then
61              closest_point = pointb
62              return
63           endif
64     
65     ! Check if P in edge region of AB, if so return projection of P onto AB
66           vc = d1*d4 - d3*d2;
67           if (vc <= zero .and. d1 >= zero .and. d3 <= zero) then
68              v = d1 / (d1 - d3);
69              closest_point =  pointa + v * ab;
70              return
71           end if
72     
73     ! Check if P in vertex region outside C
74           cp = pointp - pointc
75           d5 = DOT_PRODUCT(ab, cp)
76           d6 = DOT_PRODUCT(ac, cp)
77           if (d6 >= zero .and. d5 <= d6) then
78              closest_point  = pointc
79              return
80           endif
81     
82     ! Check if P in edge region of AC, if so return projection of P onto AC
83           vb = d5*d2 - d1*d6
84     
85           if (vb <= zero .and. d2 >= zero .and. d6 <= zero) then
86              w = d2 / (d2 - d6)
87              closest_point = pointa + w * ac
88              return
89           end if
90     
91     ! Check if P in edge region of BC, if so return projection of P onto BC
92           va = d3*d6 - d5*d4
93           if (va <= zero .and.(d4 - d3) >= zero .and. (d5 - d6) >= zero) then
94              w = (d4 - d3) / ((d4 - d3) + (d5 - d6))
95              closest_point = pointb + w * (pointc - pointb)
96              return
97           end if
98     
99     
100     ! P inside face region. Compute Q through its barycentric coordinates (u,v,w)
101           denom = one / (va + vb + vc)
102           v = vb * denom
103           w = vc * denom
104           closest_point = pointa + ab * v + ac * w;
105           return
106           end Subroutine ClosestPtPointTriangle
107     
108     
109     
110     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
111     !                                                                      !
112     !  Subroutine: intersectLnPlane                                        !
113     !  Author: Rahul Garg                                 Date: 24-Oct-13  !
114     !                                                                      !
115     !  Purpose:                                                            !
116     !                                                                      !
117     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
118           Subroutine intersectLnPlane(ref_line, dir_line, ref_plane,       &
119              norm_plane, line_param)
120     
121           USE discretelement, only: dimn
122           USE param1, only: zero
123     
124           IMPLICIT NONE
125     
126     ! Reference point and direction of the line
127           DOUBLE PRECISION, INTENT(IN) :: REF_LINE(3),  DIR_LINE(3)
128     ! reference point and normal of the plane
129           DOUBLE PRECISION, INTENT(IN) :: REF_PLANE(3), NORM_PLANE(3)
130     
131     ! line is parameterized as p = p_ref + t * dir_line, t is line_param
132           double precision, intent(out) :: line_param
133     
134           !local vars
135           double precision :: denom
136     
137           denom = DOT_PRODUCT(dir_line, norm_plane)
138     
139           if(denom*denom.gt.zero) then
140              line_param = DOT_PRODUCT(ref_plane-ref_line, norm_plane)
141              line_param = line_param/denom
142           endif
143     
144           return
145           end subroutine intersectLnPlane
146     
147     !......................................................................!
148     ! Subroutine TRI_BOX_OVERLAP                                           !
149     ! Author: J.Musser                                   Date: 10-22-2015  !
150     !                                                                      !
151     ! Purpose: Determine if a box (DES grid cell) intersects the triangle  !
152     !    (SLT). Note that the DES grid size is slightly increased to       !
153     !    capture STLs near the boarder of the cell. Otherwise, some        !
154     !    collisions could ve over looked.                                  !
155     !                                                                      !
156     ! Author: Tomas Akenine-Moller                   Accessed: 10-22-2015  !
157     ! REF: http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/    !
158     !         code/tribox2.txt                                             !
159     !......................................................................!
160           SUBROUTINE TRI_BOX_OVERLAP(pCENTER, pHALFSIZE, pVERTS, pOVERLAP)
161     
162           IMPLICIT NONE
163     
164           DOUBLE PRECISION, INTENT(IN) :: pCENTER(3), pHALFSIZE(3)
165           DOUBLE PRECISION, INTENT(IN) :: pVERTS(3,3)
166           LOGICAL, INTENT(OUT) :: pOVERLAP
167     
168           DOUBLE PRECISION :: v0(3), v1(3), v2(3)
169           DOUBLE PRECISION :: fex, fey, fez
170           DOUBLE PRECISION :: normal(3), e0(3), e1(3), e2(3)
171     
172           pOVERLAP = .FALSE.
173     
174           v0 = pVERTS(1,:) - pCENTER
175           v1 = pVERTS(2,:) - pCENTER
176           v2 = pVERTS(3,:) - pCENTER
177     
178           e0 = v1-v0
179           e1 = v2-v1
180           e2 = v0-v2
181     
182           fex = abs(e0(1))
183           fey = abs(e0(2))
184           fez = abs(e0(3))
185     
186           if(ATEST_X01(e0(3),e0(2),fez,fey)) return
187           if(ATEST_Y02(e0(3),e0(1),fez,fex)) return
188           if(ATEST_Z12(e0(2),e0(1),fey,fex)) return
189     
190           fex = abs(e1(1))
191           fey = abs(e1(2))
192           fez = abs(e1(3))
193     
194           if(ATEST_X01(e1(3),e1(2),fez,fey)) return
195           if(ATEST_Y02(e1(3),e1(1),fez,fex)) return
196           if(ATEST_Z0 (e1(2),e1(1),fey,fex)) return
197     
198           fex = abs(e2(1))
199           fey = abs(e2(2))
200           fez = abs(e2(3))
201     
202           if(ATEST_X2 (e2(3),e2(2),fez,fey)) return
203           if(ATEST_Y1 (e2(3),e2(1),fez,fex)) return
204           if(ATEST_Z12(e2(2),e2(1),fey,fex)) return
205     
206           if(findMin(v0(1),v1(1),v2(1)) > phalfsize(1) .OR. &
207              findMax(v0(1),v1(1),v2(1)) <-phalfsize(1)) return
208     
209           if(findMin(v0(2),v1(2),v2(2)) > phalfsize(2) .OR. &
210              findMax(v0(2),v1(2),v2(2)) <-phalfsize(2)) return
211     
212           if(findMin(v0(3),v1(3),v2(3)) > phalfsize(3) .OR. &
213              findMax(v0(3),v1(3),v2(3)) <-phalfsize(3)) return
214     
215     
216           normal(1) = e0(2)*e1(3)-e0(3)*e1(2)
217           normal(2) = e0(3)*e1(1)-e0(1)*e1(3)
218           normal(3) = e0(1)*e1(2)-e0(2)*e1(1)
219     
220           if(.NOT.planeBoxOverlap(normal,v0,phalfsize)) return
221     
222           pOVERLAP = .TRUE.
223     
224           RETURN
225     
226           CONTAINS
227     
228     !``````````````````````````````````````````````````````````````````````!
229     ! Function: planeBoxOverlap                                            !
230     ! Support function for TRI_BOX_OVERLAP                                 !
231     !``````````````````````````````````````````````````````````````````````!
232           LOGICAL FUNCTION planeBoxOverlap(norm, vert, maxbox)
233     
234           double precision :: norm(3), vert(3), maxbox(3)
235     
236           integer :: lc
237           double precision :: vmin(3), vmax(3), v
238     
239           do lc=1,3
240              v=vert(lc)
241              if(norm(lc) > 0.0d0) then
242                 vmin(lc) = -maxbox(lc) - v
243                 vmax(lc) =  maxbox(lc) - v
244              else
245                 vmin(lc) = maxbox(lc) - v
246                 vmax(lc) =-maxbox(lc) - v
247              endif
248           enddo
249     
250           if(dot_product(norm,vmin) > 0.0d0) then
251              planeBoxOverlap = .false.
252              return
253           elseif(dot_product(norm,vmax) >= 0.0d0) then
254              planeBoxOverlap = .true.
255              return
256           endif
257     
258           planeBoxOverlap = .false.
259           return
260     
261           RETURN
262           END FUNCTION planeBoxOverlap
263     
264     !``````````````````````````````````````````````````````````````````````!
265     ! Function: findMin                                                    !
266     ! Support function for TRI_BOX_OVERLAP                                 !
267     !``````````````````````````````````````````````````````````````````````!
268           DOUBLE PRECISION FUNCTION findMin(x0,x1,x2)
269     
270           double precision :: x0,x1,x2
271     
272           findMin = x0
273     
274           if(x1<findMin) findMin=x1
275           if(x2<findMin) findMin=x2
276     
277           RETURN
278           END FUNCTION findMin
279     
280     !``````````````````````````````````````````````````````````````````````!
281     ! Function: findMax                                                    !
282     ! Support function for TRI_BOX_OVERLAP                                 !
283     !``````````````````````````````````````````````````````````````````````!
284           DOUBLE PRECISION FUNCTION findMax(x0,x1,x2)
285     
286           double precision :: x0,x1,x2
287     
288           findMax = x0
289     
290           if(x1>findMax) findMax=x1
291           if(x2>findMax) findMax=x2
292     
293           RETURN
294           END FUNCTION findMax
295     
296     !``````````````````````````````````````````````````````````````````````!
297     ! Function: ATEST_X01                                                  !
298     ! Support function for TRI_BOX_OVERLAP                                 !
299     !``````````````````````````````````````````````````````````````````````!
300           LOGICAL FUNCTION ATEST_X01(a,b,fa,fb)
301     
302           double precision :: a, b, fa, fb
303           double precision :: lMin, lMax, p0, p2, rad
304     
305           p0 = a*v0(2) - b*v0(3)
306           p2 = a*v2(2) - b*v2(3)
307     
308           if(p0<p2) then; lMIN=p0; lMAX=p2
309           else; lMIN=p2; lMAX=p0; endif
310     
311           rad=fa*phalfsize(2) + fb*phalfsize(3)
312           ATEST_X01=(lmin>rad .OR. lmax<-rad)
313     
314           END FUNCTION ATEST_X01
315     
316     !``````````````````````````````````````````````````````````````````````!
317     ! Function: ATEST_X2                                                   !
318     ! Support function for TRI_BOX_OVERLAP                                 !
319     !``````````````````````````````````````````````````````````````````````!
320           LOGICAL FUNCTION ATEST_X2(a,b,fa,fb)
321     
322           double precision :: a, b, fa, fb
323           double precision :: lMin, lMax, p0, p1, rad
324     
325           p0 = a*v0(2) - b*v0(3)
326           p1 = a*v1(2) - b*v1(3)
327     
328           if(p0<p1) then; lMIN=p0; lMAX=p1
329           else; lMIN=p1; lMAX=p0; endif
330     
331           rad=fa*phalfsize(2) + fb*phalfsize(3)
332           ATEST_X2=(lmin>rad .OR. lmax<-rad)
333     
334           END FUNCTION ATEST_X2
335     
336     !``````````````````````````````````````````````````````````````````````!
337     ! Function: ATEST_Y02                                                  !
338     ! Support function for TRI_BOX_OVERLAP                                 !
339     !``````````````````````````````````````````````````````````````````````!
340           LOGICAL FUNCTION ATEST_Y02(a,b,fa,fb)
341     
342           double precision :: a, b, fa, fb
343           double precision :: lMin, lMax, p0, p2, rad
344     
345           p0 = -a*v0(1) + b*v0(3)
346           p2 = -a*v2(1) + b*v2(3)
347     
348           if(p0<p2) then; lMIN=p0; lMAX=p2
349           else; lMIN=p2; lMAX=p0; endif
350     
351           rad=fa*phalfsize(1) + fb*phalfsize(3)
352           ATEST_Y02=(lmin>rad .OR. lmax<-rad)
353     
354           END FUNCTION ATEST_Y02
355     
356     !``````````````````````````````````````````````````````````````````````!
357     ! Function: ATEST_Y1                                                   !
358     ! Support function for TRI_BOX_OVERLAP                                 !
359     !``````````````````````````````````````````````````````````````````````!
360           LOGICAL FUNCTION ATEST_Y1(a,b,fa,fb)
361     
362           double precision :: a, b, fa, fb
363           double precision :: lMin, lMax, p0, p1, rad
364     
365           p0 = -a*v0(1) + b*v0(3)
366           p1 = -a*v1(1) + b*v1(3)
367     
368           if(p0<p1) then; lMIN=p0; lMAX=p1
369           else; lMIN=p1; lMAX=p0; endif
370     
371           rad=fa*phalfsize(1) + fb*phalfsize(3)
372           ATEST_Y1=(lmin>rad .OR. lmax<-rad)
373     
374           END FUNCTION ATEST_Y1
375     
376     !``````````````````````````````````````````````````````````````````````!
377     ! Function: ATEST_Z12                                                  !
378     ! Support function for TRI_BOX_OVERLAP                                 !
379     !``````````````````````````````````````````````````````````````````````!
380           LOGICAL FUNCTION ATEST_Z12(a,b,fa,fb)
381     
382           double precision :: a, b, fa, fb
383           double precision :: lMin, lMax, p1, p2, rad
384     
385           p1 = a*v1(1) - b*v1(2)
386           p2 = a*v2(1) - b*v2(2)
387     
388           if(p2<p1) then; lMIN=p2; lMAX=p1
389           else; lMIN=p1; lMAX=p2; endif
390     
391           rad=fa*phalfsize(1) + fb*phalfsize(2)
392           ATEST_Z12=(lmin>rad .OR. lmax<-rad)
393     
394           END FUNCTION ATEST_Z12
395     
396     !``````````````````````````````````````````````````````````````````````!
397     ! Function: ATEST_Z0                                                   !
398     ! Support function for TRI_BOX_OVERLAP                                 !
399     !``````````````````````````````````````````````````````````````````````!
400           LOGICAL FUNCTION ATEST_Z0(a,b,fa,fb)
401     
402           double precision :: a, b, fa, fb
403           double precision :: lMin, lMax, p0, p1, rad
404     
405           p0 = a*v0(1) - b*v0(2)
406           p1 = a*v1(1) - b*v1(2)
407     
408           if(p0<p1) then; lMIN=p0; lMAX=p1
409           else; lMIN=p1; lMAX=p0; endif
410     
411           rad=fa*phalfsize(1) + fb*phalfsize(2)
412           ATEST_Z0=(lmin>rad .OR. lmax<-rad)
413     
414           END FUNCTION ATEST_Z0
415     
416           END SUBROUTINE TRI_BOX_OVERLAP
417     
418     
419     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
420     !                                                                      !
421     !  SUBROUTINE: CHECK_IF_PARTICLE_OVERLAPS_STL                          !
422     !  Authors: Rahul Garg                               Date: 21-Mar-2014 !
423     !                                                                      !
424     !  Purpose: This subroutine is special written to check if a particle  !
425     !          overlaps any of the STL faces. The routine exits on         !
426     !          detecting an overlap. It is called after initial            !
427     !          generation of lattice configuration to remove out of domain !
428     !          particles                                                   !
429     !                                                                      !
430     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
431           SUBROUTINE CHECK_IF_PARTICLE_OVERLAPS_STL(POS, fI, fJ, fK, REMOVE)
432     
433           USE compar
434           USE constant
435           USE cutcell
436           USE desgrid
437           USE discretelement, ONLY: dimn, xe, yn, zt, max_radius
438           USE functions
439           USE geometry
440           USE indices
441           USE param1
442           USE run
443           USE stl
444     
445           Implicit none
446     
447           DOUBLE PRECISION, INTENT(IN) :: POS(DIMN)
448           INTEGER, INTENT(IN) :: fI, fJ, fK
449           LOGICAL, INTENT(OUT) :: REMOVE
450     
451     ! Integers mapping the fluid cell corners to DES Grid indices.
452           INTEGER :: I1, I2, J1, J2, K1, K2
453     
454           INTEGER I, J, K, IJK, NF, LC
455     
456           DOUBLE PRECISION :: LINE_T
457           DOUBLE PRECISION :: RADSQ, DIST(3)
458     
459           double precision :: vv(3)
460     
461           REMOVE = .TRUE.
462     
463           I1 = IofPOS(XE(fI-1))
464           I2 = IofPOS(XE( fI ))
465     
466           J1 = JofPOS(YN(fJ-1))
467           J2 = JofPOS(YN( fJ ))
468     
469           K1 = KofPOS(ZT(fK-1))
470           K2 = KofPOS(ZT( fK ))
471     
472           RADSQ = (1.05d0*MAX_RADIUS)**2
473     
474           DO K = K1, K2
475           DO J = J1, J2
476           DO I = I1, I2
477     
478              IF(.NOT.DG_is_ON_myPE_plus1layers(I,J,K)) CYCLE
479     
480              IJK = DG_FUNIJK(I,J,K)
481     
482     ! The point is on the non-fluid side of the plane if t>0
483              DO LC = 1, FACETS_AT_DG(IJK)%COUNT
484                 NF = FACETS_AT_DG(IJK)%ID(LC)
485     
486                 vv = VERTEX(1,:,NF)
487                 CALL INTERSECTLNPLANE(POS, NORM_FACE(:,NF), &
488                    vv, NORM_FACE(:,NF), LINE_T)
489     
490     ! Orthogonal projection puts the point outside of the domain or less than
491     ! one particle radius to the facet.
492                 DIST = LINE_T*NORM_FACE(:,NF)
493                 IF(LINE_T > ZERO .OR. dot_product(DIST,DIST)<=RADSQ) RETURN
494     
495              ENDDO
496           ENDDO
497           ENDDO
498           ENDDO
499     
500           REMOVE = .FALSE.
501     
502           RETURN
503           END SUBROUTINE CHECK_IF_PARTICLE_OVERLAPS_STL
504     
505           END MODULE STL_FUNCTIONS_DES
506