CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Macros | Functions
checkconvergence.c File Reference
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "CalculiX.h"
Include dependency graph for checkconvergence.c:

Go to the source code of this file.

Macros

#define max(a, b)   ((a) >= (b) ? (a) : (b))
 

Functions

void checkconvergence (double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, double *stn, ITG *nmethod, ITG *kode, char *filab, double *een, double *t1act, double *time, double *epn, ITG *ielmat, char *matname, double *enern, double *xstaten, ITG *nstate_, ITG *istep, ITG *iinc, ITG *iperturb, double *ener, ITG *mi, char *output, ITG *ithermal, double *qfn, ITG *mode, ITG *noddiam, double *trab, ITG *inotr, ITG *ntrans, double *orab, ITG *ielorien, ITG *norien, char *description, double *sti, ITG *icutb, ITG *iit, double *dtime, double *qa, double *vold, double *qam, double *ram1, double *ram2, double *ram, double *cam, double *uam, ITG *ntg, double *ttime, ITG *icntrl, double *theta, double *dtheta, double *veold, double *vini, ITG *idrct, double *tper, ITG *istab, double *tmax, ITG *nactdof, double *b, double *tmin, double *ctrl, double *amta, ITG *namta, ITG *itpamp, ITG *inext, double *dthetaref, ITG *itp, ITG *jprint, ITG *jout, ITG *uncoupled, double *t1, ITG *iitterm, ITG *nelemload, ITG *nload, ITG *nodeboun, ITG *nboun, ITG *itg, ITG *ndirboun, double *deltmx, ITG *iflagact, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, double *emn, double *thicke, char *jobnamec, ITG *mortar, ITG *nmat, ITG *ielprop, double *prop, ITG *ialeatoric, ITG *kscale, double *energy, double *allwk, double *energyref, double *emax, double *r_abs, double *enetoll, double *energyini, double *allwkini, double *temax, double *sizemaxinc, ITG *ne0, ITG *neini, double *dampwk, double *dampwkini, double *energystartstep)
 

Macro Definition Documentation

◆ max

#define max (   a,
  b 
)    ((a) >= (b) ? (a) : (b))

Function Documentation

◆ checkconvergence()

void checkconvergence ( double *  co,
ITG *  nk,
ITG *  kon,
ITG *  ipkon,
char *  lakon,
ITG *  ne,
double *  stn,
ITG *  nmethod,
ITG *  kode,
char *  filab,
double *  een,
double *  t1act,
double *  time,
double *  epn,
ITG *  ielmat,
char *  matname,
double *  enern,
double *  xstaten,
ITG *  nstate_,
ITG *  istep,
ITG *  iinc,
ITG *  iperturb,
double *  ener,
ITG *  mi,
char *  output,
ITG *  ithermal,
double *  qfn,
ITG *  mode,
ITG *  noddiam,
double *  trab,
ITG *  inotr,
ITG *  ntrans,
double *  orab,
ITG *  ielorien,
ITG *  norien,
char *  description,
double *  sti,
ITG *  icutb,
ITG *  iit,
double *  dtime,
double *  qa,
double *  vold,
double *  qam,
double *  ram1,
double *  ram2,
double *  ram,
double *  cam,
double *  uam,
ITG *  ntg,
double *  ttime,
ITG *  icntrl,
double *  theta,
double *  dtheta,
double *  veold,
double *  vini,
ITG *  idrct,
double *  tper,
ITG *  istab,
double *  tmax,
ITG *  nactdof,
double *  b,
double *  tmin,
double *  ctrl,
double *  amta,
ITG *  namta,
ITG *  itpamp,
ITG *  inext,
double *  dthetaref,
ITG *  itp,
ITG *  jprint,
ITG *  jout,
ITG *  uncoupled,
double *  t1,
ITG *  iitterm,
ITG *  nelemload,
ITG *  nload,
ITG *  nodeboun,
ITG *  nboun,
ITG *  itg,
ITG *  ndirboun,
double *  deltmx,
ITG *  iflagact,
char *  set,
ITG *  nset,
ITG *  istartset,
ITG *  iendset,
ITG *  ialset,
double *  emn,
double *  thicke,
char *  jobnamec,
ITG *  mortar,
ITG *  nmat,
ITG *  ielprop,
double *  prop,
ITG *  ialeatoric,
ITG *  kscale,
double *  energy,
double *  allwk,
double *  energyref,
double *  emax,
double *  r_abs,
double *  enetoll,
double *  energyini,
double *  allwkini,
double *  temax,
double *  sizemaxinc,
ITG *  ne0,
ITG *  neini,
double *  dampwk,
double *  dampwkini,
double *  energystartstep 
)
59  {
60 
61  ITG i0,ir,ip,ic,il,ig,ia,iest,iest1=0,iest2=0,iconvergence,idivergence,
62  ngraph=1,k,*ipneigh=NULL,*neigh=NULL,*inum=NULL,id,istart,iend,inew,
63  i,j,mt=mi[1]+1,iexceed,
64  iforceincsize=0;
65 
66  double df,dc,db,dd,ran,can,rap,ea,cae,ral,da,*vr=NULL,*vi=NULL,*stnr=NULL,
67  *stni=NULL,*vmax=NULL,*stnmax=NULL,*cs=NULL,c1[2],c2[2],reftime,
68  *fn=NULL,*eenmax=NULL,*fnr=NULL,*fni=NULL,*qfx=NULL,*cdn=NULL,
69  *cdnr=NULL,*cdni=NULL,tmp, maxdecay=0.0, r_rel,cetol;
70 
71  /* reset ialeatoric to zero */
72 
73  *ialeatoric=0;
74 
75  cetol=ctrl[39];
76 
77  /* next lines are active if the number of contact elements was
78  changed in the present increment */
79 
80  if ((*iflagact==1)&&(*mortar!=1)){
81  if(ctrl[0]<*iit+4)ctrl[0]=*iit+4;
82  if(ctrl[1]<*iit+8)ctrl[1]=*iit+8;
83  ctrl[3]+=1;
84  }
85 
86  i0=ctrl[0];ir=ctrl[1];ip=ctrl[2];ic=ctrl[3];il=ctrl[4];ig=ctrl[5];
87  ia=ctrl[7];df=ctrl[10];dc=ctrl[11];db=ctrl[12];da=ctrl[13];dd=ctrl[16];
88  ran=ctrl[18];can=ctrl[19];rap=ctrl[22];
89  ea=ctrl[23];cae=ctrl[24];ral=ctrl[25];
90 
91  /* for face-to-face penalty contact: increase the number of iterations
92  in two subsequent increments in order to increase the increment size */
93 
94  if(*mortar==1){ig+=12;il+=12;}
95 
96  /* if iconvergence=0 the increment did not yet converge, iterations are
97  continued
98  if idivergence=1 the increment diverged and has to be reiterated
99  with a smaller size */
100 
101  idivergence=0;
102 
103  /* check for forced divergence (due to divergence of a user material
104  routine */
105 
106  if(qa[2]>0.){idivergence=1;}
107 
108  if(*ithermal!=2){
109  if(qa[0]>ea*qam[0]){
110  if(*iit<=ip){c1[0]=ran;}
111  else{c1[0]=rap;}
112  c2[0]=can;
113  }
114  else{
115  c1[0]=ea;
116  c2[0]=cae;
117  }
118  if(ram1[0]<ram2[0]){ram2[0]=ram1[0];}
119  }
120  if(*ithermal>1){
121  if(qa[1]>ea*qam[1]){
122  if(*iit<=ip){c1[1]=ran;}
123  else{c1[1]=rap;}
124  c2[1]=can;
125  }
126  else{
127  c1[1]=ea;
128  c2[1]=cae;
129  }
130  if(ram1[1]<ram2[1]){ram2[1]=ram1[1];}
131  }
132 
133  iconvergence=0;
134 
135  /* mechanical */
136 
137  if(*ithermal<2){
138  if((*iit>1)&&(ram[0]<=c1[0]*qam[0])&&(*iflagact==0)&&
139  ((*nmethod!=-1)||(qa[3]<=cetol))&&
140  ((cam[0]<=c2[0]*uam[0])||
141  (((ram[0]*cam[0]<c2[0]*uam[0]*ram2[0])||(ram[0]<=ral*qam[0])||
142  (qa[0]<=ea*qam[0]))&&(*ntg==0))||
143  (cam[0]<1.e-8))) iconvergence=1;
144  }
145 
146  /* thermal */
147 
148  if(*ithermal==2){
149  if((ram[1]<=c1[1]*qam[1])&&
150  (cam[2]<*deltmx)&&
151  ((cam[1]<=c2[1]*uam[1])||
152  (((ram[1]*cam[1]<c2[1]*uam[1]*ram2[1])||(ram[1]<=ral*qam[1])||
153  (qa[1]<=ea*qam[1]))&&(*ntg==0))||
154  (cam[1]<1.e-8)))iconvergence=1;
155  }
156 
157  /* thermomechanical */
158 
159  if(*ithermal==3){
160  if(((*iit>1)&&(ram[0]<=c1[0]*qam[0])&&
161  ((*nmethod!=-1)||(qa[3]<=cetol))&&
162  ((cam[0]<=c2[0]*uam[0])||
163  (((ram[0]*cam[0]<c2[0]*uam[0]*ram2[0])||(ram[0]<=ral*qam[0])||
164  (qa[0]<=ea*qam[0]))&&(*ntg==0))||
165  (cam[0]<1.e-8)))&&
166  ((ram[1]<=c1[1]*qam[1])&&
167  (cam[2]<*deltmx)&&
168  ((cam[1]<=c2[1]*uam[1])||
169  (((ram[1]*cam[1]<c2[1]*uam[1]*ram2[1])||(ram[1]<=ral*qam[1])||
170  (qa[1]<=ea*qam[1]))&&(*ntg==0))||
171  (cam[1]<1.e-8))))iconvergence=1;
172  }
173 
174  /* reset kscale */
175 
176  if(iconvergence==1){
177  if(*kscale>1){
178  *kscale=1;
179  iconvergence=0;
180  printf("\n restoring the elastic contact stifnesses to their original values \n\n");
181  }
182  }
183 
184 // # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
185 // MPADD start
186  /*
187  Energy conservation convergence: only for implicit dynamic calculations
188  (convergence is not checked in explicit dynamics)
189 
190  Main variables and meaning
191 
192  r_rel : modified energy conservation criterion
193  emax : maximum value of the energy over the time history
194  r : energy residual (Eint + K + Econt - Wext - Wdamp -Eref)
195  enetoll : energy conservation tolerance
196  tmp : auxiliary temporary variable
197  maxdecay : \hat(r)^{max}(\theta) -> value of the decay boundary
198  for r_hat.
199 
200  Proposed by Matteo Pacher */
201 
202  if((*nmethod==4)&&(*ithermal<2)&&(*uncoupled==0)&&(iconvergence==1)&&(*ne==*ne0)&&(*neini==*ne0)&&(*idrct==0)) {
203 
204  /* Update the value of the maximum energy of the system emax
205  (contact energy is not taken into account because small) */
206 
207  *emax=max(*emax,fabs(energy[0]-energystartstep[0]));
208  *emax=max(*emax,fabs(energy[1]));
209  *emax=max(*emax,fabs(*allwk));
210 
211  // energy residual (only calculated in the absence of contact);
212  // if <=0: loss of energy
213 
214  *r_abs=energy[0]+energy[1]+energy[2]+energy[3]-*energyref-*allwk-*dampwk;
215 
216  // Absolute tolerance check (when the error is really small --> beginning of simulation)
217 
218  if(fabs(*r_abs)>=*enetoll/4) {
219 
220  // Normal strategy: Relative error
221 
222  /* Compute admissible decay*/
223 
224  maxdecay=*enetoll/2*(1+sqrt(*theta));
225 
226  /* modified r_hat criterion */
227 
228  r_rel=*r_abs/(*emax);
229  if(r_rel<=-maxdecay) {
230  idivergence=1;
231  }else{
232 
233  /* Check if the residual is too close to the boundary */
234 
235  if(r_rel<=-0.9*maxdecay) {
236  *istab=0; // keep the increment size
237  }
238  }
239  }
240  }
241 
242  /* Contact Strategy: limit jumps and time increment during contact based
243  on the natural frequency of oscillation of contact elements
244  Implicit dynamic calculations only */
245 
246  if((*nmethod==4)&&(*ithermal<2)&&(*uncoupled==0)&&(iconvergence==1)&&((*ne!=*ne0)||(*neini!=*ne0))){
247 
248  /* store temporarly the value of emax: in case of forced divergence
249  emax has to be reset. */
250 
251  tmp=*emax;
252 
253  /* Update the value of the maximum energy of the system emax
254  (contact energy is not taken into account because small) */
255 
256  *emax=max(*emax,fabs(energy[0]-energystartstep[0]));
257  *emax=max(*emax,fabs(energy[1]));
258  *emax=max(*emax,fabs(*allwk));
259 
260  /* maximum decay boundary */
261 
262  maxdecay=*enetoll/2*(1+sqrt(*theta));
263 
264  FORTRAN(checkimpacts,(ne,neini,temax,sizemaxinc,energyref,
265  tmin,tper,&idivergence,
266  &iforceincsize,istab,dtheta,r_abs,energy,energyini,
267  allwk,allwkini,dampwk,dampwkini,emax,mortar,
268  &maxdecay,enetoll));
269 
270  /* reset emax in case of forced divergence */
271 
272  if(idivergence==1){
273  *emax=tmp;
274  }
275 
276  /* Adaption of the energy tolerance in case of violation due to
277  contact jumps (rebounds).
278  The user is aware of it via the output string. */
279 
280  if((*ne==*ne0)&&(*neini>*ne0)&&(idivergence==0)){
281  *r_abs=fabs(energy[0]+energy[1]+energy[2]+energy[3]-*energyref-*allwk-*dampwk);
282  tmp=1.3*(2.0*(*r_abs)/(*emax))/(1.0+sqrt(*theta));
283  *enetoll=max(*enetoll,tmp);
284  printf("\n Adaption of the max-decay boundary, enetoll = %f \n",*enetoll);
285  }
286 
287  /*
288  Adaption of the energy residual during long periods of persistent contact.
289  Take care of the general (increasing) trend of the residual to avoid code stucks.
290  */
291 
292  if((iconvergence==1)&&(*ne<=*neini)){
293  tmp=energy[0]+energy[1]+energy[2]+energy[3]-*energyref-*allwk-*dampwk;
294  if(tmp>*r_abs){
295  *r_abs=tmp;
296  printf("\n Adaption of the energy residual in persistent contact, \n");
297  printf(" an increasing trend has been detected.\n");
298  }
299  }
300  } else {
301  *sizemaxinc=*tmax;
302  }
303 // MPADD end
304 // # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
305 
306  /* increment convergence reached */
307 
308  if((iconvergence==1)&&(idivergence==0)){
309 
310  FORTRAN(writesummary,(istep,iinc,icutb,iit,ttime,time,dtime));
311  if(*uncoupled){
312  if(*ithermal==2){
313  *iitterm=*iit;
314  *ithermal=1;
315  for(k=0;k<*nk;++k){t1[k]=vold[mt*k];}
316  *iit=1;
317  (ctrl[0])*=4;
318  printf(" thermal convergence\n\n");
319  *iflagact=0;
320  return;
321  }else{
322  *ithermal=3;
323  *iit=*iitterm;
324  (ctrl[0])/=4;
325  }
326  }
327 
328  *icntrl=1;
329  *icutb=0;
330  *theta=*theta+*dtheta;
331 
332  /* defining a mean "velocity" for static calculations: is used to
333  extrapolate the present results for next increment */
334 
335  if(*nmethod != 4){
336  for(i=0;i<*nk;i++){
337  for(j=1;j<mt;j++){
338  veold[mt*i+j]=(vold[mt*i+j]-vini[mt*i+j])/(*dtime);
339  }
340  }
341  }
342 
343  /* check whether size is to be set to a fixed value */
344 
345  if(iforceincsize==1){
346  *dtheta=*sizemaxinc;
347  *dthetaref=*sizemaxinc;
348  printf(" convergence; new increment size is forced to %e\n\n",*dtheta**tper);
349  }
350 
351  /* check whether next increment size must be decreased */
352 
353  else if((*iit>il)&&(*idrct==0)){
354  if(*mortar==0){
355  *dtheta=*dthetaref*db;
356  *dthetaref=*dtheta;
357  printf(" convergence; the increment size is decreased to %e\n\n",*dtheta**tper);
358  if(*dtheta<*tmin){
359  printf("\n *ERROR: increment size smaller than minimum\n");
360  printf(" best solution and residuals are in the frd file\n\n");
361  NNEW(fn,double,mt**nk);
362  NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
363  FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
364  nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
365  ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
366  ++*kode;
367 
368  (*ttime)+=(*time);
369  frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
370  kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
371  xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
372  trab,inotr,ntrans,orab,ielorien,norien,description,
373  ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
374  &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
375  ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,
376  cdn,mortar,cdnr,cdni,nmat);
377 
378  FORTRAN(stop,());
379  }
380  }
381  else{
382  printf("convergence\n\n");}
383  }
384 
385  /* check whether next increment size can be increased */
386 
387  else if(*iit<=ig){
388  if((*istab==1)&&(*idrct==0)){
389  *dtheta=*dthetaref*dd;
390  *dthetaref=*dtheta;
391  printf(" convergence; the increment size is increased to %e\n\n",*dtheta**tper);
392  }
393  else{
394  *istab=1;
395  printf(" convergence\n\n");
396  *dtheta=*dthetaref;
397  }
398  }
399  else{
400  *istab=0;
401  printf(" convergence\n\n");
402  *dtheta=*dthetaref;
403  }
404 
405  /* check whether new increment size exceeds maximum increment
406  size allowed (by the user) */
407 
408  if((*dtheta>*sizemaxinc)&&(*idrct==0)){
409  *dtheta=*sizemaxinc;
410  *dthetaref=*dtheta;
411  printf(" the increment size exceeds thetamax and is decreased to %e\n\n",*dtheta**tper);
412  }
413 
414  /* check whether new time point exceeds end of step */
415 
416  if(*dtheta>=1.-*theta){
417  if(*dtheta>1.-*theta){iexceed=1;}else{iexceed=0;}
418  *dtheta=1.-*theta;
419  *dthetaref=*dtheta;
420  if(iexceed==1)
421  printf(" the increment size exceeds the remainder of the step and is decreased to %e\n\n",*dtheta**tper);
422  }
423 
424  /* check whether the end of the new increment exceeds a time point;
425  if itp=1 the increment just finished ends at a time point */
426 
427  if((*itpamp>0)&&(*idrct==0)){
428  if(*itp==1){
429  *jprint=*jout;
430  }else{
431  *jprint=*jout+1;
432  }
433  if(namta[3**itpamp-1]<0){
434  reftime=*ttime+*time+*dtheta**tper;
435  }else{
436  reftime=*time+*dtheta**tper;
437  }
438  istart=namta[3**itpamp-3];
439  iend=namta[3**itpamp-2];
440  FORTRAN(identamta,(amta,&reftime,&istart,&iend,&id));
441  if(id<istart){
442  inew=istart;
443  }else{
444  inew=id+1;
445  }
446 
447  /* if the end of the new increment is less than a time
448  point by less than 1.e-6 (theta-value) dtheta is
449  enlarged up to this time point */
450 
451  if((*inext==inew)&&(inew<=iend)){
452  if(amta[2*inew-2]-reftime<1.e-6**tper){inew++;}
453  }
454 
455  /* inew: smallest time point exceeding time+dtheta*tper
456  inext: smallest time point exceeding time */
457 
458  if(*inext<inew){
459  if(namta[3**itpamp-1]<0){
460  *dtheta=(amta[2**inext-2]-*ttime-*time)/(*tper);
461  }else{
462  *dtheta=(amta[2**inext-2]-*time)/(*tper);
463  }
464  (*inext)++;
465  *itp=1;
466  printf(" the increment size exceeds a time point and is decreased to %e\n\n",*dtheta**tper);
467  }else{*itp=0;}
468  }
469  }
470  else{
471 
472  /* no convergence */
473 
474  /* check for the amount of iterations */
475 
476  if(((*iit>ic)&&(*mortar==0))||((*mortar>1)&&(*iit>200))){
477  printf("\n *ERROR: too many iterations needed\n");
478  printf(" best solution and residuals are in the frd file\n\n");
479 
480  FORTRAN(writesummarydiv,(istep,iinc,icutb,iit,ttime,time,dtime));
481 
482  NNEW(fn,double,mt**nk);
483  NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
484  FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,nk,sti,stn,
485  ipkon,inum,kon,lakon,ne,mi,orab,ielorien,co,itg,ntg,vold,
486  ielmat,thicke,ielprop,prop));
487  ++*kode;
488 
489  (*ttime)+=(*time);
490  frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
491  kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
492  xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
493  trab,inotr,ntrans,orab,ielorien,norien,description,
494  ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
495  &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
496  ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
497  mortar,cdnr,cdni,nmat);
498 
499  FORTRAN(stop,());
500  }
501 
502  /* check for diverging residuals */
503 
504  /* if the user has not defined deltmx on the *HEAT
505  TRANSFER card it is set to a large value (1.e30,
506  cf. CalculiX.c); therefore, a comparison of cam[2]
507  with deltmx only makes sense for cam[2]<1.e30 */
508 
509  if((*iit>=i0)||(fabs(ram[0])>1.e20)||(fabs(cam[0])>1.e20)||
510  (fabs(ram[1])>1.e20)||(fabs(cam[1])>1.e20)||
511  ((cam[2]<1.e30)&&(cam[2]>*deltmx))||(idivergence==1)||
512  (iforceincsize==1)){
513  if((*ithermal!=2)&&(*mortar!=1)){
514  if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
515  idivergence=1;
516  }
517 
518  if((*ithermal!=2)&&(*mortar==1)){
519 
520  if(ram[0]>1.e9){
521  printf("divergence allowed: residual force too large\n");
522  if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
523  idivergence=1;
524  }
525 
526  /* number of contact elements does not change */
527 
528  if(*iflagact==0){
529  printf("divergence allowed: number of contact elements stabilized\n");
530  if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0])){
531  if(cam[0]<=c2[0]*uam[0]){
532  *ialeatoric=1;
533  }
534  }
535  }
536 
537  /* rate of number of contact elements is increasing */
538 
539  if(((ITG)ram[6]*(ITG)ram1[6]<0)&&((ITG)ram1[6]*(ITG)ram2[6]<0)){
540 
541  if(((ram[4]>0.98*ram1[4])&&(ram[4]<1.02*ram1[4]))&&
542  ((ram[4]>0.98*ram2[4])&&(ram[4]<1.02*ram2[4]))){
543  printf("divergence allowed: repetitive pattern detected\n");
544  if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
545  idivergence=1;
546  }
547  }
548  }
549 
550  /* check whether in a viscous step the allowable increase in viscous
551  strain has been exceeded */
552 
553  if((idivergence==0)&&((*nmethod==-1)&&(qa[3]>cetol))) idivergence=2;
554 
555 
556  /* for thermal calculations the maximum temperature change
557  is checked as well */
558 
559  if(*ithermal>1){
560  if((ram1[1]>ram2[1])&&(ram[1]>ram2[1])&&(ram[1]>c1[1]*qam[1]))
561  idivergence=1;
562 
563  /* if the user has not defined deltmx on the *HEAT
564  TRANSFER card it is set to a large value (1.e30,
565  cf. CalculiX.c); therefore, a comparison of cam[2]
566  with deltmx only makes sense for cam[2]<1.e30 */
567 
568  if((cam[2]<1.e30)&&(cam[2]>*deltmx)) idivergence=2;
569  }
570 
571  if(idivergence>0){
572  if(*idrct==1){
573  if((*mortar<=1)||((*mortar>1)&&(*iit>200))) {
574 
575  /* fixed time increments */
576 
577  printf("\n *ERROR: solution seems to diverge; please try \n");
578  printf(" automatic incrementation; program stops\n");
579  printf(" best solution and residuals are in the frd file\n\n");
580 
581  FORTRAN(writesummarydiv,(istep,iinc,icutb,iit,ttime,time,
582  dtime));
583 
584  NNEW(fn,double,mt**nk);
585  NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
586  FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,nk,
587  sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
588  ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
589  ++*kode;
590 
591  (*ttime)+=(*time);
592  frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
593  kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
594  xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
595  trab,inotr,ntrans,orab,ielorien,norien,description,
596  ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
597  &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
598  ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
599  mortar,cdnr,cdni,nmat);
600 
601  FORTRAN(stop,());
602  }
603  }
604  else {
605 
606  /* variable time increments */
607 
608  if(qa[2]>0.){
609  *dtheta=*dtheta*qa[2];
610  printf("increment size decrease requested by a material user routine (through pnewdt)\n\n");
611  }else{
612  if(idivergence==1){
613  if((*mortar!=1)||(*icutb!=0)){ // MPADD
614  if(iforceincsize != 1){ // MPADD
615  *dtheta=*dtheta*df; // MPADD
616  }else{ // MPADD
617  *dtheta=*sizemaxinc; // MPADD
618  } // MPADD
619  } // MPADD
620  }else{
621  if(*nmethod==-1){
622  *dtheta=*dtheta*cetol/qa[3]*da;
623  }else{
624  *dtheta=*dtheta**deltmx/cam[2]*da;
625  }
626  }
627  }
628  *dthetaref=*dtheta;
629  printf(" divergence; the increment size is decreased to %e\n",*dtheta**tper);
630  printf(" the increment is reattempted\n\n");
631 
632  FORTRAN(writesummarydiv,(istep,iinc,icutb,iit,ttime,time,
633  dtime));
634 
635  *istab=0;
636  if(*itp==1){
637  *itp=0;
638  (*inext)--;
639  }
640 
641  /* check whether new increment size is smaller than minimum */
642 
643  if(*dtheta<*tmin){
644  printf("\n *ERROR: increment size smaller than minimum\n");
645  printf(" best solution and residuals are in the frd file\n\n");
646  NNEW(fn,double,mt**nk);
647  NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
648  FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
649  nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
650  ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
651  ++*kode;
652 
653  (*ttime)+=(*time);
654  frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
655  kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
656  xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
657  trab,inotr,ntrans,orab,ielorien,norien,description,
658  ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
659  &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
660  ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
661  mortar,cdnr,cdni,nmat);
662 
663  FORTRAN(stop,());
664  }
665  *icntrl=1;
666  (*icutb)++;
667  if(*mortar==1){
668  *kscale=100;
669  printf("\n reducing the constant stiffnesses by a factor of 100 \n\n");
670  }
671 
672  /* check whether too many cutbacks */
673 
674  if(*icutb>ia){
675  printf("\n *ERROR: too many cutbacks\n");
676  printf(" best solution and residuals are in the frd file\n\n");
677  NNEW(fn,double,mt**nk);
678  NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
679  FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
680  nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
681  ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
682  ++*kode;
683 
684  (*ttime)+=(*time);
685  frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
686  kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
687  xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
688  trab,inotr,ntrans,orab,ielorien,norien,description,
689  ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
690  &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
691  ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
692  mortar,cdnr,cdni,nmat);
693 
694  FORTRAN(stop,());
695  }
696  if(*uncoupled){
697  if(*ithermal==1){
698  (ctrl[0])/=4;
699  }
700  *ithermal=3;
701  }
702 
703  /* default value for qa[2] */
704 
705  qa[2]=-1.;
706 
707  *iflagact=0;
708  return;
709  }
710  }
711  }
712 
713  /* check for too slow convergence */
714 
715  if(*iit>=ir){
716  if(*ithermal!=2){
717  iest1=(ITG)ceil(*iit+log(ran*qam[0]/(ram[0]))/
718  log(ram[0]/(ram1[0])));
719  }
720  if(*ithermal>1){
721  iest2=(ITG)ceil(*iit+log(ran*qam[1]/(ram[1]))/
722  log(ram[1]/(ram1[1])));
723  }
724  if(iest1>iest2){iest=iest1;}else{iest=iest2;}
725  if((iest>0)&&(*mortar!=1)){
726  printf(" estimated number of iterations till convergence = %" ITGFORMAT "\n",
727  iest);
728  }
729  if((((iest>ic)||(*iit==ic))&&(*mortar!=1))||((*mortar==1)&&(*iit==60))){
730 
731  if(*idrct!=1){
732  if((*mortar!=1)||(*icutb!=0)) *dtheta=*dtheta*dc;
733  *dthetaref=*dtheta;
734  printf(" too slow convergence; the increment size is decreased to %e\n",*dtheta**tper);
735  printf(" the increment is reattempted\n\n");
736 
737  FORTRAN(writesummarydiv,(istep,iinc,icutb,iit,ttime,
738  time,dtime));
739  *istab=0;
740  if(*itp==1){
741  *itp=0;
742  (*inext)--;
743  }
744 
745  /* check whether new increment size is smaller than minimum */
746 
747  if(*dtheta<*tmin){
748  printf("\n *ERROR: increment size smaller than minimum\n");
749  printf(" best solution and residuals are in the frd file\n\n");
750  NNEW(fn,double,mt**nk);
751  NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
752  FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
753  nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
754  ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
755  ++*kode;
756 
757  (*ttime)+=(*time);
758  frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
759  kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
760  xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
761  trab,inotr,ntrans,orab,ielorien,norien,description,
762  ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
763  &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
764  ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
765  mortar,cdnr,cdni,nmat);
766 
767  FORTRAN(stop,());
768  }
769  *icntrl=1;
770  (*icutb)++;
771  if(*mortar==1){
772  *kscale=100;
773  printf("\n reducing the constant stiffnesses by a factor of 100 \n\n");
774  }
775 // if(*mortar==1) *kscale=100;
776 
777  if(*icutb>ia){
778  printf("\n *ERROR: too many cutbacks\n");
779  printf(" best solution and residuals are in the frd file\n\n");
780  NNEW(fn,double,mt**nk);
781  NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
782  FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
783  nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
784  ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
785  ++*kode;
786 
787  (*ttime)+=(*time);
788  frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
789  kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
790  xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
791  trab,inotr,ntrans,orab,ielorien,norien,description,
792  ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
793  &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
794  ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
795  mortar,cdnr,cdni,nmat);
796 
797  FORTRAN(stop,());
798  }
799  if(*uncoupled){
800  if(*ithermal==1){
801  (ctrl[0])/=4;
802  }
803  *ithermal=3;
804  }
805 
806  /* default value for qa[2] */
807 
808  qa[2]=-1.;
809 
810  *iflagact=0;
811  return;
812  }
813  }
814  }
815 
816  printf(" no convergence\n\n");
817 
818  (*iit)++;
819 
820  }
821 
822  *iflagact=0;
823  return;
824 }
checkimpacts
subroutine checkimpacts(ne, neini, temax, sizemaxinc, energyref, tmin, tper, idivergence, iforceincsize, istab, dtheta, r_abs, energy, energyini, allwk, allwkini, dampwk, dampwkini, emax, mortar, maxdecay, enetoll)
Definition: checkimpacts.f:23
ITGFORMAT
#define ITGFORMAT
Definition: CalculiX.h:52
frd
void frd(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne0, double *v, double *stn, ITG *inum, ITG *nmethod, ITG *kode, char *filab, double *een, double *t1, double *fn, double *time, double *epn, ITG *ielmat, char *matname, double *enern, double *xstaten, ITG *nstate_, ITG *istep, ITG *iinc, ITG *ithermal, double *qfn, ITG *mode, ITG *noddiam, double *trab, ITG *inotr, ITG *ntrans, double *orab, ITG *ielorien, ITG *norien, char *description, ITG *ipneigh, ITG *neigh, ITG *mi, double *stx, double *vr, double *vi, double *stnr, double *stni, double *vmax, double *stnmax, ITG *ngraph, double *veold, double *ener, ITG *ne, double *cs, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, double *eenmax, double *fnr, double *fni, double *emn, double *thicke, char *jobnamec, char *output, double *qfx, double *cdn, ITG *mortar, double *cdnr, double *cdni, ITG *nmat)
Definition: frd.c:32
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
writesummary
subroutine writesummary(istep, j, icutb, l, ttime, time, dtime)
Definition: writesummary.f:20
storeresidual
subroutine storeresidual(nactdof, b, fn, filab, ithermal, nk, sti, stn, ipkon, inum, kon, lakon, ne, mi, orab, ielorien, co, itg, ntg, vold, ielmat, thicke, ielprop, prop)
Definition: storeresidual.f:22
FORTRAN
void FORTRAN(actideacti,(char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, char *objectset, ITG *ipkon, ITG *ibject, ITG *ne))
c1
static double * c1
Definition: mafillvcompmain.c:30
stop
subroutine stop()
Definition: stop.f:20
writesummarydiv
subroutine writesummarydiv(istep, j, icutb, l, ttime, time, dtime)
Definition: writesummarydiv.f:20
identamta
subroutine identamta(amta, reftime, istart, iend, id)
Definition: identamta.f:26
ITG
#define ITG
Definition: CalculiX.h:51
max
#define max(a, b)
Definition: checkconvergence.c:32
NNEW
#define NNEW(a, b, c)
Definition: CalculiX.h:39
Hosted by OpenAircraft.com, (Michigan UAV, LLC)