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

Go to the source code of this file.

Functions/Subroutines

subroutine gen3dforc (ikforc, ilforc, nforc, nforc_, nodeforc, ndirforc, xforc, iamforc, ntrans, inotr, trab, rig, ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, labmpc, iponoel, inoel, iponoelmax, kon, ipkon, lakon, ne, iponor, xnor, knor, nam, nk, nk_, co, thicke, nodeboun, ndirboun, ikboun, ilboun, nboun, nboun_, iamboun, typeboun, xboun, nmethod, iperturb, istep, vold, mi, idefforc)
 

Function/Subroutine Documentation

◆ gen3dforc()

subroutine gen3dforc ( integer, dimension(*)  ikforc,
integer, dimension(*)  ilforc,
integer  nforc,
integer  nforc_,
integer, dimension(2,*)  nodeforc,
integer, dimension(*)  ndirforc,
real*8, dimension(*)  xforc,
integer, dimension(*)  iamforc,
integer  ntrans,
integer, dimension(2,*)  inotr,
real*8, dimension(7,*)  trab,
integer, dimension(*)  rig,
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(*)  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  nam,
integer  nk,
integer  nk_,
real*8, dimension(3,*)  co,
real*8, dimension(mi(3),*)  thicke,
integer, dimension(*)  nodeboun,
integer, dimension(*)  ndirboun,
integer, dimension(*)  ikboun,
integer, dimension(*)  ilboun,
integer  nboun,
integer  nboun_,
integer, dimension(*)  iamboun,
character*1, dimension(*)  typeboun,
real*8, dimension(*)  xboun,
integer  nmethod,
integer  iperturb,
integer  istep,
real*8, dimension(0:mi(2),*)  vold,
integer, dimension(*)  mi,
integer, dimension(*)  idefforc 
)
25 !
26 ! connects nodes of 1-D and 2-D elements, for which
27 ! concentrated forces were
28 ! defined, to the nodes of their expanded counterparts
29 !
30  implicit none
31 !
32  logical add,fixed,user,quadratic
33 !
34  character*1 type,typeboun(*)
35  character*8 lakon(*)
36  character*20 labmpc(*),label
37 !
38  integer ikforc(*),ilforc(*),nodeforc(2,*),ndirforc(*),itr,idirref,
39  & iamforc(*),idim,ier,matz,nodeact,linc,lend,lstart,nnodes,
40  & nforc,nforc_,ntrans,inotr(2,*),rig(*),ipompc(*),nodempc(3,*),
41  & nmpc,nmpc_,mpcfree,ikmpc(*),ilmpc(*),iponoel(*),inoel(3,*),
42  & iponoelmax,kon(*),ipkon(*),ne,iponor(2,*),knor(*),nforcold,
43  & i,node,index,ielem,j,indexe,indexk,nam,iamplitude,idir,
44  & irotnode,nk,nk_,newnode,idof,id,mpcfreenew,k,isector,ndepnodes,
45  & idepnodes(80),l,iexpnode,indexx,irefnode,imax,isol,mpcfreeold,
46  & nod,impc,istep,nodeboun(*),ndirboun(*),ikboun(*),ilboun(*),
47  & nboun,nboun_,iamboun(*),nmethod,iperturb,nrhs,ipiv(3),info,m,
48  & mi(*),idefforc(*),nedge
49 !
50  real*8 xforc(*),trab(7,*),coefmpc(*),xnor(*),val,co(3,*),dot,
51  & thicke(mi(3),*),pi,xboun(*),xnoref(3),dmax,d(3,3),e(3,3,3),
52  & alpha,q(3),w(3),xn(3),a(3,3),a1(3),a2(3),dd,c1,c2,c3,ww,c(3,3),
53  & vold(0:mi(2),*),e1(3),e2(3),t1(3),b(3,3),x(3),y(3),fv1(3),
54  & fv2(3),z(3,3),xi1,xi2,xi3,u(3,3),r(3,3)
55 !
56  data d /1.,0.,0.,0.,1.,0.,0.,0.,1./
57  data e /0.,0.,0.,0.,0.,-1.,0.,1.,0.,
58  & 0.,0.,1.,0.,0.,0.,-1.,0.,0.,
59  & 0.,-1.,0.,1.,0.,0.,0.,0.,0./
60 !
61  label=' '
62  fixed=.false.
63 !
64  add=.false.
65  user=.false.
66  pi=4.d0*datan(1.d0)
67  isector=0
68 !
69  nforcold=nforc
70  do i=1,nforcold
71  node=nodeforc(1,i)
72  if(node.gt.iponoelmax) then
73  if(ndirforc(i).gt.3) then
74  write(*,*) '*WARNING: in gen3dforc: node ',i,
75  & ' does not'
76  write(*,*) ' belong to a beam nor shell'
77  write(*,*) ' element and consequently has no'
78  write(*,*) ' rotational degrees of freedom'
79  endif
80  cycle
81  endif
82  index=iponoel(node)
83  if(index.eq.0) then
84  if(ndirforc(i).gt.3) then
85  write(*,*) '*WARNING: in gen3dforc: node ',i,
86  & ' does not'
87  write(*,*) ' belong to a beam nor shell'
88  write(*,*) ' element and consequently has no'
89  write(*,*) ' rotational degrees of freedom'
90  endif
91  cycle
92  endif
93  ielem=inoel(1,index)
94 !
95 ! checking whether element is linear or quardratic
96 !
97  if((lakon(ielem)(4:4).eq.'6').or.
98  & (lakon(ielem)(4:4).eq.'8')) then
99  quadratic=.false.
100  else
101  quadratic=.true.
102  endif
103 !
104 ! checking whether element is 3-sided or 4-sided
105 !
106  if((lakon(ielem)(4:4).eq.'6').or.
107  & (lakon(ielem)(4:5).eq.'15')) then
108  nedge=3
109  else
110  nedge=4
111  endif
112 !
113  j=inoel(2,index)
114  indexe=ipkon(ielem)
115  indexk=iponor(2,indexe+j)
116  if(nam.gt.0) iamplitude=iamforc(i)
117  idir=ndirforc(i)
118  val=xforc(i)
119 !
120  if(rig(node).ne.0) then
121  if(idir.gt.3) then
122  if(rig(node).lt.0) then
123  write(*,*) '*ERROR in gen3dforc: in node ',node
124  write(*,*) ' a rotational DOF is loaded;'
125  write(*,*) ' however, the elements to which'
126  write(*,*) ' this node belongs do not have'
127  write(*,*) ' rotational DOFs'
128  call exit(201)
129  endif
130 c val=xforc(i)
131  k=idir-3
132  irotnode=rig(node)
133  call forcadd(irotnode,k,val,nodeforc,
134  & ndirforc,xforc,nforc,nforc_,iamforc,
135  & iamplitude,nam,ntrans,trab,inotr,co,
136  & ikforc,ilforc,isector,add,user,idefforc,
137  & ipompc,nodempc,nmpc,ikmpc,ilmpc,labmpc)
138  endif
139  else
140 !
141 ! check for moments defined in any but the first step
142 ! nonlinear dynamic case: creation of knots
143 ! knots (expandable rigid bodies) can take rotational
144 ! values arbitrarily exceeding 90 degrees
145 !
146  if((idir.gt.3).and.((nmethod.eq.4).and.(iperturb.gt.1)))then
147 !
148 ! create a knot: determine the knot
149 !
150  ndepnodes=0
151  if(lakon(ielem)(7:7).eq.'L') then
152  do k=1,3
153  ndepnodes=ndepnodes+1
154  idepnodes(ndepnodes)=knor(indexk+k)
155  enddo
156  idim=1
157  elseif(lakon(ielem)(7:7).eq.'B') then
158  do k=1,8
159  ndepnodes=ndepnodes+1
160  idepnodes(ndepnodes)=knor(indexk+k)
161  enddo
162  idim=3
163  else
164  write(*,*)
165  & '*ERROR in gen3dboun: a rotational DOF was applied'
166  write(*,*)
167  & '* to node',node,' without rotational DOFs'
168  call exit(201)
169  endif
170 !
171 ! remove all MPC's in which the knot nodes are
172 ! dependent nodes
173 !
174  do k=1,ndepnodes
175  nod=idepnodes(k)
176  do l=1,3
177  idof=8*(nod-1)+l
178  call nident(ikmpc,idof,nmpc,id)
179  if(id.gt.0) then
180  if(ikmpc(id).eq.idof) then
181  impc=ilmpc(id)
182  call mpcrem(impc,mpcfree,nodempc,nmpc,
183  & ikmpc,ilmpc,labmpc,coefmpc,ipompc)
184  endif
185  endif
186  enddo
187  enddo
188 !
189 ! generate a rigid body knot
190 !
191  irefnode=node
192  nk=nk+1
193  if(nk.gt.nk_) then
194  write(*,*) '*ERROR in rigidbodies: increase nk_'
195  call exit(201)
196  endif
197  irotnode=nk
198  rig(node)=irotnode
199  nk=nk+1
200  if(nk.gt.nk_) then
201  write(*,*) '*ERROR in rigidbodies: increase nk_'
202  call exit(201)
203  endif
204  iexpnode=nk
205  do k=1,ndepnodes
206  call knotmpc(ipompc,nodempc,coefmpc,irefnode,
207  & irotnode,iexpnode,
208  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,
209  & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,
210  & idepnodes,typeboun,co,xboun,istep,k,ndepnodes,
211  & idim,e1,e2,t1)
212  enddo
213 !
214 ! determine the location of the center of gravity of
215 ! the section and its displacements
216 !
217  do l=1,3
218  q(l)=0.d0
219  w(l)=0.d0
220  enddo
221  if(ndepnodes.eq.3) then
222  do k=1,ndepnodes,2
223  nod=idepnodes(k)
224  do l=1,3
225  q(l)=q(l)+co(l,nod)
226  w(l)=w(l)+vold(l,nod)
227  enddo
228  enddo
229  do l=1,3
230  q(l)=q(l)/2.d0
231  w(l)=w(l)/2.d0
232  enddo
233  else
234  do k=1,ndepnodes
235  nod=idepnodes(k)
236  do l=1,3
237  q(l)=q(l)+co(l,nod)
238  w(l)=w(l)+vold(l,nod)
239  enddo
240  enddo
241  do l=1,3
242  q(l)=q(l)/ndepnodes
243  w(l)=w(l)/ndepnodes
244  enddo
245  endif
246 !
247 ! determine the uniform expansion
248 !
249  alpha=0.d0
250  do k=1,ndepnodes
251  nod=idepnodes(k)
252  dd=(co(1,nod)-q(1))**2
253  & +(co(2,nod)-q(2))**2
254  & +(co(3,nod)-q(3))**2
255  if(dd.lt.1.d-20) cycle
256  alpha=alpha+dsqrt(
257  & ((co(1,nod)+vold(1,nod)-q(1)-w(1))**2
258  & +(co(2,nod)+vold(2,nod)-q(2)-w(2))**2
259  & +(co(3,nod)+vold(3,nod)-q(3)-w(3))**2)/dd)
260  enddo
261  alpha=alpha/ndepnodes
262 !
263 ! determine the displacements of irotnodes
264 !
265  do l=1,3
266  do m=1,3
267  a(l,m)=0.d0
268  enddo
269  xn(l)=0.d0
270  enddo
271  do k=1,ndepnodes
272  nod=idepnodes(k)
273  dd=0.d0
274  do l=1,3
275  a1(l)=co(l,nod)-q(l)
276  a2(l)=vold(l,nod)-w(l)
277  dd=dd+a1(l)*a1(l)
278  enddo
279  dd=dsqrt(dd)
280  if(dd.lt.1.d-10) cycle
281  do l=1,3
282  a1(l)=a1(l)/dd
283  a2(l)=a2(l)/dd
284  enddo
285  xn(1)=xn(1)+(a1(2)*a2(3)-a1(3)*a2(2))
286  xn(2)=xn(2)+(a1(3)*a2(1)-a1(1)*a2(3))
287  xn(3)=xn(3)+(a1(1)*a2(2)-a1(2)*a2(1))
288  do l=1,3
289  do m=1,3
290  a(l,m)=a(l,m)+a1(l)*a1(m)
291  enddo
292  enddo
293  enddo
294 !
295  do l=1,3
296  do m=1,3
297  a(l,m)=a(l,m)/ndepnodes
298  enddo
299  xn(l)=xn(l)/ndepnodes
300  a(l,l)=1.d0-a(l,l)
301  enddo
302 !
303  m=3
304  nrhs=1
305  call dgesv(m,nrhs,a,m,ipiv,xn,m,info)
306  if(info.ne.0) then
307  write(*,*) '*ERROR in gen3dforc:'
308  write(*,*) ' singular system of equations'
309  call exit(201)
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 r
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 c vold(1,iexpnode)=alpha
354  vold(1,iexpnode)=alpha-1.d0
355 !
356 ! correction of the expansion values for beam sections
357 !
358  if(idim.eq.2) then
359 !
360 ! initializing matrices b and c
361 !
362  do l=1,3
363  do m=1,3
364  b(l,m)=0.d0
365  c(l,m)=0.d0
366  enddo
367  enddo
368 !
369 ! solving a least squares problem to determine
370 ! the transpose of the deformation gradient:
371 ! c.F^T=b
372 !
373  do k=1,ndepnodes
374  nod=idepnodes(k)
375  do l=1,3
376  x(l)=co(l,nod)-q(l)
377  y(l)=x(l)+vold(l,nod)-w(l)
378  enddo
379  do l=1,3
380  do m=1,3
381  c(l,m)=c(l,m)+x(l)*x(m)
382  b(l,m)=b(l,m)+x(l)*y(m)
383  enddo
384  enddo
385  enddo
386 !
387 ! solving the linear equation system
388 !
389  m=3
390  nrhs=3
391  call dgesv(m,nrhs,c,m,ipiv,b,m,info)
392  if(info.ne.0) then
393  write(*,*) '*ERROR in gen3dforc:'
394  write(*,*) ' singular system of equations'
395  call exit(201)
396  endif
397 !
398 ! now b=F^T
399 !
400 ! constructing the right stretch tensor
401 ! U=F^T.R
402 !
403  do l=1,3
404  do m=l,3
405  u(l,m)=b(l,1)*r(1,m)+b(l,2)*r(2,m)+
406  & b(l,3)*r(3,m)
407  enddo
408  enddo
409  u(2,1)=u(1,2)
410  u(3,1)=u(1,3)
411  u(3,2)=u(2,3)
412 !
413 ! determining the eigenvalues and eigenvectors of U
414 !
415  m=3
416  matz=1
417  ier=0
418  call rs(m,m,u,w,matz,z,fv1,fv2,ier)
419  if(ier.ne.0) then
420  write(*,*)
421  & '*ERROR in knotmpc while calculating the'
422  write(*,*) ' eigenvalues/eigenvectors'
423  call exit(201)
424  endif
425 !
426  if((dabs(w(1)-1.d0).lt.dabs(w(2)-1.d0)).and.
427  & (dabs(w(1)-1.d0).lt.dabs(w(3)-1.d0))) then
428  l=2
429  m=3
430  elseif((dabs(w(2)-1.d0).lt.dabs(w(1)-1.d0)).and.
431  & (dabs(w(2)-1.d0).lt.dabs(w(3)-1.d0))) then
432  l=1
433  m=3
434  else
435  l=1
436  m=2
437  endif
438  xi1=datan2((z(1,l)*e2(1)+z(2,l)*e2(2)+z(3,l)*e2(2)),
439  & (z(1,l)*e1(1)+z(2,l)*e1(2)+z(3,l)*e1(2)))
440  xi2=w(l)-1.d0
441  xi3=w(m)-1.d0
442 !
443  vold(1,iexpnode)=xi1
444  vold(2,iexpnode)=xi2
445  vold(3,iexpnode)=xi3
446  endif
447 !
448 ! apply the moment
449 !
450  idir=idir-3
451  val=xforc(i)
452  call forcadd(irotnode,idir,val,nodeforc,
453  & ndirforc,xforc,nforc,nforc_,iamforc,
454  & iamplitude,nam,ntrans,trab,inotr,co,
455  & ikforc,ilforc,isector,add,user,idefforc,
456  & ipompc,nodempc,nmpc,ikmpc,ilmpc,labmpc)
457 !
458 ! check for shells whether the rotation about the normal
459 ! on the shell has been eliminated
460 !
461  if(lakon(ielem)(7:7).eq.'L') then
462  indexx=iponor(1,indexe+j)
463  do k=1,3
464  xnoref(k)=xnor(indexx+k)
465  enddo
466  dmax=0.d0
467  imax=0
468  do k=1,3
469  if(dabs(xnoref(k)).gt.dmax) then
470  dmax=dabs(xnoref(k))
471  imax=k
472  endif
473  enddo
474 !
475 ! check whether a SPC suffices
476 !
477  if(dabs(1.d0-dmax).lt.1.d-3) then
478  val=0.d0
479  if(nam.gt.0) iamplitude=0
480  type='R'
481  call bounadd(irotnode,imax,imax,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  else
489 !
490 ! check for an unused rotational DOF
491 !
492  isol=0
493  do l=1,3
494  idof=8*(node-1)+3+imax
495  call nident(ikboun,idof,nboun,id)
496  if((id.gt.0).and.(ikboun(id).eq.idof)) then
497  imax=imax+1
498  if(imax.gt.3) imax=imax-3
499  cycle
500  endif
501  isol=1
502  exit
503  enddo
504 !
505 ! if one of the rotational dofs was not used so far,
506 ! it can be taken as dependent side for fixing the
507 ! rotation about the normal. If all dofs were used,
508 ! no additional equation is needed.
509 !
510  if(isol.eq.1) then
511  idof=8*(irotnode-1)+imax
512  call nident(ikmpc,idof,nmpc,id)
513  nmpc=nmpc+1
514  if(nmpc.gt.nmpc_) then
515  write(*,*)
516  & '*ERROR in gen3dnor: increase nmpc_'
517  call exit(201)
518  endif
519 !
520  ipompc(nmpc)=mpcfree
521  labmpc(nmpc)=' '
522 !
523  do l=nmpc,id+2,-1
524  ikmpc(l)=ikmpc(l-1)
525  ilmpc(l)=ilmpc(l-1)
526  enddo
527  ikmpc(id+1)=idof
528  ilmpc(id+1)=nmpc
529 !
530  nodempc(1,mpcfree)=irotnode
531  nodempc(2,mpcfree)=imax
532  coefmpc(mpcfree)=xnoref(imax)
533  mpcfree=nodempc(3,mpcfree)
534  imax=imax+1
535  if(imax.gt.3) imax=imax-3
536  nodempc(1,mpcfree)=irotnode
537  nodempc(2,mpcfree)=imax
538  coefmpc(mpcfree)=xnoref(imax)
539  mpcfree=nodempc(3,mpcfree)
540  imax=imax+1
541  if(imax.gt.3) imax=imax-3
542  nodempc(1,mpcfree)=irotnode
543  nodempc(2,mpcfree)=imax
544  coefmpc(mpcfree)=xnoref(imax)
545  mpcfreeold=mpcfree
546  mpcfree=nodempc(3,mpcfree)
547  nodempc(3,mpcfreeold)=0
548  endif
549  endif
550  endif
551  cycle
552  endif
553 !
554 ! all cases except nonlinear dynamic case: creation
555 ! of mean rotation MPC's
556 !
557  if((idir.gt.3).and.((nmethod.ne.4).or.(iperturb.le.1)))then
558 !
559 ! create a mean rotation MPC
560 ! advantage: more accurate since less constraining
561 ! disadvantage: cannot exceed 90 degrees rotation
562 !
563 ! if a mean rotation MPC has already been created
564 ! for idof, ilforc(id) contains the index of the
565 ! SPC created for the angle value
566 !
567  idof=8*(node-1)+idir
568  call nident(ikforc,idof,nforc,id)
569  if(ilforc(id).ne.i) cycle
570 
571  idirref=idir-3
572 !
573  if(lakon(ielem)(7:7).eq.'L') then
574  lstart=3
575  lend=1
576  linc=-2
577 !
578 ! calculating the normal vector =
579 ! vector along the drilling direction
580 !
581  do j=1,3
582  xn(j)=co(j,knor(indexk+3))-co(j,knor(indexk+1))
583  enddo
584  dd=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3))
585  do j=1,3
586  xn(j)=xn(j)/dd
587  enddo
588 !
589  elseif(lakon(ielem)(7:7).eq.'B') then
590  lstart=4
591  lend=1
592  linc=-1
593  endif
594 !
595 ! check for transformations
596 !
597  if(ntrans.le.0) then
598  itr=0
599  elseif(inotr(1,node).eq.0) then
600  itr=0
601  else
602  itr=inotr(1,node)
603  endif
604 !
605 ! determine a unit vector on the rotation axis
606 !
607  if(itr.eq.0) then
608  do j=1,3
609  do k=1,3
610  a(j,k)=0.d0
611  enddo
612  a(j,j)=1.d0
613  enddo
614  else
615  call transformatrix(trab(1,itr),co(1,node),a)
616  endif
617 !
618 ! check whether the rotation vector does not
619 ! have a component along the drilling direction
620 !
621  if(lakon(ielem)(7:7).eq.'L') then
622  dot=a(1,idirref)*xn(1)+a(2,idirref)*xn(2)+
623  & a(3,idirref)*xn(3)
624  if(dot.gt.0.05) then
625  write(*,*) '*ERROR in gen3dforc: applied'
626  write(*,*) ' moment in node ',node
627  write(*,*) ' and direction ',idir-1
628  write(*,*) ' has a significant'
629  write(*,*)
630  & ' component along the drilling'
631  write(*,*) ' direction; this is not'
632  write(*,*) ' allowed'
633  call exit(201)
634  endif
635  endif
636 !
637 ! specific label for mean rotations for beams and
638 ! shells
639 !
640  label='MEANROTBS '
641  nnodes=0
642  do j=lstart,lend,linc
643  nodeact=knor(indexk+j)
644  do k=1,3
645  nnodes=nnodes+1
646  call usermpc(ipompc,nodempc,coefmpc,
647  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
648  & nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
649  & nboun,nboun_,nnodes,nodeact,co,label,
650  & typeboun,iperturb,node,idirref,xboun)
651  enddo
652  enddo
653 !
654 ! rotation value term
655 !
656  nodeact=nk+1
657  do k=1,3
658  co(k,nodeact)=a(k,idirref)
659  enddo
660  nnodes=nnodes+1
661  call usermpc(ipompc,nodempc,coefmpc,
662  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
663  & nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
664  & nboun,nboun_,nnodes,nodeact,co,label,
665  & typeboun,iperturb,node,idirref,xboun)
666 !
667 ! inhomogeneous term
668 !
669  nodeact=0
670  call usermpc(ipompc,nodempc,coefmpc,
671  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
672  & nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
673  & nboun,nboun_,nnodes,nodeact,co,label,
674  & typeboun,iperturb,node,idirref,xboun)
675 !
676 ! end meanrotationmpc
677 !
678 ! SPC angle term
679 !
680  if(nodeact.ne.-1) then
681  idir=1
682  call forcadd(nk,idir,val,nodeforc,
683  & ndirforc,xforc,nforc,nforc_,iamforc,
684  & iamplitude,nam,ntrans,trab,inotr,co,
685  & ikforc,ilforc,isector,add,user,idefforc,
686  & ipompc,nodempc,nmpc,ikmpc,ilmpc,labmpc)
687 !
688 ! storing the index of the SPC with the angle
689 ! value in ilboun(id)
690 !
691  ilforc(id)=nforc
692  endif
693 !
694  cycle
695  endif
696 !
697 ! 2d shell element: generate MPC's
698 !
699  if(lakon(ielem)(7:7).eq.'L') then
700  newnode=knor(indexk+1)
701  idof=8*(newnode-1)+idir
702  call nident(ikmpc,idof,nmpc,id)
703  if((id.le.0).or.(ikmpc(id).ne.idof)) then
704  nmpc=nmpc+1
705  if(nmpc.gt.nmpc_) then
706  write(*,*)
707  & '*ERROR in gen3dforc: increase nmpc_'
708  call exit(201)
709  endif
710  labmpc(nmpc)=' '
711  ipompc(nmpc)=mpcfree
712  do k=nmpc,id+2,-1
713  ikmpc(k)=ikmpc(k-1)
714  ilmpc(k)=ilmpc(k-1)
715  enddo
716  ikmpc(id+1)=idof
717  ilmpc(id+1)=nmpc
718 !
719 ! for middle nodes: u_1+u_3-2*u_node=0
720 ! for end nodes: -u_1+4*u_2-u_3-2*u_node=0
721 !
722 ! u_1 corresponds to knor(indexk+1)....
723 !
724  nodempc(1,mpcfree)=newnode
725  nodempc(2,mpcfree)=idir
726  if((j.gt.nedge).or.(.not.quadratic)) then
727  coefmpc(mpcfree)=1.d0
728  else
729  coefmpc(mpcfree)=-1.d0
730  endif
731  mpcfree=nodempc(3,mpcfree)
732  if(mpcfree.eq.0) then
733  write(*,*)
734  & '*ERROR in gen3dforc: increase memmpc_'
735  call exit(201)
736  endif
737 !
738  if((j.le.nedge).and.(quadratic)) then
739  nodempc(1,mpcfree)=knor(indexk+2)
740  nodempc(2,mpcfree)=idir
741  coefmpc(mpcfree)=4.d0
742  mpcfree=nodempc(3,mpcfree)
743  if(mpcfree.eq.0) then
744  write(*,*)
745  & '*ERROR in gen3dforc: increase memmpc_'
746  call exit(201)
747  endif
748  endif
749 !
750  nodempc(1,mpcfree)=knor(indexk+3)
751  nodempc(2,mpcfree)=idir
752  if((j.gt.nedge).or.(.not.quadratic)) then
753  coefmpc(mpcfree)=1.d0
754  else
755  coefmpc(mpcfree)=-1.d0
756  endif
757  mpcfree=nodempc(3,mpcfree)
758  if(mpcfree.eq.0) then
759  write(*,*)
760  & '*ERROR in gen3dforc: increase memmpc_'
761  call exit(201)
762  endif
763 !
764  nodempc(1,mpcfree)=node
765  nodempc(2,mpcfree)=idir
766  coefmpc(mpcfree)=-2.d0
767  mpcfreenew=nodempc(3,mpcfree)
768  if(mpcfreenew.eq.0) then
769  write(*,*)
770  & '*ERROR in gen3dforc: increase memmpc_'
771  call exit(201)
772  endif
773 !
774  nodempc(3,mpcfree)=0
775  mpcfree=mpcfreenew
776  endif
777  elseif(lakon(ielem)(7:7).eq.'B') then
778 !
779 ! 1d beam element: generate MPC's
780 !
781  newnode=knor(indexk+1)
782  idof=8*(newnode-1)+idir
783  call nident(ikmpc,idof,nmpc,id)
784  if((id.le.0).or.(ikmpc(id).ne.idof)) then
785  nmpc=nmpc+1
786  if(nmpc.gt.nmpc_) then
787  write(*,*)
788  & '*ERROR in gen3dforc: increase nmpc_'
789  call exit(201)
790  endif
791  labmpc(nmpc)=' '
792  ipompc(nmpc)=mpcfree
793  do k=nmpc,id+2,-1
794  ikmpc(k)=ikmpc(k-1)
795  ilmpc(k)=ilmpc(k-1)
796  enddo
797  ikmpc(id+1)=idof
798  ilmpc(id+1)=nmpc
799  do k=1,4
800  nodempc(1,mpcfree)=knor(indexk+k)
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 gen3dforc: increase memmpc_'
807  call exit(201)
808  endif
809  enddo
810  nodempc(1,mpcfree)=node
811  nodempc(2,mpcfree)=idir
812  coefmpc(mpcfree)=-4.d0
813  mpcfreenew=nodempc(3,mpcfree)
814  if(mpcfreenew.eq.0) then
815  write(*,*)
816  & '*ERROR in gen3dforc: increase memmpc_'
817  call exit(201)
818  endif
819  nodempc(3,mpcfree)=0
820  mpcfree=mpcfreenew
821  endif
822  else
823 !
824 ! 2d plane stress, plane strain or axisymmetric
825 ! element: MPC in all but z-direction
826 !
827  newnode=knor(indexk+2)
828  idof=8*(newnode-1)+idir
829  call nident(ikmpc,idof,nmpc,id)
830  if(((id.le.0).or.(ikmpc(id).ne.idof)).and.
831  & (idir.ne.3)) then
832  nmpc=nmpc+1
833  if(nmpc.gt.nmpc_) then
834  write(*,*)
835  & '*ERROR in gen3dmpc: increase nmpc_'
836  call exit(201)
837  endif
838  labmpc(nmpc)=' '
839  ipompc(nmpc)=mpcfree
840  do j=nmpc,id+2,-1
841  ikmpc(j)=ikmpc(j-1)
842  ilmpc(j)=ilmpc(j-1)
843  enddo
844  ikmpc(id+1)=idof
845  ilmpc(id+1)=nmpc
846  nodempc(1,mpcfree)=newnode
847  nodempc(2,mpcfree)=idir
848  coefmpc(mpcfree)=1.d0
849  mpcfree=nodempc(3,mpcfree)
850  if(mpcfree.eq.0) then
851  write(*,*)
852  & '*ERROR in gen3dmpc: increase memmpc_'
853  call exit(201)
854  endif
855  nodempc(1,mpcfree)=node
856  nodempc(2,mpcfree)=idir
857  coefmpc(mpcfree)=-1.d0
858  mpcfreenew=nodempc(3,mpcfree)
859  if(mpcfreenew.eq.0) then
860  write(*,*)
861  & '*ERROR in gen3dmpc: increase memmpc_'
862  call exit(201)
863  endif
864  nodempc(3,mpcfree)=0
865  mpcfree=mpcfreenew
866  endif
867  endif
868  endif
869  enddo
870 !
871  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
forcadd
subroutine forcadd(node, i, val, nodeforc, ndirforc, xforc, nforc, nforc_, iamforc, iamplitude, nam, ntrans, trab, inotr, co, ikforc, ilforc, isector, add, user, idefforc, ipompc, nodempc, nmpc, ikmpc, ilmpc, labmpc)
Definition: forcadd.f:23
Hosted by OpenAircraft.com, (Michigan UAV, LLC)