30
31
32
33
34
35
36
38 use element_mod , only : nixq
39
40
41
42#include "implicit_f.inc"
43
44
45
46#include "mvsiz_p.inc"
47
48
49
50#include "com01_c.inc"
51#include "com04_c.inc"
52#include "vect01_c.inc"
53
54
55
56 INTEGER,INTENT(IN) :: IXQ(7,NUMELQ),NEL
57 my_real,
INTENT(IN) :: x(3,numnod)
58 my_real,
INTENT(INOUT) :: grad(nel,4)
59 TYPE(t_ale_connectivity), INTENT(INOUT) :: ALE_CONNECTIVITY
60
61
62
63 INTEGER I, II, IE, IV1, IV2, IV3, IV4, IAD1
64 my_real y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
65 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
66 . yc(mvsiz), zc(mvsiz),
67 . n1y(mvsiz), n2y(mvsiz), n3y(mvsiz), n4y(mvsiz),
68 . n1z(mvsiz), n2z(mvsiz), n3z(mvsiz), n4z(mvsiz),
69 . dd1(mvsiz), dd2(mvsiz), dd3(mvsiz), dd4(mvsiz),
70 . d1y(mvsiz), d2y(mvsiz), d3y(mvsiz), d4y(mvsiz),
71 . d1z(mvsiz), d2z(mvsiz), d3z(mvsiz), d4z(mvsiz)
72
73
74
75
76
77 DO i=lft,llt
78 ii=i+nft
79
80 y1(i)=x(2,ixq(2,ii))
81 z1(i)=x(3,ixq(2,ii))
82
83 y2(i)=x(2,ixq(3,ii))
84 z2(i)=x(3,ixq(3,ii))
85
86 y3(i)=x(2,ixq(4,ii))
87 z3(i)=x(3,ixq(4,ii))
88
89 y4(i)=x(2,ixq(5,ii))
90 z4(i)=x(3,ixq(5,ii))
91 END DO
92
93
94 DO i=lft,llt
95 n1y(i)= (z2(i)-z1(i))
96 n1z(i)=-(y2(i)-y1(i))
97
98 n2y(i)= (z3(i)-z2(i))
99 n2z(i)=-(y3(i)-y2(i))
100
101 n3y(i)= (z4(i)-z3(i))
102 n3z(i)=-(y4(i)-y3(i))
103
104 n4y(i)= (z1(i)-z4(i))
105 n4z(i)=-(y1(i)-y4(i))
106
107 yc(i) = (y1(i)+y2(i)+y3(i)+y4(i))
108 zc(i) = (z1(i)+z2(i)+z3(i)+z4(i))
109 END DO
110
111 IF(n2d == 1)THEN
112 DO i=lft,llt
113 n1y(i)= n1y(i)*(y1(i)+y2(i))*0.5
114 n1z(i)= n1z(i)*(y1(i)+y2(i))*0.5
115 n2y(i)= n2y(i)*(y2(i)+y3(i))*0.5
116 n2z(i)= n2z(i)*(y2(i)+y3(i))*0.5
117 n3y(i)= n3y(i)*(y3(i)+y4(i))*0.5
118 n3z(i)= n3z(i)*(y3(i)+y4(i))*0.5
119 n4y(i)= n4y(i)*(y1(i)+y4(i))*0.5
120 n4z(i)= n4z(i)*(y1(i)+y4(i))*0.5
121 END DO
122 ENDIF
123
124
125 DO i=lft,llt
126 ie =nft+i
127 iad1 = ale_connectivity%ee_connect%iad_connect(ie)
128 iv1 = ale_connectivity%ee_connect%connected(iad1 + 1 - 1)
129 iv2 = ale_connectivity%ee_connect%connected(iad1 + 2 - 1)
130 iv3 = ale_connectivity%ee_connect%connected(iad1 + 3 - 1)
131 iv4 = ale_connectivity%ee_connect%connected(iad1 + 4 - 1)
132
133 IF(iv1 <= 0) iv1=ie
134 IF(iv2 <= 0) iv2=ie
135 IF(iv3 <= 0) iv3=ie
136 IF(iv4 <= 0) iv4=ie
137
138 d1y(i) = - yc(i)+x(2,ixq(2,iv1))+x(2,ixq(3,iv1))+x(2,ixq(4,iv1))+x(2,ixq(5,iv1))
139 d1z(i) = - zc(i)+x(3,ixq(2,iv1))+x(3,ixq(3,iv1))+x(3,ixq(4,iv1))+x(3,ixq(5,iv1))
140
141 d2y(i) = - yc(i)+x(2,ixq(2,iv2))+x(2,ixq(3,iv2))+x(2,ixq(4,iv2))+x(2,ixq(5,iv2))
142 d2z(i) = - zc(i)+x(3,ixq(2,iv2))+x(3,ixq(3,iv2))+x(3,ixq(4,iv2))+x(3,ixq(5,iv2))
143
144 d3y(i) = - yc(i)+x(2,ixq(2,iv3))+x(2,ixq(3,iv3))+x(2,ixq(4,iv3))+x(2,ixq(5,iv3))
145 d3z(i) = - zc(i)+x(3,ixq(2,iv3))+x(3,ixq(3,iv3))+x(3,ixq(4,iv3))+x(3,ixq(5,iv3))
146
147 d4y(i) = - yc(i)+x(2,ixq(2,iv4))+x(2,ixq(3,iv4))+x(2,ixq(4,iv4))+x(2,ixq(5,iv4))
148 d4z(i) = - zc(i)+x(3,ixq(2,iv4))+x(3,ixq(3,iv4))+x(3,ixq(4,iv4))+x(3,ixq(5,iv4))
149 END DO
150
151 DO i=lft,llt
152 dd1(i)=d1y(i)**2+d1z(i)**2
153 dd2(i)=d2y(i)**2+d2z(i)**2
154 dd3(i)=d3y(i)**2+d3z(i)**2
155 dd4(i)=d4y(i)**2+d4z(i)**2
156 END DO
157
158
159 DO i=lft,llt
160 grad(i,1)= four*(d1y(i)*n1y(i)+d1z(i)*n1z(i)) /
max(em15,dd1(i))
161 grad(i,2)= four*(d2y(i)*n2y(i)+d2z(i)*n2z(i)) /
max(em15,dd2(i))
162 grad(i,3)= four*(d3y(i)*n3y(i)+d3z(i)*n3z(i)) /
max(em15,dd3(i))
163 grad(i,4)= four*(d4y(i)*n4y(i)+d4z(i)*n4z(i)) /
max(em15,dd4(i))
164 END DO
165
166 RETURN