← Back to branch summary

~gabriel1984sibiu/calculix/ccx

« back to all changes in this revision

Viewing changes to CalculiX/ccx_2.11/src/createmddof.f

  • Committer: Grevutiu Gabriel
  • Date: 2016-12-30 12:06:41 UTC
  • Revision ID: gabriel1984sibiu@gmail.com-20161230120641-kzmhfy8mn00w3mhg
New upstream version

Show diffs side-by-side

added added

removed removed

Lines of Context:
CalculiX/ccx_2.11/src/createmddof.f
 
1
!
 
2
!     CalculiX - A 3-dimensional finite element program
 
3
!              Copyright (C) 1998-2015 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,
 
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 /47/
 
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 for which *NODE FILE or *EL FILE was selected
 
103
!
 
104
      do i=1,nlabel
 
105
!
 
106
!        CDIS,CSTR und CELS are not taken into account:
 
107
!        contact area is treated separately (no set can
 
108
!        be specified for CDIS, CSTR und CELS)
 
109
!
 
110
         if((i.eq.26).or.(i.eq.27)) cycle
 
111
!
 
112
         if(filab(i)(1:1).ne.' ') then
 
113
            read(filab(i)(7:87),'(a81)') noset
 
114
            nrset=0
 
115
            do k=1,nset
 
116
               if(set(k).eq.noset) then
 
117
                  nrset=k
 
118
                  exit
 
119
               endif
 
120
            enddo
 
121
!
 
122
!           if output for all nodes is selected, use
 
123
!           of imdnode is deactivated
 
124
!
 
125
            if(nrset.eq.0) then
 
126
               if(cyclicsymmetry.eq.1) then
 
127
                  write(*,*) '*ERROR in createmddof: in a cylic'
 
128
                  write(*,*) '       symmetric modal dynamic or'
 
129
                  write(*,*) '       steady static dynamics calculation'
 
130
                  write(*,*) '       a node set MUST be defined on each'
 
131
                  write(*,*) '       *NODE FILE, *NODE OUTPUT, *EL FILE'
 
132
                  write(*,*) '       or *ELEMENT OUTPUT card.'
 
133
                  write(*,*) '       Justification: in a steady state'
 
134
                  write(*,*) '       dynamics calculation with cyclic'
 
135
                  write(*,*) '       symmetry the segment is expanded'
 
136
                  write(*,*) '       into 360 °. Storing results for'
 
137
                  write(*,*) '       this expansion may lead to huge'
 
138
                  write(*,*) '       frd-files. Specifying a set can'
 
139
                  write(*,*) '       reduce this output.'
 
140
                  call exit(201)
 
141
               endif
 
142
               nmdnode=0
 
143
               nmddof=0
 
144
               nmdboun=0
 
145
               nmdmpc=0
 
146
               return
 
147
            endif
 
148
!
 
149
!           adding the nodes belonging to nrset
 
150
!
 
151
            do j=istartset(nrset),iendset(nrset)
 
152
               if(ialset(j).gt.0) then
 
153
                  node=ialset(j)
 
154
                  call addimd(imdnode,nmdnode,node)
 
155
                  if(ithermal.ne.2) then
 
156
                     do k=1,3
 
157
                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
158
     &                   nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
159
     &                   nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
 
160
     &                   ikboun,nboun,ilboun)
 
161
                     enddo
 
162
                  else
 
163
                     k=0
 
164
                     call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
165
     &                nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
166
     &                nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,ikboun,
 
167
     &                nboun,ilboun)
 
168
                  endif
 
169
               else
 
170
                  node=ialset(j-2)
 
171
                  do
 
172
                     node=node-ialset(j)
 
173
                     if(node.ge.ialset(j-1)) exit
 
174
                     call addimd(imdnode,nmdnode,node)
 
175
                     if(ithermal.ne.2) then
 
176
                        do k=1,3
 
177
                           call addimdnodedof(node,k,ikmpc,ilmpc,
 
178
     &                      ipompc,nodempc,nmpc,imdnode,nmdnode,imddof,
 
179
     &                      nmddof,nactdof,mi,imdmpc,nmdmpc,imdboun,
 
180
     &                      nmdboun,ikboun,nboun,ilboun)
 
181
                        enddo
 
182
                     else
 
183
                        k=0
 
184
                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
185
     &                   nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
186
     &                   nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
 
187
     &                   ikboun,nboun,ilboun)
 
188
                     endif
 
189
                  enddo
 
190
               endif
 
191
            enddo
 
192
!
 
193
         endif
 
194
      enddo
 
195
!
 
196
!     storing the nodes for which *NODE PRINT was selected
 
197
!
 
198
      do i=1,nprint
 
199
         if((prlab(i)(1:4).eq.'U   ').or.
 
200
     &        (prlab(i)(1:4).eq.'NT  ').or.
 
201
     &        (prlab(i)(1:4).eq.'RF  ').or.
 
202
     &        (prlab(i)(1:4).eq.'RFL ').or.
 
203
     &        (prlab(i)(1:4).eq.'PS  ').or.
 
204
     &        (prlab(i)(1:4).eq.'PN  ').or.
 
205
     &        (prlab(i)(1:4).eq.'MF  ').or.
 
206
     &        (prlab(i)(1:4).eq.'V   ')) then
 
207
            noset=prset(i)
 
208
            nrset=0
 
209
            do k=1,nset
 
210
               if(set(k).eq.noset) then
 
211
                  nrset=k
 
212
                  exit
 
213
               endif
 
214
            enddo
 
215
!
 
216
!           adding the nodes belonging to nrset
 
217
!
 
218
            do j=istartset(nrset),iendset(nrset)
 
219
               if(ialset(j).gt.0) then
 
220
                  node=ialset(j)
 
221
                  call addimd(imdnode,nmdnode,node)
 
222
                  if(ithermal.ne.2) then
 
223
                     do k=1,3
 
224
                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
225
     &                   nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
226
     &                   nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
 
227
     &                   ikboun,nboun,ilboun)
 
228
                     enddo
 
229
                  else
 
230
                     k=0
 
231
                     call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
232
     &                nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
233
     &                nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,ikboun,
 
234
     &                nboun,ilboun)
 
235
                  endif
 
236
               else
 
237
                  node=ialset(j-2)
 
238
                  do
 
239
                     node=node-ialset(j)
 
240
                     if(node.ge.ialset(j-1)) exit
 
241
                     call addimd(imdnode,nmdnode,node)
 
242
                     if(ithermal.ne.2) then
 
243
                        do k=1,3
 
244
                           call addimdnodedof(node,k,ikmpc,ilmpc,
 
245
     &                      ipompc,nodempc,nmpc,imdnode,nmdnode,imddof,
 
246
     &                      nmddof,nactdof,mi,imdmpc,nmdmpc,imdboun,
 
247
     &                      nmdboun,ikboun,nboun,ilboun)
 
248
                        enddo
 
249
                     else
 
250
                        k=0
 
251
                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
252
     &                   nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
253
     &                   nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
 
254
     &                   ikboun,nboun,ilboun)
 
255
                     endif
 
256
                  enddo
 
257
               endif
 
258
            enddo
 
259
         endif
 
260
      enddo
 
261
!
 
262
!     check whether all contact slave and master nodes were selected
 
263
!
 
264
      do i=1,ntie
 
265
!     
 
266
!     check for contact conditions
 
267
!     'C' are active contact conditions
 
268
!     '-' are temporarily deactivated contact conditions
 
269
!     
 
270
         if((tieset(1,i)(81:81).eq.'C').or.
 
271
     &      (tieset(1,i)(81:81).eq.'-')) then
 
272
            rightset=tieset(3,i)
 
273
!     
 
274
!     determining the master surface
 
275
!     
 
276
            do j=1,nset
 
277
               if(set(j).eq.rightset) exit
 
278
            enddo
 
279
            if(j.gt.nset) then
 
280
               write(*,*) '*ERROR in createmddof: master surface',
 
281
     &              rightset
 
282
               write(*,*) '       does not exist'
 
283
               call exit(201)
 
284
            endif
 
285
            imast=j
 
286
!     
 
287
            do j=istartset(imast),iendset(imast)
 
288
!     
 
289
               nelem=int(ialset(j)/10.d0)
 
290
               jface=ialset(j)-10*nelem
 
291
!     
 
292
               indexe=ipkon(nelem)
 
293
!     
 
294
               if(lakon(nelem)(4:4).eq.'2') then
 
295
                  nnodelem=8
 
296
                  nface=6
 
297
               elseif(lakon(nelem)(4:4).eq.'8') then
 
298
                  nnodelem=4
 
299
                  nface=6
 
300
               elseif(lakon(nelem)(4:5).eq.'10') then
 
301
                  nnodelem=6
 
302
                  nface=4
 
303
               elseif(lakon(nelem)(4:4).eq.'4') then
 
304
                  nnodelem=3
 
305
                  nface=4
 
306
               elseif(lakon(nelem)(4:5).eq.'15') then
 
307
                  if(jface.le.2) then
 
308
                     nnodelem=6
 
309
                  else
 
310
                     nnodelem=8
 
311
                  endif
 
312
                  nface=5
 
313
                  nope=15
 
314
               elseif(lakon(nelem)(4:4).eq.'6') then
 
315
                  if(jface.le.2) then
 
316
                     nnodelem=3
 
317
                  else
 
318
                     nnodelem=4
 
319
                  endif
 
320
                  nface=5
 
321
                  nope=6
 
322
               else
 
323
                  cycle
 
324
               endif
 
325
!     
 
326
!     determining the master nodes 
 
327
!     
 
328
               if(nface.eq.4) then
 
329
                  do k=1,nnodelem
 
330
                     nodef(k)=kon(indexe+ifacet(k,jface))
 
331
                  enddo
 
332
               elseif(nface.eq.5) then
 
333
                  if(nope.eq.6) then
 
334
                     do k=1,nnodelem
 
335
                        nodef(k)=kon(indexe+ifacew1(k,jface))
 
336
                     enddo
 
337
                  elseif(nope.eq.15) then
 
338
                     do k=1,nnodelem
 
339
                        nodef(k)=kon(indexe+ifacew2(k,jface))
 
340
                     enddo
 
341
                  endif
 
342
               elseif(nface.eq.6) then
 
343
                  do k=1,nnodelem
 
344
                     nodef(k)=kon(indexe+ifaceq(k,jface))
 
345
                  enddo
 
346
               endif
 
347
!
 
348
               do l=1,nnodelem
 
349
                  node=nodef(l)
 
350
                  call addimd(imdnode,nmdnode,node)
 
351
                  if(ithermal.ne.2) then
 
352
                     do k=1,3
 
353
                        call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
354
     &                       nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
355
     &                       nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
 
356
     &                       ikboun,nboun,ilboun)
 
357
                     enddo
 
358
                  endif
 
359
               enddo
 
360
            enddo
 
361
!
 
362
!           determining the slave nodes 
 
363
!
 
364
            slavset=tieset(2,i)
 
365
!
 
366
!           check whether facial slave surface; 
 
367
!
 
368
            ipos=index(slavset,' ')-1
 
369
! 
 
370
!           default for node-to-surface contact is
 
371
!           a nodal slave surface
 
372
!
 
373
            if(slavset(ipos:ipos).eq.'S') then
 
374
               nodeslavsurf=.true.
 
375
            endif
 
376
!
 
377
!           determining the slave surface 
 
378
!
 
379
            do j=1,nset
 
380
               if(set(j).eq.slavset) exit
 
381
            enddo
 
382
            if(j.gt.nset) then
 
383
               do j=1,nset
 
384
                  if((set(j)(1:ipos-1).eq.slavset(1:ipos-1)).and.
 
385
     &                 (set(j)(ipos:ipos).eq.'T')) then
 
386
                     nodeslavsurf=.false.
 
387
                     exit
 
388
                  endif
 
389
               enddo
 
390
            endif
 
391
!
 
392
            islav=j
 
393
!
 
394
            if(nodeslavsurf) then
 
395
!
 
396
!              nodal slave surface
 
397
!
 
398
               do j=istartset(islav),iendset(islav)
 
399
                  if(ialset(j).gt.0) then
 
400
                     node=ialset(j)
 
401
                     call addimd(imdnode,nmdnode,node)
 
402
                     if(ithermal.ne.2) then
 
403
                        do k=1,3
 
404
                           call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
405
     &                       nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
406
     &                       nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
 
407
     &                       ikboun,nboun,ilboun)
 
408
                        enddo
 
409
                     endif
 
410
                  else
 
411
                     k=ialset(j-2)
 
412
                     do
 
413
                        k=k-ialset(j)
 
414
                        if(k.ge.ialset(j-1)) exit
 
415
                        node=k
 
416
                        call addimd(imdnode,nmdnode,node)
 
417
                        if(ithermal.ne.2) then
 
418
                           do k=1,3
 
419
                              call addimdnodedof(node,k,ikmpc,ilmpc,
 
420
     &                          ipompc,nodempc,nmpc,imdnode,nmdnode,
 
421
     &                          imddof,nmddof,nactdof,mi,imdmpc,nmdmpc,
 
422
     &                          imdboun,nmdboun,ikboun,nboun,ilboun)
 
423
                           enddo
 
424
                        endif
 
425
                     enddo
 
426
                  endif
 
427
               enddo
 
428
            else
 
429
!     
 
430
!             facial slave surface
 
431
!
 
432
               do j=istartset(islav),iendset(islav)
 
433
!     
 
434
                  nelem=int(ialset(j)/10.d0)
 
435
                  jface=ialset(j)-10*nelem
 
436
!     
 
437
                  indexe=ipkon(nelem)
 
438
!     
 
439
                  if(lakon(nelem)(4:4).eq.'2') then
 
440
                     nnodelem=8
 
441
                     nface=6
 
442
                  elseif(lakon(nelem)(4:4).eq.'8') then
 
443
                     nnodelem=4
 
444
                     nface=6
 
445
                  elseif(lakon(nelem)(4:5).eq.'10') then
 
446
                     nnodelem=6
 
447
                     nface=4
 
448
                  elseif(lakon(nelem)(4:4).eq.'4') then
 
449
                     nnodelem=3
 
450
                     nface=4
 
451
                  elseif(lakon(nelem)(4:5).eq.'15') then
 
452
                     if(jface.le.2) then
 
453
                        nnodelem=6
 
454
                     else
 
455
                        nnodelem=8
 
456
                     endif
 
457
                     nface=5
 
458
                     nope=15
 
459
                  elseif(lakon(nelem)(4:4).eq.'6') then
 
460
                     if(jface.le.2) then
 
461
                        nnodelem=3
 
462
                     else
 
463
                        nnodelem=4
 
464
                     endif
 
465
                     nface=5
 
466
                     nope=6
 
467
                  else
 
468
                     cycle
 
469
                  endif
 
470
!     
 
471
!     determining the slave nodes 
 
472
!     
 
473
                  if(nface.eq.4) then
 
474
                     do k=1,nnodelem
 
475
                        nodef(k)=kon(indexe+ifacet(k,jface))
 
476
                     enddo
 
477
                  elseif(nface.eq.5) then
 
478
                     if(nope.eq.6) then
 
479
                        do k=1,nnodelem
 
480
                           nodef(k)=kon(indexe+ifacew1(k,jface))
 
481
                        enddo
 
482
                     elseif(nope.eq.15) then
 
483
                        do k=1,nnodelem
 
484
                           nodef(k)=kon(indexe+ifacew2(k,jface))
 
485
                        enddo
 
486
                     endif
 
487
                  elseif(nface.eq.6) then
 
488
                     do k=1,nnodelem
 
489
                        nodef(k)=kon(indexe+ifaceq(k,jface))
 
490
                     enddo
 
491
                  endif
 
492
!     
 
493
                  do l=1,nnodelem
 
494
                     node=nodef(l)
 
495
                     call addimd(imdnode,nmdnode,node)
 
496
                     if(ithermal.ne.2) then
 
497
                        do k=1,3
 
498
                           call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
499
     &                       nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
500
     &                       nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
 
501
     &                       ikboun,nboun,ilboun)
 
502
                        enddo
 
503
                     endif
 
504
                  enddo
 
505
               enddo
 
506
            endif
 
507
!     
 
508
         endif
 
509
      enddo
 
510
!
 
511
!     adding nodes belonging to nonlinear MPC's (why only dependent nodes?)
 
512
!      
 
513
      do i=1,nmpc
 
514
         if((labmpc(i)(1:20).ne.'                    ').and.
 
515
     &          (labmpc(i)(1:7).ne.'CONTACT').and.
 
516
     &          (labmpc(i)(1:6).ne.'CYCLIC').and.
 
517
     &          (labmpc(i)(1:9).ne.'SUBCYCLIC')) then
 
518
            indexe1=ipompc(i)
 
519
            if(indexe1.eq.0) cycle
 
520
            node=nodempc(1,indexe1)
 
521
            call addimd(imdnode,nmdnode,node)
 
522
            if(ithermal.ne.2) then
 
523
               do k=1,3
 
524
                  call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
525
     &                 nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
526
     &                 nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,
 
527
     &                 ikboun,nboun,ilboun)
 
528
               enddo
 
529
            else
 
530
               k=0
 
531
               call addimdnodedof(node,k,ikmpc,ilmpc,ipompc,
 
532
     &              nodempc,nmpc,imdnode,nmdnode,imddof,nmddof,
 
533
     &              nactdof,mi,imdmpc,nmdmpc,imdboun,nmdboun,ikboun,
 
534
     &              nboun,ilboun)
 
535
            endif
 
536
         endif
 
537
      enddo
 
538
!
 
539
!     subtracting 1 to comply with the C-convention
 
540
!
 
541
c      do j=1,nmddof
 
542
c         imddof(j)=imddof(j)-1
 
543
c      enddo
 
544
!
 
545
c      write (*,*) 'nmddof, nmdnode',nmddof,nmdnode
 
546
!
 
547
      return
 
548
      end
 
549
 
 
550
 
 
551
 
 
552