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

Go to the source code of this file.

Functions/Subroutines

subroutine nonlinmpc (co, vold, ipompc, nodempc, coefmpc, labmpc, nmpc, ikboun, ilboun, nboun, xbounact, aux, iaux, maxlenmpc, ikmpc, ilmpc, icascade, kon, ipkon, lakon, ne, reltime, newstep, xboun, fmpc, iit, idiscon, ncont, trab, ntrans, ithermal, mi)
 

Function/Subroutine Documentation

◆ nonlinmpc()

subroutine nonlinmpc ( real*8, dimension(3,*)  co,
real*8, dimension(0:mi(2),*)  vold,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
character*20, dimension(*)  labmpc,
integer  nmpc,
integer, dimension(*)  ikboun,
integer, dimension(*)  ilboun,
integer  nboun,
real*8, dimension(*)  xbounact,
real*8, dimension(*)  aux,
integer, dimension(*)  iaux,
integer  maxlenmpc,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
integer  icascade,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer  ne,
real*8  reltime,
integer  newstep,
real*8, dimension(*)  xboun,
real*8, dimension(*)  fmpc,
integer  iit,
integer  idiscon,
integer  ncont,
real*8, dimension(7,*)  trab,
integer  ntrans,
integer  ithermal,
integer, dimension(*)  mi 
)
23 !
24 ! updates the coefficients in nonlinear MPC's
25 !
26  implicit none
27 !
28  character*8 lakon(*)
29  character*20 labmpc(*),label
30 !
31  integer ipompc(*),nodempc(3,*),irefnode,irotnode,idir,
32  & nmpc,index,ii,inode,node,id,ikboun(*),ilboun(*),nboun,
33  & i,j,k,idof,na,nb,nc,np,i1,i2,i3,iaux(*),maxlenmpc,n,
34  & l,m,lmax,mmax,ikmpc(*),ilmpc(*),icascade,neigh(7,8),
35  & mpc,kon(*),ipkon(*),indexe,ne,idofrem,idofins,nmpc0,nmpc01,
36  & newstep,iit,idiscon,ncont,iexpnode,indexexp,nmpcdif,ntrans,
37  & nodei,noded,lathyp(3,6),inum,ndir,number,ithermal,mi(*),
38  & newknot,indexexp1,indexexp2,indexexp3,idim,nodea,nodeb,
39  & jmax,idira
40 !
41  real*8 co(3,*),coefmpc(*),vold(0:mi(2),*),c(3,3),dc(3,3,3),ww,
42  & e(3,3,3),d(3,3),w(3),f(3,3),c1,c2,c3,c4,c5,c6,xbounact(*),
43  & xboun(*),fmpc(*),expan,dd,a11,a12,a13,a21,a22,a23,a31,a32,a33,
44  & b11,b12,b13,b21,b22,b23,b31,b32,b33,aux(*),const,e1(3),e2(3),
45  & ddmax,a(3,3),b(3,3),xj,xi,et,ze,xlag(3,20),xeul(3,20),t1(3),
46  & coloc(3,8),reltime,csab(7),trab(7,*),pd(3),pi(3),e11(3,3),
47  & ad(3,3),ai(3,3),e22(3,3),e12(3,3),ru(3,3),ru1(3,3),
48  & ru2(3,3),ru3(3,3),u1(3,3),u2(3,3),u3(3,3),dcu(3,3,3),u(3,3),
49  & xi1,xi2,xi3,dco,dsi,dco2,dsi2,cj,ck,cl,dmax
50 !
51  data d /1.,0.,0.,0.,1.,0.,0.,0.,1./
52  data e /0.,0.,0.,0.,0.,-1.,0.,1.,0.,
53  & 0.,0.,1.,0.,0.,0.,-1.,0.,0.,
54  & 0.,-1.,0.,1.,0.,0.,0.,0.,0./
55  data neigh /1,9,2,12,4,17,5,2,9,1,10,3,18,6,
56  & 3,11,4,10,2,19,7,4,11,3,12,1,20,8,
57  & 5,13,6,16,8,17,1,6,13,5,14,7,18,2,
58  & 7,15,8,14,6,19,3,8,15,7,16,5,20,4/
59  data coloc /-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.,-1.,1.,-1.,
60  & -1.,-1.,1.,1.,-1.,1.,1.,1.,1.,-1.,1.,1./
61 !
62 ! latin hypercube positions in a 3 x 3 matrix
63 !
64  data lathyp /1,2,3,1,3,2,2,1,3,2,3,1,3,1,2,3,2,1/
65 !
66  irotnode=0
67  irefnode=0
68  if((icascade.eq.1).and.(newstep.ne.1).and.(ncont.eq.0)) icascade=0
69 !
70  ii=0
71  loop: do
72  ii=ii+1
73  if(ii.gt.nmpc) exit
74  if(labmpc(ii)(1:5).eq.'RIGID') then
75 !
76  index=ipompc(ii)
77  inode=nodempc(1,index)
78  idir=nodempc(2,index)
79  coefmpc(index)=1.d0
80 !
81  index=nodempc(3,index)
82  irefnode=nodempc(1,index)
83  coefmpc(index)=-1.d0
84 !
85  index=nodempc(3,index)
86  node=nodempc(1,index)
87 !
88 ! check whether the rotational node is the same as in
89 ! the last rigid body MPC
90 !
91  if(node.ne.irotnode) then
92  irotnode=node
93  w(1)=vold(1,node)
94  w(2)=vold(2,node)
95  w(3)=vold(3,node)
96  ww=dsqrt(w(1)*w(1)+w(2)*w(2)+w(3)*w(3))
97 !
98  c1=dcos(ww)
99  if(ww.gt.1.d-10) then
100  c2=dsin(ww)/ww
101  else
102  c2=1.d0
103  endif
104  if(ww.gt.1.d-5) then
105  c3=(1.d0-c1)/ww**2
106  else
107  c3=0.5d0
108  endif
109 !
110 ! rotation matrix c
111 !
112  do i=1,3
113  do j=1,3
114  c(i,j)=c1*d(i,j)+
115  & c2*(e(i,1,j)*w(1)+e(i,2,j)*w(2)+e(i,3,j)*w(3))+
116  & c3*w(i)*w(j)
117  enddo
118  enddo
119 !
120  c4=-c2
121  if(ww.gt.0.00464159) then
122  c5=(ww*dcos(ww)-dsin(ww))/ww**3
123  else
124  c5=-1.d0/3.d0
125  endif
126  if(ww.gt.0.0031623) then
127  c6=(ww*dsin(ww)-2.d0+2.d0*dcos(ww))/ww**4
128  else
129  c6=-1.d0/12.d0
130  endif
131 !
132 ! derivative of the rotation matrix c with respect to
133 ! the rotation vector w
134 !
135  do i=1,3
136  do j=1,3
137  do k=1,3
138  dc(i,j,k)=c4*w(k)*d(i,j)+
139  & c5*w(k)*(e(i,1,j)*w(1)+
140  & e(i,2,j)*w(2)+e(i,3,j)*w(3))+
141  & c2*e(i,k,j)+
142  & c6*w(k)*w(i)*w(j)+
143  & c3*(d(i,k)*w(j)+d(j,k)*w(i))
144  enddo
145  enddo
146  enddo
147 !
148 ! dummy variable
149 !
150  do i=1,3
151  do j=1,3
152 c f(i,j)=c(i,j)-d(i,j)-dc(i,j,1)*w(1)-dc(i,j,2)*w(2)-
153 c & dc(i,j,3)*w(3)
154  f(i,j)=c(i,j)-d(i,j)
155  enddo
156  enddo
157  endif
158 !
159 ! determining the coefficients of the rotational degrees
160 ! of freedom
161 !
162  coefmpc(index)=dc(idir,1,1)*(co(1,irefnode)-co(1,inode))+
163  & dc(idir,2,1)*(co(2,irefnode)-co(2,inode))+
164  & dc(idir,3,1)*(co(3,irefnode)-co(3,inode))
165 !
166  index=nodempc(3,index)
167  coefmpc(index)=dc(idir,1,2)*(co(1,irefnode)-co(1,inode))+
168  & dc(idir,2,2)*(co(2,irefnode)-co(2,inode))+
169  & dc(idir,3,2)*(co(3,irefnode)-co(3,inode))
170 !
171  index=nodempc(3,index)
172  coefmpc(index)=dc(idir,1,3)*(co(1,irefnode)-co(1,inode))+
173  & dc(idir,2,3)*(co(2,irefnode)-co(2,inode))+
174  & dc(idir,3,3)*(co(3,irefnode)-co(3,inode))
175 !
176 ! determining the nonhomogeneous part
177 !
178  index=nodempc(3,index)
179  coefmpc(index)=1.d0
180 !
181 ! old value of the nonhomogeneous term must be zero
182 !
183  vold(nodempc(2,index),nodempc(1,index))=0.d0
184  idof=8*(nodempc(1,index)-1)+nodempc(2,index)
185  call nident(ikboun,idof,nboun,id)
186  xbounact(ilboun(id))=f(idir,1)*(co(1,irefnode)-co(1,inode))+
187  & f(idir,2)*(co(2,irefnode)-co(2,inode))+
188  & f(idir,3)*(co(3,irefnode)-co(3,inode))-
189  & vold(idir,irefnode)+vold(idir,inode)
190 !
191  elseif(labmpc(ii)(1:4).eq.'KNOT') then
192 !
193 ! dependent node
194 !
195  index=ipompc(ii)
196  inode=nodempc(1,index)
197  idir=nodempc(2,index)
198  coefmpc(index)=1.d0
199 !
200 ! translation node
201 !
202  index=nodempc(3,index)
203  node=nodempc(1,index)
204  coefmpc(index)=-1.d0
205 !
206 ! check whether knot is the same as in the previous MPC
207 !
208  if(node.ne.irefnode) then
209  newknot=1
210  irefnode=node
211  else
212  newknot=0
213  endif
214 !
215  read(labmpc(ii)(5:5),'(i1)') idim
216 !
217 ! expansion node
218 !
219  index=nodempc(3,index)
220  iexpnode=nodempc(1,index)
221 !
222  if((idim.eq.1).or.(idim.eq.3)) then
223 !
224 ! nodes of knot lie on a straight line (1 term in MPC)
225 !
226  indexexp=index
227  elseif(idim.eq.2) then
228 !
229 ! node of knot lie in a plane (3 terms in MPC)
230 !
231  indexexp1=index
232  index=nodempc(3,index)
233 !
234  indexexp2=index
235  index=nodempc(3,index)
236 !
237  indexexp3=index
238 !
239  endif
240 !
241  if(newknot.eq.1) then
242  if((idim.eq.1).or.(idim.eq.3)) then
243  expan=1.d0+vold(1,iexpnode)
244  elseif(idim.eq.2) then
245  xi1=vold(1,iexpnode)
246  xi2=1.d0+vold(2,iexpnode)
247  xi3=1.d0+vold(3,iexpnode)
248  dco=dcos(xi1)
249  dsi=dsin(xi1)
250  dco2=dcos(2.d0*xi1)
251  dsi2=dsin(2.d0*xi1)
252  dd=xi2**2-xi3**2
253  endif
254  endif
255 !
256 ! rotation node
257 !
258  index=nodempc(3,index)
259  irotnode=nodempc(1,index)
260 !
261  if(newknot.eq.1) then
262  w(1)=vold(1,irotnode)
263  w(2)=vold(2,irotnode)
264  w(3)=vold(3,irotnode)
265  ww=dsqrt(w(1)*w(1)+w(2)*w(2)+w(3)*w(3))
266 !
267  c1=dcos(ww)
268  if(ww.gt.1.d-10) then
269  c2=dsin(ww)/ww
270  else
271  c2=1.d0
272  endif
273  if(ww.gt.1.d-5) then
274  c3=(1.d0-c1)/ww**2
275  else
276  c3=0.5d0
277  endif
278 !
279 ! rotation matrix c
280 !
281  do i=1,3
282  do j=1,3
283  c(i,j)=c1*d(i,j)+
284  & c2*(e(i,1,j)*w(1)+e(i,2,j)*w(2)+e(i,3,j)*w(3))+
285  & c3*w(i)*w(j)
286  enddo
287  enddo
288 !
289  c4=-c2
290  if(ww.gt.0.00464159) then
291  c5=(ww*dcos(ww)-dsin(ww))/ww**3
292  else
293  c5=-1.d0/3.d0
294  endif
295  if(ww.gt.0.0031623) then
296  c6=(ww*dsin(ww)-2.d0+2.d0*dcos(ww))/ww**4
297  else
298  c6=-1.d0/12.d0
299  endif
300 !
301 ! derivative of the rotation matrix c with respect to
302 ! the rotation vector w
303 !
304  do i=1,3
305  do j=1,3
306  do k=1,3
307  dc(i,j,k)=c4*w(k)*d(i,j)+
308  & c5*w(k)*(e(i,1,j)*w(1)+
309  & e(i,2,j)*w(2)+e(i,3,j)*w(3))+
310  & c2*e(i,k,j)+
311  & c6*w(k)*w(i)*w(j)+
312  & c3*(d(i,k)*w(j)+d(j,k)*w(i))
313  enddo
314  enddo
315  enddo
316 !
317  if((idim.eq.1).or.(idim.eq.3)) then
318 !
319 ! dummy variable for constant term
320 !
321  do i=1,3
322  do j=1,3
323  f(i,j)=expan*c(i,j)-d(i,j)
324  enddo
325  enddo
326 !
327 ! derivative of the rotation matrix w.r.t the rotation
328 ! vector multiplied by the expansion coefficient
329 !
330  do i=1,3
331  do j=1,3
332  do k=1,3
333  dc(i,j,k)=dc(i,j,k)*expan
334  enddo
335  enddo
336  enddo
337  elseif(idim.eq.2) then
338 !
339 ! local unit vectors
340 !
341  do i=1,3
342  t1(i)=co(i,irotnode)
343  e1(i)=co(i,iexpnode)
344  enddo
345  e2(1)=t1(2)*e1(3)-t1(3)*e1(2)
346  e2(2)=t1(3)*e1(1)-t1(1)*e1(3)
347  e2(3)=t1(1)*e1(2)-t1(2)*e1(1)
348 !
349  do i=1,3
350  do j=1,3
351  e11(i,j)=e1(i)*e1(j)
352  e22(i,j)=e2(i)*e2(j)
353  e12(i,j)=e1(i)*e2(j)+e2(i)*e1(j)
354 !
355  u(i,j)=t1(i)*t1(j)
356  & +((xi2*dco)**2+(xi3*dsi)**2)*e11(i,j)
357  & +((xi2*dsi)**2+(xi3*dco)**2)*e22(i,j)
358  & +dd*dco*dsi*e12(i,j)
359  u1(i,j)=dd*(dco2*e12(i,j)
360  & -dsi2*(e11(i,j)-e22(i,j)))
361  u2(i,j)=2.d0*xi2*(dco*dco*e11(i,j)
362  & +dsi*dsi*e22(i,j))+xi2*dsi2*e12(i,j)
363  u3(i,j)=2.d0*xi3*(dsi*dsi*e11(i,j)
364  & +dco*dco*e22(i,j))-xi3*dsi2*e12(i,j)
365  enddo
366  enddo
367 !
368 ! calculating r.u, r.u1, r.u2 and r.u3
369 !
370  do i=1,3
371  do j=1,3
372  ru(i,j)=0.d0
373  ru1(i,j)=0.d0
374  ru2(i,j)=0.d0
375  ru3(i,j)=0.d0
376 !
377  do k=1,3
378  ru(i,j)=ru(i,j)+c(i,k)*u(k,j)
379  ru1(i,j)=ru1(i,j)+c(i,k)*u1(k,j)
380  ru2(i,j)=ru2(i,j)+c(i,k)*u2(k,j)
381  ru3(i,j)=ru3(i,j)+c(i,k)*u3(k,j)
382  enddo
383  f(i,j)=ru(i,j)-d(i,j)
384  enddo
385  enddo
386 !
387 ! calculating dc.u
388 !
389  do i=1,3
390  do j=1,3
391  do k=1,3
392  dcu(i,j,k)=0.d0
393  do l=1,3
394  dcu(i,j,k)=dcu(i,j,k)+dc(i,l,k)*u(l,j)
395  enddo
396  dc(i,j,k)=dcu(i,j,k)
397  enddo
398  enddo
399  enddo
400  endif
401 !
402  endif
403 !
404 ! determining the coefficients of the expansion degrees
405 ! of freedom
406 !
407  if((idim.eq.1).or.(idim.eq.3)) then
408  coefmpc(indexexp)=c(idir,1)*(co(1,irefnode)-co(1,inode))+
409  & c(idir,2)*(co(2,irefnode)-co(2,inode))+
410  & c(idir,3)*(co(3,irefnode)-co(3,inode))
411  elseif(idim.eq.2) then
412 !
413 ! if xi2=xi3 xi1 cannot be determined, since its coefficient
414 ! is always zero (cf. definition of u1)
415 !
416  if(dabs(xi2-xi3).lt.1.d-10) then
417  coefmpc(indexexp1)=0.d0
418 c if(nodempc(2,indexexp1).ne.2) then
419  if(icascade.lt.1) icascade=1
420  nodempc(2,indexexp1)=2
421 c endif
422  else
423  coefmpc(indexexp1)=ru1(idir,1)*
424  & (co(1,irefnode)-co(1,inode))+
425  & ru1(idir,2)*(co(2,irefnode)-co(2,inode))+
426  & ru1(idir,3)*(co(3,irefnode)-co(3,inode))
427 c if(nodempc(2,indexexp1).ne.1) then
428  if(icascade.lt.1) icascade=1
429  nodempc(2,indexexp1)=1
430 c endif
431  endif
432  coefmpc(indexexp2)=ru2(idir,1)*
433  & (co(1,irefnode)-co(1,inode))+
434  & ru2(idir,2)*(co(2,irefnode)-co(2,inode))+
435  & ru2(idir,3)*(co(3,irefnode)-co(3,inode))
436  coefmpc(indexexp3)=ru3(idir,1)*
437  & (co(1,irefnode)-co(1,inode))+
438  & ru3(idir,2)*(co(2,irefnode)-co(2,inode))+
439  & ru3(idir,3)*(co(3,irefnode)-co(3,inode))
440  endif
441 !
442 ! determining the coefficients of the rotational degrees
443 ! of freedom
444 !
445  coefmpc(index)=(dc(idir,1,1)*(co(1,irefnode)-co(1,inode))+
446  & dc(idir,2,1)*(co(2,irefnode)-co(2,inode))+
447  & dc(idir,3,1)*(co(3,irefnode)-co(3,inode)))
448 !
449  index=nodempc(3,index)
450  coefmpc(index)=(dc(idir,1,2)*(co(1,irefnode)-co(1,inode))+
451  & dc(idir,2,2)*(co(2,irefnode)-co(2,inode))+
452  & dc(idir,3,2)*(co(3,irefnode)-co(3,inode)))
453 !
454  index=nodempc(3,index)
455  coefmpc(index)=(dc(idir,1,3)*(co(1,irefnode)-co(1,inode))+
456  & dc(idir,2,3)*(co(2,irefnode)-co(2,inode))+
457  & dc(idir,3,3)*(co(3,irefnode)-co(3,inode)))
458 !
459 ! determining the nonhomogeneous part
460 !
461  index=nodempc(3,index)
462  coefmpc(index)=1.d0
463 !
464 ! old value of the nonhomogeneous term must be zero
465 !
466  vold(nodempc(2,index),nodempc(1,index))=0.d0
467  idof=8*(nodempc(1,index)-1)+nodempc(2,index)
468  call nident(ikboun,idof,nboun,id)
469  xbounact(ilboun(id))=f(idir,1)*(co(1,irefnode)-co(1,inode))+
470  & f(idir,2)*(co(2,irefnode)-co(2,inode))+
471  & f(idir,3)*(co(3,irefnode)-co(3,inode))-
472  & vold(idir,irefnode)+vold(idir,inode)
473 !
474  elseif(labmpc(ii)(1:8).eq.'STRAIGHT') then
475 !
476 ! determining nodes and directions involved in MPC
477 !
478  index=ipompc(ii)
479  np=nodempc(1,index)
480  j=nodempc(2,index)
481  index=nodempc(3,index)
482  i=nodempc(2,index)
483  index=nodempc(3,index)
484  na=nodempc(1,index)
485  index=nodempc(3,nodempc(3,index))
486  nb=nodempc(1,index)
487 !
488 ! determining the coefficients
489 !
490  index=ipompc(ii)
491  c2=co(i,na)+vold(i,na)-co(i,nb)-vold(i,nb)
492  if(dabs(c2).lt.1.d-5) then
493  write(*,*) '*WARNING in nonlinmpc: coefficient of'
494  write(*,*)
495  & ' dependent node in STRAIGHT MPC is zero'
496  idofrem=8*(np-1)+j
497 !
498 ! determining a new dependent term
499 !
500  ddmax=abs(c2)
501  l=i
502  m=j
503  do k=1,2
504  l=l+1
505  m=m+1
506  if(l.gt.3) l=l-3
507  if(m.gt.3) m=m-3
508  dd=dabs(co(l,na)+vold(l,na)-co(l,nb)-vold(l,nb))
509  if(dd.gt.ddmax) then
510  ddmax=dd
511  lmax=l
512  mmax=m
513  endif
514  enddo
515  i=lmax
516  j=mmax
517  idofins=8*(np-1)+j
518 !
519  call changedepterm(ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
520 !
521  index=ipompc(ii)
522  nodempc(2,index)=j
523  index=nodempc(3,index)
524  nodempc(2,index)=i
525  index=nodempc(3,index)
526  nodempc(2,index)=j
527  index=nodempc(3,index)
528  nodempc(2,index)=i
529  index=nodempc(3,index)
530  nodempc(2,index)=j
531  index=nodempc(3,index)
532  nodempc(2,index)=i
533  index=nodempc(3,index)
534  nodempc(2,index)=j
535  if(icascade.eq.0) icascade=1
536  c2=co(i,na)+vold(i,na)-co(i,nb)-vold(i,nb)
537  endif
538  coefmpc(index)=c2
539  index=nodempc(3,index)
540  c3=co(j,nb)+vold(j,nb)-co(j,na)-vold(j,na)
541  coefmpc(index)=c3
542  index=nodempc(3,index)
543  c5=co(i,nb)+vold(i,nb)-co(i,np)-vold(i,np)
544  coefmpc(index)=c5
545  index=nodempc(3,index)
546  c6=co(j,np)+vold(j,np)-co(j,nb)-vold(j,nb)
547  coefmpc(index)=c6
548  index=nodempc(3,index)
549  c4=co(i,np)+vold(i,np)-co(i,na)-vold(i,na)
550  coefmpc(index)=c4
551  index=nodempc(3,index)
552  c1=co(j,na)+vold(j,na)-co(j,np)-vold(j,np)
553  coefmpc(index)=c1
554  index=nodempc(3,index)
555 !
556 ! nonhomogeneous term
557 !
558  coefmpc(index)=1.d0
559 !
560 ! old value of the nonhomogeneous term must be zero
561 !
562  idof=8*(nodempc(1,index)-1)+nodempc(2,index)
563  call nident(ikboun,idof,nboun,id)
564  xbounact(ilboun(id))=-c1*c2+c3*c4
565  if(newstep.eq.1) xboun(ilboun(id))=xbounact(ilboun(id))
566  vold(nodempc(2,index),nodempc(1,index))=
567  & (1.d0-reltime)*xboun(ilboun(id))
568  elseif(labmpc(ii)(1:5).eq.'PLANE') then
569 !
570 ! determining nodes and directions involved in MPC
571 !
572  index=ipompc(ii)
573  np=nodempc(1,index)
574  i1=nodempc(2,index)
575  index=nodempc(3,index)
576  i2=nodempc(2,index)
577  index=nodempc(3,index)
578  i3=nodempc(2,index)
579  index=nodempc(3,index)
580  na=nodempc(1,index)
581  index=nodempc(3,nodempc(3,nodempc(3,index)))
582  nb=nodempc(1,index)
583  index=nodempc(3,nodempc(3,nodempc(3,index)))
584  nc=nodempc(1,index)
585 !
586 ! determining the coefficients
587 !
588  a11=co(i1,np)+vold(i1,np)-co(i1,nc)-vold(i1,nc)
589  a12=co(i2,np)+vold(i2,np)-co(i2,nc)-vold(i2,nc)
590  a13=co(i3,np)+vold(i3,np)-co(i3,nc)-vold(i3,nc)
591  a21=co(i1,na)+vold(i1,na)-co(i1,nc)-vold(i1,nc)
592  a22=co(i2,na)+vold(i2,na)-co(i2,nc)-vold(i2,nc)
593  a23=co(i3,na)+vold(i3,na)-co(i3,nc)-vold(i3,nc)
594  a31=co(i1,nb)+vold(i1,nb)-co(i1,nc)-vold(i1,nc)
595  a32=co(i2,nb)+vold(i2,nb)-co(i2,nc)-vold(i2,nc)
596  a33=co(i3,nb)+vold(i3,nb)-co(i3,nc)-vold(i3,nc)
597 !
598  b11=a22*a33-a23*a32
599  b12=a31*a23-a21*a33
600  b13=a21*a32-a31*a22
601  b21=a32*a13-a12*a33
602  b22=a11*a33-a31*a13
603  b23=a31*a12-a11*a32
604  b31=a12*a23-a22*a13
605  b32=a21*a13-a11*a23
606  b33=a11*a22-a12*a21
607 !
608  index=ipompc(ii)
609  if(dabs(b11).lt.1.d-5) then
610  write(*,*) '*WARNING in nonlinmpc: coefficient of'
611  write(*,*) ' dependent node in PLANE MPC is zero'
612 !
613  idofrem=8*(nodempc(1,index)-1)+i1
614 !
615  if(dabs(b12).gt.dabs(b13)) then
616  idofins=8*(nodempc(1,index)-1)+i2
617  call changedepterm
618  & (ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
619  coefmpc(index)=b12
620  nodempc(2,index)=i2
621  index=nodempc(3,index)
622  coefmpc(index)=b11
623  nodempc(2,index)=i1
624  index=nodempc(3,index)
625  coefmpc(index)=b13
626  index=nodempc(3,index)
627  coefmpc(index)=b22
628  nodempc(2,index)=i2
629  index=nodempc(3,index)
630  coefmpc(index)=b21
631  nodempc(2,index)=i1
632  index=nodempc(3,index)
633  coefmpc(index)=b23
634  index=nodempc(3,index)
635  coefmpc(index)=b32
636  nodempc(2,index)=i2
637  index=nodempc(3,index)
638  coefmpc(index)=b31
639  nodempc(2,index)=i1
640  index=nodempc(3,index)
641  coefmpc(index)=b33
642  index=nodempc(3,index)
643  coefmpc(index)=-b12-b22-b32
644  nodempc(2,index)=i2
645  index=nodempc(3,index)
646  coefmpc(index)=-b11-b21-b31
647  nodempc(2,index)=i1
648  index=nodempc(3,index)
649  coefmpc(index)=-b13-b23-b33
650  if(icascade.eq.0) icascade=1
651  else
652  idofins=8*(nodempc(1,index)-1)+i3
653  call changedepterm
654  & (ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
655  coefmpc(index)=b13
656  nodempc(2,index)=i3
657  index=nodempc(3,index)
658  coefmpc(index)=b12
659  index=nodempc(3,index)
660  coefmpc(index)=b11
661  nodempc(2,index)=i1
662  index=nodempc(3,index)
663  coefmpc(index)=b23
664  nodempc(2,index)=i3
665  index=nodempc(3,index)
666  coefmpc(index)=b22
667  index=nodempc(3,index)
668  coefmpc(index)=b21
669  nodempc(2,index)=i1
670  index=nodempc(3,index)
671  coefmpc(index)=b33
672  nodempc(2,index)=i3
673  index=nodempc(3,index)
674  coefmpc(index)=b32
675  index=nodempc(3,index)
676  coefmpc(index)=b31
677  nodempc(2,index)=i1
678  index=nodempc(3,index)
679  coefmpc(index)=-b13-b23-b33
680  nodempc(2,index)=i3
681  index=nodempc(3,index)
682  coefmpc(index)=-b12-b22-b32
683  index=nodempc(3,index)
684  coefmpc(index)=-b11-b21-b31
685  nodempc(2,index)=i1
686  if(icascade.eq.0) icascade=1
687  endif
688  else
689  coefmpc(index)=b11
690  index=nodempc(3,index)
691  coefmpc(index)=b12
692  index=nodempc(3,index)
693  coefmpc(index)=b13
694  index=nodempc(3,index)
695  coefmpc(index)=b21
696  index=nodempc(3,index)
697  coefmpc(index)=b22
698  index=nodempc(3,index)
699  coefmpc(index)=b23
700  index=nodempc(3,index)
701  coefmpc(index)=b31
702  index=nodempc(3,index)
703  coefmpc(index)=b32
704  index=nodempc(3,index)
705  coefmpc(index)=b33
706  index=nodempc(3,index)
707  coefmpc(index)=-b11-b21-b31
708  index=nodempc(3,index)
709  coefmpc(index)=-b12-b22-b32
710  index=nodempc(3,index)
711  coefmpc(index)=-b13-b23-b33
712  endif
713  index=nodempc(3,index)
714  coefmpc(index)=1.d0
715  idof=8*(nodempc(1,index)-1)+nodempc(2,index)
716 !
717 ! old value of the nonhomogeneous term must be zero
718 !
719  call nident(ikboun,idof,nboun,id)
720  xbounact(ilboun(id))=a11*b11+a12*b12+a13*b13
721  if(newstep.eq.1) xboun(ilboun(id))=xbounact(ilboun(id))
722  vold(nodempc(2,index),nodempc(1,index))=0.d0
723  elseif(labmpc(ii)(1:4).eq.'BEAM') then
724  index=ipompc(ii)
725  nodea=nodempc(1,index)
726  idira=nodempc(2,index)
727  nodeb=nodempc(1,nodempc(3,nodempc(3,nodempc(3,index))))
728 !
729 ! check for degree of freedom with largest coefficient in the MPC
730 !
731  dmax=dabs(co(1,nodea)+vold(1,nodea)
732  & -co(1,nodeb)-vold(1,nodeb))
733  jmax=1
734  do j=2,3
735  dd=dabs(co(j,nodea)+vold(j,nodea)
736  & -co(j,nodeb)-vold(j,nodeb))
737  if(dd.gt.dmax) then
738  dmax=dd
739  jmax=j
740  endif
741  enddo
742 !
743 ! change degree of freedom if necessary
744 !
745  if(jmax.ne.idira) then
746  write(*,*) '*WARNING in nonlinmpc: dependent dof',idira
747  write(*,*) ' in a BEAM MPC is not maximum; it'
748  write(*,*) ' will be changed to dof ',jmax
749  write(*,*) ' which corresponds to the maximum'
750  write(*,*) ' coefficient'
751  idofrem=8*(nodea-1)+idira
752  idofins=8*(nodea-1)+jmax
753  call changedepterm(ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
754 !
755  j=jmax
756  k=j+1
757  if(k.gt.3) k=1
758  l=k+1
759  if(l.gt.3) l=1
760 !
761  nodempc(2,index)=j
762  index=nodempc(3,index)
763  nodempc(2,index)=k
764  index=nodempc(3,index)
765  nodempc(2,index)=l
766  index=nodempc(3,index)
767  nodempc(2,index)=j
768  index=nodempc(3,index)
769  nodempc(2,index)=k
770  index=nodempc(3,index)
771  nodempc(2,index)=l
772  index=nodempc(3,index)
773  index=ipompc(ii)
774  endif
775 !
776  j=nodempc(2,index)
777  cj=2.d0*(co(j,nodea)+vold(j,nodea)
778  & -co(j,nodeb)-vold(j,nodeb))
779  coefmpc(index)=cj
780  index=nodempc(3,index)
781 !
782  k=nodempc(2,index)
783  ck=2.d0*(co(k,nodea)+vold(k,nodea)
784  & -co(k,nodeb)-vold(k,nodeb))
785  coefmpc(index)=ck
786  index=nodempc(3,index)
787 !
788  l=nodempc(2,index)
789  cl=2.d0*(co(l,nodea)+vold(l,nodea)
790  & -co(l,nodeb)-vold(l,nodeb))
791  coefmpc(index)=cl
792  index=nodempc(3,index)
793 !
794  coefmpc(index)=-cj
795  index=nodempc(3,index)
796 !
797  coefmpc(index)=-ck
798  index=nodempc(3,index)
799 !
800  coefmpc(index)=-cl
801  index=nodempc(3,index)
802 !
803  coefmpc(index)=1.d0
804 !
805  idof=8*(nodempc(1,index)-1)+nodempc(2,index)
806  call nident(ikboun,idof,nboun,id)
807  xbounact(ilboun(id))=(cj**2+ck**2+cl**2)/4.d0
808  & -(co(1,nodea)-co(1,nodeb))**2
809  & -(co(2,nodea)-co(2,nodeb))**2
810  & -(co(3,nodea)-co(3,nodeb))**2
811 !
812  if(newstep.eq.1) xboun(ilboun(id))=xbounact(ilboun(id))
813  vold(nodempc(2,index),nodempc(1,index))=0.d0
814  elseif((labmpc(ii)(1:20).ne.' ').and.
815  & (labmpc(ii)(1:10).ne.'PRETENSION').and.
816  & (labmpc(ii)(1:11).ne.'THERMALPRET').and.
817 c & (labmpc(ii)(1:7).ne.'CONTACT').and.
818  & (labmpc(ii)(1:7).ne.'NETWORK').and.
819  & (labmpc(ii)(1:5).ne.'FLUID').and.
820  & (labmpc(ii)(1:6).ne.'CYCLIC').and.
821  & (labmpc(ii)(1:14).ne.'ROTTRACOUPLING').and.
822  & (labmpc(ii)(1:9).ne.'SUBCYCLIC')) then
823  index=ipompc(ii)
824  i=0
825  do
826  if(index.eq.0) exit
827  node=nodempc(1,index)
828  i=i+1
829  iaux(i)=nodempc(2,index)
830  iaux(maxlenmpc+i)=node
831  aux(6*maxlenmpc+i)=coefmpc(index)
832  do j=1,3
833  aux(3*(i-1)+j)=co(j,node)
834  aux(3*(maxlenmpc+i-1)+j)=vold(j,node)
835  enddo
836  index=nodempc(3,index)
837  enddo
838  n=i-1
839  if((labmpc(ii)(1:7).eq.'MEANROT').or.
840  & (labmpc(ii)(1:1).eq.'1')) then
841  call umpc_mean_rot(aux,aux(3*maxlenmpc+1),const,
842  & aux(6*maxlenmpc+1),iaux,n,fmpc(ii),iit,idiscon,
843  & iaux(maxlenmpc+1),ikmpc,nmpc,ikboun,nboun,
844  & labmpc(ii))
845  elseif(labmpc(ii)(1:4).eq.'DIST') then
846  call umpc_dist(aux,aux(3*maxlenmpc+1),const,
847  & aux(6*maxlenmpc+1),iaux,n,fmpc(ii),iit,idiscon)
848  elseif(labmpc(ii)(1:4).eq.'USER') then
849  call umpc_user(aux,aux(3*maxlenmpc+1),const,
850  & aux(6*maxlenmpc+1),iaux,n,fmpc(ii),iit,idiscon)
851  else
852  write(*,*) '*ERROR in nonlinmpc: mpc of type ',labmpc(ii)
853  write(*,*) ' is unknown'
854  call exit(201)
855  endif
856  index=ipompc(ii)
857 !
858  if((iaux(1).ne.nodempc(2,index)).or.
859  & (iaux(maxlenmpc+1).ne.nodempc(1,index))) then
860 !
861 ! dependent MPC has changed
862 !
863  idofrem=8*(nodempc(1,index)-1)+nodempc(2,index)
864  idofins=8*(nodempc(1,index)-1)+iaux(1)
865  call changedepterm(ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
866  if(icascade.eq.0) icascade=1
867  endif
868 !
869  i=0
870  do
871  if(index.eq.0) exit
872  i=i+1
873  if(i.le.n) then
874 !
875 ! check whether any directions have changed:
876 ! necessitates calling of remastruct
877 !
878  if(iaux(i).ne.nodempc(2,index)) then
879  if(icascade.eq.0) icascade=1
880  endif
881  nodempc(2,index)=iaux(i)
882  coefmpc(index)=aux(6*maxlenmpc+i)
883  else
884  coefmpc(index)=1.d0
885 !
886 ! old value of the nonhomogeneous term must be zero
887 !
888  vold(nodempc(2,index),nodempc(1,index))=0.d0
889  idof=8*(nodempc(1,index)-1)+nodempc(2,index)
890  call nident(ikboun,idof,nboun,id)
891  xbounact(ilboun(id))=const
892  endif
893  index=nodempc(3,index)
894  enddo
895  endif
896  enddo loop
897 !
898  return
umpc_user
subroutine umpc_user(x, u, f, a, jdof, n, force, iit, idiscon)
Definition: umpc_user.f:20
changedepterm
subroutine changedepterm(ikmpc, ilmpc, nmpc, mpc, idofrem, idofins)
Definition: changedepterm.f:20
umpc_mean_rot
subroutine umpc_mean_rot(x, u, f, a, jdof, n, force, iit, idiscon, jnode, ikmpc, nmpc, ikboun, nboun, label)
Definition: umpc_mean_rot.f:21
a31
static double * a31
Definition: mafillkmain.c:30
c1
static double * c1
Definition: mafillvcompmain.c:30
a21
static double * a21
Definition: mafillkmain.c:30
umpc_dist
subroutine umpc_dist(x, u, f, a, jdof, n, force, iit, idiscon)
Definition: umpc_dist.f:20
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
a11
static double * a11
Definition: mafillkmain.c:30
e11
static double * e11
Definition: radflowload.c:42
Hosted by OpenAircraft.com, (Michigan UAV, LLC)