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 createmddof(imddof,nmddof,istartset,iendset,
20     &            ialset,nactdof,ithermal,mi,imdnode,nmdnode,ikmpc,
21     &            ilmpc,ipompc,nodempc,nmpc,
22     &            imdmpc,nmdmpc,imdboun,nmdboun,ikboun,nboun,
23     &            nset,ntie,tieset,set,lakon,kon,ipkon,labmpc,
24     &            ilboun,filab,prlab,prset,nprint,ne,cyclicsymmetry)
25!
26!     creating a set imddof containing the degrees of freedom
27!     selected by the user for modal dynamic calculations. The
28!     solution will be calculated for these dof's only in order
29!     to speed up the calculation.
30!
31      implicit none
32!
33      logical nodeslavsurf
34!
35      character*6 prlab(*)
36      character*8 lakon(*)
37      character*20 labmpc(*)
38      character*81 tieset(3,*),rightset,set(*),slavset,noset,prset(*)
39      character*87 filab(*)
40!
41      integer imddof(*),nmddof,nrset,istartset(*),iendset(*),mi(*),
42     &  ialset(*),nactdof(0:mi(2),*),node,ithermal(*),j,k,l,
43     &  ikmpc(*),ilmpc(*),ipompc(*),nodempc(3,*),nmpc,id,
44     &  imdnode(*),nmdnode,imdmpc(*),nmdmpc,nprint,ipos,
45     &  imdboun(*),nmdboun,ikboun(*),nboun,indexe1,indexe,islav,
46     &  jface,nset,ntie,nnodelem,nope,nodef(8),nelem,nface,imast,
47     &  ifaceq(8,6),ifacet(6,4),ifacew1(4,5),ifacew2(8,5),kon(*),
48     &  ipkon(*),i,ilboun(*),nlabel,ne,cyclicsymmetry
49!
50!     nodes per face for hex elements
51!
52      data ifaceq /4,3,2,1,11,10,9,12,
53     &            5,6,7,8,13,14,15,16,
54     &            1,2,6,5,9,18,13,17,
55     &            2,3,7,6,10,19,14,18,
56     &            3,4,8,7,11,20,15,19,
57     &            4,1,5,8,12,17,16,20/
58!
59!     nodes per face for tet elements
60!
61      data ifacet /1,3,2,7,6,5,
62     &             1,2,4,5,9,8,
63     &             2,3,4,6,10,9,
64     &             1,4,3,8,10,7/
65!
66!     nodes per face for linear wedge elements
67!
68      data ifacew1 /1,3,2,0,
69     &             4,5,6,0,
70     &             1,2,5,4,
71     &             2,3,6,5,
72     &             3,1,4,6/
73!
74!     nodes per face for quadratic wedge elements
75!
76      data ifacew2 /1,3,2,9,8,7,0,0,
77     &             4,5,6,10,11,12,0,0,
78     &             1,2,5,4,7,14,10,13,
79     &             2,3,6,5,8,15,11,14,
80     &             3,1,4,6,9,13,12,15/
81!
82      data nlabel /55/
83!
84!     if 1d/2d elements are part of the mesh, no node selection
85!     is performed (because of the renumbering due to the
86!     expansion node selection is excessively difficult)
87!
88      do i=1,ne
89         if((lakon(i)(7:7).eq.'E').or.
90     &      (lakon(i)(7:7).eq.'S').or.
91     &      ((lakon(i)(7:7).eq.'A').and.(lakon(i)(1:1).eq.'C')).or.
92     &      (lakon(i)(7:7).eq.'L').or.
93     &      (lakon(i)(7:7).eq.'B')) then
94            nmdnode=0
95            nmddof=0
96            nmdboun=0
97            nmdmpc=0
98            return
99         endif
100      enddo
101!
102!     storing the nodes, dofs, spcs and mpcs for which *NODE FILE
103!     or *EL FILE was selected
104!
105      do i=1,nlabel
106!
107!        CDIS,CSTR und CELS are not taken into account:
108!        contact area is treated separately (no set can
109!        be specified for CDIS, CSTR und CELS)
110!
111         if((i.eq.26).or.(i.eq.27)) cycle
112!
113         if(filab(i)(1:1).ne.' ') then
114            read(filab(i)(7:87),'(a81)') noset
115c            nrset=0
116c            do k=1,nset
117c               if(set(k).eq.noset) then
118c                  nrset=k
119c                  exit
120c               endif
121c            enddo
122            call cident81(set,noset,nset,id)
123            nrset=0
124            if(id.gt.0) then
125              if(noset.eq.set(id)) then
126                nrset=id
127              endif
128            endif
129!
130!           if output for all nodes is selected, use
131!           of imdnode is deactivated
132!
133            if(nrset.eq.0) then
134               if(cyclicsymmetry.eq.1) then
135                  write(*,*) '*ERROR in createmddof: in a cylic'
136                  write(*,*) '       symmetric modal dynamic or'
137                  write(*,*) '       steady static dynamics calculation'
138                  write(*,*) '       a node set MUST be defined on each'
139                  write(*,*) '       *NODE FILE, *NODE OUTPUT, *EL FILE'
140                  write(*,*) '       or *ELEMENT OUTPUT card.'
141                  write(*,*) '       Justification: in a steady state'
142                  write(*,*) '       dynamics calculation with cyclic'
143                  write(*,*) '       symmetry the segment is expanded'
144                  write(*,*) '       into 360 °. Storing results for'
145                  write(*,*) '       this expansion may lead to huge'
146                  write(*,*) '       frd-files. Specifying a set can'
147                  write(*,*) '       reduce this output.'
148                  call exit(201)
149               endif
150               nmdnode=0
151               nmddof=0
152               nmdboun=0
153               nmdmpc=0
154               return
155            endif
156!
157!           adding the nodes belonging to nrset
158!
159            do j=istartset(nrset),iendset(nrset)
160               if(ialset(j).gt.0) then
161                  node=ialset(j)
162                  call addimd(imdnode,nmdnode,node)
163                  if(ithermal(1).ne.2) then
164                     do k=1,3
165                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
166     &                   nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
167     &                   nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
168     &                   ikboun,nboun,ilboun)
169                     enddo
170                  else
171                     k=0
172                     call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
173     &                nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
174     &                nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,ikboun,
175     &                nboun,ilboun)
176                  endif
177               else
178                  node=ialset(j-2)
179                  do
180                     node=node-ialset(j)
181                     if(node.ge.ialset(j-1)) exit
182                     call addimd(imdnode,nmdnode,node)
183                     if(ithermal(1).ne.2) then
184                        do k=1,3
185                           call addimdnodedof(node,k,ikmpc,ilmpc,
186     &                      ipompc,nodempc,nmpc,imdnode,nmdnode,imddof,
187     &                      nmddof,nactdof,mi,imdmpc,nmdmpc,imdboun,
188     &                      nmdboun,ikboun,nboun,ilboun)
189                        enddo
190                     else
191                        k=0
192                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
193     &                   nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
194     &                   nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
195     &                   ikboun,nboun,ilboun)
196                     endif
197                  enddo
198               endif
199            enddo
200!
201         endif
202      enddo
203!
204!     storing the nodes, dofs, spcs and mpcs for which *NODE PRINT
205!     was selected
206!
207      do i=1,nprint
208         if((prlab(i)(1:4).eq.'U   ').or.
209     &        (prlab(i)(1:4).eq.'NT  ').or.
210     &        (prlab(i)(1:4).eq.'RF  ').or.
211     &        (prlab(i)(1:4).eq.'RFL ').or.
212     &        (prlab(i)(1:4).eq.'PS  ').or.
213     &        (prlab(i)(1:4).eq.'PN  ').or.
214     &        (prlab(i)(1:4).eq.'MF  ').or.
215     &        (prlab(i)(1:4).eq.'V   ')) then
216            noset=prset(i)
217c            nrset=0
218c            do k=1,nset
219c               if(set(k).eq.noset) then
220c                  nrset=k
221c                  exit
222c               endif
223c            enddo
224            call cident81(set,noset,nset,id)
225            nrset=0
226            if(id.gt.0) then
227              if(noset.eq.set(id)) then
228                nrset=id
229              endif
230            endif
231!
232!           adding the nodes belonging to nrset
233!
234            do j=istartset(nrset),iendset(nrset)
235               if(ialset(j).gt.0) then
236                  node=ialset(j)
237                  call addimd(imdnode,nmdnode,node)
238                  if(ithermal(1).ne.2) then
239                     do k=1,3
240                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
241     &                   nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
242     &                   nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
243     &                   ikboun,nboun,ilboun)
244                     enddo
245                  else
246                     k=0
247                     call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
248     &                nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
249     &                nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,ikboun,
250     &                nboun,ilboun)
251                  endif
252               else
253                  node=ialset(j-2)
254                  do
255                     node=node-ialset(j)
256                     if(node.ge.ialset(j-1)) exit
257                     call addimd(imdnode,nmdnode,node)
258                     if(ithermal(1).ne.2) then
259                        do k=1,3
260                           call addimdnodedof(node,k,ikmpc,ilmpc,
261     &                      ipompc,nodempc,nmpc,imdnode,nmdnode,imddof,
262     &                      nmddof,nactdof,mi,imdmpc,nmdmpc,imdboun,
263     &                      nmdboun,ikboun,nboun,ilboun)
264                        enddo
265                     else
266                        k=0
267                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
268     &                   nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
269     &                   nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
270     &                   ikboun,nboun,ilboun)
271                     endif
272                  enddo
273               endif
274            enddo
275         endif
276      enddo
277!
278!     check whether all contact slave and master nodes (and corresponding
279!     dofs, spcs and mpcs were selected
280!
281      do i=1,ntie
282!
283!     check for contact conditions
284!     'C' are active contact conditions
285!     '-' are temporarily deactivated contact conditions
286!
287         if((tieset(1,i)(81:81).eq.'C').or.
288     &      (tieset(1,i)(81:81).eq.'-')) then
289            rightset=tieset(3,i)
290!
291!     determining the master surface
292!
293c            do j=1,nset
294c               if(set(j).eq.rightset) exit
295c            enddo
296            call cident81(set,rightset,nset,id)
297            j=nset+1
298            if(id.gt.0) then
299              if(rightset.eq.set(id)) then
300                j=id
301              endif
302            endif
303            if(j.gt.nset) then
304               write(*,*) '*ERROR in createmddof: master surface',
305     &              rightset
306               write(*,*) '       does not exist'
307               call exit(201)
308            endif
309            imast=j
310!
311            do j=istartset(imast),iendset(imast)
312!
313               nelem=int(ialset(j)/10.d0)
314               jface=ialset(j)-10*nelem
315!
316               indexe=ipkon(nelem)
317!
318               if(lakon(nelem)(4:4).eq.'2') then
319                  nnodelem=8
320                  nface=6
321               elseif(lakon(nelem)(4:4).eq.'8') then
322                  nnodelem=4
323                  nface=6
324               elseif(lakon(nelem)(4:5).eq.'10') then
325                  nnodelem=6
326                  nface=4
327               elseif(lakon(nelem)(4:4).eq.'4') then
328                  nnodelem=3
329                  nface=4
330               elseif(lakon(nelem)(4:5).eq.'15') then
331                  if(jface.le.2) then
332                     nnodelem=6
333                  else
334                     nnodelem=8
335                  endif
336                  nface=5
337                  nope=15
338               elseif(lakon(nelem)(4:4).eq.'6') then
339                  if(jface.le.2) then
340                     nnodelem=3
341                  else
342                     nnodelem=4
343                  endif
344                  nface=5
345                  nope=6
346               else
347                  cycle
348               endif
349!
350!     determining the master nodes
351!
352               if(nface.eq.4) then
353                  do k=1,nnodelem
354                     nodef(k)=kon(indexe+ifacet(k,jface))
355                  enddo
356               elseif(nface.eq.5) then
357                  if(nope.eq.6) then
358                     do k=1,nnodelem
359                        nodef(k)=kon(indexe+ifacew1(k,jface))
360                     enddo
361                  elseif(nope.eq.15) then
362                     do k=1,nnodelem
363                        nodef(k)=kon(indexe+ifacew2(k,jface))
364                     enddo
365                  endif
366               elseif(nface.eq.6) then
367                  do k=1,nnodelem
368                     nodef(k)=kon(indexe+ifaceq(k,jface))
369                  enddo
370               endif
371!
372               do l=1,nnodelem
373                  node=nodef(l)
374                  call addimd(imdnode,nmdnode,node)
375                  if(ithermal(1).ne.2) then
376                     do k=1,3
377                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
378     &                       nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
379     &                       nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
380     &                       ikboun,nboun,ilboun)
381                     enddo
382                  endif
383               enddo
384            enddo
385!
386!           determining the slave nodes
387!
388            slavset=tieset(2,i)
389!
390!           check whether facial slave surface;
391!
392            ipos=index(slavset,' ')-1
393!
394!           default for node-to-surface contact is
395!           a nodal slave surface
396!
397            if(slavset(ipos:ipos).eq.'S') then
398!
399!              'S' means node-to-face contact, it does not mean
400!              that the slave surface is node-based;
401!              start with assuming a nodal surface (default),
402!              if not present, check for a facial surface
403!
404               nodeslavsurf=.true.
405            else
406               nodeslavsurf=.false.
407            endif
408!
409!           determining the slave surface
410!
411c            do j=1,nset
412c               if(set(j).eq.slavset) exit
413c            enddo
414            call cident81(set,slavset,nset,id)
415            j=nset+1
416            if(id.gt.0) then
417              if(slavset.eq.set(id)) then
418                j=id
419              endif
420            endif
421            if(j.gt.nset) then
422               do j=1,nset
423                  if((set(j)(1:ipos-1).eq.slavset(1:ipos-1)).and.
424     &                 (set(j)(ipos:ipos).eq.'T')) then
425                     nodeslavsurf=.false.
426                     exit
427                  endif
428               enddo
429            endif
430!
431            islav=j
432!
433            if(nodeslavsurf) then
434!
435!              nodal slave surface
436!
437               do j=istartset(islav),iendset(islav)
438                  if(ialset(j).gt.0) then
439                     node=ialset(j)
440                     call addimd(imdnode,nmdnode,node)
441                     if(ithermal(1).ne.2) then
442                        do k=1,3
443                           call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
444     &                       nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
445     &                       nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
446     &                       ikboun,nboun,ilboun)
447                        enddo
448                     endif
449                  else
450                     k=ialset(j-2)
451                     do
452                        k=k-ialset(j)
453                        if(k.ge.ialset(j-1)) exit
454                        node=k
455                        call addimd(imdnode,nmdnode,node)
456                        if(ithermal(1).ne.2) then
457                           do k=1,3
458                              call addimdnodedof(node,k,ikmpc,ilmpc,
459     &                          ipompc,nodempc,nmpc,imdnode,nmdnode,
460     &                          imddof,nmddof,nactdof,mi,imdmpc,nmdmpc,
461     &                          imdboun,nmdboun,ikboun,nboun,ilboun)
462                           enddo
463                        endif
464                     enddo
465                  endif
466               enddo
467            else
468!
469!             facial slave surface
470!
471               do j=istartset(islav),iendset(islav)
472!
473                  nelem=int(ialset(j)/10.d0)
474                  jface=ialset(j)-10*nelem
475!
476                  indexe=ipkon(nelem)
477!
478                  if(lakon(nelem)(4:4).eq.'2') then
479                     nnodelem=8
480                     nface=6
481                  elseif(lakon(nelem)(4:4).eq.'8') then
482                     nnodelem=4
483                     nface=6
484                  elseif(lakon(nelem)(4:5).eq.'10') then
485                     nnodelem=6
486                     nface=4
487                  elseif(lakon(nelem)(4:4).eq.'4') then
488                     nnodelem=3
489                     nface=4
490                  elseif(lakon(nelem)(4:5).eq.'15') then
491                     if(jface.le.2) then
492                        nnodelem=6
493                     else
494                        nnodelem=8
495                     endif
496                     nface=5
497                     nope=15
498                  elseif(lakon(nelem)(4:4).eq.'6') then
499                     if(jface.le.2) then
500                        nnodelem=3
501                     else
502                        nnodelem=4
503                     endif
504                     nface=5
505                     nope=6
506                  else
507                     cycle
508                  endif
509!
510!     determining the slave nodes
511!
512                  if(nface.eq.4) then
513                     do k=1,nnodelem
514                        nodef(k)=kon(indexe+ifacet(k,jface))
515                     enddo
516                  elseif(nface.eq.5) then
517                     if(nope.eq.6) then
518                        do k=1,nnodelem
519                           nodef(k)=kon(indexe+ifacew1(k,jface))
520                        enddo
521                     elseif(nope.eq.15) then
522                        do k=1,nnodelem
523                           nodef(k)=kon(indexe+ifacew2(k,jface))
524                        enddo
525                     endif
526                  elseif(nface.eq.6) then
527                     do k=1,nnodelem
528                        nodef(k)=kon(indexe+ifaceq(k,jface))
529                     enddo
530                  endif
531!
532                  do l=1,nnodelem
533                     node=nodef(l)
534                     call addimd(imdnode,nmdnode,node)
535                     if(ithermal(1).ne.2) then
536                        do k=1,3
537                           call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
538     &                       nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
539     &                       nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
540     &                       ikboun,nboun,ilboun)
541                        enddo
542                     endif
543                  enddo
544               enddo
545            endif
546!
547         endif
548      enddo
549!
550!     adding nodes, dofs, spcs and mpcs belonging to nonlinear MPC's
551!     (why only dependent nodes?)
552!
553      do i=1,nmpc
554         if((labmpc(i)(1:20).ne.'                    ').and.
555     &          (labmpc(i)(1:6).ne.'CYCLIC').and.
556     &          (labmpc(i)(1:9).ne.'SUBCYCLIC')) then
557            indexe1=ipompc(i)
558            if(indexe1.eq.0) cycle
559            node=nodempc(1,indexe1)
560            call addimd(imdnode,nmdnode,node)
561            if(ithermal(1).ne.2) then
562               do k=1,3
563                  call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
564     &                 nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
565     &                 nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
566     &                 ikboun,nboun,ilboun)
567               enddo
568            else
569               k=0
570               call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
571     &              nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
572     &              nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,ikboun,
573     &              nboun,ilboun)
574            endif
575         endif
576      enddo
577!
578      return
579      end
580
581
582
583
584