1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
SUBROUTINE Q2BCD (EST,PLANAR,RMAT,ET,IERROR)
C
C BASIC CALCULATIONS ARE PERFORMED FOR THE QDMEM2 ELEMENT IN THIS
C ROUTINE (DOUBLE-PRECISION VERSION)
C
LOGICAL PLANAR
REAL EST(1)
DOUBLE PRECISION MAG ,D12(3) ,G1(3) ,IAREA ,D13(3),GRID(3,5),
1 G2(3) ,ITWOH ,D24(3),VEC(3) ,G3(3) ,ET(3,3) ,
2 G5(3) ,G4(3) ,DADOTB,RMAT(3,5)
EQUIVALENCE (GRID(1,1),G1(1)),(GRID(1,2),G2(1)),
1 (GRID(1,3),G3(1)),(GRID(1,4),G4(1)),
2 (GRID(1,5),G5(1))
C
C MOVE GRID COORDINATES AND MAKE DOUBLE-PRECISION IF THIS IS THE
C DOUBLE-PRECISION VERSION.
C
DO 10 I = 1,3
G1(I) = EST(I+10)
G2(I) = EST(I+14)
G3(I) = EST(I+18)
G4(I) = EST(I+22)
10 CONTINUE
C
C FORM D , D AND D VECTORS
C 13 24 12
C
DO 20 I = 1,3
D12(I) = G2(I) - G1(I)
D13(I) = G3(I) - G1(I)
D24(I) = G4(I) - G2(I)
20 CONTINUE
C
C NVEC = D13 CROSS D24 = K-VECTOR (UN-NORMALIZED)
C
CALL DAXB (D13,D24,VEC)
MAG = DSQRT(DADOTB(VEC,VEC))
IAREA = 0.5D0*MAG
C
C NORMALIZE K-VECTOR
C
IF (MAG) 100,100,30
30 ET(1,3) = VEC(1)/MAG
ET(2,3) = VEC(2)/MAG
ET(3,3) = VEC(3)/MAG
C
C H = .5 * (D DOT K-VEC)
C 12
C
ITWOH = DADOTB(D12,ET(1,3))
C
C I-VECTOR (UN-NORMALIZED) = (D ) - 2 H (K-VECTOR)
C 12
C
DO 40 I = 1,3
VEC(I) = D12(I) - ITWOH*ET(I,3)
40 CONTINUE
MAG = DSQRT(DADOTB(VEC,VEC))
C
C NORMALIZE I-VECTOR
C
IF (MAG) 100,100,50
50 ET(1,1) = VEC(1)/MAG
ET(2,1) = VEC(2)/MAG
ET(3,1) = VEC(3)/MAG
C
C JVEC = KVEC CROSS IVEC
C
CALL DAXB (ET(1,3),ET(1,1),ET(1,2))
C
C FILL THE SUB-TRIANGLE ELEMENT COORDINATE MATRIX
C
DO 60 I = 1,3
G5(I) = 0.25D0*(G1(I) + G2(I) + G3(I) + G4(I))
60 CONTINUE
RMAT(1,1) = 0.0D0
RMAT(2,1) = 0.0D0
RMAT(3,1) =-ITWOH/2.0D0
DO 70 I = 2,5
VEC(1) = GRID(1,I) - G1(1)
VEC(2) = GRID(2,I) - G1(2)
VEC(3) = GRID(3,I) - G1(3)
CALL GMMATD (ET,3,3,0, VEC,3,1,0, RMAT(1,I))
RMAT(1,I) = RMAT(1,I) + RMAT(1,1)
RMAT(2,I) = RMAT(2,I) + RMAT(2,1)
RMAT(3,I) = RMAT(3,I) + RMAT(3,1)
70 CONTINUE
C
C SET PLANAR FLAG .TRUE. OR .FALSE.
C
IF ((ITWOH/2.0D0)**2/IAREA .LE. 0.01D0) GO TO 80
PLANAR = .FALSE.
GO TO 90
80 PLANAR = .TRUE.
C
C ALL BASIC CALCULATIONS NOW COMPLETE
C
90 IERROR = 0
RETURN
C
C ERROR CONDITION, BAD ELEMENT GEOMETRY.
C
100 IERROR = 1
RETURN
END
|
