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

Go to the source code of this file.

Functions/Subroutines

subroutine rhs (co, nk, kon, ipkon, lakon, ne, ipompc, nodempc, coefmpc, nmpc, nodeforc, ndirforc, xforc, nforc, nelemload, sideload, xload, nload, xbody, ipobody, nbody, cgr, fext, nactdof, neq, nmethod, ikmpc, ilmpc, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, iprestr, vold, iperturb, iexpl, plicon, nplicon, plkcon, nplkcon, npmat_, ttime, time, istep, iinc, dtime, physcon, ibody, xloadold, reltime, veold, matname, mi, ikactmech, nactmech, ielprop, prop, sti, xstateini, xstate, nstate_)
 

Function/Subroutine Documentation

◆ rhs()

subroutine rhs ( real*8, dimension(3,*)  co,
integer  nk,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer  ne,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer, dimension(2,*)  nodeforc,
integer, dimension(*)  ndirforc,
real*8, dimension(*)  xforc,
integer  nforc,
integer, dimension(2,*)  nelemload,
character*20, dimension(*)  sideload,
real*8, dimension(2,*)  xload,
integer  nload,
real*8, dimension(7,*)  xbody,
integer, dimension(2,*)  ipobody,
integer  nbody,
real*8, dimension(4,*)  cgr,
real*8, dimension(*)  fext,
integer, dimension(0:mi(2),*)  nactdof,
integer  neq,
integer  nmethod,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
real*8, dimension(0:21,ntmat_,*)  elcon,
integer, dimension(2,*)  nelcon,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
real*8, dimension(0:6,ntmat_,*)  alcon,
integer, dimension(2,*)  nalcon,
real*8, dimension(*)  alzero,
integer, dimension(mi(3),*)  ielmat,
integer, dimension(mi(3),*)  ielorien,
integer  norien,
real*8, dimension(7,*)  orab,
integer  ntmat_,
real*8, dimension(*)  t0,
real*8, dimension(*)  t1,
integer  ithermal,
integer  iprestr,
real*8, dimension(0:mi(2),*)  vold,
integer  iperturb,
integer  iexpl,
real*8, dimension(0:2*npmat_,ntmat_,*)  plicon,
integer, dimension(0:ntmat_,*)  nplicon,
real*8, dimension(0:2*npmat_,ntmat_,*)  plkcon,
integer, dimension(0:ntmat_,*)  nplkcon,
integer  npmat_,
real*8  ttime,
real*8  time,
integer  istep,
integer  iinc,
real*8  dtime,
real*8, dimension(*)  physcon,
integer, dimension(3,*)  ibody,
real*8, dimension(2,*)  xloadold,
real*8  reltime,
real*8, dimension(0:mi(2),*)  veold,
character*80, dimension(*)  matname,
integer, dimension(*)  mi,
integer, dimension(*)  ikactmech,
integer  nactmech,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop,
real*8, dimension(6,mi(1),*)  sti,
real*8, dimension(nstate_,mi(1),*)  xstateini,
real*8, dimension(nstate_,mi(1),*)  xstate,
integer  nstate_ 
)
29 !
30 ! filling the right hand side load vector b
31 !
32 ! b contains the contributions due to mechanical forces only
33 !
34  implicit none
35 !
36  character*8 lakon(*)
37  character*20 sideload(*)
38  character*80 matname(*)
39 !
40  integer kon(*),ipompc(*),nodempc(3,*),ipobody(2,*),nbody,
41  & nodeforc(2,*),ndirforc(*),nelemload(2,*),ikmpc(*),mi(*),
42  & ilmpc(*),nactdof(0:mi(2),*),konl(20),nelcon(2,*),ibody(3,*),
43  & nrhcon(*),nalcon(2,*),ielmat(mi(3),*),ielorien(mi(3),*),
44  & ipkon(*),ielprop(*),nstate_,
45  & nk,ne,nmpc,nforc,nload,neq,nmethod,nom,m,idm,
46  & ithermal,iprestr,iperturb,i,j,k,idist,jj,
47  & id,ist,index,jdof1,jdof,node1,ntmat_,indexe,nope,norien,
48  & iexpl,idof1,iinc,istep,icalccg,nplicon(0:ntmat_,*),
49  & nplkcon(0:ntmat_,*),npmat_,ikactmech(*),nactmech
50 !
51  real*8 co(3,*),coefmpc(*),xforc(*),xload(2,*),p1(3,2),
52  & p2(3,2),fext(*),bodyf(3),elcon(0:21,ntmat_,*),
53  & rhcon(0:1,ntmat_,*),xloadold(2,*),reltime,prop(*),
54  & alcon(0:6,ntmat_,*),alzero(*),orab(7,*),xbody(7,*),cgr(4,*),
55  & t0(*),t1(*),vold(0:mi(2),*),ff(60),time,ttime,dtime,
56  & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
57  & om(2),physcon(*),veold(0:mi(2),*),sti(6,mi(1),*),
58  & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
59 !
60  icalccg=0
61 !
62  if((nmethod.ge.4).and.(iperturb.lt.2).and.(nactmech.lt.neq/2))then
63 !
64 ! modal dynamics and steady state dynamics:
65 ! only nonzeros are reset to zero
66 !
67  do i=1,nactmech
68  fext(ikactmech(i)+1)=0.d0
69  enddo
70  else
71  do i=1,neq
72  fext(i)=0.d0
73  enddo
74  endif
75  nactmech=0
76 !
77 ! distributed forces (body forces or thermal loads or
78 ! residual stresses or distributed face loads)
79 !
80 c if((nbody.ne.0).or.(ithermal.ne.0).or.
81 c & (iprestr.ne.0).or.(nload.ne.0)) then
82  if((nbody.ne.0).or.(nload.ne.0)) then
83  idist=1
84  else
85  idist=0
86  endif
87 !
88 c if((ithermal.le.1).or.(ithermal.eq.3)) then
89  if(((ithermal.le.1).or.(ithermal.eq.3)).and.(idist.ne.0)) then
90 !
91 ! mechanical analysis: loop over all elements
92 !
93  do i=1,ne
94 !
95  if(ipkon(i).lt.0) cycle
96  indexe=ipkon(i)
97  if(lakon(i)(4:4).eq.'2') then
98  nope=20
99  elseif(lakon(i)(4:4).eq.'8') then
100  nope=8
101  elseif(lakon(i)(4:5).eq.'10') then
102  nope=10
103  elseif(lakon(i)(4:4).eq.'4') then
104  nope=4
105  elseif(lakon(i)(4:5).eq.'15') then
106  nope=15
107  elseif(lakon(i)(4:4).eq.'6') then
108  nope=6
109  else
110  cycle
111  endif
112 !
113  do j=1,nope
114  konl(j)=kon(indexe+j)
115  enddo
116 !
117 ! assigning centrifugal forces
118 !
119  if(nbody.gt.0) then
120  nom=0
121  om(1)=0.d0
122  om(2)=0.d0
123  bodyf(1)=0.d0
124  bodyf(2)=0.d0
125  bodyf(3)=0.d0
126 !
127  index=i
128  do
129  j=ipobody(1,index)
130  if(j.eq.0) exit
131  if(ibody(1,j).eq.1) then
132  nom=nom+1
133  if(nom.gt.2) then
134  write(*,*)'*ERROR in rhs: no more than two centri-'
135  write(*,*)' fugal loading cards allowed'
136  call exit(201)
137  endif
138  om(nom)=xbody(1,j)
139  p1(1,nom)=xbody(2,j)
140  p1(2,nom)=xbody(3,j)
141  p1(3,nom)=xbody(4,j)
142  p2(1,nom)=xbody(5,j)
143  p2(2,nom)=xbody(6,j)
144  p2(3,nom)=xbody(7,j)
145 !
146 ! assigning gravity forces
147 !
148  elseif(ibody(1,j).eq.2) then
149  bodyf(1)=bodyf(1)+xbody(1,j)*xbody(2,j)
150  bodyf(2)=bodyf(2)+xbody(1,j)*xbody(3,j)
151  bodyf(3)=bodyf(3)+xbody(1,j)*xbody(4,j)
152 !
153 ! assigning newton gravity forces
154 !
155  elseif(ibody(1,j).eq.3) then
156  call newton(icalccg,ne,ipkon,lakon,kon,t0,co,rhcon,
157  & nrhcon,ntmat_,physcon,i,cgr,bodyf,ielmat,ithermal,
158  & vold)
159  endif
160  index=ipobody(2,index)
161  if(index.eq.0) exit
162  enddo
163  endif
164 !
165  if(idist.ne.0)
166  & call e_c3d_rhs(co,nk,konl,lakon(i),p1,p2,om,bodyf,nbody,
167  & ff,i,nmethod,rhcon,ielmat,ntmat_,vold,iperturb,
168  & nelemload,sideload,xload,nload,idist,ttime,time,istep,
169  & iinc,dtime,xloadold,reltime,ipompc,nodempc,coefmpc,nmpc,
170  & ikmpc,ilmpc,veold,matname,mi,ielprop,prop)
171 !
172 ! modal dynamics and steady state dynamics:
173 ! location of nonzeros is stored
174 !
175  if((nmethod.ge.4).and.(iperturb.lt.2)) then
176  do jj=1,3*nope
177 !
178  j=(jj-1)/3+1
179  k=jj-3*(j-1)
180 !
181  node1=kon(indexe+j)
182  jdof1=nactdof(k,node1)
183 !
184 ! distributed forces
185 !
186  if(idist.ne.0) then
187  if(dabs(ff(jj)).lt.1.d-30) cycle
188  if(jdof1.le.0) then
189  if(nmpc.ne.0) then
190  idof1=(node1-1)*8+k
191  call nident(ikmpc,idof1,nmpc,id)
192  if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
193  id=ilmpc(id)
194  ist=ipompc(id)
195  index=nodempc(3,ist)
196  do
197  jdof1=nactdof(nodempc(2,index),
198  & nodempc(1,index))
199  if(jdof1.gt.0) then
200  fext(jdof1)=fext(jdof1)
201  & -coefmpc(index)*ff(jj)/coefmpc(ist)
202  call nident(ikactmech,jdof1-1,nactmech,
203  & idm)
204  do
205  if(idm.gt.0) then
206  if(ikactmech(idm).eq.jdof1-1) exit
207  endif
208  nactmech=nactmech+1
209  do m=nactmech,idm+2,-1
210  ikactmech(m)=ikactmech(m-1)
211  enddo
212  ikactmech(idm+1)=jdof1-1
213  exit
214  enddo
215  endif
216  index=nodempc(3,index)
217  if(index.eq.0) exit
218  enddo
219  endif
220  endif
221  cycle
222  endif
223  fext(jdof1)=fext(jdof1)+ff(jj)
224  call nident(ikactmech,jdof1-1,nactmech,
225  & idm)
226  do
227  if(idm.gt.0) then
228  if(ikactmech(idm).eq.jdof1-1) exit
229  endif
230  nactmech=nactmech+1
231  do m=nactmech,idm+2,-1
232  ikactmech(m)=ikactmech(m-1)
233  enddo
234  ikactmech(idm+1)=jdof1-1
235  exit
236  enddo
237  endif
238 !
239  enddo
240 !
241 ! other procedures
242 !
243  else
244  do jj=1,3*nope
245 !
246  j=(jj-1)/3+1
247  k=jj-3*(j-1)
248 !
249  node1=kon(indexe+j)
250  jdof1=nactdof(k,node1)
251 !
252 ! distributed forces
253 !
254  if(idist.ne.0) then
255  if(jdof1.le.0) then
256  if(nmpc.ne.0) then
257  idof1=(node1-1)*8+k
258  call nident(ikmpc,idof1,nmpc,id)
259  if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
260  id=ilmpc(id)
261  ist=ipompc(id)
262  index=nodempc(3,ist)
263  do
264  jdof1=nactdof(nodempc(2,index),
265  & nodempc(1,index))
266  if(jdof1.gt.0) then
267  fext(jdof1)=fext(jdof1)
268  & -coefmpc(index)*ff(jj)/coefmpc(ist)
269  endif
270  index=nodempc(3,index)
271  if(index.eq.0) exit
272  enddo
273  endif
274  endif
275  cycle
276  endif
277  fext(jdof1)=fext(jdof1)+ff(jj)
278  endif
279 !
280  enddo
281  endif
282  enddo
283 !
284 c else
285  elseif((ithermal.eq.2).and.(nload.gt.0)) then
286 !
287 ! thermal analysis: loop over all elements
288 !
289  do i=1,ne
290 !
291  if(ipkon(i).lt.0) cycle
292  indexe=ipkon(i)
293  if(lakon(i)(4:4).eq.'2') then
294  nope=20
295  elseif(lakon(i)(4:4).eq.'8') then
296  nope=8
297  elseif(lakon(i)(4:5).eq.'10') then
298  nope=10
299  elseif(lakon(i)(4:4).eq.'4') then
300  nope=4
301  elseif(lakon(i)(4:5).eq.'15') then
302  nope=15
303  elseif(lakon(i)(4:4).eq.'6') then
304  nope=6
305  else
306  cycle
307  endif
308 !
309  do j=1,nope
310  konl(j)=kon(indexe+j)
311  enddo
312 !
313  if(nload.gt.0)
314  & call e_c3d_rhs_th(co,nk,konl,lakon(i),
315  & ff,i,nmethod,t0,t1,vold,nelemload,
316  & sideload,xload,nload,idist,dtime,
317  & ttime,time,istep,iinc,xloadold,reltime,
318  & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,mi,
319  & ielprop,prop,sti,xstateini,xstate,nstate_)
320 !
321 ! modal dynamics and steady state dynamics:
322 ! location of nonzeros is stored
323 !
324  if((nmethod.ge.4.and.(iperturb.lt.2))) then
325  do jj=1,nope
326 !
327  j=jj
328 !
329  node1=kon(indexe+j)
330  jdof1=nactdof(0,node1)
331 !
332 ! distributed forces
333 !
334  if(idist.ne.0) then
335  if(dabs(ff(jj)).lt.1.d-30) cycle
336  if(jdof1.le.0) then
337  if(nmpc.ne.0) then
338  idof1=(node1-1)*8
339  call nident(ikmpc,idof1,nmpc,id)
340  if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
341  id=ilmpc(id)
342  ist=ipompc(id)
343  index=nodempc(3,ist)
344  do
345  jdof1=nactdof(nodempc(2,index),
346  & nodempc(1,index))
347  if(jdof1.gt.0) then
348  fext(jdof1)=fext(jdof1)
349  & -coefmpc(index)*ff(jj)/coefmpc(ist)
350  call nident(ikactmech,jdof1-1,nactmech,
351  & idm)
352  do
353  if(idm.gt.0) then
354  if(ikactmech(idm).eq.jdof1-1) exit
355  endif
356  nactmech=nactmech+1
357  do m=nactmech,idm+2,-1
358  ikactmech(m)=ikactmech(m-1)
359  enddo
360  ikactmech(idm+1)=jdof1-1
361  exit
362  enddo
363  endif
364  index=nodempc(3,index)
365  if(index.eq.0) exit
366  enddo
367  endif
368  endif
369  cycle
370  endif
371  fext(jdof1)=fext(jdof1)+ff(jj)
372  call nident(ikactmech,jdof1-1,nactmech,
373  & idm)
374  do
375  if(idm.gt.0) then
376  if(ikactmech(idm).eq.jdof1-1) exit
377  endif
378  nactmech=nactmech+1
379  do m=nactmech,idm+2,-1
380  ikactmech(m)=ikactmech(m-1)
381  enddo
382  ikactmech(idm+1)=jdof1-1
383  exit
384  enddo
385  endif
386 !
387  enddo
388 !
389 !
390 ! other procedures
391 !
392  else
393  do jj=1,nope
394 !
395  j=jj
396 !
397  node1=kon(indexe+j)
398  jdof1=nactdof(0,node1)
399 !
400 ! distributed forces
401 !
402  if(idist.ne.0) then
403  if(jdof1.le.0) then
404  if(nmpc.ne.0) then
405  idof1=(node1-1)*8
406  call nident(ikmpc,idof1,nmpc,id)
407  if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
408  id=ilmpc(id)
409  ist=ipompc(id)
410  index=nodempc(3,ist)
411  do
412  jdof1=nactdof(nodempc(2,index),
413  & nodempc(1,index))
414  if(jdof1.gt.0) then
415  fext(jdof1)=fext(jdof1)
416  & -coefmpc(index)*ff(jj)/coefmpc(ist)
417  endif
418  index=nodempc(3,index)
419  if(index.eq.0) exit
420  enddo
421  endif
422  endif
423  cycle
424  endif
425  fext(jdof1)=fext(jdof1)+ff(jj)
426  endif
427 !
428  enddo
429  endif
430  enddo
431 !
432  endif
433 !
434 ! point forces
435 !
436 ! modal dynamics and steady state dynamics:
437 ! location of nonzeros is stored
438 !
439  if((nmethod.ge.4).and.(iperturb.lt.2)) then
440  do i=1,nforc
441  if(ndirforc(i).gt.3) cycle
442  if(dabs(xforc(i)).lt.1.d-30) cycle
443  jdof=nactdof(ndirforc(i),nodeforc(1,i))
444  if(jdof.gt.0) then
445  fext(jdof)=fext(jdof)+xforc(i)
446  call nident(ikactmech,jdof-1,nactmech,
447  & idm)
448  do
449  if(idm.gt.0) then
450  if(ikactmech(idm).eq.jdof-1) exit
451  endif
452  nactmech=nactmech+1
453  do m=nactmech,idm+2,-1
454  ikactmech(m)=ikactmech(m-1)
455  enddo
456  ikactmech(idm+1)=jdof-1
457  exit
458  enddo
459  else
460 !
461 ! node is a dependent node of a MPC: distribute
462 ! the forces among the independent nodes
463 ! (proportional to their coefficients)
464 !
465  jdof=8*(nodeforc(1,i)-1)+ndirforc(i)
466  call nident(ikmpc,jdof,nmpc,id)
467  if(id.gt.0) then
468  if(ikmpc(id).eq.jdof) then
469  ist=ipompc(id)
470  index=nodempc(3,ist)
471  if(index.eq.0) cycle
472  do
473  jdof=nactdof(nodempc(2,index),nodempc(1,index))
474  if(jdof.gt.0) then
475  fext(jdof)=fext(jdof)-
476  & coefmpc(index)*xforc(i)/coefmpc(ist)
477  call nident(ikactmech,jdof-1,nactmech,
478  & idm)
479  do
480  if(idm.gt.0) then
481  if(ikactmech(idm).eq.jdof-1) exit
482  endif
483  nactmech=nactmech+1
484  do m=nactmech,idm+2,-1
485  ikactmech(m)=ikactmech(m-1)
486  enddo
487  ikactmech(idm+1)=jdof-1
488  exit
489  enddo
490  endif
491  index=nodempc(3,index)
492  if(index.eq.0) exit
493  enddo
494  endif
495  endif
496  endif
497  enddo
498  else
499 !
500 ! other procedures
501 !
502  do i=1,nforc
503  if(ndirforc(i).gt.3) cycle
504  jdof=nactdof(ndirforc(i),nodeforc(1,i))
505  if(jdof.gt.0) then
506  fext(jdof)=fext(jdof)+xforc(i)
507  else
508 !
509 ! node is a dependent node of a MPC: distribute
510 ! the forces among the independent nodes
511 ! (proportional to their coefficients)
512 !
513  jdof=8*(nodeforc(1,i)-1)+ndirforc(i)
514  call nident(ikmpc,jdof,nmpc,id)
515  if(id.gt.0) then
516  if(ikmpc(id).eq.jdof) then
517  ist=ipompc(id)
518  index=nodempc(3,ist)
519  if(index.eq.0) cycle
520  do
521  jdof=nactdof(nodempc(2,index),nodempc(1,index))
522  if(jdof.gt.0) then
523  fext(jdof)=fext(jdof)-
524  & coefmpc(index)*xforc(i)/coefmpc(ist)
525  endif
526  index=nodempc(3,index)
527  if(index.eq.0) exit
528  enddo
529  endif
530  endif
531  endif
532  enddo
533  endif
534 c write(*,*) 'rhs '
535 c write(*,'(6(1x,e11.4))') (fext(i),i=1,neq)
536 !
537  return
newton
subroutine newton(icalccg, ne, ipkon, lakon, kon, t0, co, rhcon, nrhcon, ntmat_, physcon, nelem, cgr, bodyf, ielmat, ithermal, vold, mi)
Definition: newton.f:22
idist
static ITG * idist
Definition: radflowload.c:39
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
e_c3d_rhs_th
subroutine e_c3d_rhs_th(co, nk, konl, lakonl, ff, nelem, nmethod, t0, t1, vold, nelemload, sideload, xload, nload, idist, dtime, ttime, time, istep, iinc, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, mi, ielprop, prop, sti, xstateini, xstate, nstate_)
Definition: e_c3d_rhs_th.f:25
e_c3d_rhs
subroutine e_c3d_rhs(co, nk, konl, lakonl, p1, p2, omx, bodyfx, nbody, ff, nelem, nmethod, rhcon, ielmat, ntmat_, vold, iperturb, nelemload, sideload, xload, nload, idist, ttime, time, istep, iinc, dtime, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, veold, matname, mi, ielprop, prop)
Definition: e_c3d_rhs.f:24
Hosted by OpenAircraft.com, (Michigan UAV, LLC)