1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine resultsmech_u1(co,kon,ipkon,lakon,ne,v, 20 & stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero, 21 & ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr, 22 & iprestr,eme,iperturb,fn,iout,qa,vold,nmethod, 23 & veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon, 24 & xstateini,xstiff,xstate,npmat_,matname,mi,ielas,icmd, 25 & ncmat_,nstate_,stiini,vini,ener,eei,enerini,istep,iinc, 26 & reltime,calcul_fn,calcul_qa,calcul_cauchy,nener, 27 & ikin,nal,ne0,thicke,emeini,i,ielprop,prop,t0g,t1g) 28! 29! calculates nal,qa,fn,xstiff,ener,eme,eei,stx for user element 1 30! 31! This is a beam type element. Reference: 32! Yunhua Luo, An Efficient 3D Timoshenko Beam Element with 33! Consistent Shape Functions, Adv. Theor. Appl. Mech., Vol. 1, 34! 2008, no. 3, 95-106 35! 36! 37! special case for which the beam axis goes through the 38! center of gravity of the cross section and the 1-direction 39! corresponds with a principal axis 40! 41! note that the strain components are Lagrange strain components, 42! the stress components are Piola-Kirchhoff components of the 43! second kind. For linear geometric calculations the strains 44! reduce to the infinitesimal strains 45! 46! 47! INPUT: 48! 49! co(1..3,i) coordinates of node i 50! kon(*) contains the topology of all elements. The 51! topology of element i starts at kon(ipkon(i)+1) 52! and continues until all nodes are covered. The 53! number of nodes depends on the element label 54! ipkon(i) points to the location in field kon preceding 55! the topology of element i 56! lakon(i) label of element i (character*8) 57! ne highest element number in the mesh 58! v(0..mi(2),i) value of the variables in node i at the end 59! of the present iteration 60! elcon elastic constants (cf. List of variables and 61! their meaning in the User's Manual) 62! nelcon integers describing the elastic constant fields 63! (cf. User's Manual) 64! rhcon density constants (cf. User's Manual) 65! nrhcon integers describing the density constant fields 66! (cf. User's Manual) 67! alcon thermal expansion constants (cf. User's Manual) 68! nalcon integers describing the thermal expansion 69! constants (cf. User's Manual) 70! alzero thermal expansion reference values (cf. User's Manual) 71! ielmat(i) material for element i 72! ielorien(i) orientation for element i 73! norien number of orientations 74! orab(7,*) description of all local coordinate systems. 75! (cf. List of variables and their meaning in the 76! User's manual) 77! ntmat_ maximum number of material temperature points 78! t0(i) temperature in node i at start of calculation 79! t1(i) temperature in node i at the end of the current 80! increment 81! ithermal(1..2) cf. List of variables and 82! their meaning in the User's Manual 83! prestr(i,j,k) residual stress component i in integration point j 84! of element k 85! iprestr if 0: no residual stresses 86! else: residual stresses 87! iperturb(*) describes the kind of nonlinearity of the 88! calculation, cf. User's Manual 89! iout if -2: v is assumed to be known and is used to 90! calculate strains, stresses..., no result output 91! corresponds to iout=-1 with in addition the 92! calculation of the internal energy density 93! if -1: v is assumed to be known and is used to 94! calculate strains, stresses..., no result 95! output; is used to take changes in SPC's and 96! MPC's at the start of a new increment or 97! iteration into account 98! if 0: v is calculated from the system solution 99! and strains, stresses.. are calculated, 100! no result output 101! if 1: v is calculated from the system solution and 102! strains, stresses.. are calculated, requested 103! results output 104! if 2: v is assumed to be known and is used to 105! calculate strains, stresses..., requested 106! results output 107! vold(0..mi(2),i) value of the variables in node i at the end 108! of the previous iteration 109! nmethod procedure: 110! 1: static analysis 111! 2: frequency analysis 112! 3: buckling analysis 113! 4: (linear or nonlinear) dynamic analysis 114! 5: steady state dynamics analysis 115! 6: Coriolis frequency calculation 116! 7: flutter frequency calculation 117! 8: magnetostatics 118! 9: magnetodynamics 119! 10: electromagnetic eigenvalue problems 120! 11: superelement creation or Green function 121! calculation 122! 12: sensitivity analysis 123! veold(j,i) time rate of variable j in node i at the end 124! of the previous iteration 125! dtime length of present time increment 126! time step time at the end of the present increment 127! ttime total time at the start of the present increment 128! plicon,nplicon fields describing isotropic hardening of 129! a plastic material or spring constants of 130! a nonlinear spring (cf. User's Manual) 131! plkcon,nplkcon fields describing kinematic hardening of 132! a plastic material or gap conductance 133! constants (cf. User's Manual) 134! xstateini(i,j,k) state variable component i in integration point j 135! of element k at the start of the present increment 136! xstate(i,j,k) state variable component i in integration point j 137! of element k at the end of the present increment 138! npmat_ maximum number of plastic constants 139! matname(i) name of material i 140! mi(1) max # of integration points per element (max 141! over all elements) 142! mi(2) max degree of freedom per node (max over all 143! nodes) in fields like v(0:mi(2))... 144! mi(3) max number of layers in the structure 145! ielas 0: no elastic iteration: irreversible effects 146! are allowed 147! 1: elastic iteration, i.e. no irreversible 148! deformation allowed 149! icmd not equal to 3: calculate stress and stiffness 150! 3: calculate only stress 151! ncmat_ max number of elastic constants 152! nstate_ max number of state variables in any integration 153! point 154! stiini(i,j,k) stress component i in integration point j 155! of element k at the start of the present 156! increment (= end of last increment) 157! vini(0..mi(2),i) value of the variables in node i at the start 158! of the present increment 159! enerini(j,k) internal energy density at integration point j 160! of element k at the start of the present increment 161! istep current step number 162! iinc current increment number within the actual step 163! reltime relative step time (between 0 and 1) 164! calcul_fn if 0: no nodal forces have to be calculated 165! else: nodal forces are required on output 166! calcul_qa if 0: no mean forces have to be calculated 167! else: mean forces are required on output 168! calcul_cauchy if 0: no Cauchy stresses are required 169! else: Cauchy stresses are required on output: have 170! to be calculated from the PK2 stresses 171! nener if 0: internal energy calculation is not required 172! else: internal energy is required on output 173! ikin if 0: kinetic energy calculation is not requred 174! else: kinetic energy is required on output 175! ne0 largest element number without contact elements (are 176! stored after all other elements) 177! thicke(j,i) layer thickness for layer j in element i 178! emeini(i,j,k) mechanical strain component i in integration point j 179! of element k at the start of the present increment 180! i actual element at stake 181! ielprop(i) points to the location in field prop preceding 182! the properties of element i 183! prop(*) contains the properties and some beam 184! elements (cf. User's Manual) 185! t0g(1..2,i) temperature gradient in node i at start of calculation 186! t1g(1..2,i) temperature gradient in node i at the end of the 187! current increment 188! 189! 190! OUTPUT: 191! 192! stx(i,j,k) stress component i in integration point j 193! of element k at the end of the present 194! iteration 195! eme(i,j,k) mechanical strain component i in integration point j 196! of element k at the end of the present iteration 197! fn(j,i) nodal force component j in node i 198! qa(1..4) qa(1): average force 199! qa(2): average flux 200! ... cf. User's Manual 201! xstiff(i,j,k) stiffness (i=1...21) and conductivity 202! (i=22..27) constants in integration point j 203! of element k 204! ener(j,k) internal energy density at integration point j 205! of element k at the end of the present increment 206! eei(i,j,k) total strain component i in integration point j 207! of element k at the end of the present iteration 208! nal number of nodal force contributions 209! 210 implicit none 211! 212 character*8 lakon(*) 213 character*80 amat,matname(*) 214! 215 integer kon(*),konl(26),mi(*),kal(2,6),j1,j2,j3,j4, 216 & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*), 217 & ielorien(mi(3),*),ntmat_,ipkon(*),ne0, 218 & istep,iinc,ne,mattyp,ithermal(*),iprestr,i,j,k,m1,m2,jj, 219 & i1,kk,nener,indexe,nope,norien,iperturb(*),iout, 220 & nal,icmd,ihyper,nmethod,kode,imat,iorien,ielas, 221 & istiff,ncmat_,nstate_,ikin,ielprop(*), 222 & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_fn, 223 & calcul_cauchy,calcul_qa,index,node 224! 225 real*8 co(3,*),v(0:mi(2),*),stiini(6,mi(1),*),xa(3,3), 226 & stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*), 227 & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*), 228 & alcon(0:6,ntmat_,*),vini(0:mi(2),*), 229 & alzero(*),orab(7,*),elas(21),rho,fn(0:mi(2),*), 230 & q(0:mi(2),26),t0(*),t1(*),prestr(6,mi(1),*),eme(6,mi(1),*), 231 & vold(0:mi(2),*),eloc(9),elconloc(21),eth(6),coords(3), 232 & ener(mi(1),*),emec(6),eei(6,mi(1),*),enerini(mi(1),*), 233 & veold(0:mi(2),*),e,un,um,tt,dl,qa(*),t0l,t1l,dtime,time,ttime, 234 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 235 & xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802), 236 & xstateini(nstate_,mi(1),*),tm(3,3),a,reltime, 237 & thicke(mi(3),*),emeini(6,mi(1),*),aly,alz,bey,bez,xi(2), 238 & vlp(6,2),xi11,xi12,xi22,xk,offset1,offset2,e1(3),e2(3),e3(3), 239 & t0g(2,*),t1g(2,*) 240! 241 kal=reshape((/1,1,2,2,3,3,1,2,1,3,2,3/),(/2,6/)) 242! 243 nope=2 244! 245 indexe=ipkon(i) 246! 247! material and orientation 248! 249 imat=ielmat(1,i) 250 amat=matname(imat) 251 if(norien.gt.0) then 252 iorien=max(0,ielorien(1,i)) 253 else 254 iorien=0 255 endif 256 if(iorien.gt.0) then 257 write(*,*) '*ERROR in resultsmech_u1: no orientation' 258 write(*,*) ' calculation for this type of element' 259 call exit(201) 260 endif 261! 262 if(nelcon(1,imat).lt.0) then 263 ihyper=1 264 else 265 ihyper=0 266 endif 267! 268! properties of the cross section 269! 270 index=ielprop(i) 271 a=prop(index+1) 272 xi11=prop(index+2) 273 xi12=prop(index+3) 274 xi22=prop(index+4) 275 xk=prop(index+5) 276 e2(1)=prop(index+6) 277 e2(2)=prop(index+7) 278 e2(3)=prop(index+8) 279 offset1=prop(index+9) 280 offset2=prop(index+10) 281! 282 if(dabs(xi12).gt.0.d0) then 283 write(*,*) '*ERROR in resultsmech_u1: no nonzero cross moment' 284 write(*,*) ' of inertia for this type of element' 285 call exit(201) 286 endif 287! 288 if(dabs(offset1).gt.0.d0) then 289 write(*,*) '*ERROR in resultsmech_u1: no offset in 1-direction' 290 write(*,*) ' for this type of element' 291 call exit(201) 292 endif 293! 294 if(dabs(offset2).gt.0.d0) then 295 write(*,*) '*ERROR in resultsmech_u1: no offset in 2-direction' 296 write(*,*) ' for this type of element' 297 call exit(201) 298 endif 299! 300! computation of the coordinates and displacements of the local nodes 301! 302 do k=1,nope 303 konl(k)=kon(indexe+k) 304 do j=1,3 305 xl(j,k)=co(j,konl(k)) 306 vl(j,k)=v(j,konl(k)) 307 enddo 308 enddo 309! 310! q contains the nodal forces per element; initialization of q 311! 312 if((iperturb(1).ge.2).or.((iperturb(1).le.0).and.(iout.lt.1))) 313 & then 314 do m1=1,nope 315 do m2=0,mi(2) 316 q(m2,m1)=fn(m2,konl(m1)) 317 enddo 318 enddo 319 endif 320! 321! local axes e1-e2-e3 (e2 is given by the user (in the input deck 322! this is called e1), e1 is parallel to the beam axis) 323! 324 do j=1,3 325 e1(j)=xl(j,2)-xl(j,1) 326 enddo 327! 328 dl=dsqrt(e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3)) 329! 330 do j=1,3 331 e1(j)=e1(j)/dl 332 enddo 333! 334! e3 = e1 x e2 335! 336 e3(1)=e1(2)*e2(3)-e1(3)*e2(2) 337 e3(2)=e1(3)*e2(1)-e1(1)*e2(3) 338 e3(3)=e1(1)*e2(2)-e1(2)*e2(1) 339! 340! transformation matrix from the global into the local system 341! 342 do j=1,3 343 tm(1,j)=e1(j) 344 tm(2,j)=e2(j) 345 tm(3,j)=e3(j) 346 enddo 347! 348! calculating the temperature in the integration 349! point 350! 351 t0l=0.d0 352 t1l=0.d0 353 if(ithermal(1).eq.1) then 354 do i1=1,nope 355 t0l=t0l+t0(konl(i1))/2.d0 356 t1l=t1l+t1(konl(i1))/2.d0 357 enddo 358 elseif(ithermal(1).ge.2) then 359 write(*,*) '*ERROR in resultsmech_u1: no thermal' 360 write(*,*) ' calculation for this type of element' 361 call exit(201) 362 endif 363 tt=t1l-t0l 364! 365 kode=nelcon(1,imat) 366! 367 if(kode.eq.2) then 368 mattyp=1 369 else 370 mattyp=0 371 endif 372! 373! material data and local stiffness matrix 374! 375 istiff=0 376 call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 377 & imat,amat,iorien,coords,orab,ntmat_,elas,rho, 378 & i,ithermal,alzero,mattyp,t0l,t1l, 379 & ihyper,istiff,elconloc,eth,kode,plicon, 380 & nplicon,plkcon,nplkcon,npmat_, 381 & plconloc,mi(1),dtime,kk, 382 & xstiff,ncmat_) 383! 384! correcting the thermal strains 385! 386 do k=2,6 387 eth(k)=0.d0 388 enddo 389! 390 if(mattyp.eq.1) then 391 e=elconloc(1) 392 un=elconloc(2) 393 um=e/(1.d0+un) 394 um=um/2.d0 395 else 396 write(*,*) '*ERROR in resultsmech_u1: no anisotropic material' 397 write(*,*) ' calculation for this type of element' 398 call exit(201) 399 endif 400! 401! calculating the local displacement and rotation values 402! 403! values 1..6 of vl are u,v,w,phi,psi,theta from Luo 404! 405 do k=1,nope 406 node=kon(indexe+k) 407 do j=1,3 408 vl(j,k)=tm(j,1)*v(1,node) 409 & +tm(j,2)*v(2,node) 410 & +tm(j,3)*v(3,node) 411 vl(j+3,k)=tm(j,1)*v(4,node) 412 & +tm(j,2)*v(5,node) 413 & +tm(j,3)*v(6,node) 414 enddo 415 enddo 416! 417 xi(1)=0.d0 418 xi(2)=1.d0 419! 420 aly=12.d0*e*xi11/(xk*um*a*dl*dl) 421 alz=12.d0*e*xi22/(xk*um*a*dl*dl) 422 bey=1.d0/(1.d0-aly) 423 bez=1.d0/(1.d0-alz) 424! 425 do jj=1,nope 426! 427! calculating the derivative of the local displacement and 428! rotation values (from Eqn. (19) and (20) in Luo) 429! 430! Equation (19) does not seem to be correct in Luo: 431! - in the third equation v_1 shoud be w_1 and the sign in front 432! of the psi-terms should be negative 433! - in the last equation the sign before the w-terms should be 434! negative 435! 436 vlp(1,jj)=(-vl(1,1)+vl(1,2))/dl 437! 438 vlp(2,jj)= 439 & bey*(6.d0*xi(jj)**2-6.d0*xi(jj)+aly)/dl*vl(2,1)+ 440 & bey*(3.d0*xi(jj)**2+(aly-4.d0)*xi(jj)+(1.d0-aly/2.d0))*vl(6,1) 441 & +bey*(-6.d0*xi(jj)**2+6.d0*xi(jj)-aly)/dl*vl(2,2) 442 & +bey*(3.d0*xi(jj)**2-(2.d0+aly)*xi(jj)+aly/2.d0)*vl(6,2) 443! 444 vlp(3,jj)=bez*(6.d0*xi(jj)*xi(jj)-6.d0*xi(jj)+alz)/dl*vl(3,1)- 445 & bez*(3.d0*xi(jj)**2+(alz-4.d0)*xi(jj)+(1.d0-alz/2.d0))*vl(5,1) 446 & +bez*(-6.d0*xi(jj)**2+6.d0*xi(jj)-alz)/dl*vl(3,2) 447 & -bez*(3.d0*xi(jj)**2-(2.d0+alz)*xi(jj)+alz/2.d0)*vl(5,2) 448! 449 vlp(4,jj)=(vl(4,2)-vl(4,1))/dl 450! 451 vlp(5,jj)=-6.d0*bez*(2.d0*xi(jj)-1.d0)/(dl*dl)*vl(3,1) 452 & +bez*(6.d0*xi(jj)+(alz-4.d0))/dl*vl(5,1) 453 & -6.d0*bez*(-2.d0*xi(jj)+1.d0)/(dl*dl)*vl(3,2) 454 & +bez*(6.d0*xi(jj)-(alz+2))/dl*vl(5,2) 455! 456 vlp(6,jj)=6.d0*bey*(2.d0*xi(jj)-1.d0)/(dl*dl)*vl(2,1) 457 & +bey*(6.d0*xi(jj)+(aly-4.d0))/dl*vl(6,1) 458 & +6.d0*bey*(-2.d0*xi(jj)+1.d0)/(dl*dl)*vl(2,2) 459 & +bey*(6.d0*xi(jj)-(aly+2.d0))/dl*vl(6,2) 460! 461! calculation of the strains (Eqn. (8) in Luo) 462! 463 eloc(1)=vlp(1,jj) 464 eloc(2)=-vlp(6,jj) 465 eloc(3)=vlp(5,jj) 466 eloc(4)=vlp(2,jj)-vl(6,jj) 467 eloc(5)=vlp(3,jj)+vl(5,jj) 468 eloc(6)=vlp(4,jj) 469! 470! determining the mechanical strain 471! 472 if(ithermal(1).ne.0) then 473 do m1=2,6 474 emec(m1)=eloc(m1)-eth(1) 475 enddo 476 else 477 do m1=1,6 478 emec(m1)=eloc(m1) 479 enddo 480 endif 481! 482! subtracting initial strains 483! 484 if(iprestr.ne.0) then 485 write(*,*) '*ERROR in resultsmech_u1:' 486 write(*,*) ' no initial strains allowed for' 487 write(*,*) ' this user element' 488 call exit(201) 489 endif 490! 491! calculating the section forces 492! simplified version of Eqn. (11) in Luo (symmetric case 493! for which Ay=Az=J=0) 494! 495 stre(1)=e*a*emec(1) 496 stre(2)=e*xi11*emec(2) 497 stre(3)=e*xi22*emec(3) 498 stre(4)=xk*um*a*emec(4) 499 stre(5)=xk*um*a*emec(5) 500 stre(6)=xk*um*(xi11+xi22)*emec(6) 501! 502! rotating the strain into the global system 503! 504 xa(1,1)=eloc(1) 505 xa(1,2)=eloc(4) 506 xa(1,3)=eloc(5) 507 xa(2,1)=eloc(4) 508 xa(2,2)=eloc(2) 509 xa(2,3)=eloc(6) 510 xa(3,1)=eloc(5) 511 xa(3,2)=eloc(6) 512 xa(3,3)=eloc(3) 513! 514 do m1=1,6 515 eloc(m1)=0.d0 516 j1=kal(1,m1) 517 j2=kal(2,m1) 518 do j3=1,3 519 do j4=1,3 520 eloc(m1)=eloc(m1)+ 521 & xa(j3,j4)*tm(j3,j1)*tm(j4,j2) 522 enddo 523 enddo 524 enddo 525! 526! rotating the stress into the global system 527! 528 xa(1,1)=stre(1) 529 xa(1,2)=stre(4) 530 xa(1,3)=stre(5) 531 xa(2,1)=stre(4) 532 xa(2,2)=stre(2) 533 xa(2,3)=stre(6) 534 xa(3,1)=stre(5) 535 xa(3,2)=stre(6) 536 xa(3,3)=stre(3) 537! 538 do m1=1,6 539 stre(m1)=0.d0 540 j1=kal(1,m1) 541 j2=kal(2,m1) 542 do j3=1,3 543 do j4=1,3 544 stre(m1)=stre(m1)+ 545 & xa(j3,j4)*tm(j3,j1)*tm(j4,j2) 546 enddo 547 enddo 548 enddo 549! 550! updating the internal energy and mechanical strain 551! 552 if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100).or. 553 & ((nmethod.eq.4).and.(iperturb(1).gt.1).and. 554 & (ithermal(1).le.1))) then 555! 556 if(nener.eq.1) then 557 ener(jj,i)=enerini(jj,i)+ 558 & ((emec(1)-emeini(1,jj,i))* 559 & (stre(1)+stiini(1,jj,i))+ 560 & (emec(2)-emeini(2,jj,i))* 561 & (stre(2)+stiini(2,jj,i))+ 562 & (emec(3)-emeini(3,jj,i))* 563 & (stre(3)+stiini(3,jj,i)))/2.d0+ 564 & (emec(4)-emeini(4,jj,i))*(stre(4)+stiini(4,jj,i))+ 565 & (emec(5)-emeini(5,jj,i))*(stre(5)+stiini(5,jj,i))+ 566 & (emec(6)-emeini(6,jj,i))*(stre(6)+stiini(6,jj,i)) 567 endif 568 eme(1,jj,i)=emec(1) 569 eme(2,jj,i)=emec(2) 570 eme(3,jj,i)=emec(3) 571 eme(4,jj,i)=emec(4) 572 eme(5,jj,i)=emec(5) 573 eme(6,jj,i)=emec(6) 574 endif 575! 576 if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100)) then 577! 578 eei(1,jj,i)=eloc(1) 579 eei(2,jj,i)=eloc(2) 580 eei(3,jj,i)=eloc(3) 581 eei(4,jj,i)=eloc(4) 582 eei(5,jj,i)=eloc(5) 583 eei(6,jj,i)=eloc(6) 584 endif 585! 586! updating the kinetic energy 587! 588 if(ikin.eq.1) then 589 write(*,*) '*ERROR in resultsmech_u1:' 590 write(*,*) ' no initial strains allowed for' 591 write(*,*) ' this user element' 592 call exit(201) 593 endif 594! 595 stx(1,jj,i)=stre(1) 596 stx(2,jj,i)=stre(2) 597 stx(3,jj,i)=stre(3) 598 stx(4,jj,i)=stre(4) 599 stx(5,jj,i)=stre(5) 600 stx(6,jj,i)=stre(6) 601! 602! calculation of the global nodal forces 603! 604 if(calcul_fn.eq.1)then 605 node=kon(indexe+jj) 606 do j=1,3 607 fn(j,node)=fn(j,node)+tm(1,j)*stre(1) 608 & +tm(2,j)*stre(2) 609 & +tm(3,j)*stre(3) 610 fn(j+3,node)=fn(j+3,node)+tm(1,j)*stre(4) 611 & +tm(2,j)*stre(5) 612 & +tm(3,j)*stre(6) 613 enddo 614 endif 615 enddo 616! 617! q contains the contributions to the nodal force in the nodes 618! belonging to the element at stake from other elements (elements 619! already treated). These contributions have to be 620! subtracted to get the contributions attributable to the element 621! at stake only 622! 623 if(calcul_qa.eq.1) then 624 do m1=1,nope 625 do m2=1,3 626 qa(1)=qa(1)+dabs(fn(m2,konl(m1))-q(m2,m1)) 627 enddo 628 enddo 629 nal=nal+3*nope 630 endif 631! 632 return 633 end 634