00001
00002
00003
00004
00005
00006 #include "elemhead.h"
00007
00008
00009
00010
00011 #include "gadaptivity.h"
00012
00013 #include "element.h"
00014 #include "intpoints.h"
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "adaptivity.h"
00026 #include "z2_smoothing.h"
00027
00028 #include "global.h"
00029 #include "elemswitch.h"
00030
00031
00032
00033 adaptivity :: adaptivity ()
00034 {
00035 tad = Mp->adaptivityflag;
00036 printflags = 0;
00037 corr = 1.0;
00038 otherflags = -1;
00039
00040 path = filename = suffix = NULL;
00041 ni = new char[1]; ni[0] = '\0';
00042
00043 refined_ders_in_nodes = NULL;
00044
00045 answer = 0;
00046 }
00047
00048 adaptivity :: ~adaptivity ()
00049 {
00050 delete [] ni;
00051 delete [] path;
00052 delete [] filename;
00053 delete [] suffix;
00054 }
00055
00056
00057 void adaptivity :: readinit (XFILE *in)
00058 {
00059 xfscanf (in, "%k%ld %k%lf %k%lf", "printflag", &printflags, "accuracy", &adapt_accuracy, "correction", &corr);
00060
00061 if (Mp->adaptivityflag<1 || Mp->adaptivityflag>2)
00062 { print_err("unknown type of adaptivity is required",__FILE__,__LINE__,__func__); exit (0); }
00063
00064 filename_decomposition (in->fname, path, filename, suffix);
00065 }
00066
00067
00068 void adaptivity :: printinit (FILE *out)
00069 {
00070 fprintf (out,"%ld %ld %.3lf %.3lf\n", tad, printflags, adapt_accuracy, corr);
00071 }
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 long adaptivity :: run (long of, int width, long nincr)
00083 {
00084 long hod,min;
00085 double sec;
00086
00087 if (Mespr==1){
00088 fprintf (stdout,"\n\n *********************************");
00089 fprintf (stdout, "\n *** Adaptivity is started ***");
00090
00091 sec = clock();
00092 }
00093
00094
00095 dim = Mt->give_dimension (0);
00096 ord = Mt->give_degree (0);
00097 ncomp = Mt->give_ncomp (0,0);
00098 nn = Mt->nn;
00099 ne = Mt->ne;
00100
00101
00102 for (long i=0; i<ne; i++) {
00103 if (ord != Mt->give_degree (i)) { print_err("mismatch in order", __FILE__, __LINE__, __func__); exit (1); }
00104 if (dim != Mt->give_dimension (i)) { print_err("mismatch in dimension", __FILE__, __LINE__, __func__); exit (1); }
00105 Gtm->gelements[i].auxinf = 100*Mt->give_nne(i) + 10*ord + dim;
00106 }
00107
00108 otherflags = of;
00109 answer = 0;
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120 refined_ders_in_nodes = new vector[nn];
00121 for (long i=0; i<nn; i++) allocv (ncomp, refined_ders_in_nodes[i]);
00122
00123 prepare_ni (width, nincr);
00124 this->check_consistency ();
00125
00126 switch (tad){
00127 case 1:{
00128 spr ();
00129 compute_error ();
00130 break;
00131 }
00132 case 2:{
00133 z2_smoothing z2_est(ncomp);
00134 z2_est.run (refined_ders_in_nodes);
00135 compute_error ();
00136 break;
00137 }
00138
00139
00140
00141
00142
00143 default: print_err("unknown type of error smoothing is required in function", __FILE__, __LINE__, __func__); exit (1);
00144 }
00145
00146
00147
00148
00149
00150 if (printflags & 16)
00151 print_test ();
00152
00153
00154 delete [] refined_ders_in_nodes;
00155
00156 if (Mespr==1){
00157 sec = (clock() - sec) / (double)CLOCKS_PER_SEC;
00158 hod = (long)sec/3600; sec -= hod*3600;
00159 min = (long)sec/60; sec -= min*60;
00160 fprintf(stdout,"\n Consumed time by Ada %2ld:%02ld:%05.2f",hod,min,sec);
00161
00162 fprintf (stdout,"\n *** Adaptivity is finished ***");
00163 fprintf (stdout,"\n *********************************");
00164 }
00165 return (answer);
00166 }
00167
00168
00169
00170 void adaptivity :: prepare_ni (int width, long nincr)
00171 {
00172 ni = new char[1+width+1];
00173
00174 int n = atol(suffix+1);
00175
00176 if (n) sprintf (ni,".%0*d", width, n);
00177 else ni[0] = '\0';
00178
00179
00180 if (strncmp(filename, "this.sifel", 10) == 0)
00181 filename[4] = '\0';
00182 }
00183
00184
00185
00186
00187
00188
00189
00190 void adaptivity :: check_consistency (void) const
00191 {
00192 int te = Mt->give_elem_type (0);
00193 for (long i=1; i<ne; i++)
00194 if (Mt->give_elem_type(i) != te) {
00195 if (te == planeelementlt && Mt->give_elem_type(i) == planeelementlq) continue;
00196 if (te == planeelementlq && Mt->give_elem_type(i) == planeelementlt) continue;
00197 print_err("Uniform type of elements is required, element %ld differs", __FILE__, __LINE__, __func__, i+1);
00198 exit (1);
00199 }
00200
00201
00202 switch (te){
00203 case planeelementlt:
00204 case planeelementlq:
00205 if (Pelt && (Pelt->nb!=1 || Pelt->nip[0][0]!=1 || Pelt->ncomp[0]!=3) ) { print_err("Pelt has nonstandard intpoints", __FILE__, __LINE__, __func__); exit (1); }
00206 if (Pelq && (Pelq->nb!=2 || Pelq->nip[0][0]!=4 || Pelq->nip[1][1]!=1 || Pelq->nip[1][0]!=0 || Pelq->nip[0][1]!=0 || Pelq->ncomp[0]!=2 || Pelq->ncomp[1]!=1)) { print_err("Pelq has nonstandard intpoints", __FILE__, __LINE__, __func__); exit (1); }
00207 break;
00208 case planeelementqt:
00209 case planeelementqq:
00210 if (Peqt && (Peqt->nb!=2 || Peqt->nip[0][0]!=3 || Peqt->nip[1][1]!=3 || Peqt->nip[1][0]!=0 || Peqt->nip[0][1]!=0 || Peqt->ncomp[0]!=2 || Peqt->ncomp[1]!=1)) { print_err("Peqt has nonstandard intpoints", __FILE__, __LINE__, __func__); exit (1); }
00211 if (Peqq && (Peqq->nb!=2 || Peqq->nip[0][0]!=4 || Peqq->nip[1][1]!=4 || Peqq->nip[1][0]!=0 || Peqq->nip[0][1]!=0 || Peqq->ncomp[0]!=2 || Peqq->ncomp[1]!=1)) { print_err("Peqq has nonstandard intpoints", __FILE__, __LINE__, __func__); exit (1); }
00212 break;
00213 case cctel:
00214 if (Cct && (Cct ->nb!=2 || Cct ->nip[0][0]!=1 || Cct ->nip[1][1]!=1 || Cct ->nip[1][0]!=0 || Cct ->nip[0][1]!=0 || Cct ->ncomp[0]!=3 || Cct ->ncomp[1]!=2)) { print_err("Cct has nonstandard intpoints", __FILE__, __LINE__, __func__); exit (1); }
00215 break;
00216 case lineartet:
00217 if (Ltet && (Ltet->nb!=1 || Ltet->nip[0][0]!=1 || Ltet->ncomp[0]!=6 ) ) { print_err("Ltet has nonstandard intpoints", __FILE__, __LINE__, __func__); exit (1); }
00218 break;
00219 default: print_err("unknown element type is required", __FILE__, __LINE__, __func__); exit (1);
00220 }
00221
00222
00223
00224
00225
00226
00227
00228 }
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 void adaptivity :: spr (void)
00239 {
00240
00241
00242
00243 matrix *spder = new matrix[ne];
00244
00245 int nip;
00246 long i, j, ipp;
00247 vector vect;
00248 vect.n = ncomp;
00249 for (i=0; i<ne; i++) {
00250
00251
00252 nip = 0;
00253 switch (Mt->elements[i].te) {
00254 case planeelementlt: nip = 1; break;
00255 case planeelementlq: { Pelq->elchar (i,spder[i]); break;}
00256 case planeelementqt: { Peqt->elchar (i,spder[i]); break;}
00257 case planeelementqq: { Peqq->elchar (i,spder[i]); break;}
00258 case lineartet: { Ltet->elchar (i,spder[i]); break;}
00259 default: print_err("unknown element type", __FILE__, __LINE__, __func__); exit (0);
00260 }
00261
00262 if (nip) {
00263 allocm(nip, ncomp, spder[i]);
00264 ipp = Mt->elements[i].ipp[0][0];
00265
00266 for (j=0; j<nip; j++) {
00267 vect.a = spder[i].a + j*ncomp;
00268
00269 if (!(otherflags & 2)) Mm->givestrain (0, ipp, vect);
00270 else Mm->givestress (0,ipp,vect);
00271 ipp++;
00272 }
00273 }
00274 }
00275 vect.a = NULL;
00276
00277
00278 patch_averaging spr (Gtm, dim, ncomp, 0);
00279
00280 spr.solve (Out, spder, refined_ders_in_nodes);
00281
00282 delete [] spder;
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292 void adaptivity :: compute_error (void)
00293 {
00294 long i;
00295 double e2 ,u2 ,volume ,corr1,corr2;
00296 vector ei2(ne),ui2(ne),sizel(ne);
00297
00298 e2 = u2 = volume = 0.0;
00299
00300 for (i=0; i<ne; i++) {
00301
00302 switch (Mt->give_elem_type (i)) {
00303 case planeelementlt: volume += Pelt->compute_error (i, ei2[i], ui2[i], sizel[i], refined_ders_in_nodes,otherflags); corr1 = 1.00; corr2 = 1.11; break;
00304 case planeelementlq: volume += Pelq->compute_error (i, ei2[i], ui2[i], sizel[i], refined_ders_in_nodes); corr1 = 1.00; corr2 = 1.20; break;
00305 case planeelementqt: volume += Peqt->compute_error (i, ei2[i], ui2[i], sizel[i], refined_ders_in_nodes); corr1 = 1.07; corr2 = 1.55; break;
00306 case planeelementqq: volume += Peqq->compute_error (i, ei2[i], ui2[i], sizel[i], refined_ders_in_nodes); corr1 = 1.00; corr2 = 1.60; break;
00307 case lineartet: volume += Ltet->compute_error (i, ei2[i], ui2[i], sizel[i], refined_ders_in_nodes); corr1 = 1.00; corr2 = 1.00; break;
00308 default: print_err("wrong mesh - element %ld", __FILE__, __LINE__, __func__, i+1); exit (1);
00309 }
00310
00311 if (tad == 1) ei2[i] *= corr1*corr1;
00312 if (tad == 2) ei2[i] *= corr2*corr2;
00313
00314 e2 += ei2[i];
00315 u2 += ui2[i];
00316 }
00317 u2 += e2;
00318
00319
00320 compute_refsizel_lin (sizel.a, ei2.a, e2, u2);
00321
00322
00323 double er = sqrt(e2/u2);
00324 if (er > adapt_accuracy) answer = 1;
00325
00326 if (Mespr==1)
00327 fprintf (stdout, "\n global error = %.1lf%% (sqrt(e2/u2) == sqrt(%.3e / %.3e)) at volume = %f", 100*er, e2, u2, volume);
00328
00329
00330 if (printflags & 2 && answer){
00331 char name[255];
00332 sprintf (name,"%s%s.T3d.bgm%s", path, filename, suffix);
00333 print_valel (Out, Gtm, path, name, "chyba", refsizel.a, 'e');
00334 }
00335
00336
00337 allocv (ne, elem_error_pct);
00338 for (i=0; i<ne; i++)
00339 elem_error_pct[i] = 100.0 * sqrt( ei2[i] / (ui2[i]+ei2[i]) );
00340
00341 }
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352 void adaptivity :: compute_refsizel_lin (const double *sizel, const double *ei2, double e2, double u2)
00353 {
00354 long i;
00355 double em2,specacc;
00356 vector epsil;
00357
00358 allocv (ne, epsil);
00359
00360 specacc = this->adapt_accuracy;
00361 em2 = specacc*specacc * u2/ne;
00362
00363
00364
00365 allocv (ne, refsizel);
00366 allocv (ne, refsizelrel);
00367
00368 for (i=0; i<ne; i++) {
00369 epsil[i] = sqrt(ei2[i] / em2);
00370 if (epsil[i]<0.05) epsil[i] = 0.05;
00371
00372 refsizelrel[i] = 1.0 / pow(epsil[i] , 1.0/(double)ord * corr);
00373 refsizel[i] = sizel[i] / pow(epsil[i] , 1.0/(double)ord * corr);
00374 }
00375
00376
00377
00378
00379
00380
00381 }
00382
00383
00384
00385
00386
00387
00388
00389 void adaptivity :: print_test (void) const
00390 {
00391 long i,j,ipp;
00392 char f_name[255];
00393 FILE *f_test;
00394
00395 sprintf (f_name,"%s%s.sifel.test",path,filename);
00396 f_test = fopen (f_name,"a"); if (f_test==NULL) fprintf (stderr,"\n .test file has not been opened.\n\n");
00397
00398
00399
00400 fprintf (f_test,"\n\n LEFT HAND SIDE \n\n");
00401
00402 for (i=0;i<Mb->nlc*Ndofm;i++){
00403 fprintf (f_test," %20.10f ",Lsrs->lhs[i]);
00404 if ((i+1)%7==0) fprintf (f_test,"\n");
00405 }
00406
00407
00408
00409 fprintf (f_test,"\n\n STRAINS AND STRESSES ON INTEGRATION POINTS \n\n");
00410
00411 for (i=0;i<ne;i++){
00412 ipp = Mt->elements[i].ipp[0][0];
00413
00414 switch ((long)Mt->elements[i].te){
00415 case planeelementlq:{
00416 for (j=0;j<4;j++){
00417 fprintf (f_test," %20.10f " ,Mm->ip[ipp+j].strain[0]);
00418 fprintf (f_test," %20.10f " ,Mm->ip[ipp+j].strain[1]);
00419 fprintf (f_test," %20.10f \n",Mm->ip[ipp+4].strain[2]);
00420 }
00421 break;
00422 }
00423 case planeelementqq:{
00424 for (j=0;j<4;j++){
00425 fprintf (f_test," %20.10f " ,Mm->ip[ipp+j ].strain[0]);
00426 fprintf (f_test," %20.10f " ,Mm->ip[ipp+j ].strain[1]);
00427 fprintf (f_test," %20.10f \n",Mm->ip[ipp+j+4].strain[2]);
00428 }
00429 break;
00430 }
00431 case planeelementlt:{
00432 fprintf (f_test," %20.10f " ,Mm->ip[ipp].strain[0]);
00433 fprintf (f_test," %20.10f " ,Mm->ip[ipp].strain[1]);
00434 fprintf (f_test," %20.10f \n",Mm->ip[ipp].strain[2]);
00435 break;
00436 }
00437 case planeelementqt:
00438 case planeelementsubqt:{
00439 for (j=0;j<3;j++){
00440 fprintf (f_test," %20.10f " ,Mm->ip[ipp+j ].strain[0]);
00441 fprintf (f_test," %20.10f " ,Mm->ip[ipp+j ].strain[1]);
00442 fprintf (f_test," %20.10f \n",Mm->ip[ipp+j+3].strain[2]);
00443 }
00444 break;
00445 }
00446 case lineartet:{
00447 for (j=0;j<6;j++)
00448 fprintf (f_test," %20.10f ",Mm->ip[ipp].strain[j]);
00449
00450 fprintf (f_test,"\n");
00451 break;
00452 }
00453 default:{
00454 fprintf (stderr,"\n\n unknown element type is required in function");
00455 fprintf (stderr," adaptivity::print_test (%s, line %d)",__FILE__,__LINE__);
00456 }
00457 }
00458
00459 }
00460
00461
00462
00463 fprintf (f_test,"\n\n REFINED_DERS_IN_NODES \n\n");
00464
00465 long l=0;
00466 for (i=0; i<ncomp; i++)
00467 for (j=0; j<nn; j++) {
00468 fprintf (f_test," %20.10f ", refined_ders_in_nodes[j].a[i]);
00469 if ((l+1)%7==0) fprintf (f_test,"\n");
00470 l++;
00471 }
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483 fprintf (f_test,"\n\n REFSIZEL \n\n");
00484
00485 for (i=0;i<Mt->ne;i++){
00486 fprintf (f_test," %20.10f ",refsizel[i]);
00487 if ((i+1)%7==0) fprintf (f_test,"\n");
00488 }
00489
00490
00491 if (Mm->ip[0].eqother!=NULL){
00492 fprintf (f_test,"\n\n OTHER \n\n");
00493
00494 for (i=0;i<Mm->tnip;i++){
00495 for (j=0;j<6;j++)
00496 fprintf (f_test," %20.10f " ,Mm->ip[i].eqother[j]);
00497
00498 fprintf (f_test,"\n");
00499 }
00500 }
00501
00502 fclose (f_test);
00503 }