OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cinit3.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "vect01_c.inc"
#include "param_c.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "scr03_c.inc"
#include "scr17_c.inc"
#include "scry_c.inc"
#include "com_xfem1.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine cinit3 (elbuf_str, ixc, pm, x, geo, xmas, in, nvc, dtelem, igrsh4n, xrefc, nel, ithk, ihbe, igrsh3n, thk, isigsh, sigsh, stifn, stifr, partsav, v, ipart, msc, inc, skew, iparg, i8mi, nsigsh, igeo, iuser, etnod, nshnod, stc, ptshel, ipm, bufmat, sh4tree, mcp, mcps, temp, cpt_eltens, part_area, itagn, itage, ixfem, npf, tf, xfem_str, isubstack, stack, rnoise, drape, sh4ang, iddlevel, geo_stack, igeo_stack, strc, perturb, iyldini, ele_area, ng, group_param, nloc_dmg, idrape, drapeg, mat_param, fail_fractal, fail_brokmann, glob_therm)
subroutine csms11_ini (jft, jlt, ixc, x, x2s, y2s, x3s, y3s, x4s, y4s)

Function/Subroutine Documentation

◆ cinit3()

subroutine cinit3 ( type(elbuf_struct_), target elbuf_str,
integer, dimension(nixc,*) ixc,
pm,
x,
geo,
xmas,
in,
integer nvc,
dtelem,
type (group_), dimension(ngrshel) igrsh4n,
xrefc,
integer nel,
integer ithk,
integer ihbe,
type (group_), dimension(ngrsh3n) igrsh3n,
thk,
integer isigsh,
sigsh,
stifn,
stifr,
partsav,
v,
integer, dimension(*) ipart,
msc,
inc,
skew,
integer, dimension(*) iparg,
integer *8, dimension(6,*) i8mi,
integer nsigsh,
integer, dimension(npropgi,*) igeo,
integer iuser,
etnod,
integer, dimension(*) nshnod,
stc,
integer, dimension(*) ptshel,
integer, dimension(npropmi,*) ipm,
bufmat,
integer, dimension(*) sh4tree,
mcp,
mcps,
temp,
integer cpt_eltens,
part_area,
integer, dimension(*) itagn,
integer, dimension(*) itage,
integer ixfem,
integer, dimension(*) npf,
tf,
type(elbuf_struct_), dimension(ngroup,*), target xfem_str,
integer isubstack,
type (stack_ply) stack,
rnoise,
type (drape_), dimension(numelc_drape + numeltg_drape) drape,
sh4ang,
integer, intent(in) iddlevel,
geo_stack,
integer, dimension(*) igeo_stack,
strc,
integer, dimension(nperturb) perturb,
integer iyldini,
ele_area,
integer ng,
type (group_param_) group_param,
type (nlocal_str_) nloc_dmg,
integer idrape,
type (drapeg_) drapeg,
type (matparam_struct_), dimension(nummat), intent(inout) mat_param,
type (fail_fractal_), intent(in) fail_fractal,
type (fail_brokmann_), intent(in) fail_brokmann,
type (glob_therm_), intent(in) glob_therm )

Definition at line 66 of file cinit3.F.

81C-----------------------------------------------
82C M o d u l e s
83C-----------------------------------------------
84 USE my_alloc_mod
85 USE elbufdef_mod
86 USE message_mod
87 USE stack_mod
88 USE groupdef_mod
89 USE group_param_mod
90 USE nlocal_reg_mod
91 USE drape_mod
92 USE matparam_def_mod
93 USE random_walk_def_mod
94 USE fractal_dmg_init_mod
95 USE names_and_titles_mod , ONLY : nchartitle
96 use brokmann_random_def_mod
97 use glob_therm_mod
98 use initemp_shell_mod
99 use element_mod , only : nixc
100 use law42c_ini_mod
101C-----------------------------------------------
102C I m p l i c i t T y p e s
103C-----------------------------------------------
104#include "implicit_f.inc"
105C-----------------------------------------------
106C G l o b a l P a r a m e t e r s
107C-----------------------------------------------
108#include "mvsiz_p.inc"
109C-----------------------------------------------
110C C o m m o n B l o c k s
111C-----------------------------------------------
112#include "vect01_c.inc"
113#include "param_c.inc"
114#include "com01_c.inc"
115#include "com04_c.inc"
116#include "scr03_c.inc"
117#include "scr17_c.inc"
118#include "scry_c.inc"
119#include "com_xfem1.inc"
120C-----------------------------------------------
121C D u m m y A r g u m e n t s
122C-----------------------------------------------
123 INTEGER NVC,NEL,ITHK,IHBE,ISIGSH,IXFEM,NSIGSH,IUSER,IYLDINI
124 INTEGER IXC(NIXC,*),IPART(*),IPARG(*),IGEO(NPROPGI,*), NSHNOD(*),
125 . PTSHEL(*),IPM(NPROPMI,*), SH4TREE(*),ITAGN(*),ITAGE(*),NPF(*),
126 . ISUBSTACK,IGEO_STACK(*),PERTURB(NPERTURB),NG,IDRAPE
127 INTEGER *8 I8MI(6,*)
128 INTEGER ,INTENT(IN) :: IDDLEVEL
129 my_real
130 . pm(npropm,*), x(3,*), geo(npropg,*), xmas(*), in(*),
131 . dtelem(*), xrefc(4,3,*),thk(*), sigsh(nsigsh,*),
132 . stifn(*),stifr(*),partsav(20,*), v(*) ,msc(*) ,inc(*),
133 . skew(lskew,*), etnod(*), stc(*),bufmat(*),mcp(*),mcps(*),
134 . temp(*),part_area(*),tf(*),rnoise(*),
135 . sh4ang(*),geo_stack(*),strc(*),ele_area(*)
136 TYPE(ELBUF_STRUCT_), TARGET :: ELBUF_STR
137 TYPE(ELBUF_STRUCT_), TARGET ,DIMENSION(NGROUP,*):: XFEM_STR
138 ! when XFEM is ON, XFEM_STR's dimension = NGROUP,NXEL
139 TYPE (STACK_PLY) :: STACK
140 TYPE (GROUP_PARAM_) :: GROUP_PARAM
141 TYPE (NLOCAL_STR_) :: NLOC_DMG
142 TYPE (DRAPE_) :: DRAPE(NUMELC_DRAPE + NUMELTG_DRAPE)
143 TYPE (DRAPEG_) :: DRAPEG
144 TYPE (MATPARAM_STRUCT_) ,DIMENSION(NUMMAT) ,INTENT(INOUT) :: MAT_PARAM
145 TYPE (FAIL_FRACTAL_) ,INTENT(IN) :: FAIL_FRACTAL
146 TYPE (FAIL_BROKMANN_) ,INTENT(IN) :: FAIL_BROKMANN
147 TYPE (glob_therm_) ,intent(in) :: glob_therm
148C-----------------------------------------------
149 TYPE (GROUP_) , DIMENSION(NGRSHEL) :: IGRSH4N
150 TYPE (GROUP_) , DIMENSION(NGRSH3N) :: IGRSH3N
151C-----------------------------------------------
152C L o c a l V a r i a b l e s
153C-----------------------------------------------
154 CHARACTER(LEN=NCHARTITLE) :: TITR,TITR1
155 INTEGER I,J,N,IP,NDEPAR,IGTYP,NUVAR,ID,NLAY,II,IREP,IPROP,NUPARAM,
156 . IL,IR,IS,IT,CPT_ELTENS,IUN,NPTR,NPTS,NPTT,IXEL,ILAW,IMAT,IFAIL,
157 . IGMAT,JJ(9),NPT_ALL,MPT,LAYNPT_MAX,LAY_MAX
158 INTEGER, DIMENSION(MVSIZ) :: IX1,IX2,IX3,IX4,IORTHLOC,MAT,PID,NGL
159 my_real :: bid
160 my_real, DIMENSION(MVSIZ) :: px1g,px2g,py1g,py2g,area,aldt,
161 . dt,vx,vy,vz,
162 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
163 . e1x,e2x,e3x,e1y,e2y,e3y,e1z,e2z,e3z,
164 . x2s,y2s,x3s,y3s,x4s,y4s,
165 . x2l,x3l,x4l,y2l,y3l,y4l
166 my_real, DIMENSION(NEL) :: tempel
167 my_real, DIMENSION(:) ,POINTER :: uvar,dir1,dir2
168 my_real, ALLOCATABLE, DIMENSION(:) :: dir_a,dir_b
169 my_real, DIMENSION(:), ALLOCATABLE :: phi1,phi2,coor1,coor2,coor3,coor4
170 parameter(laynpt_max = 10)
171 parameter(lay_max = 100)
172 INTEGER, DIMENSION(:),ALLOCATABLE :: MATLY
173 my_real, DIMENSION(:,:),ALLOCATABLE :: posly
174
175 TYPE(G_BUFEL_) ,POINTER :: GBUF
176 TYPE(L_BUFEL_) ,POINTER :: LBUF
177C=======================================================================
178 CALL my_alloc(matly,mvsiz*lay_max)
179 CALL my_alloc(posly,mvsiz,lay_max*laynpt_max)
180 gbuf => elbuf_str%GBUF
181c
182 imat = ixc(1,1+nft) ! mat N
183 iprop = ixc(nixc-1,1+nft) ! property N
184 igtyp = igeo(11,iprop)
185 id = igeo(1,iprop)
186 igmat = igeo(98,ixc(6,1+nft))
187 irep = iparg(35)
188 ifail = iparg(43)
189c
190 CALL fretitl2(titr,igeo(npropgi-ltitr+1,iprop),ltitr)
191 vx = zero
192 vy = zero
193 vz = zero
194 iorthloc = 0
195 bid = zero
196C
197 iun = 1
198 ir = 1
199 is = 1
200 nlay = elbuf_str%NLAY
201 nxel = elbuf_str%NXEL
202 nptt = elbuf_str%NPTT
203 idrape = elbuf_str%IDRAPE
204 npt_all = 0
205 DO il=1,nlay
206 npt_all = npt_all + elbuf_str%BUFLY(il)%NPTT
207 ENDDO
208 mpt = max(1,npt_all)
209 IF(npt_all == 0 ) npt_all = nlay
210C----- NPT=0 for some cases
211 IF((igtyp == 51 .OR. igtyp == 52) .AND. idrape > 0) THEN
212 ALLOCATE(phi1(mvsiz*npt_all))
213 ALLOCATE(phi2(nvsiz*npt_all))
214 ALLOCATE(dir_a(npt_all*nel*2))
215 ALLOCATE(dir_b(npt_all*nel*2))
216 phi1 = zero
217 phi2 = zero
218 dir_a = zero
219 dir_b = zero
220 ALLOCATE(coor1(npt_all*mvsiz))
221 ALLOCATE(coor2(npt_all*mvsiz))
222 ALLOCATE(coor3(npt_all*mvsiz))
223 ALLOCATE(coor4(npt_all*mvsiz))
224 coor1 = zero
225 coor2 = zero
226 coor3 = zero
227 coor4 = zero
228 ELSE
229 ALLOCATE(phi1(nlay*mvsiz))
230 ALLOCATE(phi2(nlay*mvsiz))
231 ALLOCATE(dir_a(nlay*nel*2))
232 ALLOCATE(dir_b(nlay*nel*2))
233 phi1 = zero
234 phi2 = zero
235 dir_a = zero
236 dir_b = zero
237 ALLOCATE(coor1(nlay*mvsiz))
238 ALLOCATE(coor2(nlay*mvsiz))
239 ALLOCATE(coor3(nlay*mvsiz))
240 ALLOCATE(coor4(nlay*mvsiz))
241 coor1 = zero
242 coor2 = zero
243 coor3 = zero
244 coor4 = zero
245 npt_all = nlay
246 ENDIF
247C
248 IF (iparg(6) == 0.OR.npt==0) mpt=0
249C
250 DO j=1,9
251 jj(j) = nel*(j-1)
252 ENDDO
253C
254 DO i=lft,llt
255 n = i+nft
256 mat(i) = ixc(1,n)
257 pid(i) = ixc(6,n)
258 ENDDO
259C
260 IF (ixfem > 0) THEN
261 DO i=lft,llt
262 n = i+nft
263 itagn(ixc(2,n)) =1
264 itagn(ixc(3,n)) =1
265 itagn(ixc(4,n)) =1
266 itagn(ixc(5,n)) =1
267 itage(n) = 1
268 ENDDO
269 ENDIF
270C
271 CALL ccoori(x,xrefc(1,1,nft+1),ixc(1,nft+1),
272 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
273 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
274 . ix1 ,ix2 ,ix3 ,ix4 ,ngl )
275c
276 CALL cveok3(nvc,4,ix1,ix2,ix3,ix4)
277c--------------------------------------
278 CALL ceveci(lft,llt ,area,
279 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
280 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
281 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
282c
283C------------
284C Tags total area of the part (needed in /ADMAS for shells)
285C------------
286 IF ((imasadd > 0).OR.(nloc_dmg%IMOD > 0)) THEN
287 DO i=lft,llt
288 ip = ipart(i+nft)
289C PART_AREA(IP) = PART_AREA(IP) + AREA(I)
290 ele_area(i+nft) = area(i)
291 IF (gbuf%G_AREA > 0) gbuf%AREA(i) = area(i)
292 ENDDO
293 ENDIF
294!
295! Initialize element temperature from /initemp
296!
297 IF (jthe == 0 .and. glob_therm%NINTEMP > 0) THEN
298 CALL initemp_shell(elbuf_str,temp,nel,numnod,numelc,4,nixc,ixc)
299 END IF
300C------------
301 CALL cinmas(x ,xrefc(1,1,nft+1),ixc ,geo ,pm,
302 . xmas ,in ,thk ,ihbe ,partsav,
303 . v ,ipart(nft+1) ,msc(nft+1),inc(nft+1) ,area ,
304 . i8mi ,igeo ,etnod ,imat ,iprop ,
305 . nshnod ,stc(nft+1) ,sh4tree ,mcp ,mcps(nft+1),
306 . temp ,bid ,bid ,bid ,bid,
307 . bid ,bid ,isubstack ,nlay ,elbuf_str,
308 . stack ,gbuf%THK_I ,rnoise ,drape ,glob_therm%NINTEMP,
309 . perturb,ix1 ,ix2 ,ix3 ,ix4 ,
310 . idrape ,drapeg%INDX)
311C-----------------------------------------------
312 CALL cderii(px1g,px2g,py1g,py2g,
313 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
314 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
315 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z ,
316 . x2l ,x3l ,x4l ,y2l ,y3l ,y4l )
317 CALL cdleni(pm ,geo ,stifn ,stifr ,ixc(1,nft+1),
318 . px1g ,px2g ,py1g ,py2g ,thk ,
319 . igeo ,dt ,sh4tree ,aldt ,bufmat ,
320 . ipm ,nlay ,stack%PM,isubstack,strc(nft+1),
321 . area ,imat ,iprop ,
322 . x2l ,x3l ,x4l ,y2l ,y3l ,y4l ,
323 . stack%IGEO,group_param)
324 CALL c1buf3(geo,gbuf%THK,gbuf%OFF,thk,ksh4tree,sh4tree)
325C-----------------------------------------------
326 IF (ixfem > 0) THEN
327 DO ixel=1,nxel
328 DO i=lft,llt
329 xfem_str(ng,ixel)%GBUF%THK(i) = thk(i)
330 xfem_str(ng,ixel)%GBUF%OFF(i) = -one
331 END DO
332 ENDDO
333 ENDIF
334C-----------------------------------------------------------------------
335c PHI, COOR used only with dimension(NLAY,MVSIZ); !! now corrected
336 CALL corthini(
337 . lft ,llt ,nft ,nlay ,numelc ,
338 . nsigsh ,nixc ,ixc(1,nft+1),igeo ,geo ,
339 . skew ,sigsh ,ptshel ,phi1 ,phi2 ,
340 . vx ,vy ,vz ,coor1 ,coor2 ,
341 . coor3 ,coor4 ,iorthloc ,isubstack ,stack ,
342 . irep ,elbuf_str ,drape ,sh4ang(nft+1),x ,
343 . geo_stack ,e3x ,e3y ,e3z ,
344 . gbuf%BETAORTH,x1 ,x2 ,y1 ,y2 ,
345 . z1 ,z2 ,nel ,gbuf%G_ADD_NODE,gbuf%ADD_NODE,
346 . npt_all ,idrape ,drapeg%INDX)
347c---
348C-----
349 IF(igtyp == 51 .OR. igtyp == 52 .AND. igmat > 0) THEN
350c
351 CALL corthdir(elbuf_str,
352 . igeo ,geo ,vx ,vy ,vz ,
353 . phi1 ,phi2 ,coor1 ,coor2 ,coor3 ,
354 . coor4 ,iorthloc ,nlay ,irep ,isubstack,
355 . stack ,geo_stack ,igeo_stack ,ir ,is ,
356 . nel ,imat ,iprop ,
357 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
358 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
359 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z ,
360 . npt_all ,idrape)
361c
362 ELSEIF (mtn == 27) THEN
363 CALL cm27in3(elbuf_str,
364 . geo ,igeo ,pm ,ipm ,ixc(1,1+nft) ,nixc,
365 . nlay,ir ,is ,imat )
366 ELSEIF (mtn == 35) THEN
367 nptr = elbuf_str%NPTR
368 npts = elbuf_str%NPTS
369 nptt = elbuf_str%NPTT
370 CALL cm35in3(elbuf_str,thk,area,nel,nlay,
371 . nptr,npts,nptt,igtyp)
372 ELSEIF (mtn==15 .or. mtn==19 .or. mtn==25 .or. mtn>=28) THEN
373 IF (mtn == 19 .AND. igtyp /= 9) THEN
374 CALL ancmsg(msgid=5,
375 . anmode=aninfo,
376 . msgtype=msgerror,
377 . i1=igeo(1,ixc(nixc-1,nft+1)))
378 ENDIF
379c
380 CALL corthdir(elbuf_str,
381 . igeo ,geo ,vx ,vy ,vz ,
382 . phi1 ,phi2 ,coor1 ,coor2 ,coor3 ,
383 . coor4 ,iorthloc ,nlay ,irep ,isubstack,
384 . stack ,geo_stack ,igeo_stack ,ir ,is ,
385 . nel ,imat ,iprop ,
386 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
387 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
388 . e1x, e2x, e3x, e1y, e2y, e3y ,e1z, e2z, e3z ,
389 . npt_all ,idrape)
390 ENDIF
391c
392 IF ((mtn == 58 .or. mtn == 158) .AND.
393 . igtyp /= 16 .AND. igtyp /= 51 .AND. igtyp /= 52) THEN
394 CALL ancmsg(msgid=658,
395 . msgtype=msgerror,
396 . anmode=aninfo_blind_1,
397 . i1=id,
398 . c1=titr,
399 . i2=mtn,
400 . i3=igtyp)
401 ELSEIF (mtn == 58 .or. mtn == 158 .OR. igtyp == 51 .OR. igtyp == 52) THEN
402C
403 IF (idrape == 0) THEN
404 DO il = 1,nlay
405 nptt = elbuf_str%BUFLY(il)%NPTT
406 imat = elbuf_str%BUFLY(il)%IMAT
407 ilaw = elbuf_str%BUFLY(il)%ILAW
408 nuvar = elbuf_str%BUFLY(il)%NVAR_MAT
409 dir1 => elbuf_str%BUFLY(il)%DIRA
410 dir2 => elbuf_str%BUFLY(il)%DIRB
411 nuparam = mat_param(imat)%NUPARAM
412!
413 IF (ilaw == 58) THEN
414 DO it=1,nptt
415 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
416 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
417 CALL cm58in3(
418 . irep ,dir1 ,dir2 ,mat_param(imat)%UPARAM,
419 . uvar ,aldt ,nel ,nuvar ,lbuf%ANG ,
420 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
421 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
422 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
423 ENDDO
424 ELSE IF (ilaw == 158) THEN
425 DO it=1,nptt
426 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
427 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
428 CALL law158_init(dir1 ,dir2 ,
429 . uvar ,aldt ,nel ,nuvar ,lbuf%ANG ,
430 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
431 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
432 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
433 ENDDO
434 ENDIF ! IF (ILAW == 58) THEN
435 ENDDO ! DO IL = 1,NLAY
436 ELSE
437 DO il = 1,nlay
438 nptt = elbuf_str%BUFLY(il)%NPTT
439 imat = elbuf_str%BUFLY(il)%IMAT
440 ilaw = elbuf_str%BUFLY(il)%ILAW
441 nuvar = elbuf_str%BUFLY(il)%NVAR_MAT
442 nuparam = mat_param(imat)%NUPARAM
443!
444 IF (ilaw == 58) THEN
445 DO it=1,nptt
446 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
447 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
448 dir1 => elbuf_str%BUFLY(il)%LBUF_DIR(it)%DIRA
449 dir2 => elbuf_str%BUFLY(il)%LBUF_DIR(it)%DIRB
450 CALL cm58in3(
451 . irep ,dir1 ,dir2 ,mat_param(imat)%UPARAM,
452 . uvar ,aldt ,nel ,nuvar ,lbuf%ANG ,
453 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
454 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
455 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
456 ENDDO
457 ELSE IF (ilaw == 158) THEN
458 DO it=1,nptt
459 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
460 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
461 dir1 => elbuf_str%BUFLY(il)%LBUF_DIR(it)%DIRA
462 dir2 => elbuf_str%BUFLY(il)%LBUF_DIR(it)%DIRB
463 CALL law158_init(dir1 ,dir2 ,
464 . uvar ,aldt ,nel ,nuvar ,lbuf%ANG ,
465 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
466 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
467 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
468 ENDDO
469 ENDIF ! IF (ILAW == 58) THEN
470 ENDDO ! DO IL = 1,NLAY
471 ENDIF
472 END IF !(IGTYP == 51 .OR. IGTYP == 52 .AND. IGMAT > 0) THEN
473 IF (mtn == 42 .OR. mtn == 69 .OR. igtyp == 51 .OR. igtyp == 52) THEN
474 DO il = 1,nlay
475 nptt = elbuf_str%BUFLY(il)%NPTT
476 ilaw = elbuf_str%BUFLY(il)%ILAW
477 nuvar = elbuf_str%BUFLY(il)%NVAR_MAT
478!
479 IF (ilaw == 42 .OR. mtn == 69) THEN
480 DO it=1,nptt
481 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
482 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
483 CALL law42c_ini(nuvar,uvar,llt)
484 ENDDO
485 ENDIF ! IF (ILAW == 42) THEN
486 ENDDO ! DO IL = 1,NLAY
487 ENDIF
488C-----------------------------------------------------------------------
489C calculation of initial constraints
490C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
491 IF (isigsh/=0 .OR. ithkshel == 2) THEN
492C to be checked for IGTYP = 51 : ok
493 IF (mpt > 0) THEN
494 CALL layini1(
495 . elbuf_str ,lft ,llt ,geo ,igeo ,
496 . mat ,pid ,matly ,posly ,igtyp ,
497 . nlay ,mpt ,isubstack ,stack ,drape ,
498 . nft ,gbuf%THK ,nel ,idrape ,scdrape ,
499 . drapeg%INDX)
500 CALL corth3(elbuf_str,dir_a ,dir_b ,lft ,llt ,
501 . nlay ,irep ,nel ,
502 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
503 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
504 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z ,
505 . idrape, igtyp )
506 END IF
507
508 CALL csigini(elbuf_str ,
509 1 lft ,llt ,nft ,mpt ,istrain ,
510 2 gbuf%THK,gbuf%EINT,gbuf%STRA,gbuf%HOURG,gbuf%PLA ,
511 3 gbuf%FOR,gbuf%MOM ,sigsh ,nlay ,gbuf%G_HOURG,
512 4 numelc ,ixc ,nixc ,nsigsh ,numshel ,
513 5 ptshel ,igeo ,thk ,nel ,e1x ,
514 6 e2x ,e3x ,e1y ,e2y ,e3y ,
515 7 e1z ,e2z ,e3z ,isigsh ,dir_a ,
516 8 dir_b ,posly ,igtyp )
517 IF (mpt == 0) THEN
518 DO i=lft,llt
519 gbuf%FOR_G(i+jj(1:5))=gbuf%FOR(i+jj(1:5))
520 END DO
521 END IF
522 ELSEIF ( ithkshel == 1 ) THEN
523 CALL thickini(lft ,llt ,nft ,ptshel,numelc,
524 2 gbuf%THK,thk ,ixc ,nixc ,nsigsh,
525 3 sigsh )
526 ENDIF
527C-----------------------------------------------------------------------
528C CALCUL DES DEFORMATIONS INITIALES (MEMBRANE)
529C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
530 IF (istrain == 1 .AND. nxref > 0) THEN
531 uvar => elbuf_str%BUFLY(1)%MAT(1,1,1)%VAR
532 imat = elbuf_str%BUFLY(1)%IMAT
533 CALL cepsini(elbuf_str ,mat_param(imat),
534 . lft ,llt ,ismstr ,mtn ,ithk ,
535 . pm ,geo ,ixc(1,nft+1),x ,xrefc(1,1,nft+1),
536 . gbuf%FOR ,gbuf%THK ,gbuf%EINT ,gbuf%STRA ,nlay ,
537 . px1g ,px2g ,py1g ,py2g ,x2s ,
538 . y2s ,x3s ,y3s ,x4s ,y4s ,
539 . uvar ,ipm ,igeo ,imat ,
540 . skew ,nel ,dir_a ,dir_b ,gbuf%SIGI,
541 . npf ,tf ,irep )
542c
543 CALL cepschk(lft, llt,nft, pm, geo,ixc(1,nft+1),gbuf%STRA,thk,
544 . nel,cpt_eltens)
545 IF (ismstr == 1 .AND. mtn==19) iparg(9)=11
546c
547 ELSEIF (ismstr == 11 .OR.(ismstr==1 .AND. mtn==19)) THEN
548C to be checked for IGTYP = 51
549 CALL csms11_ini(lft ,llt ,ixc(1,nft+1),x ,
550 . x2s ,y2s ,x3s ,y3s ,x4s ,y4s )
551 ENDIF
552c----------------------------
553 IF (ismstr == 10 ) THEN
554 DO i=lft,llt
555 ii = nft + i
556 elbuf_str%GBUF%SMSTR(jj(1)+i) = x(1,ixc(3,ii))-x(1,ixc(2,ii))
557 elbuf_str%GBUF%SMSTR(jj(2)+i) = x(2,ixc(3,ii))-x(2,ixc(2,ii))
558 elbuf_str%GBUF%SMSTR(jj(3)+i) = x(3,ixc(3,ii))-x(3,ixc(2,ii))
559 elbuf_str%GBUF%SMSTR(jj(4)+i) = x(1,ixc(4,ii))-x(1,ixc(2,ii))
560 elbuf_str%GBUF%SMSTR(jj(5)+i) = x(2,ixc(4,ii))-x(2,ixc(2,ii))
561 elbuf_str%GBUF%SMSTR(jj(6)+i) = x(3,ixc(4,ii))-x(3,ixc(2,ii))
562 elbuf_str%GBUF%SMSTR(jj(7)+i) = x(1,ixc(5,ii))-x(1,ixc(2,ii))
563 elbuf_str%GBUF%SMSTR(jj(8)+i) = x(2,ixc(5,ii))-x(2,ixc(2,ii))
564 elbuf_str%GBUF%SMSTR(jj(9)+i) = x(3,ixc(5,ii))-x(3,ixc(2,ii))
565 ENDDO
566 ELSEIF (ismstr == 11 .OR.(ismstr==1 .AND. mtn==19)) THEN
567 DO i=lft,llt
568 elbuf_str%GBUF%SMSTR(jj(1)+i) = x2s(i)
569 elbuf_str%GBUF%SMSTR(jj(2)+i) = y2s(i)
570 elbuf_str%GBUF%SMSTR(jj(3)+i) = x3s(i)
571 elbuf_str%GBUF%SMSTR(jj(4)+i) = y3s(i)
572 elbuf_str%GBUF%SMSTR(jj(5)+i) = x4s(i)
573 elbuf_str%GBUF%SMSTR(jj(6)+i) = y4s(i)
574 ENDDO
575 ENDIF
576C-----
577 IF (iuser == 1 .AND. mtn > 28) THEN
578C to be checked for IGTYP = 51
579 CALL cuserini(elbuf_str,
580 1 lft ,llt ,nft ,nel ,npt ,
581 2 istrain,sigsh ,numelc ,ixc ,nixc ,
582 3 nsigsh ,numshel,ptshel ,iun ,iun ,
583 4 nlay )
584 ENDIF
585
586 IF (iyldini == 1 .AND. (mtn== 36.OR. mtn==87))THEN
587 CALL cmatini(elbuf_str,
588 1 lft ,llt ,nft ,nel ,npt ,
589 2 istrain,sigsh ,numelc ,ixc ,nixc ,
590 3 nsigsh ,numshel,ptshel ,iun ,iun ,
591 4 nlay )
592 ENDIF
593
594c-----------------------------------------------------------------------
595c Failure model initialisation
596c-----------------------------------------------------------------------
597c set element damage initialized in /FAIL/FRACTAL_DMG
598 IF (fail_fractal%NFAIL > 0) THEN
599 CALL fractal_dmg_init(elbuf_str,mat_param,fail_fractal,
600 . nummat ,numelc ,nel ,nft ,ngl ,ity )
601 ENDIF
602c tag edge elements in local UVAR for /FAIL/ALTER (XFEM)
603 IF (ifail > 0 .and. iddlevel == 1) THEN
604 CALL fail_windshield_init(elbuf_str,mat_param,fail_brokmann,
605 . nel ,nft ,ity ,igrsh4n ,igrsh3n ,
606 . aldt ,thk ,ngl )
607 ENDIF
608C to be checked for IGTYP = 51
609 CALL cfailini(elbuf_str,mat_param,
610 . nptt ,nlay ,sigsh ,nsigsh ,ptshel ,
611 . rnoise ,perturb ,aldt ,thk )
612C-------------------------------------------
613C calculation of elementary dt
614C-------------------------------------------
615c IGTYP=GEO(12,IXC(6,I+NFT))
616 IF (igtyp /= 0 .AND. igtyp /= 1 .AND.
617 . igtyp /= 9 .AND. igtyp /= 10 .AND.
618 . igtyp /= 11 .AND. igtyp /= 16 .AND.
619 . igtyp /= 17 .AND. igtyp /= 51 .AND.
620 . igtyp /= 52) THEN
621 CALL ancmsg(msgid=25,
622 . anmode=aninfo,
623 . msgtype=msgerror,
624 . i1=id,
625 . c1=titr,
626 . i2=iprop)
627 ENDIF
628 ndepar=numels+nft
629 DO i=lft,llt
630 dtelem(ndepar+i)=dt(i)
631 ENDDO
632C---
633 IF (ixfem > 0) THEN
634 CALL cbufxfe(elbuf_str,xfem_str,isubstack,stack ,
635 . igeo ,geo ,lft ,llt ,mat,
636 . pid ,npt ,nptt ,nlay,ir ,
637 . is ,ixfem,mtn ,ng)
638 ENDIF
639C------------
640 ! Compute the initial volume
641 DO i=lft,llt
642 IF (gbuf%G_VOL > 0) gbuf%VOL(i) = area(i)*gbuf%THK(i)
643 ENDDO
644 IF (ixfem > 0) THEN
645 DO ixel=1,nxel
646 DO i=lft,llt
647 IF (xfem_str(ng,ixel)%GBUF%G_VOL > 0)
648 . xfem_str(ng,ixel)%GBUF%VOL(i) = area(i)*gbuf%THK(i)
649 END DO
650 ENDDO
651 ENDIF
652C---
653 IF (ALLOCATED(dir_b)) DEALLOCATE(dir_b)
654 IF (ALLOCATED(dir_a)) DEALLOCATE(dir_a)
655 DEALLOCATE(phi1,phi2,coor1,coor2,coor3,coor4)
656 DEALLOCATE(matly)
657 DEALLOCATE(posly)
658C-----------------------------
659 RETURN
c1buf3
subroutine c1buf3(geo, thk, off, thke, kshtree, shtree)
Definition c1buf3.F:32
cbufxfe
subroutine cbufxfe(elbuf_str, xfem_str, isubstack, stack, igeo, geo, lft, llt, mat, pid, npt, nptt, nlay, ir, is, ixfem, mtn, ng)
Definition cbufxfe.F:39
ccoori
subroutine ccoori(x, xrefc, ixc, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, ix1, ix2, ix3, ix4, ngl)
Definition ccoori.F:40
cderii
subroutine cderii(px1, px2, py1, py2, x1g, x2g, x3g, x4g, y1g, y2g, y3g, y4g, z1g, z2g, z3g, z4g, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, x2l, x3l, x4l, y2l, y3l, y4l)
Definition cderii.F:35
cdleni
subroutine cdleni(pm, geo, stifn, stifr, ixc, px1, px2, py1, py2, thk, igeo, dt, sh4tree, aldt, uparam, ipm, nlay, pm_stack, isubstack, strc, area, imat, iprop, x2l, x3l, x4l, y2l, y3l, y4l, igeo_stack, group_param)
Definition cdleni.F:36
cepsini
subroutine cepsini(elbuf_str, mat_param, jft, jlt, ismstr, ilaw, ithk, pm, geo, ixc, x, xrefc, for, thk, eint, gstr, nlay, px1g, px2g, py1g, py2g, x2s, y2s, x3s, y3s, x4s, y4s, uvar, ipm, igeo, imat, skew, nel, dir_a, dir_b, sigi, npf, tf, irep)
Definition cepsini.F:44
cepschk
subroutine cepschk(jft, jlt, nft, pm, geo, ixc, gstr, thk, nel, cpt_eltens)
Definition cepsini.F:829
ceveci
subroutine ceveci(jft, jlt, area, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z)
Definition ceveci.F:37
cfailini
subroutine cfailini(elbuf_str, mat_param, nptt, nlay, sigsh, nsigsh, ptsh, rnoise, perturb, aldt, thk)
Definition cfailini.F:36
csms11_ini
subroutine csms11_ini(jft, jlt, ixc, x, x2s, y2s, x3s, y3s, x4s, y4s)
Definition cinit3.F:671
cinmas
subroutine cinmas(x, xrefc, ix, geo, pm, ms, tiner, thke, ihbe, partsav, v, ipart, msc, inc, area, i8mi, igeo, etnod, imid, iprop, nshnod, stc, sh4tree, mcp, mcps, temp, ms_layer, zi_layer, ms_layerc, zi_layerc, msz2c, zply, isubstack, nlay, elbuf_str, stack, thki, rnoise, drape, nintemp, perturb, ix1, ix2, ix3, ix4, idrape, indx)
Definition cinmas.F:95
cm27in3
subroutine cm27in3(elbuf_str, geo, igeo, pm, ipm, ix, nix, nlay, ir, is, imat)
Definition cm27in3.F:38
cm35in3
subroutine cm35in3(elbuf_str, thk, area, nel, nlay, nptr, npts, nptt, igtyp)
Definition cm35in3.F:34
cm58in3
subroutine cm58in3(irep, dir1, dir2, uparam, uvar, aldt, nel, nuvar, tan_phi, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z)
Definition cm58in3.F:35
cmatini
subroutine cmatini(elbuf_str, jft, jlt, nft, nel, npt, istrain, sigsh, numel, ix, nix, nsigsh, numsh, ptsh, ir, is, nlay)
Definition cmatini.F:36
corth3
subroutine corth3(elbuf_str, dir_a, dir_b, jft, jlt, nlay, irep, nel, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, idrape, igtyp)
Definition corth3.F:42
corthdir
subroutine corthdir(elbuf_str, igeo, geo, vx, vy, vz, phi1, phi2, coor1, coor2, coor3, coor4, iorthloc, nlay, irep, isubstack, stack, geo_stack, igeo_stack, ir, is, nel, imat, iprop, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, npt_all, idrape)
Definition corthdir.F:42
corthini
subroutine corthini(jft, jlt, nft, nlay, numel, nsigsh, nix, ix, igeo, geo, skew, sigsh, ptsh, phi1, phi2, vx, vy, vz, coor1, coor2, coor3, coor4, iorthloc, isubstack, stack, irep, elbuf_str, drape, angle, x, geo_stack, e3x, e3y, e3z, betaorth, x1, x2, y1, y2, z1, z2, nel, g_add_node, add_node, npt_all, idrape, indx)
Definition corthini.F:48
my_real
#define my_real
Definition cppsort.cpp:32
csigini
subroutine csigini(elbuf_str, jft, jlt, nft, npt, istrain, thk, eint, gstr, hh, plas, for, mom, sigsh, nlay, g_hourg, numel, ix, nix, nsigsh, numsh, ptsh, igeo, thke, nel, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, isigsh, dir_a, dir_b, posly, igtyp)
Definition csigini.F:46
cuserini
subroutine cuserini(elbuf_str, jft, jlt, nft, nel, npt, istrain, sigsh, numel, ix, nix, nsigsh, numsh, ptsh, ir, is, nlay)
Definition cuserini.F:36
cveok3
subroutine cveok3(nvc, nod, ix1, ix2, ix3, ix4)
Definition cveok3.F:35
fail_brokmann
subroutine fail_brokmann(nel, nuparam, nuvar, time, timestep, uparam, ngl, signxx, signyy, signxy, uvar, off, ipt, nindxf, indxf, tdel)
Definition fail_brokmann.F:37
fail_windshield_init
subroutine fail_windshield_init(elbuf_str, mat_param, fail_brokmann, nel, nft, ity, igrsh4n, igrsh3n, aldt, thk, ngl)
Definition fail_windshield_init.F:40
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
law158_init
subroutine law158_init(dir1, dir2, uvar, aldt, nel, nuvar, tan_phi, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z)
Definition law158_init.F:35
layini1
subroutine layini1(elbuf_str, jft, jlt, geo, igeo, mat, pid, matly, posly, igtyp, nlay, npt, isubstack, stack, drape, nft, thk, nel, idrape, numel_drape, indx)
Definition layini1.F:40
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
drape_mod
Definition drape_mod.F:54
drape_mod::scdrape
integer scdrape
Definition drape_mod.F:92
groupdef_mod
Definition groupdef_mod.F:668
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
nlocal_reg_mod
Definition nlocal_reg_mod.F:119
stack_mod
Definition stack_mod.F:75
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
thickini
subroutine thickini(jft, jlt, nft, ptsh, numel, thk, thke, ix, nix, nsigsh, sigsh)
Definition thickini.F:36

◆ csms11_ini()

subroutine csms11_ini ( integer jft,
integer jlt,
integer, dimension(nixc,*) ixc,
x,
x2s,
y2s,
x3s,
y3s,
x4s,
y4s )

Definition at line 668 of file cinit3.F.

671C-----------------------------------------------
672C M o d u l e s
673C-----------------------------------------------
674 use element_mod , only : nixc
675C-----------------------------------------------
676C I m p l i c i t T y p e s
677C-----------------------------------------------
678#include "implicit_f.inc"
679C-----------------------------------------------
680C G l o b a l P a r a m e t e r s
681C-----------------------------------------------
682#include "mvsiz_p.inc"
683C-----------------------------------------------
684C D u m m y A r g u m e n t s
685C-----------------------------------------------
686 INTEGER JFT, JLT
687 INTEGER IXC(NIXC,*)
688 my_real
689 . x(3,*), x2s(*), y2s(*), x3s(*), y3s(*), x4s(*), y4s(*)
690C------------------------------------------------------
691C L o c a l V a r i a b l e s
692C-----------------------------------------------
693 INTEGER I, NIT, NT, I1, I2, I3, I4
694 my_real
695 . x1(mvsiz) , x2(mvsiz) , x3(mvsiz) , x4(mvsiz) ,
696 . y1(mvsiz) , y2(mvsiz) , y3(mvsiz) , y4(mvsiz) ,
697 . z1(mvsiz) , z2(mvsiz) , z3(mvsiz) , z4(mvsiz)
698 my_real
699 . suma,s1,s2
700
701 my_real
702 . e1x(mvsiz), e1y(mvsiz), e1z(mvsiz),
703 . e2x(mvsiz), e2y(mvsiz), e2z(mvsiz),
704 . e3x(mvsiz), e3y(mvsiz), e3z(mvsiz)
705C---------------------------------------------------------------
706 DO i=jft,jlt
707 x1(i)=zero
708 y1(i)=zero
709 z1(i)=zero
710 x2(i)=x(1,ixc(3,i))-x(1,ixc(2,i))
711 y2(i)=x(2,ixc(3,i))-x(2,ixc(2,i))
712 z2(i)=x(3,ixc(3,i))-x(3,ixc(2,i))
713 x3(i)=x(1,ixc(4,i))-x(1,ixc(2,i))
714 y3(i)=x(2,ixc(4,i))-x(2,ixc(2,i))
715 z3(i)=x(3,ixc(4,i))-x(3,ixc(2,i))
716 x4(i)=x(1,ixc(5,i))-x(1,ixc(2,i))
717 y4(i)=x(2,ixc(5,i))-x(2,ixc(2,i))
718 z4(i)=x(3,ixc(5,i))-x(3,ixc(2,i))
719 ENDDO
720 DO i=jft,jlt
721 e1x(i) = x2(i)+x3(i)-x1(i)-x4(i)
722 e1y(i) = y2(i)+y3(i)-y1(i)-y4(i)
723 e1z(i) = z2(i)+z3(i)-z1(i)-z4(i)
724C
725 e2x(i) = x3(i)+x4(i)-x1(i)-x2(i)
726 e2y(i) = y3(i)+y4(i)-y1(i)-y2(i)
727 e2z(i) = z3(i)+z4(i)-z1(i)-z2(i)
728C
729 e3x(i) = e1y(i)*e2z(i)-e1z(i)*e2y(i)
730 e3y(i) = e1z(i)*e2x(i)-e1x(i)*e2z(i)
731 e3z(i) = e1x(i)*e2y(i)-e1y(i)*e2x(i)
732 ENDDO
733C---
734 DO i=jft,jlt
735 suma = e3x(i)*e3x(i)+e3y(i)*e3y(i)+e3z(i)*e3z(i)
736 suma = one/max(sqrt(suma),em20)
737 e3x(i) = e3x(i)*suma
738 e3y(i) = e3y(i)*suma
739 e3z(i) = e3z(i)*suma
740C
741 s1 = e1x(i)*e1x(i)+e1y(i)*e1y(i)+e1z(i)*e1z(i)
742 s2 = e2x(i)*e2x(i)+e2y(i)*e2y(i)+e2z(i)*e2z(i)
743 suma = sqrt(s1/s2)
744 e1x(i) = e1x(i) + (e2y(i)*e3z(i)-e2z(i)*e3y(i))*suma
745 e1y(i) = e1y(i) + (e2z(i)*e3x(i)-e2x(i)*e3z(i))*suma
746 e1z(i) = e1z(i) + (e2x(i)*e3y(i)-e2y(i)*e3x(i))*suma
747C
748 suma = e1x(i)*e1x(i)+e1y(i)*e1y(i)+e1z(i)*e1z(i)
749 suma = one/max(sqrt(suma),em20)
750 e1x(i) = e1x(i)*suma
751 e1y(i) = e1y(i)*suma
752 e1z(i) = e1z(i)*suma
753C
754 e2x(i) = e3y(i) * e1z(i) - e3z(i) * e1y(i)
755 e2y(i) = e3z(i) * e1x(i) - e3x(i) * e1z(i)
756 e2z(i) = e3x(i) * e1y(i) - e3y(i) * e1x(i)
757 ENDDO
758 DO i=jft,jlt
759 x2s(i) = e1x(i)*x2(i) + e1y(i)*y2(i) + e1z(i)*z2(i)
760 y2s(i) = e2x(i)*x2(i) + e2y(i)*y2(i) + e2z(i)*z2(i)
761 x3s(i) = e1x(i)*x3(i) + e1y(i)*y3(i) + e1z(i)*z3(i)
762 y3s(i) = e2x(i)*x3(i) + e2y(i)*y3(i) + e2z(i)*z3(i)
763 x4s(i) = e1x(i)*x4(i) + e1y(i)*y4(i) + e1z(i)*z4(i)
764 y4s(i) = e2x(i)*x4(i) + e2y(i)*y4(i) + e2z(i)*z4(i)
765 ENDDO
766C-----------
767 RETURN