q4ke2.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    q4ke2          ../engine/source/elements/solid_2d/quad4/q4ke2.f
!||--- called by ------------------------------------------------------
!||    imp_glob_k0    ../engine/source/implicit/imp_glob_k.F
!||--- calls      -----------------------------------------------------
!||    assem_q4       ../engine/source/implicit/assem_q4.f
!||    mmat_h1        ../engine/source/elements/solid/solide8z/mmat_h1.F
!||    mmats          ../engine/source/elements/solid/solide8z/mmats.F
!||    mmstifs        ../engine/source/elements/solid/solide8z/mmstifs.F
!||    q4deri2        ../engine/source/elements/solid_2d/quad4/q4deri2.F
!||    q4deric2       ../engine/source/elements/solid_2d/quad4/q4deric2.F
!||    q4kega2        ../engine/source/elements/solid_2d/quad4/q4kega2.F
!||    q4kel2         ../engine/source/elements/solid_2d/quad4/q4kel2.F
!||    q4kep2         ../engine/source/elements/solid_2d/quad4/q4kep2.f
!||    q4kerot2       ../engine/source/elements/solid_2d/quad4/q4kerot2.F
!||    qcoor2         ../engine/source/elements/solid_2d/quad/qcoor2.F
!||    qrcoor2        ../engine/source/elements/solid_2d/quad/qrcoor2.F
!||    qvolu2         ../engine/source/elements/solid_2d/quad/qvolu2.F
!||    s8zsigp3       ../engine/source/elements/solid/solide8z/s8zsigp3.F
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod   ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    element_mod    ../common_source/modules/elements/element_mod.F90
!||====================================================================
      SUBROUTINE q4ke2(
     1   PM,       GEO,      IXQ,      X,
     2   ELBUF_STR,NEL,      LIAD,     ICP,
     3   ICSIG,    ETAG,     IDDL,     NDOF,
     4   K_DIAG,   K_LT,     IADK,     JDIK,
     5   NPG,      IPM,      IGEO,     IKGEO,
     6   BUFMAT,   NFT,      MTN,      JMULT,
     7   JHBE,     JCVT,     IGTYP,    ISORTH,
     8   ISMSTR)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE elbufdef_mod            
      use element_mod , only : nixq
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: ISMSTR
      INTEGER, INTENT(IN) :: NFT
      INTEGER, INTENT(IN) :: MTN
      INTEGER, INTENT(IN) :: JMULT
      INTEGER, INTENT(IN) :: JHBE
      INTEGER, INTENT(IN) :: JCVT
      INTEGER, INTENT(IN) :: IGTYP
      INTEGER, INTENT(IN) :: ISORTH
      INTEGER IXQ(NIXQ,*), ICP, ICSIG, IKGEO
      INTEGER NEL, LIAD, NPG,
     .   ipm(npropmi,*), igeo(npropgi,*), etag(*), iddl(*),
     .   ndof(*), iadk(*), jdik(*)
      my_real
     .   pm(npropm,*), geo(npropg,*), x(3,*),
     .   bufmat(*), k_diag(*), k_lt(*)
      TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_STR
C-----------------------------------------------
c FUNCTION: elemental stiffness of Q4 element
c ARGUMENTS:  (I: input, O: output, IO: input & output, W: workspace)
c TYPE NAME                FUNCTION
c  I   PM(),GEO()          material and geometrical property data
c  I   IXQ(7,NUM_QUAD)     connectivity and mid,pid integer data
c  I   X(3,NUMNOD)         coordinates
c  I   NEL                 number of quad element in this group
c  I   ICP                 flag for constant pressure
c  I   ^ICSG                flag for solid shell usage
c  I   ^ETAG()              for eigein value computation usage
c  I   IDDL()              position for nodal DOFs
c  I   NDOF()              number of nodal DOFs
c  I   K_DIAG()            diagonal K
c  I   K_LT()              off-diagonal K
c  I   IADK()              position for K
c  I   JDIK()              position for K
c  I   NPG                 flag for number of integration points
c  I   IPM(NPROPMI,*)      material property data (integer)
c  I   IGEO(NPROPGI,*)     geometrical property data (integer)
c  I   IKGEO               flag for geometric stiffness computation method
c  I   BUFMAT()            buffer for material constitution computation
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J, EP
C
      INTEGER IKORTH,IADBUF,ICPG,IPREDU
C
      INTEGER IAD0
C     INDEX OF ELEMENT INFORMATION IN GLOBAL TABLE "IXQ"
      INTEGER LCO
C     SN OF THE FIRST ELEMENT OF THE GROUP IN GLOBAL STORAGE
      INTEGER NF1
C     MATERIAL SN, CONNECTIVITY, ID, PROPERTY SN OF ELEMENTS
      INTEGER MXT(MVSIZ),
     +    NC1(MVSIZ),NC2(MVSIZ),NC3(MVSIZ),NC4(MVSIZ),
     +    ngl(mvsiz),ngeo(mvsiz)
C     NUMBER AND INDEXES OF INTEGRATION POINTS
      INTEGER NNPT,NPTR,NPTS,IR,IS,IT,IP
C
C     ELEMENT D/A FLAG
C     NODAL COORDINATES (t+i)
C     DIFFERENCES OF "Y", "Z"
C     SUMMERIES OF "Y"
C     SHAPE DERIVATIVES (dNi/dY, dNi/dZ, Ni/r) AT CENTER
C     AREA, VOLUME/(THICKNESS .OR. 2*PI)
C     SHAPE DERIVATIVES (dNi/dY, dNi/dZ) AT INTEGRATION POINT
C     (W*|J| .OR. r'* W*|J|) AT INTEGRATION POINT
C     STIFFNESS MATRIX
C     TRANSFORMATION MATRIX [R] FOR CO-ROTATIONAL CASE
C     {X}=[R]{X'} <=> {X'}=T([R]){X}
C     [K]=[R][K']T([R])
      my_real
     +    offg(mvsiz),off(mvsiz),gama(mvsiz,6),
     +    y1(mvsiz),y2(mvsiz),y3(mvsiz),y4(mvsiz),
     +    z1(mvsiz),z2(mvsiz),z3(mvsiz),z4(mvsiz),
     +    y12(mvsiz),y34(mvsiz),y13(mvsiz),y24(mvsiz),
     +    y14(mvsiz),y23(mvsiz),
     +    z12(mvsiz),z34(mvsiz),z13(mvsiz),z24(mvsiz),
     +    z14(mvsiz),z23(mvsiz),
     +    y234(mvsiz),y124(mvsiz),yavg(mvsiz),
     +    pyc1(mvsiz),pyc2(mvsiz),pzc1(mvsiz),pzc2(mvsiz),
     + ay(mvsiz),
     +    aire(mvsiz),volu(mvsiz),
     +    py1(mvsiz),py2(mvsiz),py3(mvsiz),py4(mvsiz),
     +    pz1(mvsiz),pz2(mvsiz),pz3(mvsiz),pz4(mvsiz),
     +    airn(mvsiz),voln(mvsiz),
     +    k11(2,2,mvsiz),k12(2,2,mvsiz),k13(2,2,mvsiz),k14(2,2,mvsiz),
     +    k22(2,2,mvsiz),k23(2,2,mvsiz),k24(2,2,mvsiz),
     +    k33(2,2,mvsiz),k34(2,2,mvsiz),k44(2,2,mvsiz),
     + k11u(2,2,mvsiz),k12u(2,2,mvsiz),k13u(2,2,mvsiz),k14u(2,2,mvsiz),
     + k22u(2,2,mvsiz),k23u(2,2,mvsiz),k24u(2,2,mvsiz),
     + k33u(2,2,mvsiz),k34u(2,2,mvsiz),k44u(2,2,mvsiz),
     + k11l(2,2,mvsiz),k12l(2,2,mvsiz),k13l(2,2,mvsiz),k14l(2,2,mvsiz),
     + k22l(2,2,mvsiz),k23l(2,2,mvsiz),k24l(2,2,mvsiz),
     + k33l(2,2,mvsiz),k34l(2,2,mvsiz),k44l(2,2,mvsiz),
     +    r11(mvsiz),r12(mvsiz),r13(mvsiz),
     +    r21(mvsiz),r22(mvsiz),r23(mvsiz),
     +    r31(mvsiz),r32(mvsiz),r33(mvsiz)
C
C     MATERIAL PARAMETERS
      my_real
     .    nu(mvsiz),c1,e0(mvsiz),fac(mvsiz),bid(1),
     .    hh(2,mvsiz),hh1(2,mvsiz),
     .    dm(9,mvsiz),dgm(9,mvsiz),gm(9,mvsiz),
     .    dd(9,mvsiz),dg(9,mvsiz),gg(mvsiz),g33(9,mvsiz),
     .    byz1(mvsiz),byz2(mvsiz),byz3(mvsiz),byz4(mvsiz),
     .    bzy1(mvsiz),bzy2(mvsiz),bzy3(mvsiz),bzy4(mvsiz),
     +    nuu(mvsiz)
C
C     USED ONLY FOR SUBROUTINE INTERFACE
      my_real
     +    vd2(mvsiz),vis(mvsiz),
     +    rx(mvsiz),ry(mvsiz),rz(mvsiz),
     +    sx(mvsiz),sy(mvsiz),sz(mvsiz)
C
      my_real
     +  wi,ksi,eta
      my_real
     +  w_gauss(9,9),a_gauss(9,9)
      TYPE(g_bufel_)  ,POINTER :: GBUF     
      TYPE(L_BUFEL_)  ,POINTER :: LBUF     
c
      DATA W_GAUSS / 
     1 2.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     2 1.               ,1.               ,0.               ,
     2 0.               ,0.               ,0.               ,
     2 0.               ,0.               ,0.               ,
     3 0.555555555555556,0.888888888888889,0.555555555555556,
     3 0.               ,0.               ,0.               ,
     3 0.               ,0.               ,0.               ,
     4 0.347854845137454,0.652145154862546,0.652145154862546,
     4 0.347854845137454,0.               ,0.               ,
     4 0.               ,0.               ,0.               ,
     5 0.236926885056189,0.478628670499366,0.568888888888889,
     5 0.478628670499366,0.236926885056189,0.               ,
     5 0.               ,0.               ,0.               ,
     6 0.171324492379170,0.360761573048139,0.467913934572691,
     6 0.467913934572691,0.360761573048139,0.171324492379170,
     6 0.               ,0.               ,0.               ,
     7 0.129484966168870,0.279705391489277,0.381830050505119,
     7 0.417959183673469,0.381830050505119,0.279705391489277,
     7 0.129484966168870,0.               ,0.               ,
     8 0.101228536290376,0.222381034453374,0.313706645877887,
     8 0.362683783378362,0.362683783378362,0.313706645877887,
     8 0.222381034453374,0.101228536290376,0.               ,
     9 0.081274388361574,0.180648160694857,0.260610696402935,
     9 0.312347077040003,0.330239355001260,0.312347077040003,
     9 0.260610696402935,0.180648160694857,0.081274388361574/
      DATA a_gauss / 
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     2 -.577350269189626,0.577350269189626,0.               ,
     2 0.               ,0.               ,0.               ,
     2 0.               ,0.               ,0.               ,
     3 -.774596669241483,0.               ,0.774596669241483,
     3 0.               ,0.               ,0.               ,
     3 0.               ,0.               ,0.               ,
     4 -.861136311594053,-.339981043584856,0.339981043584856,
     4 0.861136311594053,0.               ,0.               ,
     4 0.               ,0.               ,0.               ,
     5 -.906179845938664,-.538469310105683,0.               ,
     5 0.538469310105683,0.906179845938664,0.               ,
     5 0.               ,0.               ,0.               ,
     6 -.932469514203152,-.661209386466265,-.238619186083197,
     6 0.238619186083197,0.661209386466265,0.932469514203152,
     6 0.               ,0.               ,0.               ,
     7 -.949107912342759,-.741531185599394,-.405845151377397,
     7 0.               ,0.405845151377397,0.741531185599394,
     7 0.949107912342759,0.               ,0.               ,
     8 -.960289856497536,-.796666477413627,-.525532409916329,
     8 -.183434642495650,0.183434642495650,0.525532409916329,
     8 0.796666477413627,0.960289856497536,0.               ,
     9 -.968160239507626,-.836031107326636,-.613371432700590,
     9 -.324253423403809,0.               ,0.324253423403809,
     9 0.613371432700590,0.836031107326636,0.968160239507626/
C-----------------------------------------------
C   S o u r c e  L i n e s
C-----------------------------------------------
      gbuf => elbuf_str%GBUF
      IF (isorth == 0) THEN
        DO i=1,nel
          gama(i,1) = one
          gama(i,2) = zero
          gama(i,3) = zero
          gama(i,4) = zero
          gama(i,5) = one
          gama(i,6) = zero
        ENDDO
      ELSE
        DO i=1,nel                                            
          gama(i,1) = gbuf%GAMA(i        )
          gama(i,2) = gbuf%GAMA(i +   nel)
          gama(i,3) = gbuf%GAMA(i + 2*nel)
          gama(i,4) = gbuf%GAMA(i + 3*nel)
          gama(i,5) = gbuf%GAMA(i + 4*nel)
          gama(i,6) = gbuf%GAMA(i + 5*nel)
        ENDDO
      ENDIF
      iad0 = 1
      IF (isorth > 0) iad0 = 1 + 6*nel
      IF (igtyp == 21.OR.igtyp == 22) THEN
       ikorth=2
      ELSEIF (isorth>0) THEN
       ikorth=1
      ELSE
       ikorth=0
      ENDIF
C
      lco = 1 + nixq*nft
      nf1 = 1 + nft
C
C     GATHER NODAL INFORMATION & 
C     COMPUTE INTRINSIC ROTATION FOR CO-ROTATIONAL CASE & 
C     PROJECT NODAL COORDINATES INTO CO-ROTATIONAL SYSTEM
      IF (jcvt==0) THEN
        CALL qcoor2(
     1   x,         ixq(1,nf1),y1,        y2,
     2   y3,        y4,        z1,        z2,
     3   z3,        z4,        nc1,       nc2,
     4   nc3,       nc4,       ngl,       mxt,
     5   ngeo,      vd2,       vis,       nel)
      ELSE
        CALL qrcoor2(
     1   x,         ixq(1,nf1),y1,        y2,
     2   y3,        y4,        z1,        z2,
     3   z3,        z4,        nc1,       nc2,
     4   nc3,       nc4,       ngl,       mxt,
     5   ngeo,      vd2,       r11,       r12,
     6   r13,       r21,       r22,       r23,
     7   r31,       r32,       r33,       gama,
     8   y234,      y124,      vis,       nel,
     9   isorth)
      ENDIF
C
C     COMPUTE FAC(*) FOR PRESSURE ACCORDING TO PLASTICITY STATE
C---- now assumed strain is used, no effect for Isolid17
        DO i=1,nel
          nu(i)=min(half,pm(21,mxt(i)))
          c1 =pm(32,mxt(i))
          e0(i) =three*(one-two*nu(i))*c1
        ENDDO
      IF(icp==2) THEN
        CALL s8zsigp3(1  ,nel       ,gbuf%SIG,e0  ,gbuf%PLA,
     2                fac  ,gbuf%G_PLA,nel     )
        DO i=1,nel
          nuu(i)=nu(i)+(half-nu(i))*fac(i)
        ENDDO
      ELSEIF(icp==1) THEN
        DO i=1,nel
          nuu(i)=half
        ENDDO
      ELSE
        DO i=1,nel
          nuu(i)=zero
        ENDDO
      ENDIF 
C     COMPUTE AREA & VOLUME/THICKNESS
      CALL qvolu2(
     1   gbuf%OFF,aire,    volu,    ngl,
     2   y1,      y2,      y3,      y4,
     3   z1,      z2,      z3,      z4,
     4   y234,    y124,    nel,     jmult,
     5   jcvt)
C
      IF(n2d==1) THEN
        DO i=1,nel
         yavg(i) = x(2,nc1(i))+x(2,nc2(i))+x(2,nc3(i))+x(2,nc4(i))
        ENDDO
      ENDIF
C     COMPUTE SHAPE DERIVATIVES AT ELEMENT CENTER
      CALL q4deric2(
     1   y1,      y2,      y3,      y4,
     2   z1,      z2,      z3,      z4,
     3   y12,     y34,     y13,     y24,
     4   y14,     y23,     z12,     z34,
     5   z13,     z24,     z14,     z23,
     6   pyc1,    pyc2,    pzc1,    pzc2,
     7   aire,    volu,    yavg,    rx,
     8   ry,      rz,      sx,      sy,
     9   sz,      nel,     jhbe)
C
C     COMPUTE VOLUME/(2*PI) FOR AXISYMMETRIC CASE
C     
C
      icpg = 0
      IF(icpg==2) icpg = 1
C
C
      nptr = 2
      npts = 2
      nnpt = nptr*npts
C
      IF (mtn>=28) THEN
        iadbuf = ipm(7,mxt(1))
      ELSE
        iadbuf = 1
      ENDIF
      CALL mmats(1    ,nel     ,pm    ,mxt    ,hh    ,
     .           mtn    ,ikorth  ,ipm   ,igeo   ,gama  ,
     .           bufmat(iadbuf)  ,dm    ,dgm    ,gm    ,
     .           jhbe  ,gbuf%SIG ,bid   ,nnpt   ,nel   )
      CALL mmat_h1(
     1   hh,      hh1,     fac,     icpg,
     2   ipredu,  nel,     mtn,     ismstr,
     3   jhbe)
C
      DO i=1,nel
        offg(i) = gbuf%OFF(i)
      ENDDO
C
      DO ep=1,nel
      DO j=1,2
      DO i=1,2
        k11(i,j,ep)=zero
        k12(i,j,ep)=zero
        k13(i,j,ep)=zero
        k14(i,j,ep)=zero
        k22(i,j,ep)=zero
        k23(i,j,ep)=zero
        k24(i,j,ep)=zero
        k33(i,j,ep)=zero
        k34(i,j,ep)=zero
        k44(i,j,ep)=zero
        k11u(i,j,ep)=zero
        k12u(i,j,ep)=zero
        k13u(i,j,ep)=zero
        k14u(i,j,ep)=zero
        k22u(i,j,ep)=zero
        k23u(i,j,ep)=zero
        k24u(i,j,ep)=zero
        k33u(i,j,ep)=zero
        k34u(i,j,ep)=zero
        k44u(i,j,ep)=zero
        k11l(i,j,ep)=zero
        k12l(i,j,ep)=zero
        k13l(i,j,ep)=zero
        k14l(i,j,ep)=zero
        k22l(i,j,ep)=zero
        k23l(i,j,ep)=zero
        k24l(i,j,ep)=zero
        k33l(i,j,ep)=zero
        k34l(i,j,ep)=zero
        k44l(i,j,ep)=zero
      ENDDO
      ENDDO
      ENDDO
C
C     ENTER THE INTEGRATION POINTS LOOP -->
      it = 1
      DO 100 ir=1,nptr
      DO 200 is=1,npts
        lbuf => elbuf_str%BUFLY(1)%LBUF(ir,is,1)
C
C     INITIALIZE WEIGHTING FACTORS
      ip = ir + (is-1)*nptr
      ksi = a_gauss(ir,nptr)
      eta = a_gauss(is,npts)
      wi = w_gauss(ir,nptr)*w_gauss(is,npts)
C
C     COMPUTE JACOBIAN & SHAPE DERIVATIVES AT INTEGRATION POINT
        CALL q4deri2(
     1   offg,    off,     ksi,     eta,
     2   wi,      yavg,    y12,     y34,
     3   y13,     y24,     y14,     y23,
     4   z12,     z34,     z13,     z24,
     5   z14,     z23,     py1,     py2,
     6   py3,     py4,     pz1,     pz2,
     7   pz3,     pz4,     pyc1,    pyc2,
     8   pzc1,    pzc2,    byz1,    byz2,
     9   byz3,    byz4,    bzy1,    bzy2,
     a   bzy3,    bzy4,    airn,    voln,
     b   nuu,     nel,     jhbe)
C-------------add BYZj... later
C
      CALL mmstifs(
     1   pm,      mxt,     hh1,     voln,
     2   icsig,   dd,      gg,      dg,
     3   g33,     dm,      gm,      dgm,
     4   ikorth,  lbuf%SIG,ir,      is,
     5   it,      nel,     jhbe,    mtn)
C
      CALL q4kel2(
     1   py1,     py2,     py3,     py4,
     2   pz1,     pz2,     pz3,     pz4,
     3   pyc1,    pyc2,    pzc1,    pzc2,
     4   ay,      r22,     r23,     k11,
     5   k12,     k13,     k14,     k22,
     6   k23,     k24,     k33,     k34,
     7   k44,     k11u,    k12u,    k13u,
     8   k14u,    k22u,    k23u,    k24u,
     9   k33u,    k34u,    k44u,    k11l,
     a   k12l,    k13l,    k14l,    k22l,
     b   k23l,    k24l,    k33l,    k34l,
     c   k44l,    dd,      gg,      dg,
     d   g33,     ikorth,  icpg,    offg,
     e   nel,     jcvt)
C
c      IF (IKGEO<0) THEN
c        CALL Q4KEG2(
c     1              PY1,PY2,PY3,PY4,
c     2              PZ1,PZ2,PZ3,PZ4,AY,
c     3              K11,K12,K13,K14,K22,
c     4              K23,K24,K33,K34,K44,
c     5              EV(NB2),VOLN,OFFG)
c      ENDIF
C
200   CONTINUE
100   CONTINUE
C     EXIT THE INTEGRATION POINTS LOOP <--
C
      IF (ipredu > 0) THEN
        CALL q4kep2(
     1   pyc1,    pyc2,    pzc1,    pzc2,
     2   ay,      r22,     r23,     k11,
     3   k12,     k13,     k14,     k22,
     4   k23,     k24,     k33,     k34,
     5   k44,     hh,      volu,    fac,
     6   icpg,    offg,    nel,     jcvt)
      ENDIF
C
c      IF (IKGEO>0) THEN
      IF (ikgeo/=0) THEN
        CALL q4kega2(
     1   pyc1,    pyc2,    pzc1,    pzc2,
     2   ay,      k11,     k12,     k13,
     3   k14,     k22,     k23,     k24,
     4   k33,     k34,     k44,     gbuf%SIG,
     5   volu,    offg,    nel)
      ENDIF
C
c      IF (N2D==1 .AND. JHBE==17) THEN
c        CALL Q4ASYM(1, NEL, K11, K11U, K11L)
c        CALL Q4ASYM(1, NEL, K22, K22U, K22L)
c        CALL Q4ASYM(1, NEL, K33, K33U, K33L)
c        CALL Q4ASYM(1, NEL, K44, K44U, K44L)
c        CALL Q4ASYM(1, NEL, K12, K12U, K12L)
c        CALL Q4ASYM(1, NEL, K13, K13U, K13L)
c        CALL Q4ASYM(1, NEL, K14, K14U, K14L)
c        CALL Q4ASYM(1, NEL, K23, K23U, K23L)
c        CALL Q4ASYM(1, NEL, K24, K24U, K24L)
c        CALL Q4ASYM(1, NEL, K34, K34U, K34L)
c      ENDIF
C
      IF (jcvt/=0) THEN
        CALL q4kerot2(
     1   r22,     r32,     r23,     r33,
     2   k11,     k12,     k13,     k14,
     3   k22,     k23,     k24,     k33,
     4   k34,     k44,     nel)
      ENDIF
C
C     IF (NEIG>0) THEN
C       CALL Q4EOFF(1, NEL, IXQ(1,NF1), ETAG, OFFG)
C     ENDIF
C
      CALL assem_q4(
     1    ixq(1,nf1),nel   ,iddl  ,ndof  ,k_diag,
     2    k_lt  ,iadk  ,jdik  ,k11   ,k12   ,
     3    k13   ,k14   ,k22   ,k23   ,k24   ,
     4    k33   ,k34   ,k44   ,offg  )
C
      RETURN
      END