33 use element_mod , only : nixq
34
35
36
37#include "implicit_f.inc"
38#include "comlock.inc"
39
40
41
42#include "com04_c.inc"
43#include "vect01_c.inc"
44#include "param_c.inc"
45#include "tabsiz_c.inc"
46
47
48
49 INTEGER, INTENT(IN) :: IXQ(NIXQ,NUMELQ)
50 INTEGER, INTENT(INOUT) :: LGAUGE(3,NBGAUGE)
51 my_real,
INTENT(IN) :: phi(sphi),x(3,numnod)
52 my_real,
INTENT(INOUT) :: gauge(llgauge,nbgauge)
53 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
54
55
56
57 INTEGER I,II,J,N,IG,
59 . det_12_1p,det_23_2p,det_34_3p,det_41_4p,
60 . y1,z1,y2,z2,y3,z3,y4,z4
61
62
63
64
65
66
67 DO i=1,llt
68 ii=i+nft
69 iad2 = ale_connect%ee_connect%iad_connect(ii)
70 DO j=1,4
71 n= ale_connect%ee_connect%connected(iad2 + j - 1)
72 IF(n<=0)cycle
73 ig=nint(phi(n))
74 IF(ig==0)cycle
75
76
77 xg = gauge(2,ig)
78 yg = gauge(3,ig)
79 zg = gauge(4,ig)
80
81
82 y1 = x(2,ixq(2,ii))
83 z1 = x(3,ixq(2,ii))
84 y2 = x(2,ixq(3,ii))
85 z2 = x(3,ixq(3,ii))
86 y3 = x(2,ixq(4,ii))
87 z3 = x(3,ixq(4,ii))
88 y4 = x(2,ixq(5,ii))
89 z4 = x(3,ixq(5,ii))
90
91 det_12_1p = (y2-y1)*(zg-z1)-(z2-z1)*(yg-y1)
92 det_23_2p = (y3-y2)*(zg-z2)-(z3-z2)*(yg-y2)
93 det_34_3p = (y4-y3)*(zg-z3)-(z4-z3)*(yg-y3)
94 det_41_4p = (y1-y4)*(zg-z4)-(z1-z4)*(yg-y4)
95
97 IF ((det_12_1p >=0 .AND. det_23_2p >=0 .AND.
98 . det_34_3p >=0 .AND. det_41_4p >=0) .OR.
99 . (det_12_1p <=0 .AND. det_23_2p <=0 .AND.
100 . det_34_3p <=0 .AND. det_41_4p <=0))
alpha = one
101
102#include "lockon.inc"
103 IF(
alpha == one .AND.
alpha >= gauge(5,ig))
THEN
104
105
106 lgauge(1,ig)=-ii
108 ENDIF
109#include "lockoff.inc"
110 ENDDO
111 ENDDO
112
113
114 RETURN