CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
gen3dnor.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine gen3dnor (nk, nk_, co, iponoel, inoel, iponoelmax, kon, ipkon, lakon, ne, thicke, offset, iponor, xnor, knor, rig, iperturb, tinc, tper, tmin, tmax, ctrl, ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, labmpc, ikboun, ilboun, nboun, nboun_, nodeboun, ndirboun, xboun, iamboun, typeboun, nam, ntrans, inotr, trab, ikfree, ixfree, nmethod, ithermal, istep, mi, icomposite, ielmat, vold)
 

Function/Subroutine Documentation

◆ gen3dnor()

subroutine gen3dnor ( integer  nk,
integer  nk_,
real*8, dimension(3,*)  co,
integer, dimension(*)  iponoel,
integer, dimension(3,*)  inoel,
integer  iponoelmax,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer  ne,
real*8, dimension(mi(3),*)  thicke,
real*8, dimension(2,*)  offset,
integer, dimension(2,*)  iponor,
real*8, dimension(*)  xnor,
integer, dimension(*)  knor,
integer, dimension(*)  rig,
integer  iperturb,
real*8  tinc,
real*8  tper,
real*8  tmin,
real*8  tmax,
real*8, dimension(*)  ctrl,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer  nmpc_,
integer  mpcfree,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
character*20, dimension(*)  labmpc,
integer, dimension(*)  ikboun,
integer, dimension(*)  ilboun,
integer  nboun,
integer  nboun_,
integer, dimension(*)  nodeboun,
integer, dimension(*)  ndirboun,
real*8, dimension(*)  xboun,
integer, dimension(*)  iamboun,
character*1, dimension(*)  typeboun,
integer  nam,
integer  ntrans,
integer, dimension(2,*)  inotr,
real*8, dimension(7,*)  trab,
integer  ikfree,
integer  ixfree,
integer  nmethod,
integer, dimension(2)  ithermal,
integer  istep,
integer, dimension(*)  mi,
integer  icomposite,
integer, dimension(mi(3),*)  ielmat,
real*8, dimension(0:mi(2),*)  vold 
)
25 !
26 ! calculates normals on 1-D and 2-D elements
27 !
28 ! degrees of freedom:
29 ! 0: Temperature
30 ! 1: tra_in_x
31 ! 2: tra_in_y
32 ! 3: tra_in_z
33 ! 4: pressure (only CFD) or rot_about_x (else)
34 ! 5: rot_about_y
35 ! 6: rot_about_z
36 !
37  implicit none
38 !
39  logical fixed,composite,beam
40 !
41  character*1 type,typeboun(*)
42  character*8 lakon(*)
43  character*20 labmpc(*),label
44 !
45  integer nk,nk_,iponoel(*),inoel(3,*),iponoelmax,kon(*),ipkon(*),
46  & ne,iponor(2,*),knor(*),rig(*),iperturb,ipompc(*),nodempc(3,*),
47  & nmpc,nmpc_,mpcfree,ikmpc(*),ilmpc(*),ikboun(*),ilboun(*),nboun,
48  & nboun_,nodeboun(*),ndirboun(*),iamboun(*),nam,ntrans,inotr(2,*),
49  & isol,istep,idummy,mi(*),icomposite,ielmat(mi(3),*),nkold,
50  & i,ndepnodes,index,nexp,nnor,nel,ielem,indexe,j,iel(100),idmpc,
51  & jl(100),ial(100),ifi(100),idepnodes(800),indexx,k,l,ifix,nemin,
52  & jact,ixfree,ikfree,node,nelshell,irefnode,idof,id,mpcfreeold,
53  & irotnode,imax,iamplitude,nmethod,ithermal(2),iexpnode,idim
54 !
55  real*8 co(3,*),thicke(mi(3),*),offset(2,*),xnor(*),tinc,tper,tmin,
56  & tmax,ctrl(*),coefmpc(*),xboun(*),trab(7,*),vold(0:mi(2),*),
57  & xno(3,100),xta(3,100),xn1(3,100),thl1(100),thl2(100),
58  & off1(100),off2(100),xi,et,coloc6(2,6),coloc8(2,8),xl(3,8),
59  & dd,xnoref(3),dot,coloc3(3),dot1,dot2,dmax,val,coloc2(2),
60  & e1(3),e2(3),t1(3)
61 !
62  data coloc2 /-1.d0,1.d0/
63  data coloc3 /-1.d0,0.d0,1.d0/
64  data coloc6 /0.d0,0.d0,1.d0,0.d0,0.d0,1.d0,0.5d0,0.d0,
65  & 0.5d0,0.5d0,0.d0,0.5d0/
66  data coloc8 /-1.d0,-1.d0,1.d0,-1.d0,1.d0,1.d0,-1.d0,1.d0,
67  & 0.d0,-1.d0,1.d0,0.d0,0.d0,1.d0,-1.d0,0.d0/
68 !
69  label=' '
70  fixed=.false.
71 !
72 ! calculating the normals in nodes belonging to shells/beams
73 !
74  nkold=nk
75 !
76  do i=1,nkold
77  ndepnodes=0
78  idim=0
79  index=iponoel(i)
80  if(index.eq.0) cycle
81 !
82 ! nexp indicates how many times the node was expanded
83 !
84  nexp=0
85 !
86 ! nnor indicates whether the expanded nodes lie on a point
87 ! (nnor=0, only for plane stress, plane strain or axisymmetric
88 ! elements), on a line (nnor=1) or in a plane (nnor=2)
89 !
90  nnor=0
91 !
92 ! locating the shell elements to which node i belongs
93 !
94  nel=0
95  do
96  if(index.eq.0) exit
97  ielem=inoel(1,index)
98  if((lakon(ielem)(1:1).ne.'B').and.
99  & (lakon(ielem)(1:1).ne.'T')) then
100  if(lakon(ielem)(1:1).eq.'S') nnor=1
101  indexe=ipkon(ielem)
102  nel=nel+1
103  if(nel.gt.100) then
104  write(*,*) '*ERROR in gen3dnor: more than 100'
105  write(*,*) ' shell elements share the'
106  write(*,*) ' same node'
107  call exit(201)
108  endif
109  j=inoel(2,index)
110  jl(nel)=j
111  iel(nel)=ielem
112  thl1(nel)=thicke(1,indexe+j)
113 !
114 ! correcting the thickness for more than one layer (composite)
115 !
116  do k=2,mi(3)
117  thl1(nel)=thl1(nel)+thicke(k,indexe+j)
118  enddo
119 !
120  off1(nel)=offset(1,ielem)
121  endif
122  index=inoel(3,index)
123  enddo
124 !
125  if(nel.gt.0) then
126  do j=1,nel
127  ial(j)=0
128  enddo
129 !
130 ! estimate the normal
131 !
132  do j=1,nel
133  indexe=ipkon(iel(j))
134  indexx=iponor(1,indexe+jl(j))
135  if(indexx.ge.0) then
136  do k=1,3
137  xno(k,j)=xnor(indexx+k)
138  enddo
139  ifi(j)=1
140  cycle
141  else
142  ifi(j)=0
143  endif
144 !
145 ! local normal on the element (Jacobian)
146 !
147  if((lakon(iel(j))(2:2).eq.'3').or.
148  & (lakon(iel(j))(4:4).eq.'3')) then
149  xi=coloc6(1,jl(j))
150  et=coloc6(2,jl(j))
151  do k=1,3
152  node=kon(indexe+k)
153  do l=1,3
154  xl(l,k)=co(l,node)
155  enddo
156  enddo
157  call norshell3(xi,et,xl,xno(1,j))
158  elseif((lakon(iel(j))(2:2).eq.'4').or.
159  & (lakon(iel(j))(4:4).eq.'4')) then
160  xi=coloc8(1,jl(j))
161  et=coloc8(2,jl(j))
162  do k=1,4
163  node=kon(indexe+k)
164  do l=1,3
165  xl(l,k)=co(l,node)
166  enddo
167  enddo
168  call norshell4(xi,et,xl,xno(1,j))
169  elseif((lakon(iel(j))(2:2).eq.'6').or.
170  & (lakon(iel(j))(4:4).eq.'6')) then
171  xi=coloc6(1,jl(j))
172  et=coloc6(2,jl(j))
173  do k=1,6
174  node=kon(indexe+k)
175  do l=1,3
176  xl(l,k)=co(l,node)
177  enddo
178  enddo
179  call norshell6(xi,et,xl,xno(1,j))
180  elseif((lakon(iel(j))(2:2).eq.'8').or.
181  & (lakon(iel(j))(4:4).eq.'8')) then
182  xi=coloc8(1,jl(j))
183  et=coloc8(2,jl(j))
184  do k=1,8
185  node=kon(indexe+k)
186  do l=1,3
187  xl(l,k)=co(l,node)
188  enddo
189  enddo
190  call norshell8(xi,et,xl,xno(1,j))
191  endif
192 !
193  dd=dsqrt(xno(1,j)**2+xno(2,j)**2+xno(3,j)**2)
194  if(dd.lt.1.d-10) then
195  write(*,*) '*ERROR in gen3dnor: size of estimated'
196  write(*,*) ' shell normal in node ',i,
197  & ' element ',iel(j)
198  write(*,*) ' is smaller than 1.e-10'
199  call exit(201)
200  endif
201  do k=1,3
202  xno(k,j)=xno(k,j)/dd
203  enddo
204  enddo
205 !
206  do
207 !
208 ! determining a fixed normal which was not treated yet,
209 ! or, if none is left, the minimum element number of all
210 ! elements containing node i and for which no normal was
211 ! determined yet
212 !
213 ! if ial(j)=0: the normal on this element has not been
214 ! treated yet
215 ! if ial(j)=2: normal has been treated
216 !
217  ifix=0
218  nemin=ne+1
219  do j=1,nel
220  if(ial(j).ne.0) cycle
221  if(ifi(j).eq.1) then
222  jact=j
223  ifix=1
224  exit
225  endif
226  enddo
227  if(ifix.eq.0) then
228  do j=1,nel
229  if(ial(j).eq.0) then
230  if(iel(j).lt.nemin) then
231  nemin=iel(j)
232  jact=j
233  endif
234  endif
235  enddo
236  if(nemin.eq.ne+1) exit
237  endif
238 !
239  do j=1,3
240  xnoref(j)=xno(j,jact)
241  enddo
242 !
243 ! determining all elements whose normal in node i makes an
244 ! angle smaller than 0.5 or 20 degrees with the reference normal,
245 ! depending whether the reference normal was given by the
246 ! user or is being calculated; the thickness and offset must
247 ! also fit.
248 !
249 ! if ial(j)=1: normal on element is being treated now
250 !
251  do j=1,nel
252  if(ial(j).eq.2) cycle
253  if(j.eq.jact) then
254  ial(jact)=1
255  else
256  dot=xno(1,j)*xnoref(1)+xno(2,j)*xnoref(2)+
257  & xno(3,j)*xnoref(3)
258  if(ifix.eq.0) then
259  if(dot.gt.0.939693d0)then
260  if((dabs(thl1(j)-thl1(jact)).lt.1.d-10)
261  & .and.
262  & (dabs(off1(j)-off1(jact)).lt.1.d-10)
263  & .and.
264  & (lakon(iel(j))(1:1).ne.'M').and.
265  & (lakon(iel(jact))(1:1).ne.'M').and.
266  & ((lakon(iel(j))(1:3).eq.lakon(iel(jact))(1:3))
267  & .or.
268  & ((lakon(iel(j))(1:1).eq.'S').and.
269  & (lakon(iel(jact))(1:1).eq.'S'))))
270  & ial(j)=1
271 c
272  if(dot.lt.0.999962d0) nnor=2
273 c
274  else
275  if((lakon(iel(j))(1:1).eq.'S').and.
276  & (lakon(iel(jact))(1:1).eq.'S')) then
277 !
278 ! if the normals have the opposite
279 ! direction, the expanded nodes are on a
280 ! straight line
281 !
282  if(dot.gt.-0.999962) then
283  nnor=2
284  else
285  write(*,*) '*INFO in gen3dnor: in some
286  & nodes opposite normals are defined'
287  endif
288  endif
289  endif
290  else
291  if(dot.gt.0.999962d0) then
292  if((dabs(thl1(j)-thl1(jact)).lt.1.d-10)
293  & .and.
294  & (dabs(off1(j)-off1(jact)).lt.1.d-10)
295  & .and.
296  & (lakon(iel(j))(1:1).ne.'M').and.
297  & (lakon(iel(jact))(1:1).ne.'M').and.
298  & ((lakon(iel(j))(1:3).eq.lakon(iel(jact))(1:3))
299  & .or.
300  & ((lakon(iel(j))(1:1).eq.'S').and.
301  & (lakon(iel(jact))(1:1).eq.'S'))))
302  & ial(j)=1
303 c
304 c if(dot.lt.0.999962) nnor=2
305 c
306  else
307  if((lakon(iel(j))(1:1).eq.'S').and.
308  & (lakon(iel(jact))(1:1).eq.'S')) then
309 !
310 ! if the normals have the opposite
311 ! direction, the expanded nodes are on a
312 ! straight line
313 !
314  if(dot.gt.-0.999962) then
315  nnor=2
316  else
317  write(*,*) '*INFO in gen3dnor: in some
318  & nodes opposite normals are defined'
319  endif
320  endif
321  endif
322  endif
323  endif
324  enddo
325 !
326 ! determining the mean normal for the selected elements
327 !
328  if(ifix.eq.0) then
329  do j=1,3
330  xnoref(j)=0.d0
331  enddo
332  do j=1,nel
333  if(ial(j).eq.1) then
334  do k=1,3
335  xnoref(k)=xnoref(k)+xno(k,j)
336  enddo
337  endif
338  enddo
339  dd=dsqrt(xnoref(1)**2+xnoref(2)**2+xnoref(3)**2)
340  if(dd.lt.1.d-10) then
341  write(*,*) '*ERROR in gen3dnor: size of'
342  write(*,*) ' estimated shell normal is'
343  write(*,*) ' smaller than 1.e-10'
344  call exit(201)
345  endif
346  do j=1,3
347  xnoref(j)=xnoref(j)/dd
348  enddo
349  endif
350 !
351 ! check whether any elements are composite elements
352 !
353  composite=.false.
354  if(icomposite.eq.1) then
355  do j=1,nel
356  if(ial(j).eq.1) then
357  ielem=iel(j)
358  if(ielmat(2,ielem).ne.0) then
359  composite=.true.
360  exit
361  endif
362  endif
363  enddo
364  endif
365 !
366 ! updating the pointers iponor
367 !
368  nexp=nexp+1
369  do j=1,nel
370  if(ial(j).eq.1) then
371  ial(j)=2
372  if(ifix.eq.0) then
373  iponor(1,ipkon(iel(j))+jl(j))=ixfree
374  elseif(j.ne.jact) then
375  iponor(1,ipkon(iel(j))+jl(j))=
376  & iponor(1,ipkon(iel(jact))+jl(jact))
377  endif
378  iponor(2,ipkon(iel(j))+jl(j))=ikfree
379  endif
380  enddo
381 !
382 ! storing the normal in xnor and generating 3 nodes
383 ! for knor
384 !
385  if(ifix.eq.0) then
386  do j=1,3
387  xnor(ixfree+j)=xnoref(j)
388  enddo
389  ixfree=ixfree+3
390  endif
391 !
392  do k=1,3
393  nk=nk+1
394  if(nk.gt.nk_) then
395  write(*,*) '*ERROR in gen3dnor: increase nk_'
396  call exit(201)
397  endif
398  knor(ikfree+k)=nk
399 !
400 ! for plane stress, plane strain and axisymmetric
401 ! elements only the middle node is included in the
402 ! rigid body definition
403 !
404  if((lakon(iel(jact))(2:2).ne.'P').and.
405  & (lakon(iel(jact))(2:2).ne.'A').and.
406  & (lakon(iel(jact))(1:1).ne.'M')) then
407  idepnodes(ndepnodes+1)=nk
408  ndepnodes=ndepnodes+1
409  idim=max(idim,1)
410  elseif(k.eq.2) then
411  if (lakon(iel(jact))(1:1).ne.'M') then
412 !
413 ! plane stress/strain axisymmetric
414 !
415  idepnodes(ndepnodes+1)=nk
416  ndepnodes=ndepnodes+1
417  idim=max(idim,1)
418  else
419 !
420 ! membrane: insert hinge
421 !
422  call gen3membrane(ipompc,nodempc,coefmpc,nmpc,
423  & nmpc_,mpcfree,ikmpc,ilmpc,labmpc,nk,
424  & ithermal,i)
425  endif
426  endif
427  enddo
428  ikfree=ikfree+3
429 !
430 ! extra nodes for composite shells
431 !
432  if(composite) then
433  do l=1,mi(3)
434  do k=1,3
435  nk=nk+1
436  if(nk.gt.nk_) then
437  write(*,*) '*ERROR in gen3dnor: increase nk_'
438  call exit(201)
439  endif
440  knor(ikfree+k)=nk
441  enddo
442  ikfree=ikfree+3
443  enddo
444  endif
445 !
446  enddo
447  endif
448 !
449  nelshell=nel+1
450 !
451 ! locating the beam elements to which node i belongs
452 !
453  index=iponoel(i)
454  do
455  if(index.eq.0) exit
456  ielem=inoel(1,index)
457  if((lakon(ielem)(1:1).eq.'B').or.
458  & (lakon(ielem)(1:1).eq.'T')) then
459  if(lakon(ielem)(1:1).eq.'B') then
460  beam=.true.
461  else
462  beam=.false.
463  endif
464  indexe=ipkon(ielem)
465  nel=nel+1
466  if(nel.gt.100) then
467  write(*,*) '*ERROR in gen3dnor: more than 100'
468  write(*,*) ' beam/shell elements share'
469  write(*,*) ' the same node'
470  call exit(201)
471  endif
472  j=inoel(2,index)
473  jl(nel)=j
474  iel(nel)=ielem
475  thl1(nel)=thicke(1,indexe+j)
476  thl2(nel)=thicke(2,indexe+j)
477  off1(nel)=offset(1,ielem)
478  off2(nel)=offset(2,ielem)
479  endif
480  index=inoel(3,index)
481  enddo
482 !
483  if(nel.ge.nelshell) then
484  if(beam) nnor=2
485  do j=nelshell,nel
486  ial(j)=0
487  enddo
488 !
489 ! estimate the normal
490 !
491  do j=nelshell,nel
492  if((lakon(iel(j))(3:3).eq.'1').or.
493  & (lakon(iel(j))(4:4).eq.'2')) then
494 !
495 ! linear beam or truss element
496 !
497  xi=coloc2(jl(j))
498  indexe=ipkon(iel(j))
499  do k=1,2
500  node=kon(indexe+k)
501  do l=1,3
502  xl(l,k)=co(l,node)
503  enddo
504  enddo
505 !
506 ! determining the tangent vector xta
507 !
508  do k=1,3
509  xta(k,j)=(xl(k,2)-xl(k,1))/2.d0
510  enddo
511  else
512 !
513 ! quadratic beam element
514 !
515  xi=coloc3(jl(j))
516  indexe=ipkon(iel(j))
517  do k=1,3
518  node=kon(indexe+k)
519  do l=1,3
520  xl(l,k)=co(l,node)
521  enddo
522  enddo
523 !
524 ! determining the tangent vector xta
525 !
526  do k=1,3
527  xta(k,j)=(xi-0.5d0)*xl(k,1)-2.d0*xi*xl(k,2)+
528  & (xi+0.5d0)*xl(k,3)
529  enddo
530  endif
531 !
532  dd=dsqrt(xta(1,j)**2+xta(2,j)**2+xta(3,j)**2)
533  if(dd.lt.1.d-10) then
534  write(*,*) '*ERROR in gen3dnor: size of estimated'
535  write(*,*)' beam tangent in node ',i,' element '
536  &,iel(j)
537  write(*,*) ' is smaller than 1.e-10'
538  call exit(201)
539  endif
540  do k=1,3
541  xta(k,j)=xta(k,j)/dd
542  enddo
543 !
544 ! check whether normal was defined with *NORMAL and
545 ! determine vector n1
546 !
547  if(iponor(1,indexe+jl(j)).ge.0) then
548  indexx=iponor(1,indexe+jl(j))
549  if(dabs(xnor(indexx+4)**2+xnor(indexx+5)**2+
550  & xnor(indexx+6)**2-1.d0).lt.1.d-5) then
551  do k=1,3
552  xno(k,j)=xnor(indexx+3+k)
553  enddo
554  ifi(j)=1
555  cycle
556  endif
557  ifi(j)=0
558  do k=1,3
559  xn1(k,j)=xnor(indexx+k)
560  enddo
561  else
562  ifi(j)=0
563  xn1(1,j)=0.d0
564  xn1(2,j)=0.d0
565  xn1(3,j)=-1.d0
566  endif
567 !
568 ! normal (=n2) = xta x xn1
569 !
570  xno(1,j)=xta(2,j)*xn1(3,j)-xta(3,j)*xn1(2,j)
571  xno(2,j)=xta(3,j)*xn1(1,j)-xta(1,j)*xn1(3,j)
572  xno(3,j)=xta(1,j)*xn1(2,j)-xta(2,j)*xn1(1,j)
573  dd=dsqrt(xno(1,j)**2+xno(2,j)**2+xno(3,j)**2)
574  if(dd.lt.1.d-10) then
575  write(*,*) '*ERROR in gen3dnor: size of estimated'
576  write(*,*)' beam normal in 2-direction in node '
577  &,i,' element ',iel(j)
578  write(*,*) ' is smaller than 1.e-10'
579  call exit(201)
580  endif
581  do k=1,3
582  xno(k,j)=xno(k,j)/dd
583  enddo
584  enddo
585 !
586  do
587 !
588 ! determining a fixed normal which was not treated yet,
589 ! or, if none is left, the minimum element number of all
590 ! elements containing node i and for which no normal was
591 ! determined yet
592 !
593  ifix=0
594  nemin=ne+1
595  do j=nelshell,nel
596  if(ial(j).ne.0) cycle
597  if(ifi(j).eq.1) then
598  jact=j
599  ifix=1
600  exit
601  endif
602  enddo
603  if(ifix.eq.0) then
604  do j=nelshell,nel
605  if(ial(j).eq.0) then
606  if(iel(j).lt.nemin) then
607  nemin=iel(j)
608  jact=j
609  endif
610  endif
611  enddo
612  if(nemin.eq.ne+1) exit
613  endif
614 !
615 ! the reference normal is the one on the element with the
616 ! smallest element number
617 !
618  do j=1,3
619  xnoref(j)=xno(j,jact)
620  enddo
621 !
622 ! determining all elements whose normal in node i makes an
623 ! angle smaller than 0.5 or 20 degrees with the reference normal,
624 ! depending whether the reference normal was given by the
625 ! user or is being calculated; the thickness and offset must
626 ! also fit.
627 !
628  do j=nelshell,nel
629  if(j.eq.jact) then
630  ial(jact)=1
631  else
632  dot1=xno(1,j)*xnoref(1)+xno(2,j)*xnoref(2)+
633  & xno(3,j)*xnoref(3)
634  dot2=xta(1,j)*xta(1,jact)+xta(2,j)*xta(2,jact)+
635  & xta(3,j)*xta(3,jact)
636  if(ifix.eq.0) then
637  if((dot1.gt.0.939693d0).and.
638  & (dot2.gt.0.939693d0)) then
639  if((dabs(thl1(j)-thl1(jact)).lt.1.d-10)
640  & .and.
641  & (dabs(thl2(j)-thl2(jact)).lt.1.d-10)
642  & .and.
643  & (dabs(off1(j)-off1(jact)).lt.1.d-10)
644  & .and.
645  & (dabs(off2(j)-off2(jact)).lt.1.d-10)
646  & .and.
647  & (lakon(iel(j))(1:1).ne.'T').and.
648  & (lakon(iel(jact))(1:1).ne.'T').and.
649  & (lakon(iel(j))(8:8).eq.lakon(iel(jact))(8:8)))
650  & ial(j)=1
651  endif
652  else
653  if((dot1.gt.0.999962d0).and.
654  & (dot2.gt.0.999962d0)) then
655  if((dabs(thl1(j)-thl1(jact)).lt.1.d-10)
656  & .and.
657  & (dabs(thl2(j)-thl2(jact)).lt.1.d-10)
658  & .and.
659  & (dabs(off1(j)-off1(jact)).lt.1.d-10)
660  & .and.
661  & (dabs(off2(j)-off2(jact)).lt.1.d-10)
662  & .and.
663  & (lakon(iel(j))(1:1).ne.'T').and.
664  & (lakon(iel(jact))(1:1).ne.'T').and.
665  & (lakon(iel(j))(8:8).eq.lakon(iel(jact))(8:8)))
666  & ial(j)=1
667  endif
668  endif
669  endif
670  enddo
671 !
672 ! determining the mean normal for the selected elements
673 !
674  if(ifix.eq.0) then
675  do j=1,3
676  xnoref(j)=0.d0
677  enddo
678  do j=nelshell,nel
679  if(ial(j).eq.1) then
680  do k=1,3
681  xnoref(k)=xnoref(k)+xno(k,j)
682  enddo
683  endif
684  enddo
685  endif
686 !
687 ! calculating the mean tangent
688 !
689  do j=nelshell,nel
690  if((ial(j).eq.1).and.(j.ne.jact)) then
691  do k=1,3
692  xta(k,jact)=xta(k,jact)+xta(k,j)
693  enddo
694  endif
695  enddo
696  dd=dsqrt(xta(1,jact)**2+xta(2,jact)**2+xta(3,jact)**2)
697  if(dd.lt.1.d-10) then
698  write(*,*) '*ERROR in gen3dnor: size of mean'
699  write(*,*)' beam tangent is smaller than 1.e-10'
700  call exit(201)
701  endif
702  do k=1,3
703  xta(k,jact)=xta(k,jact)/dd
704  enddo
705 !
706 ! taking care that the mean normal is orthogonal towards
707 ! the mean tangent
708 !
709  dd=xnoref(1)*xta(1,jact)+xnoref(2)*xta(2,jact)+
710  & xnoref(3)*xta(3,jact)
711  do j=1,3
712  xnoref(j)=xnoref(j)-dd*xta(j,jact)
713  enddo
714  dd=dsqrt(xnoref(1)**2+xnoref(2)**2+xnoref(3)**2)
715  if(dd.lt.1.d-10) then
716  write(*,*) '*ERROR in gen3dnor: size of'
717  write(*,*) ' estimated beam normal is'
718  write(*,*) ' smaller than 1.e-10'
719  call exit(201)
720  endif
721  do j=1,3
722  xnoref(j)=xnoref(j)/dd
723  enddo
724 !
725 ! calculating xn1 = xn2 x tangent
726 !
727  xn1(1,jact)=xnoref(2)*xta(3,jact)-xnoref(3)*xta(2,jact)
728  xn1(2,jact)=xnoref(3)*xta(1,jact)-xnoref(1)*xta(3,jact)
729  xn1(3,jact)=xnoref(1)*xta(2,jact)-xnoref(2)*xta(1,jact)
730 !
731 ! storing the normals in xnor and generating 8 nodes
732 ! for knor
733 !
734  nexp=nexp+1
735  do j=nelshell,nel
736  if(ial(j).eq.1) then
737  ial(j)=2
738  if(ifix.eq.0) then
739  iponor(1,ipkon(iel(j))+jl(j))=ixfree
740  else
741  iponor(1,ipkon(iel(j))+jl(j))=
742  & iponor(1,ipkon(iel(jact))+jl(jact))
743  endif
744  iponor(2,ipkon(iel(j))+jl(j))=ikfree
745  endif
746  enddo
747 !
748  do j=1,3
749  xnor(ixfree+j)=xn1(j,jact)
750  enddo
751  do j=1,3
752  xnor(ixfree+3+j)=xnoref(j)
753  enddo
754  ixfree=ixfree+6
755  if(lakon(iel(jact))(1:1).ne.'T') then
756 !
757 ! beam
758 !
759  do k=1,8
760  nk=nk+1
761  if(nk.gt.nk_) then
762  write(*,*) '*ERROR in gen3dnor: increase nk_'
763  call exit(201)
764  endif
765  knor(ikfree+k)=nk
766  idepnodes(ndepnodes+k)=nk
767  enddo
768  ndepnodes=ndepnodes+8
769  idim=max(idim,3)
770  else
771 !
772 ! truss
773 !
774  do k=1,8
775  nk=nk+1
776  if(nk.gt.nk_) then
777  write(*,*) '*ERROR in gen3dnor: increase nk_'
778  call exit(201)
779  endif
780  knor(ikfree+k)=nk
781 !
782 ! assigning coordinates (necessary for the
783 ! no-rotation-about-the-truss-axis-MPC created
784 ! in gen3dtruss
785 !
786  if(k.eq.1) then
787  do j=1,3
788  co(j,nk)=co(j,i)
789  & -thl1(jact)*xn1(j,jact)*(.5d0+off1(jact))
790  & +thl2(jact)*xnoref(j)*(.5d0-off2(jact))
791  enddo
792  elseif(k.eq.2) then
793  do j=1,3
794  co(j,nk)=co(j,i)
795  & -thl1(jact)*xn1(j,jact)*(.5d0+off1(jact))
796  & -thl2(jact)*xnoref(j)*(.5d0+off2(jact))
797  enddo
798  elseif(k.eq.3) then
799  do j=1,3
800  co(j,nk)=co(j,i)
801  & +thl1(jact)*xn1(j,jact)*(.5d0-off1(jact))
802  & -thl2(jact)*xnoref(j)*(.5d0+off2(jact))
803  enddo
804  elseif(k.eq.4) then
805  do j=1,3
806  co(j,nk)=co(j,i)
807  & +thl1(jact)*xn1(j,jact)*(.5d0-off1(jact))
808  & +thl2(jact)*xnoref(j)*(.5d0-off2(jact))
809  enddo
810  endif
811  enddo
812 !
813  call gen3dtruss(ipompc,nodempc,coefmpc,nmpc,
814  & nmpc_,mpcfree,ikmpc,ilmpc,labmpc,nk,
815  & ithermal,i,nodeboun,ndirboun,ikboun,ilboun,
816  & nboun,nboun_,typeboun,xboun,xta,jact,co,
817  & knor,ntrans,inotr,trab,vold,mi,nmethod,nk_,
818  & nam,iperturb,ikfree,iamboun)
819  endif
820  ikfree=ikfree+8
821  enddo
822  endif
823 !
824 ! check whether the user has specified rotational degrees
825 ! of freedom (in that case rig(i)=-1 was assigned in
826 ! subroutine gen3delem); if so, a rigid MPC must be defined
827 !
828  if(rig(i).ne.0) then
829  rig(i)=0
830  if(nexp.le.1) then
831  nexp=2
832  endif
833  endif
834 !
835 ! storing the expanded nodes
836 !
837 !
838 ! generate rigid MPC's if necessary
839 !
840  if(nexp.gt.1) then
841  irefnode=i
842 !
843  rig(i)=-1
844 !
845  if(ithermal(2).ne.2) then
846  if(nnor.eq.0) then
847 !
848 ! the node belongs to plane stress, plane strain
849 ! or axisymmetric elements only. These are only linked
850 ! through the node in the midplane: the nodes
851 ! coincide; only DOF1 and DOF2 are linked.
852 ! rig(i)=-1 to indicate that a knot has
853 ! been generated without rotational node
854 !
855  do k=1,ndepnodes
856  node=idepnodes(k)
857  do j=1,2
858  idof=8*(node-1)+j
859  call nident(ikmpc,idof,nmpc,id)
860  nmpc=nmpc+1
861  if(nmpc.gt.nmpc_) then
862  write(*,*)
863  & '*ERROR in rigidmpc: increase nmpc_'
864  call exit(201)
865  endif
866 !
867  ipompc(nmpc)=mpcfree
868  labmpc(nmpc)=' '
869 !
870  do l=nmpc,id+2,-1
871  ikmpc(l)=ikmpc(l-1)
872  ilmpc(l)=ilmpc(l-1)
873  enddo
874  ikmpc(id+1)=idof
875  ilmpc(id+1)=nmpc
876 !
877  nodempc(1,mpcfree)=node
878  nodempc(2,mpcfree)=j
879  coefmpc(mpcfree)=1.d0
880  mpcfree=nodempc(3,mpcfree)
881  nodempc(1,mpcfree)=irefnode
882  nodempc(2,mpcfree)=j
883  coefmpc(mpcfree)=-1.d0
884  mpcfreeold=mpcfree
885  mpcfree=nodempc(3,mpcfree)
886  nodempc(3,mpcfreeold)=0
887  enddo
888  enddo
889  else
890 !
891 ! generate a rigid body knot; rig(i) contains the
892 ! rotational node of the knot
893 !
894  nk=nk+1
895  if(nk.gt.nk_) then
896  write(*,*) '*ERROR in rigidbodies: increase nk_'
897  call exit(201)
898  endif
899  irotnode=nk
900  rig(i)=irotnode
901  nk=nk+1
902  if(nk.gt.nk_) then
903  write(*,*) '*ERROR in rigidbodies: increase nk_'
904  call exit(201)
905  endif
906  iexpnode=nk
907  do k=1,ndepnodes
908  call knotmpc(ipompc,nodempc,coefmpc,irefnode,
909  & irotnode,iexpnode,
910  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,
911  & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,
912  & idepnodes,typeboun,co,xboun,istep,k,
913  & ndepnodes,idim,e1,e2,t1)
914  enddo
915  endif
916  endif
917 !
918 ! MPC's for the temperature DOFs
919 !
920  if(ithermal(2).ge.2) then
921  do k=1,ndepnodes
922  node=idepnodes(k)
923  idof=8*(node-1)
924  call nident(ikmpc,idof,nmpc,id)
925  nmpc=nmpc+1
926  if(nmpc.gt.nmpc_) then
927  write(*,*)
928  & '*ERROR in gen3dnor: increase nmpc_'
929  call exit(201)
930  endif
931 !
932  ipompc(nmpc)=mpcfree
933  labmpc(nmpc)=' '
934 !
935  do l=nmpc,id+2,-1
936  ikmpc(l)=ikmpc(l-1)
937  ilmpc(l)=ilmpc(l-1)
938  enddo
939  ikmpc(id+1)=idof
940  ilmpc(id+1)=nmpc
941 !
942  nodempc(1,mpcfree)=node
943  nodempc(2,mpcfree)=0
944  coefmpc(mpcfree)=1.d0
945  mpcfree=nodempc(3,mpcfree)
946  nodempc(1,mpcfree)=irefnode
947  nodempc(2,mpcfree)=0
948  coefmpc(mpcfree)=-1.d0
949  mpcfreeold=mpcfree
950  mpcfree=nodempc(3,mpcfree)
951  nodempc(3,mpcfreeold)=0
952  enddo
953  endif
954 !
955  if((nnor.eq.1).and.(ithermal(2).ne.2)) then
956 !
957 ! generate an additional SPC or MPC for rigid body nodes
958 ! lying on a line to prevent rotation about the
959 ! line
960 !
961  dmax=0.d0
962  imax=0
963  do j=1,3
964  if(dabs(xnoref(j)).gt.dmax) then
965  dmax=dabs(xnoref(j))
966  imax=j
967  endif
968  enddo
969 !
970 ! check whether a SPC suffices
971 !
972  if(dabs(1.d0-dmax).lt.1.d-3) then
973  val=0.d0
974  if(nam.gt.0) iamplitude=0
975  type='R'
976 !
977 ! check whether the dof has been used as dependent
978 ! term of a MPC
979 !
980  idof=8*(i-1)+3+imax
981  call nident(ikmpc,idof,nmpc,idmpc)
982 !
983  if((idmpc.le.0).or.(ikmpc(idmpc).ne.idof)) then
984  call bounadd(irotnode,imax,imax,val,nodeboun,
985  & ndirboun,xboun,nboun,nboun_,iamboun,
986  & iamplitude,nam,ipompc,nodempc,coefmpc,
987  & nmpc,nmpc_,mpcfree,inotr,trab,ntrans,
988  & ikboun,ilboun,ikmpc,ilmpc,co,nk,nk_,labmpc,
989  & type,typeboun,nmethod,iperturb,fixed,vold,
990  & idummy,mi,label)
991  endif
992  else
993 !
994 ! check whether the rotational degree of freedom
995 ! imax is fixed through a SPC
996 !
997  isol=0
998  do l=1,3
999  idof=8*(i-1)+3+imax
1000  call nident(ikboun,idof,nboun,id)
1001  call nident(ikmpc,idof,nmpc,idmpc)
1002  if(((idmpc.gt.0).and.(ikmpc(idmpc).eq.idof)).or.
1003  & ((id.gt.0).and.(ikboun(id).eq.idof)).or.
1004  & (dabs(xnoref(imax)).lt.1.d-2)) then
1005 c changed 24.11.2014
1006 c & (dabs(xnoref(imax)).lt.1.d-10)) then
1007  imax=imax+1
1008  if(imax.gt.3) imax=imax-3
1009  cycle
1010  endif
1011  isol=1
1012  exit
1013  enddo
1014 !
1015 ! if one of the rotational dofs was not used so far,
1016 ! it can be taken as dependent side for fixing the
1017 ! rotation about the normal. If all dofs were used,
1018 ! no additional equation is needed.
1019 !
1020  if(isol.eq.1) then
1021  idof=8*(irotnode-1)+imax
1022  call nident(ikmpc,idof,nmpc,id)
1023  nmpc=nmpc+1
1024  if(nmpc.gt.nmpc_) then
1025  write(*,*)
1026  & '*ERROR in gen3dnor: increase nmpc_'
1027  call exit(201)
1028  endif
1029 !
1030  ipompc(nmpc)=mpcfree
1031  labmpc(nmpc)=' '
1032 !
1033  do l=nmpc,id+2,-1
1034  ikmpc(l)=ikmpc(l-1)
1035  ilmpc(l)=ilmpc(l-1)
1036  enddo
1037  ikmpc(id+1)=idof
1038  ilmpc(id+1)=nmpc
1039 !
1040  nodempc(1,mpcfree)=irotnode
1041  nodempc(2,mpcfree)=imax
1042  coefmpc(mpcfree)=xnoref(imax)
1043  mpcfree=nodempc(3,mpcfree)
1044  imax=imax+1
1045  if(imax.gt.3) imax=imax-3
1046  nodempc(1,mpcfree)=irotnode
1047  nodempc(2,mpcfree)=imax
1048  coefmpc(mpcfree)=xnoref(imax)
1049  mpcfree=nodempc(3,mpcfree)
1050  imax=imax+1
1051  if(imax.gt.3) imax=imax-3
1052  nodempc(1,mpcfree)=irotnode
1053  nodempc(2,mpcfree)=imax
1054  coefmpc(mpcfree)=xnoref(imax)
1055  mpcfreeold=mpcfree
1056  mpcfree=nodempc(3,mpcfree)
1057  nodempc(3,mpcfreeold)=0
1058  endif
1059  endif
1060  endif
1061  endif
1062  enddo
1063 !
1064  return
knotmpc
subroutine knotmpc(ipompc, nodempc, coefmpc, irefnode, irotnode, iexpnode, labmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, nk, nk_, nodeboun, ndirboun, ikboun, ilboun, nboun, nboun_, idepnodes, typeboun, co, xboun, istep, k, ndepnodes, idim, e1, e2, t1)
Definition: knotmpc.f:24
max
#define max(a, b)
Definition: cascade.c:32
gen3dtruss
subroutine gen3dtruss(ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, labmpc, nk, ithermal, i, nodeboun, ndirboun, ikboun, ilboun, nboun, nboun_, typeboun, xboun, xta, jact, co, knor, ntrans, inotr, trab, vold, mi, nmethod, nk_, nam, iperturb, indexk, iamboun)
Definition: gen3dtruss.f:24
norshell8
subroutine norshell8(xi, et, xl, xnor)
Definition: norshell8.f:20
norshell4
subroutine norshell4(xi, et, xl, xnor)
Definition: norshell4.f:20
gen3membrane
subroutine gen3membrane(ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, labmpc, nk, ithermal, i)
Definition: gen3dmembrane.f:21
norshell6
subroutine norshell6(xi, et, xl, xnor)
Definition: norshell6.f:20
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
bounadd
subroutine bounadd(node, is, ie, val, nodeboun, ndirboun, xboun, nboun, nboun_, iamboun, iamplitude, nam, ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, inotr, trab, ntrans, ikboun, ilboun, ikmpc, ilmpc, co, nk, nk_, labmpc, type, typeboun, nmethod, iperturb, fixed, vold, nodetrue, mi, label)
Definition: bounadd.f:24
norshell3
subroutine norshell3(xi, et, xl, xnor)
Definition: norshell3.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)