00001 #include <string.h>
00002 #include "outdiagm.h"
00003 #include "iotools.h"
00004 #include "gtopology.h"
00005 #include "global.h"
00006 #include "probdesc.h"
00007 #include "meshtransfer.h"
00008 #include "mechtop.h"
00009 #include "element.h"
00010 #include "node.h"
00011 #include "intpoints.h"
00012 #include "globmat.h"
00013
00014
00015
00016
00017
00018
00019
00020
00021 outdiagm::outdiagm()
00022 {
00023 pid = NULL; eid = NULL; ipeid = NULL; nif = NULL;
00024 pu = NULL;
00025 ipu = NULL;
00026 }
00027
00028
00029
00030
00031
00032
00033
00034 outdiagm::~outdiagm()
00035 {
00036 delete [] nif;
00037 delete [] pid;
00038 delete [] eid;
00039 delete [] ipeid;
00040 delete [] pu;
00041 delete [] ipu;
00042 }
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 long outdiagm::read(XFILE *in)
00058 {
00059 long i;
00060 if (xfscanf (in, "%k%ld", "numunknowns", &npun) != 2)
00061 {
00062 print_err("cannot read number of printed unknowns", __FILE__,__LINE__, __func__);
00063 return 1;
00064 }
00065 xfscanf(in, "%k", "sel_diagstep");
00066 dstep.read(in);
00067 if (dstep.st == sel_no)
00068 return 0;
00069 nif = new nodip[npun];
00070 pid = new long [npun];
00071 eid = new long [npun];
00072 ipeid = new long [npun];
00073 pu = new prunk[npun];
00074 ipu = new long [npun];
00075 memset(nif, 0, sizeof(*nif)*npun);
00076 memset(pu, 0, sizeof(*pu)*npun);
00077 memset(ipu, 0, sizeof(*ipu)*npun);
00078 for (i=0; i<npun; i++)
00079 {
00080 pid[i] = eid[i] = ipeid[i] = -1;
00081 if (xfscanf (in, "%k%m", "point", &nodip_kwdset, (int *)(nif+i)) != 2)
00082 {
00083 print_err("cannot read type of point", __FILE__,__LINE__, __func__);
00084 return 1;
00085 }
00086 switch (nif[i])
00087 {
00088 case atnode:
00089 if (xfscanf (in, "%k%ld", "node", pid+i) != 2)
00090 {
00091 print_err("cannot read node id", __FILE__,__LINE__, __func__);
00092 return 1;
00093 }
00094 pid[i]--;
00095 break;
00096 case atip:
00097 if (xfscanf (in, "%k %ld %k %ld", "elem", eid+i, "ip", ipeid+i) != 4)
00098 {
00099 print_err("cannot read element or ip number", __FILE__,__LINE__, __func__);
00100 return 1;
00101 }
00102 eid[i]--; ipeid[i]--;
00103 break;
00104 case atxyz:
00105 if (xfscanf (in, "%k %le %k %le %k %le", "x", &x, "y", &y, "z", &z) != 6)
00106 {
00107 print_err("cannot read point coordinates", __FILE__,__LINE__, __func__);
00108 return 1;
00109 }
00110 pid[i] = adjacnode (Gtm,x,y,z);
00111 break;
00112 default:
00113 print_err("unknown type of point is required", __FILE__, __LINE__, __func__);
00114 return 1;
00115 }
00116
00117 if (xfscanf (in, "%k%m", "quant_type", &prunk_kwdset, (int *)(pu+i)) != 2)
00118 {
00119 print_err("cannot read type of printed unknown", __FILE__,__LINE__, __func__);
00120 return 1;
00121 }
00122 if ((pu[i] == pr_stepid) || (pu[i] == pr_appload) || (pu[i] == pr_time))
00123 continue;
00124 if (xfscanf (in, "%k%ld", "compid", ipu+i) != 2)
00125 {
00126 print_err("cannot read type of printed unknown", __FILE__,__LINE__, __func__);
00127 return 1;
00128 }
00129 ipu[i]--;
00130 if (nif[i] == atnode)
00131 {
00132 switch (pu[i])
00133 {
00134 case pr_strains:
00135 if(Mp->straincomp == 0)
00136 Mp->straincomp = 1;
00137 if(Mp->strainaver == 0)
00138 Mp->strainaver = 1;
00139 if(Mp->strainpos == 0)
00140 Mp->strainpos = 2;
00141 break;
00142 case pr_stresses:
00143 if(Mp->straincomp == 0)
00144 Mp->straincomp = 1;
00145 if(Mp->strainaver == 0)
00146 Mp->strainaver = 1;
00147 if(Mp->strainpos == 0)
00148 Mp->strainpos = 2;
00149 if(Mp->stresscomp == 0)
00150 Mp->stresscomp = 1;
00151 if(Mp->stressaver == 0)
00152 Mp->stressaver = 1;
00153 if(Mp->stresspos == 0)
00154 Mp->stresspos = 2;
00155 break;
00156 case pr_react:
00157 Mp->reactcomp = 1;
00158 break;
00159 case pr_other:
00160 if(Mp->othercomp == 0)
00161 Mp->othercomp = 1;
00162 if(Mp->otheraver == 0)
00163 Mp->otheraver = 1;
00164 if(Mp->otherpos == 0)
00165 Mp->otherpos = 2;
00166 break;
00167 default:
00168 break;
00169 }
00170 }
00171
00172 if ((nif[i] == atip) || (nif[i] == atxyz))
00173 {
00174 switch (pu[i])
00175 {
00176 case pr_strains:
00177 if(Mp->straincomp == 0)
00178 Mp->straincomp = 1;
00179 if(Mp->strainaver == 0)
00180 Mp->strainaver = 1;
00181 if(Mp->strainpos == 0)
00182 Mp->strainpos = 1;
00183 break;
00184 case pr_stresses:
00185 if(Mp->straincomp == 0)
00186 Mp->straincomp = 1;
00187 if(Mp->strainaver == 0)
00188 Mp->strainaver = 1;
00189 if(Mp->strainpos == 0)
00190 Mp->strainpos = 1;
00191 if(Mp->stresscomp == 0)
00192 Mp->stresscomp = 1;
00193 if(Mp->stressaver == 0)
00194 Mp->stressaver = 1;
00195 if(Mp->stresspos == 0)
00196 Mp->stresspos = 1;
00197 break;
00198 case pr_react:
00199 Mp->reactcomp = 1;
00200 break;
00201 case pr_other:
00202 if(Mp->othercomp == 0)
00203 Mp->othercomp = 1;
00204 if(Mp->otheraver == 0)
00205 Mp->otheraver = 1;
00206 if(Mp->otherpos == 0)
00207 Mp->otherpos = 1;
00208 break;
00209 default:
00210 break;
00211 }
00212 }
00213 }
00214 return 0;
00215 }
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 long outdiagm::print(FILE *out)
00229 {
00230 long i;
00231 fprintf(out, "%ld ", npun);
00232 dstep.print(out);
00233 if (dstep.st == sel_no)
00234 return 0;
00235
00236 for (i=0; i<npun; i++)
00237 {
00238 fprintf(out, "%d ", int(nif[i]));
00239 switch (nif[i])
00240 {
00241 case atnode:
00242 fprintf (out, "%ld ", pid[i]+1);
00243 break;
00244 case atip:
00245 fprintf (out, "%ld %ld ", eid[i]+1, ipeid[i]+1);
00246 break;
00247 case atxyz:
00248 fprintf (out, "%e %e %e ", x, y, z);
00249 break;
00250 default:
00251 print_err("unknown type of point is required", __FILE__, __LINE__, __func__);
00252 return 1;
00253 }
00254 fprintf (out, "%d", int(pu[i]));
00255
00256 if ((pu[i] == pr_stepid) || (pu[i] == pr_appload))
00257 {
00258 fprintf (out, "\n");
00259 continue;
00260 }
00261 fprintf (out, " %ld\n", ipu[i]+1);
00262 }
00263
00264 return 0;
00265 }
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279 long outdiagm::print_header(FILE *out)
00280 {
00281 long i;
00282 for (i=0; i<npun; i++)
00283 {
00284 switch (pu[i])
00285 {
00286 case pr_displ:
00287 fprintf(out, "r_%ld ", ipu[i]);
00288 break;
00289 case pr_strains:
00290 fprintf(out, "eps_%ld ", ipu[i]);
00291 break;
00292 case pr_stresses:
00293 fprintf(out, "sig_%ld ", ipu[i]);
00294 break;
00295 case pr_forces:
00296 fprintf(out, "f_%ld ", ipu[i]);
00297 break;
00298 case pr_react:
00299 fprintf(out, "R_%ld ", ipu[i]);
00300 break;
00301 case pr_stepid:
00302 fprintf(out, "step ");
00303 break;
00304 case pr_appload:
00305 fprintf(out, "lambda ");
00306 break;
00307 case pr_time:
00308 fprintf(out, "time ");
00309 break;
00310 case pr_other:
00311 fprintf(out, "other_%ld ", ipu[i]);
00312 break;
00313 default:
00314 print_err("unknown type of value for diagram is required", __FILE__, __LINE__, __func__, pid[i]+1);
00315 return 1;
00316 }
00317 }
00318 fprintf(out, "\n");
00319 return 0;
00320 }
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341 long outdiagm::printval(FILE *out, long lcid, double lambda, long istep, double *fi)
00342 {
00343 long i;
00344
00345
00346 if (dstep.presence_id(istep, lambda, Mp->timecon) == 0)
00347 return 0;
00348
00349 for (i=0; i<npun; i++)
00350 {
00351 if ((pid[i] < 0) && (nif[i] == atip))
00352 {
00353 if ((eid[i] < Mt->ne) && (ipeid[i] >= 0) && (ipeid[i] < Mt->give_tnip(eid[i])))
00354 pid[i] = Mt->elements[eid[i]].ipp[0][0]+ipeid[i];
00355 else
00356 {
00357 print_err("invalid element number or element ip number is required\n"
00358 "(element %ld, ip %ld)", __FILE__, __LINE__, __func__, eid[i]+1, ipeid[i]+1);
00359 return 1;
00360 }
00361 }
00362 switch (pu[i])
00363 {
00364 case pr_displ:
00365 print_displacements(out, lcid, i);
00366 break;
00367 case pr_strains:
00368 print_strains(out, lcid, i);
00369 break;
00370 case pr_stresses:
00371 print_stresses(out, lcid, i);
00372 break;
00373 case pr_macrostr:
00374 print_macrostr(out, lcid, i);
00375 break;
00376 case pr_forces:
00377 print_forces(out, i, fi);
00378 break;
00379 case pr_react:
00380 print_reactions(out, lcid, i);
00381 break;
00382 case pr_stepid :
00383 fprintf(out, "%ld ", istep);
00384 break;
00385 case pr_time:
00386 if(Mp->tpr == seconds)
00387 fprintf(out, "% e ", lambda);
00388 if(Mp->tpr == minutes)
00389 fprintf(out, "% e ", lambda/60.0);
00390 if(Mp->tpr == hours)
00391 fprintf(out, "% e ", lambda/3600.0);
00392 if(Mp->tpr == days)
00393 fprintf(out, "% e ", lambda/86400.0);
00394 break;
00395 case pr_appload:
00396 case pr_eigval:
00397 fprintf(out, "% e ", lambda);
00398 break;
00399 case pr_other:
00400 print_others(out, lcid, i);
00401 break;
00402 default:
00403 print_err("unknown type of value for diagram is required at node or ip number %ld",
00404 __FILE__, __LINE__, __func__, pid[i]+1);
00405 return 2;
00406 }
00407 }
00408 fprintf(out, "\n");
00409
00410 return 0;
00411 }
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 long outdiagm::printval_forced(FILE *out, long lcid, double lambda, long istep, double *fi)
00434 {
00435 long i;
00436
00437
00438 for (i=0; i<npun; i++)
00439 {
00440 if ((pid[i] < 0) && (nif[i] == atip))
00441 {
00442 if ((eid[i] < Mt->ne) && (ipeid[i] >= 0) && (ipeid[i] < Mt->give_tnip(eid[i])))
00443 pid[i] = Mt->elements[eid[i]].ipp[0][0]+ipeid[i];
00444 else
00445 {
00446 print_err("invalid element number or element ip number is required\n"
00447 "(element %ld, ip %ld)", __FILE__, __LINE__, __func__, eid[i]+1, ipeid[i]+1);
00448 return 1;
00449 }
00450 }
00451 switch (pu[i])
00452 {
00453 case pr_displ:
00454 print_displacements(out, lcid, i);
00455 break;
00456 case pr_strains:
00457 print_strains(out, lcid, i);
00458 break;
00459 case pr_stresses:
00460 print_stresses(out, lcid, i);
00461 break;
00462 case pr_forces:
00463 print_forces(out, i, fi);
00464 break;
00465 case pr_react:
00466 print_reactions(out, lcid, i);
00467 break;
00468 case pr_stepid:
00469 fprintf(out, "%ld ", istep);
00470 break;
00471 case pr_time:
00472 if(Mp->tpr == seconds)
00473 fprintf(out, "% e ", lambda);
00474 if(Mp->tpr == minutes)
00475 fprintf(out, "% e ", lambda/60.0);
00476 if(Mp->tpr == hours)
00477 fprintf(out, "% e ", lambda/3600.0);
00478 if(Mp->tpr == days)
00479 fprintf(out, "% e ", lambda/86400.0);
00480 break;
00481 case pr_appload:
00482 case pr_eigval:
00483 fprintf(out, "% e ", lambda);
00484 break;
00485 case pr_other:
00486 print_others(out, lcid, i);
00487 break;
00488 default:
00489 print_err("unknown type of value for diagram is required at node or ip number %ld",
00490 __FILE__, __LINE__, __func__, pid[i]+1);
00491 return 2;
00492 }
00493 }
00494 fprintf(out, "\n");
00495
00496 return 0;
00497 }
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514 long outdiagm::print_displacements(FILE *out, long lcid, long idp)
00515 {
00516 double *r;
00517 switch (nif[idp])
00518 {
00519 case atxyz:
00520 case atnode:
00521 r = new double [Gtm->gnodes[pid[idp]].ndofn];
00522 noddispl (lcid, r, pid[idp]);
00523 fprintf(out, "% .15e ", r[ipu[idp]]);
00524 delete [] r;
00525 break;
00526 case atip:
00527 print_err("unsupported combination of printed values is required -\n"
00528 " displacements at ip number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00529 return 1;
00530 default:
00531 print_err("unknown type of point is required (point number %ld)", __FILE__, __LINE__, __func__, pid[idp]+1);
00532 return 1;
00533 }
00534 return 0;
00535 }
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553 long outdiagm::print_strains(FILE *out, long lcid, long idp)
00554 {
00555 switch (nif[idp])
00556 {
00557 case atxyz:
00558 case atnode:
00559 if (ipu[idp] < Mt->nodes[pid[idp]].ncompstr)
00560 fprintf(out, "% .15e ", Mt->nodes[pid[idp]].strain[ipu[idp]]);
00561 else
00562 {
00563 print_err("invalid component number is required at node number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00564 return 2;
00565 }
00566 break;
00567 case atip:
00568 if (ipu[idp] < Mm->ip[pid[idp]].ncompstr)
00569 fprintf(out, "% .15e ", Mm->ip[pid[idp]].strain[ipu[idp]]);
00570 else
00571 {
00572 print_err("invalid component number is required at integration point number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00573 return 2;
00574 }
00575 break;
00576 default:
00577 print_err("unknown type of point is required (point number %ld)", __FILE__, __LINE__, __func__, pid[idp]+1);
00578 return 1;
00579 }
00580 return 0;
00581 }
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599 long outdiagm::print_stresses(FILE *out, long lcid, long idp)
00600 {
00601 switch (nif[idp])
00602 {
00603 case atxyz:
00604 case atnode:
00605 if (ipu[idp] < Mt->nodes[pid[idp]].ncompstr)
00606 fprintf(out, "% .15e ", Mt->nodes[pid[idp]].stress[ipu[idp]]);
00607 else
00608 {
00609 print_err("invalid component number is required at node number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00610 return 2;
00611 }
00612 break;
00613 case atip:
00614 if (ipu[idp] < Mm->ip[pid[idp]].ncompstr)
00615 fprintf(out, "% .15e ", Mm->ip[pid[idp]].stress[ipu[idp]]);
00616 else
00617 {
00618 print_err("invalid component number is required at integration point number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00619 return 2;
00620 }
00621 break;
00622 default:
00623 print_err("unknown type of point is required (point number %ld)", __FILE__, __LINE__, __func__, pid[idp]+1);
00624 return 1;
00625 }
00626 return 0;
00627 }
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646 long outdiagm::print_macrostr(FILE *out, long lcid, long idp)
00647 {
00648 if ((Mp->homog == 3) || (Mp->homog == 4))
00649 {
00650 long j,ncomp = Mt->give_ncomp(0, 0);
00651 double *lhs = Lsrs->give_lhs(lcid);
00652 fprintf(out, "** Calculated macro-strains for the homogenization problem:\n\n");
00653 j=Ndofm-ncomp;
00654 switch (nif[idp])
00655 {
00656 case atxyz:
00657 case atnode:
00658 if (ipu[idp] < ncomp)
00659 fprintf(out, "% .15e ", lhs[j+ipu[idp]]);
00660 else
00661 {
00662 print_err("invalid component number is required at node number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00663 return 2;
00664 }
00665 break;
00666 case atip:
00667 if (ipu[idp] < ncomp)
00668 fprintf(out, "% .15e ", lhs[j+ipu[idp]]);
00669 else
00670 {
00671 print_err("invalid component number is required at integration point number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00672 return 2;
00673 }
00674 break;
00675 default:
00676 print_err("unknown type of point is required (point number %ld)", __FILE__, __LINE__, __func__, pid[idp]+1);
00677 return 1;
00678 }
00679 }
00680 else
00681 {
00682 print_err("macro-stress/strain is required but they are not defined\n in the problem (point number %ld)",
00683 __FILE__, __LINE__, __func__, pid[idp]+1);
00684 return 3;
00685 }
00686
00687 return 0;
00688 }
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705 long outdiagm::print_forces(FILE *out, long idp, double *fi)
00706 {
00707 long ii;
00708 switch (nif[idp])
00709 {
00710 case atxyz:
00711 case atnode:
00712 ii=Mt->give_dof(pid[idp], ipu[idp]);
00713 if (ii > 0)
00714 fprintf(out, "% .15e ", fi[ii-1]);
00715 else
00716 {
00717 print_err("invalid component number is required at node number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00718 return 2;
00719 }
00720 break;
00721 case atip:
00722 print_err("unsupported combination of printed values is required -\n"
00723 " forces in ip number %ld",
00724 __FILE__, __LINE__, __func__, pid[idp]+1);
00725 return 1;
00726 default:
00727 print_err("unknown type of point is required (point number %ld)", __FILE__, __LINE__, __func__, pid[idp]+1);
00728 return 1;
00729 }
00730 return 0;
00731 }
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748 long outdiagm::print_reactions(FILE *out, long lcid, long idp)
00749 {
00750 switch (nif[idp])
00751 {
00752 case atxyz:
00753 case atnode:
00754 if (Mt->nodes[pid[idp]].react)
00755 fprintf(out, "% .15e ", Mt->nodes[pid[idp]].r[ipu[idp]]);
00756 break;
00757 case atip:
00758 print_err("unsupported combination of printed values is required -\n"
00759 " forces in ip number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00760 return 1;
00761 default:
00762 print_err("unknown type of point is required (point number %ld)", __FILE__, __LINE__, __func__, pid[idp]+1);
00763 return 1;
00764 }
00765 return 0;
00766 }
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784 long outdiagm::print_others(FILE *out, long lcid, long idp)
00785 {
00786 switch (nif[idp])
00787 {
00788 case atxyz:
00789 case atnode:
00790 if (ipu[idp] < Mt->nodes[pid[idp]].ncompother)
00791 fprintf(out, "% .15e ", Mt->nodes[pid[idp]].other[ipu[idp]]);
00792 else
00793 {
00794 print_err("invalid component number is required at node number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00795 return 2;
00796 }
00797 break;
00798 case atip:
00799 if (ipu[idp] < Mm->ip[pid[idp]].ncompeqother)
00800 fprintf(out, "% .15e ", Mm->ip[pid[idp]].eqother[ipu[idp]]);
00801 else
00802 {
00803 print_err("invalid component number is required at integration point number %ld", __FILE__, __LINE__, __func__, pid[idp]+1);
00804 return 2;
00805 }
00806 break;
00807 default:
00808 print_err("unknown type of point is required (point number %ld)", __FILE__, __LINE__, __func__, pid[idp]+1);
00809 return 1;
00810 }
00811 return 0;
00812 }