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