OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
smorth3.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| smorth3 ../starter/source/elements/solid/solide/smorth3.F
25!||--- called by ------------------------------------------------------
26!|| multifluid_init3 ../starter/source/multifluid/multifluid_init3.F
27!|| s10init3 ../starter/source/elements/solid/solide10/s10init3.F
28!|| s4init3 ../starter/source/elements/solid/solide4/s4init3.F
29!|| s8zinit3 ../starter/source/elements/solid/solide8z/s8zinit3.F
30!|| sinit3 ../starter/source/elements/solid/solide/sinit3.F
31!||--- calls -----------------------------------------------------
32!|| ancmsg ../starter/source/output/message/message.F
33!|| fretitl2 ../starter/source/starter/freform.F
34!|| inv3 ../starter/source/tools/univ/inv3.F
35!||--- uses -----------------------------------------------------
36!|| message_mod ../starter/share/message_module/message_mod.F
37!||====================================================================
38 SUBROUTINE smorth3(PID ,GEO ,IGEO ,SKEW ,IREP ,GAMA ,
39 . RX ,RY ,RZ ,SX ,SY ,SZ ,TX ,TY ,TZ ,
40 . E1X ,E1Y ,E1Z ,E2X ,E2Y ,E2Z ,E3X ,E3Y ,E3Z ,
41 . F1X ,F1Y ,F1Z ,F2X ,F2Y ,F2Z ,NSIGI,SIGSP,NSIGS,
42 . SIGI ,IXS ,X ,JHBE ,PT ,NEL ,ISOLNOD)
43C-----------------------------------------------
44C M o d u l e s
45C-----------------------------------------------
46 USE message_mod
47 USE names_and_titles_mod , ONLY : nchartitle
48C-----------------------------------------------
49C I m p l i c i t T y p e s
50C-----------------------------------------------
51#include "implicit_f.inc"
52C-----------------------------------------------
53C A n a l y s e M o d u l e
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56C-----------------------------------------------
57C C o m m o n B l o c k s
58C-----------------------------------------------
59#include "com01_c.inc"
60#include "vect01_c.inc"
61#include "param_c.inc"
62#include "scr17_c.inc"
63C-----------------------------------------------
64C D u m m y A r g u m e n t s
65C-----------------------------------------------
66 INTEGER JHBE,IREP,NSIGI,NSIGS,NEL,ISOLNOD
67 INTEGER PID(*),IGEO(NPROPGI,*),IXS(NIXS,*),PT(*)
68C REAL
69 my_real
70 . geo(npropg,*),skew(lskew,*),gama(nel,6),
71 . rx(*) ,ry(*) ,rz(*) ,sx(*) ,sy(*) ,sz(*) ,tx(*) ,ty(*) ,tz(*),
72 . e1x(*),e1y(*),e1z(*),e2x(*),e2y(*),e2z(*),e3x(*),e3y(*),e3z(*),
73 . f1x(*),f1y(*),f1z(*),f2x(*),f2y(*),f2z(*),sigsp(nsigi,*),
74 . sigi(nsigs,*),x(3,*)
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER I,ISK,IPNUM,IG,IIS,II,J,JJ,N,IFLAGINI,N1,N2,N4,NNOD,INIORTH(MVSIZ)
79C REAL
80 my_real
81 . XL,YL,ZL,SUM,HX,HY,HZ,KX,KY,KZ,LX,LY,LZ,PHI,CP,SP,VX,VY,VZ,VN,
82 . f3x,f3y,f3z,
83 . g11,g22,g33,g12,g21,g23,g32,g13,g31,pts(3)
84 my_real
85 . sk(6),a(9),b(9)
86 INTEGER ID
87 CHARACTER(LEN=NCHARTITLE)::TITR
88C-----------------------------------------------------------------------
89C Repere orthotrope
90C Stockage de Transpose(G) tq Xortho = Transpose(G) Xcvt
91C GAMA(1)= TG11 , GAMA(2) = TG12, TG13..., TG21..., TG22..., TG23...
92C=======================================================================
93C---- tag elm /w /INIBRI/ORTHO
94 iniorth(lft:llt)=0
95 IF (nvsolid3 /= 0) THEN
96 iis= nvsolid1 + nvsolid2 + 4 +nusolid
97 DO i=lft,llt
98 jj=pt(nft+i)
99 IF(jj ==0 ) cycle
100 IF(
101 . sigsp(iis+1,jj) /= zero .OR. sigsp(iis+2,jj) /=zero .OR.
102 . sigsp(iis+3,jj) /= zero .OR. sigsp(iis+4,jj) /=zero .OR.
103 . sigsp(iis+5,jj) /= zero .OR. sigsp(iis+6,jj) /=zero )THEN
104 iniorth(i) = 1
105 ENDIF
106 ENDDO
107 ENDIF
108 DO i=lft,llt
109 IF(iniorth(i) ==1 ) cycle
110 ig = pid(i)
111 id=igeo(1,ig)
112 CALL fretitl2(titr,
113 . igeo(npropgi-ltitr+1,ig),ltitr)
114 ipnum = igeo(2,ig)
115 phi = geo(1,ig) * pi/hundred80
116 vx = geo(7,ig)
117 vy = geo(8,ig)
118 vz = geo(9,ig)
119 cp = cos(phi)
120 sp = sin(phi)
121 IF (ipnum > 20) THEN
122 pts(1:3) = geo(33:35,ig)
123 END IF
124C
125 IF (ipnum < 0) THEN
126 isk = -ipnum
127 gama(i,1)=
128 . skew(1,isk)*e1x(i)+skew(2,isk)*e1y(i)+skew(3,isk)*e1z(i)
129 gama(i,2)=
130 . skew(1,isk)*e2x(i)+skew(2,isk)*e2y(i)+skew(3,isk)*e2z(i)
131 gama(i,3)=
132 . skew(1,isk)*e3x(i)+skew(2,isk)*e3y(i)+skew(3,isk)*e3z(i)
133 gama(i,4)=
134 . skew(4,isk)*e1x(i)+skew(5,isk)*e1y(i)+skew(6,isk)*e1z(i)
135 gama(i,5)=
136 . skew(4,isk)*e2x(i)+skew(5,isk)*e2y(i)+skew(6,isk)*e2z(i)
137 gama(i,6)=
138 . skew(4,isk)*e3x(i)+skew(5,isk)*e3y(i)+skew(6,isk)*e3z(i)
139 ELSEIF (ipnum == 4) THEN
140C DIR1 = projection vect V on bottom face
141C DIR2 = orthogonal to DIR1, on bottom face
142C DIR3 = orthogonal to bottom face
143 f3x = f1y(i)*f2z(i) - f1z(i)*f2y(i)
144 f3y = f1z(i)*f2x(i) - f1x(i)*f2z(i)
145 f3z = f1x(i)*f2y(i) - f1y(i)*f2x(i)
146 sum = one / max(sqrt(f3x*f3x+f3y*f3y+f3z*f3z),em20)
147 f3x = f3x * sum
148 f3y = f3y * sum
149 f3z = f3z * sum
150C
151 vx = geo(7,ig)
152 vy = geo(8,ig)
153 vz = geo(9,ig)
154 vn = vx*f3x + vy*f3y + vz*f3z
155 vx = vx - vn*f3x
156 vy = vy - vn*f3y
157 vz = vz - vn*f3z
158 sum= sqrt(vx*vx+vy*vy+vz*vz)
159 IF (sum < em20) THEN
160 vx = f1x(i)
161 vy = f1y(i)
162 vz = f1z(i)
163 sum = one
164 ELSE
165 sum = one / max(sqrt(vx*vx+vy*vy+vz*vz),em20)
166 ENDIF
167C repere orthogononalise facette
168C E3 = E1 x E2, E1' = V, E2' = E3 x E1'
169 f1x(i) = vx * sum
170 f1y(i) = vy * sum
171 f1z(i) = vz * sum
172 f2x(i) = f3y*f1z(i) - f3z*f1y(i)
173 f2y(i) = f3z*f1x(i) - f3x*f1z(i)
174 f2z(i) = f3x*f1y(i) - f3y*f1x(i)
175C dir ortho ds repere global
176C S1 = vect F1 + rot(phi) dans le plan de facette
177C S2 = vect F2 + rot(phi) dans le plan de facette
178 sk(1) = cp*f1x(i) + sp*f2x(i)
179 sk(2) = cp*f1y(i) + sp*f2y(i)
180 sk(3) = cp*f1z(i) + sp*f2z(i)
181 sk(4) =-sp*f1x(i) + cp*f2x(i)
182 sk(5) =-sp*f1y(i) + cp*f2y(i)
183 sk(6) =-sp*f1z(i) + cp*f2z(i)
184C dir ortho ds repere elementaire
185 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
186 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
187 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
188 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
189 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
190 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
191 ELSEIF (ipnum == 5) THEN
192C DIR1 = orthogonal to bottom face
193C DIR2 = projection vect V on bottom face
194 f3x = f1y(i)*f2z(i) - f1z(i)*f2y(i)
195 f3y = f1z(i)*f2x(i) - f1x(i)*f2z(i)
196 f3z = f1x(i)*f2y(i) - f1y(i)*f2x(i)
197 sum = one / max(sqrt(f3x*f3x+f3y*f3y+f3z*f3z),em20)
198 f3x = f3x * sum
199 f3y = f3y * sum
200 f3z = f3z * sum
201C
202 vx = geo(7,ig)
203 vy = geo(8,ig)
204 vz = geo(9,ig)
205 vn = vx*f3x + vy*f3y + vz*f3z
206 vx = vx - vn*f3x
207 vy = vy - vn*f3y
208 vz = vz - vn*f3z
209 sum= sqrt(vx*vx+vy*vy+vz*vz)
210 IF (sum < em20) THEN
211 vx = f1x(i)
212 vy = f1y(i)
213 vz = f1z(i)
214 sum = one
215 ELSE
216 sum= one / max(sqrt(vx*vx+vy*vy+vz*vz),em20)
217 ENDIF
218C repere orthogononalise facette
219C E3 = E1 x E2, E1' = V, E2' = E3 x E1'
220 f1x(i) = vx * sum
221 f1y(i) = vy * sum
222 f1z(i) = vz * sum
223 f2x(i) = f3y*f1z(i) - f3z*f1y(i)
224 f2y(i) = f3z*f1x(i) - f3x*f1z(i)
225 f2z(i) = f3x*f1y(i) - f3y*f1x(i)
226C dir ortho ds repere global
227C S1 = vect F1 + rot(phi) dans le plan de facette
228C S2 = vect F2 + rot(phi) dans le plan de facette
229 sk(1) = f3x
230 sk(2) = f3y
231 sk(3) = f3z
232 sk(4) = cp*f1x(i) + sp*f2x(i)
233 sk(5) = cp*f1y(i) + sp*f2y(i)
234 sk(6) = cp*f1z(i) + sp*f2z(i)
235C dir ortho ds repere elementaire
236 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
237 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
238 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
239 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
240 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
241 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
242 ELSEIF (ipnum == 20) THEN
243C -- g1: 12 (SK(1:3)); g2: 14(SK(4:6))
244 ii=nft+i
245 IF (isolnod == 4 .OR. isolnod == 10) THEN
246 sk(1) = tx(i)
247 sk(2) = ty(i)
248 sk(3) = tz(i)
249 sk(4) = rx(i)
250 sk(5) = ry(i)
251 sk(6) = rz(i)
252 ELSE
253 n1=ixs(2,ii)
254 n2=ixs(3,ii)
255 n4=ixs(5,ii)
256 sk(1) = x(1,n2)-x(1,n1)
257 sk(2) = x(2,n2)-x(2,n1)
258 sk(3) = x(3,n2)-x(3,n1)
259 sk(4) = x(1,n4)-x(1,n1)
260 sk(5) = x(2,n4)-x(2,n1)
261 sk(6) = x(3,n4)-x(3,n1)
262 END IF
263 sum = one / max(sqrt(sk(1)*sk(1)+sk(2)*sk(2)+sk(3)*sk(3)),em20)
264 sk(1) = sk(1) * sum
265 sk(2) = sk(2) * sum
266 sk(3) = sk(3) * sum
267 sum = sk(1)*sk(4)+sk(2)*sk(5)+sk(3)*sk(6)
268 sk(4) = sk(4) - sum*sk(1)
269 sk(5) = sk(5) - sum*sk(2)
270 sk(6) = sk(6) - sum*sk(3)
271 sum = one / max(sqrt(sk(4)*sk(4)+sk(5)*sk(5)+sk(6)*sk(6)),em20)
272 sk(4) = sk(4) * sum
273 sk(5) = sk(5) * sum
274 sk(6) = sk(6) * sum
275 IF (jcvt > 0) THEN
276C dir ortho v.s. local elem
277 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
278 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
279 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
280 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
281 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
282 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
283 ELSE
284 gama(i,1:6) = sk(1:6)
285 END IF
286 ELSEIF (ipnum == 21) THEN
287C -- g1: o-pt(o element centrid); g2:Z^g1
288 ii=nft+i
289 nnod = isolnod
290 IF (nnod==10) nnod=4
291 lx=zero
292 ly=zero
293 lz=zero
294 DO j = 1,nnod
295 n = ixs(j+1,ii)
296 lx=lx+x(1,n)
297 ly=ly+x(2,n)
298 lz=lz+x(3,n)
299 END DO
300 lx=lx/nnod
301 ly=ly/nnod
302 lz=lz/nnod
303 sk(1) = lx - pts(1)
304 sk(2) = ly - pts(2)
305 sk(3) = lz - pts(3)
306 sum = sqrt(sk(1)*sk(1)+sk(2)*sk(2)+sk(3)*sk(3))
307C---- check if SUM>0 wrong input pt position same than elem id center
308 IF (sum < em20) THEN
309 CALL ancmsg(msgid=1919,
310 . msgtype=msgerror,
311 . anmode=aninfo_blind_1,
312 . i1=id,
313 . c1=titr,
314 . i2=ipnum)
315 ELSE
316 sk(1) = sk(1) / sum
317 sk(2) = sk(2) / sum
318 sk(3) = sk(3) / sum
319 sk(4) = - sk(2)
320 sk(5) = sk(1)
321 sk(6) = zero
322 IF (jcvt > 0) THEN
323C dir ortho v.s. local elem
324 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
325 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
326 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
327 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
328 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
329 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
330 ELSE
331 gama(i,1:6) = sk(1:6)
332 END IF
333 END IF
334 ELSEIF (ipnum == 23) THEN
335C -- Vj+phi
336 sum = one / max(sqrt(sx(i)*sx(i)+sy(i)*sy(i)+sz(i)*sz(i)),em20)
337 f3x = sx(i) * sum
338 f3y = sy(i) * sum
339 f3z = sz(i) * sum
340 vx = geo(7,ig)
341 vy = geo(8,ig)
342 vz = geo(9,ig)
343 sum = vx*f3x+vy*f3y+vz*f3z
344 sk(4) = vx - sum*f3x
345 sk(5) = vy - sum*f3y
346 sk(6) = vz - sum*f3z
347 sk(1) = sk(5)*f3z - sk(6)*f3y
348 sk(2) = sk(6)*f3x - sk(4)*f3z
349 sk(3) = sk(4)*f3y - sk(5)*f3x
350 sum = sqrt(sk(1)*sk(1)+sk(2)*sk(2)+sk(3)*sk(3))
351C---- check if SUM>0 wrong input Vj, same than s
352 IF (sum < em20) THEN
353 CALL ancmsg(msgid=1920,
354 . msgtype=msgerror,
355 . anmode=aninfo_blind_1,
356 . i1=id,
357 . c1=titr,
358 . i2=ipnum)
359 ELSE
360 f1x(i) = sk(1) / sum
361 f1y(i) = sk(2) / sum
362 f1z(i) = sk(3) / sum
363 f2x(i) = f3y*f1z(i) - f3z*f1y(i)
364 f2y(i) = f3z*f1x(i) - f3x*f1z(i)
365 f2z(i) = f3x*f1y(i) - f3y*f1x(i)
366C
367 sk(1) = cp*f1x(i) + sp*f2x(i)
368 sk(2) = cp*f1y(i) + sp*f2y(i)
369 sk(3) = cp*f1z(i) + sp*f2z(i)
370 sk(4) =-sp*f1x(i) + cp*f2x(i)
371 sk(5) =-sp*f1y(i) + cp*f2y(i)
372 sk(6) =-sp*f1z(i) + cp*f2z(i)
373 IF (jcvt > 0) THEN
374C dir ortho v.s. local elem
375 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
376 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
377 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
378 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
379 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
380 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
381 ELSE
382 gama(i,1:6) = sk(1:6)
383 END IF
384 END IF
385 ELSEIF (ipnum == 24) THEN
386C -- Vj+Pt
387 ii=nft+i
388 nnod = isolnod
389 IF (nnod==10) nnod=4
390 lx=zero
391 ly=zero
392 lz=zero
393 DO j = 1,nnod
394 n = ixs(j+1,ii)
395 lx=lx+x(1,n)
396 ly=ly+x(2,n)
397 lz=lz+x(3,n)
398 END DO
399 lx=lx/nnod
400 ly=ly/nnod
401 lz=lz/nnod
402 vx = geo(7,ig)
403 vy = geo(8,ig)
404 vz = geo(9,ig)
405 sum = one / max(sqrt(vx*vx+vy*vy+vz*vz),em20)
406 vx = vx*sum
407 vy = vy*sum
408 vz = vz*sum
409 sk(1) = lx - pts(1)
410 sk(2) = ly - pts(2)
411 sk(3) = lz - pts(3)
412 sum = vx*sk(1)+vy*sk(2)+vz*sk(3)
413 f3x = sk(1) - sum*vx
414 f3y = sk(2) - sum*vy
415 f3z = sk(3) - sum*vz
416 sk(1) = vy*f3z - vz*f3y
417 sk(2) = vz*f3x - vx*f3z
418 sk(3) = vx*f3y - vy*f3x
419 sum = sqrt(sk(1)*sk(1)+sk(2)*sk(2)+sk(3)*sk(3))
420C---- check if SUM>0
421 IF (sum < em20) THEN
422 CALL ancmsg(msgid=1920,
423 . msgtype=msgerror,
424 . anmode=aninfo_blind_1,
425 . i1=id,
426 . c1=titr,
427 . i2=ipnum)
428 ELSE
429 sk(1) = sk(1) / sum
430 sk(2) = sk(2) / sum
431 sk(3) = sk(3) / sum
432 sk(4) = vx
433 sk(5) = vy
434 sk(6) = vz
435 IF (jcvt > 0) THEN
436C dir ortho v.s. local elem
437 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
438 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
439 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
440 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
441 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
442 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
443 ELSE
444 gama(i,1:6) = sk(1:6)
445 END IF
446 END IF
447 ELSEIF (jcvt == 0) THEN
448 SELECT CASE (ipnum)
449 CASE (1)
450 gama(i,1)= cp
451 gama(i,2)= sp
452 gama(i,3)= zero
453 gama(i,4)=-sp
454 gama(i,5)= cp
455 gama(i,6)= zero
456 CASE (2)
457 gama(i,1)= zero
458 gama(i,2)= cp
459 gama(i,3)= sp
460 gama(i,4)= zero
461 gama(i,5)=-sp
462 gama(i,6)= cp
463 CASE (3)
464 gama(i,1)= sp
465 gama(i,2)= zero
466 gama(i,3)= cp
467 gama(i,4)= cp
468 gama(i,5)= zero
469 gama(i,6)=-sp
470 CASE (11)
471 vn = vx*e3x(i) + vy*e3y(i) + vz*e3z(i)
472 vx = vx - vn*e3x(i)
473 vy = vy - vn*e3y(i)
474 vz = vz - vn*e3z(i)
475 sum = sqrt(vx*vx+vy*vy+vz*vz)
476 IF (sum < em10) THEN
477 CALL ancmsg(msgid=811,
478 . msgtype=msgwarning,
479 . anmode=aninfo_blind_1,
480 . i1=id,
481 . c1=titr)
482 sk(1) = e1x(i)
483 sk(2) = e1y(i)
484 sk(3) = e1z(i)
485 ELSE
486 sum = one / sum
487 sk(1) = vx * sum
488 sk(2) = vy * sum
489 sk(3) = vz * sum
490 ENDIF
491 sk(4) = e3y(i)* sk(3) - e3z(i)* sk(2)
492 sk(5) = e3z(i)* sk(1) - e3x(i)* sk(3)
493 sk(6) = e3x(i)* sk(2) - e3y(i)* sk(1)
494 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
495 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
496 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
497 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
498 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
499 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
500 CASE (12)
501 vn = vx*e1x(i) + vy*e1y(i) + vz*e1z(i)
502 vx = vx - vn*e1x(i)
503 vy = vy - vn*e1y(i)
504 vz = vz - vn*e1z(i)
505 sum = sqrt(vx*vx+vy*vy+vz*vz)
506 IF (sum < em10) THEN
507 CALL ancmsg(msgid=811,
508 . msgtype=msgwarning,
509 . anmode=aninfo_blind_1,
510 . i1=id,
511 . c1=titr)
512 sk(1) = e2x(i)
513 sk(2) = e2y(i)
514 sk(3) = e2z(i)
515 ELSE
516 sum = one / sum
517 sk(1) = vx * sum
518 sk(2) = vy * sum
519 sk(3) = vz * sum
520 ENDIF
521 sk(4) = e1y(i)* sk(3) - e1z(i)* sk(2)
522 sk(5) = e1z(i)* sk(1) - e1x(i)* sk(3)
523 sk(6) = e1x(i)* sk(2) - e1y(i)* sk(1)
524 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
525 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
526 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
527 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
528 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
529 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
530 CASE (13)
531 vn = vx*e2x(i) + vy*e2y(i) + vz*e2z(i)
532 vx = vx - vn*e2x(i)
533 vy = vy - vn*e2y(i)
534 vz = vz - vn*e2z(i)
535 sum = sqrt(vx*vx+vy*vy+vz*vz)
536 IF (sum < em10) THEN
537 CALL ancmsg(msgid=811,
538 . msgtype=msgwarning,
539 . anmode=aninfo_blind_1,
540 . i1=id,
541 . c1=titr)
542 sk(1) = e3x(i)
543 sk(2) = e3y(i)
544 sk(3) = e3z(i)
545 ELSE
546 sum= one / sum
547 sk(1) = vx * sum
548 sk(2) = vy * sum
549 sk(3) = vz * sum
550 ENDIF
551 sk(4) = e2y(i)* sk(3) - e2z(i)* sk(2)
552 sk(5) = e2z(i)* sk(1) - e2x(i)* sk(3)
553 sk(6) = e2x(i)* sk(2) - e2y(i)* sk(1)
554 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
555 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
556 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
557 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
558 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
559 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
560 END SELECT
561 ELSEIF (jcvt > 0) THEN
562C R'=R, T'=R^S, S'=T'^R
563 sum=sqrt(rx(i)**2+ry(i)**2+rz(i)**2)
564 IF (sum > zero) sum=one/sum
565 hx=rx(i)*sum
566 hy=ry(i)*sum
567 hz=rz(i)*sum
568 lx=hy*sz(i)-hz*sy(i)
569 ly=hz*sx(i)-hx*sz(i)
570 lz=hx*sy(i)-hy*sx(i)
571 sum = sqrt(lx**2+ly**2+lz**2)
572 IF (sum > zero) sum=one/sum
573 lx=lx*sum
574 ly=ly*sum
575 lz=lz*sum
576 kx=ly*hz-lz*hy
577 ky=lz*hx-lx*hz
578 kz=lx*hy-ly*hx
579 sum = sqrt(kx**2+ky**2+kz**2)
580 IF (sum > zero) sum=one/sum
581 kx=kx*sum
582 ky=ky*sum
583 kz=kz*sum
584 SELECT CASE (ipnum)
585 CASE (1)
586 sk(1)= cp*hx+sp*kx
587 sk(2)= cp*hy+sp*ky
588 sk(3)= cp*hz+sp*kz
589 sk(4)=-sp*hx+cp*kx
590 sk(5)=-sp*hy+cp*ky
591 sk(6)=-sp*hz+cp*kz
592 CASE (2)
593 sk(1)= cp*kx+sp*lx
594 sk(2)= cp*ky+sp*ly
595 sk(3)= cp*kz+sp*lz
596 sk(4)=-sp*kx+cp*lx
597 sk(5)=-sp*ky+cp*ly
598 sk(6)=-sp*kz+cp*lz
599 CASE (3)
600 sk(1)= cp*lx+sp*hx
601 sk(2)= cp*ly+sp*hy
602 sk(3)= cp*lz+sp*hz
603 sk(4)=-sp*lx+cp*hx
604 sk(5)=-sp*ly+cp*hy
605 sk(6)=-sp*lz+cp*hz
606 CASE (11)
607 vn = vx*lx + vy*ly + vz*lz
608 vx = vx - vn*lx
609 vy = vy - vn*ly
610 vz = vz - vn*lz
611 sum = sqrt(vx*vx+vy*vy+vz*vz)
612 IF (sum < em10) THEN
613 CALL ancmsg(msgid=811,
614 . msgtype=msgwarning,
615 . anmode=aninfo_blind_1,
616 . i1=id,
617 . c1=titr)
618 sk(1) = hx
619 sk(2) = hy
620 sk(3) = hz
621 ELSE
622 sum = one / sum
623 sk(1) = vx * sum
624 sk(2) = vy * sum
625 sk(3) = vz * sum
626 ENDIF
627 sk(4) = ly* sk(3) - lz* sk(2)
628 sk(5) = lz* sk(1) - lx* sk(3)
629 sk(6) = lx* sk(2) - ly* sk(1)
630 CASE (12)
631 vn = vx*hx + vy*hy + vz*hz
632 vx = vx - vn*hx
633 vy = vy - vn*hy
634 vz = vz - vn*hz
635 sum = sqrt(vx*vx+vy*vy+vz*vz)
636 IF (sum < em10) THEN
637 CALL ancmsg(msgid=811,
638 . msgtype=msgwarning,
639 . anmode=aninfo_blind_1,
640 . i1=id,
641 . c1=titr)
642 sk(1) = kx
643 sk(2) = ky
644 sk(3) = kz
645 ELSE
646 sum = one / sum
647 sk(1) = vx * sum
648 sk(2) = vy * sum
649 sk(3) = vz * sum
650 ENDIF
651 sk(4) = hy* sk(3) - hz* sk(2)
652 sk(5) = hz* sk(1) - hx* sk(3)
653 sk(6) = hx* sk(2) - hy* sk(1)
654 CASE (13)
655 vn = vx*kx + vy*ky + vz*kz
656 vx = vx - vn*kx
657 vy = vy - vn*ky
658 vz = vz - vn*kz
659 sum = sqrt(vx*vx+vy*vy+vz*vz)
660 IF (sum < em10) THEN
661 CALL ancmsg(msgid=811,
662 . msgtype=msgwarning,
663 . anmode=aninfo_blind_1,
664 . i1=id,
665 . c1=titr)
666 sk(1) = lx
667 sk(2) = ly
668 sk(3) = lz
669 ELSE
670 sum = one / sum
671 sk(1) = vx * sum
672 sk(2) = vy * sum
673 sk(3) = vz * sum
674 ENDIF
675 sk(4) = ky* sk(3) - kz* sk(2)
676 sk(5) = kz* sk(1) - kx* sk(3)
677 sk(6) = kx* sk(2) - ky* sk(1)
678 END SELECT
679 gama(i,1) = sk(1)*e1x(i) + sk(2)*e1y(i) + sk(3)*e1z(i)
680 gama(i,2) = sk(1)*e2x(i) + sk(2)*e2y(i) + sk(3)*e2z(i)
681 gama(i,3) = sk(1)*e3x(i) + sk(2)*e3y(i) + sk(3)*e3z(i)
682 gama(i,4) = sk(4)*e1x(i) + sk(5)*e1y(i) + sk(6)*e1z(i)
683 gama(i,5) = sk(4)*e2x(i) + sk(5)*e2y(i) + sk(6)*e2z(i)
684 gama(i,6) = sk(4)*e3x(i) + sk(5)*e3y(i) + sk(6)*e3z(i)
685 ENDIF
686 ENDDO
687C-----
688 IF (irep > 0) THEN
689C dir orthotrope attache au repere isoparametrique
690 DO i=lft,llt
691 a(1) = rx(i)*e1x(i) + ry(i)*e1y(i) + rz(i)*e1z(i)
692 a(2) = rx(i)*e2x(i) + ry(i)*e2y(i) + rz(i)*e2z(i)
693 a(3) = rx(i)*e3x(i) + ry(i)*e3y(i) + rz(i)*e3z(i)
694 a(4) = sx(i)*e1x(i) + sy(i)*e1y(i) + sz(i)*e1z(i)
695 a(5) = sx(i)*e2x(i) + sy(i)*e2y(i) + sz(i)*e2z(i)
696 a(6) = sx(i)*e3x(i) + sy(i)*e3y(i) + sz(i)*e3z(i)
697 a(7) = tx(i)*e1x(i) + ty(i)*e1y(i) + tz(i)*e1z(i)
698 a(8) = tx(i)*e2x(i) + ty(i)*e2y(i) + tz(i)*e2z(i)
699 a(9) = tx(i)*e3x(i) + ty(i)*e3y(i) + tz(i)*e3z(i)
700 CALL inv3(a, b)
701 a(1) = gama(i,1)
702 a(2) = gama(i,2)
703 a(3) = gama(i,3)
704 a(4) = gama(i,4)
705 a(5) = gama(i,5)
706 a(6) = gama(i,6)
707C
708 gama(i,1) = b(1)*a(1) + b(4)*a(2) + b(7)*a(3)
709 gama(i,2) = b(2)*a(1) + b(5)*a(2) + b(8)*a(3)
710 gama(i,3) = b(3)*a(1) + b(6)*a(2) + b(9)*a(3)
711 sum = one / max(em20,sqrt(gama(i,1)**2 + gama(i,2)**2 + gama(i,3)**2))
712 gama(i,1) = gama(i,1) * sum
713 gama(i,2) = gama(i,2) * sum
714 gama(i,3) = gama(i,3) * sum
715C
716 gama(i,4) = b(1)*a(4) + b(4)*a(5) + b(7)*a(6)
717 gama(i,5) = b(2)*a(4) + b(5)*a(5) + b(8)*a(6)
718 gama(i,6) = b(3)*a(4) + b(6)*a(5) + b(9)*a(6)
719 sum = one / max(em20,sqrt(gama(i,4)**2 + gama(i,5)**2 + gama(i,6)**2))
720 gama(i,4) = gama(i,4) * sum
721 gama(i,5) = gama(i,5) * sum
722 gama(i,6) = gama(i,6) * sum
723 ENDDO
724 ENDIF
725C
726C---
727 IF (nvsolid3 /= 0) THEN
728 iis= nvsolid1 + nvsolid2 + 4 + nusolid
729 DO i=lft,llt
730 ii=nft+i
731 jj=pt(ii)
732 iflagini = 1
733 IF(jj==0)iflagini = 0
734 IF(iflagini == 1 .AND.
735 . ( sigsp(iis+1,jj) /= zero.OR.sigsp(iis+2,jj)/=zero.OR.
736 . sigsp(iis+3,jj) /= zero .OR. sigsp(iis+4,jj) /= zero .OR.
737 . sigsp(iis+5,jj) /= zero .OR. sigsp(iis+6,jj) /= zero) )THEN
738 g11 = sigsp(iis+1,jj)
739 g21 = sigsp(iis+2,jj)
740 g31 = sigsp(iis+3,jj)
741 g12 = sigsp(iis+4,jj)
742 g22 = sigsp(iis+5,jj)
743 g32 = sigsp(iis+6,jj)
744 g13 = g21*g32-g31*g22
745 g23 = g31*g12-g11*g32
746 g33 = g11*g22-g21*g12
747C Change of basis : GLOBAL -> LOCAL.
748 gama(i,1)=e1x(i)*g11+e1y(i)*g21+e1z(i)*g31
749 gama(i,2)=e2x(i)*g11+e2y(i)*g21+e2z(i)*g31
750 gama(i,3)=e3x(i)*g11+e3y(i)*g21+e3z(i)*g31
751 gama(i,4)=e1x(i)*g12+e1y(i)*g22+e1z(i)*g32
752 gama(i,5)=e2x(i)*g12+e2y(i)*g22+e2z(i)*g32
753 gama(i,6)=e3x(i)*g12+e3y(i)*g22+e3z(i)*g32
754 ENDIF
755 ENDDO
756 ENDIF
757C
758C---
759 RETURN
760 END SUBROUTINE smorth3
761
inv3
subroutine inv3(a, b)
Definition inv3.F:29
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1249
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003
smorth3
subroutine smorth3(pid, geo, igeo, skew, irep, gama, rx, ry, rz, sx, sy, sz, tx, ty, tz, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, f1x, f1y, f1z, f2x, f2y, f2z, nsigi, sigsp, nsigs, sigi, ixs, x, jhbe, pt, nel, isolnod)
Definition smorth3.F:43
ancmsg
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)
Definition message.F:889
fretitl2
subroutine fretitl2(titr, iasc, l)
Definition freform.F:804