00001 #include "adaptivityt.h"
00002
00003 #include "globalt.h"
00004 #include "globmatt.h"
00005 #include "elemswitcht.h"
00006
00007
00008
00009 adaptivityt :: adaptivityt ()
00010 {
00011 dim = ord = ncomp = nn = ne = ntm = 0;
00012
00013 tad = Tp->adaptivityflag;
00014 printflags = 3;
00015 corr = 1.0;
00016 adapt_accuracy = 1.0;
00017
00018
00019 path = filename = suffix = NULL;
00020 interadaloop = false;
00021 step = 0;
00022 niwidth = 2;
00023 ni = NULL;
00024
00025 refined_ders_in_nodes = elem_error_pct = refsizelrel = refsizelabs = NULL;
00026
00027 r = rdr = NULL;
00028 tctrl = NULL;
00029
00030 answer = 0;
00031 }
00032
00033 adaptivityt :: ~adaptivityt ()
00034 {
00035 delete [] ni;
00036 delete [] path;
00037 delete [] filename;
00038 delete [] suffix;
00039 delete [] r;
00040 delete [] rdr;
00041 delete tctrl;
00042 }
00043
00044
00045 void adaptivityt :: set_step (long s)
00046 {
00047 step = s;
00048 if (ni == NULL) ni = new char[1+niwidth+1];
00049 if (step >= 0) sprintf (ni, ".%0*ld", niwidth, step);
00050 else ni[0] = '\0';
00051 }
00052
00053
00054 void adaptivityt :: readinit (XFILE *in)
00055 {
00056 xfscanf (in, "%k%lf %k%lf %k%lf", "accuracy", &adapt_accuracy, "enlargement", &enlarg, "reduction", &reduct);
00057
00058 if (Tp->adaptivityflag<1 || Tp->adaptivityflag>2)
00059 { print_err("unknown type of adaptivity is required",__FILE__,__LINE__,__func__); exit (0); }
00060
00061 filename_decomposition (in->fname, path, filename, suffix);
00062 }
00063
00064
00065 void adaptivityt :: printinit (FILE *out)
00066 {
00067 fprintf (out,"%ld %.3lf %lf %.2lf\n", tad, adapt_accuracy, enlarg, reduct);
00068 }
00069
00070
00071 void adaptivityt :: initialize (long s)
00072 {
00073
00074 dim = Tt->give_dimension (0);
00075 ord = Tt->give_degree (0);
00076 ncomp = Tt->give_ncomp (0);
00077 nn = Tt->nn;
00078 ne = Tt->ne;
00079 ntm = Tp->ntm;
00080
00081 for (long i=0; i<ne; i++)
00082 Gtt->gelements[i].auxinf = 100*Tt->give_nne(i) + 10*ord + dim;
00083
00084 answer = 0;
00085
00086 step = s;
00087
00088 prepare_ni ();
00089
00090
00091 this->check_consistency ();
00092 }
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 long adaptivityt :: run (long of, bool ial)
00105 {
00106 long hod,min;
00107 double sec;
00108
00109 if (Mesprt>0){
00110 fprintf (stdout,"\n\n *********************************");
00111 fprintf (stdout, "\n *** Adaptivity is started ***");
00112
00113 sec = clock();
00114 }
00115
00116 interadaloop = ial;
00117
00118
00119
00120 compute_ipgrads ();
00121
00122
00123
00124 if (dim == 0)
00125 this->initialize(0);
00126
00127
00128 refined_ders_in_nodes = new vector[nn];
00129 for (long i=0; i<nn; i++) allocv (ncomp*ntm, refined_ders_in_nodes[i]);
00130
00131 elem_error_pct = new vector[ntm+1];
00132 refsizelrel = new vector[ntm+1];
00133 refsizelabs = new vector[ntm+1];
00134
00135
00136 for (int i=0; i<ntm; i++) {
00137
00138 if (ntm==2 && i==1)
00139 for (long j=0; j<nn; j++)
00140 for (long k=0; k<ncomp; k++)
00141 refined_ders_in_nodes[j].a[k+i*ncomp] = refined_ders_in_nodes[j].a[k];
00142
00143 switch (tad){
00144 case 1:{
00145 spr (i);
00146 compute_error (i);
00147 break;
00148 }
00149 default: print_err("unknown type of error smoothing is required in function", __FILE__, __LINE__, __func__); exit (1);
00150 }
00151
00152
00153 }
00154
00155
00156 if (ntm > 1) {
00157 allocv (ne, elem_error_pct[ntm]);
00158 allocv (ne, refsizelrel[ntm]);
00159 allocv (ne, refsizelabs[ntm]);
00160
00161 for (long i=0; i<ne; i++) {
00162 elem_error_pct[ntm][i] = elem_error_pct[0][i];
00163 refsizelrel [ntm][i] = refsizelrel [0][i];
00164 refsizelabs [ntm][i] = refsizelabs [0][i];
00165
00166 if (ntm > 1) {
00167 if (elem_error_pct[ntm][i] < elem_error_pct[1][i]) elem_error_pct[ntm][i] = elem_error_pct[1][i];
00168 if (refsizelrel [ntm][i] > refsizelrel [1][i]) refsizelrel [ntm][i] = refsizelrel [1][i];
00169 if (refsizelabs [ntm][i] > refsizelabs [1][i]) refsizelabs [ntm][i] = refsizelabs [1][i];
00170 }
00171
00172 if (ntm > 2) {
00173 print_err("number of transported matters greater then 2", __FILE__, __LINE__, __func__); exit (1);
00174 }
00175 }
00176 }
00177
00178
00179 if (printflags & 2 && answer){
00180 char name[255];
00181 sprintf (name,"%s%s.T3d.bgm%s", path, filename, ni);
00182 print_valel (Outt, Gtt, path, name, "chyba", refsizelabs[ntm==1 ? 0 : ntm].a, 'e');
00183 }
00184
00185 this->print_addat_vtk ();
00186
00187
00188 delete [] refined_ders_in_nodes;
00189 delete [] elem_error_pct;
00190 delete [] refsizelrel;
00191 delete [] refsizelabs;
00192
00193 if (Mesprt>0){
00194 sec = (clock() - sec) / (double)CLOCKS_PER_SEC;
00195 hod = (long)sec/3600; sec -= hod*3600;
00196 min = (long)sec/60; sec -= min*60;
00197 fprintf(stdout,"\n Consumed time by Ada %2ld:%02ld:%05.2f",hod,min,sec);
00198
00199 fprintf (stdout,"\n *** Adaptivity is finished ***");
00200 fprintf (stdout,"\n *********************************");
00201 }
00202 return (answer);
00203 }
00204
00205
00206
00207 void adaptivityt :: prepare_ni (void)
00208 {
00209 if (! interadaloop) {
00210 step = atol(suffix+1);
00211
00212 if ( step ||
00213 (suffix[1] == '0' && suffix[2] == '\0') ||
00214 (suffix[1] == '0' && suffix[2] == '0' && suffix[3] == '\0') ||
00215 (suffix[1] == '0' && suffix[2] == '0' && suffix[3] == '0' && suffix[4] == '\0') )
00216 ;
00217 else
00218 step = -1;
00219 }
00220
00221 this->set_step(this->step);
00222
00223
00224 if (strncmp(filename, "this.sifel", 10) == 0)
00225 filename[4] = '\0';
00226 }
00227
00228
00229
00230
00231
00232
00233 void adaptivityt :: check_consistency (void) const
00234 {
00235 int te = Tt->give_elem_type (0);
00236 for (long i=0; i<nn; i++) {
00237 if (ntm != Tt->give_ndofn (i)) { print_err("mismatch in ndofn, node %ld differs", __FILE__, __LINE__, __func__, i+1); exit (1); }
00238 }
00239 for (long i=0; i<ne; i++) {
00240 if (te != Tt->give_elem_type(i)) { print_err("Uniform type of elements is required, element %ld differs", __FILE__, __LINE__, __func__, i+1); exit (1); }
00241 if (ord != Tt->give_degree (i)) { print_err("mismatch in order", __FILE__, __LINE__, __func__); exit (1); }
00242 if (dim != Tt->give_dimension(i)) { print_err("mismatch in dimension", __FILE__, __LINE__, __func__); exit (1); }
00243 }
00244
00245 switch (te){
00246 case trlint:
00247 case quadlint:
00248 if (Ltt && (Ltt->intordkm[0][0]!=1 || Ltt->ncomp!=2)) { print_err("Ltt has nonstandard intpoints", __FILE__, __LINE__, __func__); exit (1); }
00249 if (Lqt && (Lqt->intordkm[0][0]!=2 || Lqt->ncomp!=2)) { print_err("Lqt has nonstandard intpoints", __FILE__, __LINE__, __func__); exit (1); }
00250 break;
00251 default: print_err("unknown element type is required", __FILE__, __LINE__, __func__); exit (1);
00252 }
00253
00254
00255
00256
00257
00258
00259
00260 }
00261
00262
00263
00264
00265
00266
00267
00268
00269 void adaptivityt :: spr (int mattid)
00270 {
00271
00272
00273
00274 matrix *spder = new matrix[ne];
00275
00276
00277 if (Tp->savemode!=1) { print_err("savemode is NOT TRUE", __FILE__, __LINE__, __func__); exit (0); }
00278
00279 int nip;
00280 long i, j, ipp;
00281 vector vect;
00282 vect.n = ncomp;
00283 for (i=0; i<ne; i++) {
00284
00285
00286 switch (Tt->elements[i].te) {
00287 case quadlint: nip = 4; break;
00288 case trlint: nip = 1; break;
00289 default: print_err("unknown element type", __FILE__, __LINE__, __func__); exit (0);
00290 }
00291
00292 allocm(nip, ncomp, spder[i]);
00293 ipp = Tt->elements[i].ipp[0][0];
00294
00295 for (j=0; j<nip; j++) {
00296 vect.a = spder[i].a + j*ncomp;
00297 Tm->givegrad (mattid, ipp, vect);
00298 ipp++;
00299 }
00300 }
00301 vect.a = NULL;
00302
00303
00304 patch_averaging spr (Gtt, dim, ncomp, 0);
00305
00306 spr.solve (Outt, spder, refined_ders_in_nodes);
00307
00308 delete [] spder;
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318 void adaptivityt :: compute_error (int mattid)
00319 {
00320 long i;
00321 double e2 ,u2 ,volume ;
00322 vector ei2(ne),ui2(ne),sizel(ne);
00323
00324 e2 = u2 = volume = 0.0;
00325
00326 for (i=0; i<ne; i++) {
00327
00328 switch (Tt->give_elem_type (i)) {
00329 case trlint: volume += Ltt->compute_error (i, refined_ders_in_nodes, mattid, ei2[i], ui2[i], sizel[i]); break;
00330 case quadlint: volume += Lqt->compute_error (i, refined_ders_in_nodes, mattid, ei2[i], ui2[i], sizel[i]); break;
00331 default: print_err("wrong mesh - element %ld", __FILE__, __LINE__, __func__, i+1); exit (1);
00332 }
00333
00334 e2 += ei2[i];
00335 u2 += ui2[i];
00336 }
00337 u2 += e2;
00338
00339
00340 compute_refsizel_lin (mattid, sizel.a, ei2.a, e2, u2);
00341
00342
00343 double er = sqrt(e2/u2);
00344 if (er > adapt_accuracy || er < 0.02) answer = 1;
00345
00346
00347 if (Mesprt>0)
00348 fprintf (stdout, "\n Matter %d: global error = %.1lf%% (sqrt(e2/u2) == sqrt(%.3e / %.3e)) at volume = %f", mattid+1, 100*er, e2, u2, volume);
00349
00350
00351 allocv (ne, elem_error_pct[mattid]);
00352 for (i=0; i<ne; i++)
00353 elem_error_pct[mattid][i] = 100.0 * sqrt( ei2[i] / (ui2[i]+ei2[i]) );
00354
00355 }
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366 void adaptivityt :: compute_refsizel_lin (int mattid, const double *sizel, const double *ei2, double e2, double u2)
00367 {
00368 long i;
00369 double em2,specacc;
00370 vector epsil;
00371
00372 allocv (ne, epsil);
00373
00374 specacc = this->adapt_accuracy;
00375 em2 = specacc*specacc * u2/ne;
00376
00377 corr = 1.0;
00378
00379 allocv (ne, refsizelrel[mattid]);
00380 allocv (ne, refsizelabs[mattid]);
00381
00382 for (i=0; i<ne; i++) {
00383 epsil[i] = sqrt(ei2[i] / em2);
00384 if (epsil[i] < 1.0 / enlarg) epsil[i] = 1.0 / enlarg;
00385 if (epsil[i] > 1.0 / reduct) epsil[i] = 1.0 / reduct;
00386
00387 refsizelrel[mattid][i] = 1.0 / pow(epsil[i] , 1.0/(double)ord * corr);
00388 refsizelabs[mattid][i] = sizel[i] * refsizelrel[mattid][i];
00389 }
00390
00391
00392
00393
00394
00395
00396 }
00397
00398
00399
00400 void adaptivityt :: print_addat_vtk () const
00401 {
00402 long i,j;
00403 char name[255];
00404 FILE *out;
00405
00406
00407 sprintf (name,"%s%s.addat.vtk%s", path, filename, ni);
00408 out = fopen (name,"w"); if (out==NULL) fprintf (stderr,"\n .test file has not been opened.\n\n");
00409
00410 fprintf (out,"#\n");
00411
00412 fprintf (out,"POINT_DATA %ld\n", nn);
00413
00414
00415 fprintf (out,"SCALARS ada_refined_derivatives float %ld\n", ncomp*ntm);
00416 fprintf (out,"LOOKUP_TABLE default\n");
00417
00418 for (i=0; i<nn; i++) {
00419 for (j=0; j<ncomp*ntm; j++)
00420 fprintf (out, " %13.6f", refined_ders_in_nodes[i].a[j]);
00421
00422 fprintf (out,"\n");
00423 }
00424
00425
00426
00427 fprintf (out,"CELL_DATA %ld\n", ne);
00428
00429
00430 fprintf (out,"SCALARS ada_error_pct float %ld\n", ntm==1 ? 1 : ntm + 1);
00431 fprintf (out,"LOOKUP_TABLE default\n");
00432
00433 for (i=0; i<ne; i++) {
00434 ; fprintf (out, "%13.6f", elem_error_pct[0] [i]);
00435 if (ntm > 1) fprintf (out, "%13.6f", elem_error_pct[1] [i]);
00436 if (ntm !=1) fprintf (out, "%13.6f", elem_error_pct[ntm][i]);
00437
00438 fprintf (out,"\n");
00439 }
00440
00441
00442 fprintf (out,"SCALARS ada_refsizel_relative float %ld\n", ntm==1 ? 1 : ntm + 1);
00443 fprintf (out,"LOOKUP_TABLE default\n");
00444
00445 for (i=0; i<ne; i++) {
00446 ; fprintf (out, "%13.6f", refsizelrel[0] [i]);
00447 if (ntm > 1) fprintf (out, "%13.6f", refsizelrel[1] [i]);
00448 if (ntm !=1) fprintf (out, "%13.6f", refsizelrel[ntm][i]);
00449
00450 fprintf (out,"\n");
00451 }
00452
00453 fclose (out);
00454 }
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474 long closest_node (gtopology *gt, long dim, double x, double y, double z)
00475 {
00476 long i, ans;
00477 double aux, dist=1e100;
00478 double coord[3];
00479
00480 coord[0] = x;
00481 coord[1] = y;
00482 coord[2] = z;
00483
00484 for (i=0; i<gt->nn; i++)
00485 if ( (aux=gt->gnodes[i].distance2(dim, coord)) < dist ) {
00486 dist = aux;
00487 ans = i;
00488 }
00489
00490 return ans;
00491 }
00492
00493
00494
00495
00496
00497
00498
00499 bool isnodonlhsofline (gnode &n1, gnode &n2, double xx, double yy)
00500 {
00501 if (n1.x != n2.x)
00502 if ( (n2.x>n1.x ? 1.0:-1.0)*((xx-n1.x)*(n1.y-n2.y)/(n1.x-n2.x)+(n1.y-yy)) <= 1.0e-13 ) return true;
00503 else return false;
00504 else
00505 if ( (n1.y-n2.y)*(xx-n1.x) > -1.0e-13 ) return true;
00506 else return false;
00507 }
00508
00509
00510
00511
00512
00513
00514
00515 bool isnodonlhsof3pcurve (gnode &n1, gnode &n2, gnode &n3, double xx, double yy)
00516 {
00517 double c;
00518
00519 c = (n1.x*(n3.y-n2.y)+n2.x*(n1.y-n3.y)+n3.x*(n2.y-n1.y))/(fabs(n1.x-n3.x)+fabs(n1.y-n3.y));
00520 if (-1.0e-13<c && c<1.0e-13)
00521 return (isnodonlhsofline (n1,n3,xx,yy));
00522 else {
00523
00524
00525
00526 double s,x[3],y[2];
00527
00528 c = sqrt((n3.x-n1.x)*(n3.x-n1.x)+(n3.y-n1.y)*(n3.y-n1.y));
00529 s = (n3.y-n1.y)/c;
00530 c = (n3.x-n1.x)/c;
00531
00532 x[0] = c*( xx-n1.x) + s*( yy-n1.y); y[0] = -s*( xx-n1.x) + c*( yy-n1.y);
00533 x[1] = c*(n2.x-n1.x) + s*(n2.y-n1.y); y[1] = -s*(n2.x-n1.x) + c*(n2.y-n1.y);
00534 x[2] = c*(n3.x-n1.x) + s*(n3.y-n1.y);
00535
00536 if ( y[1]*x[0]*x[0]/x[1]/(x[1]-x[2]) + y[1]*x[2]*x[0]/x[1]/(x[2]-x[1]) -y[0] <= 1.0e-13 )
00537 return true;
00538
00539 return false;
00540 }
00541 }
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556 long ispointinel (gnode *nod, long *elnod, double x, double y, long ndd)
00557 {
00558 switch (ndd) {
00559 case 312:{
00560 if ( isnodonlhsofline (nod[elnod[0]],nod[elnod[1]],x,y) &&
00561 isnodonlhsofline (nod[elnod[1]],nod[elnod[2]],x,y) &&
00562 isnodonlhsofline (nod[elnod[2]],nod[elnod[0]],x,y) )
00563 return true;
00564 break;
00565 }
00566 case 412:{
00567 if ( isnodonlhsofline (nod[elnod[0]],nod[elnod[1]],x,y) &&
00568 isnodonlhsofline (nod[elnod[1]],nod[elnod[2]],x,y) &&
00569 isnodonlhsofline (nod[elnod[2]],nod[elnod[3]],x,y) &&
00570 isnodonlhsofline (nod[elnod[3]],nod[elnod[0]],x,y) )
00571 return true;
00572 break;
00573 }
00574 case 622:{
00575 if ( isnodonlhsof3pcurve (nod[elnod[0]],nod[elnod[3]],nod[elnod[1]],x,y) &&
00576 isnodonlhsof3pcurve (nod[elnod[1]],nod[elnod[4]],nod[elnod[2]],x,y) &&
00577 isnodonlhsof3pcurve (nod[elnod[2]],nod[elnod[5]],nod[elnod[0]],x,y) )
00578 return true;
00579 break;
00580 }
00581 case 822:{
00582 if ( isnodonlhsof3pcurve (nod[elnod[0]],nod[elnod[4]],nod[elnod[1]],x,y) &&
00583 isnodonlhsof3pcurve (nod[elnod[1]],nod[elnod[5]],nod[elnod[2]],x,y) &&
00584 isnodonlhsof3pcurve (nod[elnod[2]],nod[elnod[6]],nod[elnod[3]],x,y) &&
00585 isnodonlhsof3pcurve (nod[elnod[3]],nod[elnod[7]],nod[elnod[0]],x,y) )
00586 return true;
00587 break;
00588 }
00589
00590 default:{
00591 print_err("wrong dimdegnne",__FILE__,__LINE__,__func__);
00592 exit (0);
00593 break;
00594 }}
00595
00596 return false;
00597 }
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613 long whereispoint (gtopology *gt, long nsusel, long *susel, long *newsusel, long *passedel, double x, double y, long rep)
00614 {
00615 long i,j,el,nnewsusel;
00616
00617
00618 if (rep==0)
00619 for (i=0;i<nsusel;i++)
00620 if ( !passedel[susel[i]]++ && (ispointinel (gt->gnodes,gt->gelements[susel[i]].nodes,x,y,gt->gelements[susel[i]].auxinf)) )
00621 return (susel[i]);
00622
00623 nnewsusel = 0;
00624 for (i=0; i<nsusel; i++) {
00625 for (j=0; j<gt->nadjelel[susel[i]]; j++) {
00626 el = gt->adjelel[susel[i]][j];
00627 if (!passedel[el]){
00628 if (ispointinel (gt->gnodes,gt->gelements[el].nodes,x,y,gt->gelements[el].auxinf))
00629 return (el);
00630
00631 passedel[el] = 1;
00632 newsusel[nnewsusel++] = el;
00633 }
00634 }
00635 }
00636
00637 if (!nnewsusel || rep==5) {
00638 print_err("Point does not lay in domain",__FILE__,__LINE__,__func__);
00639 return (-1);
00640 }
00641
00642 return ( whereispoint (gt,nnewsusel,newsusel,susel,passedel,x,y,rep+1) );
00643 }
00644
00645 double distance2 (double *a,double *b,long dim)
00646 {
00647 if (dim == 1) return ((a[0]-b[0])*(a[0]-b[0]));
00648 else if (dim == 2) return ((a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]));
00649 else if (dim == 3) return ((a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]) + (a[2]-b[2])*(a[2]-b[2]));
00650 else { print_err("Wrong dimension",__FILE__,__LINE__,__func__); return (-1); }
00651 return 0;
00652 }
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664 long whereispoint_outpoint (gtopology *gt, long dim, double x, double y, double z)
00665 {
00666 long i,cl,nadjel,*adjel;
00667 double aux, centr[3], point[3],dist=1e100;
00668
00669 point[0]=x;
00670 point[1]=y;
00671 point[1]=z;
00672
00673
00674 cl = closest_node (gt, dim, x,y,z);
00675 nadjel = gt->nadjelnod[cl];
00676 adjel = gt->adjelnod[cl];
00677
00678 for (i=0;i<nadjel;i++)
00679 if ( ispointinel (gt->gnodes,gt->gelements[adjel[i]].nodes,x,y,gt->gelements[adjel[i]].auxinf) )
00680 return (adjel[i]);
00681
00682 for (i=0;i<nadjel;i++){
00683 gt->gelements[adjel[i]].centroid (dim,gt->gnodes,centr);
00684 aux = distance2 (point,centr,dim);
00685 if (aux<dist){
00686 dist=aux;
00687 cl=i;
00688 }
00689 }
00690
00691 return adjel[cl];
00692 }
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704 void findout_parentel_nod (gtopology *gt_old, gtopology *gt_new, long dim, long *parentel)
00705 {
00706 long i,j,k, n, nod, endheap;
00707 ivector oldsusel(gt_old->ne), newsusel(gt_old->ne), passedel(gt_old->ne), nodheap(gt_new->nn), nsusel(gt_new->nn);
00708 imatrix susel;
00709
00710 if (dim==2) allocm (gt_new->nn,18,susel);
00711 if (dim==3) allocm (gt_new->nn,120,susel);
00712
00713 fillv (-1, parentel, gt_new->nn);
00714
00715 nodheap[0] = 0;
00716 endheap = 1;
00717
00718
00719 nod = closest_node (gt_old, dim, gt_new->gnodes[nodheap[0]].x, gt_new->gnodes[nodheap[0]].y, gt_new->gnodes[nodheap[0]].z);
00720
00721 nsusel[nodheap[0]] = gt_old->nadjelnod[nod];
00722 for (i=0; i<nsusel[nodheap[0]]; i++)
00723 susel[nodheap[0]][i] = gt_old->adjelnod[nod][i];
00724
00725 for (i=0; i<gt_new->nn; i++) {
00726 nod = nodheap[i];
00727
00728 nullv (passedel);
00729 copyv (susel[nod],oldsusel.a,nsusel[nod]);
00730
00731 parentel[nod] = whereispoint (gt_old, nsusel[nod], oldsusel.a, newsusel.a, passedel.a, gt_new->gnodes[nod].x, gt_new->gnodes[nod].y, 0);
00732 if (parentel[nod]==-1)
00733 parentel[nod] = whereispoint_outpoint (gt_old, dim, gt_new->gnodes[nod].x, gt_new->gnodes[nod].y, gt_new->gnodes[nod].z);
00734
00735 for (j=0;j<gt_new->nadjelnod[nod];j++)
00736 for (k=0;k<gt_new->gelements[gt_new->adjelnod[nod][j]].nne;k++){
00737 n = gt_new->gelements[gt_new->adjelnod[nod][j]].nodes[k];
00738 if (parentel[n] == -1){
00739 parentel[n] = -2;
00740 nodheap[endheap++] = n;
00741 }
00742 if (parentel[n] == -2)
00743 susel[n][nsusel[n]++] = parentel[nod];
00744 }
00745
00746 }
00747 }
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765 void give_valuesinpoints (gtopology *gt, long eid, long npoints, double *xx, double *yy, long nval, const double *nodvalues, double **pointval, bool pvmtrx)
00766 {
00767 long i,j,k,nne;
00768 double xi,eta,value;
00769 vector x, y, bf, nodval;
00770 ivector nodes;
00771
00772 nne = gt->gelements[eid].nne;
00773
00774 allocv (nne,x);
00775 allocv (nne,y);
00776 allocv (nne,bf);
00777 allocv (nne,nodes);
00778 allocv (nne,nodval);
00779
00780 gt->give_nodes (eid,nodes);
00781 gt->give_node_coord2d (x,y,eid);
00782
00783 for (i=0; i<npoints; i++) {
00784
00785 switch (gt->gelements[eid].auxinf) {
00786 case 312: nc_lin_3_2d ( xx[i],yy[i],x.a,y.a, xi,eta); bf_lin_3_2d (bf.a,xi,eta); break;
00787 case 622: nc_quad_3_2d (1e-8,xx[i],yy[i],x.a,y.a, xi,eta); bf_quad_3_2d (bf.a,xi,eta); break;
00788 case 412: nc_lin_4_2d (1e-5,xx[i],yy[i],x.a,y.a, xi,eta); bf_lin_4_2d (bf.a,xi,eta); break;
00789 case 822: nc_quad_4_2d (1e-5,xx[i],yy[i],x.a,y.a, xi,eta); bf_quad_4_2d (bf.a,xi,eta); break;
00790 default: {
00791 print_err("unknown nnedegdim", __FILE__, __LINE__, __func__); exit (1);
00792 }}
00793
00794 for (j=0; j<nval; j++) {
00795 for (k=0; k<nne; k++)
00796 nodval[k] = nodvalues[nodes[k]+j*gt->nn];
00797
00798 scprd (bf,nodval,value);
00799
00800 if (pvmtrx) pointval[i][j] = value;
00801 else pointval[0][i+j*npoints] = value;
00802 }
00803 }
00804 }
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819 void transfvalues_nod (gtopology *gt_old, gtopology *gt_new, const double *r_old, double *r_new, long *parentel, long dim, long ndofn, bool spr)
00820 {
00821 long i,j,k,nid,maxnchild;
00822 long *nchilds,**childnod;
00823 double **spcoord, *r_loc;
00824 vector x,y;
00825
00826 if (dim != 2) { print_err("dim != 2", __FILE__, __LINE__, __func__); exit (1); }
00827
00828
00829 nchilds = new long [gt_old->ne];
00830 memset (nchilds, 0, gt_old->ne*sizeof(long));
00831
00832 for (i=0; i<gt_new->nn; i++)
00833 nchilds[parentel[i]]++;
00834
00835 childnod = new long* [gt_old->ne];
00836 for (i=0; i<gt_old->ne; i++)
00837 childnod[i] = new long [nchilds[i]];
00838
00839 memset (nchilds,0,gt_old->ne*sizeof(long));
00840 for (i=0; i<gt_new->nn; i++)
00841 childnod[parentel[i]][nchilds[parentel[i]]++] = i;
00842
00843
00844
00845 if (spr) {
00846 spcoord = new double* [gt_old->ne];
00847 for (i=0; i<gt_old->ne; i++) spcoord[i] = new double [nchilds[i]*dim];
00848 }
00849 else {
00850 maxnchild = 0;
00851 for (i=0; i<gt_old->ne; i++) if (maxnchild<nchilds[i]) maxnchild = nchilds[i];
00852 allocv (maxnchild, x);
00853 allocv (maxnchild, y);
00854 r_loc = new double [maxnchild*ndofn];
00855 }
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876 if (spr) {
00877 { print_err("nezkontrolovano", __FILE__, __LINE__, __func__); exit (1); }
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888 }
00889 else {
00890 for (i=0; i<gt_old->ne; i++) {
00891 for (j=0; j<nchilds[i]; j++) {
00892 x[j] = gt_new->gnodes[childnod[i][j]].x;
00893 y[j] = gt_new->gnodes[childnod[i][j]].y;
00894 }
00895
00896 give_valuesinpoints (gt_old, i, nchilds[i], x.a, y.a, ndofn, r_old, &r_loc, false);
00897
00898
00899 for (j=0; j<nchilds[i]; j++) {
00900 nid = childnod[i][j];
00901 for (k=0; k<ndofn; k++)
00902 r_new[nid+k*gt_new->nn] = r_loc[j+k*nchilds[i]];
00903 }
00904 }
00905 }
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952 delete [] r_loc;
00953 delete [] nchilds;
00954 destrm (childnod,gt_old->ne);
00955 if (spr) destrm (spcoord,gt_old->ne);
00956 }
00957
00958
00959
00960
00961 void adaptivityt :: statedata_backup (void)
00962 {
00963 if (Tb->nlc != ntm) { print_err("Tb->nlc != dim", __FILE__, __LINE__, __func__); exit (1); }
00964
00965
00966 delete [] r; this->r = new double [ntm * Gtt->nn];
00967 delete [] rdr; this->rdr = new double [ntm * Gtt->nn];
00968
00969 long i, j, ii;
00970 for (i=0; i<Gtt->nn; i++)
00971 for (j=0; j<ntm; j++) {
00972 ii = Gtt->give_dof(i,j);
00973 if (ii>0) {
00974 this->r [i+j*Gtt->nn] = Lsrst-> lhs[ii-1];
00975 this->rdr[i+j*Gtt->nn] = Lsrst->tdlhs[ii-1];
00976 }
00977 else {
00978 this->r [i+j*Gtt->nn] = 0.0;
00979 this->rdr[i+j*Gtt->nn] = 0.0;
00980 }
00981 }
00982
00983
00984 delete tctrl; this->tctrl = new timecontr; this->tctrl->be_copy_of (Tp->timecont);
00985 ; this->istep = Tp->istep;
00986 }
00987
00988
00989 void adaptivityt :: statedata_transfer (adaptivityt *Adat_old, gtopology *Gtt_old)
00990 {
00991
00992 long *parentel = new long [Gtt->nn];
00993 findout_parentel_nod (Gtt_old, Gtt, dim, parentel);
00994
00995 bool nod_spr = false;
00996
00997 r = new double [ntm * Gtt->nn];
00998 rdr = new double [ntm * Gtt->nn];
00999 transfvalues_nod (Gtt_old, Gtt, Adat_old->r, r, parentel, dim, ntm, nod_spr);
01000 transfvalues_nod (Gtt_old, Gtt, Adat_old->rdr, rdr, parentel, dim, ntm, nod_spr);
01001
01002 delete [] parentel;
01003
01004
01005
01006 if (tctrl) { print_err("tctrl is not NULL",__FILE__,__LINE__,__func__); exit (0); }
01007 this->tctrl = Adat_old->tctrl; Adat_old->tctrl = NULL;
01008
01009
01010 istep = Adat_old->istep;
01011 }
01012
01013
01014 void adaptivityt :: statedata_restore (void)
01015 {
01016 long i, j, ii;
01017 for (i=0; i<Gtt->nn; i++)
01018 for (j=0; j<ntm; j++) {
01019 ii = Gtt->give_dof(i,j);
01020 if (ii>0) {
01021 Lsrst-> lhs[ii-1] = this->r [i+j*Gtt->nn];
01022 Lsrst->tdlhs[ii-1] = this->rdr[i+j*Gtt->nn];
01023 }
01024 }
01025
01026
01027 Tp->istep = this->istep;
01028 }