OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i24pen3.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "vect07_c.inc"
#include "scr05_c.inc"
#include "com04_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i24pen3 (x, irect, gapv, cand_e, cand_n, nsv, inacti, itab, tag, iwpene, nsn, irtlm, msegtyp, iwpene0, pmin, gap_n, mvoisn, ixs, ixs10, ixs16, ixs20, penmax, penmin, id, titr, ilev, pen_old, knod2els, nod2els, ipartns, ipen0, icont_i, xfic, nrtm, irtse, is2se)
subroutine ini_st3 (xx, yy, zz, xi, yi, zi, nn, ssc, ttc, ier, alp, xc, yc, zc)
subroutine i24penmax (pen, penmax, etyp, el, ns, ixs, ixs10, ixs16, ixs20, ielim)
subroutine iconnet (irect, ixs, knod2els, nod2els, ixs10, ixs16, ixs20, ns, iconn)
subroutine i24cand (cand_e, cand_n, nsn, irtlm, ii_stok, msegtyp)

Function/Subroutine Documentation

◆ i24cand()

subroutine i24cand ( integer, dimension(*) cand_e,
integer, dimension(*) cand_n,
integer nsn,
integer, dimension(2,*) irtlm,
integer ii_stok,
integer, dimension(*) msegtyp )

Definition at line 731 of file i24pen3.F.

733C
734C-----------------------------------------------
735C I m p l i c i t T y p e s
736C-----------------------------------------------
737#include "implicit_f.inc"
738C-----------------------------------------------
739C C o m m o n B l o c k s
740C-----------------------------------------------
741 INTEGER CAND_E(*),CAND_N(*),NSN,IRTLM(2,*),II_STOK,
742 * MSEGTYP(*)
743C-----------------------------------------------
744C L o c a l V a r i a b l e s
745C-----------------------------------------------
746 INTEGER E, I,ISH
747 .
748C-----------------------------------------------
749C E x t e r n a l F u n c t i o n s
750C-----------------------------------------------
751 ii_stok = 0
752 DO i=1,nsn
753 e = irtlm(1,i)
754 IF (e >0) THEN
755 ii_stok =ii_stok + 1
756 cand_n(ii_stok) = i
757 cand_e(ii_stok) = e
758
759 ish = abs(msegtyp(e))
760 IF (ish > 0)THEN
761 ii_stok =ii_stok + 1
762 cand_n(ii_stok) = i
763 cand_e(ii_stok) = ish
764 ENDIF
765
766 END IF
767 END DO
768C
769 RETURN

◆ i24pen3()

subroutine i24pen3 ( x,
integer, dimension(4,*) irect,
gapv,
integer, dimension(*) cand_e,
integer, dimension(*) cand_n,
integer, dimension(*) nsv,
integer inacti,
integer, dimension(*) itab,
integer, dimension(*) tag,
integer iwpene,
integer nsn,
integer, dimension(2,*) irtlm,
integer, dimension(*) msegtyp,
integer iwpene0,
pmin,
gap_n,
integer, dimension(4,*) mvoisn,
integer, dimension(nixs,*) ixs,
integer, dimension(6,*) ixs10,
integer, dimension(8,*) ixs16,
integer, dimension(12,*) ixs20,
penmax,
penmin,
integer id,
character(len=nchartitle) titr,
integer ilev,
pen_old,
integer, dimension(*) knod2els,
integer, dimension(*) nod2els,
integer, dimension(*) ipartns,
integer ipen0,
integer, dimension(*) icont_i,
xfic,
integer nrtm,
integer, dimension(*) irtse,
integer, dimension(*) is2se )

Definition at line 38 of file i24pen3.F.

46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49#ifndef HYPERMESH_LIB
50 USE message_mod
51#endif
53 USE format_mod , ONLY : fmt_i_3f
54 use element_mod , only :nixs
55C
56C-----------------------------------------------
57C I m p l i c i t T y p e s
58C-----------------------------------------------
59#include "implicit_f.inc"
60C-----------------------------------------------
61C C o m m o n B l o c k s
62C-----------------------------------------------
63#include "units_c.inc"
64#include "vect07_c.inc"
65#include "scr05_c.inc"
66#include "com04_c.inc"
67C-----------------------------------------------
68 INTEGER IWPENE,TAG(*),INACTI,NSV(*),NSN,MSEGTYP(*),IWPENE0,
69 . MVOISN(4,*),ILEV,KNOD2ELS(*),NOD2ELS(*),IPARTNS(*),NRTM
70C---- IWPENE0 : nb of tiny initial pene; IWPENE nb of initial pene
71 my_real gapv(*)
72 INTEGER IRECT(4,*), ITAB(*),CAND_E(*),CAND_N(*),IRTLM(2,*)
73 INTEGER IXS(NIXS,*),IXS10(6,*), IXS16(8,*), IXS20(12,*),ICONT_I(*),
74 . IRTSE(*) ,IS2SE(*)
75 my_real x(3,*),pmin(*),gap_n(12,*),penmax,penmin,pen_old(5,nsn),xfic(3,*)
76C--------GAP_N(1,*) stock temporarily characteristic length for Pen_max
77 INTEGER ID,IPEN0
78 CHARACTER(LEN=NCHARTITLE) :: TITR
79C-----------------------------------------------
80C L o c a l V a r i a b l e s
81C-----------------------------------------------
82 INTEGER II, I, J, K, L, JJ, NJ, IER,NS,IC,I0,IELIM,NI,ICONN,ip,NS1,
83 . IDEL,NN1,NN2,IE
84C REAL
86 . pen, alp,xx(4),yy(4),zz(4),ssc,ttc,dist,dist0,
87 . xi,yi,zi,xc,yc,zc,nn(3),tol,pen0,dpen,norm,maxpen
88C-----------------------------------------------
89C E x t e r n a l F u n c t i o n s
90C-----------------------------------------------
91C--- MVOISN,IRTLM(2,*) is used temporarily for Pen_ini MVOISN(1,*) -> MTYPE(solid),MVOISN(2,*) -> E_id
92C-----MVOISN(3,*) -> part_id, PEN_OLD(1-3,*)->SECONDARY nodal normal, reput PEN_OLD(1-2,*)=0, 3,)=5,)after
93C-----to be consistent with engine for PENMIN
94 IF (iresp==1.AND.penmin<=em06) penmin = two*em06
95 tol = penmin
96 alp = two*em06
97 IF (iresp==1) alp = two*em05
98 DO i=lft,llt
99 l = cand_e(i)
100 ni = cand_n(i)
101 ns = nsv(ni)
102 IF (ns >numnod) THEN
103 ns1 = ns -numnod
104 xi=xfic(1,ns1)
105 yi=xfic(2,ns1)
106 zi=xfic(3,ns1)
107 ELSE
108 xi=x(1,ns)
109 yi=x(2,ns)
110 zi=x(3,ns)
111 END IF
112 DO jj=1,4
113 nj=irect(jj,l)
114 xx(jj)=x(1,nj)
115 yy(jj)=x(2,nj)
116 zz(jj)=x(3,nj)
117 END DO
118C--------
119 CALL ini_st3(xx,yy,zz,xi,yi,zi,nn,ssc,ttc,ier,alp,
120 2 xc,yc,zc)
121 IF(ier==-1)THEN
122#ifndef HYPERMESH_LIB
123 CALL ancmsg(msgid=85,
124 . msgtype=msgerror,
125 . anmode=aninfo,
126 . i1=id,
127 . c1=titr,
128 . i2=itab(ns),
129 . i3=l,
130 . i4=l,
131 . i5=itab(irect(1,l)),
132 . i6=itab(irect(2,l)),
133 . i7=itab(irect(3,l)),
134 . i8=itab(irect(4,l)))
135#endif
136C
137 ELSE IF(ier==1.AND.(msegtyp(l)/=0.AND.msegtyp(l)<=nrtm))THEN
138C shells except coating shells
139C---------outside
140c IF(IPRI>=1)WRITE(IOUT,FMT=FMT_6I)ITAB(NS),L,
141c . (ITAB(IRECT(JJ,L)),JJ=1,4)
142 ELSE
143C------initial penetration case 1) |PEN|<TOL : on, 2) Inacti<0 on,3)Inacti=3,4
144C-------warnning out for 2),3),4)
145C --------MSEGTYP /=0 ->shell ---
146 pen0=nn(1)*(xi-xc)+nn(2)*(yi-yc)+nn(3)*(zi-zc)
147 IF(ier==1) THEN
148 dist = sqrt((xi-xc)*(xi-xc)+(yi-yc)*(yi-yc)+(zi-zc)*(zi-zc))
149 ELSE
150 dist = abs(pen0)
151 END IF
152C--------for exception of elimination---
153 idel = 1
154C----------coating shell is like solid now----
155 IF (msegtyp(l)/=0.AND.msegtyp(l)<=nrtm) THEN
156 pen=gapv(i)-abs(pen0)
157 IF (pen > penmax ) idel = 0
158C-debug sandwish shell : avoid elimination with high thickness outside
159 IF (pen > zero) dist = abs(gapv(i)-pen0)
160C----------give up the the wrong one (normal direction, and)
161 IF (pen0 < zero .OR. pen > penmax) pen=-abs(pen)-tol
162C----------distance after shifted
163 ELSE
164 pen=gapv(i)-pen0
165C------------used only for eliminating wrong contact w/ smaller distance
166 IF(ier==1) pen=-abs(pen)-tol
167 IF (pen > zero .OR. abs(pen) < tol) THEN
168 maxpen = gap_n(1,l)
169 IF (inacti /= 0) maxpen = penmax
170 CALL i24penmax(pen,maxpen ,mvoisn(1,l),mvoisn(2,l),
171 + ns ,ixs, ixs10, ixs16, ixs20 ,
172 + ielim)
173 iconn = 0
174 IF (ns>numnod) THEN
175 CALL i24fic_getn(ns1 ,irtse ,is2se ,ie ,nn1 ,
176 4 nn2 )
177 CALL iconnet(irect(1,l),ixs ,knod2els,nod2els,
178 . ixs10 ,ixs16 ,ixs20 ,nn1 ,iconn )
179 IF (iconn == 0)
180 . CALL iconnet(irect(1,l),ixs ,knod2els,nod2els,
181 . ixs10 ,ixs16 ,ixs20 ,nn2 ,iconn )
182 ELSE
183 CALL iconnet(irect(1,l),ixs ,knod2els,nod2els,
184 . ixs10 ,ixs16 ,ixs20 ,ns ,iconn )
185 END IF
186 IF ((ielim+iconn) > 0) pen = -abs(pen)-tol
187 IF (pen < zero ) idel = 0
188 END IF
189C------Elimine the impact take into account to SECONDARY nodal normal
190 IF (inacti/=0.AND.(pen > zero .OR. abs(pen) < tol).AND.ilev/=3) THEN
191 norm = nn(1)*pen_old(1,ni)+nn(2)*pen_old(2,ni)
192 + +nn(3)*pen_old(3,ni)
193 IF (norm >= zero) THEN
194 pen = -abs(pen)-tol
195c print *,'impact pair eliminated due to N_SECONDARY'
196 idel = 0
197 END IF
198 END IF
199 END IF !(MSEGTYP(L)/=0 ) THEN
200C------Elimine the impact between the same part
201 IF (ipen0==0) THEN
202 IF (inacti/=0.AND.(pen > zero .OR. abs(pen) < tol)) THEN
203 IF (ipartns(ni) == mvoisn(3,l)) THEN
204 pen = -abs(pen)-tol
205 END IF
206 END IF
207 END IF !(IPEN0==0) THEN
208C------don't take into account auto-impact case for elimination
209 IF (ipartns(ni) == mvoisn(3,l)) idel = 0
210C--------exception for SECONDARY shell (test TV)----
211 IF (gapv(i)>zero.AND.(msegtyp(l)==0.OR.msegtyp(l)>nrtm))idel=0
212C--------PMIN() has been changed from PENE to dist excepting (INACTI ==0,1)
213C--------PMIN() = -dist for ABS(PEN) < TOL .OR. PEN<ZERO
214C------------ cas 1 this part is removed in Engine at T=0 for consisting
215 IF(abs(pen) < tol .OR. (pen<zero.AND.idel>0)) THEN
216C---------only used to calculate Dist_min and to eliminate wrong contact (higher)
217 IF (tag(ns)==0) THEN
218 pmin(ni)=-dist
219 tag(ns)=ni
220 ELSE
221 i0=tag(ns)
222 pen0=pmin(i0)
223 IF (dist <abs(pen0)) THEN
224C----------only update dist_min
225 pmin(ni)=-dist
226 tag(ns)=ni
227 IF (pen0 > zero) THEN
228C----------elimine wrong contact
229 irtlm(1,i0)=0
230 irtlm(2,i0)=0
231 pen_old(5,i0)=zero
232 END IF
233 END IF
234 END IF
235 ELSEIF(pen > penmax) THEN
236C----------warning w/o treatment
237#ifndef HYPERMESH_LIB
238 WRITE(iout,1200)pen
239#endif
240 ELSEIF(pen > zero) THEN
241C------Warning anyway-------------
242 IF (tag(ns)==0) iwpene=iwpene+1
243C------------PMIN has been changed from PENE to dist excepting (INACTI ==0,1)
244 IF(inacti ==0 .OR. inacti ==1) THEN
245C------------use IRTLM(2,NI)-> ICONT_I<0 for initial penetration
246 IF (tag(ns)>0) THEN
247 i0=tag(ns)
248 pen0=pmin(i0)
249C----------exclude case of PMIN(I0)<0 : -dist
250 IF (pen < pen0) THEN
251 icont_i(ni)=-l
252 pmin(ni)=pen
253 tag(ns) = ni
254#ifdef HYPERMESH_LIB
255 pen_old(1:3,ni) = nn(1:3)
256#endif
257 ENDIF
258 ELSE
259 icont_i(ni)=-l
260 pmin(ni)=pen
261 tag(ns) = ni
262#ifdef HYPERMESH_LIB
263 pen_old(1:3,ni) = nn(1:3)
264#endif
265 END IF
266 ELSEIF(inacti ==-1) THEN
267C------------multi-cont-> single by overwriting with min_pene
268 IF (tag(ns)>0) THEN
269 i0=tag(ns)
270 pen0=pmin(i0)
271 dist0 = abs(pmin(i0))
272 IF (dist < dist0) THEN
273 irtlm(1,ni)=l
274 irtlm(2,ni)=1
275 pmin(ni)=dist
276 pen_old(5,ni)=pen
277 tag(ns) = ni
278#ifdef HYPERMESH_LIB
279 pen_old(1:3,ni) = nn(1:3)
280#endif
281 ENDIF
282 ELSE
283 irtlm(1,ni)=l
284 irtlm(2,ni)=1
285 pmin(ni)=dist
286 pen_old(5,ni)=pen
287 tag(ns) = ni
288#ifdef HYPERMESH_LIB
289 pen_old(1:3,ni) = nn(1:3)
290#endif
291 END IF
292C--------hide option------
293 ELSEIF(inacti ==3 ) THEN
294 IF (ilev ==3) THEN
295 dpen = pen + tol
296 ELSE
297 dpen = half*(pen + tol)
298 END IF
299C-------change SECONDARY node
300 IF (tag(ns)==0) THEN
301 irtlm(1,ni)=l
302 irtlm(2,ni)=1
303 iwpene=iwpene+1
304 tag(ns)=ni
305#ifndef HYPERMESH_LIB
306 WRITE(iout,1000)pen
307#endif
308 IF (ns >numnod) THEN
309 ns1 = ns -numnod
310 xfic(1,ns1) = xi + dpen*nn(1)
311 xfic(2,ns1) = yi + dpen*nn(2)
312 xfic(3,ns1) = zi + dpen*nn(3)
313#ifndef HYPERMESH_LIB
314 WRITE(iout,fmt=fmt_i_3f)(itab(numnod)+ns1),xfic(1,ns1),xfic(2,ns1),xfic(3,ns1)
315#endif
316 ELSE
317 x(1,ns) = xi + dpen*nn(1)
318 x(2,ns) = yi + dpen*nn(2)
319 x(3,ns) = zi + dpen*nn(3)
320#ifndef HYPERMESH_LIB
321 WRITE(iout,fmt=fmt_i_3f)itab(ns),x(1,ns),x(2,ns),x(3,ns)
322#endif
323 END IF !(NS >NUMNOD) THEN
324 END IF
325 ELSEIF(inacti ==5) THEN
326C------------multi-cont-> single by overwriting with min_pene
327 IF (tag(ns)>0) THEN
328 i0=tag(ns)
329 pen0=pen_old(5,i0)
330 dist0 = abs(pmin(i0))
331 IF (dist < dist0) THEN
332 irtlm(1,ni)=l
333 irtlm(2,ni)=1
334 pen_old(5,ni)=pen
335 pmin(ni)=dist
336 tag(ns) = ni
337#ifdef HYPERMESH_LIB
338 pen_old(1:3,ni) = nn(1:3)
339#endif
340 ENDIF
341 ELSE
342 irtlm(1,ni)=l
343 irtlm(2,ni)=1
344 pen_old(5,ni)=pen
345 pmin(ni)=dist
346 tag(ns) = ni
347#ifdef HYPERMESH_LIB
348 pen_old(1:3,ni) = nn(1:3)
349#endif
350 END IF
351 END IF !IF(INACTI ==0 .OR. INACTI ==1) THEN
352 END IF !(PEN > ZERO) THEN
353 END IF !(IER==-1)THEN
354 END DO !I=LFT,LLT
355C
356 RETURN
357 1000 FORMAT(2x,'** INITIAL PENETRATION =',1pg20.13,
358 . ' CHANGE COORDINATES OF SECONDARY NODE TO:')
359 1100 FORMAT(2x,'** INITIAL PENETRATION =',1pg20.13,
360 . ' CHANGE COORDINATES OF MAIN NODE TO:')
361 1200 FORMAT(2x,'** TOO HIGH INITIAL PENETRATION=, WILL BE IGNORED',
362 . 1pg20.13)
363C
#define my_real
Definition cppsort.cpp:32
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
subroutine i24fic_getn(ns, irtse, is2se, ie, ns1, ns2)
Definition i24surfi.F:1923
initmumps id
integer, parameter nchartitle
subroutine ini_st3(xx, yy, zz, xi, yi, zi, nn, ssc, ttc, ier, alp, xc, yc, zc)
Definition i24pen3.F:372
subroutine i24penmax(pen, penmax, etyp, el, ns, ixs, ixs10, ixs16, ixs20, ielim)
Definition i24pen3.F:596
subroutine iconnet(irect, ixs, knod2els, nod2els, ixs10, ixs16, ixs20, ns, iconn)
Definition i24pen3.F:655
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

◆ i24penmax()

subroutine i24penmax ( pen,
penmax,
integer etyp,
integer el,
integer ns,
integer, dimension(nixs,*) ixs,
integer, dimension(6,*) ixs10,
integer, dimension(8,*) ixs16,
integer, dimension(12,*) ixs20,
integer ielim )

Definition at line 594 of file i24pen3.F.

596 use element_mod , only : nixs
597C-----------------------------------------------
598C I m p l i c i t T y p e s
599C-----------------------------------------------
600#include "implicit_f.inc"
601#include "com04_c.inc"
602C-----------------------------------------------
603C D u m m y A r g u m e n t s
604C-----------------------------------------------
605 INTEGER IXS(NIXS,*),IXS10(6,*), IXS16(8,*), IXS20(12,*)
606 INTEGER ETYP ,EL ,NS,IELIM
607C REAL
608 my_real
609 . pen ,penmax
610C-----------------------------------------------
611C External function
612C-----------------------------------------------
613 LOGICAL INTAB
614 EXTERNAL intab
615C-----------------------------------------------
616C L o c a l V a r i a b l e s
617C-----------------------------------------------
618C REAL
619 INTEGER EL2
620 my_real
621 . s
622C---- Common Add Id--And End --------------------------------------------------------------------
623C--------eliminier some initial penetrations---------
624C------automatic for self_contact
625 ielim=0
626 SELECT CASE (etyp)
627 CASE(1)
628 IF (intab(8,ixs(2,el),ns)) ielim=1
629 CASE(10)
630 el2=el-numels8
631 IF (intab(8,ixs(2,el),ns).OR.intab(6,ixs10(1,el2),ns))
632 + ielim=1
633 CASE(20)
634 el2=el-numels8-numels10
635 IF (intab(8,ixs(2,el),ns).OR.intab(12,ixs20(1,el2),ns))
636 + ielim=1
637 CASE(16)
638 el2=el-numels8-numels10-numels20
639 IF (intab(8,ixs(2,el),ns).OR.intab(8,ixs16(1,el2),ns))
640 + ielim=1
641 END SELECT
642C-------
643 IF (pen >= penmax ) ielim = 1
644
645 RETURN
logical function intab(nic, ic, n)
Definition i24tools.F:95

◆ iconnet()

subroutine iconnet ( integer, dimension(4) irect,
integer, dimension(nixs,*) ixs,
integer, dimension(*) knod2els,
integer, dimension(*) nod2els,
integer, dimension(6,*) ixs10,
integer, dimension(8,*) ixs16,
integer, dimension(12,*) ixs20,
integer ns,
integer iconn )

Definition at line 653 of file i24pen3.F.

655 use element_mod , only : nixs
656C-----------------------------------------------
657C I m p l i c i t T y p e s
658C-----------------------------------------------
659#include "implicit_f.inc"
660C-----------------------------------------------
661C C o m m o n B l o c k s
662C-----------------------------------------------
663#include "com04_c.inc"
664C-----------------------------------------------
665C D u m m y A r g u m e n t s
666C-----------------------------------------------
667C REAL
668 INTEGER IRECT(4), IXS(NIXS,*), KNOD2ELS(*), NOD2ELS(*),
669 . IXS10(6,*), IXS16(8,*), IXS20(12,*),ICONN,NS
670C REAL
671C-----------------------------------------------
672C L o c a l V a r i a b l e s
673C-----------------------------------------------
674 INTEGER N, JJ, II, K, NN, KK, IC, IAD
675C REAL
676C-----------------------------------------------
677C E x t e r n a l F u n c t i o n s
678C-----------------------------------------------
679C
680 iconn = 0
681 IF(numels==0) RETURN
682 DO 230 iad=knod2els(ns)+1,knod2els(ns+1)
683 n = nod2els(iad)
684 IF(n <= numels8)THEN
685 DO jj=1,4
686 ii=irect(jj)
687 DO k=1,8
688 IF(ixs(k+1,n)==ii) iconn = 1
689 ENDDO
690 ENDDO
691 ELSEIF(n <= numels8+numels10)THEN
692 DO jj=1,4
693 ii=irect(jj)
694 DO k=1,8
695 IF(ixs(k+1,n)==ii) iconn = 1
696 ENDDO
697 DO k=1,6
698 IF(ixs10(k,n-numels8)==ii) iconn = 1
699 ENDDO
700 ENDDO
701 ELSEIF(n <= numels8+numels10+numels20)THEN
702 DO jj=1,4
703 ii=irect(jj)
704 DO k=1,8
705 IF(ixs(k+1,n)==ii) iconn = 1
706 ENDDO
707 DO k=1,12
708 IF(ixs20(k,n-numels8-numels10)==ii) iconn = 1
709 ENDDO
710 ENDDO
711 ELSEIF(n <= numels8+numels10+numels20+numels16)THEN
712 DO jj=1,4
713 ii=irect(jj)
714 DO k=1,8
715 IF(ixs(k+1,n)==ii) iconn = 1
716 ENDDO
717 DO k=1,8
718 IF(ixs16(k,n-numels8-numels10-numels20)==ii) iconn = 1
719 ENDDO
720 ENDDO
721 END IF
722 230 CONTINUE
723 RETURN

◆ ini_st3()

subroutine ini_st3 ( xx,
yy,
zz,
xi,
yi,
zi,
nn,
ssc,
ttc,
integer ier,
alp,
xc,
yc,
zc )

Definition at line 370 of file i24pen3.F.

372C
373C-----------------------------------------------
374C I m p l i c i t T y p e s
375C-----------------------------------------------
376#include "implicit_f.inc"
377C-----------------------------------------------
378C D u m m y A r g u m e n t s
379C-----------------------------------------------
380 INTEGER IER
381 my_real
382 . xx(4),yy(4),zz(4),nn(3), ssc, ttc, alp,xi,yi,zi,xc,yc,zc
383C-----------------------------------------------
384C L o c a l V a r i a b l e s
385C-----------------------------------------------
386 my_real
387 . h(4), x0, y0, z0, xl1, xl2, xl3, xl4, yy1, yy2, yy3, yy4,
388 . zz1, zz2, zz3, zz4, xi1, xi2, xi3, xi4, yi1, yi2, yi3, yi4,
389 . zi1, zi2, zi3, zi4, xn1, yn1, zn1, xn2, yn2, zn2, xn3, yn3,
390 . zn3, xn4, yn4, zn4, an, area, a12, a23, a34, a41, b12, b23,
391 . b34, b41, ab1, ab2, tp, tm, sp, sm, x1,x2,x3,x4,
392 . y1,y2,y3,y4,z1,z2,z3,z4,n1,n2,n3,la,lb,lc,lbs,lcs,tt1,ss1
393C
394 x1=xx(1)
395 x2=xx(2)
396 x3=xx(3)
397 x4=xx(4)
398 y1=yy(1)
399 y2=yy(2)
400 y3=yy(3)
401 y4=yy(4)
402 z1=zz(1)
403 z2=zz(2)
404 z3=zz(3)
405 z4=zz(4)
406C
407 x0 = fourth*(x1+x2+x3+x4)
408 y0 = fourth*(y1+y2+y3+y4)
409 z0 = fourth*(z1+z2+z3+z4)
410C
411 xl1 = x1-x0
412 xl2 = x2-x0
413 xl3 = x3-x0
414 xl4 = x4-x0
415 yy1 = y1-y0
416 yy2 = y2-y0
417 yy3 = y3-y0
418 yy4 = y4-y0
419 zz1 = z1-z0
420 zz2 = z2-z0
421 zz3 = z3-z0
422 zz4 = z4-z0
423C
424 xi1 = x1-xi
425 xi2 = x2-xi
426 xi3 = x3-xi
427 xi4 = x4-xi
428 yi1 = y1-yi
429 yi2 = y2-yi
430 yi3 = y3-yi
431 yi4 = y4-yi
432 zi1 = z1-zi
433 zi2 = z2-zi
434 zi3 = z3-zi
435 zi4 = z4-zi
436C
437 xn1 = yy1*zz2 - yy2*zz1
438 yn1 = zz1*xl2 - zz2*xl1
439 zn1 = xl1*yy2 - xl2*yy1
440 n1=xn1
441 n2=yn1
442 n3=zn1
443C
444 xn2 = yy2*zz3 - yy3*zz2
445 yn2 = zz2*xl3 - zz3*xl2
446 zn2 = xl2*yy3 - xl3*yy2
447 n1=n1+xn2
448 n2=n2+yn2
449 n3=n3+zn2
450C
451 xn3 = yy3*zz4 - yy4*zz3
452 yn3 = zz3*xl4 - zz4*xl3
453 zn3 = xl3*yy4 - xl4*yy3
454 n1=n1+xn3
455 n2=n2+yn3
456 n3=n3+zn3
457C
458 xn4 = yy4*zz1 - yy1*zz4
459 yn4 = zz4*xl1 - zz1*xl4
460 zn4 = xl4*yy1 - xl1*yy4
461 n1=n1+xn4
462 n2=n2+yn4
463 n3=n3+zn4
464C
465 an= max(em20,sqrt(n1*n1+n2*n2+n3*n3))
466 n1=n1/an
467 n2=n2/an
468 n3=n3/an
469 nn(1)=n1
470 nn(2)=n2
471 nn(3)=n3
472 IF(an<=em19) THEN
473 ier=-1
474 RETURN
475 ENDIF
476 area=half*an
477C
478 a12=(n1*xn1+n2*yn1+n3*zn1)
479 a23=(n1*xn2+n2*yn2+n3*zn2)
480 a34=(n1*xn3+n2*yn3+n3*zn3)
481 a41=(n1*xn4+n2*yn4+n3*zn4)
482C
483 xn1 = yi1*zi2 - yi2*zi1
484 yn1 = zi1*xi2 - zi2*xi1
485 zn1 = xi1*yi2 - xi2*yi1
486 b12=(n1*xn1+n2*yn1+n3*zn1)
487C
488 xn2 = yi2*zi3 - yi3*zi2
489 yn2 = zi2*xi3 - zi3*xi2
490 zn2 = xi2*yi3 - xi3*yi2
491 b23=(n1*xn2+n2*yn2+n3*zn2)
492C
493 xn3 = yi3*zi4 - yi4*zi3
494 yn3 = zi3*xi4 - zi4*xi3
495 zn3 = xi3*yi4 - xi4*yi3
496 b34=(n1*xn3+n2*yn3+n3*zn3)
497C
498 xn4 = yi4*zi1 - yi1*zi4
499 yn4 = zi4*xi1 - zi1*xi4
500 zn4 = xi4*yi1 - xi1*yi4
501 b41=(n1*xn4+n2*yn4+n3*zn4)
502C
503 ab1=a23*b41
504 ab2=b23*a41
505C
506 IF(abs(ab1+ab2)/area>em10)THEN
507 ssc=(ab1-ab2)/(ab1+ab2)
508 ELSE
509 ssc=zero
510 ENDIF
511 IF(abs(a34/area)>em10)THEN
512 ab1=b12*a34
513 ab2=b34*a12
514 IF(abs(ab1+ab2)/area>em10)THEN
515 ttc=(ab1-ab2)/(ab1+ab2)
516 ELSE
517 ttc=zero
518 END IF
519 ELSE
520 ttc=b12/a12-one
521 IF(b23<=zero.AND.b41<=zero)THEN
522 IF(-b23/a12<=alp.AND.-b41/a12<=alp)ssc=zero
523 ELSEIF(b23<=zero)THEN
524 IF(-b23/a12<=alp) THEN
525 ssc=one
526 ELSE
527 ssc=two
528 END IF
529 ELSEIF(b41<=zero)THEN
530 IF(-b41/a12<=alp) THEN
531 ssc=-one
532 ELSE
533 ssc=-two
534 END IF
535 ENDIF
536 ENDIF
537C-------------out of seg
538 IF(abs(ssc)>one+alp.OR.abs(ttc)>one+alp) THEN
539 ier=1
540C------case tria re-compute
541 IF (a34==zero.AND.ttc< one) THEN
542 lb=fourth*(one - ttc)*(one - ssc)
543 lc=fourth*(one - ttc)*(one + ssc)
544 la = one - lb - lc
545 IF (la>=zero) THEN
546 lb= min(one,max(zero,lb))
547 lc= min(one,max(zero,lc))
548 ELSEIF(lc>lb.AND.lc >= one) THEN
549 lc = one
550 lb = zero
551 ELSEIF(lb >= one) THEN
552 lc = zero
553 lb = one
554 ELSE
555 lbs = half*(one+lb-lc)
556 lcs = half*(one-lb+lc)
557 lb= min(one,max(zero,lbs))
558 lc= min(one,max(zero,lcs))
559 ENDIF
560 ssc= (lc-lb)/(lc+lb)
561 ttc= one - two*lb - two*lc
562 END IF
563 IF(abs(ssc)>one)ssc=ssc/abs(ssc)
564 IF(abs(ttc)>one)ttc=ttc/abs(ttc)
565 ELSE
566 ier=0
567 IF(abs(ssc)>one)ssc=ssc/abs(ssc)
568 IF(abs(ttc)>one)ttc=ttc/abs(ttc)
569 ENDIF
570C
571 tp=fourth*(one+ttc)
572 tm=fourth*(one-ttc)
573C
574 sp=one+ssc
575 sm=one-ssc
576 h(1)=tm*sm
577 h(2)=tm*sp
578 h(3)=tp*sp
579 h(4)=tp*sm
580C
581 xc =h(1)*x1+h(2)*x2+h(3)*x3+h(4)*x4
582 yc =h(1)*y1+h(2)*y2+h(3)*y3+h(4)*y4
583 zc =h(1)*z1+h(2)*z2+h(3)*z3+h(4)*z4
584 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
#define min(a, b)
Definition macros.h:20
#define max(a, b)
Definition macros.h:21