33
34
35
38
39
40
41#include "implicit_f.inc"
42
43
44
45
46
47
48 INTEGER IBUFN(*),ITY,ID,ADMID
49 my_real ms(*),amasu,x(3,*),addmas
50 CHARACTER(LEN=NCHARTITLE)::TITR
51
52
53
55 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
56 . x21,y21,z21,x31,y31,z31,x42,y42,z42,
57 . x32,y32,z32,e3x,e3y,e3z,smass,sum,
area,
58 . a2,b2,c2,aa,bb,cc,ang1,ang2,ang3
59
60 IF(ity /= 7)THEN
61
62
63
64 x1=x(1,ibufn(1))
65 y1=x(2,ibufn(1))
66 z1=x(3,ibufn(1))
67 x2=x(1,ibufn(2))
68 y2=x(2,ibufn(2))
69 z2=x(3,ibufn(2))
70 x3=x(1,ibufn(3))
71 y3=x(2,ibufn(3))
72 z3=x(3,ibufn(3))
73 x4=x(1,ibufn(4))
74 y4=x(2,ibufn(4))
75 z4=x(3,ibufn(4))
76
77 x21=x2-x1
78 y21=y2-y1
79 z21=z2-z1
80 x31=x3-x1
81 y31=y3-y1
82 z31=z3-z1
83 x42=x4-x2
84 y42=y4-y2
85 z42=z4-z2
86
87 e3x=y31*z42-z31*y42
88 e3y=z31*x42-x31*z42
89 e3z=x31*y42-y31*x42
90 sum=sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
92
93
94
96
97 ms(ibufn(1))=ms(ibufn(1)) + smass * fourth
98 ms(ibufn(2))=ms(ibufn(2)) + smass * fourth
99 ms(ibufn(3))=ms(ibufn(3)) + smass * fourth
100 ms(ibufn(4))=ms(ibufn(4)) + smass * fourth
101
102 addmas = addmas + smass
103
104 ELSE IF(ity == 7)THEN
105
106
107
108 x1=x(1,ibufn(1))
109 y1=x(2,ibufn(1))
110 z1=x(3,ibufn(1))
111 x2=x(1,ibufn(2))
112 y2=x(2,ibufn(2))
113 z2=x(3,ibufn(2))
114 x3=x(1,ibufn(3))
115 y3=x(2,ibufn(3))
116 z3=x(3,ibufn(3))
117
118 x21=x2-x1
119 y21=y2-y1
120 z21=z2-z1
121 x31=x3-x1
122 y31=y3-y1
123 z31=z3-z1
124 x32=x3-x2
125 y32=y3-y2
126 z32=z3-z2
127
128 a2 = x21**2 + y21**2 + z21**2
129 b2 = x32**2 + y32**2 + z32**2
130 c2 = x31**2 + y31**2 + z31**2
131 aa = sqrt(a2)
132 bb = sqrt(b2)
133 cc = sqrt(c2)
134
135 ang1 = acos((a2 + c2 - b2)/(two * aa * cc)) / pi
136 ang2 = acos((a2 + b2 - c2)/(two * aa * bb)) / pi
137 ang3 = acos((b2 + c2 - a2)/(two * bb * cc)) / pi
138
139 IF ( ( (a2 + c2 - b2)/(2. * aa * cc) <= -one ) .OR.
140 . ( (a2 + c2 - b2)/(2. * aa * cc) >= one ) .OR.
141 . ( (a2 + b2 - c2)/(2. * aa * bb) <= -one ) .OR.
142 . ( (a2 + b2 - c2)/(2. * aa * bb) >= one ) .OR.
143 . ( (b2 + c2 - a2)/(2. * bb * cc) <= -one ) .OR.
144 . ( (b2 + c2 - a2)/(2. * bb * cc) >= one ) ) THEN
146 . msgtype=msgerror,
147 . anmode=aninfo,
148 . i1=admid,
150 ENDIF
151
152 e3x=y21*z31-z21*y31
153 e3y=z21*x31-x21*z31
154 e3z=x21*y31-y21*x31
155 sum=sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
157
159
160 ms(ibufn(1))=ms(ibufn(1)) + smass * ang1
161 ms(ibufn(2))=ms(ibufn(2)) + smass * ang2
162 ms(ibufn(3))=ms(ibufn(3)) + smass * ang3
163
164
165
166 addmas = addmas + smass
167 END IF
168
169 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
integer, parameter nchartitle
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)