00001 #include <errno.h>
00002 #include "globmat.h"
00003 #include "genfile.h"
00004 #include "mathem.h"
00005 #include "global.h"
00006 #include "elemhead.h"
00007 #include "loadcase.h"
00008 #include "dloadcase.h"
00009 #include "dloadpd.h"
00010 #include "node.h"
00011 #include "element.h"
00012 #include "intpoints.h"
00013 #include "problem.h"
00014 #include "mechmat.h"
00015 #include "creepb.h"
00016 #include "creepbs.h"
00017 #include "mechcrsec.h"
00018 #include "elemswitch.h"
00019 #include "vecttens.h"
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 void stiffness_matrix (long lcid)
00033 {
00034 long i,j,nl,ndofe,ndofemn,ndofn,nmult,nid;
00035 ivector cn,cna;
00036 matrix lm,lmt;
00037 matrix alm,bd,db,dd;
00038 time_t bt,et;
00039 int err;
00040
00041 bt = time (NULL);
00042
00043 if (Smat==NULL)
00044 Smat = new gmatrix;
00045
00046 Smat->setval (Mp->ssle);
00047 Smat->initiate (Gtm,Ndofm,Mp->tstorsm,Mespr,Out);
00048
00049
00050
00051
00052
00053 if (Mp->lbtype==1)
00054 abort ();
00055
00056
00057
00058
00059
00060
00061
00062 for (i=0;i<Mt->ne;i++){
00063 if (Gtm->leso[i]==1){
00064
00065
00066
00067
00068
00069 ndofe = Mt->give_ndofe (i);
00070
00071 reallocm (ndofe,ndofe,lm);
00072
00073 stiffmat (lcid,i,lm);
00074
00075 err =check_math_errel(i);
00076
00077
00078
00079 ndofemn = Gtm->give_ndofe (i);
00080 reallocv (ndofemn,cn);
00081 Mt->give_code_numbers (i,cn.a);
00082
00083 if (ndofe != ndofemn){
00084 if (Mp->homog!=3){
00085
00086
00087
00088
00089 reallocm (ndofe,ndofemn,alm);
00090 mxm (lm,*Mt->elements[i].tmat,alm);
00091 reallocm (ndofemn,ndofemn,lm);
00092 mtxm (*Mt->elements[i].tmat,alm,lm);
00093 }else{
00094
00095
00096 long ii,jj,kk;
00097 reallocm (ndofe,ndofe,alm);
00098 copym (lm,alm);
00099 reallocm (ndofemn,ndofemn,lm);
00100 for (ii=0;ii<ndofe;ii++){
00101 for (jj=0;jj<ndofe;jj++){
00102 lm(ii, jj)=alm(ii, jj);
00103 }
00104 }
00105 reallocm (ndofe,ndofemn-ndofe,bd);
00106 Ltet->bd_matrix (i,0,0,bd);
00107 for (ii=0;ii<ndofe;ii++){
00108 kk=ndofe;
00109 for (jj=0;jj<ndofemn-ndofe;jj++){
00110 lm(ii, kk)=bd(ii, jj);
00111 kk++;
00112 }
00113 }
00114 kk=ndofe;
00115 for (ii=0;ii<ndofemn-ndofe;ii++){
00116 for (jj=0;jj<ndofe;jj++){
00117 lm(kk, jj)=bd(jj, ii);
00118 }
00119 kk++;
00120 }
00121
00122 reallocm (ndofemn-ndofe,ndofemn-ndofe,dd);
00123 Ltet->dd_matrix (i,0,0,dd);
00124 for (ii=0;ii<ndofemn-ndofe;ii++){
00125 for (jj=0;jj<ndofemn-ndofe;jj++){
00126 lm(ii+ndofe, jj+ndofe)=dd(ii, jj);
00127 }
00128 }
00129
00130 }
00131 }
00132
00133 Smat->localize (lm,cn,i);
00134 }
00135 }
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165 if (Mp->tprob==layered_linear_statics){
00166 for (i=0;i<Mt->nln;i++){
00167 nl=Gtm->lgnodes[i].nl;
00168 ndofn=Gtm->lgnodes[i].ndofn;
00169 nmult=Gtm->lgnodes[i].nmult;
00170
00171 reallocm (ndofn*2,nmult,lm);
00172 reallocm (nmult,ndofn*2,lmt);
00173
00174 reallocv (ndofn*2,cn);
00175
00176 reallocv (nmult,cna);
00177
00178 for (j=1;j<nl;j++){
00179 constr_matrix (i,j,lm);
00180
00181 nid=Gtm->lgnodes[i].nodes[j-1];
00182 Mt->give_node_code_numbers (nid,cn.a);
00183 nid=Gtm->lgnodes[i].nodes[j];
00184 Mt->give_node_code_numbers (nid,cn.a+ndofn);
00185 Mt->give_mult_code_numbers (i,j,cna.a);
00186
00187 tranm (lm,lmt);
00188
00189 Smat->glocalize (lmt,cna,cn);
00190 Smat->glocalize (lm,cn,cna);
00191
00192
00193 }
00194 }
00195 }
00196
00197 Smat->prepmat (0.0,Mespr);
00198
00199
00200
00201
00202 et = time (NULL);
00203 fprintf (stdout,"\n stiffness matrix assembling time %ld",et-bt);
00204 fflush(stdout);
00205
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220 void mass_matrix (long lcid)
00221 {
00222 long i,ndofe,ndofemn;
00223 ivector cn;
00224 matrix lm;
00225 matrix alm;
00226
00227 if (Mmat==NULL) Mmat = new gmatrix;
00228
00229 Mmat->setval (Mp->ssle);
00230 Mmat->initiate (Gtm,Ndofm,Mp->tstormm,Mespr,Out);
00231
00232
00233
00234 for (i=0;i<Mt->ne;i++){
00235 if (Gtm->leso[i]==1){
00236
00237 ndofe = Mt->give_ndofe (i);
00238
00239 reallocm (ndofe,ndofe,lm);
00240
00241 massmat (lcid,i,lm);
00242
00243
00244
00245
00246
00247 ndofemn = Gtm->give_ndofe (i);
00248 reallocv (ndofemn,cn);
00249 Mt->give_code_numbers (i,cn.a);
00250
00251 if (ndofe != ndofemn){
00252
00253
00254 reallocm (ndofe,ndofemn,alm);
00255 mxm (lm,*Mt->elements[i].tmat,alm);
00256 reallocm (ndofemn,ndofemn,lm);
00257 mtxm (*Mt->elements[i].tmat,alm,lm);
00258 }
00259
00260
00261 Mmat->localize (lm,cn,i);
00262 }
00263 }
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 Mmat->prepmat (0.0,Mespr);
00286 }
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 void initial_stiffness_matrix (long lcid)
00304 {
00305 if (Ismat==NULL) Ismat = new gmatrix;
00306
00307 Ismat->setval (Mp->ssle);
00308 Ismat->initiate (Gtm,Ndofm,Mp->tstorsm,Mespr,Out);
00309
00310
00311
00312 long i,ndofe,ndofemn;
00313 ivector cn;
00314 matrix lm;
00315 matrix alm;
00316
00317 for (i=0;i<Mt->ne;i++){
00318 if (Gtm->leso[i]==1){
00319
00320 ndofe = Mt->give_ndofe (i);
00321
00322 reallocm (ndofe,ndofe,lm);
00323
00324
00325 initstiffmat (lcid,i,lm);
00326
00327
00328
00329 ndofemn = Gtm->give_ndofe (i);
00330 reallocv (ndofemn,cn);
00331
00332 Mt->give_code_numbers (i,cn.a);
00333
00334 if (ndofe != ndofemn){
00335
00336
00337 reallocm (ndofe,ndofemn,alm);
00338 mxm (lm,*Mt->elements[i].tmat,alm);
00339 reallocm (ndofemn,ndofemn,lm);
00340 mtxm (*Mt->elements[i].tmat,alm,lm);
00341 }
00342
00343
00344 Ismat->localize (lm,cn,i);
00345 }
00346 }
00347
00348 Ismat->prepmat (0.0,Mespr);
00349 }
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360 void damping_matrix ()
00361 {
00362 if (Dmat==NULL) Dmat = new gmatrix;
00363
00364 Dmat->setval (Mp->ssle);
00365 Dmat->initiate (Gtm,Ndofm,Mp->tstorsm,Mespr,Out);
00366
00367
00368 switch (Mp->damp){
00369 case nodamp:{
00370 break;
00371 }
00372 case massdamp:{
00373 Dmat->addgm (Mp->dmass,*Mmat);
00374 break;
00375 }
00376 case stiffdamp:{
00377 Dmat->addgm (Mp->dstiff,*Smat);
00378 break;
00379 }
00380 case rayleighdamp:{
00381 Dmat->addgm (Mp->dmass,*Mmat);
00382 Dmat->addgm (Mp->dstiff,*Smat);
00383 break;
00384 }
00385 default:{
00386 print_err("unknown type of damping is required",__FILE__,__LINE__,__func__);
00387 }
00388 }
00389
00390 Dmat->prepmat (0.0,Mespr);
00391 }
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 void internal_forces (long lcid,double *intfor)
00407 {
00408 if (Mp->matmodel==local){
00409 loc_internal_forces (lcid,intfor);
00410 }
00411 if (Mp->matmodel==nonlocal){
00412 nonloc_internal_forces (lcid,intfor);
00413 }
00414
00415
00416
00417
00418 if (Mp->tprob == lin_floating_subdomain || Mp->tprob == nonlin_floating_subdomain){
00419
00420 }
00421
00422
00423
00424
00425 Mp->strainstate=1;
00426
00427
00428 Mp->stressstate=1;
00429
00430
00431 Mp->otherstate=1;
00432 }
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448 void loc_internal_forces (long lcid,double *intfor)
00449 {
00450 long i,j,ne,ndofe,ndofemn, ncomp=0L;
00451 ivector cn;
00452 int err;
00453 vector ifor, aifor;
00454 vector r, mstr;
00455 matrix bd, dd;
00456
00457
00458 if (Mp->strcomp)
00459 compute_ipstrains (lcid);
00460
00461
00462 ne=Mt->ne;
00463
00464
00465 nullv (intfor,Ndofm);
00466
00467 for (i=0;i<ne;i++){
00468 if (Gtm->leso[i]==1){
00469
00470
00471
00472 ndofe=Mt->give_ndofe (i);
00473 reallocv (ndofe,ifor);
00474
00475 elem_internal_forces (i,lcid,ifor);
00476
00477 err = check_math_errel(i);
00478
00479 if (err)
00480 {
00481 print_err("standard library error No.%d detected on element %ld\n", __FILE__, __LINE__, __func__, err, i+1);
00482 abort();
00483 }
00484
00485
00486
00487 ndofemn = Gtm->give_ndofe (i);
00488 reallocv (ndofemn,cn);
00489
00490 Mt->give_code_numbers (i,cn.a);
00491
00492 if (ndofe != ndofemn)
00493 {
00494 if (Mp->homog!=3)
00495 {
00496
00497
00498
00499
00500 reallocv (ndofemn,aifor);
00501 mtxv (*Mt->elements[i].tmat,ifor,aifor);
00502 locglob (intfor,aifor.a,cn.a,ndofemn);
00503 }
00504 else
00505 {
00506
00507 if (ncomp == 0L)
00508 ncomp = ndofemn-ndofe;
00509
00510
00511
00512
00513 locglob (intfor,ifor.a,cn.a,ndofe);
00514
00515
00516
00517 reallocv (ncomp, ifor);
00518 elem_integration_quant(i, locstress, lcid, ifor);
00519 for (j=0; j<ncomp; j++)
00520 intfor[Ndofm-ncomp+j] += ifor[j];
00521 }
00522 }
00523 else
00524 {
00525
00526 locglob (intfor,ifor.a,cn.a,ndofe);
00527 }
00528 }
00529 }
00530 }
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546 void nonloc_internal_forces (long lcid,double *intfor)
00547 {
00548 long i,ne,ndofe,ndofemn;
00549 ivector cn;
00550 vector ifor;
00551 vector aifor;
00552
00553
00554 if (Mp->strcomp)
00555 compute_ipstrains (lcid);
00556
00557 ne=Mt->ne;
00558 nullv (intfor,Ndofm);
00559 Mp->nonlocphase=1;
00560
00561 for (i = 0; i < ne; i++)
00562 {
00563 if (Gtm->leso[i]==1)
00564 elem_local_values (i, lcid);
00565 }
00566
00567 for (i = 0; i < Mm->tnip; i++)
00568 Mm->nonlocaverage(i);
00569 Mp->nonlocphase=2;
00570
00571 for (i = 0; i < ne; i++)
00572 {
00573 if (Gtm->leso[i]==1)
00574 {
00575
00576
00577 ndofe=Mt->give_ndofe (i);
00578 reallocv (ndofe,ifor);
00579
00580 elem_nonloc_internal_forces(i, lcid, ifor);
00581
00582 check_math_errel(i);
00583
00584
00585
00586 ndofemn = Gtm->give_ndofe (i);
00587 reallocv (ndofemn,cn);
00588
00589 Mt->give_code_numbers (i,cn.a);
00590 if (ndofe != ndofemn){
00591
00592
00593 reallocv (ndofemn,aifor);
00594 mtxv (*Mt->elements[i].tmat,ifor,aifor);
00595 locglob (intfor,aifor.a,cn.a,ndofemn);
00596 }else{
00597
00598 locglob (intfor,ifor.a,cn.a,ndofe);
00599 }
00600 }
00601 }
00602 }
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617 void incr_internal_forces (long lcid,double *intfor)
00618 {
00619 long i,ne,ndofe,ndofemn;
00620 ivector cn;
00621 vector ifor;
00622 vector aifor;
00623
00624
00625 compute_ipstrains (lcid);
00626
00627
00628 ne=Mt->ne;
00629
00630
00631 nullv (intfor,Ndofm);
00632
00633 for (i=0;i<ne;i++){
00634 if (Gtm->leso[i]==1){
00635
00636
00637
00638 ndofe=Mt->give_ndofe (i);
00639 reallocv (ndofe,ifor);
00640
00641 elem_incr_internal_forces (i,lcid,ifor);
00642
00643 check_math_errel(i);
00644
00645
00646
00647 ndofemn = Gtm->give_ndofe (i);
00648 reallocv (ndofemn,cn);
00649
00650 Mt->give_code_numbers (i,cn.a);
00651
00652 if (ndofe != ndofemn){
00653
00654
00655 reallocv (ndofemn,aifor);
00656 mtxv (*Mt->elements[i].tmat,ifor,aifor);
00657 locglob (intfor,aifor.a,cn.a,ndofemn);
00658 }else{
00659
00660 locglob (intfor,ifor.a,cn.a,ndofe);
00661 }
00662 }
00663 }
00664 }
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
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
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750 void nodal_eigstrain_forces (long lcid,double *nodfor,double time)
00751 {
00752 long i,ne,ndofe,ndofemn;
00753 ivector cn;
00754 vector nfor;
00755 vector anfor;
00756
00757
00758
00759
00760
00761 if (Mp->eigstrains==1 || Mp->eigstrains==2)
00762 Mb->eigstrain_computation (time);
00763
00764 if (Mp->eigstrains==4 || Mp->eigstrains==5)
00765 {
00766 Mb->eigstrain_computation (time);
00767 Mp->eigstrcomp = 0;
00768 }
00769
00770
00771 ne=Mt->ne;
00772
00773 for (i=0;i<ne;i++){
00774 if (Gtm->leso[i]==1){
00775
00776
00777
00778 ndofe=Mt->give_ndofe (i);
00779 reallocv (ndofe,nfor);
00780
00781
00782 elem_eigstrain_forces (lcid,i,nfor);
00783
00784 check_math_errel(i);
00785
00786
00787
00788 ndofemn = Gtm->give_ndofe (i);
00789 reallocv (ndofemn,cn);
00790
00791 Mt->give_code_numbers (i,cn.a);
00792
00793 if (ndofe != ndofemn){
00794
00795
00796 reallocv (ndofemn,anfor);
00797 mtxv (*Mt->elements[i].tmat,nfor,anfor);
00798 locglob (nodfor,anfor.a,cn.a,ndofemn);
00799 }else{
00800
00801 locglob (nodfor,nfor.a,cn.a,ndofe);
00802 }
00803 }
00804 }
00805 }
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822 void nodal_pore_press_forces (long lcid,double *nodfor,double time)
00823 {
00824 long i,ne,ndofe,ndofemn, besc;
00825 ivector cn;
00826 vector nfor;
00827 vector anfor;
00828 double pp;
00829 matrix sigt(3,3);
00830 vector sig;
00831
00832
00833
00834
00835
00836 for(i=0; i<Mm->tnip; i++)
00837 {
00838 pp = Mm->givenonmechq(pore_pressure, i);
00839 reallocv(Mm->ip[i].ncompstr, sig);
00840 sigt[0][0] = sigt[1][1] = sigt[2][2] = -pp;
00841 tensor_vector(sig, sigt, Mm->ip[i].ssst, stress);
00842 Mm->storeeigstress(i, sig);
00843 }
00844
00845
00846 besc = Mp->eigstrcomp;
00847 Mp->eigstrcomp = 0;
00848
00849
00850 ne=Mt->ne;
00851
00852 for (i=0;i<ne;i++){
00853 if (Gtm->leso[i]==1){
00854
00855
00856
00857 ndofe=Mt->give_ndofe (i);
00858 reallocv (ndofe,nfor);
00859
00860
00861 elem_eigstrain_forces (lcid,i,nfor);
00862
00863 check_math_errel(i);
00864
00865
00866
00867 ndofemn = Gtm->give_ndofe (i);
00868 reallocv (ndofemn,cn);
00869
00870 Mt->give_code_numbers (i,cn.a);
00871
00872 if (ndofe != ndofemn){
00873
00874
00875 reallocv (ndofemn,anfor);
00876 mtxv (*Mt->elements[i].tmat,nfor,anfor);
00877 locglob (nodfor,anfor.a,cn.a,ndofemn);
00878 }else{
00879
00880 locglob (nodfor,nfor.a,cn.a,ndofe);
00881 }
00882 }
00883 }
00884
00885
00886 Mp->eigstrcomp = besc;
00887 }
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902 void eldispl (long lcid,long eid,double *r)
00903 {
00904 long i,j,ii,ndofe,ndofemn,nne,ndofn,acn;
00905 ivector cn,nod,ncn;
00906 vector rr;
00907
00908
00909
00910 ndofe = Mt->give_ndofe (eid);
00911
00912
00913 ndofemn = Gtm->give_ndofe (eid);
00914
00915 if (ndofe != ndofemn){
00916
00917 reallocv (ndofemn,cn);
00918 reallocv (ndofemn,rr);
00919 ndofe=ndofemn;
00920 }
00921 else{
00922
00923 reallocv (ndofe,cn);
00924 reallocv (ndofe,rr);
00925 }
00926
00927 Mt->give_code_numbers (eid,cn.a);
00928
00929
00930 switch (Mp->tprob){
00931 case forced_dynamics:
00932 case mech_timedependent_prob:{
00933
00934 for (i=0;i<ndofe;i++){
00935 ii=cn[i];
00936 if (ii<0) rr[i]=Mb->dlc[lcid].get_pd(Mp->time, ii);
00937 if (ii==0) rr[i]=0.0;
00938 if (ii>0) rr[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
00939 }
00940 break;
00941 }
00942
00943 case mat_nonlinear_statics:{
00944 for (i=0;i<ndofe;i++){
00945 ii=cn[i];
00946 if (ii<0)
00947 {
00948 rr[i] = 0.0;
00949 if (Mb->lc[lcid].pd)
00950 rr[i] += Mb->lc[lcid].pd[0-ii-1] * Mp->lambda;
00951 if (Mb->lc[lcid+1].pd)
00952 rr[i]+= Mb->lc[lcid+1].pd[0-ii-1];
00953 }
00954 if (ii==0) rr[i]=0.0;
00955 if (ii>0) rr[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
00956 }
00957 break;
00958 }
00959
00960 case growing_mech_structure:{
00961
00962
00963
00964
00965 nne = Gtm->give_nne(eid);
00966
00967 reallocv (nne,nod);
00968
00969 Gtm->give_nodes (eid,nod);
00970
00971 acn=0;
00972
00973 for (i=0;i<nne;i++){
00974
00975
00976 ndofn = Gtm->give_ndofn (nod[i]);
00977
00978 reallocv (ndofn,ncn);
00979
00980 Gtm->give_node_code_numbers (nod[i],ncn.a);
00981
00982
00983 for (j=0;j<ndofn;j++){
00984 ii=ncn[j];
00985 if (ii<0) rr[acn+j]=Mb->dlc[lcid].get_pd(Mp->time, ii);
00986 if (ii==0) rr[acn+j]=0.0;
00987 if (ii>0) rr[acn+j]=Lsrs->lhs[lcid*Ndofm+ii-1];
00988
00989 if (Mt->nodedispl[nod[i]] && (ii < 1))
00990
00991
00992 rr[acn+j] += Mt->nodedispl[nod[i]][j];
00993 }
00994 acn += ndofn;
00995 }
00996
00997 break;
00998 }
00999
01000 default:{
01001 for (i=0;i<ndofe;i++){
01002 ii=cn[i];
01003 if (ii<0) rr[i]=Mb->lc[lcid].pd[0-ii-1];
01004 if (ii==0) rr[i]=0.0;
01005 if (ii>0) rr[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
01006 }
01007 }
01008 }
01009
01010
01011
01012
01013
01014
01015 ndofe = Mt->give_ndofe (eid);
01016
01017
01018 ndofemn = Gtm->give_ndofe (eid);
01019
01020 if (ndofe != ndofemn){
01021 if (Mp->homog!=3){
01022
01023
01024
01025 mxv (Mt->elements[eid].tmat->a,rr.a,r,ndofe,ndofemn);
01026 }
01027 else{
01028
01029
01030
01031 copyv (rr.a,r,ndofe);
01032 }
01033 }
01034 else{
01035
01036 copyv (rr,r);
01037 }
01038
01039 if (Mp->tprob == growing_mech_structure){
01040
01041 Mt->elements[eid].subtrinitdispl (r, ndofe);
01042 }
01043
01044
01045 destrv (rr);
01046 destrv (cn);
01047 }
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062 void elprdispl (long lcid,long eid,double *r)
01063 {
01064 long i,ii,ndofe,ndofemn;
01065 ivector cn,nod,ncn;
01066 vector rr;
01067
01068
01069
01070 ndofe = Mt->give_ndofe (eid);
01071
01072
01073 ndofemn = Gtm->give_ndofe (eid);
01074
01075 if (ndofe != ndofemn){
01076
01077 reallocv (ndofemn,cn);
01078 reallocv (ndofemn,rr);
01079 ndofe=ndofemn;
01080 }
01081 else{
01082
01083 reallocv (ndofe,cn);
01084 reallocv (ndofe,rr);
01085 }
01086
01087 Mt->give_code_numbers (eid,cn.a);
01088
01089 switch (Mp->tprob){
01090 case forced_dynamics:
01091 case mech_timedependent_prob:{
01092
01093 for (i=0;i<ndofe;i++){
01094 ii=cn[i];
01095 if (ii<0)
01096 rr[i]=Mb->dlc[lcid].get_pd(Mp->time, ii);
01097 }
01098 break;
01099 }
01100
01101 case mat_nonlinear_statics:{
01102 for (i=0;i<ndofe;i++){
01103 ii=cn[i];
01104 if (ii<0)
01105 {
01106 rr[i]=Mb->lc[lcid].pd[0-ii-1];
01107 if ((lcid % 2) == 0)
01108 rr[i]*=Mp->lambda;
01109 }
01110 }
01111 break;
01112 }
01113
01114 case growing_mech_structure:{
01115
01116
01117
01118 for (i=0;i<ndofe;i++){
01119 ii=cn[i];
01120 if (ii<0)
01121 rr[i]=Mb->dlc[lcid].get_pd(Mp->time, ii);
01122 }
01123
01124
01125
01126
01127
01128 break;
01129 }
01130
01131 default:{
01132 for (i=0;i<ndofe;i++){
01133 ii=cn[i];
01134 if (ii<0)
01135 rr[i]=Mb->lc[lcid].pd[0-ii-1];
01136 }
01137 }
01138 }
01139
01140
01141
01142
01143
01144
01145 ndofe = Mt->give_ndofe (eid);
01146
01147
01148 ndofemn = Gtm->give_ndofe (eid);
01149
01150 if (ndofe != ndofemn){
01151
01152 mxv (Mt->elements[eid].tmat->a,rr.a,r,ndofe,ndofemn);
01153 }
01154 else{
01155
01156 copyv (rr,r);
01157 }
01158
01159
01160 if (Mp->tprob == growing_mech_structure){
01161
01162 Mt->elements[eid].subtrinitdispl (r,ndofe);
01163 }
01164
01165 }
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181 void noddispl (long lcid,double *r, long nid)
01182 {
01183 long i,j,ndofn,ndofnm,nmn;
01184 double **rr;
01185 vector w;
01186
01187
01188 ndofn = Gtm->give_ndofn(nid);
01189
01190 if (ndofn<0){
01191
01192
01193
01194 nmn=0-ndofn;
01195
01196 rr = new double* [nmn];
01197
01198 for (i=0;i<nmn;i++){
01199
01200 ndofnm = Gtm->give_ndofn (Gtm->gnodes[nid].mnodes[i]);
01201 rr[i] = new double [ndofnm];
01202 select_noddispl (lcid,rr[i],Gtm->gnodes[nid].mnodes[i]);
01203 }
01204
01205
01206 allocv (nmn,w);
01207 Gtm->approx_weights (Gtm->gnodes[nid].masentity,nid,w,Out);
01208
01209
01210
01211 for (i=0;i<ndofnm;i++){
01212 r[i]=0.0;
01213 for (j=0;j<nmn;j++){
01214 r[i]+=rr[j][i]*w[j];
01215 }
01216 }
01217
01218 for (i=0;i<0-ndofn;i++){
01219 delete [] rr[i];
01220 }
01221 delete [] rr;
01222
01223 ndofn = ndofnm;
01224 }else{
01225
01226
01227
01228
01229 select_noddispl (lcid,r,nid);
01230 }
01231
01232
01233 }
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 void select_noddispl (long lcid,double *r,long nid)
01247 {
01248 long i,ii;
01249
01250 switch (Mp->tprob){
01251 case forced_dynamics:
01252 case mech_timedependent_prob:
01253 {
01254 for (i=0;i<Mt->give_ndofn (nid);i++)
01255 {
01256 ii=Mt->give_dof(nid,i);
01257 if (ii<0) r[i]=Mb->dlc[lcid].get_pd(Mp->time, ii);
01258 if (ii==0) r[i]=0.0;
01259 if (ii>0) r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
01260 }
01261 break;
01262 }
01263 case growing_mech_structure:
01264 for (i=0;i<Mt->give_ndofn (nid);i++)
01265 {
01266 ii=Mt->give_dof(nid,i);
01267 if (ii<0) r[i]=Mb->dlc[lcid].get_pd(Mp->time, ii);
01268 if (ii==0) r[i]=0.0;
01269 if (ii>0) r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
01270
01271 if (Mt->nodedispl[nid] && (ii < 1))
01272
01273
01274 r[i] += Mt->nodedispl[nid][i];
01275 }
01276 break;
01277 case mat_nonlinear_statics:{
01278 for (i=0;i<Mt->give_ndofn (nid);i++){
01279 ii=Mt->give_dof(nid,i);
01280 if (ii<0)
01281 {
01282 r[i] = 0.0;
01283 if (Mb->lc[lcid].pd)
01284 r[i] += Mb->lc[lcid].pd[0-ii-1] * Mp->lambda;
01285 if (Mb->lc[lcid+1].pd)
01286 r[i]+= Mb->lc[lcid+1].pd[0-ii-1];
01287 }
01288 if (ii==0) r[i]=0.0;
01289 if (ii>0) r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
01290 }
01291 break;
01292 }
01293 default:{
01294 for (i=0;i<Mt->give_ndofn (nid);i++){
01295 ii=Mt->give_dof(nid,i);
01296 if (ii<0) r[i]=Mb->lc[lcid].pd[0-ii-1];
01297 if (ii==0) r[i]=0.0;
01298 if (ii>0) r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
01299 }
01300 }
01301 }
01302 }
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318 void macrodispl(long lcid, long ncomp, double *r)
01319 {
01320 double *d;
01321
01322 d = Lsrs->give_lhs(lcid);
01323 copyv(d+Ndofm-ncomp, r, ncomp);
01324 }
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339 void nodforce (double *f, long nid, vector &nf)
01340 {
01341 long i,j,ndofn,ndofnm,nmn;
01342 vector *ff;
01343 vector w;
01344
01345
01346 ndofn = Gtm->give_ndofn(nid);
01347
01348 if (ndofn<0)
01349 {
01350
01351
01352 nmn=0-ndofn;
01353
01354 ff = new vector[nmn];
01355
01356
01357 for (i=0;i<nmn;i++)
01358 {
01359
01360 ndofnm = Gtm->give_ndofn(Gtm->gnodes[nid].mnodes[i]);
01361 reallocv(ndofnm, ff[i]);
01362 select_nodforce (f, Gtm->gnodes[nid].mnodes[i], ff[i]);
01363 }
01364
01365
01366 allocv(nmn, w);
01367 Gtm->approx_weights(Gtm->gnodes[nid].masentity, nid, w, Out);
01368
01369
01370 for (i=0;i<ndofnm;i++)
01371 {
01372 nf[i]=0.0;
01373 for (j=0; j<nmn; j++)
01374 nf[i] += ff[j][i]*w[j];
01375 }
01376 delete [] ff;
01377
01378 ndofn = ndofnm;
01379 }
01380 else
01381 {
01382
01383 select_nodforce (f, nid, nf);
01384 }
01385 }
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399 void select_nodforce(double *f, long nid, vector &nf)
01400 {
01401 long i,ii;
01402
01403 for (i=0;i<Mt->give_ndofn(nid);i++)
01404 {
01405 ii=Mt->give_dof(nid,i);
01406 if (ii < 1)
01407 {
01408 if (Mp->reactcomp)
01409 nf[i]=Mt->nodes[nid].r[i];
01410 else
01411 nf[i]=0.0;
01412 }
01413 else
01414 {
01415 if (f != NULL)
01416 nf[i]=f[ii-1];
01417 else
01418 nf[i]=0.0;
01419 }
01420 }
01421 }
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492 void constr_matrix (long nid,long cid,matrix &lcm)
01493 {
01494 long i,idcs;
01495 crsectype crst;
01496 double tl,tu;
01497
01498 i=Gtm->lgnodes[nid].nodes[cid-1];
01499 crst = Mt->nodes[i].crst;
01500 idcs = Mt->nodes[i].idcs;
01501 tl = Mc->give_onethickness (crst,idcs);
01502 i=Gtm->lgnodes[nid].nodes[cid];
01503 crst = Mt->nodes[i].crst;
01504 idcs = Mt->nodes[i].idcs;
01505 tu = Mc->give_onethickness (crst,idcs);
01506
01507 fillm (0.0,lcm);
01508
01509 if (Mp->tlm==1){
01510 lcm[0][0] = -1.0;
01511 lcm[1][1] = -1.0;
01512 lcm[2][0] = tl/(-2.0);
01513
01514 lcm[3][0] = 1.0;
01515 lcm[4][1] = 1.0;
01516 lcm[5][0] = tu/(-2.0);
01517 }
01518
01519 if (Mp->tlm==2){
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539 lcm[0][0] = -1.0;
01540 lcm[1][1] = -1.0;
01541 lcm[2][2] = -1.0;
01542 lcm[3][1] = tl/2.0;
01543 lcm[4][0] = tl/(-2.0);
01544 lcm[5][3] = -1.0;
01545
01546 lcm[6][0] = 1.0;
01547 lcm[7][1] = 1.0;
01548 lcm[8][2] = 1.0;
01549 lcm[9][1] = tu/2.0;
01550 lcm[10][0] = tu/(-2.0);
01551 lcm[11][3] = 1.0;
01552 }
01553 }
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575 void mefel_right_hand_side (long lcid,double *rhs,double *flv)
01576 {
01577 if (flv)
01578 nullv(flv,Ndofm);
01579
01580 switch (Mp->tprob){
01581 case linear_statics:{
01582
01583
01584 nullv (rhs+lcid*Ndofm,Ndofm);
01585 Mb->lc[lcid].assemble (lcid,rhs+lcid*Ndofm,flv,1.0);
01586 break;
01587 }
01588 case eigen_dynamics:{
01589 break;
01590 }
01591 case forced_dynamics:{
01592
01593
01594 nullv (rhs+lcid*Ndofm,Ndofm);
01595
01596 Mb->dlc[lcid].assemble (lcid,rhs,flv,Ndofm,Mp->time);
01597 break;
01598 }
01599 case linear_stability:{
01600 break;
01601 }
01602 case mat_nonlinear_statics:{
01603
01604
01605 nullv (rhs+(lcid*2+0)*Ndofm, Ndofm);
01606 nullv (rhs+(lcid*2+1)*Ndofm, Ndofm);
01607 nullv (flv+(lcid*2+1)*Ndofm, Ndofm);
01608 Mb->lc[lcid*2+0].assemble (lcid*2+0,rhs+(lcid*2+0)*Ndofm,flv+(lcid*2+0)*Ndofm,1.0);
01609 Mb->lc[lcid*2+1].assemble (lcid*2+1,rhs+(lcid*2+1)*Ndofm,flv+(lcid*2+1)*Ndofm,1.0);
01610 break;
01611 }
01612 case geom_nonlinear_statics:{
01613
01614
01615 nullv (rhs+(lcid*2+0)*Ndofm, Ndofm);
01616 nullv (rhs+(lcid*2+1)*Ndofm, Ndofm);
01617 Mb->lc[lcid*2+0].assemble (lcid*2+0,rhs+(lcid*2+0)*Ndofm,flv,1.0);
01618 Mb->lc[lcid*2+1].assemble (lcid*2+1,rhs+(lcid*2+1)*Ndofm,flv,1.0);
01619 break;
01620 }
01621 case mech_timedependent_prob:{
01622
01623
01624 nullv (rhs+lcid*Ndofm, Ndofm);
01625 Mb->dlc[lcid].assemble (lcid,rhs+lcid*Ndofm,flv,Ndofm,Mp->time);
01626 break;
01627 }
01628 case growing_mech_structure:{
01629
01630
01631
01632
01633 nullv (rhs, Ndofm);
01634 Mb->dlc[lcid].assemble (0,rhs,flv,Ndofm,Mp->time);
01635 break;
01636 }
01637 case earth_pressure:{
01638
01639 break;
01640 }
01641 case layered_linear_statics:{
01642
01643 break;
01644 }
01645 case lin_floating_subdomain:{
01646
01647
01648 nullv (rhs+lcid*Ndofm,Ndofm);
01649 Mb->lc[lcid].assemble (lcid,rhs+lcid*Ndofm,flv,1.0);
01650 break;
01651 }
01652 case nonlin_floating_subdomain:{
01653
01654
01655 nullv (rhs+(lcid*2+0)*Ndofm, Ndofm);
01656 nullv (rhs+(lcid*2+1)*Ndofm, Ndofm);
01657 Mb->lc[lcid*2+0].assemble (lcid*2+0,rhs+(lcid*2+0)*Ndofm,flv,1.0);
01658 Mb->lc[lcid*2+1].assemble (lcid*2+1,rhs+(lcid*2+1)*Ndofm,flv,1.0);
01659 break;
01660 }
01661
01662 case hemivar_inequal:{
01663
01664
01665 nullv (rhs+lcid*Ndofm,Ndofm);
01666 Mb->lc[lcid].assemble (lcid,rhs+lcid*Ndofm,flv,1.0);
01667 break;
01668 }
01669 case load_balancing:{
01670
01671
01672 nullv (rhs+lcid*Ndofm,Ndofm);
01673 Mb->lc[lcid].assemble (lcid,rhs+lcid*Ndofm,flv,1.0);
01674 break;
01675 }
01676
01677 default:{
01678 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
01679 }
01680 }
01681
01682 if (Mp->eigstrains > 0){
01683
01684
01685
01686
01687 Mm->est=eigstrain;
01688
01689
01690 nodal_eigstrain_forces (lcid,rhs+lcid*Ndofm,Mp->time);
01691 }
01692
01693 if ((Mp->pore_press > 0) && (Mp->pore_press < 3))
01694 {
01695
01696
01697 nodal_pore_press_forces (lcid, rhs+lcid*Ndofm, Mp->time);
01698 }
01699 }
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714 void compute_req_val (long lcid)
01715 {
01716 long aux;
01717
01718
01719 if (Mp->straincomp == 1)
01720 {
01721 if (Mp->strainstate==0)
01722 {
01723 compute_ipstrains (lcid);
01724
01725
01726 Mp->strainstate=1;
01727 }
01728 if (Mp->strainpos == 2)
01729 compute_nodestrains (lcid);
01730 }
01731
01732
01733 if (Mp->stresscomp == 1)
01734 {
01735
01736 if (Mp->strainstate==0)
01737 {
01738 compute_ipstrains (lcid);
01739
01740
01741 Mp->strainstate=1;
01742 }
01743 if (Mp->stressstate==0)
01744 {
01745 compute_ipstresses (lcid);
01746
01747
01748 Mp->stressstate=1;
01749
01750
01751 Mp->otherstate=1;
01752 }
01753 if (Mp->stresspos == 2)
01754 compute_nodestresses (lcid);
01755 }
01756
01757
01758 if ((Mp->othercomp == 1) && (Mp->otherpos == 2) && (Mp->otherstate==1))
01759 compute_nodeothers ();
01760
01761
01762 if (Mp->reactcomp==1){
01763
01764 if (Mp->strainstate==0)
01765 {
01766 compute_ipstrains (lcid);
01767
01768
01769 Mp->strainstate=1;
01770 }
01771 if (Mp->stressstate==0)
01772 {
01773 compute_ipstresses (lcid);
01774
01775
01776 Mp->stressstate=1;
01777
01778
01779 Mp->otherstate=1;
01780 }
01781 aux = Mp->strcomp;
01782 Mp->strcomp = 0;
01783
01784 switch (Mp->tprob){
01785 case linear_statics:
01786 case load_balancing:
01787 case mat_nonlinear_statics:
01788 case earth_pressure:
01789 case lin_floating_subdomain:
01790 case nonlin_floating_subdomain:
01791 Mb->lc[lcid].compute_reactions (lcid);
01792 break;
01793 case forced_dynamics:
01794 case mech_timedependent_prob:
01795 case growing_mech_structure:
01796 Mb->dlc[lcid].compute_reactions(lcid);
01797 break;
01798 case eigen_dynamics:
01799 break;
01800 default:
01801 print_err("unknown type of problem is required",__FILE__,__LINE__,__func__);
01802 }
01803 Mp->strcomp = aux;
01804 }
01805
01806
01807
01808 if (Mp->straincomp == 1){
01809 if (Mp->strainstate==0)
01810 {
01811 compute_ipstrains (lcid);
01812
01813
01814 Mp->strainstate=1;
01815 }
01816
01817
01818
01819 Mm->totstrains ();
01820 if (Mp->strainpos == 2)
01821 compute_nodestrains (lcid);
01822 }
01823 }
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837 void local_global_displ_transf (long lcid)
01838 {
01839 long i,j,k,ndofn;
01840 long *cn;
01841 vector l,g;
01842 matrix tm(3,3);
01843
01844 for (i=0;i<Mt->nn;i++){
01845 if (Mt->nodes[i].transf>0){
01846
01847
01848
01849 ndofn = Mt->give_ndofn (i);
01850 allocv (ndofn,l);
01851 allocv (ndofn,g);
01852
01853 noddispl (lcid,l.a,i);
01854
01855 tm[0][0]=Mt->nodes[i].e1[0];
01856 tm[1][0]=Mt->nodes[i].e1[1];
01857 tm[2][0]=Mt->nodes[i].e1[2];
01858
01859 tm[0][1]=Mt->nodes[i].e2[0];
01860 tm[1][1]=Mt->nodes[i].e2[1];
01861 tm[2][1]=Mt->nodes[i].e2[2];
01862
01863 tm[0][2]=Mt->nodes[i].e3[0];
01864 tm[1][2]=Mt->nodes[i].e3[1];
01865 tm[2][2]=Mt->nodes[i].e3[2];
01866
01867 mxv(tm,l,g);
01868
01869 cn = new long [ndofn];
01870 Mt->give_node_code_numbers (i,cn);
01871
01872 for (j=0;j<ndofn;j++){
01873 k=cn[j]-1;
01874 if (k>-1)
01875 Lsrs->lhs[k]=g[j];
01876 }
01877 delete [] cn;
01878 destrv (l);
01879 destrv (g);
01880 }
01881 }
01882 }
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899 void stress_initdispl(long lcid)
01900 {
01901 long i, j, ne, tnip, ncomp, ipp;
01902 matrix d;
01903 vector eps, sig;
01904
01905 compute_ipstrains(lcid);
01906 ne = Mt->ne;
01907 for (i=0;i<ne;i++)
01908 {
01909 if (Gtm->leso[i]==1)
01910 {
01911 tnip = Mt->give_tnip(i);
01912 ipp = Mt->elements[i].ipp[0][0];
01913 for (j=0; j<tnip; j++)
01914 {
01915 ncomp = Mm->ip[ipp].ncompstr;
01916 reallocm(ncomp, ncomp, d);
01917 reallocv(ncomp, eps);
01918 reallocv(ncomp, sig);
01919 Mm->matstiff(d, ipp);
01920 Mm->givestrain(lcid, ipp, eps);
01921 mxv(d, eps, sig);
01922 Mm->storestress(lcid, ipp, sig);
01923 ipp++;
01924 }
01925 }
01926 }
01927
01928
01929
01930 Mp->strainstate=1;
01931
01932
01933 Mp->stressstate=1;
01934
01935
01936 Mp->otherstate=1;
01937 }
01938