44
45
46
48
49
50
51#include "implicit_f.inc"
52
53
54
55#include "mvsiz_p.inc"
56
57
58
59#include "param_c.inc"
60#include "com04_c.inc"
61
62
63
64 INTEGER :: IGEO(NPROPGI,NUMGEO), NGL(NEL), (NEL)
65 INTEGER :: NEL, JEUL, NXREF, IMULTI_FVM
66
68 . vol(nel), veul(lveul,*) , geo(npropg,numgeo),
69 . jac1(nel), jac2(nel), jac3(nel), jac4(nel), jac5(nel), jac6(nel),
70 . jac12(nel), jac45(nel), jac78(nel),
71 . px1(nel), px2(nel), px3(nel), px4(nel),
72 . py1(nel), py2(nel), py3(nel), py4(nel),
73 . pz1(nel), pz2(nel), pz3(nel), pz4(nel), det(nel)
74 double precision
75 . xd1(mvsiz), xd2(mvsiz), xd3(mvsiz), xd4(mvsiz),
76 . xd5(mvsiz), xd6(mvsiz), xd7(mvsiz), xd8(mvsiz),
77 . yd1(mvsiz), yd2(mvsiz), yd3(mvsiz), yd4(mvsiz),
78 . yd5(mvsiz), yd6(mvsiz), yd7(mvsiz), yd8(mvsiz),
79 . zd1(mvsiz), zd2(mvsiz), zd3(mvsiz), zd4(mvsiz),
80 . zd5(mvsiz), zd6(mvsiz), zd7(mvsiz), zd8(mvsiz),
81 . voldp(nel)
82
83
84
85 INTEGER :: I
86
87 my_real :: jac7(mvsiz), jac8(mvsiz) , jac9(mvsiz)
89 . x_17_46, x_28_35,
90 . y_17_46, y_28_35,
91 . z_17_46, z_28_35
93 . dett(mvsiz),
94 . jaci1(mvsiz), jaci2(mvsiz), jaci3(mvsiz),
95 . jaci4(mvsiz), jaci5(mvsiz), jaci6(mvsiz),
96 . jaci7(mvsiz), jaci8(mvsiz), jaci9(mvsiz),
97 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz)
98 double precision
99 . x17, x28, x35, x46,
100 . y17, y28, y35, y46,
101 . z17, z28, z35, z46
102
103
104
105
106
107 DO i=1,nel
108 x17=xd7(i)-xd1(i)
109 x28=xd8(i)-xd2(i)
110 x35=xd5(i)-xd3(i)
111 x46=xd6(i)-xd4(i)
112
113 y17=yd7(i)-yd1(i)
114 y28=yd8(i)-yd2(i)
115 y35=yd5(i)-yd3(i)
116 y46=yd6(i)-yd4(i)
117
118 z17=zd7(i)-zd1(i)
119 z28=zd8(i)-zd2(i)
120 z35=zd5(i)-zd3(i)
121 z46=zd6(i)-zd4(i)
122
123 jac1(i)=x17+x28-x35-x46
124 jac2(i)=y17+y28-y35-y46
125 jac3(i)=z17+z28-z35-z46
126
127 x_17_46=x17+x46
128 x_28_35=x28+x35
129 y_17_46=y17+y46
130 y_28_35=y28+y35
131 z_17_46=z17+z46
132 z_28_35=z28+z35
133
134 jac4(i)=x_17_46+x_28_35
135 jac5(i)=y_17_46+y_28_35
136 jac6(i)=z_17_46+z_28_35
137
138 jac7(i)=x_17_46-x_28_35
139 jac8(i)=y_17_46-y_28_35
140 jac9(i)=z_17_46-z_28_35
141 ENDDO
142
143 DO i=1,nel
144 jac_59_68(i)=jac5(i)*jac9(i)-jac6(i)*jac8(i)
145 jac_67_49(i)=jac6(i)*jac7(i)-jac4(i)*jac9(i
146 jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
147 ENDDO
148
149 DO i=1,nel
150 voldp(i)=one_over_64*(jac1(i)*jac_59_68(i)+jac2(i)*jac_67_49(i)+jac3(i)*jac_48_57(i))
151 det(i)=voldp(i)
152 vol(i)=det(i)
153 ENDDO
154
155 IF(jeul * (1 - imulti_fvm) /= 0)THEN
156 DO i=1,nel
157 veul(32,i) = vol(i)
158 ENDDO
159 ENDIF
160
161 DO i=1,nel
162 IF (det(i) > zero) cycle
163 IF (igeo(11,ngeo(i)) /= 0 .AND. igeo(11,ngeo(i)) /= 43) THEN
165 . msgtype=msgerror,
166 . anmode=aninfo,
167 . i1=ngl(i))
168 ELSE
170 . msgtype=msgwarning,
171 . anmode=aninfo,
172 . i1=ngl(i))
173 ENDIF
174 ENDDO
175
176 IF( jeul==0 .AND. nxref==0) RETURN
177
178 DO i=1,nel
179 dett(i)=one_over_64/
max(det(i),em20)
180 ENDDO
181
182
183 DO i=1,nel
184 jaci1(i)=dett(i)*jac_59_68(i)
185 jaci4(i)=dett(i)*jac_67_49(i)
186 jaci7(i)=dett(i)*jac_48_57(i)
187 jaci2(i)=dett(i)*(-jac2(i)*jac9(i)+jac3(i)*jac8(i))
188 jaci5(i)=dett(i)*( jac1(i)*jac9(i)-jac3(i)*jac7(i))
189 jaci8(i)=dett(i)*(-jac1(i)*jac8(i)+jac2(i)*jac7(i))
190 jaci3(i)=dett(i)*( jac2(i)*jac6(i)-jac3(i)*jac5(i))
191 jaci6(i)=dett(i)*(-jac1(i)*jac6(i)+jac3(i)*jac4(i))
192 jaci9(i)=dett(i)*( jac1(i)*jac5(i)-jac2(i)*jac4(i))
193 ENDDO
194
195
196 DO i=1,nel
197 jac12(i)=jaci1(i)-jaci2(i)
198 jac45(i)=jaci4(i)-jaci5(i)
199 jac78(i)=jaci7(i)-jaci8(i)
200 ENDDO
201
202 DO i=1,nel
203 px3(i)= jac12(i)+jaci3(i)
204 py3(i)= jac45(i)+jaci6(i)
205 pz3(i)= jac78(i)+jaci9(i)
206 px4(i)= jac12(i)-jaci3(i)
207 py4(i)= jac45(i)-jaci6(i)
208 pz4(i)= jac78(i)-jaci9(i)
209 ENDDO
210
211 DO i=1,nel
212 jac12(i)=jaci1(i)+jaci2(i)
213 jac45(i)=jaci4(i)+jaci5(i)
214 jac78(i)=jaci7(i)+jaci8(i)
215 ENDDO
216
217 DO i=1,nel
218 px1(i)=-jac12(i)-jaci3(i)
219 py1(i)=-jac45(i)-jaci6(i)
220 pz1(i)=-jac78(i)-jaci9(i)
221 px2(i)=-jac12(i)+jaci3(i)
222 py2(i)=-jac45(i)+jaci6(i)
223 pz2(i)=-jac78(i)+jaci9(i)
224 ENDDO
225
226 IF(jeul * (1 - imulti_fvm) /= 0)THEN
227 DO i=1,nel
228 veul(1,i) = px1(i)
229 veul(2,i) = px2(i)
230 veul(3,i) = px3(i)
231 veul(4,i) = px4(i)
232 veul(5,i) = py1(i)
233 veul(6,i) = py2(i)
234 veul(7,i) = py3(i)
235 veul(8,i) = py4(i)
236 veul(9,i) = pz1(i)
237 veul(10,i)= pz2(i)
238 veul(11,i)= pz3(i)
239 veul(12,i)= pz4(i)
240 END DO
241
242 IF (igeo(11,ngeo(1)) == 15) THEN
243 DO i=1,nel
244 vol(i)=vol(i)*geo(1,ngeo(i))
245 END DO
246 END IF
247 END IF
248
249 RETURN
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)