00001 #include <stdlib.h>
00002 #include <string.h>
00003
00004 #include "mechbclc.h"
00005 #include "global.h"
00006 #include "alias.h"
00007 #include "loadcase.h"
00008 #include "dloadcase.h"
00009 #include "inicd.h"
00010 #include "vector.h"
00011 #include "matrix.h"
00012 #include "elemhead.h"
00013 #include "gtopology.h"
00014 #include "element.h"
00015 #include "node.h"
00016 #include "elemswitch.h"
00017 #include "intpoints.h"
00018
00019
00020
00021
00022
00023
00024
00025
00026 mechbclc::mechbclc ()
00027 {
00028
00029 nlc=0L;
00030
00031 nico = 0L;
00032
00033 lc=NULL;
00034
00035 dlc = NULL;
00036 ico = NULL;
00037 ncsum = 0L;
00038 sumcomp = NULL;
00039 reactsumcomp = NULL;
00040 pd_reactsumcomp = NULL;
00041
00042
00043 ngfes=0L;
00044
00045 eigstrfun = NULL;
00046 strastre = (strastrestate)0;
00047 }
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 mechbclc::~mechbclc ()
00058 {
00059 if (lc!=NULL)
00060 delete [] lc;
00061 if (dlc!=NULL)
00062 delete [] dlc;
00063 if (ico!=NULL)
00064 delete [] ico;
00065 if (sumcomp!=NULL)
00066 delete [] sumcomp;
00067 if (reactsumcomp!=NULL)
00068 delete [] reactsumcomp;
00069 if (pd_reactsumcomp!=NULL)
00070 delete [] pd_reactsumcomp;
00071 if (eigstrfun!=NULL)
00072 delete [] eigstrfun;
00073 }
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087 void mechbclc::read (XFILE *in)
00088 {
00089 long i;
00090 vector nodval,ipval;
00091
00092 switch (Mp->tprob){
00093 case linear_statics:{
00094 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00095 if (nlc<0){
00096 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00097 abort ();
00098 }
00099 lc = new loadcase [nlc];
00100 for (i=0;i<nlc;i++){
00101 lc[i].read (in);
00102 }
00103 break;
00104 }
00105 case mat_nonlinear_statics:{
00106 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00107 if (nlc!=1){
00108 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00109 abort ();
00110 }
00111 lc = new loadcase [2*nlc];
00112 for (i=0;i<2*nlc;i++){
00113 lc[i].read (in);
00114 }
00115 readinic(in);
00116 break;
00117 }
00118 case geom_nonlinear_statics:{
00119 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00120 if (nlc!=1){
00121 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00122 abort ();
00123 }
00124 lc = new loadcase [2*nlc];
00125 for (i=0;i<2*nlc;i++){
00126 lc[i].read (in);
00127 }
00128 readinic(in);
00129 break;
00130 }
00131 case eigen_dynamics:{
00132 if (Mp->straincomp==1 || Mp->stresscomp==1)
00133 nlc=Mp->eigsol.neigv;
00134 break;
00135 }
00136 case forced_dynamics:{
00137 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00138 if (nlc<0){
00139 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00140 abort ();
00141 }
00142 lc = new loadcase [nlc];
00143 dlc = new dloadcase [nlc];
00144 for (i=0;i<nlc;i++){
00145 lc[i].read (in);
00146 dlc[i].read (in);
00147 }
00148 break;
00149 }
00150 case mech_timedependent_prob:{
00151 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00152 if (nlc!=1){
00153 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00154 abort ();
00155 }
00156 dlc = new dloadcase [nlc];
00157 for (i=0;i<nlc;i++){
00158 dlc[i].read (in);
00159 }
00160 readinic(in);
00161 break;
00162 }
00163 case growing_mech_structure:{
00164 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00165 if (nlc!=1){
00166 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00167 abort ();
00168 }
00169 dlc = new dloadcase [nlc];
00170 for (i=0;i<nlc;i++){
00171 dlc[i].read (in);
00172 }
00173 readinic(in);
00174 break;
00175 }
00176 case earth_pressure:{
00177 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00178 if (Mp->tepsol < epplast_sol)
00179 {
00180 lc = new loadcase [nlc];
00181 dlc = new dloadcase [nlc];
00182 for (i=0;i<nlc;i++){
00183 lc[i].read (in);
00184 dlc[i].read (in);
00185 }
00186 readinic(in);
00187 }
00188 else
00189 {
00190 lc = new loadcase [2*nlc];
00191 for (i=0;i<2*nlc;i++){
00192 lc[i].read (in);
00193 }
00194 readinic(in);
00195 }
00196 break;
00197 }
00198 case layered_linear_statics:{
00199 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00200 if (nlc<0){
00201 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00202 abort ();
00203 }
00204 lc = new loadcase [nlc];
00205 for (i=0;i<nlc;i++){
00206 lc[i].read (in);
00207 }
00208 break;
00209 }
00210
00211 case lin_floating_subdomain:{
00212 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00213 if (nlc<0){
00214 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00215 abort ();
00216 }
00217 lc = new loadcase [nlc];
00218 for (i=0;i<nlc;i++){
00219 lc[i].read (in);
00220 }
00221 break;
00222 }
00223
00224 case nonlin_floating_subdomain:{
00225 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00226 if (nlc!=1){
00227 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00228 abort ();
00229 }
00230 lc = new loadcase [2*nlc];
00231 for (i=0;i<2*nlc;i++){
00232 lc[i].read (in);
00233 }
00234 break;
00235 }
00236
00237 case nonlinear_dynamics:{
00238 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00239 if (nlc!=1){
00240 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00241 abort ();
00242 }
00243 lc = new loadcase [nlc];
00244 dlc = new dloadcase [nlc];
00245 for (i=0;i<nlc;i++){
00246 lc[i].read (in);
00247 dlc[i].read (in);
00248 }
00249 break;
00250 }
00251
00252 case hemivar_inequal:{
00253 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00254 if (nlc<0){
00255 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00256 abort ();
00257 }
00258 lc = new loadcase [nlc];
00259 for (i=0;i<nlc;i++){
00260 lc[i].read (in);
00261 }
00262 break;
00263 }
00264
00265 case load_balancing:{
00266 xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
00267 if (nlc<0){
00268 print_err("wrong number of load cases", __FILE__, __LINE__, __func__);
00269 abort ();
00270 }
00271 lc = new loadcase [nlc];
00272 for (i=0;i<nlc;i++){
00273 lc[i].read (in);
00274 }
00275 break;
00276 }
00277
00278 default:{
00279 print_err("unknown problem type", __FILE__, __LINE__, __func__);
00280 }
00281 }
00282
00283
00284
00285
00286 read_eigenstrains (in);
00287
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 void mechbclc::print (FILE *out)
00301 {
00302 long i;
00303
00304 fprintf (out,"\n number of load cases:");
00305 switch (Mp->tprob){
00306 case linear_statics:{
00307 fprintf (out,"\n %ld",nlc);
00308 for (i=0;i<nlc;i++){
00309 lc[i].print (out);
00310 }
00311 break;
00312 }
00313 case mat_nonlinear_statics:{
00314 fprintf (out,"\n %ld",nlc);
00315 for (i=0;i<2*nlc;i++){
00316 lc[i].print (out);
00317 }
00318 printinic(out);
00319 break;
00320 }
00321 case geom_nonlinear_statics:{
00322 fprintf (out,"\n %ld",nlc);
00323 for (i=0;i<2*nlc;i++){
00324 lc[i].print (out);
00325 }
00326 printinic(out);
00327 break;
00328 }
00329 case eigen_dynamics:{
00330 break;
00331 }
00332 case forced_dynamics:{
00333 fprintf (out,"\n %ld",nlc);
00334 for (i=0;i<nlc;i++){
00335 lc[i].print (out);
00336 dlc[i].print (out);
00337 }
00338 break;
00339 }
00340 case mech_timedependent_prob:{
00341 fprintf (out,"\n %ld",nlc);
00342 for (i=0;i<nlc;i++){
00343 dlc[i].print (out);
00344 }
00345 printinic(out);
00346 break;
00347 }
00348 case growing_mech_structure:{
00349 fprintf (out,"\n %ld",nlc);
00350 for (i=0;i<nlc;i++){
00351 dlc[i].print (out);
00352 }
00353 printinic(out);
00354 break;
00355 }
00356 case earth_pressure:{
00357 fprintf (out,"\n %ld",nlc);
00358 if (Mp->tepsol < epplast_sol)
00359 {
00360 for (i=0;i<nlc;i++){
00361 lc[i].print (out);
00362 dlc[i].print (out);
00363 }
00364 printinic(out);
00365 }
00366 else
00367 {
00368 for (i=0;i<2*nlc;i++){
00369 lc[i].print (out);
00370 }
00371 printinic(out);
00372 }
00373 break;
00374 }
00375 case layered_linear_statics:{
00376 fprintf (out,"\n %ld",nlc);
00377 for (i=0;i<nlc;i++){
00378 lc[i].print (out);
00379 }
00380 break;
00381 }
00382
00383 case lin_floating_subdomain:{
00384 fprintf (out,"\n %ld",nlc);
00385 for (i=0;i<nlc;i++){
00386 lc[i].print (out);
00387 }
00388 break;
00389 }
00390
00391 case nonlin_floating_subdomain:{
00392 fprintf (out,"\n %ld",nlc);
00393 for (i=0;i<2*nlc;i++){
00394 lc[i].print (out);
00395 }
00396 break;
00397 }
00398
00399 case nonlinear_dynamics:{
00400 fprintf (out,"\n %ld",nlc);
00401 for (i=0;i<nlc;i++){
00402 lc[i].print (out);
00403 dlc[i].print (out);
00404 }
00405 break;
00406 }
00407
00408 case hemivar_inequal:{
00409 fprintf (out,"\n %ld",nlc);
00410 for (i=0;i<nlc;i++){
00411 lc[i].print (out);
00412 }
00413 break;
00414 }
00415
00416 case load_balancing:{
00417 fprintf (out,"\n %ld",nlc);
00418 for (i=0;i<nlc;i++){
00419 lc[i].print (out);
00420 }
00421 break;
00422 }
00423
00424 default:{
00425 print_err("unknown problem type", __FILE__, __LINE__, __func__);
00426 }
00427 }
00428
00429
00430
00431
00432 print_eigenstrains (out);
00433
00434 }
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444 void mechbclc::read_eigenstrains (XFILE *in)
00445 {
00446 long i,j,k;
00447 long ncomp,nip,ipid;
00448 long *id;
00449
00450
00451
00452
00453
00454
00455
00456 xfscanf (in,"%k%ld","eigenstrains",&(Mp->eigstrains));
00457 if (Mp->eigstrains<0 || Mp->eigstrains>5){
00458 print_err("wrong definition of eigenstrains", __FILE__, __LINE__, __func__);
00459 abort ();
00460 }
00461
00462 if (Mespr == 1){
00463 if (Mp->eigstrains==0)
00464 fprintf(stdout, "\n eigenstrains are not defined");
00465 if (Mp->eigstrains==1)
00466 fprintf(stdout, "\n eigenstrains are defined for the whole problem by a single set of values");
00467 if (Mp->eigstrains==2)
00468 fprintf(stdout, "\n eigenstrains are defined independently for each element");
00469 if (Mp->eigstrains==3)
00470 fprintf(stdout, "\n eigenstrains are defined independently in TRFEL");
00471 if (Mp->eigstrains==4)
00472 fprintf(stdout, "\n eigenstresses are defined for the whole problem by a single set of values");
00473 if (Mp->eigstrains==5)
00474 fprintf(stdout, "\n eigenstresses are defined independently for each element");
00475 }
00476
00477
00478 if (Mp->eigstrains>0){
00479
00480
00481
00482 if (Mm->eigstrid==NULL)
00483 Mm->eigstrid = new long* [Mm->tnip];
00484
00485 if (Mm->eigstrains==NULL)
00486 Mm->eigstrains = new double* [Mm->tnip];
00487
00488 if (Mm->eigstresses==NULL)
00489 Mm->eigstresses = new double* [Mm->tnip];
00490
00491 id = new long [6];
00492 }
00493
00494
00495
00496
00497
00498 if (Mp->eigstrains==1){
00499
00500
00501
00502
00503 xfscanf (in,"%m", &strastrestate_kwdset, (int *)&strastre);
00504
00505 switch (strastre){
00506 case bar:{
00507 xfscanf (in,"%ld",id);
00508 ncomp=1;
00509 break;
00510 }
00511 case planestress:{
00512 xfscanf (in,"%ld %ld %ld",id,id+1,id+2);
00513 ncomp=3;
00514 break;
00515 }
00516 case spacestress:{
00517 xfscanf (in,"%ld %ld %ld %ld %ld %ld",id,id+1,id+2,id+3,id+4,id+5);
00518 ncomp=6;
00519 break;
00520 }
00521 default:{
00522 print_err("unknown type of strain-stress state", __FILE__, __LINE__, __func__);
00523 }
00524 }
00525
00526 for (i=0;i<Mm->tnip;i++){
00527
00528 Mm->eigstrains[i] = new double [ncomp];
00529 memset(Mm->eigstrains[i], 0, sizeof(*Mm->eigstrains[i])*ncomp);
00530
00531 Mm->eigstresses[i] = new double [ncomp];
00532 memset(Mm->eigstresses[i], 0, sizeof(*Mm->eigstresses[i])*ncomp);
00533
00534 Mm->eigstrid[i] = new long [ncomp];
00535
00536 for (j=0;j<ncomp;j++){
00537 Mm->eigstrid[i][j]=id[j]-1;
00538 }
00539 }
00540
00541 }
00542
00543
00544
00545
00546
00547 if (Mp->eigstrains==2){
00548
00549
00550 for (i=0;i<Mt->ne;i++){
00551
00552 strastre = Mt->give_ssst (i,0);
00553
00554 switch (strastre){
00555 case bar:{
00556 xfscanf (in,"%ld",id);
00557 ncomp=1;
00558 break;
00559 }
00560 case planestress:{
00561 xfscanf (in,"%ld %ld %ld",id,id+1,id+2);
00562 ncomp=3;
00563 break;
00564 }
00565 case spacestress:{
00566 xfscanf (in,"%ld %ld %ld %ld %ld %ld",id,id+1,id+2,id+3,id+4,id+5);
00567 ncomp=6;
00568 break;
00569 }
00570 default:{
00571 print_err("unknown type of strain-stress state", __FILE__, __LINE__, __func__);
00572 }
00573 }
00574
00575
00576 nip = Mt->give_tnip (i);
00577
00578 ipid=Mt->elements[i].ipp[0][0];
00579
00580
00581 for (j=0;j<nip;j++){
00582 Mm->eigstrains[ipid] = new double [ncomp];
00583 Mm->eigstresses[ipid] = new double [ncomp];
00584 Mm->eigstrid[ipid] = new long [ncomp];
00585
00586 for (k=0;k<ncomp;k++){
00587 Mm->eigstrid[ipid][k] = id[k]-1;
00588 }
00589
00590 ipid++;
00591 }
00592
00593 }
00594 }
00595
00596
00597 if (Mp->eigstrains==3){
00598
00599 for (j=0;j<Mm->tnip;j++){
00600 ncomp = Mm->ip[j].ncompstr;
00601 Mm->eigstrains[j] = new double[ncomp];
00602 memset(Mm->eigstrains[j], 0, sizeof(*Mm->eigstrains[j])*ncomp);
00603 Mm->eigstresses[j] = new double [ncomp];
00604 memset(Mm->eigstresses[j], 0, sizeof(*Mm->eigstresses[j])*ncomp);
00605 }
00606
00607 }
00608
00609
00610
00611
00612
00613 if (Mp->eigstrains==4){
00614
00615
00616
00617
00618 xfscanf (in,"%m", &strastrestate_kwdset, (int *)&strastre);
00619
00620 switch (strastre){
00621 case bar:{
00622 xfscanf (in,"%ld",id);
00623 ncomp=1;
00624 break;
00625 }
00626 case planestress:{
00627 xfscanf (in,"%ld %ld %ld",id,id+1,id+2);
00628 ncomp=3;
00629 break;
00630 }
00631 case axisymm:{
00632 xfscanf(in,"%ld %ld %ld %ld", id, id+1, id+2, id+3);
00633 ncomp = 4;
00634 break;
00635 }
00636 case spacestress:{
00637 xfscanf (in,"%ld %ld %ld %ld %ld %ld",id,id+1,id+2,id+3,id+4,id+5);
00638 ncomp=6;
00639 break;
00640 }
00641 default:{
00642 print_err("unknown type of strain-stress state", __FILE__, __LINE__, __func__);
00643 }
00644 }
00645
00646 for (i=0;i<Mm->tnip;i++){
00647
00648 Mm->eigstrains[i] = new double [ncomp];
00649 memset(Mm->eigstrains[i], 0, sizeof(*Mm->eigstrains[i])*ncomp);
00650
00651 Mm->eigstresses[i] = new double [ncomp];
00652 memset(Mm->eigstresses[i], 0, sizeof(*Mm->eigstresses[i])*ncomp);
00653
00654 Mm->eigstrid[i] = new long [ncomp];
00655
00656 for (j=0;j<ncomp;j++){
00657 Mm->eigstrid[i][j]=id[j]-1;
00658 }
00659 }
00660
00661 }
00662
00663
00664
00665
00666
00667 if (Mp->eigstrains==5){
00668
00669
00670 for (i=0;i<Mt->ne;i++){
00671
00672 strastre = Mt->give_ssst (i,0);
00673
00674 switch (strastre){
00675 case bar:{
00676 xfscanf (in,"%ld",id);
00677 ncomp=1;
00678 break;
00679 }
00680 case planestress:{
00681 xfscanf (in,"%ld %ld %ld",id,id+1,id+2);
00682 ncomp=3;
00683 break;
00684 }
00685 case axisymm:{
00686 xfscanf(in,"%ld %ld %ld %ld", id, id+1, id+2, id+3);
00687 ncomp = 4;
00688 break;
00689 }
00690 case spacestress:{
00691 xfscanf (in,"%ld %ld %ld %ld %ld %ld",id,id+1,id+2,id+3,id+4,id+5);
00692 ncomp=6;
00693 break;
00694 }
00695 default:{
00696 print_err("unknown type of strain-stress state", __FILE__, __LINE__, __func__);
00697 }
00698 }
00699
00700
00701 nip = Mt->give_tnip (i);
00702
00703 ipid=Mt->elements[i].ipp[0][0];
00704
00705
00706 for (j=0;j<nip;j++){
00707
00708 Mm->eigstrains[ipid] = new double [ncomp];
00709 memset(Mm->eigstrains[ipid], 0, sizeof(*Mm->eigstrains[ipid])*ncomp);
00710 Mm->eigstresses[ipid] = new double [ncomp];
00711 Mm->eigstrid[ipid] = new long [ncomp];
00712
00713 for (k=0;k<ncomp;k++){
00714 Mm->eigstrid[ipid][k] = id[k]-1;
00715 }
00716
00717 ipid++;
00718 }
00719
00720 }
00721 }
00722
00723
00724
00725 if (Mp->eigstrains>0){
00726 delete [] id;
00727 }
00728
00729
00730
00731
00732 if (Mp->eigstrains==1 || Mp->eigstrains==2 || Mp->eigstrains==4 || Mp->eigstrains==5){
00733
00734 xfscanf (in,"%k%ld","number_of_gfunct_eigenstr",&ngfes);
00735
00736 eigstrfun = new gfunct [ngfes];
00737
00738
00739 for (i=0;i<ngfes;i++){
00740 eigstrfun[i].read (in);
00741 }
00742
00743
00744 eigstrain_computation (Mp->time);
00745
00746 }
00747
00748 }
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759 void mechbclc::print_eigenstrains (FILE *out)
00760 {
00761 long i;
00762 long ipid;
00763
00764
00765
00766
00767
00768 fprintf (out,"\n\n%ld",Mp->eigstrains);
00769
00770
00771
00772
00773 if ((Mp->eigstrains==1) || (Mp->eigstrains==4)){
00774
00775
00776
00777
00778 fprintf (out,"\n\n%d",(int)strastre);
00779
00780 switch (strastre){
00781 case bar:{
00782 fprintf (out,"\n %ld",Mm->eigstrid[0][0]+1);
00783 break;
00784 }
00785 case planestress:{
00786 fprintf (out,"\n %ld %ld %ld",Mm->eigstrid[0][0]+1,Mm->eigstrid[0][1]+1,Mm->eigstrid[0][2]+1);
00787 break;
00788 }
00789 case axisymm:{
00790 fprintf (out,"\n %ld %ld %ld %ld",Mm->eigstrid[0][0]+1,Mm->eigstrid[0][1]+1,Mm->eigstrid[0][2]+1,Mm->eigstrid[0][3]+1);
00791 break;
00792 }
00793 case spacestress:{
00794 fprintf (out,"\n %ld %ld %ld %ld %ld %ld",Mm->eigstrid[0][0]+1,Mm->eigstrid[0][1]+1,Mm->eigstrid[0][2]+1
00795 ,Mm->eigstrid[0][3]+1,Mm->eigstrid[0][4]+1,Mm->eigstrid[0][5]+1);
00796 break;
00797 }
00798 default:{
00799 print_err("unknown type of strain-stress state", __FILE__, __LINE__, __func__);
00800 }
00801 }
00802 }
00803
00804
00805
00806
00807
00808 if ((Mp->eigstrains==2) || (Mp->eigstrains==5)){
00809
00810
00811 for (i=0;i<Mt->ne;i++){
00812
00813
00814 strastre = Mt->give_ssst (i,0);
00815
00816
00817 ipid=Mt->elements[i].ipp[0][0];
00818
00819 switch (strastre){
00820 case bar:{
00821 fprintf (out,"\n %ld",Mm->eigstrid[ipid][0]+1);
00822 break;
00823 }
00824 case planestress:{
00825 fprintf (out,"\n %ld %ld %ld",Mm->eigstrid[ipid][0]+1,Mm->eigstrid[ipid][1]+1,Mm->eigstrid[ipid][2]+1);
00826 break;
00827 }
00828 case axisymm:{
00829 fprintf (out,"\n %ld %ld %ld %ld",Mm->eigstrid[ipid][0]+1,Mm->eigstrid[ipid][1]+1,Mm->eigstrid[ipid][2]+1,Mm->eigstrid[ipid][3]+1);
00830 break;
00831 }
00832 case spacestress:{
00833 fprintf (out,"\n %ld %ld %ld %ld %ld %ld",Mm->eigstrid[ipid][0]+1,Mm->eigstrid[ipid][1]+1,Mm->eigstrid[ipid][2]+1
00834 ,Mm->eigstrid[ipid][3]+1,Mm->eigstrid[ipid][4]+1,Mm->eigstrid[ipid][5]+1);
00835 break;
00836 }
00837 default:{
00838 print_err("unknown type of strain-stress state", __FILE__, __LINE__, __func__);
00839 }
00840 }
00841 }
00842 }
00843
00844
00845
00846
00847
00848 if (Mp->eigstrains==1 || Mp->eigstrains==2 || Mp->eigstrains==4 || Mp->eigstrains==5){
00849
00850 fprintf (out,"\n\n %ld",ngfes);
00851
00852
00853 for (i=0;i<ngfes;i++){
00854 eigstrfun[i].print (out);
00855 }
00856
00857 }
00858
00859 }
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869 void mechbclc::eigstrain_computation (double time)
00870 {
00871 long i,j,gfid,ncomp;
00872
00873
00874 for (i=0;i<Mm->tnip;i++){
00875
00876 ncomp = Mm->ip[i].give_ncompstr ();
00877
00878
00879 for (j=0;j<ncomp;j++){
00880
00881 gfid = Mm->eigstrid[i][j];
00882
00883 if (Mp->eigstrains < 3)
00884 Mm->eigstrains[i][j] = eigstrfun[gfid].getval (time);
00885 if (Mp->eigstrains > 3)
00886 Mm->eigstresses[i][j] = -eigstrfun[gfid].getval (time);
00887
00888 }
00889
00890 }
00891
00892 }
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906 long mechbclc::readinic (XFILE *in)
00907 {
00908 long i, n;
00909
00910 xfscanf(in, "%k%ld", "num_ini_cond", &nico);
00911
00912 if ((nico < 0) && (nico > Mt->nn))
00913 {
00914 print_err("invalid number of inital conditions", __FILE__, __LINE__, __func__);
00915 return 1;
00916 }
00917 if (nico == 0)
00918 return 0;
00919 ico = new inicd[Mt->nn];
00920 for (i = 0L; i < nico; i++)
00921 {
00922 xfscanf(in, "%ld", &n);
00923 if ((n < 0) || (n > Mt->nn))
00924 {
00925 print_err("invalid node number", __FILE__, __LINE__, __func__);
00926 return 2;
00927 }
00928 ico[n-1].read(in);
00929 if ((ico[n-1].type & inicond) && (Mm->ic == NULL))
00930 {
00931 Mm->ic = new double *[Mm->tnip];
00932 memset(Mm->ic, 0, sizeof(*Mm->ic)*Mm->tnip);
00933 }
00934 }
00935 return 0;
00936 }
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952 long mechbclc::printinic (FILE *out)
00953 {
00954 long i, n=1;
00955
00956 fprintf(out, "\n ## initial conditions:");
00957 fprintf(out, "\n %ld\n", nico);
00958
00959 for (i = 0L; i < nico; i++)
00960 {
00961 fprintf(out, "\n %ld ",n);
00962 ico[n-1].print(out);
00963 }
00964 return 0;
00965 }
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977 void mechbclc::inicipval(void)
00978 {
00979 long i, j, k, l, nne, nval, eid;
00980 elemtype et;
00981 ivector enod;
00982 matrix nodval;
00983 inictype *ictn;
00984
00985
00986 for (i = 0; i < Mt->ne; i++)
00987 {
00988 if (Gtm->leso[i]==1)
00989 {
00990 nne = Mt->give_nne(i);
00991 allocv(nne, enod);
00992
00993
00994 Mt->give_elemnodes (i, enod);
00995 nval = ico[enod[0]].nval;
00996 if (ico[enod[0]].type & inidisp)
00997 nval -= Gtm->give_ndofn(enod[0]);
00998 for (j = 0; j < nne; j++)
00999 {
01000 if (nval < ico[enod[j]].nval)
01001 {
01002 nval = ico[enod[j]].nval;
01003 if (ico[enod[j]].type & inidisp)
01004 nval -= Gtm->give_ndofn(enod[j]);
01005 }
01006 }
01007 allocm (nne, nval, nodval);
01008 fillm(0.0, nodval);
01009 ictn = new inictype[nne];
01010 for (j = 0; j < nne; j++)
01011 {
01012 ictn[j] = ico[enod[j]].type;
01013 for (k = 0; k < nval; k++)
01014 {
01015 if (ictn[j] & inidisp)
01016 l = k + Gtm->give_ndofn(enod[j]);
01017 else
01018 l = k;
01019 nodval[j][k] = ico[enod[j]].val[l];
01020 }
01021 }
01022 eid = i;
01023 et = Mt->give_elem_type(eid);
01024 switch (et)
01025 {
01026 case bar2d:{ Bar2d->inicipval(eid, 0, 0, nodval, ictn); break; }
01027 case barq2d:{ Barq2d->inicipval(eid, 0, 0, nodval, ictn); break; }
01028
01029
01030
01031 case spring_1:{ Spring->inicipval(eid, 0, 0, nodval, ictn); break; }
01032 case spring_2:{ Spring->inicipval(eid, 0, 0, nodval, ictn); break; }
01033 case spring_3:{ Spring->inicipval(eid, 0, 0, nodval, ictn); break; }
01034 case spring_4:{ Spring->inicipval(eid, 0, 0, nodval, ictn); break; }
01035 case spring_5:{ Spring->inicipval(eid, 0, 0, nodval, ictn); break; }
01036 case spring_6:{ Spring->inicipval(eid, 0, 0, nodval, ictn); break; }
01037 case planeelementlt:{ Pelt->inicipval(eid, 0, 0, nodval, ictn); break; }
01038 case planeelementqt:{ Peqt->inicipval(eid, 0, 0, nodval, ictn); break; }
01039 case planeelementrotlt:{ Perlt->inicipval(eid, 0, 0, nodval, ictn); break; }
01040 case planeelementlq:{ Pelq->inicipval(eid, 0, 0, nodval, ictn); break; }
01041 case planeelementqq:{ Peqq->inicipval(eid, 0, 0, nodval, ictn); break; }
01042 case planeelementrotlq:{ Perlq->inicipval(eid, 0, 0, nodval, ictn); break; }
01043 case planeelementsubqt:{ Pesqt->inicipval(eid, 0, 0, nodval, ictn); break; }
01044 case cctel:{ Cct->inicipval(eid, 0, 0, nodval, ictn); break; }
01045 case dktel:{ Dkt->inicipval(eid, 0, 0, nodval, ictn); break; }
01046 case dstel:{ Dst->inicipval(eid, 0, 0, nodval, ictn); break; }
01047 case q4plateel:{ Q4pl->inicipval(eid, 0, 0, nodval, ictn); break; }
01048
01049
01050 case axisymmlq:{ Asymlq->inicipval(eid, 0, 0, nodval, ictn); break; }
01051 case axisymmlt:{ Asymlt->inicipval(eid, 0, 0, nodval, ictn); break; }
01052 case shelltrelem:{ Shtr->inicipval(eid, 0, 0, nodval, ictn); break; }
01053 case shellqelem:{ Shq->inicipval(eid, 0, 0, nodval, ictn); break; }
01054 case lineartet:{ Ltet->inicipval(eid, 0, 0, nodval, ictn); break; }
01055 case quadrtet:{ Qtet->inicipval(eid, 0, 0, nodval, ictn); break; }
01056 case linearhex:{ Lhex->inicipval(eid, 0, 0, nodval, ictn); break; }
01057 case quadrhex:{ Qhex->inicipval(eid, 0, 0, nodval, ictn); break; }
01058 default:
01059 {
01060 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01061 }
01062 }
01063 destrv(enod);
01064 destrm(nodval);
01065 delete [] ictn;
01066 }
01067 }
01068 }
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080 void mechbclc::alloc_sumcomp ()
01081 {
01082 ncsum = Mt->give_ndofn (0);
01083 sumcomp = new double [ncsum];
01084 reactsumcomp = new double [ncsum];
01085 pd_reactsumcomp = new double [ncsum];
01086 memset(sumcomp, 0, sizeof(*sumcomp)*ncsum);
01087 memset(reactsumcomp, 0, sizeof(*reactsumcomp)*ncsum);
01088 memset(pd_reactsumcomp, 0, sizeof(*pd_reactsumcomp)*ncsum);
01089 }
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101 void mechbclc::comp_sum (double *rhs)
01102 {
01103 long i,j,ndofn,cn;
01104 long *aux;
01105
01106 aux=new long [Ndofm];
01107 for (i=0;i<Ndofm;i++){
01108 aux[i]=0;
01109 }
01110
01111 nullv (sumcomp,ncsum);
01112
01113 for (i=0;i<Mt->nn;i++){
01114 ndofn=Gtm->gnodes[i].ndofn;
01115 for (j=0;j<ndofn;j++){
01116 cn=Mt->give_dof (i,j);
01117 if (cn>0){
01118 if (aux[cn-1]==0){
01119 sumcomp[j]+=rhs[cn-1];
01120 aux[cn-1]=1;
01121 }
01122 }
01123 }
01124 }
01125
01126 delete [] aux;
01127 }
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139 void mechbclc::comp_sum_react()
01140 {
01141 long i,j,ndofn,cn;
01142
01143 nullv (reactsumcomp,ncsum);
01144
01145 for (i=0;i<Mt->nn;i++){
01146 ndofn=Gtm->gnodes[i].ndofn;
01147 for (j=0;j<ndofn;j++){
01148 cn=Mt->give_dof(i,j);
01149 if (cn==0)
01150 reactsumcomp[j]+=Mt->nodes[i].r[j];
01151 }
01152 }
01153 }
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165 void mechbclc::comp_sum_pdreact()
01166 {
01167 long i,j,ndofn,cn;
01168
01169 nullv (sumcomp,ncsum);
01170
01171 for (i=0;i<Mt->nn;i++){
01172 ndofn=Gtm->gnodes[i].ndofn;
01173 for (j=0;j<ndofn;j++){
01174 cn=Mt->give_dof (i,j);
01175 if (cn<0)
01176 pd_reactsumcomp[j]+=Mt->nodes[i].r[j];
01177 }
01178 }
01179 }
01180
01181
01182
01183
01184
01185
01186
01187
01188 void mechbclc::give_comp_sum (double *sum)
01189 {
01190 long i;
01191 for (i=0;i<ncsum;i++)
01192 sum[i]=sumcomp[i];
01193 }
01194