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