00001 #include <stdlib.h>
00002 #include <errno.h>
00003 #include "globmatt.h"
00004 #include "globalt.h"
00005 #include "lhsrhst.h"
00006 #include "constrel.h"
00007 #include "onemedium.h"
00008 #include "twomedia.h"
00009 #include "threemedia.h"
00010 #include "fourmedia.h"
00011 #include "elemswitcht.h"
00012
00013
00014
00015
00016
00017
00018 void conductivity_matrix (long lcid)
00019 {
00020 long i,ndofe;
00021 ivector cn;
00022 matrix lm;
00023
00024 if (Kmat==NULL) Kmat = new gmatrix;
00025
00026 Kmat->setval (Tp->ssle);
00027 Kmat->initiate (Gtt,Ndoft,Tp->tstorkm,Mesprt,Outt);
00028
00029
00030 for (i=0;i<Tt->ne;i++){
00031 if (Gtt->leso[i]==1){
00032
00033 ndofe = Gtt->give_ndofe (i);
00034
00035 allocm (ndofe,ndofe,lm);
00036
00037 conductmat (i,lcid,lm);
00038
00039 check_math_errel(i);
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 allocv (ndofe,cn);
00060
00061 Gtt->give_code_numbers (i,cn.a);
00062
00063 Kmat->localize (lm,cn,i);
00064
00065 destrm (lm);
00066
00067 destrv (cn);
00068 }
00069 }
00070
00071 Kmat->prepmat (0.0,Mesprt);
00072
00073 }
00074
00075
00076
00077
00078
00079
00080 void capacity_matrix (long lcid)
00081 {
00082 long i,ndofe;
00083 ivector cn;
00084 matrix lm;
00085
00086 if (Cmat==NULL) Cmat = new gmatrix;
00087
00088 Cmat->setval (Tp->ssle);
00089 Cmat->initiate (Gtt,Ndoft,Tp->tstorcm,Mesprt,Outt);
00090
00091
00092 for (i=0;i<Tt->ne;i++){
00093 if (Gtt->leso[i]==1){
00094
00095 ndofe = Gtt->give_ndofe (i);
00096
00097 allocm (ndofe,ndofe,lm);
00098
00099 capacmat (i,lcid,lm);
00100
00101 check_math_errel(i);
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 allocv (ndofe,cn);
00113
00114 Gtt->give_code_numbers (i,cn.a);
00115
00116 Cmat->localize (lm,cn,i);
00117
00118 destrm (lm);
00119 destrv (cn);
00120
00121 }
00122 }
00123
00124 Cmat->prepmat (0.0,Mesprt);
00125
00126 }
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137 void residuum (double *r,double *p,double *v,double dt,long n,long lcid)
00138 {
00139 capacity_matrix (0);
00140 conductivity_matrix (0);
00141 trfel_right_hand_side (0,Lsrst->rhs,Ndoft);
00142
00143 Cmat->gmxv (p,v);
00144 cmulv (1.0/dt,v,n);
00145 Kmat->gmxv (Lsrst->lhs,r);
00146 addv (r,v,n);
00147 subv (r,Lsrst->rhs,n);
00148 }
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 void elemvalues_puc (long lcid,long eid,double *r,long ndofe)
00160 {
00161 long i,j,ii,nn,ncomp,ndofn,nne,k;
00162 ivector nod;
00163 double *rd;
00164 j = 0;
00165 nn = 0;
00166
00167 ncomp = Tt->give_ncomp(eid);
00168
00169 nne = Tt->give_nne (eid);
00170 allocv (nne,nod);
00171 Tt->give_elemnodes (eid,nod);
00172
00173
00174 ii = 0;
00175 for (i=0;i<nne;i++){
00176 k = nod[i];
00177 ncomp = Tt->nodes[k].ncompgrad;
00178 ndofn=Gtt->give_ndofn (k);
00179 rd = new double [ndofn];
00180 nodalval (0, rd, k);
00181
00182 for (j=0;j<ndofn;j++){
00183 if(ncomp == 1){
00184 r[ii] = rd[j] + Tb->lc[j].masterval + Tb->lc[j].mastergrad[0]*Gtt->gnodes[k].x;
00185 }
00186 if(ncomp == 2){
00187 r[ii] = rd[j] + Tb->lc[j].masterval + Tb->lc[j].mastergrad[0]*Gtt->gnodes[k].x+Tb->lc[j].mastergrad[1]*Gtt->gnodes[k].y;
00188 }
00189 if(ncomp == 3){
00190 r[ii] = rd[j] + Tb->lc[j].masterval + Tb->lc[j].mastergrad[0]*Gtt->gnodes[k].x+Tb->lc[j].mastergrad[1]*Gtt->gnodes[k].y+Tb->lc[j].mastergrad[2]*Gtt->gnodes[k].z;
00191 }
00192 ii++;
00193 }
00194 delete [] rd;
00195 }
00196 }
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 void elemvalues (long lcid,long eid,double *r,long *cn,long ndofe)
00210 {
00211 long i,j,ii,nn,ncomp,no,nne;
00212 ivector nod;
00213 j = 0;
00214 nn = 0;
00215
00216 ncomp = Tt->give_ncomp(eid);
00217
00218 nne = Tt->give_nne (eid);
00219 allocv (nne,nod);
00220 Tt->give_elemnodes (eid,nod);
00221
00222 switch (Tp->tprob){
00223 case stationary_problem:{
00224 for (i=0;i<ndofe;i++){
00225 ii=cn[i];
00226
00227
00228
00229 if(Tp->homogt == 1){
00230 j++;
00231 if(j > Tp->ntm){
00232 j = 1;
00233 nn++;
00234 }
00235 }
00236
00237 if (ii<0){
00238 if(Tp->homogt == 1){
00239 r[i]=Tb->lc[lcid].pv[0-ii-1].getval() + Tb->lc[j-1].masterval;
00240 no = nod[nn];
00241 if(ncomp == 1){
00242 r[i]=r[i]+Tb->lc[j-1].mastergrad[0]*Gtt->gnodes[no].x;
00243 }
00244 if(ncomp == 2){
00245 r[i]=r[i]+Tb->lc[j-1].mastergrad[0]*Gtt->gnodes[no].x+Tb->lc[j-1].mastergrad[1]*Gtt->gnodes[no].y;
00246 }
00247 if(ncomp == 3){
00248 r[i]=r[i]+Tb->lc[j-1].mastergrad[0]*Gtt->gnodes[no].x+Tb->lc[j-1].mastergrad[1]*Gtt->gnodes[no].y+Tb->lc[j-1].mastergrad[2]*Gtt->gnodes[no].z;
00249 }
00250 }
00251 else
00252 r[i]=Tb->lc[lcid].pv[0-ii-1].getval();
00253 }
00254
00255 if (ii==0) r[i]=0.0;
00256 if (ii>0) r[i]=Lsrst->lhs[ii-1];
00257 }
00258 break;
00259 }
00260 case nonlinear_stationary_problem:{
00261 for (i=0;i<ndofe;i++){
00262 ii=cn[i];
00263
00264
00265
00266 if(Tp->homogt == 1){
00267 j++;
00268 if(j > Tp->ntm){
00269 j = 1;
00270 nn++;
00271 }
00272 }
00273
00274 if (ii<0){
00275 if(Tp->homogt == 1){
00276 r[i]=Tb->lc[lcid].pv[0-ii-1].getval() + Tb->lc[j-1].masterval;
00277 no = nod[nn];
00278 if(ncomp == 1){
00279 r[i]=r[i]+Tb->lc[j-1].mastergrad[0]*Gtt->gnodes[no].x;
00280 }
00281 if(ncomp == 2){
00282 r[i]=r[i]+Tb->lc[j-1].mastergrad[0]*Gtt->gnodes[no].x+Tb->lc[j-1].mastergrad[1]*Gtt->gnodes[no].y;
00283 }
00284 if(ncomp == 3){
00285 r[i]=r[i]+Tb->lc[j-1].mastergrad[0]*Gtt->gnodes[no].x+Tb->lc[j-1].mastergrad[1]*Gtt->gnodes[no].y+Tb->lc[j-1].mastergrad[2]*Gtt->gnodes[no].z;
00286 }
00287 }
00288 else
00289 r[i]=Tb->lc[lcid].pv[0-ii-1].getval();
00290 }
00291
00292 if (ii==0) r[i]=0.0;
00293 if (ii>0) r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00294 }
00295 break;
00296 }
00297 case nonstationary_problem:
00298 case nonlinear_nonstationary_problem:
00299 case discont_nonstat_problem:
00300 case discont_nonlin_nonstat_problem:
00301 case growing_np_problem:
00302 case growing_np_problem_nonlin:{
00303 for (i=0;i<ndofe;i++){
00304 ii=cn[i];
00305 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00306 if (ii==0) r[i]=0.0;
00307 if (ii>0){
00308 r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00309 }
00310 }
00311 break;
00312 }
00313 default:{
00314 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00315 }
00316 }
00317
00318 destrv (nod);
00319 }
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331 void nodalval (long lcid,double *r, long nid)
00332 {
00333 long i,ii,ndofn;
00334
00335
00336 ndofn = Tt->give_ndofn (nid);
00337
00338 switch (Tp->tprob){
00339 case stationary_problem:{
00340 for (i=0;i<ndofn;i++){
00341 ii=Tt->give_dof (nid,i);
00342 if (ii<0){
00343 r[i]=Tb->lc[0].pv[0-ii-1].getval();
00344 }
00345 if (ii==0){
00346 r[i]=0.0;
00347 }
00348 if (ii>0){
00349 r[i]=Lsrst->lhs[ii-1];
00350 }
00351 }
00352 break;
00353 }
00354 case nonlinear_stationary_problem:{
00355 for (i=0;i<ndofn;i++){
00356 ii=Tt->give_dof (nid,i);
00357 if (ii<0){
00358 r[i]=Tb->lc[0].pv[0-ii-1].getval();
00359 }
00360 if (ii==0){
00361 r[i]=0.0;
00362 }
00363 if (ii>0){
00364 if(Tp->homogt == 1)
00365 r[i]=Lsrst->lhs[ii-1];
00366 else
00367 r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00368 }
00369 }
00370 break;
00371 }
00372 case nonstationary_problem:
00373 case nonlinear_nonstationary_problem:
00374 case discont_nonstat_problem:
00375 case discont_nonlin_nonstat_problem:
00376 case growing_np_problem:
00377 case growing_np_problem_nonlin:{
00378
00379 for (i=0;i<ndofn;i++){
00380 ii=Tt->give_dof (nid,i);
00381 if (ii<0){
00382 r[i]=Tb->lc[0].pv[0-ii-1].getval();
00383 }
00384 if (ii==0){
00385 r[i]=0.0;
00386 }
00387 if (ii>0){
00388 r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00389 }
00390 }
00391 break;
00392 }
00393 default:{
00394 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00395 }
00396 }
00397 }
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408 double nodalval (long lcid,long nid,long dofid)
00409 {
00410 long ii;
00411 double v=0.0;
00412
00413 switch (Tp->tprob){
00414 case stationary_problem:{
00415 ii=Tt->give_dof (nid,dofid);
00416 if (ii<0){
00417 v=Tb->lc[0].pv[0-ii-1].getval();
00418 }
00419 if (ii==0){
00420 v=0.0;
00421 }
00422 if (ii>0){
00423 v=Lsrst->lhs[ii-1];
00424 }
00425 break;
00426 }
00427 case nonlinear_stationary_problem:{
00428 ii=Tt->give_dof (nid,dofid);
00429 if (ii<0){
00430 v=Tb->lc[0].pv[0-ii-1].getval();
00431 }
00432 if (ii==0){
00433 v=0.0;
00434 }
00435 if (ii>0){
00436 if(Tp->homogt == 1)
00437 v=Lsrst->lhs[ii-1];
00438 else
00439 v=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00440 }
00441 break;
00442 }
00443 case nonstationary_problem:
00444 case nonlinear_nonstationary_problem:
00445 case discont_nonstat_problem:
00446 case discont_nonlin_nonstat_problem:
00447 case growing_np_problem:
00448 case growing_np_problem_nonlin:{
00449
00450 ii=Tt->give_dof (nid,dofid);
00451 if (ii<0){
00452 v=Tb->lc[0].pv[0-ii-1].getval();
00453 }
00454 if (ii==0){
00455 v=0.0;
00456 }
00457 if (ii>0){
00458 v=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00459 }
00460 break;
00461 }
00462 default:{
00463 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00464 }
00465 }
00466 return v;
00467 }
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480 void initnodval (long lcid,double *r, long idn)
00481 {
00482 long i,ii;
00483
00484 switch (Tp->tprob){
00485 case nonlinear_stationary_problem:
00486 case nonstationary_problem:{
00487 for (i=0;i<Tt->give_ndofn (idn);i++){
00488 ii=Tt->give_dof(idn,i);
00489 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00490 if (ii==0) r[i]=0.0;
00491 if (ii>0){
00492 r[i]=Lsrst->lhsi[ii-1];
00493 }
00494 }
00495
00496 break;
00497 }
00498 case nonlinear_nonstationary_problem:{
00499 for (i=0;i<Tt->give_ndofn (idn);i++){
00500 ii=Tt->give_dof(idn,i);
00501 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00502 if (ii==0) r[i]=0.0;
00503 if (ii>0){
00504 r[i]=Lsrst->lhsi[ii-1];
00505 }
00506 }
00507
00508 break;
00509 }
00510 case discont_nonstat_problem:
00511 case discont_nonlin_nonstat_problem:{
00512 for (i=0;i<Tt->give_ndofn (idn);i++){
00513 ii=Tt->give_dof(idn,i);
00514 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00515 if (ii==0) r[i]=0.0;
00516 if (ii>0){
00517 r[i]=Lsrst->lhsi[ii-1];
00518 }
00519 }
00520
00521 break;
00522 }
00523
00524 case growing_np_problem:
00525 case growing_np_problem_nonlin:{
00526 for (i=0;i<Tt->give_ndofn (idn);i++){
00527 ii=Tt->give_dof(idn,i);
00528 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00529 if (ii==0) r[i]=0.0;
00530 if (ii>0){
00531 r[i]=Lsrst->lhsi[ii-1];
00532 }
00533 }
00534
00535 break;
00536 }
00537 default:{
00538 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00539 }
00540 }
00541 }
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557 void prevnodval (long lcid,double *r, long idn,double *lhs)
00558 {
00559 long i,ii;
00560
00561 switch (Tp->tprob){
00562 case stationary_problem:{
00563 for (i=0;i<Tt->give_ndofn (idn);i++){
00564 ii=Tt->give_dof(idn,i);
00565 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00566 if (ii==0) r[i]=0.0;
00567 if (ii>0) r[i]=lhs[ii-1];
00568 }
00569 break;
00570 }
00571 case nonlinear_stationary_problem:
00572 case nonstationary_problem:{
00573 for (i=0;i<Tt->give_ndofn (idn);i++){
00574 ii=Tt->give_dof(idn,i);
00575 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00576 if (ii==0) r[i]=0.0;
00577 if (ii>0){
00578 r[i]=Lsrst->lhsi[ii-1]+lhs[ii-1];
00579 }
00580 }
00581
00582 break;
00583 }
00584 case nonlinear_nonstationary_problem:{
00585 for (i=0;i<Tt->give_ndofn (idn);i++){
00586 ii=Tt->give_dof(idn,i);
00587 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00588 if (ii==0) r[i]=0.0;
00589 if (ii>0){
00590 r[i]=Lsrst->lhsi[ii-1]+lhs[ii-1];
00591 }
00592 }
00593
00594 break;
00595 }
00596 case discont_nonstat_problem:
00597 case discont_nonlin_nonstat_problem:{
00598 for (i=0;i<Tt->give_ndofn (idn);i++){
00599 ii=Tt->give_dof(idn,i);
00600 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00601 if (ii==0) r[i]=0.0;
00602 if (ii>0){
00603 r[i]=Lsrst->lhsi[ii-1]+lhs[ii-1];
00604 }
00605 }
00606
00607 break;
00608 }
00609
00610 case growing_np_problem:
00611 case growing_np_problem_nonlin:{
00612 for (i=0;i<Tt->give_ndofn (idn);i++){
00613 ii=Tt->give_dof(idn,i);
00614 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00615 if (ii==0) r[i]=0.0;
00616 if (ii>0){
00617 r[i]=Lsrst->lhsi[ii-1]+lhs[ii-1];
00618 }
00619 }
00620
00621
00622
00623 break;
00624 }
00625 default:{
00626 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00627 }
00628 }
00629 }
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642 void nodalotherval (double *r,ivector &nod)
00643 {
00644 long i,j,k,ncompo;
00645
00646 k=0;
00647 for (i=0;i<nod.n;i++){
00648 ncompo = Tt->nodes[nod[i]].ncompother;
00649 for (j=0;j<ncompo;j++){
00650 r[k]=Tt->nodes[nod[i]].other[j];
00651 k++;
00652 }
00653 }
00654
00655 }
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666 void nodalderivatives (double *r,long *cn,long ndofe)
00667 {
00668 long i,ii;
00669
00670 switch (Tp->tprob){
00671 case nonlinear_stationary_problem:
00672 case nonstationary_problem:
00673 case nonlinear_nonstationary_problem:
00674 case growing_np_problem:
00675 case discont_nonstat_problem:
00676 case discont_nonlin_nonstat_problem:
00677 case growing_np_problem_nonlin:{
00678 for (i=0;i<ndofe;i++){
00679 ii=cn[i];
00680 if (ii<0) r[i]=0.0;
00681 if (ii==0) r[i]=0.0;
00682 if (ii>0) r[i]=Lsrst->tdlhs[ii-1];
00683 }
00684 break;
00685 }
00686 default:{
00687 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00688 }
00689 }
00690 }
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701 void prescvalues (double *r,long *cn,long ndofe)
00702 {
00703 long i,ii;
00704
00705 for (i=0;i<ndofe;i++){
00706 ii=cn[i];
00707 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00708 if (ii==0) r[i]=0.0;
00709 if (ii>0){
00710 if(Tp->homogt == 1)
00711 r[i]=0.0;
00712 else
00713 r[i]=Lsrst->lhsi[ii-1];
00714 }
00715 }
00716 }
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727 void initialvalues (double *r,long *cn,long ndofe)
00728 {
00729 long i,ii;
00730
00731 for (i=0;i<ndofe;i++){
00732 ii=cn[i];
00733 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getipv();
00734 if (ii==0) r[i]=0.0;
00735 if (ii>0) r[i]=Lsrst->lhsi[ii-1];
00736 }
00737 }
00738
00739
00740
00741
00742
00743
00744
00745 void approximation_puc ()
00746 {
00747 long i;
00748
00749
00750 for (i=0;i<Tt->ne;i++){
00751 if (Gtt->leso[i]==1){
00752
00753
00754
00755 intpointvalues_puc (i);
00756
00757
00758 }
00759 }
00760
00761 Tm->aux_values (0);
00762 }
00763
00764
00765
00766
00767
00768
00769
00770
00771 void approximation ()
00772 {
00773 long i,ii,j,ipp,nip;
00774
00775 if (Tp->nvs==1){
00776
00777
00778 copy_nodval (0);
00779 }
00780
00781
00782 for (i=0;i<Tt->ne;i++){
00783 if (Gtt->leso[i]==1){
00784
00785
00786
00787 intpointvalues (i);
00788
00789 intpointgradients (i);
00790 }
00791 }
00792
00793
00794 Tm->aux_values (0);
00795
00796
00797 if (Tp->nvs==1){
00798
00799
00800 actual_previous_change ();
00801 }
00802
00803
00804
00805 if (Tm->cemhydr){
00806
00807
00808 for (i=0;i<Tm->nmt;i++){
00809 if (Tm->mattype[i] == cementhydrmat)
00810 for (j=0;j<Tm->numtype[i];j++){
00811 Tm->cemhydr[j].average_temp = 0.0;
00812 Tm->cemhydr[j].n_ipp = 0;
00813 }
00814 }
00815
00816
00817
00818 for (i=0;i<Tt->ne;i++){
00819 if (Gtt->leso[i]==1){
00820 ipp=Tt->elements[i].ipp[0][0];
00821 nip=Tt->give_tnip (i);
00822 for (j=0;j<nip;j++){
00823 ii = Tm->ip[ipp+j].idm;
00824
00825 if(Tp->tmatt == onemedium)
00826 Tm->cemhydr[ii].average_temp += Tm->ip[ipp+j].av[0];
00827 if(Tp->tmatt == twomediacoup)
00828 Tm->cemhydr[ii].average_temp += Tm->ip[ipp+j].av[1];
00829 if(Tp->tmatt == threemediacoup)
00830 Tm->cemhydr[ii].average_temp += Tm->ip[ipp+j].av[2];
00831
00832 Tm->cemhydr[ii].n_ipp++;
00833 }
00834 }
00835 }
00836
00837
00838 for (i=0;i<Tm->nmt;i++){
00839 if (Tm->mattype[i] == cementhydrmat)
00840 for (j=0;j<Tm->numtype[i];j++){
00841 if (Tm->cemhydr[j].n_ipp != 0)
00842 Tm->cemhydr[j].average_temp /= Tm->cemhydr[j].n_ipp;
00843 }
00844 }
00845
00846
00847 for (i=0;i<Tt->ne;i++){
00848 if (Gtt->leso[i]==1){
00849 ipp=Tt->elements[i].ipp[0][0];
00850 nip=Tt->give_tnip (i);
00851 for (j=0;j<nip;j++){
00852 ii = Tm->ip[ipp+j].idm;
00853 Tm->ip[ipp+j].other[0] = Tm->cemhydr[ii].average_temp;
00854 }
00855 }
00856 }
00857 }
00858 }
00859
00860
00861
00862
00863
00864
00865
00866
00867 void actual_previous_change ()
00868 {
00869 long i;
00870
00871 for (i=0;i<Tm->tnip;i++){
00872 Tm->ip[i].actual_previous_change ();
00873 }
00874 }
00875
00876
00877
00878
00879
00880
00881 void actual_previous_nodval ()
00882 {
00883 long i;
00884
00885 for (i=0;i<Tt->nn;i++){
00886 Tt->nodes[i].actual_previous_change ();
00887 }
00888 }
00889
00890
00891
00892
00893
00894
00895
00896
00897 void assemble_init (double *rhs)
00898 {
00899 long i,ndofe,*cn,lcid=0;
00900 vector r,f;
00901 matrix km;
00902 long ne = Tt->ne;
00903
00904 for (i=0;i<ne;i++){
00905
00906 if (Gtt->leso[i]==1){
00907
00908
00909
00910 ndofe=Tt->give_ndofe(i);
00911 allocm (ndofe,ndofe,km);
00912
00913
00914 conductmat (i,lcid,km);
00915
00916 check_math_errel(i);
00917
00918 allocv (ndofe,r);
00919 allocv (ndofe,f);
00920 cn = new long [ndofe];
00921
00922 Tt->give_code_numbers (i,cn);
00923
00924
00925
00926
00927
00928 prescvalues (r.a,cn,ndofe);
00929
00930
00931 mxv (km,r,f);
00932
00933 cmulv (-1.0,f);
00934
00935
00936 locglob (rhs,f.a,cn,ndofe);
00937
00938 destrm (km); destrv (f); destrv (r);
00939 delete [] cn;
00940 }
00941 }
00942 }
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953 void trfel_right_hand_side (long lcid,double *rhs,long n)
00954 {
00955 long i,ndofe;
00956 vector lv;
00957 ivector cn;
00958 double *av;
00959 av = new double [n];
00960
00961
00962 nullv (rhs,n);
00963
00964
00965
00966 assemble_init (rhs);
00967
00968
00969
00970
00971
00972
00973 if(Tp->homogt == 1){
00974 nullv (av,n);
00975 assemble_gradients_contrib(av);
00976 addv (rhs,av,n);
00977 }
00978
00979 for (i=0; i<Tp->ntm; i++)
00980 {
00981 if (Tp->tprob==nonlinear_stationary_problem){
00982
00983
00984
00985 if (Tp->ntm > 1)
00986 {
00987 print_err("support for %ld media have not been implemented\n"
00988 " in nonlinear_stationary_problem", __FILE__, __LINE__, __func__, Tp->ntm);
00989 abort();
00990 }
00991 nullv (av,n);
00992 copyv (rhs,av,n);
00993
00994 nullv (rhs+(lcid*2+0)*n,n);
00995 nullv (rhs+(lcid*2+1)*n,n);
00996
00997 Tb->lc[lcid*2+0].assemble (lcid*2+0,rhs+(lcid*2+0)*n);
00998 addv (rhs,av,n);
00999
01000
01001 }
01002 else{
01003
01004
01005
01006 nullv (av,n);
01007 Tb->lc[i].assemble(i,av);
01008 addv (rhs,av,n);
01009 }
01010 }
01011
01012
01013 nullv (av,n);
01014 for (i=0;i<Tt->ne;i++){
01015 if (Gtt->leso[i]==1){
01016
01017
01018
01019 ndofe=Tt->give_ndofe (i);
01020 reallocv (ndofe,lv);
01021 fillv (0.0,lv);
01022
01023
01024 volume_rhs_vector (lv,lcid,i);
01025
01026 check_math_errel(i);
01027
01028
01029 reallocv(ndofe, cn);
01030 Tt->give_code_numbers (i,cn.a);
01031
01032
01033 locglob (av,lv.a,cn.a,ndofe);
01034 }
01035 }
01036
01037 addv (rhs,av,n);
01038
01039 delete [] av;
01040 }
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051 void trfel_bound_flux (long lcid,double *iflux,long n)
01052 {
01053 long i;
01054 double *av;
01055 av = new double [n];
01056
01057 for(i=0; i<Tp->ntm; i++)
01058 {
01059 Tb->lc[lcid+i].assemble_flux (lcid+i,av,n);
01060 addv (iflux,av,n);
01061 nullv (av,n);
01062 }
01063
01064 delete [] av;
01065 }
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078 void compute_req_valt (long lcid)
01079 {
01080
01081 if (Tp->gradcomp == 1){
01082 approximation();
01083
01084 if (Tp->gradpos==2)
01085 compute_nodegrads ();
01086
01087
01088 }
01089
01090
01091 if (Tp->fluxcomp == 1){
01092
01093
01094
01095 if (Tp->fluxpos==1)
01096 compute_ipfluxes ();
01097 if (Tp->fluxpos==2)
01098 compute_nodefluxes ();
01099
01100
01101 }
01102
01103 if (Tp->tprob == nonlinear_nonstationary_problem ||
01104 Tp->tprob == nonstationary_problem ||
01105 Tp->tprob == discont_nonstat_problem ||
01106 Tp->tprob == discont_nonlin_nonstat_problem ||
01107 Tp->tprob == growing_np_problem ||
01108 Tp->tprob == growing_np_problem_nonlin){
01109
01110
01111 if (Tp->othercomp==1){
01112 if (Tp->otherpos==1)
01113 compute_ipotherst ();
01114 if (Tp->otherpos==2)
01115 compute_nodeotherst ();
01116 if (Tp->otherpos==3)
01117 compute_nodeotherst_comp ();
01118 }
01119
01120
01121 if (Tp->eqothercomp==1){
01122 if (Tp->eqotherpos==1)
01123 compute_ipeqotherst ();
01124 if (Tp->eqotherpos==2)
01125 compute_nodeeqotherst ();
01126
01127
01128 }
01129
01130 }
01131 }
01132
01133
01134
01135 void give_nodal_humid (double *gv,long mnt)
01136 {
01137 long i,cn,compother;
01138 double *r;
01139
01140 r = new double [Tp->ntm];
01141
01142 switch (Tp->tmatt){
01143 case onemedium:{
01144 med1 m1;
01145
01146 for (i=0;i<mnt;i++){
01147 cn=Tt->give_dof (i,0);
01148 if (cn>0) gv[i] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01149 if (cn<0) gv[i] = Tb->lc[0].pv[0-cn-1].getval();
01150
01151
01152
01153
01154
01155 }
01156 break;
01157
01158 gv[i] = m1.compute_othervalues(1,0,r);
01159 if(gv[i] >= 1.0)
01160 gv[i] = 1.0;
01161 }
01162 case twomediacoup:{
01163 med2 m2;
01164
01165 for (i=0;i<mnt;i++){
01166 cn=Tt->give_dof (i,0);
01167 if (cn>0) r[0] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01168 if (cn==0) r[0] = 0.0;
01169 if (cn<0) r[0] = Tb->lc[0].pv[0-cn-1].getval();
01170
01171 cn=Tt->give_dof (i,1);
01172 if (cn>0) r[1] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01173 if (cn==0) r[1] = 0.0;
01174 if (cn<0) r[1] = Tb->lc[0].pv[0-cn-1].getval();
01175
01176 compother = 1;
01177 if ((Tm->ip[0].tm) == kunzel)
01178 compother = 0;
01179
01180 gv[i] = m2.compute_othervalues(compother,0,r);
01181 if(gv[i] >= 1.0)
01182 gv[i] = 1.0;
01183 }
01184 break;
01185 }
01186 case threemediacoup:{
01187 med3 m3;
01188
01189 for (i=0;i<mnt;i++){
01190 cn=Tt->give_dof (i,0);
01191 if (cn>0) r[0] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01192 if (cn==0) r[0] = 0.0;
01193 if (cn<0) r[0] = Tb->lc[0].pv[0-cn-1].getval();
01194
01195 cn=Tt->give_dof (i,1);
01196 if (cn>0) r[1] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01197 if (cn==0) r[1] = 0.0;
01198 if (cn<0) r[1] = Tb->lc[0].pv[0-cn-1].getval();
01199
01200 cn=Tt->give_dof (i,2);
01201 if (cn>0) r[2] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01202 if (cn==0) r[2] = 0.0;
01203 if (cn<0) r[2] = Tb->lc[0].pv[0-cn-1].getval();
01204
01205 gv[i] = m3.compute_othervalues(3,0,r);
01206 if(gv[i] >= 1.0)
01207 gv[i] = 1.0;
01208 }
01209 break;
01210 }
01211 case fourmediacoup:{
01212 med4 m4;
01213
01214 for (i=0;i<mnt;i++){
01215 cn=Tt->give_dof (i,0);
01216 if (cn>0) r[0] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01217 if (cn==0) r[0] = 0.0;
01218 if (cn<0) r[0] = Tb->lc[0].pv[0-cn-1].getval();
01219
01220 cn=Tt->give_dof (i,1);
01221 if (cn>0) r[1] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01222 if (cn==0) r[1] = 0.0;
01223 if (cn<0) r[1] = Tb->lc[0].pv[0-cn-1].getval();
01224
01225 cn=Tt->give_dof (i,2);
01226 if (cn>0) r[2] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01227 if (cn==0) r[2] = 0.0;
01228 if (cn<0) r[2] = Tb->lc[0].pv[0-cn-1].getval();
01229
01230 cn=Tt->give_dof (i,3);
01231 if (cn>0) r[3] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01232 if (cn==0) r[3] = 0.0;
01233 if (cn<0) r[3] = Tb->lc[0].pv[0-cn-1].getval();
01234
01235 gv[i] = m4.compute_othervalues(4,0,r);
01236 if(gv[i] >= 1.0)
01237 gv[i] = 1.0;
01238 }
01239 break;
01240 }
01241 default:{
01242 print_err("unknown number of transported matters is required",__FILE__,__LINE__,__func__);
01243 }
01244 }
01245
01246 delete [] r;
01247 }
01248
01249
01250
01251
01252
01253
01254
01255
01256 void solution_correction ()
01257 {
01258 long i,j,k,cn,ne,nne,ndofn,lcid,ipp;
01259 ivector nodes;
01260 vector nv;
01261
01262
01263 ne = Tt->ne;
01264
01265
01266 lcid=0;
01267
01268
01269 for (i=0;i<ne;i++){
01270
01271 nne = Tt->give_nne (i);
01272
01273
01274 allocv (nne,nodes);
01275
01276
01277 Tt->give_elemnodes (i,nodes);
01278
01279
01280 ipp=Tt->elements[i].ipp[0][0];
01281
01282
01283 for (j=0;j<nne;j++){
01284
01285 ndofn = Tt->give_ndofn (nodes[j]);
01286
01287 allocv (ndofn,nv);
01288
01289 nodalval (lcid,nv.a,nodes[j]);
01290
01291 Tm->values_correction (nv,ipp);
01292
01293
01294 for (k=0;k<ndofn;k++){
01295 cn = Tt->give_dof (nodes[j],k);
01296 if (cn>0){
01297 Lsrst->lhs[cn-1]=nv[k]-Lsrst->lhsi[cn-1];
01298 }
01299 }
01300
01301
01302 destrv (nv);
01303 }
01304
01305 destrv (nodes);
01306 }
01307
01308 }
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321 void compute_cycles (double *plhs,double *pplhs)
01322 {
01323 long i,j,lcid,ndofn,nne,ipp,nc;
01324 double *r,*pr,*ppr;
01325 ivector nod;
01326
01327
01328 lcid=0;
01329
01330 if (Tt->ncycl == NULL){
01331 Tt->ncycl = new long [Tt->nn];
01332 Tt->aux = new long[Tt->nn];
01333
01334 for (i=0;i<Tt->nn;i++){
01335 Tt->ncycl[i]=0;
01336 }
01337 }
01338
01339 for (i=0;i<Tt->nn;i++){
01340 Tt->aux[i]=0;
01341 }
01342
01343 for (i=0;i<Tt->ne;i++){
01344
01345 nne = Tt->give_nne (i);
01346 allocv (nne,nod);
01347
01348 Tt->give_elemnodes (i,nod);
01349
01350 ipp = Tt->elements[i].ipp[0][0];
01351
01352 for (j=0;j<nne;j++){
01353 if (Tt->aux[nod[j]]>0)
01354 continue;
01355
01356
01357 ndofn = Tt->give_ndofn (nod[j]);
01358 r = new double [ndofn];
01359 pr = new double [ndofn];
01360 ppr = new double [ndofn];
01361
01362 nodalval (lcid,r,nod[j]);
01363 prevnodval (lcid,pr,nod[j],plhs);
01364 prevnodval (lcid,ppr,nod[j],pplhs);
01365
01366 nc=Tm->cycle_detection (r,pr,ppr,ipp);
01367
01368 Tt->ncycl[nod[j]]+=nc;
01369 Tt->aux[nod[j]]=1;
01370
01371 delete [] r; delete [] pr; delete [] ppr;
01372 }
01373
01374 destrv (nod);
01375 }
01376 }
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387 void copy_nodval (long lcid)
01388 {
01389 long i,ndofn;
01390 double *r;
01391
01392
01393 for (i=0;i<Tt->nn;i++){
01394
01395 ndofn = Tt->give_ndofn (i);
01396
01397 r = new double [ndofn];
01398
01399 nodalval (lcid,r,i);
01400
01401 Tt->nodes[i].save_nodval (r);
01402
01403 delete [] r;
01404 }
01405 }
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416 void assemble_gradients_contrib (double *rhs)
01417 {
01418 long i,ii,jj,k,ndofe,*cn,lcid=0,ntm,ncomp;
01419 vector grad,f;
01420 matrix lm;
01421 long ne = Tt->ne;
01422
01423
01424 for (i=0;i<ne;i++){
01425
01426 if (Gtt->leso[i]==1){
01427
01428
01429
01430 ndofe=Tt->give_ndofe(i);
01431 ntm = Tp->ntm;
01432 ncomp = Tt->give_ncomp(i);
01433 allocm (ndofe,ntm*ncomp,lm);
01434
01435
01436 ltmat (i,lcid,lm);
01437
01438 allocv (ntm*ncomp,grad);
01439 allocv (ndofe,f);
01440 cn = new long [ndofe];
01441
01442 Tt->give_code_numbers (i,cn);
01443
01444
01445 k = 0;
01446 for(ii=0;ii<ntm;ii++){
01447 for(jj=0;jj<ncomp;jj++){
01448 grad[k]=Tb->lc[ii].mastergrad[jj];
01449 k++;
01450 }
01451 }
01452
01453
01454 mxv (lm,grad,f);
01455
01456 cmulv (-1.0,f);
01457
01458
01459 locglob (rhs,f.a,cn,ndofe);
01460
01461 destrm (lm); destrv (f); destrv (grad);
01462 delete [] cn;
01463 }
01464 }
01465 }
01466
01467
01468
01469
01470
01471
01472
01473
01474 void assemble_l_matrix (matrix &lm, matrix <m)
01475 {
01476 long i,ndofe,*rcn,*ccn,lcid=0,ntm,ncomp;
01477 matrix lkm,lkmt;
01478
01479
01480 for (i=0;i<Tt->ne;i++){
01481 if (Gtt->leso[i]==1){
01482
01483
01484
01485 ndofe=Tt->give_ndofe(i);
01486 ntm = Tp->ntm;
01487 ncomp = Tt->give_ncomp(i);
01488 allocm (ndofe,ntm*ncomp,lkmt);
01489 allocm (ntm*ncomp,ndofe,lkm);
01490
01491
01492 ltmat (i,lcid,lkmt);
01493
01494
01495 lmat (i,lcid,lkm);
01496
01497 rcn = new long [ndofe];
01498 ccn = new long [ntm*ncomp];
01499
01500 Gtt->give_code_numbers (i,rcn);
01501 if((ntm*ncomp) == 1){
01502 ccn[0] = 1;
01503 }
01504 if((ntm*ncomp) == 2){
01505 ccn[0] = 1;
01506 ccn[1] = 2;
01507 }
01508 if((ntm*ncomp) == 3){
01509 ccn[0] = 1;
01510 ccn[1] = 2;
01511 ccn[2] = 3;
01512 }
01513 if((ntm*ncomp) == 4){
01514 ccn[0] = 1;
01515 ccn[1] = 2;
01516 ccn[2] = 3;
01517 ccn[3] = 4;
01518 }
01519 if((ntm*ncomp) == 6){
01520 ccn[0] = 1;
01521 ccn[1] = 2;
01522 ccn[2] = 3;
01523 ccn[3] = 4;
01524 ccn[4] = 5;
01525 ccn[5] = 6;
01526 }
01527 if((ntm*ncomp) == 9){
01528 ccn[0] = 1;
01529 ccn[1] = 2;
01530 ccn[2] = 3;
01531 ccn[3] = 4;
01532 ccn[4] = 5;
01533 ccn[5] = 6;
01534 ccn[6] = 7;
01535 ccn[7] = 8;
01536 ccn[8] = 9;
01537 }
01538
01539
01540 mat_localize (ltm,lkmt,rcn,ccn);
01541
01542
01543 mat_localize (lm,lkm,ccn,rcn);
01544
01545 destrm (lkm);
01546 destrm (lkmt);
01547 delete [] rcn;
01548 delete [] ccn;
01549 }
01550 }
01551 }
01552
01553
01554
01555
01556
01557
01558
01559 void assemble_conductivity_matrix (matrix &km)
01560 {
01561 long i,ndofe,*cn;
01562 matrix lm;
01563
01564
01565 for (i=0;i<Tt->ne;i++){
01566 if (Gtt->leso[i]==1){
01567
01568 ndofe = Gtt->give_ndofe (i);
01569
01570 allocm (ndofe,ndofe,lm);
01571
01572 conductmat (i,0,lm);
01573
01574 check_math_errel(i);
01575
01576 cn = new long [ndofe];
01577 Gtt->give_code_numbers (i,cn);
01578
01579 mat_localize (km,lm,cn,cn);
01580
01581 destrm (lm);
01582 delete [] cn;
01583 }
01584 }
01585 }
01586
01587
01588
01589
01590
01591
01592 void assemble_average_d_matrix (matrix &km, double &area)
01593 {
01594 long i,ii,jj,ntm,ncomp;
01595 double elemarea,invarea;
01596 matrix lm;
01597
01598 area = 0.0;
01599
01600 for (i=0;i<Tt->ne;i++){
01601 elemarea = 0.0;
01602 if (Gtt->leso[i]==1){
01603
01604 ntm = Tp->ntm;
01605 ncomp = Tt->give_ncomp(i);
01606
01607 allocm (ncomp*ntm,ncomp*ntm,lm);
01608
01609
01610 averdmat (i,elemarea,lm);
01611
01612 check_math_errel(i);
01613
01614 for (ii=0;ii<ntm*ncomp;ii++){
01615 for (jj=0;jj<ntm*ncomp;jj++){
01616 km[ii][jj]+=lm[ii][jj]*elemarea;
01617 }
01618 }
01619 destrm (lm);
01620
01621 area += elemarea;
01622 }
01623 }
01624
01625 invarea = 1.0/area;
01626 cmulm(invarea,km);
01627
01628 }
01629
01630
01631
01632
01633
01634
01635 void assemble_average_c_matrix (matrix &cm)
01636 {
01637 long i,ii,jj,ntm,ncomp;
01638 double elemarea,area;
01639 matrix lm;
01640
01641 area = 0.0;
01642
01643 for (i=0;i<Tt->ne;i++){
01644 elemarea = 0.0;
01645 if (Gtt->leso[i]==1){
01646
01647 ntm = Tp->ntm;
01648 ncomp = Tt->give_ncomp(i);
01649
01650 allocm (ntm,ntm,lm);
01651
01652
01653 avercmat (i,elemarea,lm);
01654
01655 check_math_errel(i);
01656
01657 for (ii=0;ii<ntm;ii++){
01658 for (jj=0;jj<ntm;jj++){
01659 cm[ii][jj]+=lm[ii][jj]*elemarea/Tb->lc[jj].masterval;
01660
01661 }
01662 }
01663 destrm (lm);
01664
01665 area += elemarea;
01666 }
01667 }
01668
01669 area = 1.0/area;
01670 cmulm(area,cm);
01671
01672 }
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682 double total_integral (long lcid)
01683 {
01684 long i,nne;
01685 double totint;
01686 ivector cn;
01687 vector nodval;
01688
01689 totint=0.0;
01690
01691
01692 for (i=0;i<Tt->ne;i++){
01693
01694
01695 nne = Tt->give_nne (i);
01696 reallocv (nne,nodval);
01697 reallocv (nne,cn);
01698
01699
01700 Tt->give_medium_code_numbers (i,lcid,cn.a);
01701
01702 elemvalues (lcid,i,nodval.a,cn.a,nne);
01703
01704 totint+=elem_total_integral (i,nodval);
01705
01706 check_math_errel(i);
01707 }
01708
01709 return totint;
01710 }
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720 double total_integral_ip (long varid)
01721 {
01722 long i;
01723 double totint;
01724
01725
01726 totint=0.0;
01727
01728
01729 for (i=0;i<Tt->ne;i++){
01730 totint+=elem_total_integral_ip (i,varid);
01731 }
01732
01733 return totint;
01734 }
01735
01736
01737
01738
01739
01740
01741 void dt_subdom ()
01742 {
01743 long i,j;
01744 long ns,nes,eid;
01745
01746
01747 ns = Gtt->stop->ns;
01748
01749
01750 for (i=0;i<ns;i++){
01751
01752 nes = Gtt->stop->ned[i];
01753
01754
01755 for (j=0;j<nes;j++){
01756
01757 eid = Gtt->stop->domel[i][j];
01758
01759
01760
01761
01762
01763 }
01764 }
01765 }
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777 void lnso_leso_setting (long *lsso)
01778 {
01779 long i,j;
01780 long ns,nes,nne,eid;
01781 ivector nodes;
01782
01783
01784 ns = Gtt->stop->ns;
01785
01786
01787
01788
01789
01790 for (i=0;i<ns;i++){
01791
01792 nes = Gtt->stop->ned[i];
01793
01794 if (lsso[i]==0){
01795
01796 for (j=0;j<nes;j++){
01797
01798 eid = Gtt->stop->domel[i][j];
01799
01800 Gtt->leso[eid]=0;
01801 }
01802 }
01803 if (lsso[i]==1){
01804
01805 for (j=0;j<nes;j++){
01806
01807 eid = Gtt->stop->domel[i][j];
01808
01809 Gtt->leso[eid]=1;
01810 }
01811 }
01812 }
01813
01814
01815
01816
01817
01818 for (i=0;i<Tt->nn;i++){
01819 Gtt->lnso[i]=1;
01820 }
01821
01822
01823 for (i=0;i<Tt->ne;i++){
01824
01825 if (Gtt->leso[i]==0){
01826
01827
01828 nne = Tt->give_nne (i);
01829 reallocv (nne,nodes);
01830
01831 Tt->give_elemnodes (i,nodes);
01832
01833
01834 for (j=0;j<nne;j++){
01835 Gtt->lnso[nodes[j]]=0;
01836 }
01837 }
01838 }
01839
01840 }