00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004
00005 #include "least_square.h"
00006 #include "gadaptivity.h"
00007 #include "vector.h"
00008 #include "matrix.h"
00009 #include "mathem.h"
00010 #include "gtopology.h"
00011 #include "basefun.h"
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 least_square::least_square (long dim,long nvals,gtopology *gt,long flag)
00027 {
00028 least_square::gt = gt;
00029 least_square::dim = dim;
00030 least_square::flag = flag;
00031 least_square::nvals = nvals;
00032 nn = gt->nn;
00033 ne = gt->ne;
00034
00035 nsp = NULL;
00036 spcoord = NULL;
00037 spvalue = NULL;
00038 nodvalue = NULL;
00039
00040 insidelem = new long [ne];
00041 insidenod = new long* [nn];
00042 maxcoord = new double* [ne];
00043
00044 nadjel = 0;
00045 adjel = NULL;
00046
00047 long i;
00048 for (i=0;i<nn;i++){
00049 insidenod[i] = new long [1];
00050 insidenod[i][0] = 1;
00051 }
00052 for (i=0;i<ne;i++)
00053 maxcoord[i] = new double[2*dim];
00054 }
00055
00056 least_square::~least_square (void)
00057 {
00058 if (insidenod != NULL) for (long i=0;i<nn;i++) delete [] insidenod[i];
00059 if (maxcoord != NULL) for (long i=0;i<ne;i++) delete [] maxcoord[i];
00060 if (spcoord != NULL) for (long i=0;i<ne;i++) delete [] spcoord[i];
00061
00062 delete [] nsp;
00063 delete [] insidenod;
00064 delete [] maxcoord;
00065 delete [] spcoord;
00066
00067 }
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 void least_square::spr_default (FILE *out,double **spvalue,double *nodvalue)
00083 {
00084 fprintf (stdout,"\n VVV=========== SPR_smoothing ===========VVV");
00085
00086 if (gt->nadjelnod == NULL)
00087 gt->adjacelem (out);
00088
00089 least_square::spvalue = spvalue;
00090 least_square::nodvalue = nodvalue;
00091
00092 nsp = new long [ne];
00093 spcoord = new double* [ne];
00094
00095 nsp_spcoord_maxcoord_assembling ('a');
00096 insidenod_assembling ();
00097 compute_patches_spr (0);
00098
00099 least_square::spvalue = NULL;
00100 least_square::nodvalue = NULL;
00101
00102 fprintf (stdout,"\n AAA=========== SPR_smoothing ===========AAA");
00103 }
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121 void least_square::spr_optional (FILE *out,long *nsp,double **spcoord,double **spvalue,double *nodvalue,long cut_flag)
00122 {
00123 fprintf (stdout,"\n VVV=========== SPR_smoothing ===========VVV");
00124
00125 if (gt->nadjelnod == NULL)
00126 gt->adjacelem (out);
00127
00128 least_square::nsp = nsp;
00129 least_square::spcoord = spcoord;
00130 least_square::spvalue = spvalue;
00131 least_square::nodvalue = nodvalue;
00132
00133 nsp_spcoord_maxcoord_assembling ('m');
00134 insidenod_assembling ();
00135 compute_patches_spr (cut_flag);
00136
00137 least_square::nsp = NULL;
00138 least_square::spcoord = NULL;
00139 least_square::spvalue = NULL;
00140 least_square::nodvalue = NULL;
00141
00142 fprintf (stdout,"\n AAA=========== SPR_smoothing ===========AAA");
00143 }
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161 void least_square::L2_nod2sp (FILE *out,long *nsp,double **spcoord,double *nodvalue,double **spvalue,char typ_patch)
00162 {
00163 fprintf (stdout,"\n VVV=========== L2_nod2sp ===========VVV");
00164
00165 if (gt->nadjelnod == NULL)
00166 gt->adjacelem (out);
00167
00168 least_square::nsp = nsp;
00169 least_square::spcoord = spcoord;
00170 least_square::spvalue = spvalue;
00171 least_square::nodvalue = nodvalue;
00172
00173 nsp_spcoord_maxcoord_assembling ('m');
00174 insidenod_assembling ();
00175 compute_patches_nod2sp (typ_patch);
00176
00177 least_square::nsp = NULL;
00178 least_square::spcoord = NULL;
00179 least_square::spvalue = NULL;
00180 least_square::nodvalue = NULL;
00181
00182 fprintf (stdout,"\n AAA=========== L2_nod2sp ===========AAA");
00183 }
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204 void least_square::L2_sp2sp (FILE *out,double **spvalue,long nap,const long *parentel,const double **apcoord,double **apvalue)
00205 {
00206
00207
00208
00209 fprintf (stdout,"\n ============== SPR_2_smoothing ==============");
00210
00211 if (gt->nadjelnod == NULL)
00212 gt->adjacelem (out);
00213
00214 least_square::spvalue = spvalue;
00215
00216 long *nadjap,**adjap;
00217
00218 nadjap = new long [nn];
00219 adjap = new long* [nn];
00220
00221 nsp = new long [ne];
00222 spcoord = new double* [ne];
00223
00224
00225 nsp_spcoord_maxcoord_assembling ('a');
00226 insidenod_assembling ();
00227 adjap_assembling (nap,parentel,apcoord,nadjap,adjap);
00228 compute_patches_sp2sp (nap,apcoord,nadjap,(const long**)adjap,apvalue);
00229
00230
00231 delete [] nadjap;
00232 for (long i=0;i<nn;i++) delete [] adjap[i];
00233 delete [] adjap;
00234
00235 least_square::spvalue = NULL;
00236 }
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249 void least_square::nsp_spcoord_maxcoord_assembling (char nsma_flag)
00250 {
00251 long i,j,k,nne;
00252 double xi,eta;
00253 vector x,y,z,bf;
00254
00255
00256 for (i=0;i<ne;i++){
00257 nne = gt->gelements[i].nne;
00258 if (dim==2) { allocv (nne,x); allocv (nne,y); gt->give_node_coord2d (x,y,i); }
00259 if (dim==3) { allocv (nne,x); allocv (nne,y); allocv (nne,z); gt->give_node_coord3d (x,y,z,i); }
00260
00261 if (nsma_flag == 'a'){
00262
00263 switch (gt->gelements[i].auxinf){
00264 case 312:{
00265 nsp[i] = 1;
00266 spcoord[i] = new double[2];
00267
00268 spcoord[i][0] = (x[0] + x[1] + x[2])/3.0;
00269 spcoord[i][1] = (y[0] + y[1] + y[2])/3.0;
00270 break;
00271 }
00272 case 412:{
00273 nsp[i] = 1;
00274 spcoord[i] = new double[2];
00275
00276 spcoord[i][0] = (x[0] + x[1] + x[2]+ x[3])/4.0;
00277 spcoord[i][1] = (y[0] + y[1] + y[2]+ y[3])/4.0;
00278 break;
00279 }
00280 case 622:{
00281 nsp[i] = 3;
00282 spcoord[i] = new double[6];
00283
00284 allocv (nne,bf);
00285 xi = eta = 1.0/6.0;
00286 for (j=0;j<3;j++){
00287 if (j==1) xi=4.0/6.0;
00288 if (j==2) {xi=1.0/6.0; eta=4.0/6.0;}
00289
00290 bf_quad_3_2d (bf.a,xi,eta);
00291 scprd (bf,x,spcoord[i][2*j ]);
00292 scprd (bf,y,spcoord[i][2*j+1]);
00293 }
00294 destrv (bf);
00295
00296 break;
00297 }
00298 case 822:{
00299 nsp[i] = 4;
00300 spcoord[i] = new double[8];
00301
00302 allocv (nne,bf);
00303 xi = eta = 0.577350269189626;
00304 for (j=0;j<2;j++){
00305 xi = -xi;
00306 for (k=0;k<2;k++){
00307 eta = -eta;
00308
00309 bf_quad_4_2d (bf.a,xi,eta);
00310 scprd (bf,x,spcoord[i][4*j+2*k ]);
00311 scprd (bf,y,spcoord[i][4*j+2*k+1]);
00312 }
00313 }
00314 destrv (bf);
00315
00316 break;
00317 }
00318 case 413:{
00319 nsp[i] = 1;
00320 spcoord[i] = new double[3];
00321
00322 spcoord[i][0] = (x[0] + x[1] + x[2] + x[3])/4.0;
00323 spcoord[i][1] = (y[0] + y[1] + y[2] + y[3])/4.0;
00324 spcoord[i][2] = (z[0] + z[1] + z[2] + z[3])/4.0;
00325
00326 break;
00327 }
00328 default:{
00329 fprintf (stderr,"\n\n wrong dimdegnne in function nsp_spcoord_maxcoord_assembling (%s, line %d)",__FILE__,__LINE__);
00330 break;
00331 }}
00332 }
00333
00334 if (nsma_flag == 'a' || nsma_flag == 'm'){
00335
00336 switch (gt->gelements[i].auxinf){
00337 case 312:
00338 case 622:{
00339 maxmin_3 (x.a,maxcoord[i][0],maxcoord[i][1]);
00340 maxmin_3 (y.a,maxcoord[i][2],maxcoord[i][3]);
00341 break;
00342 }
00343 case 412:
00344 case 822:{
00345 maxmin_4 (x.a,maxcoord[i][0],maxcoord[i][1]);
00346 maxmin_4 (y.a,maxcoord[i][2],maxcoord[i][3]);
00347 break;
00348 }
00349 case 413:{
00350 maxmin_4 (x.a,maxcoord[i][0],maxcoord[i][1]);
00351 maxmin_4 (y.a,maxcoord[i][2],maxcoord[i][3]);
00352 maxmin_4 (z.a,maxcoord[i][4],maxcoord[i][5]);
00353 break;
00354 }
00355 default:{
00356 fprintf (stderr,"\n\n wrong dimdegnne in function nsp_spcoord_maxcoord_assembling (%s, line %d)",__FILE__,__LINE__);
00357 break;
00358 }}
00359 }
00360
00361
00362 if (dim==2) { destrv (x); destrv (y); }
00363 if (dim==3) { destrv (x); destrv (y); destrv (z); }
00364 }
00365 }
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386 void least_square::insidenod_assembling (void)
00387 {
00388 long i,j,k,l,nne,*n,in[8];
00389 long jnod1,jnod2,sepnod,on,sum;
00390 long jnod,joinel,sepnodes[3];
00391 long aux3[3],n4,aux4[4];
00392
00393
00394
00395 for (i=0;i<ne;i++){
00396 nne = gt->give_nne (i);
00397 n = gt->gelements[i].nodes;
00398 for (j=0;j<nne;j++) in[j] = 0;
00399
00400 switch (gt->gelements[i].auxinf){
00401 case 312:{
00402 if (edge_position (gt,n[0],n[1]) == 1) {in[0]--; in[1]--;}
00403 if (edge_position (gt,n[1],n[2]) == 1) {in[1]--; in[2]--;}
00404 if (edge_position (gt,n[2],n[0]) == 1) {in[2]--; in[0]--;}
00405 break;
00406 }
00407 case 412:{
00408 if (edge_position (gt,n[0],n[1]) == 1) {in[0]--; in[1]--;}
00409 if (edge_position (gt,n[1],n[2]) == 1) {in[1]--; in[2]--;}
00410 if (edge_position (gt,n[2],n[3]) == 1) {in[2]--; in[3]--;}
00411 if (edge_position (gt,n[3],n[0]) == 1) {in[3]--; in[0]--;}
00412 break;
00413 }
00414 case 622:{
00415 if (gt->nadjelnod[n[3]] == 1) {in[0]--; in[1]--; in[3]--;}
00416 if (gt->nadjelnod[n[4]] == 1) {in[1]--; in[2]--; in[4]--;}
00417 if (gt->nadjelnod[n[5]] == 1) {in[2]--; in[0]--; in[5]--;}
00418 break;
00419 }
00420 case 822:{
00421 if (gt->nadjelnod[n[4]] == 1) {in[0]--; in[1]--; in[4]--;}
00422 if (gt->nadjelnod[n[5]] == 1) {in[1]--; in[2]--; in[5]--;}
00423 if (gt->nadjelnod[n[6]] == 1) {in[2]--; in[3]--; in[6]--;}
00424 if (gt->nadjelnod[n[7]] == 1) {in[3]--; in[0]--; in[7]--;}
00425 break;
00426 }
00427 case 413:{
00428 if (surface_position (gt,n[0],n[1],n[2]) == 1) {in[0]--; in[1]--; in[2]--;}
00429 if (surface_position (gt,n[0],n[2],n[3]) == 1) {in[0]--; in[2]--; in[3]--;}
00430 if (surface_position (gt,n[0],n[3],n[1]) == 1) {in[0]--; in[3]--; in[1]--;}
00431 if (surface_position (gt,n[1],n[2],n[3]) == 1) {in[1]--; in[2]--; in[3]--;}
00432 break;
00433 }
00434 default:{
00435 fprintf (stderr,"\n\n wrong dimdegnne in function insidenod_assembling (%s, line %d)",__FILE__,__LINE__);
00436 break;
00437 }
00438 }
00439
00440 for (j=0;j<nne;j++) if (in[j] < insidenod[n[j]][0] && in[j] != 0) insidenod[n[j]][0] = in[j];
00441 }
00442
00443
00444
00445 for (i=0;i<ne;i++){
00446 nne = gt->gelements[i].nne;
00447 n = gt->gelements[i].nodes;
00448
00449 insidelem[i] = 1;
00450 for (j=0;j<nne;j++)
00451 if (insidenod[n[j]][0] < 0){
00452 insidelem[i] = 0;
00453 break;
00454 }
00455
00456 switch (gt->gelements[i].auxinf){
00457 case 312:{
00458 if (insidenod[n[0]][0] == -1 && insidenod[n[1]][0] == -1 && insidenod[n[2]][0] == -1)
00459 fprintf (stderr,"\n\n !!!!!!! wrong mesh - element %ld (%s, line %d) !!!!!!!!! \n\n",i+1,__FILE__,__LINE__);
00460
00461 if (insidenod[n[0]][0] == -2 || insidenod[n[1]][0] == -2 || insidenod[n[2]][0] == -2){
00462 if (insidenod[n[0]][0] == -2) {jnod1 = n[1]; jnod2 = n[2]; sepnod = n[0];}
00463 else if (insidenod[n[1]][0] == -2) {jnod1 = n[0]; jnod2 = n[2]; sepnod = n[1];}
00464 else {jnod1 = n[0]; jnod2 = n[1]; sepnod = n[2];}
00465
00466 for (j=0;j<gt->nadjelnod[jnod1];j++)
00467 for (k=0;k<gt->nadjelnod[jnod2];k++)
00468 if ( gt->adjelnod[jnod1][j] == gt->adjelnod[jnod2][k] && gt->adjelnod[jnod2][k] != i){
00469 on = gt->adjelnod[jnod1][j];
00470 j = k = 9999999;
00471 }
00472
00473 if (gt->gelements[on].nodes[0] != jnod1 && gt->gelements[on].nodes[0] != jnod2) on = gt->gelements[on].nodes[0];
00474 else if (gt->gelements[on].nodes[1] != jnod1 && gt->gelements[on].nodes[1] != jnod2) on = gt->gelements[on].nodes[1];
00475 else on = gt->gelements[on].nodes[2];
00476
00477 if (insidenod[on][0] < 0)
00478 fprintf (stderr,"\n\n !!!!!!! wrong mesh - element %ld (%s, line %d) !!!!!!!!! \n\n",i+1,__FILE__,__LINE__);
00479 else
00480 if (insidenod[on][0] == 1){
00481 delete [] insidenod[on];
00482 insidenod[on] = new long [10];
00483 insidenod[on][0] = 2;
00484 }
00485 else
00486 insidenod[on][0]++;
00487
00488 insidenod[on][insidenod[on][0]-1] = sepnod;
00489 }
00490 break;
00491 }
00492 case 412:{
00493 sum = 0;
00494 for (j=0;j<nne;j++)
00495 if (insidenod[n[j]][0] < 0)
00496 sum += insidenod[n[j]][0];
00497
00498 if (sum < -4 || sum == -3 || (sum == -4 && insidenod[n[0]][0] == insidenod[n[1]][0]))
00499 fprintf (stderr,"\n\n !!!!!!! wrong mesh - element %ld (%s, line %d) !!!!!!!!! \n\n",i+1,__FILE__,__LINE__);
00500
00501 break;
00502 }
00503 case 622:{
00504 for (j=3;j<6;j++)
00505 if (insidenod[n[j]][0] == 1) insidenod[n[j]][0] = 0;
00506 else if (insidenod[n[j]][0] == -1) insidenod[n[j]][0] = -7;
00507
00508 if (insidenod[n[0]][0] == -1 && insidenod[n[1]][0] == -1 && insidenod[n[2]][0] == -1)
00509 fprintf (stderr,"\n\n !!!!!!! wrong mesh - element %ld (%s, line %d) !!!!!!!!! \n\n",i+1,__FILE__,__LINE__);
00510
00511 if (insidenod[n[0]][0] == -2 || insidenod[n[1]][0] == -2 || insidenod[n[2]][0] == -2){
00512 if (insidenod[n[4]][0] == 0) {jnod = n[4]; insidenod[n[4]][0] = -7; sepnodes[0] = n[0]; sepnodes[1] = n[3]; sepnodes[2] = n[5];}
00513 else if (insidenod[n[5]][0] == 0) {jnod = n[5]; insidenod[n[5]][0] = -7; sepnodes[0] = n[1]; sepnodes[1] = n[3]; sepnodes[2] = n[4];}
00514 else {jnod = n[3]; insidenod[n[3]][0] = -7; sepnodes[0] = n[2]; sepnodes[1] = n[4]; sepnodes[2] = n[5];}
00515
00516 joinel = (gt->adjelnod[jnod][0] == i) ? gt->adjelnod[jnod][1] : gt->adjelnod[jnod][0];
00517
00518 if (gt->gelements[joinel].nodes[4] == jnod) on = gt->gelements[joinel].nodes[0];
00519 else if (gt->gelements[joinel].nodes[5] == jnod) on = gt->gelements[joinel].nodes[1];
00520 else on = gt->gelements[joinel].nodes[2];
00521
00522 if (insidenod[on][0] < 0)
00523 fprintf (stderr,"\n\n !!!!!!! wrong mesh - element %ld (%s, line %d) !!!!!!!!! \n\n",i+1,__FILE__,__LINE__);
00524 else
00525 if (insidenod[on][0] == 1){
00526 delete [] insidenod[on];
00527 insidenod[on] = new long [30];
00528 insidenod[on][0] = 4;
00529 }
00530 else
00531 insidenod[on][0] += 3;
00532
00533 for (k=0,j=insidenod[on][0]-3;j<insidenod[on][0];k++,j++)
00534 insidenod[on][j] = sepnodes[k];
00535 }
00536 break;
00537 }
00538 case 822:{
00539 for (j=4;j<8;j++)
00540 if (insidenod[n[j]][0] == 1) insidenod[n[j]][0] = 0;
00541 else if (insidenod[n[j]][0] == -1) insidenod[n[j]][0] = -7;
00542
00543 sum = 0;
00544 for (j=0;j<4;j++)
00545 if (insidenod[n[j]][0] < 0)
00546 sum += insidenod[n[j]][0];
00547
00548 if (sum < -4 || sum == -3 || (sum == -4 && insidenod[n[0]][0] == insidenod[n[1]][0]))
00549 fprintf (stderr,"\n\n !!!!!!! wrong mesh - element %ld (%s, line %d) !!!!!!!!! \n\n",i+1,__FILE__,__LINE__);
00550
00551 break;
00552 }
00553 case 413:{
00554 if (insidenod[n[0]][0] == -3 || insidenod[n[1]][0] == -3 || insidenod[n[2]][0] == -3 || insidenod[n[3]][0] == -3){
00555 for (j=0,k=0;j<4;j++){
00556 if (insidenod[n[j]][0] != -3) {
00557 aux3[k] = n[j];
00558 k++;
00559 }
00560 else n4 = n[j];
00561 }
00562 on = opposite_node (gt,aux3,i);
00563
00564 if (insidenod[on][0] < 0)
00565 fprintf (stderr,"\n\n !!!!!!! wrong mesh - element %ld (%s, line %d) !!!!!!!!! \n\n",i+1,__FILE__,__LINE__);
00566 else
00567 if (insidenod[on][0] == 1){
00568 delete [] insidenod[on];
00569 insidenod[on] = new long [40];
00570 insidenod[on][0] = 2;
00571 insidenod[on][1] = n4;
00572 }
00573 else{
00574 insidenod[on][insidenod[on][0]] = n4;
00575 insidenod[on][0]++;
00576 }
00577 }
00578
00579 else if (insidenod[n[0]][0] != 0 && insidenod[n[1]][0] != 0 && insidenod[n[2]][0] != 0 && insidenod[n[3]][0] != 0){
00580 for (j=0;j<nne;j++)
00581 aux4[j] = -2;
00582
00583 for (j=0;j<nne;j++){
00584 for (l=0,k=0;k<nne;k++)
00585 {
00586 if (k!=j)
00587 aux3[l++] = n[k];
00588 }
00589
00590 on = opposite_node (gt,aux3,i);
00591
00592 if (on>=0)
00593 {
00594 if (insidenod[on][0]>0){
00595 for (l=0;l<nne;l++)
00596 if (l!=j) aux4[l] = -1;
00597 if (aux4[j] == -2)
00598 aux4[j] = on;
00599 }
00600 }
00601 }
00602
00603 for (j=0;j<nne;j++)
00604 {
00605 if (aux4[j] >= 0)
00606 {
00607 if (insidenod[aux4[j]][0] == 1){
00608 delete [] insidenod[aux4[j]];
00609 insidenod[aux4[j]] = new long [40];
00610 insidenod[aux4[j]][0] = 2;
00611 insidenod[aux4[j]][1] = n[j];
00612 }
00613 else {
00614 insidenod[aux4[j]][insidenod[aux4[j]][0]] = n[j];
00615 insidenod[aux4[j]][0]++;
00616 }
00617 }
00618 }
00619 }
00620 break;
00621 }
00622 default:{
00623 fprintf (stderr,"\n\n wrong dimdegnne in function insidenod_assembling (%s, line %d)",__FILE__,__LINE__);
00624 break;
00625 }
00626 }
00627 }
00628 n = NULL;
00629 }
00630
00631
00632
00633
00634
00635 void least_square::adjap_assembling (long nap,const long *parentel,const double **apcoord,long *nadjap,long **adjap)
00636 {
00637 long i,j,k,nne,*n,*aux;
00638 double dist,ldist;
00639
00640 aux = new long [nap];
00641
00642 for (i=0;i<nap;i++){
00643 nne = gt->gelements[parentel[i]].nne;
00644 n = gt->gelements[parentel[i]].nodes;
00645
00646 ldist = 1.0e100;
00647 for (j=0;j<nne;j++){
00648 if (insidenod[n[j]][0] != -1){
00649 dist = gt->gnodes[n[j]].distance2 (dim,apcoord[i]);
00650 if (dist < ldist){
00651 ldist = dist;
00652 aux[i] = n[j];
00653 }
00654 }
00655 }
00656 }
00657
00658 for (i=0;i<nn;i++)
00659 if (insidenod[i][0] > 1)
00660 for (j=1;j<insidenod[i][0];j++)
00661 for (k=0;k<nap;k++)
00662 if (aux[k] == insidenod[i][j])
00663 aux[k] = i;
00664
00665
00666 for (i=0;i<nap;i++)
00667 if (insidenod[aux[i]][0] < 1)
00668 fprintf (stderr,"\n\n something is wrong in function adjap_assembling\n\n\n");
00669
00670 memset (nadjap,0,nn*sizeof(long));
00671 for (i=0;i<nap;i++)
00672 nadjap[aux[i]]++;
00673
00674 for (i=0;i<nn;i++){
00675 adjap[i] = new long [nadjap[i]];
00676 nadjap[i] = 0;
00677 }
00678
00679 for (i=0;i<nap;i++)
00680 adjap[aux[i]][nadjap[aux[i]]++] = i;
00681
00682 }
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694 void least_square::compute_patches_spr (long cut_flag)
00695 {
00696 long i,j,k,eid,auxl,ncoef;
00697 double auxd;
00698 vector coef_normcoord(2*dim),*a;
00699 ivector magnitude(nn);
00700
00701 nullv (nodvalue,nn*nvals);
00702
00703 switch (gt->gelements[0].auxinf%100){
00704 case 12:{ ncoef = 3; break; }
00705 case 22:{ ncoef = 6; break; }
00706 case 13:{ ncoef = 4; break; }
00707 default:{
00708 fprintf (stderr,"\n\n unknown element type is required in function compute_patches_spr (%s, line %d)",__FILE__,__LINE__);
00709 break;
00710 }}
00711
00712 a = new vector[nvals];
00713 for (i=0;i<nvals;i++)
00714 allocv (ncoef,a[i]);
00715
00716 for (i=0;i<nn;i++){
00717 if (insidenod[i][0] < 1) continue;
00718
00719 nadjel = gt->nadjelnod[i];
00720 adjel = gt->adjelnod[i];
00721
00722 normal_coordinates_ae (coef_normcoord);
00723 polynom_coefficients_ae (coef_normcoord,a);
00724 nodvalue_assembling_ae (coef_normcoord,a,i,magnitude);
00725 }
00726
00727 for (i=0;i<nn;i++){
00728 if (magnitude[i] > 1)
00729 for (j=0;j<nvals;j++)
00730 nodvalue[i+j*nn] /= (double)magnitude[i];
00731 else if (magnitude[i] == 0)
00732 fprintf (stderr,"\n !!!!!!! stress of node %ld was not computed !!!!!!!!! \n",i+1);
00733 }
00734
00735
00736 if (cut_flag & 1){
00737 for (i=0;i<nn;i++){
00738 auxl=0;
00739 for (j=0;j<gt->nadjelnod[i];j++){
00740 eid=gt->adjelnod[i][j];
00741 auxd=0.0;
00742 for (k=0;k<nsp[eid]*nvals;k++)
00743 auxd+=spvalue[eid][k];
00744
00745 if (auxd==0.0) auxl++;
00746 }
00747 if ((gt->nadjelnod[i]-auxl)<=auxl)
00748 for (j=0;j<nvals;j++)
00749 nodvalue[i+j*nn]=0.0;
00750
00751 }
00752 }
00753
00754 for (i=0;i<nvals;i++) destrv (a[i]);
00755 delete [] a;
00756 }
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766 void least_square::compute_patches_nod2sp (char typ_patch)
00767 {
00768 long i,j,k;
00769 double **magnitude;
00770 vector coef_normcoord(2*dim),*a;
00771
00772
00773 magnitude = new double* [ne];
00774 for (i=0;i<ne;i++){
00775 magnitude[i] = new double [nsp[i]];
00776 memset (magnitude[i],0,nsp[i]*sizeof(double));
00777 memset (spvalue[i],0,nsp[i]*nvals*sizeof(double));
00778 }
00779
00780 a = new vector[nvals];
00781
00782
00783 if (typ_patch == 'n'){
00784 for (i=0;i<nn;i++){
00785 if (insidenod[i][0] < 1) continue;
00786
00787 nadjel = gt->nadjelnod[i];
00788 adjel = gt->adjelnod[i];
00789
00790 normal_coordinates_ae (coef_normcoord);
00791 polynom_coefficients_inv_ae (coef_normcoord,a,'n');
00792 spvalue_assembling_patch_ae (coef_normcoord,a,-1,magnitude);
00793 }
00794 }
00795 else if (typ_patch == 'e'){
00796 for (i=0;i<ne;i++){
00797 if (insidelem[i]==0) continue;
00798
00799 nadjel = gt->nadjelel[i];
00800 adjel = gt->adjelel[i];
00801
00802 normal_coordinates_ae (coef_normcoord);
00803 polynom_coefficients_inv_ae (coef_normcoord,a,'e');
00804 spvalue_assembling_patch_ae (coef_normcoord,a,i,magnitude);
00805 }
00806 }
00807 else fprintf (stderr,"\n\n typ_patch is wrong in function compute_patches_nod2sp (%s, line %d).\n",__FILE__,__LINE__);
00808
00809 for (i=0;i<ne;i++)
00810 for (j=0;j<nsp[i];j++){
00811 if (magnitude[i][j] == 0.0) fprintf (stderr,"\n !!!!!!! stress of element %ld was not computed !!!!!!!!! \n",i+1);
00812 for (k=0;k<nvals;k++)
00813 spvalue[i][j*nvals+k] /= magnitude[i][j];
00814 }
00815
00816
00817 destrm (magnitude,ne);
00818 for (i=0;i<nvals;i++) destrv (a[i]);
00819 delete [] a;
00820 }
00821
00822
00823
00824
00825
00826
00827
00828 void least_square::compute_patches_sp2sp (long nap,const double **apcoord,const long *nadjap,const long **adjap,double **apvalue)
00829 {
00830 long i,ncoef;
00831 vector coef_normcoord(2*dim),*a;
00832
00833 a = new vector[nvals];
00834 for (i=0;i<nvals;i++)
00835 allocv (3,a[i]);
00836
00837 for (i=0;i<nn;i++){
00838 if (insidenod[i][0] < 1) continue;
00839
00840 nadjel = gt->nadjelnod[i];
00841 adjel = gt->adjelnod[i];
00842
00843
00844
00845
00846 ncoef = 3;
00847
00848
00849
00850
00851 normal_coordinates_ae (coef_normcoord);
00852 polynom_coefficients_ae (coef_normcoord,a);
00853 apvalue_assembling (coef_normcoord,a,nadjap[i],adjap[i],apcoord,apvalue);
00854 }
00855
00856 for (i=0;i<nvals;i++) destrv (a[i]);
00857 delete [] a;
00858 }
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869 void least_square::normal_coordinates_ae (vector &coef_normcoord)
00870 {
00871 long i,j;
00872 ivector mm(2*dim);
00873
00874 fillv (adjel[0],mm);
00875
00876 for (i=1;i<nadjel;i++)
00877 for (j=0;j<2*dim;j+=2){
00878 if (maxcoord[adjel[i]][j ] > maxcoord[mm[j ]][j ]) mm[j ] = adjel[i];
00879 if (maxcoord[adjel[i]][j+1] < maxcoord[mm[j+1]][j+1]) mm[j+1] = adjel[i];
00880 }
00881
00882 for (i=0;i<2*dim;i+=2){
00883 coef_normcoord[i ] = (maxcoord[mm[i]][i] + maxcoord[mm[i+1]][i+1])/(maxcoord[mm[i+1]][i+1] - maxcoord[mm[i]][i]);
00884 coef_normcoord[i+1] = 2.0/(maxcoord[mm[i]][i] - maxcoord[mm[i+1]][i+1]);
00885 }
00886 }
00887
00888
00889
00890
00891
00892
00893
00894 void least_square::polynom_coefficients_ae (const vector &coef_normcoord,vector *a)
00895 {
00896 long i,j,k,l,m,ncoef=a[0].n;
00897 double zero=1.0e-20;
00898 vector normcoord(dim),pol(ncoef),contrv(ncoef);
00899 matrix ptp(ncoef,ncoef),contrm(ncoef,ncoef);
00900
00901 for (i=0;i<nvals;i++)
00902 fillv (0.0,a[i]);
00903
00904 for (i=0;i<nadjel;i++){
00905 for (j=0,l=0,m=0;j<nsp[adjel[i]];j++){
00906 for (k=0;k<dim;k++){
00907 normcoord[k] = coef_normcoord[2*k] + coef_normcoord[2*k+1] * spcoord[adjel[i]][l];
00908 l++;
00909 }
00910
00911 polynom (ncoef,normcoord.a,pol.a);
00912
00913 for (k=0;k<nvals;k++){
00914 cmulv (spvalue[adjel[i]][m],pol,contrv);
00915 addv (a[k],contrv,a[k]);
00916 m++;
00917 }
00918
00919 vxv (pol,pol,contrm);
00920 addm (ptp,contrm,ptp);
00921 }
00922 }
00923
00924 for (i=0;i<nvals;i++)
00925 gause (ptp,contrm,a[i],zero);
00926 }
00927
00928
00929
00930
00931
00932
00933
00934 void least_square::polynom_coefficients_inv_ae (const vector &coef_normcoord,vector *a,char typ_patch)
00935 {
00936 long i,j,k,nid,nne,ncoef,nadjn,adjanod[51],*nodes;
00937 double zero=1.0e-20;
00938
00939 nadjn = 0;
00940 for (i=0;i<nadjel;i++){
00941 nne = gt->gelements[adjel[i]].nne;
00942 nodes = gt->gelements[adjel[i]].nodes;
00943 for (j=0;j<nne;j++){
00944 for (k=0;k<nadjn;k++)
00945 if (nodes[j]==adjanod[k]) break;
00946
00947 if (k==nadjn)
00948 adjanod[nadjn++]=nodes[j];
00949
00950 }
00951 }
00952 nodes=NULL;
00953
00954
00955 if (typ_patch=='n'){
00956 if (nadjel>4) ncoef = 6;
00957 else ncoef = 3;
00958 }
00959 else if (typ_patch=='e'){
00960 if (nadjel>10) ncoef = 10;
00961 else ncoef = 6;
00962 }
00963 else fprintf (stderr,"\n\n typ_patch is wrong in function polynom_coefficients_inv_ae (%s, line %d).\n",__FILE__,__LINE__);
00964
00965 if (ncoef!=a[0].n)
00966 for (i=0;i<nvals;i++) { destrv (a[i]); allocv (ncoef,a[i]); }
00967 else
00968 for (i=0;i<nvals;i++) nullv (a[i].a,a[i].n);
00969
00970
00971 vector normcoord(dim),pol(ncoef),contrv(ncoef);
00972 matrix ptp(ncoef,ncoef),contrm(ncoef,ncoef);
00973
00974 for (i=0;i<nadjn;i++){
00975 nid = adjanod[i];
00976
00977 normcoord[0] = coef_normcoord[0] + coef_normcoord[1] * gt->gnodes[nid].x;
00978 normcoord[1] = coef_normcoord[2] + coef_normcoord[3] * gt->gnodes[nid].y;
00979 if (dim==3)
00980 normcoord[2] = coef_normcoord[4] + coef_normcoord[5] * gt->gnodes[nid].z;
00981
00982 polynom (ncoef,normcoord.a,pol.a);
00983
00984 for (j=0;j<nvals;j++){
00985 cmulv (nodvalue[nid+j*nn],pol,contrv);
00986 addv (a[j],contrv,a[j]);
00987 }
00988
00989 vxv (pol,pol,contrm);
00990 addm (ptp,contrm,ptp);
00991 }
00992
00993 for (i=0;i<nvals;i++)
00994 gause (ptp,contrm,a[i],zero);
00995 }
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007 void least_square::nodvalue_assembling_ae (const vector &coef_normcoord,const vector *a,long nid,ivector &magnitude)
01008 {
01009 long i,j,nne,*nodes,ncoef=a[0].n;
01010 ivector aux(nn);
01011
01012 sigma (ncoef,nid,coef_normcoord,a,magnitude);
01013
01014 for (i=0;i<nadjel;i++){
01015 nne = gt->give_nne (adjel[i]);
01016 nodes = gt->gelements[adjel[i]].nodes;
01017 for (j=0;j<nne;j++)
01018 aux[nodes[j]]++;
01019 }
01020 nodes=NULL;
01021
01022 for (i=1;i<insidenod[nid][0];i++)
01023 aux[insidenod[nid][i]] += 2;
01024
01025 for (i=0;i<nn;i++)
01026 if ((aux[i] > 1 && insidenod[i][0] < 1) || (aux[i] > 0 && insidenod[i][0] < -1))
01027 sigma (ncoef,i,coef_normcoord,a,magnitude);
01028 }
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040 void least_square::spvalue_assembling_patch_ae (const vector &coef_normcoord,const vector *a,long id,double **magnitude)
01041 {
01042 long i,j,k,eid;
01043 double contr,mag;
01044 vector pol(a[0].n),normcoord(dim);
01045
01046 for (i=0;i<nadjel;i++){
01047 eid = adjel[i];
01048 for (j=0;j<nsp[eid];j++){
01049
01050 if (eid == id) mag = 10;
01051 else mag = 1.0;
01052
01053 magnitude[eid][j] += mag;
01054
01055 for (k=0;k<dim;k++)
01056 normcoord[k] = coef_normcoord[2*k] + coef_normcoord[2*k+1] * spcoord[eid][k+j*dim];
01057
01058 polynom (a[0].n,normcoord.a,pol.a);
01059
01060 for (k=0;k<nvals;k++){
01061 scprd (pol,a[k],contr);
01062 spvalue[eid][j*nvals+k] += contr*mag;
01063 }
01064 }
01065 }
01066 }
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081 void least_square::apvalue_assembling (const vector &coef_normcoord,const vector *a,long nadjap,const long *adjap,const double **apcoord,double **apvalue)
01082 {
01083 long i,j,ncoef=a[0].n;
01084 vector pol(ncoef),normcoord(dim);
01085
01086 for (i=0;i<nadjap;i++){
01087 for (j=0;j<dim;j++)
01088 normcoord[j] = coef_normcoord[2*j] + coef_normcoord[2*j+1] * apcoord[adjap[i]][j];
01089
01090 polynom (ncoef,normcoord.a,pol.a);
01091
01092 for (j=0;j<nvals;j++)
01093 scprd (pol,a[j],apvalue[adjap[i]][j]);
01094
01095 }
01096 }
01097
01098
01099
01100
01101 void least_square::sigma (long ncoef,long nid,const vector &coef_normcoord,const vector *a,ivector &magnitude)
01102 {
01103 long i;
01104 double contr;
01105 vector pol(ncoef),normcoord(dim);
01106
01107 magnitude[nid]++;
01108
01109 normcoord[0] = coef_normcoord[0] + coef_normcoord[1] * gt->gnodes[nid].x;
01110 normcoord[1] = coef_normcoord[2] + coef_normcoord[3] * gt->gnodes[nid].y;
01111 if (dim == 3)
01112 normcoord[2] = coef_normcoord[4] + coef_normcoord[5] * gt->gnodes[nid].z;
01113
01114 polynom (ncoef,normcoord.a,pol.a);
01115
01116 for (i=0;i<nvals;i++){
01117 scprd (pol,a[i],contr);
01118 nodvalue[nid+i*nn] += contr;
01119 }
01120 }
01121
01122
01123
01124
01125
01126
01127
01128 void least_square::polynom (long ncoef,const double *normcoord,double *p)
01129 {
01130 if (dim==2) {
01131 if (ncoef<= 0) { fprintf (stderr,"\n\n wrong ncoef in function polynom (%s, line %d).\n",__FILE__,__LINE__); exit (1); }
01132 if (ncoef>= 1) p[ 0] = 1.0;
01133 if (ncoef>= 2) p[ 1] = normcoord[0];
01134 if (ncoef>= 3) p[ 2] = normcoord[1];
01135 if (ncoef>= 4) p[ 3] = normcoord[0] * normcoord[1];
01136 if (ncoef>= 5) p[ 4] = normcoord[0] * normcoord[0];
01137 if (ncoef>= 6) p[ 5] = normcoord[1] * normcoord[1];
01138 if (ncoef>= 7) p[ 6] = normcoord[0] * normcoord[0] * normcoord[1];
01139 if (ncoef>= 8) p[ 7] = normcoord[0] * normcoord[1] * normcoord[1];
01140 if (ncoef>= 9) p[ 8] = normcoord[0] * normcoord[0] * normcoord[0];
01141 if (ncoef>=10) p[ 9] = normcoord[1] * normcoord[1] * normcoord[1];
01142 if (ncoef>=11) p[10] = normcoord[0] * normcoord[0] * normcoord[0] * normcoord[1];
01143 if (ncoef>=12) p[11] = normcoord[0] * normcoord[0] * normcoord[1] * normcoord[1];
01144 if (ncoef>=13) p[12] = normcoord[0] * normcoord[1] * normcoord[1] * normcoord[1];
01145 if (ncoef>=14) p[13] = normcoord[0] * normcoord[0] * normcoord[0] * normcoord[0];
01146 if (ncoef>=15) p[14] = normcoord[1] * normcoord[1] * normcoord[1] * normcoord[1];
01147 if (ncoef >15) { fprintf (stderr,"\n\n wrong ncoef in function polynom (%s, line %d).\n",__FILE__,__LINE__); exit (1); }
01148 }
01149 else if (dim==3) {
01150 if (ncoef<= 0) { fprintf (stderr,"\n\n wrong ncoef in function polynom (%s, line %d).\n",__FILE__,__LINE__); exit (1); }
01151 if (ncoef>= 1) p[0] = 1.0;
01152 if (ncoef>= 2) p[1] = normcoord[0];
01153 if (ncoef>= 3) p[2] = normcoord[1];
01154 if (ncoef>= 4) p[3] = normcoord[2];
01155 if (ncoef > 4) { fprintf (stderr,"\n\n wrong ncoef in function polynom (%s, line %d).\n",__FILE__,__LINE__); exit (1); }
01156 }
01157 else { fprintf (stderr,"\n\n wrong dim in function polynom (%s, line %d).\n",__FILE__,__LINE__); exit (1); }
01158 }
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323