00001 #include "globmatc.h"
00002 #include "globmat.h"
00003 #include "globmatt.h"
00004 #include "global.h"
00005 #include "globalt.h"
00006 #include "globalc.h"
00007 #include "loadcase.h"
00008 #include "element.h"
00009 #include "elemswitch.h"
00010 #include "elemswitcht.h"
00011 #include "elemswitchc.h"
00012 #include "intpoints.h"
00013 #include "node.h"
00014 #include "element.h"
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 void zero_order_matrix (long lcid)
00026 {
00027 long i,mndofe,tndofe;
00028 ivector rcn,ccn;
00029 matrix lm;
00030
00031 if (D0mat==NULL) D0mat = new gmatrix;
00032
00033
00034
00035
00036 D0mat->setval (Cp->ssle);
00037 D0mat->initiate (Gtu,Ndofc,Cp->tstord0,Mesprc,Outc);
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 for (i=0;i<Mt->ne;i++){
00052 mndofe = Gtm->give_ndofe (i);
00053
00054 allocm (mndofe,mndofe,lm);
00055 stiffmat (i,lcid,lm);
00056
00057 allocv (mndofe,rcn);
00058
00059 Gtm->give_code_numbers (i,rcn.a);
00060
00061 D0mat->localize (lm,rcn,i);
00062
00063 destrm (lm);
00064 destrv (rcn);
00065
00066 }
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 for (i=0;i<Tt->ne;i++){
00077 tndofe = Gtt->give_ndofe (i);
00078
00079 allocm (tndofe,tndofe,lm);
00080 conductmat (i,lcid,lm);
00081
00082 allocv (tndofe,rcn);
00083
00084 Gtt->give_code_numbers (i,rcn.a);
00085
00086 D0mat->localize (lm,rcn,i);
00087
00088 destrm (lm);
00089 destrv (rcn);
00090
00091 }
00092
00093
00094
00095
00096
00097 for (i=0;i<Ct->ne;i++){
00098 mndofe = Gtm->give_ndofe (i);
00099 tndofe = Gtt->give_ndofe (i);
00100
00101 allocm (mndofe,tndofe,lm);
00102 upper_cond_coupl_mat (i,lcid,lm);
00103
00104 allocv (mndofe,rcn);
00105
00106 allocv (tndofe,ccn);
00107
00108 Gtm->give_code_numbers (i,rcn.a);
00109 Gtt->give_code_numbers (i,ccn.a);
00110
00111 D0mat->glocalize (lm,rcn,ccn);
00112
00113 destrm (lm);
00114 destrv (rcn);
00115
00116 destrv (ccn);
00117
00118 }
00119
00120
00121
00122
00123
00124 for (i=0;i<Ct->ne;i++){
00125 mndofe = Gtm->give_ndofe (i);
00126 tndofe = Gtt->give_ndofe (i);
00127
00128 allocm (tndofe,mndofe,lm);
00129 lower_cond_coupl_mat (i,lcid,lm);
00130
00131 allocv (tndofe,rcn);
00132
00133 allocv (mndofe,ccn);
00134
00135 Gtm->give_code_numbers (i,ccn.a);
00136 Gtt->give_code_numbers (i,rcn.a);
00137
00138 D0mat->glocalize (lm,rcn,ccn);
00139
00140 destrm (lm);
00141 destrv (rcn);
00142
00143 destrv (ccn);
00144
00145 }
00146
00147
00148
00149
00150
00151 D0mat->prepmat (Cp->zero,Mesprc);
00152 }
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 void first_order_matrix (long lcid)
00164 {
00165 long i,mndofe,tndofe;
00166 ivector rcn,ccn;
00167 matrix lm;
00168
00169 if (D1mat==NULL) D1mat = new gmatrix;
00170
00171
00172
00173 D1mat->setval (Cp->ssle);
00174 D1mat->initiate (Gtu,Ndofc,Cp->tstord1,Mesprc,Outc);
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198 for (i=0;i<Tt->ne;i++){
00199 tndofe = Gtt->give_ndofe (i);
00200
00201 allocm (tndofe,tndofe,lm);
00202 capacmat (i,lcid,lm);
00203
00204 allocv (tndofe,rcn);
00205
00206 Gtt->give_code_numbers (i,rcn.a);
00207
00208 D1mat->localize (lm,rcn,i);
00209
00210 destrm (lm);
00211 destrv (rcn);
00212
00213 }
00214
00215
00216
00217
00218 for (i=0;i<Ct->ne;i++){
00219 mndofe = Gtm->give_ndofe (i);
00220 tndofe = Gtt->give_ndofe (i);
00221
00222 allocm (mndofe,tndofe,lm);
00223 upper_cap_coupl_mat (i,lcid,lm);
00224
00225 allocv (mndofe,rcn);
00226
00227 allocv (tndofe,ccn);
00228
00229 Gtm->give_code_numbers (i,rcn.a);
00230 Gtt->give_code_numbers (i,ccn.a);
00231
00232 D1mat->glocalize (lm,rcn,ccn);
00233
00234 destrm (lm);
00235 destrv (rcn);
00236
00237 destrv (ccn);
00238
00239 }
00240
00241
00242
00243 for (i=0;i<Ct->ne;i++){
00244 mndofe = Gtm->give_ndofe (i);
00245 tndofe = Gtt->give_ndofe (i);
00246
00247 allocm (tndofe,mndofe,lm);
00248 lower_cap_coupl_mat (i,lcid,lm);
00249
00250 allocv (tndofe,rcn);
00251
00252 allocv (mndofe,ccn);
00253
00254 Gtm->give_code_numbers (i,ccn.a);
00255 Gtt->give_code_numbers (i,rcn.a);
00256
00257 D1mat->glocalize (lm,rcn,ccn);
00258
00259 destrm (lm);
00260 destrv (rcn);
00261
00262 destrv (ccn);
00263
00264 }
00265
00266
00267
00268
00269 D1mat->prepmat (Cp->zero,Mesprc);
00270 }
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 void mefel_trfel_ip_mapping(long *tm_ip, long *mt_ip)
00288 {
00289 elemtypet ett;
00290 long ippt, ippm;
00291 long iokm;
00292 long i, j, ii, jj, n;
00293
00294 if (Mt->ne != Tt->ne)
00295 {
00296 print_err("number of elements of TRFEL(=%ld) and MEFEL(=%ld)\n meshes are not the same",
00297 __FILE__, __LINE__, __func__, Tt->ne, Mt->ne);
00298 abort();
00299 }
00300
00301
00302 for(i=0; i<Tm->tnip; i++)
00303 tm_ip[i] = -1;
00304
00305
00306 for(i=0; i<Mm->tnip; i++)
00307 mt_ip[i] = -1;
00308
00309 for(i=0; i<Tt->ne; i++)
00310 {
00311 ett = Tt->give_elem_type(i);
00312 switch (ett)
00313 {
00314 case barlint:
00315 for (ii=0;ii<Tp->ntm;ii++)
00316 {
00317 for (jj=0;jj<Tp->ntm;jj++)
00318 {
00319 ippt = Tt->elements[i].ipp[ii][jj];
00320 n = Tt->give_nip(i, ii, jj);
00321 for(j=0; j < n; j++)
00322 {
00323 tm_ip[ippt] = Mt->elements[i].ipp[0][0];
00324 ippt++;
00325 }
00326 }
00327 }
00328 ippm = Mt->elements[i].ipp[0][0];
00329 for(j=0; j<Mt->give_tnip(i); j++)
00330 {
00331 mt_ip[ippm] = Tt->elements[i].ipp[0][0]+1;
00332 ippm++;
00333 }
00334 break;
00335 case trlint:
00336 case trlaxisym:
00337 for (ii=0;ii<Tp->ntm;ii++)
00338 {
00339 for (jj=0;jj<Tp->ntm;jj++)
00340 {
00341 ippt = Tt->elements[i].ipp[ii][jj];
00342 n = Tt->give_nip(i, ii, jj);
00343 for(j=0; j < n; j++)
00344 {
00345 tm_ip[ippt] = Mt->elements[i].ipp[0][0];
00346 ippt++;
00347 }
00348 }
00349 }
00350 ippm = Mt->elements[i].ipp[0][0];
00351 for(j=0; j<Mt->give_tnip(i); j++)
00352 {
00353 mt_ip[ippm] = Tt->elements[i].ipp[0][0];
00354 ippm++;
00355 }
00356 break;
00357 case quadlint:
00358 for (ii=0;ii<Tp->ntm;ii++)
00359 {
00360 for (jj=0;jj<Tp->ntm;jj++)
00361 {
00362 ippt = Tt->elements[i].ipp[ii][jj];
00363 n = Tt->give_nip(i, ii, jj);
00364 iokm = Tt->give_intordkm(i, ii, jj);
00365 for(j=0; j < n; j++)
00366 {
00367 tm_ip[ippt] = Mt->elements[i].ipp[0][0]+(j%iokm);
00368 ippt++;
00369 }
00370 }
00371 }
00372 ippm = Mt->elements[i].ipp[0][0];
00373 n = Mt->give_nip(i, 0, 0);
00374 for(j=0; j<n; j++)
00375 {
00376 mt_ip[ippm] = Tt->elements[i].ipp[0][0]+j;
00377 ippm++;
00378 }
00379 break;
00380 case barquadt:
00381 case quadquadt:
00382 case quadquadtax:
00383 case lineartett:
00384 case linearhext:
00385 case quadratichext:
00386 for (ii=0;ii<Tp->ntm;ii++)
00387 {
00388 for (jj=0;jj<Tp->ntm;jj++)
00389 {
00390 ippt = Tt->elements[i].ipp[ii][jj];
00391 n = Tt->give_nip(i, ii, jj);
00392 iokm = Tt->give_intordkm(i, ii, jj);
00393 for(j=0; j < n; j++)
00394 {
00395 tm_ip[ippt] = Mt->elements[i].ipp[0][0]+(j%iokm);
00396 ippt++;
00397 }
00398 }
00399 }
00400 ippm = Mt->elements[i].ipp[0][0];
00401 for(j=0; j<Mt->give_tnip(i); j++)
00402 {
00403 mt_ip[ippm] = Tt->elements[i].ipp[0][0]+j;
00404 ippm++;
00405 }
00406 break;
00407 default:
00408 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
00409 }
00410 }
00411 }
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424 void mefel_trfel_copyip(void)
00425 {
00426 long i,j;
00427 double ipval;
00428 long antq[Tm->tnkntq];
00429
00430
00431 for(i=0; i<Tm->tnip; i++)
00432 {
00433 memset(antq, 0, sizeof(*antq)*Tm->tnkntq);
00434 Tm->give_reqntq(i, antq);
00435
00436 for (j=0; j<Tm->nntq; j++)
00437 {
00438 if (antq[int(Tm->ntqo[j])-1] == 1)
00439 {
00440 ipval = Mm->givemechq(Tm->ntqo[j], TM_ip[i]);
00441 Tm->storenontransq(Tm->ntqo[j], i, ipval);
00442 }
00443 else
00444 Tm->storenontransq(Tm->ntqo[j], i, 0.0);
00445 }
00446 }
00447 }
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460 void trfel_mefel_copyip(void)
00461 {
00462 long i,j;
00463 double ipval;
00464 long anmq[Mm->tnknmq];
00465
00466
00467 for(i=0; i<Mm->tnip; i++)
00468 {
00469 memset(anmq, 0, sizeof(*anmq)*Mm->tnknmq);
00470 for (j=0; j<Mm->ip[i].nm; j++)
00471 Mm->give_reqnmq(i, j, anmq);
00472
00473 for (j=0; j<Mm->nnmq; j++)
00474 {
00475 if (anmq[int(Mm->nmqo[j])-1] == 1)
00476 {
00477 if (MT_ip[i] < 0)
00478 {
00479 ipval = Tm->givetransq(Mm->nmqo[j], MT_ip[i-4]);
00480 ipval += Tm->givetransq(Mm->nmqo[j], MT_ip[i-3]);
00481 ipval += Tm->givetransq(Mm->nmqo[j], MT_ip[i-2]);
00482 ipval += Tm->givetransq(Mm->nmqo[j], MT_ip[i-1]);
00483 ipval /= 4.0;
00484 }
00485 else
00486 ipval = Tm->givetransq(Mm->nmqo[j], MT_ip[i]);
00487 Mm->storenonmechq(Mm->nmqo[j], i, ipval);
00488 }
00489 else
00490 Mm->storenonmechq(Mm->nmqo[j], i, 0.0);
00491 }
00492 }
00493 }
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505 void mefel_trfel_by_nodes(void)
00506 {
00507 long i, j, k, nnet, nnem, nnv, nip, ipid, nentq;
00508 ivector nodes;
00509 vector nodval, ipval, auxnv;
00510 long antq[Tm->tnkntq];
00511 nontransquant *entqo;
00512
00513 for (i=0;i<Tt->ne;i++)
00514 {
00515 if (Gtt->leso[i]==1)
00516 {
00517
00518 memset(antq, 0, sizeof(*antq)*Tm->tnkntq);
00519 ipid=Tt->elements[i].ipp[0][0];
00520 for(j=0; j<Tt->give_tnip(i); j++)
00521 {
00522 Tm->give_reqntq(ipid, antq);
00523 ipid++;
00524 }
00525
00526 for (j=0, nentq=0; j<Tm->tnkntq; j++)
00527 {
00528 if (antq[j] == 1) nentq++;
00529 }
00530
00531 entqo = new nontransquant[nentq];
00532 for (j=0, k=0; j<Tm->tnkntq; j++)
00533 {
00534 if (antq[j] == 1)
00535 entqo[k]= nontransquant(j+1);
00536 }
00537
00538 nnet = Tt->give_nne (i);
00539 nnem = Mt->give_nne (i);
00540 nnv = (nnet < nnem) ? nnet : nnem;
00541 reallocv (nnv*Tm->nntq,nodval);
00542
00543
00544 elem_mechq_nodval_comp(i, nodval, nnv, nentq, entqo);
00545
00546
00547 auxnv.n = nnet;
00548 nip = Tt->give_tnip (i);
00549 reallocv (nip,ipval);
00550 for (j=0; j<nentq; j++)
00551 {
00552 auxnv.a = nodval.a + j*nnv;
00553
00554 if (nnet > nnem)
00555 {
00556
00557
00558
00559
00560 print_err("TRFEL element approximation order is higher than on MEFEL element", __FILE__, __LINE__, __func__);
00561 abort();
00562 }
00563 else
00564 {
00565
00566
00567
00568
00569
00570
00571
00572
00573 elem_intpointvalt (i, auxnv, ipval);
00574 }
00575
00576
00577 ipid=Tt->elements[i].ipp[0][0];
00578
00579 for (k=0;k<nip;k++)
00580 {
00581 Tm->storenontransq(entqo[j], ipid, ipval[k]);
00582 ipid++;
00583 }
00584 }
00585 }
00586 }
00587 auxnv.a = NULL;
00588 }
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601 void internal_gforces (long lcid,double *intfor)
00602 {
00603 double *ifort,*iform;
00604
00605 ifort= new double [Ndofc];
00606 nullv (ifort,Ndofc);
00607 iform= new double [Ndofc];
00608 nullv (iform,Ndofc);
00609 nullv (intfor,Ndofc);
00610
00611
00612 internal_forces (lcid,iform);
00613 addv (intfor,iform,Ndofc);
00614
00615
00616 internal_fluxes (ifort,Ndofc);
00617 addv (intfor,ifort,Ndofc);
00618
00619
00620
00621
00622
00623 delete [] iform;
00624 delete [] ifort;
00625 }
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638 void incr_internal_gforces (long lcid,double *intfor)
00639 {
00640 double *ifort,*iform;
00641
00642 ifort= new double [Ndofc];
00643 nullv (ifort,Ndofc);
00644 iform= new double [Ndofc];
00645 nullv (iform,Ndofc);
00646 nullv (intfor,Ndofc);
00647
00648
00649 incr_internal_forces (lcid,iform);
00650 addv (intfor,iform,Ndofc);
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660 delete [] iform;
00661 delete [] ifort;
00662 }
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676 void elemvaluesc (long lcid,double *r,long *cn,long ndofe)
00677 {
00678 long i,ii;
00679
00680 for (i=0;i<ndofe;i++){
00681 ii=cn[i];
00682 if (ii<0) r[i]=Tb->lc[lcid].pv[0-ii-1].getval();
00683 if (ii==0) r[i]=0.0;
00684 if (ii>0) r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00685 }
00686 }
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702 void metr_right_hand_side (long lcid,double *rhs,double *flv)
00703 {
00704 double *av;
00705
00706 av = new double [Ndofc];
00707
00708 nullv (rhs,Ndofc);
00709 nullv (flv,Ndofc);
00710 nullv (av,Ndofc);
00711
00712 mefel_right_hand_side (lcid,rhs,flv);
00713 trfel_right_hand_side (lcid,av,Ndofc);
00714 addv (rhs,av,Ndofc);
00715
00716 nullv (av,Ndofc);
00717 right_hand_side (lcid,av,Ndofc);
00718 addv (rhs,av,Ndofc);
00719
00720 delete [] av;
00721 }
00722
00723
00724
00725
00726
00727
00728
00729
00730 void updateval()
00731 {
00732
00733 Tm->updateipval ();
00734 Mm->updateipval();
00735 Cmu->updateipval();
00736
00737 }
00738
00739
00740
00741
00742
00743
00744
00745
00746 void copy_data()
00747 {
00748 if (Cp->tprob==fully_coupled_mech_trans){
00749
00750
00751
00752
00753 }
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763 trfel_mefel ();
00764
00765 mefel_trfel ();
00766
00767 }
00768
00769
00770
00771
00772
00773
00774
00775
00776 void approximationc ()
00777 {
00778
00779 approximation ();
00780
00781
00782 approximationcoup ();
00783 }
00784
00785
00786
00787
00788
00789
00790
00791
00792 void approximationcoup ()
00793 {
00794 long i;
00795
00796 switch (Cp->tmatt){
00797 case mech_onemedium:{
00798 for (i=0;i<Ct->ne;i++){
00799 intpointvaluesc (i);
00800 intpointstrainsc (i);
00801 intpointgradientsc (i);
00802 }
00803 break;
00804 }
00805 case mech_twomedia:{
00806 for (i=0;i<Ct->ne;i++){
00807 intpointvaluesc (i);
00808 intpointstrainsc (i);
00809 intpointgradientsc (i);
00810 }
00811 break;
00812 }
00813 case mech_threemedia:{
00814 for (i=0;i<Ct->ne;i++){
00815 intpointvaluesc (i);
00816 intpointstrainsc (i);
00817 intpointgradientsc (i);
00818 }
00819 break;
00820 }
00821 default:{
00822 print_err("unknown transported matter is required",__FILE__,__LINE__,__func__);
00823 abort();
00824 }
00825 }
00826
00827 }
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839 void internal_coup_fluxes (double *intflux,long n)
00840 {
00841 long i,ne,ndofe,*cn;
00842 elemtypec te;
00843 vector lif;
00844
00845 ne = Ct->ne;
00846
00847 for (i=0;i<ne;i++){
00848 te = Ct->give_elem_type (i);
00849 ndofe = Gtt->give_ndofe (i);
00850 allocv (ndofe,lif);
00851 fillv (0.0,lif);
00852
00853 cn = new long [ndofe];
00854 Gtt->give_code_numbers (i,cn);
00855
00856 switch (te){
00857 case coupbar:{
00858 Cbar->res_lower_internal_fluxes (i,lif);
00859 break;
00860 }
00861 default:{
00862 fprintf (stderr,"\n\n unknown element type is required in function");
00863 fprintf (stderr,"\n internal_coup_fluxes (file %s, line %d).\n",__FILE__,__LINE__);
00864 abort ();
00865 }
00866 }
00867
00868 locglob (intflux,lif.a,cn,ndofe);
00869 destrv (lif);
00870 delete [] cn;
00871
00872 }
00873
00874
00875 }
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886 void internal_coup_forces (long lcid,double *intflux)
00887 {
00888 long i,ne,ndofe,*cn;
00889 elemtypec te;
00890 vector lif;
00891
00892 ne = Ct->ne;
00893
00894 for (i=0;i<ne;i++){
00895 te = Ct->give_elem_type (i);
00896 ndofe = Gtm->give_ndofe (i);
00897 allocv (ndofe,lif);
00898 fillv (0.0,lif);
00899
00900 cn = new long [ndofe];
00901 Gtm->give_code_numbers (i,cn);
00902
00903 switch (te){
00904 case coupbar:{
00905 Cbar->res_upper_internal_forces (i,lif);
00906 break;
00907 }
00908 default:{
00909 fprintf (stderr,"\n\n unknown element type is required in function");
00910 fprintf (stderr,"\n internal_coup_forces (file %s, line %d).\n",__FILE__,__LINE__);
00911 abort ();
00912 }
00913 }
00914
00915 locglob (intflux,lif.a,cn,ndofe);
00916 destrv (lif);
00917 delete [] cn;
00918
00919 }
00920
00921
00922
00923 }
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935 void right_hand_side (long lcid,double *rhs,long n)
00936 {
00937 long i,mndofe,*cn;
00938 vector lv;
00939 double *av;
00940 av = new double [n];
00941
00942 assemble_init_coup_upper (rhs);
00943
00944 switch (Cp->tmatt){
00945 case mech_onemedium:{
00946 nullv (av,n);
00947 assemble_coup (lcid,av,n);
00948 addv (rhs,av,n);
00949 break;
00950 }
00951 case mech_twomedia:{
00952 nullv (av,n);
00953 assemble_coup (lcid,av,n);
00954 addv (rhs,av,n);
00955 break;
00956 }
00957 case mech_threemedia:{
00958 nullv (av,n);
00959 assemble_coup (lcid,av,n);
00960 addv (rhs,av,n);
00961 break;
00962 }
00963 default:{
00964 fprintf (stderr,"\n\n unknown transported matter is required in function");
00965 fprintf (stderr,"\n right_hand_side (file %s, line %d).\n",__FILE__,__LINE__);
00966 }
00967 }
00968
00969
00970 nullv (av,n);
00971 for (i=0;i<Ct->ne;i++){
00972 if (Gtm->leso[i]==1){
00973
00974
00975
00976 mndofe=Gtm->give_ndofe (i);
00977 allocv (mndofe,lv);
00978 fillv (0.0,lv);
00979
00980
00981 volume_rhs_vectorc (lv,lcid,i);
00982
00983
00984
00985 cn = new long [mndofe];
00986 Gtm->give_code_numbers (i,cn);
00987
00988
00989 locglob (av,lv.a,cn,mndofe);
00990
00991 destrv (lv);
00992 delete [] cn;
00993 }
00994 }
00995
00996 addv (rhs,av,n);
00997 delete [] av;
00998 }
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009 void assemble_init_coup_upper (double *rhs)
01010 {
01011 long i,mndofe,tndofe,lcid=0;
01012 long *rcn,*ccn;
01013 vector r,f;
01014 matrix km;
01015
01016
01017 for (i=0;i<Ct->ne;i++){
01018 mndofe = Gtm->give_ndofe (i);
01019 tndofe = Gtt->give_ndofe (i);
01020
01021 allocv (tndofe,r);
01022 allocv (mndofe,f);
01023 allocm (mndofe,tndofe,km);
01024
01025 upper_cond_coupl_mat (i,lcid,km);
01026
01027 rcn = new long [mndofe];
01028 ccn = new long [tndofe];
01029 Gtm->give_code_numbers (i,rcn);
01030 Gtt->give_code_numbers (i,ccn);
01031
01032 elemvalues (0,i,r.a,ccn,tndofe);
01033
01034 mxv (km,r,f);
01035 cmulv (-1.0,f);
01036
01037 locglob (rhs,f.a,rcn,mndofe);
01038
01039 destrm (km); destrv (f); destrv (r);
01040 delete [] rcn;
01041 delete [] ccn;
01042 }
01043 }
01044
01045
01046
01047
01048
01049
01050
01051
01052 void approximation_temper ()
01053 {
01054 long nnm;
01055 double *gv;
01056
01057
01058 nnm = Mt->nn;
01059
01060 gv = new double [nnm];
01061 nullv (gv,nnm);
01062
01063
01064 nodal_nodal_values (gv,temperature);
01065
01066
01067 if (Cp->lbb==quad_lin){
01068 intpointval2 (gv,temperature);
01069 }
01070 if (Cp->lbb==quad_quad || Cp->lbb==lin_lin){
01071 intpointval (gv,temperature,1.0);
01072 }
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082 delete [] gv;
01083 }
01084
01085
01086
01087
01088
01089
01090
01091
01092 void approximation_humid ()
01093 {
01094 long nnm;
01095 double *gv;
01096
01097
01098 nnm = Mt->nn;
01099
01100 gv = new double [nnm];
01101 nullv (gv,nnm);
01102
01103
01104 nodal_nodal_values (gv,rel_hum);
01105
01106
01107 if (Cp->lbb==quad_lin){
01108 intpointval2 (gv,rel_hum);
01109 }
01110 if (Cp->lbb==quad_quad || Cp->lbb==lin_lin){
01111 intpointval (gv,rel_hum,1.0);
01112 }
01113
01114 delete [] gv;
01115 }
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127 void nodal_nodal_values (double *gv, nonmechquant nmq)
01128 {
01129 long i,nnt,cn, idof;
01130
01131 nnt = Tt->nn;
01132
01133 switch (nmq){
01134
01135 case temperature:{
01136 idof = Tp->give_var_dofid(trf_temperature);
01137 if (idof < 0)
01138 {
01139 print_err("temperature is not defined in the problem", __FILE__, __LINE__, __func__);
01140 abort();
01141 }
01142 for (i=0;i<nnt;i++)
01143 {
01144 cn=Tt->give_dof (i,idof);
01145 if (cn>0) gv[i] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01146 if (cn<0) gv[i] = Tb->lc[0].pv[0-cn-1].getval();
01147 }
01148 break;
01149 }
01150
01151 case initial_temperature:{
01152 idof = Tp->give_var_dofid(trf_temperature);
01153 if (idof < 0)
01154 {
01155 print_err("initial temperature is not defined in the problem", __FILE__, __LINE__, __func__);
01156 abort();
01157 }
01158 for (i=0;i<nnt;i++)
01159 {
01160 cn=Tt->give_dof (i,idof);
01161 if (cn>0) gv[i] = Lsrst->lhsi[cn-1];
01162 if (cn<0) gv[i] = Tb->lc[0].pv[0-cn-1].getval();
01163 }
01164 break;
01165 }
01166
01167 case water_pressure:{
01168 idof = Tp->give_var_dofid(trf_water_press);
01169 if (idof < 0)
01170 {
01171 print_err("water pressure is not defined in the problem", __FILE__, __LINE__, __func__);
01172 abort();
01173 }
01174 for (i=0;i<nnt;i++)
01175 {
01176 cn=Tt->give_dof (i,idof);
01177 if (cn>0) gv[i] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01178 if (cn<0) gv[i] = Tb->lc[0].pv[0-cn-1].getval();
01179 }
01180 break;
01181 }
01182 case saturation_degree:{
01183
01184 give_transq_nodval (gv,nnt,nmq);
01185 break;
01186 }
01187
01188 case suction:{
01189
01190 give_transq_nodval (gv,nnt,nmq);
01191 break;
01192 }
01193
01194 case rel_hum:{
01195 give_nodal_humid(gv,nnt);
01196 break;
01197 }
01198 case pore_pressure:{
01199
01200
01201
01202 give_transq_nodval (gv,nnt,nmq);
01203 break;
01204 }
01205 case vol_moist_cont:{
01206 give_transq_nodval(gv,nnt,nmq);
01207 break;
01208 }
01209
01210 default:{
01211 print_err("unknown type of nmq %d is required",__FILE__,__LINE__,__func__, nmq);
01212 }
01213 }
01214 }
01215
01216
01217
01218
01219
01220
01221
01222
01223 void approximation_inittemper ()
01224 {
01225 long nnm;
01226 double *gv;
01227
01228
01229 nnm = Mt->nn;
01230
01231 gv = new double [nnm];
01232 nullv (gv,nnm);
01233
01234
01235 nodal_nodal_values (gv,initial_temperature);
01236
01237 if (Cp->lbb==quad_quad || Cp->lbb==lin_lin){
01238 intpointval (gv,initial_temperature,1.0);
01239 }
01240 if (Cp->lbb==quad_lin){
01241 intpointval2 (gv,initial_temperature);
01242 }
01243
01244 delete [] gv;
01245 }
01246
01247
01248
01249
01250
01251
01252
01253
01254 void metr_mefel ()
01255 {
01256 }
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266 void mefel_trfel ()
01267 {
01268 long i, j, k, nne, nip, ipid;
01269 ivector nodes;
01270 vector nodval, ipval;
01271
01272
01273 for(i=0; i<Tm->nntq; i++)
01274 {
01275 for (j=0;j<Tt->ne;j++)
01276 {
01277 if (Gtt->leso[j]==1)
01278 {
01279 nne = Tt->give_nne (j);
01280 nip = Tt->give_tnip (j);
01281
01282 reallocv (nne,nodes);
01283 reallocv (nne,nodval);
01284 reallocv (nip,ipval);
01285
01286
01287 Tt->give_elemnodes (j, nodes);
01288
01289
01290
01291
01292 elem_mechq_nodval (j, nodval, Tm->ntqo[i]);
01293
01294
01295 elem_intpointvalt (j, nodval, ipval);
01296
01297
01298 ipid=Tt->elements[j].ipp[0][0];
01299
01300 for (k=0;k<nip;k++)
01301 {
01302 Tm->storenontransq(Tm->ntqo[i], ipid, ipval[k]);
01303 ipid++;
01304 }
01305 }
01306 }
01307 }
01308 }
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321 void trfel_mefel ()
01322 {
01323 long i, nnm;
01324 double *gv;
01325
01326
01327 nnm = Mt->nn;
01328
01329 gv = new double [nnm];
01330 for (i=0; i<Mm->nnmq; i++)
01331 {
01332 if (Mm->nmqo[i] == initial_temperature)
01333 {
01334
01335 continue;
01336 }
01337
01338 nullv (gv,nnm);
01339
01340 nodal_nodal_values (gv, Mm->nmqo[i]);
01341
01342 if (Cp->lbb==quad_quad || Cp->lbb==lin_lin){
01343 intpointval(gv,Mm->nmqo[i],1.0);
01344 }
01345 if (Cp->lbb==quad_lin){
01346 intpointval2(gv,Mm->nmqo[i]);
01347 }
01348 }
01349
01350 delete [] gv;
01351 }
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363 void trfel_mefel_by_nodes(void)
01364 {
01365 long i, j, k, nnet, nnem, nip, ipid, nenmq;
01366 ivector nodes;
01367 vector nodval, ipval, auxnv;
01368 long anmq[Mm->tnknmq];
01369 nonmechquant *enmqo;
01370
01371 for (i=0;i<Mt->ne;i++)
01372 {
01373 if (Gtm->leso[i]==1)
01374 {
01375 nnem = Mt->give_nne (i);
01376 nip = Mt->give_tnip (i);
01377
01378 reallocv (nnem,nodes);
01379 reallocv (nnem*Mm->nnmq,nodval);
01380 reallocv (nip,ipval);
01381
01382
01383 Mt->give_elemnodes (i, nodes);
01384
01385 memset(anmq, 0, sizeof(*anmq)*Mm->tnknmq);
01386
01387 ipid=Mt->elements[i].ipp[0][0];
01388 for(j=0; j<Mt->give_tnip(i); j++)
01389 {
01390 for (k=0; k<Mm->ip[ipid].nm; j++)
01391 Mm->give_reqnmq(ipid, k, anmq);
01392 ipid++;
01393 }
01394
01395 for (j=0, nenmq=0; j<Mm->tnknmq; j++)
01396 {
01397 if (anmq[j] == 1) nenmq++;
01398 }
01399
01400 enmqo = new nonmechquant[nenmq];
01401 for (j=0, k=0; j<Mm->tnknmq; j++)
01402 {
01403 if (anmq[j] == 1)
01404 enmqo[k]= nonmechquant(j+1);
01405 }
01406
01407
01408 elem_transq_nodval_comp(i, nodval, nnem, nenmq, enmqo);
01409
01410 nnet = Tt->give_nne (i);
01411
01412
01413 auxnv.n = nnem;
01414 for (j=0; j<nenmq; j++)
01415 {
01416 auxnv.a = nodval.a + j*nnem;
01417
01418 if (nnem > nnet)
01419 {
01420
01421
01422
01423 elem_intpointval2 (i, auxnv, ipval);
01424 }
01425 else
01426 {
01427
01428
01429
01430
01431
01432
01433
01434
01435 elem_intpointval (i, auxnv, ipval);
01436 }
01437
01438
01439
01440 ipid=Mt->elements[i].ipp[0][0];
01441
01442 for (k=0;k<nip;k++)
01443 {
01444 Mm->storenonmechq(enmqo[j], ipid, ipval[k]);
01445 ipid++;
01446 }
01447 }
01448 delete [] enmqo;
01449 }
01450 }
01451 auxnv.a = NULL;
01452 }
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464 void init_trfel_mefel ()
01465 {
01466 long i, nnm;
01467 double *gv;
01468
01469
01470 nnm = Mt->nn;
01471
01472 gv = new double [nnm];
01473 for (i=0; i<Mm->nnmq; i++)
01474 {
01475 if (Mm->nmqo[i] == initial_temperature)
01476 {
01477
01478 nullv (gv,nnm);
01479
01480
01481 nodal_nodal_values (gv, Mm->nmqo[i]);
01482
01483 if (Cp->lbb==quad_quad || Cp->lbb==lin_lin){
01484 intpointval (gv,Mm->nmqo[i],1.0);
01485 }
01486 if (Cp->lbb==quad_lin){
01487 intpointval2 (gv,Mm->nmqo[i]);
01488 }
01489 }
01490 }
01491
01492 delete [] gv;
01493 }