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