00001
00002
00003
00004
00005
00006
00007 #include "globalt.h"
00008 #include "quadquadrilatt.h"
00009 #include "genfile.h"
00010 #include "globmatt.h"
00011
00012 quadquadrilatt::quadquadrilatt (void)
00013 {
00014 long i;
00015
00016
00017 nne=8;
00018
00019 ned=4;
00020
00021 nned=3;
00022
00023 ncomp=2;
00024
00025 ntm=Tp->ntm;
00026
00027
00028 dofe = new long* [ntm];
00029 nip = new long* [ntm];
00030 intordkm = new long* [ntm];
00031 intordcm = new long* [ntm];
00032 ordering = new long* [ntm];
00033 for (i=0;i<ntm;i++){
00034 dofe[i] = new long [ntm];
00035 nip[i] = new long [ntm];
00036 intordkm[i] = new long [ntm];
00037 intordcm[i] = new long [ntm];
00038 ordering[i] = new long [nne];
00039 }
00040
00041
00042 switch (Tp->tmatt){
00043 case nomedium:{ break; }
00044 case onemedium:{
00045 ordering[0][0]=1; ordering[0][1]=2; ordering[0][2]=3; ordering[0][3]=4;
00046 ordering[0][4]=5; ordering[0][5]=6; ordering[0][6]=7; ordering[0][7]=8;
00047 dofe[0][0]=8; intordkm[0][0]=3; intordcm[0][0]=3;
00048 nip[0][0]=18;
00049 ndofe=8; napfun=1;
00050 break;
00051 }
00052 case twomediacoup:{
00053 ordering[0][0]=1; ordering[0][1]=3; ordering[0][2]=5; ordering[0][3]=7;
00054 ordering[0][4]=9; ordering[0][5]=11; ordering[0][6]=13; ordering[0][7]=15;
00055
00056 ordering[1][0]=2; ordering[1][1]=4; ordering[1][2]=6; ordering[1][3]=8;
00057 ordering[1][4]=10; ordering[1][5]=12; ordering[1][6]=14; ordering[1][7]=16;
00058
00059 intordkm[0][0]=3; intordkm[0][1]=3; intordkm[1][0]=3; intordkm[1][1]=3;
00060 intordcm[0][0]=3; intordcm[0][1]=3; intordcm[1][0]=3; intordcm[1][1]=3;
00061
00062 if (Tp->savemode==0){
00063 nip[0][0]=18; nip[0][1]=18; nip[1][0]=18; nip[1][1]=18;
00064 }
00065 if (Tp->savemode==1){
00066 nip[0][0]=18; nip[0][1]=0; nip[1][0]=0; nip[1][1]=0;
00067 }
00068
00069 dofe[0][0]=8; dofe[0][1]=8; dofe[1][0]=8; dofe[1][1]=8;
00070
00071 ndofe=16; napfun=2;
00072 break;
00073 }
00074 case threemediacoup:{
00075 ordering[0][0]=1; ordering[0][1] =4; ordering[0][2] =7; ordering[0][3] =10;
00076 ordering[0][4]=13; ordering[0][5] =16; ordering[0][6] =19; ordering[0][7] =22;
00077
00078 ordering[1][0]=2; ordering[1][1] =5; ordering[1][2] =8; ordering[1][3] =11;
00079 ordering[1][4]=14; ordering[1][5] =17; ordering[1][6] =20; ordering[1][7] =23;
00080
00081 ordering[2][0]=3; ordering[2][1] =6; ordering[2][2] =9; ordering[2][3] =12;
00082 ordering[2][4]=15; ordering[2][5] =18; ordering[2][6] =21; ordering[2][7] =24;
00083
00084 intordkm[0][0]=3; intordkm[0][1]=3; intordkm[0][2]=3;
00085 intordkm[1][0]=3; intordkm[1][1]=3; intordkm[1][2]=3;
00086 intordkm[2][0]=3; intordkm[2][1]=3; intordkm[2][2]=3;
00087
00088 intordcm[0][0]=3; intordcm[0][1]=3; intordcm[0][2]=3;
00089 intordcm[1][0]=3; intordcm[1][1]=3; intordcm[1][2]=3;
00090 intordcm[2][0]=3; intordcm[2][1]=3; intordcm[2][2]=3;
00091
00092 if (Tp->savemode==0){
00093 nip[0][0]=18; nip[0][1]=18; nip[0][2]=18;
00094 nip[1][0]=18; nip[1][1]=18; nip[1][2]=18;
00095 nip[2][0]=18; nip[2][1]=18; nip[2][2]=18;
00096 }
00097 if (Tp->savemode==1){
00098 nip[0][0]=18; nip[0][1]=0; nip[0][2]=0;
00099 nip[1][0]=0; nip[1][1]=0; nip[1][2]=0;
00100 nip[2][0]=0; nip[2][1]=0; nip[2][2]=0;
00101 }
00102
00103 dofe[0][0]=8; dofe[0][1]=8; dofe[0][2]=8;
00104 dofe[1][0]=8; dofe[1][1]=8; dofe[1][2]=8;
00105 dofe[2][0]=8; dofe[2][1]=8; dofe[2][2]=8;
00106
00107 ndofe=24; napfun=3;
00108 break;
00109 }
00110 default:{
00111 print_err("unknown type of transported matter is required",__FILE__,__LINE__,__func__);
00112 }
00113 }
00114 }
00115
00116 quadquadrilatt::~quadquadrilatt (void)
00117 {
00118 long i;
00119
00120 for (i=0;i<ntm;i++){
00121 delete [] dofe[i];
00122 delete [] nip[i];
00123 delete [] intordkm[i];
00124 delete [] intordcm[i];
00125 delete [] ordering[i];
00126 }
00127 delete [] dofe;
00128 delete [] nip;
00129 delete [] intordkm;
00130 delete [] intordcm;
00131 delete [] ordering;
00132 }
00133
00134 void quadquadrilatt::codnum (long *cn,long ri)
00135 {
00136 long i;
00137 for (i=0;i<nne;i++){
00138 cn[i]=ordering[ri][i];
00139 }
00140 }
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 double quadquadrilatt::approx (double xi,double eta,vector &nodval)
00151 {
00152 double f;
00153 vector bf(nne);
00154
00155 bf_quad_4_2d (bf.a,xi,eta);
00156
00157 scprd (bf,nodval,f);
00158
00159 return f;
00160 }
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 void quadquadrilatt::ipcoordblock (long eid,long ri,long ci,double **coord)
00173 {
00174 long i,j,k;
00175 double xi,eta;
00176 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]);
00177
00178 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00179 Tt->give_node_coord2d (x,y,eid);
00180
00181 k=0;
00182 for (i=0;i<intordkm[ri][ci];i++){
00183 xi=gp[i];
00184 for (j=0;j<intordkm[ri][ci];j++){
00185 eta=gp[j];
00186
00187 coord[k][0]=approx (xi,eta,x);
00188 coord[k][1]=approx (xi,eta,y);
00189 coord[k++][2]=0.0;
00190 }
00191 }
00192 }
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 void quadquadrilatt::intpointval (long eid)
00203 {
00204 long i,j,k,ii,jj,ipp;
00205 double xi,eta,val;
00206 ivector cn(ndofe);
00207 vector r(ndofe),t(nne),gp,w;
00208
00209 Tt->give_code_numbers (eid,cn.a);
00210 elemvalues (0,eid,r.a,cn.a,ndofe);
00211
00212 for (k=0;k<Tp->ntm;k++){
00213
00214
00215 for (i=0;i<dofe[k][k];i++){
00216 t[i]=r[ordering[k][i]-1];
00217 }
00218
00219 for (ii=0;ii<Tp->ntm;ii++){
00220 for (jj=0;jj<Tp->ntm;jj++){
00221
00222
00223 reallocv (intordkm[ii][jj],gp);
00224 reallocv (intordkm[ii][jj],w);
00225 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00226
00227 ipp=Tt->elements[eid].ipp[ii][jj];
00228 for (i=0;i<intordkm[ii][jj];i++){
00229 xi=gp[i];
00230 for (j=0;j<intordkm[ii][jj];j++){
00231 eta=gp[j];
00232
00233 val = approx (xi,eta,t);
00234 Tm->ip[ipp].av[k]=val;
00235 ipp++;
00236 }
00237 }
00238
00239
00240 reallocv (intordcm[ii][jj],gp);
00241 reallocv (intordcm[ii][jj],w);
00242 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00243
00244 for (i=0;i<intordcm[ii][jj];i++){
00245 xi=gp[i];
00246 for (j=0;j<intordcm[ii][jj];j++){
00247 eta=gp[j];
00248
00249 val = approx (xi,eta,t);
00250 Tm->ip[ipp].av[k]=val;
00251 ipp++;
00252 }
00253 }
00254
00255 if (Tp->savemode==1) break;
00256 }
00257 if (Tp->savemode==1) break;
00258 }
00259 }
00260 }
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270 void quadquadrilatt::intpointgrad (long eid)
00271 {
00272 long i,j,ii,jj,k,ipp;
00273 double xi,eta,jac;
00274 ivector cn(ndofe);
00275 vector x(nne),y(nne),r(ndofe),t(nne),gp,w,grad(ncomp);
00276 matrix gm(ncomp,nne);
00277
00278 Tt->give_node_coord2d (x,y,eid);
00279 Tt->give_code_numbers (eid,cn.a);
00280 elemvalues (0,eid,r.a,cn.a,ndofe);
00281
00282 for (k=0;k<Tp->ntm;k++){
00283
00284
00285 for (i=0;i<dofe[k][k];i++){
00286 t[i]=r[ordering[k][i]-1];
00287 }
00288
00289 for (ii=0;ii<Tp->ntm;ii++){
00290 for (jj=0;jj<Tp->ntm;jj++){
00291
00292
00293 reallocv (intordkm[ii][jj],gp);
00294 reallocv (intordkm[ii][jj],w);
00295 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00296
00297 ipp=Tt->elements[eid].ipp[ii][jj];
00298 for (i=0;i<intordkm[ii][jj];i++){
00299 xi=gp[i];
00300 for (j=0;j<intordkm[ii][jj];j++){
00301 eta=gp[j];
00302
00303
00304 grad_matrix (gm,x,y,xi,eta,jac);
00305 mxv (gm,t,grad);
00306 Tm->storegrad (k,ipp,grad);
00307 ipp++;
00308 }
00309 }
00310
00311
00312
00313 reallocv (intordcm[ii][jj],gp);
00314 reallocv (intordcm[ii][jj],w);
00315 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00316
00317 for (i=0;i<intordcm[ii][jj];i++){
00318 xi=gp[i];
00319 for (j=0;j<intordcm[ii][jj];j++){
00320 eta=gp[j];
00321
00322
00323 grad_matrix (gm,x,y,xi,eta,jac);
00324 mxv (gm,t,grad);
00325 Tm->storegrad (k,ipp,grad);
00326 ipp++;
00327 }
00328 }
00329
00330 if (Tp->savemode==1) break;
00331 }
00332 if (Tp->savemode==1) break;
00333 }
00334 }
00335 }
00336
00337
00338
00339
00340
00341
00342
00343
00344 void quadquadrilatt::intpointother (long eid)
00345 {
00346 long i,j,k,ii,jj,ipp,ncompo,nodid;
00347 double xi,eta,val,*r;
00348 ivector nodes(nne);
00349 vector t(nne),gp,w;
00350
00351
00352 Tt->give_elemnodes (eid,nodes);
00353
00354
00355 nodid=nodes[0];
00356
00357
00358 ncompo=Tt->nodes[nodid].ncompother;
00359 r = new double [ncompo*nne];
00360
00361
00362 nodalotherval (r,nodes);
00363
00364 for (k=0;k<ncompo;k++){
00365
00366
00367 for (i=0;i<nne;i++){
00368 t[i]=r[i*ncompo+k];
00369 }
00370
00371 for (ii=0;ii<Tp->ntm;ii++){
00372 for (jj=0;jj<Tp->ntm;jj++){
00373
00374
00375 reallocv (intordkm[ii][jj],gp);
00376 reallocv (intordkm[ii][jj],w);
00377 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00378
00379 ipp=Tt->elements[eid].ipp[ii][jj];
00380 for (i=0;i<intordkm[ii][jj];i++){
00381 xi=gp[i];
00382 for (j=0;j<intordkm[ii][jj];j++){
00383 eta=gp[j];
00384
00385 val = approx (xi,eta,t);
00386 Tm->ip[ipp].other[k]=val;
00387 ipp++;
00388 }
00389 }
00390
00391
00392
00393 reallocv (intordcm[ii][jj],gp);
00394 reallocv (intordcm[ii][jj],w);
00395 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00396
00397 for (i=0;i<intordcm[ii][jj];i++){
00398 xi=gp[i];
00399 for (j=0;j<intordcm[ii][jj];j++){
00400 eta=gp[j];
00401
00402 val = approx (xi,eta,t);
00403 Tm->ip[ipp].other[k]=val;
00404 ipp++;
00405 }
00406 }
00407
00408 if (Tp->savemode==1) break;
00409 }
00410 if (Tp->savemode==1) break;
00411 }
00412 }
00413
00414 delete [] r;
00415 }
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425 void quadquadrilatt::bf_matrix (matrix &n,double xi,double eta)
00426 {
00427 fillm (0.0,n);
00428 bf_quad_4_2d (n.a,xi,eta);
00429 }
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440 void quadquadrilatt::give_approx_fun (double &f,double xi,double eta,long i)
00441 {
00442 double n[4];
00443 bf_quad_4_2d (n,xi,eta);
00444 f=n[i];
00445 }
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458 void quadquadrilatt::grad_matrix (matrix &gm,vector &x,vector &y,double xi,double eta,double &jac)
00459 {
00460 long i;
00461 vector dx(nne),dy(nne);
00462
00463 dx_bf_quad_4_2d (dx.a,xi,eta);
00464 dy_bf_quad_4_2d (dy.a,xi,eta);
00465
00466 derivatives_2d (dx,dy,jac,x,y,xi,eta);
00467
00468 fillm (0.0,gm);
00469
00470 for (i=0;i<nne;i++){
00471 gm[0][i]=dx[i];
00472 gm[1][i]=dy[i];
00473 }
00474
00475 }
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490 void quadquadrilatt::conductivity_matrix (long lcid,long eid,long ri,long ci,matrix &km)
00491 {
00492 long i,j,ii;
00493 double thick,xi,eta,ww1,ww2,jac;
00494 ivector nodes(nne);
00495 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
00496 matrix gm(ncomp,dofe[ri][ci]),d;
00497
00498 matrix n(1,dofe[ri][ci]);
00499
00500 Tt->give_elemnodes (eid,nodes);
00501 Tc->give_thickness (eid,nodes,t);
00502
00503 Tt->give_node_coord2d (x,y,eid);
00504 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00505
00506 fillm (0.0,km);
00507
00508 if (Tp->savemode==0)
00509 ii=Tt->elements[eid].ipp[ri][ci];
00510 if (Tp->savemode==1)
00511 ii=Tt->elements[eid].ipp[0][0];
00512
00513 for (i=0;i<intordkm[ri][ci];i++){
00514 xi=gp[i]; ww1=w[i];
00515 for (j=0;j<intordkm[ri][ci];j++){
00516 eta=gp[j]; ww2=w[j];
00517
00518
00519 grad_matrix (gm,x,y,xi,eta,jac);
00520
00521
00522 reallocm(ncomp,ncomp,d);
00523 Tm->matcond (d,ii,ri,ci);
00524
00525 thick = approx (xi,eta,t);
00526
00527 if (jac<0.0){
00528 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00529 jac=fabs(jac);
00530 }
00531
00532 jac*=thick*ww1*ww2;
00533
00534
00535 bdbj (km.a,gm.a,d.a,jac,gm.m,gm.n);
00536
00537
00538 reallocm(1,ncomp,d);
00539 Tm->matcond2(d,ii,ri,ci);
00540 bf_matrix (n, xi, eta);
00541 bdbjac(km, n, d, gm, jac);
00542
00543 ii++;
00544 }
00545 }
00546
00547 if ((Tt->elements[eid].transi[lcid]==3)||(Tt->elements[eid].transi[lcid]==4)){
00548
00549 transmission_matrix (lcid,eid,ri,ci,km);
00550 }
00551 }
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565 void quadquadrilatt::capacity_matrix (long eid,long ri,long ci,matrix &cm)
00566 {
00567 long i,j,ii;
00568 double jac,xi,eta,w1,w2,thick,rho,c;
00569 ivector nodes(nne);
00570 vector x(nne),y(nne),w(intordcm[ri][ci]),gp(intordcm[ri][ci]),t(nne),dens(nne);
00571 matrix n(1,dofe[ri][ci]);
00572
00573 Tt->give_elemnodes (eid,nodes);
00574 Tc->give_thickness (eid,nodes,t);
00575 Tc->give_density (eid,nodes,dens);
00576
00577 Tt->give_node_coord2d (x,y,eid);
00578 gauss_points (gp.a,w.a,intordcm[ri][ci]);
00579
00580 fillm (0.0,cm);
00581
00582 if (Tp->savemode==0)
00583 ii=Tt->elements[eid].ipp[ri][ci]+intordkm[ri][ci]*intordkm[ri][ci];
00584 if (Tp->savemode==1)
00585 ii=Tt->elements[eid].ipp[0][0]+intordkm[0][0]*intordkm[0][0];
00586
00587 for (i=0;i<intordcm[ri][ci];i++){
00588 xi=gp[i]; w1=w[i];
00589 for (j=0;j<intordcm[ri][ci];j++){
00590 eta=gp[j]; w2=w[j];
00591 jac_2d (jac,x,y,xi,eta);
00592 bf_matrix (n,xi,eta);
00593
00594 c=Tm->capcoeff (ii,ri,ci);
00595 thick = approx (xi,eta,t);
00596 rho = approx (xi,eta,dens);
00597
00598 if (jac<0.0){
00599 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00600 jac=fabs(jac);
00601 }
00602
00603 jac*=w1*w2*thick*rho*c;
00604
00605 nnj (cm.a,n.a,jac,n.m,n.n);
00606 ii++;
00607 }
00608 }
00609 }
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624 void quadquadrilatt::quantity_source_vector (vector &sv,vector &nodval,long eid,long ri,long ci)
00625 {
00626 long i,j;
00627 double jac,xi,eta,w1,w2,thick;
00628 ivector nodes(nne);
00629 vector x(nne),y(nne),w(intordcm[ri][ci]),gp(intordcm[ri][ci]),t(nne),dens(nne),v(dofe[ri][ci]);
00630 matrix n(1,dofe[ri][ci]),nm(dofe[ri][ci],dofe[ri][ci]);
00631
00632 Tt->give_elemnodes (eid,nodes);
00633
00634 Tc->give_thickness (eid,nodes,t);
00635
00636 Tt->give_node_coord2d (x,y,eid);
00637 gauss_points (gp.a,w.a,intordcm[ri][ci]);
00638
00639 fillm (0.0,nm);
00640
00641 for (i=0;i<intordcm[ri][ci];i++){
00642 xi=gp[i]; w1=w[i];
00643 for (j=0;j<intordcm[ri][ci];j++){
00644 eta=gp[j]; w2=w[j];
00645 jac_2d (jac,x,y,xi,eta);
00646 bf_matrix (n,xi,eta);
00647
00648 thick = approx (xi,eta,t);
00649
00650 if (jac<0.0){
00651 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00652 jac=fabs(jac);
00653 }
00654
00655 jac*=w1*w2*thick;
00656
00657 nnj (nm.a,n.a,jac,n.m,n.n);
00658 }
00659 }
00660 mxv (nm,nodval,v); addv (sv,v,sv);
00661
00662 }
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675 void quadquadrilatt::internal_fluxes (long lcid,long eid,vector &ifl)
00676 {
00677 long i,j,ipp;
00678 double xi,eta,jac,thick;
00679 ivector nodes(nne);
00680 vector x(nne),y(nne),w,gp,t(nne),fl(ncomp),contr(dofe[lcid][lcid]);
00681 matrix gm(ncomp,dofe[lcid][lcid]);
00682
00683
00684 Tt->give_elemnodes (eid,nodes);
00685 Tc->give_thickness (eid,nodes,t);
00686 Tt->give_node_coord2d (x,y,eid);
00687
00688 reallocv (intordkm[lcid][lcid],gp);
00689 reallocv (intordkm[lcid][lcid],w);
00690
00691 gauss_points (gp.a,w.a,intordkm[lcid][lcid]);
00692 if (Tp->savemode==0)
00693 ipp=Tt->elements[eid].ipp[lcid][lcid];
00694 if (Tp->savemode==1)
00695 ipp=Tt->elements[eid].ipp[0][0];
00696
00697 for (i=0;i<intordkm[lcid][lcid];i++){
00698 xi=gp[i];
00699 for (j=0;j<intordkm[lcid][lcid];j++){
00700 eta=gp[j];
00701 thick = approx (xi,eta,t);
00702
00703 Tm->computenlfluxes (lcid,ipp);
00704
00705 Tm->givefluxes (lcid,ipp,fl);
00706
00707 grad_matrix (gm,x,y,xi,eta,jac);
00708 if (jac<0.0){
00709 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00710 jac=fabs(jac);
00711 }
00712
00713 mtxv (gm,fl,contr);
00714
00715 cmulv (jac*w[i]*w[j]*thick,contr);
00716
00717 addv (contr,ifl,ifl);
00718
00719 ipp++;
00720 }
00721 }
00722
00723 }
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734 void quadquadrilatt::res_conductivity_matrix (long eid,long lcid,matrix &km)
00735 {
00736 long i,j,*rcn,*ccn;
00737 matrix lkm;
00738
00739 for (i=0;i<ntm;i++){
00740 for (j=0;j<ntm;j++){
00741 rcn = new long [dofe[i][j]];
00742 ccn = new long [dofe[i][j]];
00743 reallocm (dofe[i][j],dofe[i][j],lkm);
00744 conductivity_matrix (i,eid,i,j,lkm);
00745 codnum (rcn,i);
00746 codnum (ccn,j);
00747 mat_localize (km,lkm,rcn,ccn);
00748 delete [] rcn; delete [] ccn;
00749 }
00750 }
00751 }
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761 void quadquadrilatt::res_capacity_matrix (long eid,matrix &cm)
00762 {
00763 long i,j,*rcn,*ccn,k,l;
00764 double s;
00765 matrix lcm;
00766
00767 for (i=0;i<ntm;i++){
00768 for (j=0;j<ntm;j++){
00769 rcn = new long [dofe[i][j]];
00770 ccn = new long [dofe[i][j]];
00771 reallocm (dofe[i][j],dofe[i][j],lcm);
00772 capacity_matrix (eid,i,j,lcm);
00773
00774
00775 if (Tp->diagcap == 1){
00776 for (k=0;k<dofe[i][j];k++){
00777 s=0.0;
00778 for (l=0;l<dofe[i][j];l++){
00779 s+=lcm[k][l];
00780 lcm[k][l]=0.0;
00781 }
00782 lcm[k][k]=s;
00783 }
00784 }
00785
00786 codnum (rcn,i);
00787 codnum (ccn,j);
00788 mat_localize (cm,lcm,rcn,ccn);
00789 delete [] rcn; delete [] ccn;
00790 }
00791 }
00792
00793 }
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805 void quadquadrilatt::res_convection_vector (vector &f,long lcid,long eid,long leid)
00806 {
00807 long *cn;
00808 vector lf;
00809
00810 cn = new long [dofe[lcid][lcid]];
00811 reallocv (dofe[lcid][lcid],lf);
00812
00813 convection_vector (lf,lcid,eid,leid,lcid,lcid);
00814
00815 codnum (cn,lcid);
00816 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
00817 delete [] cn;
00818
00819 cmulv (-1.0,f,f);
00820
00821 }
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833 void quadquadrilatt::res_transmission_vector (vector &f,long lcid,long eid,long leid)
00834 {
00835 long *cn,i;
00836 vector lf;
00837
00838 cn = new long [dofe[lcid][lcid]];
00839 codnum (cn,lcid);
00840
00841
00842 for (i=0;i<ntm;i++){
00843 reallocv (dofe[lcid][lcid],lf);
00844
00845 if ((Tt->elements[eid].transi[i]==3)||(Tt->elements[eid].transi[i]==4)){
00846 transmission_vector (lf,i,eid,leid,lcid,i);
00847 }
00848
00849 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
00850 }
00851
00852 delete [] cn;
00853 }
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864 void quadquadrilatt::res_quantity_source_vector (vector &sv,vector &nodval,long lcid,long eid)
00865 {
00866 long *cn;
00867 vector lsv;
00868
00869 cn = new long [dofe[lcid][lcid]];
00870 reallocv (dofe[lcid][lcid],lsv);
00871 quantity_source_vector (lsv,nodval,eid,lcid,lcid);
00872 codnum (cn,lcid);
00873 locglob (sv.a,lsv.a,cn,dofe[lcid][lcid]);
00874 delete [] cn;
00875 }
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885 void quadquadrilatt::res_internal_fluxes (long eid,vector &elemif)
00886 {
00887 long i,*cn;
00888 vector lif,tdnv;
00889 matrix cm;
00890
00891 for (i=0;i<ntm;i++){
00892 cn = new long [dofe[i][i]];
00893 reallocv (dofe[i][i],lif);
00894 internal_fluxes (i,eid,lif);
00895 codnum (cn,i);
00896 locglob (elemif.a,lif.a,cn,dofe[i][i]);
00897 delete [] cn;
00898 }
00899
00900 cn = new long [ndofe];
00901 Gtt->give_code_numbers (eid,cn);
00902 reallocv (ndofe,tdnv);
00903 reallocv (ndofe,lif);
00904 reallocm (ndofe,ndofe,cm);
00905 nodalderivatives (tdnv.a,cn,ndofe);
00906 res_capacity_matrix (eid,cm);
00907 mxv (cm,tdnv,lif);
00908
00909 subv (elemif,lif,elemif);
00910
00911 delete [] cn;
00912 }
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924 double quadquadrilatt::total_integral(long eid,vector &nodval)
00925 {
00926 long i,j;
00927 double thick,xi,eta,ww1,ww2,jac,value,f;
00928 ivector nodes(nne);
00929 vector x(nne),y(nne),w(3),gp(3),t(nne);
00930
00931 Tt->give_elemnodes (eid,nodes);
00932 Tc->give_thickness (eid,nodes,t);
00933
00934 Tt->give_node_coord2d (x,y,eid);
00935 gauss_points (gp.a,w.a,3);
00936
00937 f = 0.0;
00938
00939 for (i=0;i<3;i++){
00940 xi=gp[i]; ww1=w[i];
00941 for (j=0;j<3;j++){
00942 eta=gp[j]; ww2=w[j];
00943
00944 jac_2d (jac,x,y,xi,eta);
00945 thick = approx (xi,eta,t);
00946 value = approx (xi,eta,nodval);
00947
00948 jac*=thick*ww1*ww2*value;
00949 f = f + jac;
00950 }
00951 }
00952 return(f);
00953 }
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967 void quadquadrilatt::res_boundary_flux (vector &f,long lcid,long eid,long leid)
00968 {
00969 long *cn,i;
00970 vector lf;
00971
00972 cn = new long [dofe[lcid][lcid]];
00973 codnum (cn,lcid);
00974
00975
00976 for (i=0;i<ntm;i++){
00977 reallocv (dofe[lcid][lcid],lf);
00978
00979 if ((Tt->elements[eid].transi[i]==3)||(Tt->elements[eid].transi[i]==4)){
00980 boundary_flux (lf,i,eid,leid,lcid,i);
00981 }
00982
00983 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
00984 }
00985
00986 delete [] cn;
00987 }
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000 void quadquadrilatt::nod_others (long lcid,long eid,long ri,long ci)
01001 {
01002 long i,j,k,ipp,ncompother,ii;
01003 vector other,h,r(ndofe);
01004 ivector nod(nne),cn(ndofe);
01005
01006 Tt->give_elemnodes (eid,nod);
01007 Tt->give_code_numbers (eid,cn.a);
01008 elemvalues (lcid,eid,r.a,cn.a,ndofe);
01009 ncompother = Tm->givencompother();
01010 reallocv (ncompother,other);
01011 reallocv (ntm,h);
01012
01013 ii = 0;
01014 for (i=0;i<nne;i++){
01015 if (Tp->savemode==0)
01016 ipp=Tt->elements[eid].ipp[ri][ci];
01017 if (Tp->savemode==1)
01018 ipp=Tt->elements[eid].ipp[0][0];
01019
01020 for (j=0;j<ntm;j++){
01021 h[j] = r[ii];
01022 ii++;
01023 }
01024
01025 for(k=0;k<ncompother;k++)
01026 other[k] = Tm->givecompother(k,ipp,h.a);
01027
01028
01029 }
01030
01031 }
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062 void quadquadrilatt::convection_vector (vector &v,long lcid,long eid,long leid,long ri,long ci)
01063 {
01064 long i,ipp;
01065 bocontypet *bc;
01066 ivector nodes(nne);
01067 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
01068 vector list(nned*ned),trc(nned*ned),nodval(nne),av(nne),coef(nne);
01069 matrix km(nne,nne);
01070
01071 fillv (0.0,v);
01072
01073
01074 Tt->give_elemnodes (eid,nodes);
01075
01076 Tc->give_thickness (eid,nodes,t);
01077
01078 Tt->give_node_coord2d (x,y,eid);
01079
01080
01081 gauss_points (gp.a,w.a,intordkm[ri][ci]);
01082
01083
01084 bc = new bocontypet [ned];
01085 Tb->lc[lcid].elemload[leid].give_bc (bc);
01086
01087 Tb->lc[lcid].elemload[leid].give_nodval (lcid,list);
01088
01089
01090 for (i=0;i<nned*ned;i++){
01091 trc[i]=1.0;
01092 }
01093
01094
01095 if (Tp->savemode==0)
01096 ipp=Tt->elements[eid].ipp[ri][ci];
01097 if (Tp->savemode==1)
01098 ipp=Tt->elements[eid].ipp[0][0];
01099
01100
01101 for (i=0;i<ned;i++){
01102
01103 if (bc[i]==2 || bc[i]==3){
01104
01105 edgenodeval (i,nodval,list);
01106
01107 edgenodeval (i,coef,trc);
01108
01109 fillm (0.0,km);
01110
01111
01112 edge_integral (i,x,y,intordkm[ri][ci],gp,w,t,coef,km);
01113
01114 mxv (km,nodval,av);
01115
01116 addv (v,av,v);
01117 }
01118 }
01119
01120 delete [] bc;
01121 }
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135 void quadquadrilatt::transmission_matrix (long lcid,long eid,long ri,long ci,matrix &km)
01136 {
01137 long i,ipp,leid;
01138 bocontypet *bc;
01139 ivector nodes(nne);
01140 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]);
01141 vector trc(nned*ned),coeff(nne),t(nne);
01142
01143
01144 Tt->give_elemnodes (eid,nodes);
01145
01146 Tc->give_thickness (eid,nodes,t);
01147
01148 Tt->give_node_coord2d (x,y,eid);
01149
01150
01151 gauss_points (gp.a,w.a,intordkm[ri][ci]);
01152
01153
01154 for (i=0;i<Tb->lc[lcid].neb;i++){
01155 if (Tb->lc[lcid].elemload[i].eid==eid){
01156 leid=i;
01157 break;
01158 }
01159 }
01160
01161
01162 bc = new bocontypet [ned];
01163 Tb->lc[lcid].elemload[leid].give_bc (bc);
01164
01165 Tb->lc[lcid].elemload[leid].give_trc (lcid,ci,trc);
01166
01167
01168
01169 if (Tp->savemode==0)
01170 ipp=Tt->elements[eid].ipp[ri][ci];
01171 if (Tp->savemode==1)
01172 ipp=Tt->elements[eid].ipp[0][0];
01173
01174
01175 for (i=0;i<ned;i++){
01176
01177 if (bc[i]>10){
01178
01179 transf_coeff (i,coeff,trc,eid,ri,ci,ipp,bc);
01180
01181
01182 edge_integral (i,x,y,intordkm[ri][ci],gp,w,t,coeff,km);
01183 }
01184 }
01185
01186 delete [] bc;
01187 }
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203 void quadquadrilatt::transmission_vector (vector &v,long lcid,long eid,long leid,long ri,long ci)
01204 {
01205 long i,ipp;
01206 bocontypet *bc;
01207 ivector nodes(nne);
01208 vector x(nne),y(nne),z(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
01209 vector list(nned*ned),trc(nned*ned),trr(nned*ned),nodval(nne),av(nne),coef(nne);
01210 matrix km(nne,nne);
01211
01212 fillv (0.0,v);
01213
01214
01215 Tt->give_elemnodes (eid,nodes);
01216
01217 Tc->give_thickness (eid,nodes,t);
01218
01219 Tt->give_node_coord2d (x,y,eid);
01220
01221
01222 gauss_points (gp.a,w.a,intordkm[ri][ci]);
01223
01224
01225 bc = new bocontypet [ned];
01226 Tb->lc[lcid].elemload[leid].give_bc (bc);
01227
01228 Tb->lc[lcid].elemload[leid].give_nodval (lcid,list);
01229
01230 Tb->lc[lcid].elemload[leid].give_trc (lcid,ci,trc);
01231
01232 Tb->lc[lcid].elemload[leid].give_trr (lcid,trr);
01233
01234
01235 if (Tp->savemode==0)
01236 ipp=Tt->elements[eid].ipp[ri][ci];
01237 if (Tp->savemode==1)
01238 ipp=Tt->elements[eid].ipp[0][0];
01239
01240 for (i=0;i<ned;i++){
01241
01242 if (bc[i]>10){
01243
01244 transf_val (i,nodval,list,trc,trr,eid,ri,ci,ipp,bc);
01245
01246 edgenodeval (i,coef,trc);
01247
01248 fillm (0.0,km);
01249
01250 edge_integral (i,x,y,intordkm[ri][ci],gp,w,t,coef,km);
01251
01252 mxv (km,nodval,av);
01253
01254 addv (v,av,v);
01255 }
01256 }
01257
01258 delete [] bc;
01259 }
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275 void quadquadrilatt::boundary_flux (vector &v,long lcid,long eid,long leid,long ri,long ci)
01276 {
01277 long i,ipp;
01278 bocontypet *bc;
01279 ivector nodes(nne);
01280 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
01281 vector list(nned*ned),trc(nned*ned),trr(nned*ned),nodval(nne),av(nne),coef(nne);
01282 matrix km(nne,nne);
01283
01284 fillv (0.0,v);
01285
01286
01287 Tt->give_elemnodes (eid,nodes);
01288
01289 Tc->give_thickness (eid,nodes,t);
01290
01291 Tt->give_node_coord2d (x,y,eid);
01292
01293
01294 gauss_points (gp.a,w.a,intordkm[ri][ci]);
01295
01296
01297 bc = new bocontypet [ned];
01298 Tb->lc[lcid].elemload[leid].give_bc (bc);
01299
01300 Tb->lc[lcid].elemload[leid].give_nodval (lcid,list);
01301
01302 Tb->lc[lcid].elemload[leid].give_trc (lcid,ci,trc);
01303
01304 Tb->lc[lcid].elemload[leid].give_trr (lcid,trr);
01305
01306
01307 if (Tp->savemode==0)
01308 ipp=Tt->elements[eid].ipp[ri][ci];
01309 if (Tp->savemode==1)
01310 ipp=Tt->elements[eid].ipp[0][0];
01311
01312 for (i=0;i<ned;i++){
01313
01314 if (bc[i]>10){
01315
01316 transf_flux (i,nodval,list,trc,trr,eid,ri,ci,ipp,bc);
01317
01318 edgenodeval (i,coef,trc);
01319
01320 fillm (0.0,km);
01321
01322 edge_integral (i,x,y,intordkm[ri][ci],gp,w,t,coef,km);
01323
01324 mxv (km,nodval,av);
01325
01326 addv (v,av,v);
01327 }
01328 }
01329
01330 delete [] bc;
01331 }
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347 void quadquadrilatt::edge_integral (long edg,vector &x,vector &y,long intord,vector &gp,vector &w,
01348 vector &t,vector &coef,matrix &km)
01349 {
01350 long i;
01351 double xi,eta,jac,ipval,thick;
01352 matrix n(1,nne);
01353
01354 if (edg==0){
01355 eta=1.0;
01356 for (i=0;i<intord;i++){
01357 xi=gp[i];
01358
01359 bf_matrix (n,xi,eta);
01360 jac1d_2d (jac,x,y,xi,edg);
01361
01362 ipval=approx (xi,eta,coef);
01363 thick=approx (xi,eta,t);
01364
01365 jac*=w[i]*ipval*thick;
01366 nnj (km.a,n.a,jac,n.m,n.n);
01367 }
01368 }
01369
01370 if (edg==1){
01371 xi=-1.0;
01372 for (i=0;i<intord;i++){
01373 eta=gp[i];
01374
01375 bf_matrix (n,xi,eta);
01376 jac1d_2d (jac,x,y,eta,edg);
01377
01378 ipval=approx (xi,eta,coef);
01379 thick=approx (xi,eta,t);
01380
01381 jac*=w[i]*ipval*thick;
01382 nnj (km.a,n.a,jac,n.m,n.n);
01383 }
01384 }
01385
01386 if (edg==2){
01387 eta=-1.0;
01388 for (i=0;i<intord;i++){
01389 xi=gp[i];
01390
01391 bf_matrix (n,xi,eta);
01392 jac1d_2d (jac,x,y,xi,edg);
01393
01394 ipval=approx (xi,eta,coef);
01395 thick=approx (xi,eta,t);
01396
01397 jac*=w[i]*ipval*thick;
01398 nnj (km.a,n.a,jac,n.m,n.n);
01399 }
01400 }
01401
01402 if (edg==3){
01403 xi=1.0;
01404 for (i=0;i<intord;i++){
01405 eta=gp[i];
01406
01407 bf_matrix (n,xi,eta);
01408 jac1d_2d (jac,x,y,eta,edg);
01409
01410 ipval=approx (xi,eta,coef);
01411 thick=approx (xi,eta,t);
01412
01413 jac*=w[i]*ipval*thick;
01414 nnj (km.a,n.a,jac,n.m,n.n);
01415 }
01416 }
01417
01418 }
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433 void quadquadrilatt::transf_flux (long edg,vector &coeff,vector &list,vector &trc,vector &trr,long eid,long ri,long ci,long ipp,bocontypet *bc)
01434 {
01435 long i,j,k;
01436 double tr,new_flux;
01437 ivector nodes(nne),edgenod(nned);
01438
01439
01440 fillv (0.0,coeff);
01441
01442
01443 Tt->give_elemnodes (eid,nodes);
01444
01445
01446 quadquadrilat_edgnod (edgenod.a,edg);
01447
01448
01449 j=edg*nned;
01450
01451 for (i=0;i<nned;i++){
01452 if (bc[edg]==90){
01453 tr = trr[j+i]/trc[j+i];
01454 }
01455
01456
01457 k=edgenod[i];
01458
01459 Tm->transmission_flux(new_flux,list[j+i],tr,ri,ci,nodes[k],bc[edg],ipp);
01460
01461 coeff[k]=new_flux;
01462 }
01463
01464 }
01465
01466
01467 void quadquadrilatt::transf_coeff (long edg,vector &coeff,vector &list,long eid,long ri,long ci,long ipp,bocontypet *bc)
01468 {
01469 long i,j,k;
01470 double new_coeff;
01471 ivector nodes(nne),edgenod(nned);
01472
01473
01474 fillv (0.0,coeff);
01475
01476
01477 Tt->give_elemnodes (eid,nodes);
01478
01479
01480 quadquadrilat_edgnod (edgenod.a,edg);
01481
01482
01483 j=edg*nned;
01484
01485 for (i=0;i<nned;i++){
01486
01487 k=edgenod[i];
01488
01489 Tm->transmission_transcoeff(new_coeff,list[j+i],ri,ci,nodes[k],bc[edg],ipp);
01490
01491 coeff[k]=new_coeff;
01492 }
01493
01494 }
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509 void quadquadrilatt::transf_val (long edg,vector &nodval,vector &list,vector &trc,vector &trr,long eid,long ri,long ci,long ipp,bocontypet *bc)
01510 {
01511 long i,j,k;
01512 double tr,new_nodval;
01513 ivector nodes(nne),edgenod(nned);
01514
01515
01516 fillv (0.0,nodval);
01517
01518
01519 Tt->give_elemnodes (eid,nodes);
01520
01521
01522 quadquadrilat_edgnod (edgenod.a,edg);
01523
01524
01525 j=edg*nned;
01526
01527
01528 for (i=0;i<nned;i++){
01529 if (bc[edg]==90){
01530 tr = trr[j+i]/trc[j+i];
01531 }
01532
01533
01534 k=edgenod[i];
01535
01536 Tm->transmission_nodval(new_nodval,list[j+i],tr,ri,ci,nodes[k],bc[edg],ipp);
01537
01538
01539 nodval[k]=new_nodval;
01540 }
01541
01542 }
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553 void quadquadrilatt::edgenodeval (long edg,vector &nodval,vector &list)
01554 {
01555 long i,j;
01556 ivector edgenod(nned);
01557
01558 fillv (0.0,nodval);
01559 quadquadrilat_edgnod (edgenod.a,edg);
01560
01561 for (i=0;i<nned;i++){
01562 j=edgenod[i];
01563 nodval[j]=list[edg*nned+i];
01564 }
01565 }
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576 void quadquadrilatt::intpointflux (long eid)
01577 {
01578 long i,ii,jj,k,ipp;
01579
01580 for (k=0;k<Tp->ntm;k++){
01581
01582 for (ii=0;ii<Tp->ntm;ii++){
01583 for (jj=0;jj<Tp->ntm;jj++){
01584
01585 ipp=Tt->elements[eid].ipp[ii][jj];
01586
01587 for (i=0;i<intordkm[ii][jj];i++){
01588
01589 if (Tp->fluxcomp==1)
01590 Tm->computenlfluxes (k,ipp);
01591
01592 ipp++;
01593 }
01594
01595 for (i=0;i<intordcm[ii][jj];i++){
01596
01597 if (Tp->fluxcomp==1)
01598 Tm->computenlfluxes (k,ipp);
01599
01600 ipp++;
01601 }
01602 if (Tp->savemode==1) break;
01603 }
01604 if (Tp->savemode==1) break;
01605 }
01606 }
01607 }
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619 void quadquadrilatt::nod_grads_ip (long eid)
01620 {
01621 long i,j,k,ipp;
01622 ivector ipnum(nne),nod(nne);
01623 vector grad(ncomp);
01624
01625
01626
01627 ipp=Tt->elements[eid].ipp[0][0];
01628 nodip_planelq (ipp,intordkm[0][0],ipnum);
01629
01630
01631 Tt->give_elemnodes (eid,nod);
01632
01633 for (i=0;i<nne;i++){
01634 for (k=0;k<Tp->ntm;k++){
01635
01636 Tm->givegrad (k,ipnum[i],grad);
01637
01638
01639 j=nod[i];
01640 Tt->nodes[j].storegrad (k,ncomp,grad);
01641 }
01642 }
01643
01644 }
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654 void quadquadrilatt::nod_fluxes_ip (long eid)
01655 {
01656 long i,j,k,ipp;
01657 ivector ipnum(nne),nod(nne);
01658 vector flux(ncomp);
01659
01660
01661
01662 ipp=Tt->elements[eid].ipp[0][0];
01663 nodip_planelq (ipp,intordkm[0][0],ipnum);
01664
01665
01666 Tt->give_elemnodes (eid,nod);
01667
01668 intpointflux (eid);
01669
01670 for (i=0;i<nne;i++){
01671 for (k=0;k<Tp->ntm;k++){
01672
01673 Tm->givefluxes (k,ipnum[i],flux);
01674
01675
01676 j=nod[i];
01677 Tt->nodes[j].storeflux (k,flux);
01678 }
01679 }
01680
01681 }
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694 void quadquadrilatt::nod_others_comp (long lcid,long eid,long ri,long ci)
01695 {
01696 long i,j,k,ipp,ncompother,ii;
01697 vector h,r(ndofe);
01698 ivector nod(nne),cn(ndofe);
01699 double other;
01700
01701 Tt->give_elemnodes (eid,nod);
01702 Tt->give_code_numbers (eid,cn.a);
01703 elemvalues (lcid,eid,r.a,cn.a,ndofe);
01704 ncompother = Tm->givencompother();
01705 reallocv (ntm,h);
01706
01707 ii = 0;
01708 for (i=0;i<nne;i++){
01709 if (Tp->savemode==0)
01710 ipp=Tt->elements[eid].ipp[ri][ci];
01711 if (Tp->savemode==1)
01712 ipp=Tt->elements[eid].ipp[0][0];
01713
01714 for (j=0;j<ntm;j++){
01715 h[j] = r[ii];
01716 ii++;
01717 }
01718
01719 for(k=0;k<ncompother;k++){
01720 other = Tm->givecompother(k,ipp,h.a);
01721 Tt->nodes[nod[i]].storeother (k,other);
01722 }
01723 }
01724
01725 }