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

Go to the source code of this file.

Functions/Subroutines

subroutine gen3dboun (ikboun, ilboun, nboun, nboun_, nodeboun, ndirboun, xboun, iamboun, typeboun, iponoel, inoel, iponoelmax, kon, ipkon, lakon, ne, iponor, xnor, knor, ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, labmpc, rig, ntrans, inotr, trab, nam, nk, nk_, co, nmethod, iperturb, istep, vold, mi)
 

Function/Subroutine Documentation

◆ gen3dboun()

subroutine gen3dboun ( 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, dimension(*)  iponoel,
integer, dimension(3,*)  inoel,
integer  iponoelmax,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer  ne,
integer, dimension(2,*)  iponor,
real*8, dimension(*)  xnor,
integer, dimension(*)  knor,
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(*)  rig,
integer  ntrans,
integer, dimension(2,*)  inotr,
real*8, dimension(7,*)  trab,
integer  nam,
integer  nk,
integer  nk_,
real*8, dimension(3,*)  co,
integer  nmethod,
integer  iperturb,
integer  istep,
real*8, dimension(0:mi(2),*)  vold,
integer, dimension(*)  mi 
)
24 !
25 ! connects nodes of 1-D and 2-D elements, for which SPC's were
26 ! defined, to the nodes of their expanded counterparts
27 !
28  implicit none
29 !
30  logical fixed
31 !
32  character*1 type,typeboun(*)
33  character*8 lakon(*)
34  character*20 labmpc(*),label
35 !
36  integer ikboun(*),ilboun(*),nboun,nboun_,nodeboun(*),nodeact,
37  & ndirboun(*),idim,ier,matz,ncgnodes,lstart,lend,linc,nnodes,
38  & iamboun(*),iponoel(*),inoel(3,*),iponoelmax,kon(*),ipkon(*),ne,
39  & iponor(2,*),knor(*),ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,
40  & ikmpc(*),ilmpc(*),rig(*),ntrans,inotr(2,*),nbounold,i,node,
41  & index,ielem,j,indexe,indexk,idir,iamplitude,irotnode,nk,nk_,
42  & newnode,idof,id,mpcfreenew,k,nam,nmethod,iperturb,ndepnodes,
43  & idepnodes(80),l,iexpnode,indexx,irefnode,imax,isol,mpcfreeold,
44  & nod,impc,istep,nrhs,ipiv(3),info,m,mi(*),itr,idirref,id1,nnode
45 !
46  real*8 xboun(*),xnor(*),coefmpc(*),trab(7,*),val,co(3,*),
47  & xnoref(3),dmax,d(3,3),e(3,3,3),alpha,q(3),w(3),xn(3),
48  & a1(3),a2(3),dd,c1,c2,c3,ww,c(3,3),vold(0:mi(2),*),a(3,3),
49  & e1(3),e2(3),t1(3),b(3,3),x(3),y(3),fv1(3),dot,
50  & fv2(3),z(3,3),xi1,xi2,xi3,u(3,3),r(3,3),xnode
51 !
52  data d /1.,0.,0.,0.,1.,0.,0.,0.,1./
53  data e /0.,0.,0.,0.,0.,-1.,0.,1.,0.,
54  & 0.,0.,1.,0.,0.,0.,-1.,0.,0.,
55  & 0.,-1.,0.,1.,0.,0.,0.,0.,0./
56 !
57  label=' '
58  fixed=.false.
59 !
60  nbounold=nboun
61  do i=1,nbounold
62  node=nodeboun(i)
63  if(node.gt.iponoelmax) then
64 c if((ndirboun(i).gt.3).and.(ndirboun(i).lt.7)) then
65 c write(*,*) '*ERROR: in gen3dboun: node ',node,
66 c & ' does not'
67 c write(*,*) ' belong to a beam nor shell'
68 c write(*,*) ' element and consequently has no'
69 c write(*,*) ' rotational degrees of freedom.'
70 c write(*,*) ' This may be caused by applying'
71 c write(*,*) ' a local coordinate system to a'
72 c write(*,*) ' node in which a rotational boundary'
73 c write(*,*) ' conditions is defined. In CalculiX'
74 c write(*,*) ' this is not allowed.'
75 c call exit(201)
76 c endif
77  cycle
78  endif
79  index=iponoel(node)
80  if(index.eq.0) then
81 c if((ndirboun(i).gt.3).and.(ndirboun(i).lt.7)) then
82 c write(*,*) '*ERROR: in gen3dboun: node ',node,
83 c & ' does not'
84 c write(*,*) ' belong to a beam nor shell'
85 c write(*,*) ' element and consequently has no'
86 c write(*,*) ' rotational degrees of freedom.'
87 c write(*,*) ' This may be caused by applying'
88 c write(*,*) ' a local coordinate system to a'
89 c write(*,*) ' node in which a rotational boundary'
90 c write(*,*) ' conditions is defined. In CalculiX'
91 c write(*,*) ' this is not allowed.'
92 c call exit(201)
93 c endif
94  cycle
95  endif
96  ielem=inoel(1,index)
97  j=inoel(2,index)
98  indexe=ipkon(ielem)
99  indexk=iponor(2,indexe+j)
100  idir=ndirboun(i)
101  val=xboun(i)
102  if(nam.gt.0) iamplitude=iamboun(i)
103 !
104  if(rig(node).ne.0) then
105 !
106 ! existing knot
107 !
108  if(idir.gt.3) then
109  if(rig(node).lt.0) then
110  write(*,*) '*ERROR in gen3dboun: in node ',node
111  write(*,*) ' a rotational DOF is constrained'
112  write(*,*) ' by a SPC; however, the elements'
113  write(*,*) ' to which this node belongs do not'
114  write(*,*) ' have rotational DOFs'
115  call exit(201)
116  endif
117  j=idir-3
118  irotnode=rig(node)
119  type='B'
120  call bounadd(irotnode,j,j,val,nodeboun,
121  & ndirboun,xboun,nboun,nboun_,iamboun,
122  & iamplitude,nam,ipompc,nodempc,coefmpc,
123  & nmpc,nmpc_,mpcfree,inotr,trab,ntrans,
124  & ikboun,ilboun,ikmpc,ilmpc,co,nk,nk_,labmpc,
125  & type,typeboun,nmethod,iperturb,fixed,vold,
126  & irotnode,mi,label)
127  endif
128  else
129 !
130 ! 1. no existing knot
131 ! rotational dof
132 ! dynamics
133 ! => knot is created
134 !
135 ! check for rotational DOFs defined in any but the first step
136 ! nonlinear dynamic case: creation of knots
137 ! knots (expandable rigid bodies) can take rotational
138 ! values arbitrarily exceeding 90 degrees
139 !
140  if((idir.gt.3).and.((nmethod.eq.4).and.(iperturb.gt.1)))then
141 !
142 ! create a knot: determine the knot
143 !
144  ndepnodes=0
145  if(lakon(ielem)(7:7).eq.'L') then
146  do k=1,3
147  ndepnodes=ndepnodes+1
148  idepnodes(ndepnodes)=knor(indexk+k)
149  enddo
150  idim=1
151  elseif(lakon(ielem)(7:7).eq.'B') then
152  do k=1,8
153  ndepnodes=ndepnodes+1
154  idepnodes(ndepnodes)=knor(indexk+k)
155  enddo
156  idim=3
157  else
158  write(*,*)
159  & '*ERROR in gen3dboun: a rotational DOF was applied'
160  write(*,*)
161  & '* to node',node,' without rotational DOFs'
162  call exit(201)
163  endif
164 !
165 ! remove all MPC's in which the knot nodes are
166 ! dependent nodes
167 !
168  do k=1,ndepnodes
169  nod=idepnodes(k)
170  do l=1,3
171  idof=8*(nod-1)+l
172  call nident(ikmpc,idof,nmpc,id)
173  if(id.gt.0) then
174  if(ikmpc(id).eq.idof) then
175  impc=ilmpc(id)
176  call mpcrem(impc,mpcfree,nodempc,nmpc,
177  & ikmpc,ilmpc,labmpc,coefmpc,ipompc)
178  endif
179  endif
180  enddo
181  enddo
182 !
183 ! generate a rigid body knot
184 !
185  irefnode=node
186  nk=nk+1
187  if(nk.gt.nk_) then
188  write(*,*) '*ERROR in rigidbodies: increase nk_'
189  call exit(201)
190  endif
191  irotnode=nk
192  rig(node)=irotnode
193  nk=nk+1
194  if(nk.gt.nk_) then
195  write(*,*) '*ERROR in rigidbodies: increase nk_'
196  call exit(201)
197  endif
198  iexpnode=nk
199  do k=1,ndepnodes
200  call knotmpc(ipompc,nodempc,coefmpc,irefnode,
201  & irotnode,iexpnode,
202  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,
203  & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,
204  & idepnodes,typeboun,co,xboun,istep,k,ndepnodes,
205  & idim,e1,e2,t1)
206  enddo
207 !
208 ! determine the location of the center of gravity of
209 ! the section and its displacements
210 !
211  do l=1,3
212  q(l)=0.d0
213  w(l)=0.d0
214  enddo
215  if(ndepnodes.eq.3) then
216  do k=1,ndepnodes,2
217  nod=idepnodes(k)
218  do l=1,3
219  q(l)=q(l)+co(l,nod)
220  w(l)=w(l)+vold(l,nod)
221  enddo
222  enddo
223  do l=1,3
224  q(l)=q(l)/2.d0
225  w(l)=w(l)/2.d0
226  enddo
227  else
228  do k=1,ndepnodes
229  nod=idepnodes(k)
230  do l=1,3
231  q(l)=q(l)+co(l,nod)
232  w(l)=w(l)+vold(l,nod)
233  enddo
234  enddo
235  do l=1,3
236  q(l)=q(l)/ndepnodes
237  w(l)=w(l)/ndepnodes
238  enddo
239  endif
240 !
241 ! check whether the displacements are zero
242 !
243  dd=dsqrt(w(1)*w(1)+w(2)*w(2)+w(3)*w(3))
244  if(dd.lt.1.d-20) then
245  do l=1,3
246  vold(l,irefnode)=0.d0
247  vold(l,irotnode)=0.d0
248  vold(l,iexpnode)=0.d0
249  enddo
250  else
251 !
252 ! determine the first displacements of iexpnode
253 !
254  alpha=0.d0
255  ncgnodes=0
256  do k=1,ndepnodes
257  nod=idepnodes(k)
258  dd=(co(1,nod)-q(1))**2
259  & +(co(2,nod)-q(2))**2
260  & +(co(3,nod)-q(3))**2
261  if(dd.lt.1.d-20) then
262  ncgnodes=ncgnodes+1
263  cycle
264  endif
265  alpha=alpha+dsqrt(
266  & ((co(1,nod)+vold(1,nod)-q(1)-w(1))**2
267  & +(co(2,nod)+vold(2,nod)-q(2)-w(2))**2
268  & +(co(3,nod)+vold(3,nod)-q(3)-w(3))**2)/dd)
269  enddo
270  if(ndepnodes-ncgnodes.gt.0) then
271  alpha=alpha/(ndepnodes-ncgnodes)
272  endif
273 !
274 ! determine the displacements of irotnodes
275 !
276  do l=1,3
277 c do m=1,3
278 c a(l,m)=0.d0
279 c enddo
280  xn(l)=0.d0
281  enddo
282 !
283  ncgnodes=0
284  do k=1,ndepnodes
285  nod=idepnodes(k)
286  dd=0.d0
287  do l=1,3
288  a1(l)=co(l,nod)-q(l)
289  a2(l)=vold(l,nod)-w(l)
290  dd=dd+a1(l)*a1(l)
291  enddo
292  dd=dsqrt(dd)
293  if(dd.lt.1.d-10) then
294  ncgnodes=ncgnodes+1
295  cycle
296  endif
297  do l=1,3
298  a1(l)=a1(l)/dd
299  a2(l)=a2(l)/dd
300  enddo
301  xn(1)=xn(1)+(a1(2)*a2(3)-a1(3)*a2(2))
302  xn(2)=xn(2)+(a1(3)*a2(1)-a1(1)*a2(3))
303  xn(3)=xn(3)+(a1(1)*a2(2)-a1(2)*a2(1))
304  enddo
305 !
306  if(ndepnodes-ncgnodes.gt.0) then
307  do l=1,3
308  xn(l)=xn(l)/(ndepnodes-ncgnodes)
309  enddo
310  endif
311 !
312  dd=0.d0
313  do l=1,3
314  dd=dd+xn(l)*xn(l)
315  enddo
316  dd=dsqrt(dd)
317  do l=1,3
318  xn(l)=dasin(dd/alpha)*xn(l)/dd
319  enddo
320 !
321 ! determine the displacements of irefnode
322 !
323  ww=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3))
324 !
325  c1=dcos(ww)
326  if(ww.gt.1.d-10) then
327  c2=dsin(ww)/ww
328  else
329  c2=1.d0
330  endif
331  if(ww.gt.1.d-5) then
332  c3=(1.d0-c1)/ww**2
333  else
334  c3=0.5d0
335  endif
336 !
337 ! rotation matrix c
338 !
339  do k=1,3
340  do l=1,3
341  r(k,l)=c1*d(k,l)+
342  & c2*(e(k,1,l)*xn(1)+e(k,2,l)*xn(2)+
343  & e(k,3,l)*xn(3))+c3*xn(k)*xn(l)
344  enddo
345  enddo
346 !
347 ! copying the displacements
348 !
349  do l=1,3
350  vold(l,irefnode)=w(l)
351  vold(l,irotnode)=xn(l)
352  enddo
353  vold(1,iexpnode)=alpha-1.d0
354 !
355 ! correction of the expansion values for beam sections
356 !
357  if(idim.eq.2) then
358 !
359 ! initializing matrices b and c
360 !
361  do l=1,3
362  do m=1,3
363  b(l,m)=0.d0
364  c(l,m)=0.d0
365  enddo
366  enddo
367 !
368 ! solving a least squares problem to determine
369 !
370 ! start meanrotationmpc
371 ! change: mean rotation MPC instead of KNOT
372 !
373  idirref=idir-3
374 !
375  if(lakon(ielem)(7:7).eq.'L') then
376  lstart=3
377  lend=1
378  linc=-2
379  elseif(lakon(ielem)(7:7).eq.'B') then
380  lstart=4
381  lend=1
382  linc=-1
383  endif
384 !
385 ! check for transformations
386 !
387  if(ntrans.le.0) then
388  itr=0
389  elseif(inotr(1,node).eq.0) then
390  itr=0
391  else
392  itr=inotr(1,node)
393  endif
394 !
395 ! determine a unit vector on the rotation axis
396 !
397 ! the transpose of the deformation gradient:
398 ! c.F^T=b
399 !
400  do k=1,ndepnodes
401  nod=idepnodes(k)
402  do l=1,3
403  x(l)=co(l,nod)-q(l)
404  y(l)=x(l)+vold(l,nod)-w(l)
405  enddo
406  do l=1,3
407  do m=1,3
408  c(l,m)=c(l,m)+x(l)*x(m)
409  b(l,m)=b(l,m)+x(l)*y(m)
410  enddo
411  enddo
412  enddo
413 !
414 ! solving the linear equation system
415 !
416  m=3
417  nrhs=3
418  call dgesv(m,nrhs,c,m,ipiv,b,m,info)
419  if(info.ne.0) then
420  write(*,*) '*ERROR in gen3dforc:'
421  write(*,*) ' singular system of equations'
422  call exit(201)
423  endif
424 !
425 ! now b=F^T
426 !
427 ! constructing the right stretch tensor
428 ! U=F^T.R
429 !
430  do l=1,3
431  do m=l,3
432  u(l,m)=b(l,1)*r(1,m)+b(l,2)*r(2,m)+
433  & b(l,3)*r(3,m)
434  enddo
435  enddo
436  u(2,1)=u(1,2)
437  u(3,1)=u(1,3)
438  u(3,2)=u(2,3)
439 !
440 ! determining the eigenvalues and eigenvectors of U
441 !
442  m=3
443  matz=1
444  ier=0
445  call rs(m,m,u,w,matz,z,fv1,fv2,ier)
446  if(ier.ne.0) then
447  write(*,*)
448  & '*ERROR in knotmpc while calculating the'
449  write(*,*) ' eigenvalues/eigenvectors'
450  call exit(201)
451  endif
452 !
453  if((dabs(w(1)-1.d0).lt.dabs(w(2)-1.d0)).and.
454  & (dabs(w(1)-1.d0).lt.dabs(w(3)-1.d0))) then
455  l=2
456  m=3
457  elseif((dabs(w(2)-1.d0).lt.dabs(w(1)-1.d0)).and.
458  & (dabs(w(2)-1.d0).lt.dabs(w(3)-1.d0))) then
459  l=1
460  m=3
461  else
462  l=1
463  m=2
464  endif
465  xi1=datan2
466  & ((z(1,l)*e2(1)+z(2,l)*e2(2)+z(3,l)*e2(2)),
467  & (z(1,l)*e1(1)+z(2,l)*e1(2)+z(3,l)*e1(2)))
468  xi2=w(l)-1.d0
469  xi3=w(m)-1.d0
470 !
471  vold(1,iexpnode)=xi1
472  vold(2,iexpnode)=xi2
473  vold(3,iexpnode)=xi3
474  endif
475  endif
476 !
477 ! apply the boundary condition
478 !
479  idir=idir-3
480  type='B'
481  call bounadd(irotnode,idir,idir,val,nodeboun,
482  & ndirboun,xboun,nboun,nboun_,iamboun,
483  & iamplitude,nam,ipompc,nodempc,coefmpc,
484  & nmpc,nmpc_,mpcfree,inotr,trab,ntrans,
485  & ikboun,ilboun,ikmpc,ilmpc,co,nk,nk_,labmpc,
486  & type,typeboun,nmethod,iperturb,fixed,vold,
487  & irotnode,mi,label)
488 !
489 ! check for shells whether the rotation about the normal
490 ! on the shell has been eliminated
491 !
492  if(lakon(ielem)(7:7).eq.'L') then
493  indexx=iponor(1,indexe+j)
494  do j=1,3
495  xnoref(j)=xnor(indexx+j)
496  enddo
497  dmax=0.d0
498  imax=0
499  do j=1,3
500  if(dabs(xnoref(j)).gt.dmax) then
501  dmax=dabs(xnoref(j))
502  imax=j
503  endif
504  enddo
505 !
506 ! check whether a SPC suffices
507 !
508  if(dabs(1.d0-dmax).lt.1.d-3) then
509  val=0.d0
510  if(nam.gt.0) iamplitude=0
511  type='R'
512  call bounadd(irotnode,imax,imax,val,nodeboun,
513  & ndirboun,xboun,nboun,nboun_,iamboun,
514  & iamplitude,nam,ipompc,nodempc,coefmpc,
515  & nmpc,nmpc_,mpcfree,inotr,trab,ntrans,
516  & ikboun,ilboun,ikmpc,ilmpc,co,nk,nk_,labmpc,
517  & type,typeboun,nmethod,iperturb,fixed,vold,
518  & irotnode,mi,label)
519  else
520 !
521 ! check for an unused rotational DOF
522 !
523  isol=0
524  do l=1,3
525  idof=8*(node-1)+3+imax
526  call nident(ikboun,idof,nboun,id)
527  if((id.gt.0).and.(ikboun(id).eq.idof)) then
528  imax=imax+1
529  if(imax.gt.3) imax=imax-3
530  cycle
531  endif
532  isol=1
533  exit
534  enddo
535 !
536 ! if one of the rotational dofs was not used so far,
537 ! it can be taken as dependent side for fixing the
538 ! rotation about the normal. If all dofs were used,
539 ! no additional equation is needed.
540 !
541  if(isol.eq.1) then
542  idof=8*(irotnode-1)+imax
543  call nident(ikmpc,idof,nmpc,id)
544  nmpc=nmpc+1
545  if(nmpc.gt.nmpc_) then
546  write(*,*)
547  & '*ERROR in gen3dboun: increase nmpc_'
548  call exit(201)
549  endif
550 !
551  ipompc(nmpc)=mpcfree
552  labmpc(nmpc)=' '
553 !
554  do l=nmpc,id+2,-1
555  ikmpc(l)=ikmpc(l-1)
556  ilmpc(l)=ilmpc(l-1)
557  enddo
558  ikmpc(id+1)=idof
559  ilmpc(id+1)=nmpc
560 !
561  nodempc(1,mpcfree)=irotnode
562  nodempc(2,mpcfree)=imax
563  coefmpc(mpcfree)=xnoref(imax)
564  mpcfree=nodempc(3,mpcfree)
565  imax=imax+1
566  if(imax.gt.3) imax=imax-3
567  nodempc(1,mpcfree)=irotnode
568  nodempc(2,mpcfree)=imax
569  coefmpc(mpcfree)=xnoref(imax)
570  mpcfree=nodempc(3,mpcfree)
571  imax=imax+1
572  if(imax.gt.3) imax=imax-3
573  nodempc(1,mpcfree)=irotnode
574  nodempc(2,mpcfree)=imax
575  coefmpc(mpcfree)=xnoref(imax)
576  mpcfreeold=mpcfree
577  mpcfree=nodempc(3,mpcfree)
578  nodempc(3,mpcfreeold)=0
579  endif
580  endif
581  endif
582  cycle
583  endif
584 !
585 ! 2. no existing knot
586 ! rotational dof
587 ! no dynamics
588 ! => mean rotation MPC is created
589 !
590 ! all cases except nonlinear dynamic case: creation
591 ! of meanrotation MPC's
592 !
593  if((idir.gt.3).and.((nmethod.ne.4).or.(iperturb.le.1)))then
594 !
595 ! create a mean rotation MPC
596 ! advantage: more accurate since less constraining
597 ! disadvantage: cannot exceed 90 degrees rotation
598 !
599 ! if a mean rotation MPC has already been created
600 ! for idof, ilboun(id) contains the index of the
601 ! SPC created for the angle value
602 !
603  idof=8*(node-1)+idir
604  call nident(ikboun,idof,nboun,id)
605  if(ilboun(id).ne.i) cycle
606 
607  idirref=idir-3
608 !
609  if(lakon(ielem)(7:7).eq.'L') then
610  lstart=3
611  lend=1
612  linc=-2
613 !
614 ! calculating the normal vector =
615 ! vector along the drilling direction
616 !
617  do j=1,3
618  xn(j)=co(j,knor(indexk+3))-co(j,knor(indexk+1))
619  enddo
620  dd=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3))
621  do j=1,3
622  xn(j)=xn(j)/dd
623  enddo
624 !
625  elseif(lakon(ielem)(7:7).eq.'B') then
626  lstart=4
627  lend=1
628  linc=-1
629  endif
630 !
631 ! check for transformations
632 !
633  if(ntrans.le.0) then
634  itr=0
635  elseif(inotr(1,node).eq.0) then
636  itr=0
637  else
638  itr=inotr(1,node)
639  endif
640 !
641 ! determine a unit vector on the rotation axis
642 !
643  if(itr.eq.0) then
644  do j=1,3
645  do k=1,3
646  a(j,k)=0.d0
647  enddo
648  a(j,j)=1.d0
649  enddo
650  else
651  call transformatrix(trab(1,itr),co(1,node),a)
652  endif
653 !
654 ! check whether the rotation vector does not
655 ! have a component along the drilling direction
656 !
657  if(lakon(ielem)(7:7).eq.'L') then
658  dot=a(1,idirref)*xn(1)+a(2,idirref)*xn(2)+
659  & a(3,idirref)*xn(3)
660 !
661 ! rotation vectors with a substantial component
662 ! along the drilling direction are not allowed if
663 ! the rotation value (SPC) is nonzero; else the
664 ! rotation vector is projected on the tangential
665 ! plane
666 !
667  if(dabs(dot).gt.0.05) then
668  if(xboun(i).gt.1.d-10) then
669  write(*,*) '*ERROR in gen3dboun: rotation'
670  write(*,*) ' vector in node ',node
671  write(*,*) ' and direction ',idir-1
672  write(*,*) ' has a significant'
673  write(*,*)
674  & ' component along the drilling'
675  write(*,*) ' direction and a nonzero'
676  write(*,*) ' boundary value; this is not'
677  write(*,*) ' allowed'
678  call exit(201)
679  endif
680  endif
681 !
682 ! rotation vectors closer than 45 degrees with the normal
683 ! (= drilling direction) are not taken into account
684 !
685  if(dabs(dot).gt.0.70710678d0) cycle
686 !
687 ! projecting the rotation vector on the tangent plane
688 !
689  do k=1,3
690  a(k,idirref)=a(k,idirref)-dot*xn(k)
691  enddo
692  dd=0.d0
693  do k=1,3
694  dd=dd+a(k,idirref)**2
695  enddo
696  dd=dsqrt(dd)
697  do k=1,3
698  a(k,idirref)=a(k,idirref)/dd
699  enddo
700  endif
701 !
702 ! specific label for mean rotations for beams and
703 ! shells
704 !
705  label='MEANROTBS '
706  nnodes=0
707  do j=lstart,lend,linc
708  nodeact=knor(indexk+j)
709  do k=1,3
710  nnodes=nnodes+1
711  call usermpc(ipompc,nodempc,coefmpc,
712  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
713  & nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
714  & nboun,nboun_,nnodes,nodeact,co,label,
715  & typeboun,iperturb,node,idirref,xboun)
716  enddo
717  enddo
718 !
719 ! rotation value term
720 !
721  nodeact=nk+1
722  do k=1,3
723  co(k,nodeact)=a(k,idirref)
724  enddo
725  nnodes=nnodes+1
726  call usermpc(ipompc,nodempc,coefmpc,
727  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
728  & nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
729  & nboun,nboun_,nnodes,nodeact,co,label,
730  & typeboun,iperturb,node,idirref,xboun)
731 !
732 ! inhomogeneous term
733 !
734  nodeact=0
735  call usermpc(ipompc,nodempc,coefmpc,
736  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
737  & nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
738  & nboun,nboun_,nnodes,nodeact,co,label,
739  & typeboun,iperturb,node,idirref,xboun)
740 !
741 ! end meanrotationmpc
742 !
743 ! SPC angle term
744 !
745  if(nodeact.ne.-1) then
746  idir=1
747  type='B'
748  call bounadd(nk,idir,idir,val,nodeboun,
749  & ndirboun,xboun,nboun,nboun_,iamboun,
750  & iamplitude,nam,ipompc,nodempc,coefmpc,
751  & nmpc,nmpc_,mpcfree,inotr,trab,ntrans,
752  & ikboun,ilboun,ikmpc,ilmpc,co,nk,nk_,labmpc,
753  & type,typeboun,nmethod,iperturb,fixed,vold,
754  & nk,mi,label)
755 !
756 ! storing the index of the SPC with the angle
757 ! value in ilboun(id)
758 !
759  ilboun(id)=nboun
760  endif
761 !
762  cycle
763  endif
764 !
765 ! 3. no existing knot
766 ! translational dof
767 !
768 ! a) 2d shell element: generate MPC's
769 !
770 ! u(n_1)+u(n_3)=2*u(n)
771 !
772  if(lakon(ielem)(7:7).eq.'L') then
773  newnode=knor(indexk+1)
774  idof=8*(newnode-1)+idir
775  call nident(ikmpc,idof,nmpc,id)
776  if((id.le.0).or.(ikmpc(id).ne.idof)) then
777  nmpc=nmpc+1
778  if(nmpc.gt.nmpc_) then
779  write(*,*)
780  & '*ERROR in gen3dboun: increase nmpc_'
781  call exit(201)
782  endif
783  labmpc(nmpc)=' '
784  ipompc(nmpc)=mpcfree
785  do j=nmpc,id+2,-1
786  ikmpc(j)=ikmpc(j-1)
787  ilmpc(j)=ilmpc(j-1)
788  enddo
789  ikmpc(id+1)=idof
790  ilmpc(id+1)=nmpc
791  nodempc(1,mpcfree)=newnode
792  nodempc(2,mpcfree)=idir
793  coefmpc(mpcfree)=1.d0
794  mpcfree=nodempc(3,mpcfree)
795  if(mpcfree.eq.0) then
796  write(*,*)
797  & '*ERROR in gen3dboun: increase memmpc_'
798  call exit(201)
799  endif
800  nodempc(1,mpcfree)=knor(indexk+3)
801  nodempc(2,mpcfree)=idir
802  coefmpc(mpcfree)=1.d0
803  mpcfree=nodempc(3,mpcfree)
804  if(mpcfree.eq.0) then
805  write(*,*)
806  & '*ERROR in gen3dboun: increase memmpc_'
807  call exit(201)
808  endif
809  nodempc(1,mpcfree)=node
810  nodempc(2,mpcfree)=idir
811  coefmpc(mpcfree)=-2.d0
812  mpcfreenew=nodempc(3,mpcfree)
813  if(mpcfreenew.eq.0) then
814  write(*,*)
815  & '*ERROR in gen3dboun: increase memmpc_'
816  call exit(201)
817  endif
818  nodempc(3,mpcfree)=0
819  mpcfree=mpcfreenew
820  endif
821 !
822 ! u(n_2)=u(n)
823 !
824  newnode=knor(indexk+2)
825  idof=8*(newnode-1)+idir
826  call nident(ikmpc,idof,nmpc,id)
827  if((id.le.0).or.(ikmpc(id).ne.idof)) then
828  nmpc=nmpc+1
829  if(nmpc.gt.nmpc_) then
830  write(*,*)
831  & '*ERROR in gen3dboun: increase nmpc_'
832  call exit(201)
833  endif
834  labmpc(nmpc)=' '
835  ipompc(nmpc)=mpcfree
836  do j=nmpc,id+2,-1
837  ikmpc(j)=ikmpc(j-1)
838  ilmpc(j)=ilmpc(j-1)
839  enddo
840  ikmpc(id+1)=idof
841  ilmpc(id+1)=nmpc
842  nodempc(1,mpcfree)=newnode
843  nodempc(2,mpcfree)=idir
844  coefmpc(mpcfree)=1.d0
845  mpcfree=nodempc(3,mpcfree)
846  if(mpcfree.eq.0) then
847  write(*,*)
848  & '*ERROR in gen3dboun: increase memmpc_'
849  call exit(201)
850  endif
851  nodempc(1,mpcfree)=node
852  nodempc(2,mpcfree)=idir
853  coefmpc(mpcfree)=-1.d0
854  mpcfreenew=nodempc(3,mpcfree)
855  if(mpcfreenew.eq.0) then
856  write(*,*)
857  & '*ERROR in gen3dboun: increase memmpc_'
858  call exit(201)
859  endif
860  nodempc(3,mpcfree)=0
861  mpcfree=mpcfreenew
862  endif
863 !
864 ! fixing the temperature degrees of freedom
865 !
866  if(idir.eq.0) then
867 !
868 ! t(n_3)=t(n)
869 !
870  newnode=knor(indexk+3)
871  idof=8*(newnode-1)+idir
872  call nident(ikmpc,idof,nmpc,id)
873  if((id.le.0).or.(ikmpc(id).ne.idof)) then
874  nmpc=nmpc+1
875  if(nmpc.gt.nmpc_) then
876  write(*,*)
877  & '*ERROR in gen3dboun: increase nmpc_'
878  call exit(201)
879  endif
880  labmpc(nmpc)=' '
881  ipompc(nmpc)=mpcfree
882  do j=nmpc,id+2,-1
883  ikmpc(j)=ikmpc(j-1)
884  ilmpc(j)=ilmpc(j-1)
885  enddo
886  ikmpc(id+1)=idof
887  ilmpc(id+1)=nmpc
888  nodempc(1,mpcfree)=newnode
889  nodempc(2,mpcfree)=idir
890  coefmpc(mpcfree)=1.d0
891  mpcfree=nodempc(3,mpcfree)
892  if(mpcfree.eq.0) then
893  write(*,*)
894  & '*ERROR in gen3dboun: increase memmpc_'
895  call exit(201)
896  endif
897  nodempc(1,mpcfree)=node
898  nodempc(2,mpcfree)=idir
899  coefmpc(mpcfree)=-1.d0
900  mpcfreenew=nodempc(3,mpcfree)
901  if(mpcfreenew.eq.0) then
902  write(*,*)
903  & '*ERROR in gen3dboun: increase memmpc_'
904  call exit(201)
905  endif
906  nodempc(3,mpcfree)=0
907  mpcfree=mpcfreenew
908  endif
909  endif
910  elseif(lakon(ielem)(7:7).eq.'B') then
911 !
912 ! b) 1d beam element: generate MPC's
913 !
914 ! u(n_1)+u(n_2)+u(n_3)+u(n_4)=4*u(n)
915 !
916  nnode=4
917  xnode=4.d0
918  newnode=knor(indexk+1)
919  idof=8*(newnode-1)+idir
920  call nident(ikmpc,idof,nmpc,id)
921  if((id.le.0).or.(ikmpc(id).ne.idof)) then
922  nmpc=nmpc+1
923  if(nmpc.gt.nmpc_) then
924  write(*,*)
925  & '*ERROR in gen3dboun: increase nmpc_'
926  call exit(201)
927  endif
928  labmpc(nmpc)=' '
929  ipompc(nmpc)=mpcfree
930  do j=nmpc,id+2,-1
931  ikmpc(j)=ikmpc(j-1)
932  ilmpc(j)=ilmpc(j-1)
933  enddo
934  ikmpc(id+1)=idof
935  ilmpc(id+1)=nmpc
936  nodempc(1,mpcfree)=newnode
937  nodempc(2,mpcfree)=idir
938  coefmpc(mpcfree)=1.d0
939  mpcfree=nodempc(3,mpcfree)
940  if(mpcfree.eq.0) then
941  write(*,*)
942  & '*ERROR in gen3dboun: increase memmpc_'
943  call exit(201)
944  endif
945  do k=2,nnode
946  nodempc(1,mpcfree)=knor(indexk+k)
947  nodempc(2,mpcfree)=idir
948  coefmpc(mpcfree)=1.d0
949  mpcfree=nodempc(3,mpcfree)
950  if(mpcfree.eq.0) then
951  write(*,*)
952  & '*ERROR in gen3dboun: increase memmpc_'
953  call exit(201)
954  endif
955  enddo
956  nodempc(1,mpcfree)=node
957  nodempc(2,mpcfree)=idir
958  coefmpc(mpcfree)=-xnode
959  mpcfreenew=nodempc(3,mpcfree)
960  if(mpcfreenew.eq.0) then
961  write(*,*)
962  & '*ERROR in gen3dboun: increase memmpc_'
963  call exit(201)
964  endif
965  nodempc(3,mpcfree)=0
966  mpcfree=mpcfreenew
967  endif
968 !
969 ! fixing the temperature degrees of freedom
970 !
971  if(idir.eq.0) then
972  do k=2,4
973 !
974 ! t(n_k)=t(n), k=2,4
975 !
976  newnode=knor(indexk+k)
977  idof=8*(newnode-1)+idir
978  call nident(ikmpc,idof,nmpc,id)
979  if((id.le.0).or.(ikmpc(id).ne.idof)) then
980  nmpc=nmpc+1
981  if(nmpc.gt.nmpc_) then
982  write(*,*)
983  & '*ERROR in gen3dboun: increase nmpc_'
984  call exit(201)
985  endif
986  labmpc(nmpc)=' '
987  ipompc(nmpc)=mpcfree
988  do j=nmpc,id+2,-1
989  ikmpc(j)=ikmpc(j-1)
990  ilmpc(j)=ilmpc(j-1)
991  enddo
992  ikmpc(id+1)=idof
993  ilmpc(id+1)=nmpc
994  nodempc(1,mpcfree)=newnode
995  nodempc(2,mpcfree)=idir
996  coefmpc(mpcfree)=1.d0
997  mpcfree=nodempc(3,mpcfree)
998  if(mpcfree.eq.0) then
999  write(*,*)
1000  & '*ERROR in gen3dboun: increase memmpc_'
1001  call exit(201)
1002  endif
1003  nodempc(1,mpcfree)=node
1004  nodempc(2,mpcfree)=idir
1005  coefmpc(mpcfree)=-1.d0
1006  mpcfreenew=nodempc(3,mpcfree)
1007  if(mpcfreenew.eq.0) then
1008  write(*,*)
1009  & '*ERROR in gen3dboun: increase memmpc_'
1010  call exit(201)
1011  endif
1012  nodempc(3,mpcfree)=0
1013  mpcfree=mpcfreenew
1014  endif
1015  enddo
1016  endif
1017  else
1018 !
1019 ! c) 2d plane stress, plane strain or axisymmetric
1020 ! element: MPC in all but z-direction
1021 !
1022  newnode=knor(indexk+2)
1023  idof=8*(newnode-1)+idir
1024  call nident(ikmpc,idof,nmpc,id)
1025  if(((id.le.0).or.(ikmpc(id).ne.idof)).and.
1026  & (idir.ne.3)) then
1027  nmpc=nmpc+1
1028  if(nmpc.gt.nmpc_) then
1029  write(*,*)
1030  & '*ERROR in gen3dmpc: increase nmpc_'
1031  call exit(201)
1032  endif
1033  labmpc(nmpc)=' '
1034  ipompc(nmpc)=mpcfree
1035  do j=nmpc,id+2,-1
1036  ikmpc(j)=ikmpc(j-1)
1037  ilmpc(j)=ilmpc(j-1)
1038  enddo
1039  ikmpc(id+1)=idof
1040  ilmpc(id+1)=nmpc
1041  nodempc(1,mpcfree)=newnode
1042  nodempc(2,mpcfree)=idir
1043  coefmpc(mpcfree)=1.d0
1044  mpcfree=nodempc(3,mpcfree)
1045  if(mpcfree.eq.0) then
1046  write(*,*)
1047  & '*ERROR in gen3dmpc: increase memmpc_'
1048  call exit(201)
1049  endif
1050  nodempc(1,mpcfree)=node
1051  nodempc(2,mpcfree)=idir
1052  coefmpc(mpcfree)=-1.d0
1053  mpcfreenew=nodempc(3,mpcfree)
1054  if(mpcfreenew.eq.0) then
1055  write(*,*)
1056  & '*ERROR in gen3dmpc: increase memmpc_'
1057  call exit(201)
1058  endif
1059  nodempc(3,mpcfree)=0
1060  mpcfree=mpcfreenew
1061  endif
1062  endif
1063  endif
1064  enddo
1065 !
1066  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
rs
subroutine rs(nm, n, a, w, matz, z, fv1, fv2, ierr)
Definition: rs.f:27
c1
static double * c1
Definition: mafillvcompmain.c:30
dgesv
subroutine dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
Definition: dgesv.f:58
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
mpcrem
subroutine mpcrem(i, mpcfree, nodempc, nmpc, ikmpc, ilmpc, labmpc, coefmpc, ipompc)
Definition: mpcrem.f:21
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
usermpc
subroutine usermpc(ipompc, nodempc, coefmpc, labmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, nk, nk_, nodeboun, ndirboun, ikboun, ilboun, nboun, nboun_, nnodes, node, co, label, typeboun, iperturb, noderef, idirref, xboun)
Definition: usermpc.f:23
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
Hosted by OpenAircraft.com, (Michigan UAV, LLC)