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