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 mafilldmss(co,nk,kon,ipkon,lakon,ne,
20     &     ipompc,nodempc,coefmpc,nmpc,
21     &     nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
22     &     ad,au,nactdof,jq,irow,neq,nmethod,
23     &     ikmpc,ilmpc,elcon,nelcon,rhcon,
24     &     nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab,
25     &     ntmat_,t0,t1,ithermal,vold,iperturb,sti,stx,iexpl,plicon,
26     &     nplicon,plkcon,nplkcon,xstiff,npmat_,dtime,
27     &     matname,mi,ncmat_,
28     &     physcon,ttime,time,istep,iinc,
29     &     ibody,xloadold,reltime,veold,springarea,nstate_,
30     &     xstateini,xstate,thicke,integerglob,doubleglob,
31     &     tieset,istartset,iendset,ialset,ntie,nasym,pslavsurf,
32     &     pmastsurf,mortar,clearini,ielprop,prop,ne0,nea,neb,
33     &     freq,ndamp,dacon,set,nset)
34!
35!     filling the stiffness matrix in spare matrix format (sm)
36!
37      implicit none
38!
39      integer mass(2),stiffness(2),buckling,rhsi,coriolis
40!
41      character*8 lakon(*)
42      character*20 sideload(*)
43      character*80 matname(*)
44      character*81 tieset(3,*),set(*)
45!
46      integer kon(*),ipompc(*),nodempc(3,*),nelemload(2,*),jq(*),
47     &     ikmpc(*),ilmpc(*),mi(*),nstate_,ne0,nasym,konl(20),
48     &     nactdof(0:mi(2),*),irow(*),icolumn,ialset(*),ielprop(*),one,
49     &     nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),ntie,
50     &     jfaces,ielorien(mi(3),*),integerglob(*),istartset(*),
51     &     iendset(*),ipkon(*),intscheme,ipobody(2,*),nbody,nset,
52     &     ibody(3,*),nk,ne,nmpc,nload,neq(2),nmethod,mscalmethod,
53     &     ithermal(*),iperturb(*),i,j,k,l,m,idist,jj,
54     &     ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2,idof1,idof2,
55     &     mpc1,mpc2,index1,index2,node1,node2,kflag,icalccg,ndamp,
56     &     ntmat_,indexe,nope,norien,iexpl,i0,ncmat_,istep,iinc,imat,
57     &     nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,mortar,kode,
58     &     nea,neb,kscale,ndof,ii,igauss,islavelinv(1),irowtloc(1),
59     &     jqtloc(1),mortartrafoflag
60!
61      real*8 co(3,*),coefmpc(*),xload(2,*),p1(3),smscale(1),
62     &     p2(3),ad(*),au(*),bodyf(3),xloadold(2,*),reltime,
63     &     t0(*),t1(*),vold(0:mi(2),*),s(60,60),
64     &     ff(60),xl(3,26),voldl(0:mi(2),26),
65     &     sti(6,mi(1),*),sm(60,60),stx(6,mi(1),*),
66     &     elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),springarea(2,*),
67     &     alcon(0:6,ntmat_,*),physcon(*),prop(*),
68     &     xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
69     &     alzero(*),orab(7,*),xbody(7,*),cgr(4,*),
70     &     plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
71     &     xstiff(27,mi(1),*),veold(0:mi(2),*),om,value,dtime,ttime,
72     &     time,thicke(mi(3),*),doubleglob(*),clearini(3,9,*),damping,
73     &     pslavsurf(3,*),pmastsurf(6,*),freq,elas(21),dacon(*),
74     &     autloc(1)
75!
76      mortartrafoflag=0
77      one=1
78!
79      kflag=2
80      i0=0
81      icalccg=0
82!
83      mass(1)=0
84      stiffness(1)=1
85      buckling=0
86      rhsi=0
87      intscheme=0
88      coriolis=0
89      kscale=1
90!
91      do i=nea,neb
92!
93        if(ipkon(i).lt.0) cycle
94!
95        if(ndamp.gt.0) then
96          damping=dacon(ielmat(1,i))
97          if(mi(3).gt.1) then
98            do j=2,mi(3)
99              if(ielmat(j,i).gt.0) then
100                if(dacon(ielmat(j,i)).ne.damping) then
101                  write(*,*) '*ERROR in mafilldmss: element',i
102                  write(*,*) '       contains several layers'
103                  write(*,*) '       with different damping '
104                  write(*,*) '       coefficients:'
105                  write(*,*) '       for layer ',one,':',damping
106                  write(*,*) '       for layer ',j,':',
107     &                 dacon(ielmat(j,i))
108                  call exit(201)
109                endif
110              endif
111            enddo
112          endif
113        else
114          damping=0.d0
115        endif
116
117        if(lakon(i)(1:2).eq.'ED') then
118          indexe=ipkon(i)
119          read(lakon(i)(8:8),'(i1)') nope
120          nope=nope+1
121!
122          do j=1,nope
123            konl(j)=kon(indexe+j)
124          enddo
125!
126          call e_damp(co,nk,konl,lakon(i),p1,p2,om,bodyf,nbody,s,
127     &         sm,ff,i,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
128     &         alzero,ielmat,ielorien,norien,orab,ntmat_,
129     &         t0,t1,ithermal,vold,iperturb,
130     &         nelemload,sideload,xload,nload,idist,sti,stx,
131     &         iexpl,plicon,nplicon,plkcon,nplkcon,xstiff,npmat_,
132     &         dtime,matname,mi(1),ncmat_,ttime,time,istep,iinc,
133     &         nmethod)
134          ndof=3
135          do jj=1,ndof*nope
136            do ii=1,jj
137              s(ii,jj)=s(ii,jj)*freq
138            enddo
139          enddo
140        elseif((lakon(i)(1:2).eq.'ES').and.
141     &         (lakon(i)(7:7).eq.'C')) then
142          indexe=ipkon(i)
143          if(mortar.eq.0) then
144            read(lakon(i)(8:8),'(i1)') nope
145            nope=nope+1
146            konl(nope+1)=kon(indexe+nope+1)
147          elseif(mortar.eq.1) then
148            nope=kon(indexe)
149          endif
150          imat=ielmat(1,i)
151!
152!     computation of the coordinates of the local nodes
153!
154          do k=1,nope
155            konl(k)=kon(indexe+k)
156            do j=1,3
157              xl(j,k)=co(j,konl(k))
158              voldl(j,k)=vold(j,konl(k))
159            enddo
160          enddo
161!
162!     initialisation of s
163!
164          do k=1,3*nope
165            do j=1,3*nope
166              s(k,j)=0.d0
167            enddo
168          enddo
169!
170          kode=nelcon(1,imat)
171!
172!     as soon as the first contact element is discovered ne0 is
173!     determined and saved
174!
175c     if(ne0.eq.0) ne0=i-1
176          if(mortar.eq.0) then
177            call springdamp_n2f(xl,elas,voldl,s,imat,elcon,
178     &           ncmat_,ntmat_,nope,iperturb,
179     &           springarea(1,konl(nope+1)),nmethod,
180     &           mi,reltime,nasym)
181          elseif(mortar.eq.1) then
182            jfaces=kon(indexe+nope+2)
183            igauss=kon(indexe+nope+1)
184            call springdamp_f2f(xl,elas,voldl,s,imat,elcon,
185     &           ncmat_,ntmat_,nope,lakon(i),iperturb,
186     &           springarea(1,igauss),
187     &           nmethod,mi,reltime,nasym,jfaces,igauss,pslavsurf,
188     &           pmastsurf,clearini)
189          endif
190          ndof=3
191          do jj=1,ndof*nope
192            do ii=1,jj
193              s(ii,jj)=s(ii,jj)*freq
194            enddo
195          enddo
196        elseif(damping.ne.0.d0) then
197          indexe=ipkon(i)
198c     Bernhardi start
199          if(lakon(i)(1:5).eq.'C3D8I') then
200            nope=11
201            ndof=3
202          elseif(lakon(i)(4:5).eq.'20') then
203c     Bernhardi end
204            nope=20
205            ndof=3
206          elseif(lakon(i)(4:4).eq.'8') then
207            nope=8
208            ndof=3
209          elseif(lakon(i)(4:5).eq.'10') then
210            nope=10
211            ndof=3
212          elseif(lakon(i)(4:4).eq.'4') then
213            nope=4
214            ndof=3
215          elseif(lakon(i)(4:5).eq.'15') then
216            nope=15
217            ndof=3
218          elseif(lakon(i)(4:4).eq.'6') then
219            nope=6
220            ndof=3
221          elseif((lakon(i)(1:2).eq.'ES').and.
222     &           (lakon(i)(7:7).ne.'F')) then
223!
224!     spring and contact spring elements (NO dashpot elements
225!     = ED... elements)
226!
227            nope=ichar(lakon(i)(8:8))-47
228            ndof=3
229          elseif(lakon(i)(1:4).eq.'MASS') then
230            nope=1
231            ndof=3
232          elseif(lakon(i)(1:1).eq.'U') then
233            ndof=ichar(lakon(i)(7:7))
234            nope=ichar(lakon(i)(8:8))
235          else
236            cycle
237          endif
238!
239          om=0.d0
240!
241          if((nbody.gt.0).and.(lakon(i)(1:1).ne.'E')) then
242!
243!     assigning centrifugal forces
244!
245            bodyf(1)=0.d0
246            bodyf(2)=0.d0
247            bodyf(3)=0.d0
248!
249            index=i
250            do
251              j=ipobody(1,index)
252              if(j.eq.0) exit
253              if(ibody(1,j).eq.1) then
254                om=xbody(1,j)
255                p1(1)=xbody(2,j)
256                p1(2)=xbody(3,j)
257                p1(3)=xbody(4,j)
258                p2(1)=xbody(5,j)
259                p2(2)=xbody(6,j)
260                p2(3)=xbody(7,j)
261!
262!     assigning gravity forces
263!
264              elseif(ibody(1,j).eq.2) then
265                bodyf(1)=bodyf(1)+xbody(1,j)*xbody(2,j)
266                bodyf(2)=bodyf(2)+xbody(1,j)*xbody(3,j)
267                bodyf(3)=bodyf(3)+xbody(1,j)*xbody(4,j)
268!
269!     assigning newton gravity forces
270!
271              elseif(ibody(1,j).eq.3) then
272                call newton(icalccg,ne,ipkon,lakon,kon,t0,co,rhcon,
273     &               nrhcon,ntmat_,physcon,i,cgr,bodyf,ielmat,
274     &               ithermal,vold,mi)
275              endif
276              index=ipobody(2,index)
277              if(index.eq.0) exit
278            enddo
279          endif
280!
281          if(lakon(i)(1:1).ne.'U') then
282            call e_c3d(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,ff,
283     &           i,nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
284     &           alzero,ielmat,ielorien,norien,orab,ntmat_,
285     &           t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload,
286     &           nload,idist,sti,stx,iexpl,plicon,
287     &           nplicon,plkcon,nplkcon,xstiff,npmat_,
288     &           dtime,matname,mi(1),ncmat_,mass(1),stiffness,buckling,
289     &           rhsi,intscheme,ttime,time,istep,iinc,coriolis,xloadold,
290     &           reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
291     &           springarea,nstate_,xstateini,xstate,ne0,ipkon,thicke,
292     &           integerglob,doubleglob,tieset,istartset,
293     &           iendset,ialset,ntie,nasym,pslavsurf,pmastsurf,mortar,
294     &           clearini,ielprop,prop,kscale,smscale(1),mscalmethod,
295     &           set,nset,islavelinv,
296     &           autloc,irowtloc,jqtloc,mortartrafoflag)
297          else
298            call e_c3d_u(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,
299     &           ff,i,nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
300     &           alzero,ielmat,ielorien,norien,orab,ntmat_,
301     &           t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload,
302     &           nload,idist,sti,stx,iexpl,plicon,
303     &           nplicon,plkcon,nplkcon,xstiff,npmat_,
304     &           dtime,matname,mi(1),ncmat_,mass(1),stiffness,buckling,
305     &           rhsi,intscheme,ttime,time,istep,iinc,coriolis,xloadold,
306     &           reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
307     &           ne0,ipkon,thicke,
308     &           integerglob,doubleglob,tieset,istartset,
309     &           iendset,ialset,ntie,nasym,
310     &           ielprop,prop)
311          endif
312          do jj=1,ndof*nope
313            do ii=1,jj
314              s(ii,jj)=s(ii,jj)*damping
315            enddo
316          enddo
317        else
318          cycle
319        endif
320!
321        if(nasym.eq.0) then
322          do jj=1,ndof*nope
323!
324            j=(jj-1)/ndof+1
325            k=jj-ndof*(j-1)
326!
327            node1=kon(indexe+j)
328            jdof1=nactdof(k,node1)
329!
330            do ll=jj,ndof*nope
331!
332              l=(ll-1)/ndof+1
333              m=ll-ndof*(l-1)
334!
335              node2=kon(indexe+l)
336              jdof2=nactdof(m,node2)
337!
338!     check whether one of the DOF belongs to a SPC or MPC
339!
340              if((jdof1.gt.0).and.(jdof2.gt.0)) then
341                call add_sm_st(au,ad,jq,irow,jdof1,jdof2,
342     &               s(jj,ll),jj,ll)
343              elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
344!
345!     idof1: genuine DOF
346!     idof2: nominal DOF of the SPC/MPC
347!
348                if(jdof1.le.0) then
349                  idof1=jdof2
350                  idof2=jdof1
351                else
352                  idof1=jdof1
353                  idof2=jdof2
354                endif
355                if(nmpc.gt.0) then
356                  if(idof2.ne.2*(idof2/2)) then
357!
358!     regular DOF / MPC
359!
360                    id=(-idof2+1)/2
361                    ist=ipompc(id)
362                    index=nodempc(3,ist)
363                    if(index.eq.0) cycle
364                    do
365                      idof2=nactdof(nodempc(2,index),
366     &                     nodempc(1,index))
367                      value=-coefmpc(index)*s(jj,ll)/
368     &                     coefmpc(ist)
369                      if(idof1.eq.idof2) value=2.d0*value
370                      if(idof2.gt.0) then
371                        call add_sm_st(au,ad,jq,irow,idof1,
372     &                       idof2,value,i0,i0)
373                      endif
374                      index=nodempc(3,index)
375                      if(index.eq.0) exit
376                    enddo
377                    cycle
378                  endif
379                endif
380!
381!     regular DOF / SPC
382!
383                if(rhsi.eq.1) then
384                elseif(nmethod.eq.2) then
385                  value=s(jj,ll)
386                  icolumn=neq(2)-idof2/2
387                  call add_bo_st(au,jq,irow,idof1,icolumn,value)
388                endif
389              else
390                idof1=jdof1
391                idof2=jdof2
392                mpc1=0
393                mpc2=0
394                if(nmpc.gt.0) then
395                  if(idof1.ne.2*(idof1/2)) mpc1=1
396                  if(idof2.ne.2*(idof2/2)) mpc2=1
397                endif
398                if((mpc1.eq.1).and.(mpc2.eq.1)) then
399                  id1=(-idof1+1)/2
400                  id2=(-idof2+1)/2
401                  if(id1.eq.id2) then
402!
403!     MPC id1 / MPC id1
404!
405                    ist=ipompc(id1)
406                    index1=nodempc(3,ist)
407                    if(index1.eq.0) cycle
408                    do
409                      idof1=nactdof(nodempc(2,index1),
410     &                     nodempc(1,index1))
411                      index2=index1
412                      do
413                        idof2=nactdof(nodempc(2,index2),
414     &                       nodempc(1,index2))
415                        value=coefmpc(index1)*coefmpc(index2)*
416     &                       s(jj,ll)/coefmpc(ist)/coefmpc(ist)
417                        if((idof1.gt.0).and.(idof2.gt.0)) then
418                          call add_sm_st(au,ad,jq,irow,
419     &                         idof1,idof2,value,i0,i0)
420                        endif
421!
422                        index2=nodempc(3,index2)
423                        if(index2.eq.0) exit
424                      enddo
425                      index1=nodempc(3,index1)
426                      if(index1.eq.0) exit
427                    enddo
428                  else
429!
430!     MPC id1 / MPC id2
431!
432                    ist1=ipompc(id1)
433                    index1=nodempc(3,ist1)
434                    if(index1.eq.0) cycle
435                    do
436                      idof1=nactdof(nodempc(2,index1),
437     &                     nodempc(1,index1))
438                      ist2=ipompc(id2)
439                      index2=nodempc(3,ist2)
440                      if(index2.eq.0) then
441                        index1=nodempc(3,index1)
442                        if(index1.eq.0) then
443                          exit
444                        else
445                          cycle
446                        endif
447                      endif
448                      do
449                        idof2=nactdof(nodempc(2,index2),
450     &                       nodempc(1,index2))
451                        value=coefmpc(index1)*coefmpc(index2)*
452     &                       s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
453                        if(idof1.eq.idof2) value=2.d0*value
454                        if((idof1.gt.0).and.(idof2.gt.0)) then
455                          call add_sm_st(au,ad,jq,irow,
456     &                         idof1,idof2,value,i0,i0)
457                        endif
458!
459                        index2=nodempc(3,index2)
460                        if(index2.eq.0) exit
461                      enddo
462                      index1=nodempc(3,index1)
463                      if(index1.eq.0) exit
464                    enddo
465                  endif
466                endif
467              endif
468            enddo
469!
470!
471          enddo
472        endif
473      enddo
474!
475      return
476      end
477
478