32 use element_mod , only : nixq
33
34
35
36
37
38
39
40
41#include "implicit_f.inc"
42
43
44
45
46
47
48
49 INTEGER, INTENT(IN) :: NEL, NFT, JALE, SYM, IXQ(NIXQ, *)
51 . xgrid(3, *), wgrid(*)
52 my_real,
INTENT(OUT) :: wfac(3, 4, nel), surf(4, nel)
54
55
56
57 INTEGER :: II, NODE1, NODE2, NODE3, NODE4, KFACE
59 . x1(3), x2(3), x3(3), x4(3),
60 . w1(3), w2(3), w3(3), w4(3)
62
63
64 DO ii = 1, nel
65 node1 = ixq(2, ii + nft)
66 node2 = ixq(3, ii + nft)
67 node3 = ixq(4, ii + nft)
68 node4 = ixq(5, ii + nft)
69 x1(1:3) = xgrid(1:3, node1)
70 x2(1:3) = xgrid(1:3, node2)
71 x3(1:3) = xgrid(1:3, node3)
72 x4(1:3) = xgrid(1:3, node4)
73 IF (jale /= 0) THEN
74
75 w1(1:3) = wgrid(3 * (node1 - 1) + 1 : 3 * (node1 - 1) + 3)
76 w2(1:3) = wgrid(3 * (node2 - 1) + 1 : 3 * (node2 - 1) + 3)
77 w3(1:3) = wgrid(3 * (node3 - 1) + 1 : 3 * (node3 - 1) + 3)
78 w4(1:3) = wgrid(3 * (node4 - 1) + 1 : 3 * (node4 - 1) + 3)
79 ELSE
80 w1(1:3) = zero
81 w2(1:3) = zero
82 w3(1:3) = zero
83 w4(1:3) = zero
84 ENDIF
85
86
87 kface = 1
88 norm(1, kface, ii) = zero
89 norm(2, kface, ii) = x2(3) - x1(3)
90 norm(3, kface, ii) = -(x2(2) - x1(2))
91 ny =>
norm(2, kface, ii)
92 nz =>
norm(3, kface, ii)
93 surf(kface, ii) = sqrt(ny * ny + nz * nz)
94 ny = ny / surf(kface, ii)
95 nz = nz / surf(kface, ii)
96 IF (sym == 1) THEN
97 surf(kface, ii) = surf(kface, ii) * half * (x1(2) + x2(2))
98 ENDIF
99
100 kface = 2
101 norm(1, kface, ii) = zero
102 norm(2, kface, ii) = x3(3) - x2(3)
103 norm(3, kface, ii) = -(x3(2) - x2(2))
104 ny =>
norm(2, kface, ii)
105 nz =>
norm(3, kface, ii)
106 surf(kface, ii) = sqrt(ny * ny + nz * nz)
107 ny = ny / surf(kface, ii)
108 nz = nz / surf(kface, ii)
109 IF (sym == 1) THEN
110 surf(kface, ii) = surf(kface, ii) * half * (x2(2) + x3(2))
111 ENDIF
112
113 kface = 3
114 norm(1, kface, ii) = zero
115 norm(2, kface, ii) = x4(3) - x3(3)
116 norm(3, kface, ii) = -(x4(2) - x3(2))
117 ny =>
norm(2, kface, ii)
118 nz =>
norm(3, kface, ii)
119 surf(kface, ii) = sqrt(ny * ny + nz * nz)
120 ny = ny / surf(kface, ii)
121 nz = nz / surf(kface, ii)
122 IF (sym == 1) THEN
123 surf(kface, ii) = surf(kface, ii) * half * (x3(2) + x4(2))
124 ENDIF
125
126 kface = 4
127 norm(1, kface, ii) = zero
128 norm(2, kface, ii) = x1(3) - x4(3)
129 norm(3, kface, ii) = -(x1(2) - x4(2))
130 ny =>
norm(2, kface, ii)
131 nz =>
norm(3, kface, ii)
132 surf(kface, ii) = sqrt(ny * ny + nz * nz)
133 ny = ny / surf(kface, ii)
134 nz = nz / surf(kface, ii)
135 IF (sym == 1) THEN
136 surf(kface, ii) = surf(kface, ii) * half * (x4(2) + x1(2))
137 ENDIF
138
139 wfac(1:3, 1, ii) = half * (w1(1:3) + w2(1:3))
140
141 wfac(1:3, 2, ii) = half * (w2(1:3) + w3(1:3))
142
143 wfac(1:3, 3, ii) = half * (w3(1:3) + w4(1:3))
144
145 wfac(1:3, 4, ii) = half * (w4(1:3) + w1(1:3))
146 ENDDO
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB