1 : C CQ4LOC SOURCE CHAT 05/01/12 22:27:05 5004
2 : SUBROUTINE CQ4LOC (XE, XEL,BPSS,NOQUAL, IVRF)
3 : IMPLICIT INTEGER(I-N)
4 : IMPLICIT REAL*8(A-H,O-Z)
5 : ************************************************************************
6 : *
7 : * C Q 4 L O C
8 : * -----------
9 : *
10 : * FONCTION:
11 : * ---------
12 : *
13 : * Calcul de caract{ristiques d'un {l{ment COQ4.
14 : *
15 : * PARAMETRES: (E)=ENTREE (S)=SORTIE (+ = CONTENU DANS UN COMMUN)
16 : * -----------
17 : *
18 : * XE (E) Coordonn{es des 4 noeuds.
19 : * IVRF (E) = 1 si demande de v{rification de l'{l{ment,
20 : * = 0 sinon.
21 : * XEL (S) Coordonn{es locales des 4 noeuds.
22 : * BPSS (S) Matrice de passage.
23 : * NOQUAL (S) Indice de non-qualit{:
24 : * = 0 si OK,
25 : * = 1 si noeuds trop voisins,
26 : * = 3 SI NOEUDS NON COPLANAIRES.
27 : * (fourni si demande de v{rification d'{l{ment)
28 : *
29 : REAL*8 XE(3,4),XEL(3,4),BPSS(3,3)
30 : INTEGER NOQUAL,IVRF
31 : *
32 : * CONSTANTES:
33 : * -----------
34 : *
35 : * IND4 = indi\age circulaire de 1 @ 4.
36 : *
37 : INTEGER IND4(0:5)
38 : *
39 : * VARIABLES:
40 : * ----------
41 : *
42 : * QSI1 = vecteur norm{ de la m{diane allant de 4-1 vers 2-3.
43 : * ETA1 = vecteur norm{ de la m{diane allant de 1-2 vers 3-4.
44 : * X1, Y1 = vecteurs du rep}re local, dans le plan moyen de
45 : * l'{l{ment.
46 : * Z1 = vecteur du rep}re local, normal au plan moyen de
47 : * l'{l{ment.
48 : *
49 : REAL*8 QSI1(3),ETA1(3),X1(3),Y1(3),Z1(3)
50 : REAL*8 XD(3,4),U1(3),V1(3)
51 : *
52 : * MODE DE FONCTIONNEMENT:
53 : * -----------------------
54 : *
55 : * Pour le calcul du rep}re local et de la matrice de passage, on
56 : * fait une estimation du plan moyen.
57 : *
58 : * AUTEUR, DATE DE CREATION:
59 : * -------------------------
60 : *
61 : * PASCAL MANIGOT 09 JUILLET 1991
62 : *
63 : * LANGAGE:
64 : * --------
65 : *
66 : * FORTRAN77
67 : *
68 : ************************************************************************
69 : *
70 : DATA IND4/4,1,2,3,4,1/
71 : *
72 : NOQUAL=0
73 : C
74 : C VERIFICA DISTANZA MINIMA DEI PUNTI DELL ELEMENTO
75 : C CALIBRE PAR RAPPORT AU PERIMETRE
76 : *+* A virer ?
77 : IF (IVRF.NE.1) GO TO 6
78 : PP=0.D0
79 : DO 2 I=1,4
80 : II=I+1
81 : IF(II.EQ.5) II=1
82 : C1 = ABS(XE(1,I)-XE(1,II))
83 : C2 = ABS(XE(2,I)-XE(2,II))
84 : C3 = ABS(XE(3,I)-XE(3,II))
85 : C1 = C1*C1+C2*C2+C3*C3
86 : PP = PP + SQRT(C1)
87 : 2 CONTINUE
88 : DMIN=PP/50.D0
89 : DO 3 I=1,3
90 : I1=I+1
91 : DO 3 N=I1,4
92 : IF(ABS(XE(1,I)-XE(1,N)).LE.DMIN.AND.
93 : $ ABS(XE(2,I)-XE(2,N)).LE.DMIN.AND.
94 : $ ABS(XE(3,I)-XE(3,N)).LE.DMIN) THEN
95 : NOQUAL=1
96 : RETURN
97 : ENDIF
98 : 3 CONTINUE
99 : 6 CONTINUE
100 : *
101 : * Calcul du rep}re local
102 : * ----------------------
103 : *
104 : * Y
105 : * 4 | 3
106 : * *---|---------*
107 : * | | |
108 : * | | |
109 : * | | |
110 : * | +------------X
111 : * | |
112 : * *-------------*
113 : * 1 2
114 : *
115 : *
116 : * Calcul des m{dianes:
117 : QSI1(1) = XE(1,2)+XE(1,3) - XE(1,1)-XE(1,4)
118 : QSI1(2) = XE(2,2)+XE(2,3) - XE(2,1)-XE(2,4)
119 : QSI1(3) = XE(3,2)+XE(3,3) - XE(3,1)-XE(3,4)
119 : CALL NORMER (QSI1)
120 : ETA1(1) = XE(1,3)+XE(1,4) - XE(1,1)-XE(1,2)
121 : ETA1(2) = XE(2,3)+XE(2,4) - XE(2,1)-XE(2,2)
122 : ETA1(3) = XE(3,3)+XE(3,4) - XE(3,1)-XE(3,2)
122 : CALL NORMER (ETA1)
123 : *
124 : * Normale = Normale aux 2 m{dianes.
125 : Z1(1)= QSI1(2)*ETA1(3) - QSI1(3)*ETA1(2)
126 : Z1(2)= QSI1(3)*ETA1(1) - QSI1(1)*ETA1(3)
127 : Z1(3)= QSI1(1)*ETA1(2) - QSI1(2)*ETA1(1)
127 : CALL NORMER (Z1)
128 : *
129 : * Axes dans le Plan = bissectrices des bissectrices des m{dianes
130 : * = m{dianes pour un {l{ment rectangulaire
131 : U1(1) = QSI1(1) - ETA1(1)
132 : U1(2) = QSI1(2) - ETA1(2)
133 : U1(3) = QSI1(3) - ETA1(3)
133 : CALL NORMER (U1)
134 : V1(1) = QSI1(1) + ETA1(1)
135 : V1(2) = QSI1(2) + ETA1(2)
136 : V1(3) = QSI1(3) + ETA1(3)
136 : CALL NORMER (V1)
137 : *
138 : X1(1) = U1(1) + V1(1)
139 : X1(2) = U1(2) + V1(2)
140 : X1(3) = U1(3) + V1(3)
140 : CALL NORMER (X1)
141 : *
142 : Y1(1)=X1(3)*Z1(2)-X1(2)*Z1(3)
143 : Y1(2)=X1(1)*Z1(3)-X1(3)*Z1(1)
144 : Y1(3)=X1(2)*Z1(1)-X1(1)*Z1(2)
144 : CALL NORMER (Y1)
145 : *
146 : * Coordonn{es locales
147 : * -------------------
148 : *
149 : DO 5 J=1,4
150 : DO 5 I=1,3
151 : XD(I,J)=XE(I,J)-XE(I,1)
152 : 5 CONTINUE
153 : *
154 : DO 10 J=1,4
155 : XEL(1,J) = XD(1,J)*X1(1)+XD(2,J)*X1(2)+XD(3,J)*X1(3)
156 : XEL(2,J) = XD(1,J)*Y1(1)+XD(2,J)*Y1(2)+XD(3,J)*Y1(3)
157 : XEL(3,J) = 0.D0
158 : 10 CONTINUE
159 : *
160 : * Matrice de passage
161 : * ------------------
162 : *
163 : DO 15 I=1,3
164 : BPSS(1,I)=X1(I)
165 : BPSS(2,I)=Y1(I)
166 : BPSS(3,I)=Z1(I)
167 : 15 CONTINUE
168 : *
169 : IF(IVRF.NE.1) RETURN
170 : *
171 : * Test de plan{it{
172 : * ----------------
173 : *
174 : * Calcul des 4 "normales" locales:
175 : DO 102 K=1,4
176 : KP1 = IND4(K+1)
177 : KM1 = IND4(K-1)
178 : DO 100 J=1,3
179 : U1(J) = XE(J,KP1) - XE(J,K)
180 : V1(J) = XE(J,KM1) - XE(J,K)
181 : 100 CONTINUE
182 : XD(1,K) = U1(2)*V1(3) - U1(3)*V1(2)
183 : XD(2,K) = U1(3)*V1(1) - U1(1)*V1(3)
184 : XD(3,K) = U1(1)*V1(2) - U1(2)*V1(1)
185 : XXD = (XD(1,K)**2) + (XD(2,K)**2) + (XD(3,K)**2)
186 : XXD = SQRT(XXD)
187 : DO 101 J=1,3
188 : XD(J,K) = XD(J,K) / XXD
189 : 101 CONTINUE
190 : 102 CONTINUE
191 : *
192 : * Calcul de la non-plan{it{:
193 : COS13 = XD(1,3)*XD(1,1) + XD(2,3)*XD(2,1) + XD(3,3)*XD(3,1)
194 : COS24 = XD(1,4)*XD(1,2) + XD(2,4)*XD(2,2) + XD(3,4)*XD(3,2)
195 : IF (MIN(COS13,COS24) .LT. 0.99999) THEN
196 : * Non-plan{it{ de 0.25 degr{ ou plus:
197 : NOQUAL = 3
198 : * Rq: 0.9999 , qui {quivaut @ 1 degr{, est insuffisant.
199 : END IF
200 : *
201 : END
202 :
© Cast3M 2003 - All rights reserved.
Disclaimer