OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cbavit_ply.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!|| cbavit_ply ../engine/source/properties/composite_options/stack/cbavit_ply.f
25!||--- called by ------------------------------------------------------
26!|| cbaforc3 ../engine/source/elements/shell/coqueba/cbaforc3.F
27!||--- uses -----------------------------------------------------
28!|| plyxfem_mod ../engine/share/modules/plyxfem_mod.F
29!||====================================================================
30 SUBROUTINE cbavit_ply(JFT,JLT,IXC,OFFG,OFF,NPLAT,IPLAT,NPT,
31 1 VCORE,DI,ZL,VQ , VXYZ,X13_T ,X24_T ,
32 2 Y13_T,Y24_T,AREA,INOD,DEL_PLY,
33 . VNI,ISTACK,VR)
34C-----------------------------------------------
35C M o d u l e s
36C-----------------------------------------------
37 USE plyxfem_mod
38C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
39C transformation au repere local est absolue
40C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
41C-----------------------------------------------
42C I M P L I C I T T Y P E S
43C-----------------------------------------------
44#include "implicit_f.inc"
45c-----------------------------------------------
46c g l o b a l p a r a m e t e r s
47c-----------------------------------------------
48#include "mvsiz_p.inc"
49#include "scr17_c.inc"
50#include "com08_c.inc"
51#include "impl1_c.inc"
52C-----------------------------------------------
53C D U M M Y A R G U M E N T S
54C-----------------------------------------------
55 INTEGER IXC(NIXC,*),JFT,JLT,NNOD,NPLAT,IPLAT(*),NPT,
56 . ISTACK(MVSIZ,NPT),INOD(*)
57 my_real
58 . offg(*),off(*)
59 parameter(nnod = 4)
60 my_real
61 . vcore(mvsiz,3,nnod),vxyz(mvsiz,3*nnod,npt),
62 . vq(mvsiz,3,3),zl(*),di(mvsiz,6),
63 . y24_t(*),x13_t(*),x24_t(*),y13_t(*),area(*),
64 . vni(4,4), del_ply(mvsiz,12,npt), vr(3,*)
65C-----------------------------------------------
66C L O C A L V A R I A B L E S
67C-----------------------------------------------
68 INTEGER J,I,K,EP,NN(4), IPLY,IP,NPLAT0
69 INTEGER L,M,J1,J2
70 my_real
71 . pg,z1,z2,mx13,my13,mx23,my23,mx34,my34,gama1,gama2, x21,
72 . x34,y21,y34 ,z21,z34,l12,l34,x41,x32,y41,y32,z41,z32,xx1
73 . yy,xy,xz1,yz ,zz,y24,x24,y13,x13,corel(3,4),xx1,yy,off_l,
74 . d1,d2,dt05,dt025,exz,eyz,ddry,v13x,v24x,vhix,ddrz1,ddrz2,
75 . ddrz,ddrx,a1,a2,d3
76 my_real
77 . vg13(3),vg24(3),vghi(3),v13(mvsiz,3), v24(mvsiz,3),
78 . vhi(mvsiz,3), ar(3),d(6),alr(3),rr(3,nnod),
79 . area_i(mvsiz), del_iply(mvsiz,3,npt-1),dn_iply(mvsiz,12,npt-1),
80 . dn_ply(mvsiz,12,npt)
81 DATA pg/.577350269189626/
82C-----------------------------------------------
83!! NPLAT0 = JLT
84 nplat0 = nplat
85 a1 = one - pg
86 a2 = one + pg
87C (-PG, -PG)
88 vni(1,1)= fourth*a2*a2
89 vni(2,1)= fourth*a1*a2
90 vni(3,1)= fourth*a1*a1
91 vni(4,1)= vni(2,1)
92C (xi = PG,eta=-PG)
93 vni(1,2)= vni(2,1)
94 vni(2,2)= vni(1,1)
95 vni(3,2)= vni(2,1)
96 vni(4,2)= vni(3,1)
97C xi = PG, eta = PG
98 vni(1,3)= vni(3,1)
99 vni(2,3)= vni(2,1)
100 vni(3,3)= vni(1,1)
101 vni(4,3)= vni(2,1)
102Cxi= -PG, eta= pG
103 vni(1,4)= vni(2,1)
104 vni(2,4)= vni(3,1)
105 vni(3,4)= vni(2,1)
106 vni(4,4)= vni(1,1)
107C
108 DO ep=jft,jlt
109 area_i(ep)=one/max(em20,area(ep))
110 ENDDO
111C
112 DO j=1, npt
113 DO ep=jft,jlt
114 ip = istack(ep,j)
115
116 nn(1) = inod(ixc(2,ep))
117 nn(2) = inod(ixc(3,ep))
118 nn(3) = inod(ixc(4,ep))
119 nn(4) = inod(ixc(5,ep))
120C
121 vg13(1)=ply(ip)%V(1,nn(1)) - ply(ip)%V(1,nn(3))
122 vg24(1)=ply(ip)%V(1,nn(2)) - ply(ip)%V(1,nn(4))
123 vghi(1)=ply(ip)%V(1,nn(1)) - ply(ip)%V(1,nn(2))
124 . + ply(ip)%V(1,nn(3)) - ply(ip)%V(1,nn(4))
125C
126 vg13(2)=ply(ip)%V(2,nn(1)) - ply(ip)%V(2,nn(3))
127 vg24(2)=ply(ip)%V(2,nn(2)) - ply(ip)%V(2,nn(4))
128 vghi(2)=ply(ip)%V(2,nn(1)) - ply(ip)%V(2,nn(2))
129 . +ply(ip)%V(2,nn(3)) - ply(ip)%V(2,nn(4))
130C
131 vg13(3)= ply(ip)%V(3,nn(1)) - ply(ip)%V(3,nn(3))
132 vg24(3)= ply(ip)%V(3,nn(2)) - ply(ip)%V(3,nn(4))
133 vghi(3)= ply(ip)%V(3,nn(1)) - ply(ip)%V(3,nn(2))
134 . + ply(ip)%V(3,nn(3)) - ply(ip)%V(3,nn(4))
135C
136 v13(ep,1) =(vq(ep,1,1)*vg13(1)+vq(ep,2,1)*vg13(2)
137 1 +vq(ep,3,1)*vg13(3))
138 v24(ep,1)=(vq(ep,1,1)*vg24(1)+vq(ep,2,1)*vg24(2)
139 1 +vq(ep,3,1)*vg24(3))
140 vhi(ep,1)=(vq(ep,1,1)*vghi(1)+vq(ep,2,1)*vghi(2)
141 1 +vq(ep,3,1)*vghi(3))
142 v13(ep,2)=(vq(ep,1,2)*vg13(1)+vq(ep,2,2)*vg13(2)
143 1 +vq(ep,3,2)*vg13(3))
144 v24(ep,2)=(vq(ep,1,2)*vg24(1)+vq(ep,2,2)*vg24(2)
145 1 +vq(ep,3,2)*vg24(3))
146 vhi(ep,2)=(vq(ep,1,2)*vghi(1)+vq(ep,2,2)*vghi(2)
147 1 +vq(ep,3,2)*vghi(3))
148 v13(ep,3)=(vq(ep,1,3)*vg13(1)+vq(ep,2,3)*vg13(2)
149 1 +vq(ep,3,3)*vg13(3))
150 v24(ep,3)=(vq(ep,1,3)*vg24(1)+vq(ep,2,3)*vg24(2)
151 1 +vq(ep,3,3)*vg24(3))
152 vhi(ep,3)=(vq(ep,1,3)*vghi(1)+vq(ep,2,3)*vghi(2)
153 1 +vq(ep,3,3)*vghi(3))
154
155 ENDDO
156C ---
157 IF (impl_s==0) THEN
158 dt05 = half*dt1
159 dt025 =fourth*dt1
160 DO i=jft,jlt
161 exz = y24_t(i)*v13(i,3)-y13_t(i)*v24(i,3)
162 eyz = -x24_t(i)*v13(i,3)+x13_t(i)*v24(i,3)
163 ddry=dt05*exz*area_i(i)
164 ddrx=dt05*eyz*area_i(i)
165 v13x = v13(i,1)
166 v24x = v24(i,1)
167 vhix = vhi(i,1)
168 ddrz1=dt025*(v13(i,2)-v24(i,2))/(x13_t(i)-x24_t(i))
169 IF (abs(x13_t(i)-x24_t(i))<em10) ddrz1 = zero
170 v13(i,1) = v13(i,1)-ddry*v13(i,3)-ddrz1*v13(i,2)
171 v24(i,1) = v24(i,1)-ddry*v24(i,3)-ddrz1*v24(i,2)
172 vhi(i,1) = vhi(i,1)-ddry*vhi(i,3)-ddrz1*vhi(i,2)
173 ddrz2=dt025*(v13x+v24x)/(y13_t(i)+y24_t(i))
174 IF (abs(y13_t(i)+y24_t(i))<em10) ddrz2 = zero
175 v13(i,2) = v13(i,2)-ddrx*v13(i,3)-ddrz2*v13x
176 v24(i,2) = v24(i,2)-ddrx*v24(i,3)-ddrz2*v24x
177 vhi(i,2) = vhi(i,2)-ddrx*vhi(i,3)-ddrz2*vhix
178 ENDDO
179 ENDIF
180C ---
181#include "vectorize.inc"
182 DO i=jft,nplat0
183 ep =iplat(i)
184 vxyz(ep,1,j)=v13(ep,1)
185 vxyz(ep,2,j)=v13(ep,2)
186 vxyz(ep,3,j)=v13(ep,3)
187C
188 vxyz(ep,4,j)=v24(ep,1)
189 vxyz(ep,5,j)=v24(ep,2)
190 vxyz(ep,6,j)=v24(ep,3)
191C
192 vxyz(ep,7,j)=vhi(ep,1)
193 vxyz(ep,8,j)=vhi(ep,2)
194 vxyz(ep,9,j)=vhi(ep,3)
195C
196 ENDDO
197#include "vectorize.inc"
198 DO i=nplat0+1,jlt
199 ep =iplat(i)
200 z1=zl(ep)
201 z2=z1*z1
202C
203
204 x13 =(vcore(ep,1,1)-vcore(ep,1,3))*half
205 x24 =(vcore(ep,1,2)-vcore(ep,1,4))*half
206 y13 =(vcore(ep,2,1)-vcore(ep,2,3))*half
207 y24 =(vcore(ep,2,2)-vcore(ep,2,4))*half
208 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
209 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
210
211C--------------------------
212 ar(1)=-z1*vhi(ep,2)+y13*v13(ep,3)+y24*v24(ep,3)+my13*vhi(ep,3)
213 ar(2)= z1*vhi(ep,1)-x13*v13(ep,3)-x24*v24(ep,3)-mx13*vhi(ep,3)
214 ar(3)= x13*v13(ep,2)+x24*v24(ep,2)+mx13*vhi(ep,2)
215 1 -y13*v13(ep,1)-y24*v24(ep,1)-my13*vhi(ep,1)
216C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
217 alr(1) =di(ep,1)*ar(1)+di(ep,4)*ar(2)+di(ep,5)*ar(3)
218 alr(2) =di(ep,4)*ar(1)+di(ep,2)*ar(2)+di(ep,6)*ar(3)
219 alr(3) =di(ep,5)*ar(1)+di(ep,6)*ar(2)+di(ep,3)*ar(3)
220C
221 v13(ep,1)= half*v13(ep,1)+alr(3)*y13
222 v24(ep,1)= half*v24(ep,1)+alr(3)*y24
223 vhi(ep,1)= fourth*vhi(ep,1)+(alr(3)*my13-z1*alr(2))
224 v13(ep,2)= half*v13(ep,2)-alr(3)*x13
225 v24(ep,2)= half*v24(ep,2)-alr(3)*x24
226 vhi(ep,2)= fourth*vhi(ep,2)-(alr(3)*mx13-z1*alr(1))
227 v13(ep,3)= half*v13(ep,3)-(y13*alr(1)-x13*alr(2))
228 v24(ep,3)= half*v24(ep,3)-(y24*alr(1)-x24*alr(2))
229 vhi(ep,3)= fourth*vhi(ep,3)+(mx13*alr(2)-my13*alr(1))
230C
231 vxyz(ep,1 ,j) = v13(ep,1) + vhi(ep,1)
232 vxyz(ep,4 ,j) = v24(ep,1) - vhi(ep,1)
233 vxyz(ep,7 ,j) = -v13(ep,1) + vhi(ep,1)
234 vxyz(ep,10,j) = -v24(ep,1) - vhi(ep,1)
235C
236 vxyz(ep,2 ,j) = v13(ep,2) + vhi(ep,2)
237 vxyz(ep,5 ,j) = v24(ep,2) - vhi(ep,2)
238 vxyz(ep,8 ,j) = -v13(ep,2) + vhi(ep,2)
239 vxyz(ep,11,j) = -v24(ep,2) - vhi(ep,2)
240C
241 vxyz(ep,3 ,j) = v13(ep,3) + vhi(ep,3)
242 vxyz(ep,6 ,j) = v24(ep,3) - vhi(ep,3)
243 vxyz(ep,9 ,j) = -v13(ep,3) + vhi(ep,3)
244 vxyz(ep,12,j) = -v24(ep,3) - vhi(ep,3)
245 ENDDO
246 ENDDO
247C----------------------------
248C---------Facteur pour longueur caracteristique---
249C----------------------------
250 off_l = 0.
251 DO ep=jft,jlt
252 off(ep) = min(one,abs(offg(ep)))
253 off_l = min(off_l,offg(ep))
254 ENDDO
255 IF(off_l < zero)THEN
256 DO j=1,npt
257 DO ep=jft,jlt
258 IF(offg(ep) < zero)THEN
259 vxyz(ep,1 ,j)= zero
260 vxyz(ep,2 ,j)= zero
261 vxyz(ep,3 ,j)= zero
262 vxyz(ep,4 ,j)= zero
263 vxyz(ep,5 ,j)= zero
264 vxyz(ep,6 ,j)= zero
265 vxyz(ep,7 ,j)= zero
266 vxyz(ep,8 ,j)= zero
267 vxyz(ep,9 ,j)= zero
268 vxyz(ep,10,j)= zero
269 vxyz(ep,11,j)= zero
270 vxyz(ep,12,j)= zero
271 ENDIF
272 ENDDO
273 ENDDO
274 ENDIF
275C ---
276C les deplacements nodaux dans le rep re local
277 DO j=1, npt
278 DO ep=jft,jlt
279 ip = istack(ep,j)
280 nn(1) = inod(ixc(2,ep))
281 nn(2) = inod(ixc(3,ep))
282 nn(3) = inod(ixc(4,ep))
283 nn(4) = inod(ixc(5,ep))
284C node 1 layer J
285 d1 = ply(ip)%U(1,nn(1))
286 d2 = ply(ip)%U(2,nn(1))
287 d3 = ply(ip)%U(3,nn(1))
288 dn_ply(ep,1,j)=vq(ep,1,1)*d1 + vq(ep,2,1)*d2 + vq(ep,3,1)*d3
289 dn_ply(ep,2,j)=vq(ep,1,2)*d1 + vq(ep,2,2)*d2 + vq(ep,3,2)*d3
290 dn_ply(ep,3,j)=vq(ep,1,3)*d1 + vq(ep,2,3)*d2 + vq(ep,3,3)*d3
291C node 2 layer J
292 d1 = ply(ip)%U(1,nn(2))
293 d2 = ply(ip)%U(2,nn(2))
294 d3 = ply(ip)%U(3,nn(2))
295 dn_ply(ep,4,j)=vq(ep,1,1)*d1 + vq(ep,2,1)*d2 + vq(ep,3,1)*d3
296 dn_ply(ep,5,j)=vq(ep,1,2)*d1 + vq(ep,2,2)*d2 + vq(ep,3,2)*d3
297 dn_ply(ep,6,j)=vq(ep,1,3)*d1 + vq(ep,2,3)*d2 + vq(ep,3,3)*d3
298C node 3 layer J
299 d1 = ply(ip)%U(1,nn(3))
300 d2 = ply(ip)%U(2,nn(3))
301 d3 = ply(ip)%U(3,nn(3))
302 dn_ply(ep,7,j)=vq(ep,1,1)*d1 + vq(ep,2,1)*d2 + vq(ep,3,1)*d3
303 dn_ply(ep,8,j)=vq(ep,1,2)*d1 + vq(ep,2,2)*d2 + vq(ep,3,2)*d3
304 dn_ply(ep,9,j)=vq(ep,1,3)*d1 + vq(ep,2,3)*d2 + vq(ep,3,3)*d3
305C node 4 layer J
306 d1 = ply(ip)%U(1,nn(4))
307 d2 = ply(ip)%U(2,nn(4))
308 d3 = ply(ip)%U(3,nn(4))
309 dn_ply(ep,10,j)=vq(ep,1,1)*d1 + vq(ep,2,1)*d2 + vq(ep,3,1)*d3
310 dn_ply(ep,11,j)=vq(ep,1,2)*d1 + vq(ep,2,2)*d2 + vq(ep,3,2)*d3
311 dn_ply(ep,12,j)=vq(ep,1,3)*d1 + vq(ep,2,3)*d2 + vq(ep,3,3)*d3
312 ENDDO
313 ENDDO
314C glissement elementaire
315 DO j=1, npt
316 DO ep=jft,jlt
317 del_ply(ep,1,j) = (
318 . dn_ply(ep,1,j)*vni(1,1) + dn_ply(ep,4,j)*vni(2,1) +
319 . dn_ply(ep,7,j)*vni(3,1) + dn_ply(ep,10,j)*vni(4,1) )
320 del_ply(ep,2,j) = (
321 . dn_ply(ep,2,j)*vni(1,1) + dn_ply(ep,5,j)*vni(2,1) +
322 . dn_ply(ep,8,j)*vni(3,1) + dn_ply(ep,11,j)*vni(4,1) )
323 del_ply(ep,3,j) = (
324 . dn_ply(ep,3,j)*vni(1,1) + dn_ply(ep,6,j)*vni(2,1) +
325 . dn_ply(ep,9,j)*vni(3,1) + dn_ply(ep,12,j)*vni(4,1) )
326C 2nd gauss p
327 del_ply(ep,4,j) = (
328 . dn_ply(ep,1,j)*vni(1,2) + dn_ply(ep,4,j)*vni(2,2) +
329 . dn_ply(ep,7,j)*vni(3,2) + dn_ply(ep,10,j)*vni(4,2) )
330 del_ply(ep,5,j) = (
331 . dn_ply(ep,2,j)*vni(1,2) + dn_ply(ep,5,j)*vni(2,2) +
332 . dn_ply(ep,8,j)*vni(3,2) + dn_ply(ep,11,j)*vni(4,2) )
333 del_ply(ep,6,j) = (
334 . dn_ply(ep,3,j)*vni(1,2) + dn_ply(ep,6,j)*vni(2,2) +
335 . dn_ply(ep,9,j)*vni(3,2) + dn_ply(ep,12,j)*vni(4,2) )
336c 3rd
337 del_ply(ep,7,j) = (
338 . dn_ply(ep,1,j)*vni(1,3) + dn_ply(ep,4,j)*vni(2,3) +
339 . dn_ply(ep,7,j)*vni(3,3) + dn_ply(ep,10,j)*vni(4,3) )
340 del_ply(ep,8,j) = (
341 . dn_ply(ep,2,j)*vni(1,3) + dn_ply(ep,5,j)*vni(2,3) +
342 . dn_ply(ep,8,j)*vni(3,3) + dn_ply(ep,11,j)*vni(4,3) )
343 del_ply(ep,9,j) = (
344 . dn_ply(ep,3,j)*vni(1,3) + dn_ply(ep,6,j)*vni(2,3) +
345 . dn_ply(ep,9,j)*vni(3,3) + dn_ply(ep,12,j)*vni(4,3) )
346C 4th
347 del_ply(ep,10,j) = (
348 . dn_ply(ep,1,j)*vni(1,4) + dn_ply(ep,4,j)*vni(2,4) +
349 . dn_ply(ep,7,j)*vni(3,4) + dn_ply(ep,10,j)*vni(4,4) )
350 del_ply(ep,11,j) = (
351 . dn_ply(ep,2,j)*vni(1,4) + dn_ply(ep,5,j)*vni(2,4) +
352 . dn_ply(ep,8,j)*vni(3,4) + dn_ply(ep,11,j)*vni(4,4) )
353 del_ply(ep,12,j) = (
354 . dn_ply(ep,3,j)*vni(1,4) + dn_ply(ep,6,j)*vni(2,4) +
355 . dn_ply(ep,9,j)*vni(3,4) + dn_ply(ep,12,j)*vni(4,4) )
356 ENDDO
357 ENDDO
358C ---
359 RETURN
360 END
cbavit_ply
subroutine cbavit_ply(jft, jlt, ixc, offg, off, nplat, iplat, npt, vcore, di, zl, vq, vxyz, x13_t, x24_t, y13_t, y24_t, area, inod, del_ply, vni, istack, vr)
Definition cbavit_ply.F:34
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
plyxfem_mod
Definition plyxfem_mod.F:31
plyxfem_mod::ply
type(ply_data), dimension(:), allocatable ply
Definition plyxfem_mod.F:91