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 mafillsmcs(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
20     &     xboun,nboun,
21     &     ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
22     &     nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
23     &     ad,au,bb,nactdof,icol,jq,irow,neq,nzl,nmethod,
24     &     ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,
25     &     nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab,
26     &     ntmat_,t0,t1,ithermal,prestr,
27     &     iprestr,vold,iperturb,sti,nzs,stx,adb,aub,iexpl,plicon,
28     &     nplicon,plkcon,nplkcon,xstiff,npmat_,dtime,
29     &     matname,mi,ics,cs,nm,ncmat_,labmpc,mass,stiffness,buckling,
30     &     rhsi,intscheme,mcs,coriolis,ibody,xloadold,reltime,ielcs,
31     &     veold,springarea,thicke,integerglob,doubleglob,tieset,
32     &     istartset,iendset,ialset,ntie,nasym,pslavsurf,pmastsurf,
33     &     mortar,clearini,ielprop,prop,ne0,kscale,xstateini,xstate,
34     &     nstate_,set,nset)
35!
36!     filling the stiffness matrix in spare matrix format (sm)
37!     for cyclic symmetry calculations
38!
39      implicit none
40!
41      logical mass,stiffness,buckling,rhsi,coriolis
42!
43      character*8 lakon(*)
44      character*20 labmpc(*),sideload(*)
45      character*80 matname(*)
46      character*81 tieset(3,*),set(*)
47!
48      integer kon(*),nodeboun(*),ndirboun(*),ipompc(*),nodempc(3,*),
49     &     nodeforc(2,*),ndirforc(*),nelemload(2,*),icol(*),jq(*),
50     &     ilmpc(*),ikboun(*),ilboun(*),mi(*),nstate_,ne0,ielprop(*),
51     &     nactdof(0:mi(2),*),irow(*),istartset(*),iendset(*),kscale,
52     &     nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),nset,
53     &     ielorien(mi(3),*),integerglob(*),ialset(*),ntie,ikmpc(*),
54     &     ipkon(*),ics(*),ij,ilength,lprev,ipobody(2,*),nbody,
55     &     ibody(3,*),nk,ne,nboun,nmpc,nforc,nload,neq,nzl,nmethod,
56     &     ithermal(*),iprestr,iperturb(*),nzs,i,j,k,l,m,idist,jj,
57     &     ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2,idof1,idof2,
58     &     mpc1,mpc2,index1,index2,node1,node2,kflag,nasym,mortar,
59     &     ntmat_,indexe,nope,norien,iexpl,i0,nm,inode,icomplex,
60     &     inode1,icomplex1,inode2,icomplex2,ner,ncmat_,intscheme,istep,
61     &     iinc,mcs,ielcs(*),nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),
62     &     npmat_,islavelinv(1),irowtloc(1),jqtloc(1),mortartrafoflag,
63     &     mscalmethod
64!
65      real*8 co(3,*),xboun(*),coefmpc(*),xforc(*),xload(2,*),p1(3),
66     &     p2(3),ad(*),au(*),bodyf(3),bb(*),xbody(7,*),cgr(4,*),prop(*),
67     &     t0(*),t1(*),prestr(6,mi(1),*),vold(0:mi(2),*),s(60,60),
68     &     xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),ff(60),
69     &     sti(6,mi(1),*),sm(60,60),stx(6,mi(1),*),adb(*),aub(*),
70     &     elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),xloadold(2,*),
71     &     alcon(0:6,ntmat_,*),cs(17,*),alzero(*),orab(7,*),reltime,
72     &     springarea(2,*),plicon(0:2*npmat_,ntmat_,*),smscale(1),
73     &     plkcon(0:2*npmat_,ntmat_,*),thicke(mi(3),*),doubleglob(*),
74     &     xstiff(27,mi(1),*),pi,theta,ti,tr,veold(0:mi(2),*),om,valu2,
75     &     value,dtime,walue,walu2,time,ttime,clearini(3,9,*),
76     &     pslavsurf(3,*),pmastsurf(6,*),autloc(1)
77!
78      mortartrafoflag=0
79!
80!     calculating the scaling factors for the cyclic symmetry calculation
81!
82      pi=4.d0*datan(1.d0)
83!
84      do i=1,mcs
85        write(*,*) '*INFO in mafillsmcs: calculating nodal diameter',nm
86        write(*,*) '      for',cs(1,i),'sectors'
87        write(*,*) '      (cyclic symmetry definition number',i,')'
88        write(*,*)
89        theta=nm*2.d0*pi/cs(1,i)
90        cs(15,i)=dcos(theta)
91        cs(16,i)=dsin(theta)
92      enddo
93!
94      kflag=2
95      i0=0
96!
97!     determining nzl
98!
99      nzl=0
100      do i=neq,1,-1
101        if(icol(i).gt.0) then
102          nzl=i
103          exit
104        endif
105      enddo
106!
107!     initializing the matrices
108!
109      do i=1,neq
110        ad(i)=0.d0
111      enddo
112      do i=1,nzs
113        au(i)=0.d0
114      enddo
115!
116      do i=1,neq
117        adb(i)=0.d0
118      enddo
119      do i=1,nzs
120        aub(i)=0.d0
121      enddo
122!
123      ner=neq/2
124!
125!     loop over all elements
126!
127!     initialisation of the error parameter
128!
129c     ne0=0
130      do i=1,ne
131!
132        if(ipkon(i).lt.0) cycle
133        indexe=ipkon(i)
134c     Bernhardi start
135        if(lakon(i)(4:5).eq.'8I') then
136          nope=11
137c     Bernhardi end
138        elseif(lakon(i)(4:5).eq.'20') then
139          nope=20
140        elseif(lakon(i)(4:4).eq.'2') then
141          nope=26
142        elseif(lakon(i)(4:4).eq.'8') then
143          nope=8
144        elseif(lakon(i)(4:5).eq.'10') then
145          nope=10
146        elseif(lakon(i)(4:4).eq.'4') then
147          nope=4
148        elseif(lakon(i)(4:5).eq.'15') then
149          nope=15
150        elseif(lakon(i)(4:4).eq.'6') then
151          nope=6
152        elseif(lakon(i)(1:2).eq.'ES') then
153c     begin 16.07.2020
154          nope=ichar(lakon(i)(8:8))-47
155c     read(lakon(i)(8:8),'(i1)') nope
156c     nope=nope+1
157c     end 16.07.2020
158!
159!     local contact spring number
160!
161c     write(*,*) 'nope before= ',nope
162          if(lakon(i)(7:7).eq.'C') then
163            if(nasym.eq.1) cycle
164c     begin 16.07.2020
165            if(mortar.eq.1) nope=kon(indexe)
166c     end 16.07.2020
167          endif
168c     write(*,*) 'nope after= ',nope
169        elseif(lakon(i)(1:4).eq.'MASS') then
170          nope=1
171        else
172          cycle
173        endif
174!
175        om=0.d0
176!
177        if((nbody.gt.0).and.(lakon(i)(1:1).ne.'E')) then
178!
179!     assigning centrifugal forces
180!
181          index=i
182          do
183            j=ipobody(1,index)
184            if(j.eq.0) exit
185            if(ibody(1,j).eq.1) then
186              om=xbody(1,j)
187              p1(1)=xbody(2,j)
188              p1(2)=xbody(3,j)
189              p1(3)=xbody(4,j)
190              p2(1)=xbody(5,j)
191              p2(2)=xbody(6,j)
192              p2(3)=xbody(7,j)
193            endif
194            index=ipobody(2,index)
195            if(index.eq.0) exit
196          enddo
197        endif
198!
199        call e_c3d(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,ff,i,
200     &       nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
201     &       alzero,ielmat,ielorien,norien,orab,ntmat_,
202     &       t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload,
203     &       nload,idist,sti,stx,iexpl,plicon,
204     &       nplicon,plkcon,nplkcon,xstiff,npmat_,
205     &       dtime,matname,mi(1),ncmat_,mass,stiffness,buckling,rhsi,
206     &       intscheme,ttime,time,istep,iinc,coriolis,xloadold,
207     &       reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
208     &       springarea,nstate_,xstateini,xstate,ne0,ipkon,thicke,
209     &       integerglob,doubleglob,tieset,istartset,
210     &       iendset,ialset,ntie,nasym,pslavsurf,pmastsurf,mortar,
211     &       clearini,ielprop,prop,kscale,smscale(1),mscalmethod,
212     &       set,nset,islavelinv,autloc,
213     &       irowtloc,jqtloc,mortartrafoflag)
214!
215        do jj=1,3*nope
216!
217          j=(jj-1)/3+1
218          k=jj-3*(j-1)
219!
220          node1=kon(indexe+j)
221          jdof1=nactdof(k,node1)
222!
223          do ll=jj,3*nope
224            if (mcs.gt.1)then
225              if(ielcs(i).gt.0) then
226                s(jj,ll)=(cs(1,(ielcs(i)+1))/cs(1,1))*s(jj,ll)
227                sm(jj,ll)=(cs(1,(ielcs(i)+1))/cs(1,1))*sm(jj,ll)
228              endif
229            endif
230!
231            l=(ll-1)/3+1
232            m=ll-3*(l-1)
233!
234            node2=kon(indexe+l)
235            jdof2=nactdof(m,node2)
236!
237!     check whether one of the DOF belongs to a SPC or MPC
238!
239            if((jdof1.gt.0).and.(jdof2.gt.0)) then
240              call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1,jdof2,
241     &             s(jj,ll),sm(jj,ll),jj,ll)
242              call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1+ner,jdof2+ner,
243     &             s(jj,ll),sm(jj,ll),jj,ll)
244            elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
245!
246!     idof1: genuine DOF
247!     idof2: nominal DOF of the SPC/MPC
248!
249              if(jdof1.le.0) then
250                idof1=jdof2
251                idof2=jdof1
252              else
253                idof1=jdof1
254                idof2=jdof2
255              endif
256!
257              if(nmpc.gt.0) then
258                if(idof2.ne.2*(idof2/2)) then
259!
260!     regular DOF / MPC
261!
262                  id1=(-idof2+1)/2
263                  ist=ipompc(id1)
264                  index=nodempc(3,ist)
265                  if(index.eq.0) cycle
266                  do
267                    inode=nodempc(1,index)
268                    icomplex=0
269                    if(labmpc(id1)(1:6).eq.'CYCLIC') then
270                      read(labmpc(id1)(7:20),'(i14)') icomplex
271                    elseif(labmpc(id1)(1:9).eq.'SUBCYCLIC') then
272                      do ij=1,mcs
273                        ilength=int(cs(4,ij))
274                        lprev=int(cs(14,ij))
275                        call nident(ics(lprev+1),inode,ilength,id)
276                        if(id.gt.0) then
277                          if(ics(lprev+id).eq.inode) then
278                            icomplex=ij
279                            exit
280                          endif
281                        endif
282                      enddo
283                    endif
284                    idof2=nactdof(nodempc(2,index),inode)
285                    if(idof2.gt.0) then
286                      value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
287                      valu2=-coefmpc(index)*sm(jj,ll)/
288     &                     coefmpc(ist)
289                      if(idof1.eq.idof2) then
290                        value=2.d0*value
291                        valu2=2.d0*valu2
292                      endif
293                      if(icomplex.eq.0) then
294                        call add_sm_ei(au,ad,aub,adb,jq,irow,
295     &                       idof1,idof2,value,valu2,i0,i0)
296                        call add_sm_ei(au,ad,aub,adb,jq,irow,
297     &                       idof1+ner,idof2+ner,value,valu2,i0,i0)
298                      else
299                        walue=value*cs(15,icomplex)
300                        walu2=valu2*cs(15,icomplex)
301                        call add_sm_ei(au,ad,aub,adb,jq,irow,
302     &                       idof1,idof2,walue,walu2,i0,i0)
303                        call add_sm_ei(au,ad,aub,adb,jq,irow,
304     &                       idof1+ner,idof2+ner,walue,walu2,i0,i0)
305                        if(idof1.ne.idof2) then
306                          walue=value*cs(16,icomplex)
307                          walu2=valu2*cs(16,icomplex)
308                          call add_sm_ei(au,ad,aub,adb,jq,irow,
309     &                         idof1,idof2+ner,walue,walu2,i0,i0)
310                          walue=-walue
311                          walu2=-walu2
312                          call add_sm_ei(au,ad,aub,adb,jq,irow,
313     &                         idof1+ner,idof2,walue,walu2,i0,i0)
314                        endif
315                      endif
316                    endif
317                    index=nodempc(3,index)
318                    if(index.eq.0) exit
319                  enddo
320                  cycle
321                endif
322              endif
323!
324            else
325              idof1=jdof1
326              idof2=jdof2
327!
328              mpc1=0
329              mpc2=0
330              if(nmpc.gt.0) then
331                if(idof1.ne.2*(idof1/2)) mpc1=1
332                if(idof2.ne.2*(idof2/2)) mpc2=1
333              endif
334              if((mpc1.eq.1).and.(mpc2.eq.1)) then
335                id1=(-idof1+1)/2
336                id2=(-idof2+1)/2
337                if(id1.eq.id2) then
338!
339!     MPC id1 / MPC id1
340!
341                  ist=ipompc(id1)
342                  index1=nodempc(3,ist)
343                  if(index1.eq.0) cycle
344                  do
345                    inode1=nodempc(1,index1)
346                    icomplex1=0
347                    if(labmpc(id1)(1:6).eq.'CYCLIC') then
348                      read(labmpc(id1)(7:20),'(i14)') icomplex1
349                    elseif(labmpc(id1)(1:9).eq.'SUBCYCLIC') then
350                      do ij=1,mcs
351                        ilength=int(cs(4,ij))
352                        lprev=int(cs(14,ij))
353                        call nident(ics(lprev+1),inode1,
354     &                       ilength,id)
355                        if(id.gt.0) then
356                          if(ics(lprev+id).eq.inode1) then
357                            icomplex1=ij
358                            exit
359                          endif
360                        endif
361                      enddo
362                    endif
363                    idof1=nactdof(nodempc(2,index1),inode1)
364                    index2=index1
365                    do
366                      inode2=nodempc(1,index2)
367                      icomplex2=0
368                      if(labmpc(id1)(1:6).eq.'CYCLIC') then
369                        read(labmpc(id1)(7:20),'(i14)') icomplex2
370                      elseif(labmpc(id1)(1:9).eq.'SUBCYCLIC') then
371                        do ij=1,mcs
372                          ilength=int(cs(4,ij))
373                          lprev=int(cs(14,ij))
374                          call nident(ics(lprev+1),inode2,
375     &                         ilength,id)
376                          if(id.gt.0) then
377                            if(ics(lprev+id).eq.inode2) then
378                              icomplex2=ij
379                              exit
380                            endif
381                          endif
382                        enddo
383                      endif
384                      idof2=nactdof(nodempc(2,index2),inode2)
385                      if((idof1.gt.0).and.(idof2.gt.0)) then
386                        value=coefmpc(index1)*coefmpc(index2)*
387     &                       s(jj,ll)/coefmpc(ist)/coefmpc(ist)
388                        valu2=coefmpc(index1)*coefmpc(index2)*
389     &                       sm(jj,ll)/coefmpc(ist)/coefmpc(ist)
390                        if((icomplex1.eq.0).and.
391     &                       (icomplex2.eq.0)) then
392                          call add_sm_ei(au,ad,aub,adb,jq,
393     &                         irow,idof1,idof2,value,valu2,i0,i0)
394                          call add_sm_ei(au,ad,aub,adb,jq,
395     &                         irow,idof1+ner,idof2+ner,value,
396     &                         valu2,i0,i0)
397                        elseif((icomplex1.ne.0).and.
398     &                         (icomplex2.ne.0)) then
399                          if(icomplex1.eq.icomplex2) then
400                            call add_sm_ei(au,ad,aub,adb,jq,
401     &                           irow,idof1,idof2,value,valu2,i0,i0)
402                            call add_sm_ei(au,ad,aub,adb,jq,
403     &                           irow,idof1+ner,idof2+ner,value,
404     &                           valu2,i0,i0)
405                          else
406                            tr=cs(15,icomplex1)*cs(15,icomplex2)
407     &                           +cs(16,icomplex1)*cs(16,icomplex2)
408                            walue=value*tr
409                            walu2=valu2*tr
410                            call add_sm_ei(au,ad,aub,adb,jq,
411     &                           irow,idof1,idof2,walue,walu2,i0,i0)
412                            call add_sm_ei(au,ad,aub,adb,jq,
413     &                           irow,idof1+ner,idof2+ner,walue,
414     &                           walu2,i0,i0)
415                            ti=cs(15,icomplex1)*cs(16,icomplex2)
416     &                           -cs(15,icomplex2)*cs(16,icomplex1)
417                            walue=value*ti
418                            walu2=valu2*ti
419                            call add_sm_ei(au,ad,aub,adb,jq,irow
420     &                           ,idof1,idof2+ner,walue,walu2,i0,i0)
421                            walue=-walue
422                            walu2=-walu2
423                            call add_sm_ei(au,ad,aub,adb,jq,irow
424     &                           ,idof1+ner,idof2,walue,walu2,i0,i0)
425                          endif
426                        elseif((icomplex1.eq.0).or.
427     &                         (icomplex2.eq.0)) then
428                          if(icomplex2.ne.0) then
429                            walue=value*cs(15,icomplex2)
430                            walu2=valu2*cs(15,icomplex2)
431                          else
432                            walue=value*cs(15,icomplex1)
433                            walu2=valu2*cs(15,icomplex1)
434                          endif
435                          call add_sm_ei(au,ad,aub,adb,jq,irow,
436     &                         idof1,idof2,walue,walu2,i0,i0)
437                          call add_sm_ei(au,ad,aub,adb,jq,irow,
438     &                         idof1+ner,idof2+ner,walue,walu2,i0,i0)
439                          if(icomplex2.ne.0) then
440                            walue=value*cs(16,icomplex2)
441                            walu2=valu2*cs(16,icomplex2)
442                          else
443                            walue=-value*cs(16,icomplex1)
444                            walu2=-valu2*cs(16,icomplex1)
445                          endif
446                          call add_sm_ei(au,ad,aub,adb,jq,irow,
447     &                         idof1,idof2+ner,walue,walu2,i0,i0)
448                          walue=-walue
449                          walu2=-walu2
450                          call add_sm_ei(au,ad,aub,adb,jq,irow,
451     &                         idof1+ner,idof2,walue,walu2,i0,i0)
452                        endif
453                      endif
454                      index2=nodempc(3,index2)
455                      if(index2.eq.0) exit
456                    enddo
457                    index1=nodempc(3,index1)
458                    if(index1.eq.0) exit
459                  enddo
460                else
461!
462!     MPC id1 / MPC id2
463!
464                  ist1=ipompc(id1)
465                  index1=nodempc(3,ist1)
466                  if(index1.eq.0) cycle
467                  do
468                    inode1=nodempc(1,index1)
469                    icomplex1=0
470                    if(labmpc(id1)(1:6).eq.'CYCLIC') then
471                      read(labmpc(id1)(7:20),'(i14)') icomplex1
472                    elseif(labmpc(id1)(1:9).eq.'SUBCYCLIC') then
473                      do ij=1,mcs
474                        ilength=int(cs(4,ij))
475                        lprev=int(cs(14,ij))
476                        call nident(ics(lprev+1),inode1,
477     &                       ilength,id)
478                        if(id.gt.0) then
479                          if(ics(lprev+id).eq.inode1) then
480                            icomplex1=ij
481                            exit
482                          endif
483                        endif
484                      enddo
485                    endif
486                    idof1=nactdof(nodempc(2,index1),inode1)
487                    ist2=ipompc(id2)
488                    index2=nodempc(3,ist2)
489                    if(index2.eq.0) then
490                      index1=nodempc(3,index1)
491                      if(index1.eq.0) then
492                        exit
493                      else
494                        cycle
495                      endif
496                    endif
497                    do
498                      inode2=nodempc(1,index2)
499                      icomplex2=0
500                      if(labmpc(id2)(1:6).eq.'CYCLIC') then
501                        read(labmpc(id2)(7:20),'(i14)') icomplex2
502                      elseif(labmpc(id2)(1:9).eq.'SUBCYCLIC') then
503                        do ij=1,mcs
504                          ilength=int(cs(4,ij))
505                          lprev=int(cs(14,ij))
506                          call nident(ics(lprev+1),inode2,
507     &                         ilength,id)
508                          if(id.gt.0) then
509                            if(ics(lprev+id).eq.inode2) then
510                              icomplex2=ij
511                              exit
512                            endif
513                          endif
514                        enddo
515                      endif
516                      idof2=nactdof(nodempc(2,index2),inode2)
517                      if((idof1.gt.0).and.(idof2.gt.0)) then
518                        value=coefmpc(index1)*coefmpc(index2)*
519     &                       s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
520                        valu2=coefmpc(index1)*coefmpc(index2)*
521     &                       sm(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
522                        if(idof1.eq.idof2) then
523                          value=2.d0*value
524                          valu2=2.d0*valu2
525                        endif
526                        if((icomplex1.eq.0).and.
527     &                       (icomplex2.eq.0)) then
528                          call add_sm_ei(au,ad,aub,adb,jq,
529     &                         irow,idof1,idof2,value,valu2,i0,i0)
530                          call add_sm_ei(au,ad,aub,adb,jq,
531     &                         irow,idof1+ner,idof2+ner,value,
532     &                         valu2,i0,i0)
533                        elseif((icomplex1.ne.0).and.
534     &                         (icomplex2.ne.0)) then
535                          if(icomplex1.eq.icomplex2) then
536                            call add_sm_ei(au,ad,aub,adb,jq,
537     &                           irow,idof1,idof2,value,valu2,i0,i0)
538                            call add_sm_ei(au,ad,aub,adb,jq,
539     &                           irow,idof1+ner,idof2+ner,value,
540     &                           valu2,i0,i0)
541                          else
542                            tr=cs(15,icomplex1)*cs(15,icomplex2)
543     &                           +cs(16,icomplex1)*cs(16,icomplex2)
544c     write(*,*) 'tr= ',tr
545                            walue=value*tr
546                            walu2=valu2*tr
547                            call add_sm_ei(au,ad,aub,adb,jq,
548     &                           irow,idof1,idof2,walue,walu2,i0,i0)
549                            call add_sm_ei(au,ad,aub,adb,jq,
550     &                           irow,idof1+ner,idof2+ner,walue,
551     &                           walu2,i0,i0)
552                            ti=cs(15,icomplex1)*cs(16,icomplex2)
553     &                           -cs(15,icomplex2)*cs(16,icomplex1)
554                            walue=value*ti
555                            walu2=valu2*ti
556                            call add_sm_ei(au,ad,aub,adb,jq,irow
557     &                           ,idof1,idof2+ner,walue,walu2,i0,i0)
558                            walue=-walue
559                            walu2=-walu2
560                            call add_sm_ei(au,ad,aub,adb,jq,irow
561     &                           ,idof1+ner,idof2,walue,walu2,i0,i0)
562                          endif
563                        elseif((icomplex1.eq.0).or.
564     &                         (icomplex2.eq.0)) then
565                          if(icomplex2.ne.0) then
566                            walue=value*cs(15,icomplex2)
567                            walu2=valu2*cs(15,icomplex2)
568                          else
569                            walue=value*cs(15,icomplex1)
570                            walu2=valu2*cs(15,icomplex1)
571                          endif
572                          call add_sm_ei(au,ad,aub,adb,jq,irow,
573     &                         idof1,idof2,walue,walu2,i0,i0)
574                          call add_sm_ei(au,ad,aub,adb,jq,irow,
575     &                         idof1+ner,idof2+ner,walue,walu2,i0,i0)
576                          if(idof1.ne.idof2) then
577                            if(icomplex2.ne.0) then
578                              walue=value*cs(16,icomplex2)
579                              walu2=valu2*cs(16,icomplex2)
580                            else
581                              walue=-value*cs(16,icomplex1)
582                              walu2=-valu2*cs(16,icomplex1)
583                            endif
584                            call add_sm_ei(au,ad,aub,adb,jq,
585     &                           irow,idof1,idof2+ner,walue,
586     &                           walu2,i0,i0)
587                            walue=-walue
588                            walu2=-walu2
589                            call add_sm_ei(au,ad,aub,adb,jq,
590     &                           irow,idof1+ner,idof2,walue,
591     &                           walu2,i0,i0)
592                          endif
593                        endif
594                      endif
595                      index2=nodempc(3,index2)
596                      if(index2.eq.0) exit
597                    enddo
598                    index1=nodempc(3,index1)
599                    if(index1.eq.0) exit
600                  enddo
601                endif
602              endif
603            endif
604          enddo
605!
606        enddo
607      enddo
608!
609      return
610      end
611