00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include <limits.h>
00005 #include <math.h>
00006 #include "mechtop.h"
00007 #include "global.h"
00008 #include "elemhead.h"
00009 #include "sequent.h"
00010 #include "gtopology.h"
00011 #include "difcalc.h"
00012 #include "mathem.h"
00013 #include "element.h"
00014 #include "node.h"
00015 #include "endnodem.h"
00016 #include "edgem.h"
00017 #include "siftop_element_types.h"
00018 #include "loadcase.h"
00019 #include "intpoints.h"
00020 #include "globmat.h"
00021 #include "siftop.h"
00022 #include "saddpoint.h"
00023 #include "vecttens.h"
00024
00025
00026
00027
00028
00029
00030
00031
00032 mechtop::mechtop ()
00033 {
00034
00035 nn=0;
00036
00037 ncn=0;
00038
00039 ne=0;
00040
00041 nln=0;
00042
00043
00044 max_ndofn = -1;
00045
00046
00047 nodes=NULL;
00048
00049 elements=NULL;
00050
00051 dist=NULL;
00052 nadjip=NULL;
00053 adjip=NULL;
00054 nodedispl=NULL;
00055 nodeforce=NULL;
00056 domvol=0.0;
00057 }
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 mechtop::~mechtop ()
00068 {
00069 long i;
00070
00071 delete [] nodes;
00072 delete [] elements;
00073 delete [] nadjip;
00074 if (adjip)
00075 {
00076 for (i=0; i<Mm->tnip; i++)
00077 delete [] adjip[i];
00078 }
00079 if (dist)
00080 {
00081 for (i=0; i<Mm->tnip; i++)
00082 delete [] dist[i];
00083 }
00084 delete [] adjip;
00085 delete [] dist;
00086
00087 if (nodedispl != NULL){
00088 for (i=0;i<nn;i++){
00089 delete [] nodedispl[i];
00090 }
00091 delete [] nodedispl;
00092 }
00093 if (nodeforce != NULL){
00094 for (i=0;i<nn;i++){
00095 delete [] nodeforce[i];
00096 }
00097 delete [] nodeforce;
00098 }
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 void mechtop::read (XFILE *in)
00113 {
00114 long i,j,k,l,ndofn;
00115
00116
00117
00118
00119 xfscanf (in,"%k%ld","number_of_nodes",&nn);
00120 Gtm->alloc_nodes (nn);
00121
00122 if (Mespr==1) fprintf (stdout,"\n number of nodes %ld",nn);
00123 nodes = new node [nn];
00124
00125 switch (Mp->tprob){
00126 case linear_statics:
00127 case eigen_dynamics:
00128 case forced_dynamics:
00129 case linear_stability:
00130 case mat_nonlinear_statics:
00131 case geom_nonlinear_statics:
00132 case mech_timedependent_prob:
00133 case growing_mech_structure:
00134 case load_balancing:
00135 case lin_floating_subdomain:
00136 case hemivar_inequal:
00137 case earth_pressure:{
00138 for (i=0;i<nn;i++){
00139 j=nn+1;
00140 xfscanf (in,"%ld",&j);
00141 if (j<1)
00142 print_err("node number in function read is less than 1",__FILE__,__LINE__,__func__);
00143 if (j>nn)
00144 print_err("node number in function read is greater than number of nodes",__FILE__,__LINE__,__func__);
00145 Gtm->gnodes[j-1].read (in);
00146 ndofn = Gtm->gnodes[j-1].ndofn;
00147 nodes[j-1].read (in,ndofn);
00148 }
00149 break;
00150 }
00151 case layered_linear_statics:{
00152 xfscanf (in,"%k%ld","number_of_layered_nodes",&nln);
00153 Gtm->alloc_lnodes (nln);
00154 for (i=0;i<nln;i++){
00155 xfscanf (in,"%ld",&l);
00156 xfscanf (in,"%ld",&(Gtm->lgnodes[l-1].nl));
00157 Gtm->lgnodes[l-1].nodes = new long [Gtm->lgnodes[l-1].nl];
00158 for (k=0;k<Gtm->lgnodes[l-1].nl;k++){
00159 j=nn+1;
00160 xfscanf (in,"%ld",&j);
00161 if (j<1)
00162 print_err("node number in function read is less than 1",__FILE__,__LINE__,__func__);
00163 if (j>nn)
00164 print_err("node number in function read is greater than number of nodes",__FILE__,__LINE__,__func__);
00165 Gtm->gnodes[j-1].read (in);
00166 ndofn = Gtm->gnodes[j-1].ndofn;
00167 nodes[j-1].read (in,ndofn);
00168 Gtm->lgnodes[l-1].nodes[k]=j-1;
00169 }
00170 }
00171 Gtm->unodelnode ();
00172 break;
00173 }
00174
00175 case nonlin_floating_subdomain:{
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 for (i=0;i<nn;i++){
00193 j=nn+1;
00194 xfscanf (in,"%ld",&j);
00195 if (j<1)
00196 print_err("node number in function read is less than 1",__FILE__,__LINE__,__func__);
00197 if (j>nn)
00198 print_err("node number in function read is greater than number of nodes",__FILE__,__LINE__,__func__);
00199 Gtm->gnodes[j-1].read (in);
00200 ndofn = Gtm->gnodes[j-1].ndofn;
00201 nodes[j-1].read (in,ndofn);
00202 }
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224 break;
00225 }
00226
00227 default:{
00228 print_err("unknown type of problem is required",__FILE__,__LINE__,__func__);
00229 }
00230 }
00231
00232
00233 xfscanf (in,"%k%ld","number_of_constraints",&ncn);
00234 if (Mespr==1) fprintf (stdout,"\n number of constrained nodes %ld",ncn);
00235 if (Gtm->dofcontr)
00236 {
00237
00238
00239 nodedispl = new double* [nn];
00240 memset(nodedispl, 0, sizeof(*nodedispl)*nn);
00241 }
00242 for (i=0;i<ncn;i++){
00243 k=nn+1;
00244 xfscanf (in,"%ld",&k);
00245 if (k<1)
00246 print_err("number of constrained node in function read is less than 1",__FILE__,__LINE__,__func__);
00247 if (k>nn){
00248 print_err("number of constrained node in function read is greater than number of all nodes",__FILE__,__LINE__,__func__);
00249 }
00250 Gtm->gnodes[k-1].constr(Gtm->dofcontr,in);
00251 if (Gtm->dofcontr)
00252 {
00253 ndofn = give_ndofn(k-1);
00254 nodedispl[k-1] = new double[ndofn];
00255 memset(nodedispl[k-1], 0, sizeof(*nodedispl[k-1])*ndofn);
00256 }
00257 }
00258
00259
00260
00261
00262
00263
00264 xfscanf (in,"%k%ld","number_of_elements",&ne);
00265 Gtm->alloc_elements (ne);
00266 if (Mespr==1) fprintf (stdout,"\n number of elements %ld",ne);
00267 elements = new element [ne];
00268
00269 switch (Mp->tprob){
00270 case linear_statics:
00271 case eigen_dynamics:
00272 case forced_dynamics:
00273 case linear_stability:
00274 case mat_nonlinear_statics:
00275 case geom_nonlinear_statics:
00276 case mech_timedependent_prob:
00277 case growing_mech_structure:
00278 case earth_pressure:
00279 case lin_floating_subdomain:
00280 case hemivar_inequal:
00281 case load_balancing:
00282 case layered_linear_statics:{
00283 for (i=0;i<ne;i++){
00284 j=ne+1;
00285 xfscanf (in,"%ld",&j);
00286 if (j<1)
00287 print_err("element number in function read is less than 1",__FILE__,__LINE__,__func__);
00288 if (j>ne){
00289 print_err("element number in function read is greater than total number of elements",__FILE__,__LINE__,__func__);
00290 }
00291 elements[j-1].read (in,j-1);
00292 }
00293 break;
00294 }
00295
00296 case nonlin_floating_subdomain:{
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312 for (i=0;i<ne;i++){
00313 j=ne+1;
00314 xfscanf (in,"%ld",&j);
00315 if ((j<1) || (j>ne))
00316 print_err ("element number is out of range <1,%ld>.",__FILE__,__LINE__, __func__, ne);
00317 elements[j-1].read (in,j-1);
00318 }
00319 break;
00320 }
00321 default:{
00322 print_err("unknown type of problem is required",__FILE__,__LINE__,__func__);
00323 }
00324 }
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357 if (Mp->homog != 3 && Mp->homog !=4)
00358 {
00359 Gtm->searching_hanging_nodes (Out);
00360 searching_hanging_nodes (Out);
00361 Gtm->hang_nodes_check (Out);
00362 }
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381 Gtm->lneso_init ();
00382
00383 if (Mp->tprob == growing_mech_structure){
00384 Gtm->auxinf_init ();
00385 }
00386
00387 if (Gtm->rst==1){
00388
00389
00390
00391 Gtm->read_seq_top (in);
00392 }
00393 if (Gtm->rst==2){
00394
00395
00396 Gtm->read_seq_top (in);
00397 Mp->ssle->feti.read_booldata (in);
00398 }
00399
00400 }
00401
00402
00403
00404
00405
00406
00407 void mechtop::searching_hanging_nodes (FILE *Out)
00408 {
00409 long i,j,k,l,nmn,npn,nneo,nnen,ndofeo,ndofen,ndofn,ri,ci,cumulci;
00410 ivector oldnodes,newnodes;
00411 vector weights;
00412 matrix tmat;
00413
00414
00415 for (i=0;i<ne;i++){
00416 if (give_ndofe(i)!=Gtm->give_ndofe(i)){
00417
00418
00419
00420 nneo = give_nne (i);
00421
00422 nnen = Gtm->give_nmne (i);
00423
00424 reallocv (nneo,oldnodes);
00425 reallocv (nnen,newnodes);
00426
00427
00428 Gtm->give_original_nodes (i,oldnodes);
00429
00430 Gtm->give_master_nodes (i,newnodes);
00431
00432
00433 ndofeo = give_ndofe (i);
00434
00435 ndofen = Gtm->give_ndofe (i);
00436
00437
00438 reallocm (ndofeo,ndofen,tmat);
00439
00440
00441 npn=0;
00442
00443 ri=0;
00444 cumulci=0;
00445
00446 for (j=0;j<nneo;j++){
00447
00448 ci=cumulci;
00449
00450
00451 ndofn = give_ndofn (oldnodes[j]);
00452
00453 if (oldnodes[j]==newnodes[npn]){
00454
00455
00456
00457 for (k=0;k<ndofn;k++){
00458 tmat[ri][ci]=1.0;
00459 ri++;
00460 ci++;
00461 cumulci++;
00462 }
00463 npn++;
00464 }
00465 else{
00466
00467
00468
00469 nmn=0-Gtm->gnodes[oldnodes[j]].ndofn;
00470
00471 reallocv (nmn,weights);
00472 Gtm->approx_weights (Gtm->gnodes[oldnodes[j]].masentity,oldnodes[j],weights,Out);
00473
00474
00475 ndofn = give_ndofn (newnodes[npn]);
00476
00477
00478 for (k=0;k<ndofn;k++){
00479 ci=cumulci+k;
00480
00481
00482 for (l=0;l<nmn;l++){
00483
00484 tmat[ri][ci]=weights[l];
00485 ci+=ndofn;
00486 }
00487 ri++;
00488 }
00489 npn+=nmn;
00490 cumulci=ci+1-ndofn;
00491 }
00492
00493 }
00494
00495 elements[i].tmat = new matrix (ndofeo,ndofen);
00496 copym (tmat,*(elements[i].tmat));
00497
00498
00499
00500
00501 }
00502 }
00503 }
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514 void mechtop::elemprescdisp (void)
00515 {
00516 long i,j,k,ndofe;
00517 long nne, ndofn;
00518 ivector enod, cn;
00519
00520 for (i=0;i<ne;i++){
00521 if (Gtm->leso[i]==1){
00522 ndofe=Gtm->give_ndofe (i);
00523 reallocv(ndofe, cn);
00524 give_code_numbers (i,cn.a);
00525 for (j=0;j<ndofe;j++){
00526 if (cn[j]<0){
00527 elements[i].prescdispl=1;
00528 break;
00529 }
00530 }
00531 }
00532 if ((Mp->tprob == growing_mech_structure) && (elements[i].prescdispl==0))
00533 {
00534
00535
00536 nne = give_nne(i);
00537 reallocv(nne, enod);
00538 give_elemnodes(i, enod);
00539 for(j=0; j<nne; j++)
00540 {
00541 if (nodedispl[enod[j]])
00542 {
00543 ndofn = give_ndofn(enod[j]);
00544 for(k=0; k<ndofn; k++)
00545 {
00546 if (nodedispl[enod[j]][k] != 0.0)
00547 {
00548
00549 elements[i].prescdispl=1;
00550 break;
00551 }
00552 }
00553 }
00554 if(elements[i].prescdispl)
00555 {
00556
00557
00558 break;
00559 }
00560 }
00561 }
00562 }
00563 }
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576 void mechtop::elempresctemp (long lcid)
00577 {
00578 long i;
00579 ivector nod;
00580
00581 if (Mp->temperature > 0){
00582 for (i=0;i<ne;i++){
00583 if (Gtm->leso[i]==1){
00584 elements[i].presctemp=1;
00585 }
00586 }
00587 }
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608 }
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621 elemtype mechtop::give_elem_type (long eid)
00622 {
00623 return (elements[eid].te);
00624 }
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637 long mechtop::give_ndofe (long eid)
00638 {
00639 long ndofe=0;
00640 elemtype te;
00641
00642 te = give_elem_type (eid);
00643
00644 switch (te){
00645 case bar2d:{ ndofe=Bar2d->ndofe; break; }
00646 case bar3d:{ ndofe=Bar3d->ndofe; break; }
00647 case barq2d:{ ndofe=Barq2d->ndofe; break; }
00648 case barq3d:{ ndofe=Barq3d->ndofe; break; }
00649 case beam2d:{ ndofe=Beam2d->ndofe; break; }
00650 case beam3d:{ ndofe=Beam3d->ndofe; break; }
00651 case beamg3d:{ ndofe=Beam3dg->ndofe; break; }
00652 case subsoilbeam:{ ndofe=Sbeam->ndofe; break; }
00653 case beam2dsp:{ ndofe=Spbeam2d->ndofe; break; }
00654 case spring_1:{ ndofe=Spring->give_ndofe(eid); break; }
00655 case spring_2:{ ndofe=Spring->give_ndofe(eid); break; }
00656 case spring_3:{ ndofe=Spring->give_ndofe(eid); break; }
00657 case spring_4:{ ndofe=Spring->give_ndofe(eid); break; }
00658 case spring_5:{ ndofe=Spring->give_ndofe(eid); break; }
00659 case spring_6:{ ndofe=Spring->give_ndofe(eid); break; }
00660 case planeelementlt:{ ndofe=Pelt->ndofe; break; }
00661 case planeelementqt:{ ndofe=Peqt->ndofe; break; }
00662 case planeelementrotlt:{ ndofe=Perlt->ndofe; break; }
00663 case planeelementlq:{ ndofe=Pelq->ndofe; break; }
00664 case planeelementqq:{ ndofe=Peqq->ndofe; break; }
00665 case planeelementrotlq:{ ndofe=Perlq->ndofe; break; }
00666 case planeelementsubqt:{ ndofe=Pesqt->ndofe; break; }
00667 case planequadcontact:{ ndofe=Pqcon->ndofe; break; }
00668 case cctel:{ ndofe=Cct->ndofe; break; }
00669 case dktel:{ ndofe=Dkt->ndofe; break; }
00670 case dstel:{ ndofe=Dst->ndofe; break; }
00671 case q4plateel:{ ndofe=Q4pl->ndofe; break; }
00672
00673 case argyristr:{ ndofe=Argtrpl->ndofe; break; }
00674 case subsoilplatetr:{ ndofe=Spltr->ndofe; break; }
00675 case subsoilplateq:{ ndofe=Splq->ndofe; break; }
00676 case axisymmlt:{ ndofe=Asymlt->ndofe; break; }
00677 case axisymmlq:{ ndofe=Asymlq->ndofe; break; }
00678 case axisymmqq:{ ndofe=Asymqq->ndofe; break; }
00679 case shelltrelem:{ ndofe=Shtr->ndofe; break; }
00680 case shellqelem:{ ndofe=Shq->ndofe; break; }
00681 case lineartet:{ ndofe=Ltet->ndofe; break; }
00682 case quadrtet:{ ndofe=Qtet->ndofe; break; }
00683 case linearhex:{ ndofe=Lhex->ndofe; break; }
00684 case quadrhex:{ ndofe=Qhex->ndofe; break; }
00685 case lineartetrot:{ ndofe=Ltetrot->ndofe; break; }
00686 case linearhexrot:{ ndofe=Lhexrot->ndofe; break; }
00687 case linearwed:{ ndofe=Lwed->ndofe; break; }
00688 case quadrwed:{ ndofe=Qwed->ndofe; break; }
00689 case particleelem:{ ndofe=Pelem->ndofe; break; }
00690 default:{
00691 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00692 }
00693 }
00694 return ndofe;
00695 }
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708 long mechtop::give_ndofn (long nid)
00709 {
00710
00711
00712 long ndofn;
00713
00714 ndofn = Gtm->give_ndofn (nid);
00715 if (ndofn<0){
00716
00717 ndofn = Gtm->give_ndofn (Gtm->gnodes[nid].mnodes[0]);
00718 }
00719
00720 return ndofn;
00721 }
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735 long mechtop::give_dof (long nid,long n)
00736 {
00737 return Gtm->give_dof (nid,n);
00738 }
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752 void mechtop::give_nodal_coord (long nid, vector &c)
00753 {
00754 c[0] = Gtm->gnodes[nid].x;
00755 c[1] = Gtm->gnodes[nid].y;
00756 c[2] = Gtm->gnodes[nid].z;
00757 }
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771 long mechtop::give_nne (long eid)
00772 {
00773 long nne;
00774 elemtype te;
00775
00776 te = give_elem_type (eid);
00777
00778 switch (te){
00779 case bar2d:{ nne=Bar2d->nne; break; }
00780 case bar3d:{ nne=Bar3d->nne; break; }
00781 case barq2d:{ nne=Barq2d->nne; break; }
00782 case barq3d:{ nne=Barq3d->nne; break; }
00783 case beam2d:{ nne=Beam2d->nne; break; }
00784 case beam3d:{ nne=Beam3d->nne; break; }
00785 case beamg3d:{ nne=Beam3dg->nne; break; }
00786 case subsoilbeam:{ nne=Sbeam->nne; break; }
00787 case beam2dsp:{ nne=Spbeam2d->nne; break; }
00788 case spring_1:{ nne=Spring->nne; break; }
00789 case spring_2:{ nne=Spring->nne; break; }
00790 case spring_3:{ nne=Spring->nne; break; }
00791 case spring_4:{ nne=Spring->nne; break; }
00792 case spring_5:{ nne=Spring->nne; break; }
00793 case spring_6:{ nne=Spring->nne; break; }
00794 case planeelementlt:{ nne=Pelt->nne; break; }
00795 case planeelementqt:{ nne=Peqt->nne; break; }
00796 case planeelementrotlt:{ nne=Perlt->nne; break; }
00797 case planeelementlq:{ nne=Pelq->nne; break; }
00798 case planeelementqq:{ nne=Peqq->nne; break; }
00799 case planeelementrotlq:{ nne=Perlq->nne; break; }
00800 case planeelementsubqt:{ nne=Pesqt->nne; break; }
00801 case planequadcontact:{ nne=Pqcon->nne; break; }
00802 case cctel:{ nne=Cct->nne; break; }
00803 case dktel:{ nne=Dkt->nne; break; }
00804 case dstel:{ nne=Dst->nne; break; }
00805 case q4plateel:{ nne=Q4pl->nne; break; }
00806
00807 case argyristr:{ nne=Argtrpl->nne; break; }
00808 case subsoilplatetr:{ nne=Spltr->nne; break; }
00809 case subsoilplateq:{ nne=Splq->nne; break; }
00810 case axisymmlt:{ nne=Asymlt->nne; break; }
00811 case axisymmlq:{ nne=Asymlq->nne; break; }
00812 case axisymmqq:{ nne=Asymqq->nne; break; }
00813 case shelltrelem:{ nne=Shtr->nne; break; }
00814 case shellqelem:{ nne=Shq->nne; break; }
00815 case lineartet:{ nne=Ltet->nne; break; }
00816 case quadrtet:{ nne=Qtet->nne; break; }
00817 case linearhex:{ nne=Lhex->nne; break; }
00818 case quadrhex:{ nne=Qhex->nne; break; }
00819 case lineartetrot:{ nne=Ltetrot->nne; break; }
00820 case linearhexrot:{ nne=Lhexrot->nne; break; }
00821 case linearwed:{ nne=Lwed->nne; break; }
00822 case quadrwed:{ nne=Qwed->nne; break; }
00823 case particleelem:{ nne=Pelem->nne; break; }
00824 default:{
00825 print_err("unknown element type is required",__FILE__,__LINE__, __func__);
00826 }
00827 }
00828 return nne;
00829 }
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843 long mechtop::give_nip (long eid,long ri,long ci)
00844 {
00845 long nip=0;
00846 elemtype te;
00847
00848 te = give_elem_type (eid);
00849
00850 switch (te){
00851
00852 case bar2d:{ nip=Bar2d->nip[ri][ci]; break; }
00853 case bar3d:{ nip=Bar3d->nip[ri][ci]; break; }
00854 case barq2d:{ nip=Barq2d->nip[ri][ci]; break; }
00855 case barq3d:{ nip=Barq3d->nip[ri][ci]; break; }
00856 case beam2d:{ nip=Beam2d->nip[ri][ci]; break; }
00857 case beam3d:{ nip=Beam3d->nip[ri][ci]; break; }
00858 case beamg3d:{ nip=Beam3dg->nip[ri][ci]; break; }
00859 case subsoilbeam:{ nip=Sbeam->nip[ri][ci]; break; }
00860 case beam2dsp:{ nip=Spbeam2d->nip[ri][ci]; break; }
00861 case spring_1:{ nip=Spring->nip[ri][ci]; break; }
00862 case spring_2:{ nip=Spring->nip[ri][ci]; break; }
00863 case spring_3:{ nip=Spring->nip[ri][ci]; break; }
00864 case spring_4:{ nip=Spring->nip[ri][ci]; break; }
00865 case spring_5:{ nip=Spring->nip[ri][ci]; break; }
00866 case spring_6:{ nip=Spring->nip[ri][ci]; break; }
00867 case planeelementlt:{ nip=Pelt->nip[ri][ci]; break; }
00868 case planeelementqt:{ nip=Peqt->nip[ri][ci]; break; }
00869 case planeelementrotlt:{ nip=Perlt->nip[ri][ci]; break; }
00870 case planeelementlq:{ nip=Pelq->nip[ri][ci]; break; }
00871 case planeelementqq:{ nip=Peqq->nip[ri][ci]; break; }
00872 case planeelementrotlq:{ nip=Perlq->nip[ri][ci]; break; }
00873 case planeelementsubqt:{ nip=Pesqt->nip[ri][ci]; break; }
00874 case planequadcontact:{ nip=Pqcon->nip[ri][ci]; break; }
00875 case cctel:{ nip=Cct->nip[ri][ci]; break; }
00876 case dktel:{ nip=Dkt->nip[ri][ci]; break; }
00877 case dstel:{ nip=Dst->nip[ri][ci]; break; }
00878 case q4plateel:{ nip=Q4pl->nip[ri][ci]; break; }
00879
00880 case argyristr:{ nip=Argtrpl->nip[ri][ci]; break; }
00881 case subsoilplatetr:{ nip=Spltr->nip[ri][ci]; break; }
00882 case subsoilplateq:{ nip=Splq->nip[ri][ci]; break; }
00883 case axisymmlt:{ nip=Asymlt->nip[ri][ci]; break; }
00884 case axisymmlq:{ nip=Asymlq->nip[ri][ci]; break; }
00885 case axisymmqq:{ nip=Asymqq->nip[ri][ci]; break; }
00886 case shelltrelem:{ nip=Shtr->nip[ri][ci]; break; }
00887 case shellqelem:{ nip=Shq->nip[ri][ci]; break; }
00888 case lineartet:{ nip=Ltet->nip[ri][ci]; break; }
00889 case quadrtet:{ nip=Qtet->nip[ri][ci]; break; }
00890 case linearhex:{ nip=Lhex->nip[ri][ci]; break; }
00891 case quadrhex:{ nip=Qhex->nip[ri][ci]; break; }
00892 case lineartetrot:{ nip=Ltetrot->nip[ri][ci]; break; }
00893 case linearhexrot:{ nip=Lhexrot->nip[ri][ci]; break; }
00894 case linearwed:{ nip=Lwed->nip[ri][ci]; break; }
00895 case quadrwed:{ nip=Qwed->nip[ri][ci]; break; }
00896 case particleelem:{ nip=Pelem->nip[ri][ci]; break; }
00897 default:{
00898 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00899 }
00900 }
00901 return nip;
00902 }
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916 long mechtop::give_intordsm (long eid,long ri,long ci)
00917 {
00918 long ret=0;
00919 elemtype te;
00920
00921 te = give_elem_type (eid);
00922
00923 switch (te){
00924
00925 case bar2d:{ ret=Bar2d->intordsm[ri][ci]; break; }
00926 case bar3d:{ ret=Bar3d->intordsm[ri][ci]; break; }
00927 case barq2d:{ ret=Barq2d->intordsm[ri][ci]; break; }
00928 case barq3d:{ ret=Barq3d->intordsm[ri][ci]; break; }
00929 case beam2d:{ ret=Beam2d->intordsm[ri][ci]; break; }
00930 case beam3d:{ ret=Beam3d->intordsm[ri][ci]; break; }
00931 case beamg3d:{ ret=Beam3dg->intordsm[ri][ci]; break; }
00932 case subsoilbeam:{ ret=Sbeam->intordsm[ri][ci]; break; }
00933 case spring_1:{ ret=Spring->intordsm[ri][ci]; break; }
00934 case spring_2:{ ret=Spring->intordsm[ri][ci]; break; }
00935 case spring_3:{ ret=Spring->intordsm[ri][ci]; break; }
00936 case spring_4:{ ret=Spring->intordsm[ri][ci]; break; }
00937 case spring_5:{ ret=Spring->intordsm[ri][ci]; break; }
00938 case spring_6:{ ret=Spring->intordsm[ri][ci]; break; }
00939 case planeelementlt:{ ret=Pelt->intordsm[ri][ci]; break; }
00940 case planeelementqt:{ ret=Peqt->intordsm[ri][ci]; break; }
00941 case planeelementrotlt:{ ret=Perlt->intordsm[ri][ci]; break; }
00942 case planeelementlq:{ ret=Pelq->intordsm[ri][ci]; break; }
00943 case planeelementqq:{ ret=Peqq->intordsm[ri][ci]; break; }
00944 case planeelementrotlq:{ ret=Perlq->intordsm[ri][ci]; break; }
00945 case planeelementsubqt:{ ret=Pesqt->intordsm[ri][ci]; break; }
00946 case planequadcontact:{ ret=Pqcon->intordsm[ri][ci]; break; }
00947 case cctel:{ ret=Cct->intordsm[ri][ci]; break; }
00948 case dktel:{ ret=Dkt->intordsm[ri][ci]; break; }
00949 case dstel:{ ret=Dst->intordsm[ri][ci]; break; }
00950 case q4plateel:{ ret=Q4pl->intordsm[ri][ci]; break; }
00951
00952 case argyristr:{ ret=Argtrpl->intordsm[ri][ci]; break; }
00953 case subsoilplatetr:{ ret=Spltr->intordsm[ri][ci]; break; }
00954 case subsoilplateq:{ ret=Splq->intordsm[ri][ci]; break; }
00955 case axisymmlt:{ ret=Asymlt->intordsm[ri][ci]; break; }
00956 case axisymmlq:{ ret=Asymlq->intordsm[ri][ci]; break; }
00957 case axisymmqq:{ ret=Asymqq->intordsm[ri][ci]; break; }
00958 case shelltrelem:{ ret=Shtr->intordsm[ri][ci]; break; }
00959 case shellqelem:{ ret=Shq->intordsm[ri][ci]; break; }
00960 case lineartet:{ ret=Ltet->intordsm[ri][ci]; break; }
00961 case quadrtet:{ ret=Qtet->intordsm[ri][ci]; break; }
00962 case linearhex:{ ret=Lhex->intordsm[ri][ci]; break; }
00963 case quadrhex:{ ret=Qhex->intordsm[ri][ci]; break; }
00964 case lineartetrot:{ ret=Ltetrot->intordsm[ri][ci]; break; }
00965 case linearhexrot:{ ret=Lhexrot->intordsm[ri][ci]; break; }
00966 default:{
00967 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00968 }
00969 }
00970 return ret;
00971 }
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984 long mechtop::give_tnip (long eid)
00985 {
00986 long tnip = -1;
00987 elemtype te;
00988
00989 te = give_elem_type (eid);
00990
00991 switch (te){
00992 case bar2d: { tnip=Bar2d->tnip; break; }
00993 case bar3d: { tnip=Bar3d->tnip; break; }
00994 case barq2d: { tnip=Barq2d->tnip; break; }
00995 case barq3d: { tnip=Barq3d->tnip; break; }
00996 case beam2d:{ tnip=Beam2d->tnip; break; }
00997 case beam3d:{ tnip=Beam3d->tnip; break; }
00998 case subsoilbeam:{ tnip=Sbeam->tnip; break; }
00999 case beamg3d:{ tnip=Beam3dg->tnip; break; }
01000 case beam2dsp:{ tnip=Spbeam2d->tnip; break; }
01001 case spring_1:{ tnip=Spring->tnip; break; }
01002 case spring_2:{ tnip=Spring->tnip; break; }
01003 case spring_3:{ tnip=Spring->tnip; break; }
01004 case spring_4:{ tnip=Spring->tnip; break; }
01005 case spring_5:{ tnip=Spring->tnip; break; }
01006 case spring_6:{ tnip=Spring->tnip; break; }
01007 case planeelementlt:{ tnip=Pelt->tnip; break; }
01008 case planeelementqt:{ tnip=Peqt->tnip; break; }
01009 case planeelementrotlt:{ tnip=Perlt->tnip; break; }
01010 case planeelementlq:{ tnip=Pelq->tnip; break; }
01011 case planeelementqq:{ tnip=Peqq->tnip; break; }
01012 case planeelementrotlq:{ tnip=Perlq->tnip; break; }
01013 case planeelementsubqt:{ tnip=Pesqt->tnip; break; }
01014 case planequadcontact:{ tnip=Pqcon->tnip; break; }
01015 case cctel:{ tnip=Cct->tnip; break; }
01016 case dktel:{ tnip=Dkt->tnip; break; }
01017 case dstel:{ tnip=Dst->tnip; break; }
01018 case q4plateel:{ tnip=Q4pl->tnip; break; }
01019
01020 case argyristr:{ tnip=Argtrpl->tnip; break; }
01021 case subsoilplatetr:{ tnip=Spltr->tnip; break; }
01022 case subsoilplateq:{ tnip=Splq->tnip; break; }
01023 case axisymmlt:{ tnip=Asymlt->tnip; break; }
01024 case axisymmlq:{ tnip=Asymlq->tnip; break; }
01025 case axisymmqq:{ tnip=Asymqq->tnip; break; }
01026 case shelltrelem:{ tnip=Shtr->tnip; break; }
01027 case shellqelem:{ tnip=Shq->tnip; break; }
01028 case lineartet:{ tnip=Ltet->tnip; break; }
01029 case quadrtet:{ tnip=Qtet->tnip; break; }
01030 case linearhex:{ tnip=Lhex->tnip; break; }
01031 case quadrhex:{ tnip=Qhex->tnip; break; }
01032 case lineartetrot:{ tnip=Ltetrot->tnip; break; }
01033 case linearhexrot:{ tnip=Lhexrot->tnip; break; }
01034 case linearwed:{ tnip=Lwed->tnip; break; }
01035 case quadrwed:{ tnip=Qwed->tnip; break; }
01036 case particleelem:{ tnip=1; break; }
01037 default:{
01038 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01039 }
01040 }
01041 return tnip;
01042 }
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055 long mechtop::give_totnip (long eid)
01056 {
01057 if (eid < Mt->ne-1)
01058 return elements[eid+1].ipp[0][0] - elements[eid].ipp[0][0];
01059
01060 return Mm->tnip-elements[eid].ipp[0][0];
01061 }
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075 long mechtop::give_ncomp (long eid,long bid)
01076 {
01077 long ncomp=0;
01078 elemtype te;
01079
01080 te = give_elem_type (eid);
01081
01082 switch (te){
01083 case bar2d:{
01084
01085 ncomp=Bar2d->tncomp;
01086 break; }
01087 case bar3d:{ ncomp=Bar3d->tncomp; break; }
01088 case barq2d:{ ncomp=Barq2d->tncomp; break; }
01089 case barq3d:{ ncomp=Barq3d->tncomp; break; }
01090 case beam2d:{ ncomp=Beam2d->tncomp; break; }
01091 case beam3d:{ ncomp=Beam3d->tncomp; break; }
01092 case beamg3d:{ ncomp=Beam3dg->tncomp; break; }
01093 case beam2dsp:{ ncomp=Spbeam2d->tncomp; break; }
01094 case subsoilbeam:{ ncomp=Sbeam->tncomp; break; }
01095 case spring_1:{ ncomp=Spring->tncomp; break; }
01096 case spring_2:{ ncomp=Spring->tncomp; break; }
01097 case spring_3:{ ncomp=Spring->tncomp; break; }
01098 case spring_4:{ ncomp=Spring->tncomp; break; }
01099 case spring_5:{ ncomp=Spring->tncomp; break; }
01100 case spring_6:{ ncomp=Spring->tncomp; break; }
01101 case planeelementlt:{ ncomp=Pelt->tncomp; break; }
01102 case planeelementqt:{ ncomp=Peqt->tncomp; break; }
01103 case planeelementrotlt:{ ncomp=Perlt->tncomp; break; }
01104 case planeelementlq:{ ncomp=Pelq->tncomp; break; }
01105 case planeelementqq:{ ncomp=Peqq->tncomp; break; }
01106 case planeelementrotlq:{ ncomp=Perlq->tncomp; break; }
01107 case planeelementsubqt:{ ncomp=Pesqt->tncomp; break; }
01108 case planequadcontact:{ ncomp=Pqcon->tncomp; break; }
01109 case cctel:{ ncomp=Cct->tncomp; break; }
01110 case dktel:{ ncomp=Dkt->tncomp; break; }
01111 case dstel:{ ncomp=Dst->tncomp; break; }
01112 case q4plateel:{ ncomp=Q4pl->tncomp; break; }
01113
01114 case argyristr:{ ncomp=Argtrpl->tncomp; break; }
01115 case subsoilplatetr:{ ncomp=Spltr->tncomp; break; }
01116 case subsoilplateq:{ ncomp=Splq->tncomp; break; }
01117 case axisymmlt:{ ncomp=Asymlt->tncomp; break; }
01118 case axisymmlq:{ ncomp=Asymlq->tncomp; break; }
01119 case axisymmqq:{ ncomp=Asymqq->tncomp; break; }
01120 case shelltrelem:{
01121 if (bid==0 || bid==1) ncomp=Perlt->tncomp;
01122 if (bid==2) ncomp=Dkt->tncomp;
01123 break;
01124 }
01125 case shellqelem:{
01126 if (bid==0 || bid==1) ncomp=Perlq->tncomp;
01127 if (bid==2 || bid==3) ncomp=Q4pl->tncomp;
01128 break;
01129 }
01130 case lineartet:{ ncomp=Ltet->tncomp; break; }
01131 case quadrtet:{ ncomp=Qtet->tncomp; break; }
01132 case linearhex:{ ncomp=Lhex->tncomp; break; }
01133 case quadrhex:{ ncomp=Qhex->tncomp; break; }
01134 case lineartetrot:{ ncomp=Ltetrot->tncomp; break; }
01135 case linearhexrot:{ ncomp=Lhexrot->tncomp; break; }
01136 case linearwed:{ ncomp=Lwed->tncomp; break; }
01137 case quadrwed:{ ncomp=Qwed->tncomp; break; }
01138 case particleelem:{ ncomp=Pelem->ndofe; break; }
01139 default:{
01140 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01141 }
01142 }
01143 return ncomp;
01144 }
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158 long mechtop::give_tncomp (long eid)
01159 {
01160 long tncomp=0;
01161 elemtype te;
01162
01163 te = give_elem_type (eid);
01164
01165 switch (te){
01166 case bar2d:{ tncomp=Bar2d->tncomp; break; }
01167 case bar3d:{ tncomp=Bar3d->tncomp; break; }
01168 case barq2d:{ tncomp=Barq2d->tncomp; break; }
01169 case barq3d:{ tncomp=Barq3d->tncomp; break; }
01170 case beam2d:{ tncomp=Beam2d->tncomp; break; }
01171 case beam3d:{ tncomp=Beam3d->tncomp; break; }
01172 case beamg3d:{ tncomp=Beam3dg->tncomp; break; }
01173 case beam2dsp:{ tncomp=Spbeam2d->tncomp; break; }
01174 case subsoilbeam:{ tncomp=Sbeam->tncomp; break; }
01175 case spring_1:{ tncomp=Spring->tncomp; break; }
01176 case spring_2:{ tncomp=Spring->tncomp; break; }
01177 case spring_3:{ tncomp=Spring->tncomp; break; }
01178 case spring_4:{ tncomp=Spring->tncomp; break; }
01179 case spring_5:{ tncomp=Spring->tncomp; break; }
01180 case spring_6:{ tncomp=Spring->tncomp; break; }
01181 case planeelementlt:{ tncomp=Pelt->tncomp; break; }
01182 case planeelementqt:{ tncomp=Peqt->tncomp; break; }
01183 case planeelementrotlt:{ tncomp=Perlt->tncomp; break; }
01184 case planeelementlq:{ tncomp=Pelq->tncomp; break; }
01185 case planeelementqq:{ tncomp=Peqq->tncomp; break; }
01186 case planeelementrotlq:{ tncomp=Perlq->tncomp; break; }
01187 case planeelementsubqt:{ tncomp=Pesqt->tncomp; break; }
01188 case planequadcontact:{ tncomp=Pqcon->tncomp; break; }
01189 case cctel:{ tncomp=Cct->tncomp; break; }
01190 case dktel:{ tncomp=Dkt->tncomp; break; }
01191 case dstel:{ tncomp=Dst->tncomp; break; }
01192 case q4plateel:{ tncomp=Q4pl->tncomp; break; }
01193
01194 case argyristr:{ tncomp=Argtrpl->tncomp; break; }
01195 case subsoilplatetr:{ tncomp=Spltr->tncomp; break; }
01196 case subsoilplateq:{ tncomp=Splq->tncomp; break; }
01197 case axisymmlt:{ tncomp=Asymlt->tncomp; break; }
01198 case axisymmlq:{ tncomp=Asymlq->tncomp; break; }
01199 case axisymmqq:{ tncomp=Asymqq->tncomp; break; }
01200 case shelltrelem:{ tncomp=Shtr->tncomp; break; }
01201 case shellqelem:{ tncomp=Shq->tncomp; break; }
01202 case lineartet:{ tncomp=Ltet->tncomp; break; }
01203 case quadrtet:{ tncomp=Qtet->tncomp; break; }
01204 case linearhex:{ tncomp=Lhex->tncomp; break; }
01205 case quadrhex:{ tncomp=Qhex->tncomp; break; }
01206 case lineartetrot:{ tncomp=Ltetrot->tncomp; break; }
01207 case linearhexrot:{ tncomp=Lhexrot->tncomp; break; }
01208 case linearwed:{ tncomp=Lwed->tncomp; break; }
01209 case quadrwed:{ tncomp=Qwed->tncomp; break; }
01210 case particleelem:{ tncomp=Pelem->ndofe; break; }
01211 default:{
01212 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01213 }
01214 }
01215 return tncomp;
01216 }
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229 long mechtop::give_ned (long eid)
01230 {
01231 long ned=0;
01232 elemtype te;
01233
01234 te = give_elem_type (eid);
01235
01236 switch (te){
01237 case bar3d:{ ned=0; break; }
01238 case beam2d:{ ned=1; break; }
01239 case beam3d:{ ned=1; break; }
01240 case planeelementlt:{ ned=Pelt->ned; break; }
01241 case planeelementqt:{ ned=Peqt->ned; break; }
01242 case planeelementrotlt:{ ned=Perlt->ned; break; }
01243 case planeelementlq:{ ned=Pelq->ned; break; }
01244 case planeelementqq:{ ned=Peqq->ned; break; }
01245 case planeelementrotlq:{ ned=Perlq->ned; break; }
01246 case planeelementsubqt:{ ned=Pesqt->ned; break; }
01247 case cctel:{ ned=Cct->ned; break; }
01248 case dktel:{ ned=Dkt->ned; break; }
01249 case dstel:{ ned=Dst->ned; break; }
01250 case q4plateel:{ ned=Q4pl->ned; break; }
01251
01252 case argyristr:{ ned=Argtrpl->ned; break; }
01253 case subsoilplatetr:{ ned=Spltr->ned; break; }
01254 case subsoilplateq:{ ned=Splq->ned; break; }
01255 case axisymmlt:{ ned=Asymlt->ned; break; }
01256 case axisymmlq:{ ned=Asymlq->ned; break; }
01257 case axisymmqq:{ ned=Asymqq->ned; break; }
01258 case shelltrelem:{ ned=Shtr->ned; break; }
01259 case shellqelem:{ ned=Shq->ned; break; }
01260 case lineartet:{ ned=Ltet->ned; break; }
01261 case quadrtet:{ ned=Qtet->ned; break; }
01262 case linearhex:{ ned=Lhex->ned; break; }
01263 case quadrhex:{ ned=Qhex->ned; break; }
01264 case lineartetrot:{ ned=Ltetrot->ned; break; }
01265 case linearhexrot:{ ned=Lhexrot->ned; break; }
01266 case linearwed:{ ned=Lwed->ned; break; }
01267 case quadrwed:{ ned=Qwed->ned; break; }
01268 case particleelem:{ break; }
01269 default:{
01270 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01271 }
01272 }
01273 return ned;
01274 }
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287 long mechtop::give_nned (long eid)
01288 {
01289 long nned=0;
01290 elemtype te;
01291
01292 te = give_elem_type (eid);
01293
01294 switch (te){
01295 case bar3d:{ nned=0; break; }
01296 case beam2d:{ nned=2; break; }
01297 case beam3d:{ nned=2; break; }
01298 case planeelementlt:{ nned=Pelt->nned; break; }
01299 case planeelementqt:{ nned=Peqt->nned; break; }
01300 case planeelementrotlt:{ nned=Perlt->nned; break; }
01301 case planeelementlq:{ nned=Pelq->nned; break; }
01302 case planeelementqq:{ nned=Peqq->nned; break; }
01303 case planeelementrotlq:{ nned=Perlq->nned; break; }
01304 case planeelementsubqt:{ nned=Pesqt->nned; break; }
01305 case cctel:{ nned=Cct->nned; break; }
01306 case dktel:{ nned=Dkt->nned; break; }
01307 case dstel:{ nned=Dst->nned; break; }
01308 case q4plateel:{ nned=Q4pl->nned; break; }
01309
01310 case argyristr:{ nned=Argtrpl->nned; break; }
01311 case subsoilplatetr:{ nned=Spltr->nned; break; }
01312 case subsoilplateq:{ nned=Splq->nned; break; }
01313 case axisymmlt:{ nned=Asymlt->nned; break; }
01314 case axisymmlq:{ nned=Asymlq->nned; break; }
01315 case axisymmqq:{ nned=Asymqq->nned; break; }
01316 case shelltrelem:{ nned=Shtr->nned; break; }
01317 case shellqelem:{ nned=Shq->nned; break; }
01318 case lineartet:{ nned=Ltet->nned; break; }
01319 case quadrtet:{ nned=Qtet->nned; break; }
01320 case linearhex:{ nned=Lhex->nned; break; }
01321 case quadrhex:{ nned=Qhex->nned; break; }
01322 case lineartetrot:{ nned=Ltetrot->nned; break; }
01323 case linearhexrot:{ nned=Lhexrot->nned; break; }
01324 case linearwed:{ nned=Lwed->nned; break; }
01325 case quadrwed:{ nned=Qwed->nned; break; }
01326 case particleelem:{ break; }
01327 default:{
01328 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01329 }
01330 }
01331 return nned;
01332 }
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345 long mechtop::give_nsurf (long eid)
01346 {
01347 long nsurf;
01348 elemtype te;
01349
01350 te = give_elem_type (eid);
01351
01352 switch (te){
01353 case dktel:{ nsurf=Dkt->nsurf; break; }
01354 case lineartet:{ nsurf=Ltet->nsurf; break; }
01355 case quadrtet:{ nsurf=Qtet->nsurf; break; }
01356 case linearhex:{ nsurf=Lhex->nsurf; break; }
01357 case quadrhex:{ nsurf=Qhex->nsurf; break; }
01358 case lineartetrot:{ nsurf=Ltetrot->nsurf; break; }
01359 case linearhexrot:{ nsurf=Lhexrot->nsurf; break; }
01360 case linearwed:{ nsurf=Lwed->nsurf; break; }
01361 case quadrwed:{ nsurf=Qwed->nsurf; break; }
01362 default:{
01363 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01364 }
01365 }
01366 return nsurf;
01367 }
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380 long mechtop::give_nnsurf (long eid)
01381 {
01382 long nnsurf;
01383 elemtype te;
01384
01385 te = give_elem_type (eid);
01386
01387 switch (te){
01388 case dktel:{ nnsurf=Dkt->nnsurf; break; }
01389 case lineartet:{ nnsurf=Ltet->nnsurf; break; }
01390 case quadrtet:{ nnsurf=Qtet->nnsurf; break; }
01391 case linearhex:{ nnsurf=Lhex->nnsurf; break; }
01392 case quadrhex:{ nnsurf=Qhex->nnsurf; break; }
01393 case lineartetrot:{ nnsurf=Ltetrot->nnsurf; break; }
01394 case linearhexrot:{ nnsurf=Lhexrot->nnsurf; break; }
01395 case linearwed:{ nnsurf=Lwed->nnsurf; break; }
01396 case quadrwed:{ nnsurf=Qwed->nnsurf; break; }
01397 default:{
01398 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01399 }
01400 }
01401 return nnsurf;
01402 }
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415 long mechtop::give_napfun (long eid)
01416 {
01417 long napfun=0;
01418 elemtype te;
01419
01420 te = give_elem_type (eid);
01421
01422 switch (te){
01423 case bar2d:{ napfun=Bar2d->napfun; break; }
01424 case bar3d:{ napfun=Bar3d->napfun; break; }
01425 case barq2d:{ napfun=Barq2d->napfun; break; }
01426 case barq3d:{ napfun=Barq3d->napfun; break; }
01427 case beam2d:{ napfun=Beam2d->napfun; break; }
01428 case beam3d:{ napfun=Beam3d->napfun; break; }
01429 case beamg3d:{ napfun=Beam3dg->napfun; break; }
01430 case subsoilbeam:{ napfun=Sbeam->napfun; break; }
01431 case beam2dsp:{ napfun=Spbeam2d->napfun; break; }
01432 case spring_1:{ napfun=Spring->napfun; break; }
01433 case spring_2:{ napfun=Spring->napfun; break; }
01434 case spring_3:{ napfun=Spring->napfun; break; }
01435 case spring_4:{ napfun=Spring->napfun; break; }
01436 case spring_5:{ napfun=Spring->napfun; break; }
01437 case spring_6:{ napfun=Spring->napfun; break; }
01438 case planeelementlt:{ napfun=Pelt->napfun; break; }
01439 case planeelementqt:{ napfun=Peqt->napfun; break; }
01440 case planeelementrotlt:{ napfun=Perlt->napfun; break; }
01441 case planeelementlq:{ napfun=Pelq->napfun; break; }
01442 case planeelementqq:{ napfun=Peqq->napfun; break; }
01443 case planeelementrotlq:{ napfun=Perlq->napfun; break; }
01444 case planeelementsubqt:{ napfun=Pesqt->napfun; break; }
01445 case cctel:{ napfun=Cct->napfun; break; }
01446 case dktel:{ napfun=Dkt->napfun; break; }
01447 case dstel:{ napfun=Dst->napfun; break; }
01448 case q4plateel:{ napfun=Q4pl->napfun; break; }
01449
01450 case argyristr:{ napfun=Argtrpl->napfun; break; }
01451 case subsoilplatetr:{ napfun=Spltr->napfun; break; }
01452 case subsoilplateq:{ napfun=Splq->napfun; break; }
01453 case axisymmlt:{ napfun=Asymlt->napfun; break; }
01454 case axisymmlq:{ napfun=Asymlq->napfun; break; }
01455 case axisymmqq:{ napfun=Asymqq->napfun; break; }
01456 case shelltrelem:{ napfun=Shtr->napfun; break; }
01457 case shellqelem:{ napfun=Shq->napfun; break; }
01458 case lineartet:{ napfun=Ltet->napfun; break; }
01459 case quadrtet:{ napfun=Qtet->napfun; break; }
01460 case linearhex:{ napfun=Lhex->napfun; break; }
01461 case quadrhex:{ napfun=Qhex->napfun; break; }
01462 case lineartetrot:{ napfun=Ltetrot->napfun; break; }
01463 case linearhexrot:{ napfun=Lhexrot->napfun; break; }
01464 case linearwed:{ napfun=Lwed->napfun; break; }
01465 case quadrwed:{ napfun=Qwed->napfun; break; }
01466 case particleelem:{ break; }
01467 default:{
01468 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01469 }
01470 }
01471 return napfun;
01472 }
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487 long mechtop::give_dimension (long eid)
01488 {
01489 long dim=0;
01490 elemtype te;
01491
01492 te = give_elem_type (eid);
01493
01494 switch (te){
01495 case bar2d:{ dim=1; break; }
01496 case bar3d:{ dim=1; break; }
01497 case barq2d:{ dim=1; break; }
01498 case barq3d:{ dim=1; break; }
01499 case beam2d:{ dim=1; break; }
01500 case beam3d:{ dim=1; break; }
01501 case beamg3d:{ dim=1; break; }
01502 case subsoilbeam:{ dim=1; break; }
01503 case beam2dsp:{ dim=1; break; }
01504 case spring_1:{ dim=1; break; }
01505 case spring_2:{ dim=1; break; }
01506 case spring_3:{ dim=1; break; }
01507 case spring_4:{ dim=1; break; }
01508 case spring_5:{ dim=1; break; }
01509 case spring_6:{ dim=1; break; }
01510 case planeelementlt:{ dim=2; break; }
01511 case planeelementqt:{ dim=2; break; }
01512 case planeelementrotlt:{ dim=2; break; }
01513 case planeelementlq:{ dim=2; break; }
01514 case planeelementqq:{ dim=2; break; }
01515 case planeelementrotlq:{ dim=2; break; }
01516 case planeelementsubqt:{ dim=2; break; }
01517 case cctel:{ dim=2; break; }
01518 case dktel:{ dim=2; break; }
01519 case dstel:{ dim=2; break; }
01520 case q4plateel:{ dim=2; break; }
01521 case argyristr:{ dim=2; break; }
01522 case axisymmlt:{ dim=2; break; }
01523 case axisymmlq:{ dim=2; break; }
01524 case axisymmqq:{ dim=2; break; }
01525 case subsoilplatetr:{ dim=2; break; }
01526 case subsoilplateq:{ dim=2; break; }
01527 case shelltrelem:{ dim=2; break; }
01528 case shellqelem:{ dim=2; break; }
01529 case lineartet:{ dim=3; break; }
01530 case quadrtet:{ dim=3; break; }
01531 case linearhex:{ dim=3; break; }
01532 case quadrhex:{ dim=3; break; }
01533 case lineartetrot:{ dim=3; break; }
01534 case linearhexrot:{ dim=3; break; }
01535 case linearwed:{ dim=3; break; }
01536 case quadrwed:{ dim=3; break; }
01537 case particleelem:{ break; }
01538 default:{
01539 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01540 }
01541 }
01542 return dim;
01543 }
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556 long mechtop::give_nb (long eid)
01557 {
01558 long nb=0;
01559 elemtype te;
01560
01561 te = give_elem_type (eid);
01562
01563 switch (te){
01564
01565 case bar2d:{ nb=Bar2d->nb; break; }
01566 case bar3d:{ nb=Bar3d->nb; break; }
01567 case barq2d:{ nb=Barq2d->nb; break; }
01568 case barq3d:{ nb=Barq3d->nb; break; }
01569 case beam2d:{ nb=Beam2d->nb; break; }
01570 case beam3d:{ nb=Beam3d->nb; break; }
01571 case beamg3d:{ nb=Beam3dg->nb; break; }
01572 case subsoilbeam:{ nb=Sbeam->nb; break; }
01573 case beam2dsp:{ nb=Spbeam2d->nb; break; }
01574 case spring_1:{ nb=Spring->nb; break; }
01575 case spring_2:{ nb=Spring->nb; break; }
01576 case spring_3:{ nb=Spring->nb; break; }
01577 case spring_4:{ nb=Spring->nb; break; }
01578 case spring_5:{ nb=Spring->nb; break; }
01579 case spring_6:{ nb=Spring->nb; break; }
01580 case planeelementlt:{ nb=Pelt->nb; break; }
01581 case planeelementqt:{ nb=Peqt->nb; break; }
01582 case planeelementrotlt:{ nb=Perlt->nb; break; }
01583 case planeelementlq:{ nb=Pelq->nb; break; }
01584 case planeelementqq:{ nb=Peqq->nb; break; }
01585 case planeelementrotlq:{ nb=Perlq->nb; break; }
01586 case planeelementsubqt:{ nb=Pesqt->nb; break; }
01587 case planequadcontact:{ nb=Pqcon->nb; break; }
01588 case cctel:{ nb=Cct->nb; break; }
01589 case dktel:{ nb=Dkt->nb; break; }
01590 case dstel:{ nb=Dst->nb; break; }
01591 case q4plateel:{ nb=Q4pl->nb; break; }
01592
01593 case argyristr:{ nb=Argtrpl->nb; break; }
01594 case subsoilplatetr:{ nb=Spltr->nb; break; }
01595 case subsoilplateq:{ nb=Splq->nb; break; }
01596 case axisymmlt:{ nb=Asymlt->nb; break; }
01597 case axisymmlq:{ nb=Asymlq->nb; break; }
01598 case axisymmqq:{ nb=Asymqq->nb; break; }
01599 case shelltrelem:{ nb=Shtr->nb; break; }
01600 case shellqelem:{ nb=Shq->nb; break; }
01601 case lineartet:{ nb=Ltet->nb; break; }
01602 case quadrtet:{ nb=Qtet->nb; break; }
01603 case linearhex:{ nb=Lhex->nb; break; }
01604 case quadrhex:{ nb=Qhex->nb; break; }
01605 case lineartetrot:{ nb=Ltetrot->nb; break; }
01606 case linearhexrot:{ nb=Lhexrot->nb; break; }
01607 case linearwed:{ nb=Lwed->nb; break; }
01608 case quadrwed:{ nb=Qwed->nb; break; }
01609 case particleelem:{ nb=Pelem->nb; break; }
01610 default:{
01611 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01612 }
01613 }
01614 return nb;
01615 }
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627 long mechtop::give_nb_te (elemtype te)
01628 {
01629 long nb=0;
01630
01631 switch (te){
01632
01633 case bar2d:{ nb=Bar2d->nb; break; }
01634 case bar3d:{ nb=Bar3d->nb; break; }
01635 case barq2d:{ nb=Barq2d->nb; break; }
01636 case barq3d:{ nb=Barq3d->nb; break; }
01637 case beam2d:{ nb=Beam2d->nb; break; }
01638 case beam3d:{ nb=Beam3d->nb; break; }
01639 case beamg3d:{ nb=Beam3dg->nb; break; }
01640 case subsoilbeam:{ nb=Sbeam->nb; break; }
01641 case beam2dsp:{ nb=Spbeam2d->nb; break; }
01642 case spring_1:{ nb=Spring->nb; break; }
01643 case spring_2:{ nb=Spring->nb; break; }
01644 case spring_3:{ nb=Spring->nb; break; }
01645 case spring_4:{ nb=Spring->nb; break; }
01646 case spring_5:{ nb=Spring->nb; break; }
01647 case spring_6:{ nb=Spring->nb; break; }
01648 case planeelementlt:{ nb=Pelt->nb; break; }
01649 case planeelementqt:{ nb=Peqt->nb; break; }
01650 case planeelementrotlt:{ nb=Perlt->nb; break; }
01651 case planeelementlq:{ nb=Pelq->nb; break; }
01652 case planeelementqq:{ nb=Peqq->nb; break; }
01653 case planeelementrotlq:{ nb=Perlq->nb; break; }
01654 case planeelementsubqt:{ nb=Pesqt->nb; break; }
01655 case planequadcontact:{ nb=Pqcon->nb; break; }
01656 case cctel:{ nb=Cct->nb; break; }
01657 case dktel:{ nb=Dkt->nb; break; }
01658 case dstel:{ nb=Dst->nb; break; }
01659 case q4plateel:{ nb=Q4pl->nb; break; }
01660
01661 case argyristr:{ nb=Argtrpl->nb; break; }
01662 case subsoilplatetr:{ nb=Spltr->nb; break; }
01663 case subsoilplateq:{ nb=Splq->nb; break; }
01664 case axisymmlt:{ nb=Asymlt->nb; break; }
01665 case axisymmlq:{ nb=Asymlq->nb; break; }
01666 case axisymmqq:{ nb=Asymqq->nb; break; }
01667 case shelltrelem:{ nb=Shtr->nb; break; }
01668 case shellqelem:{ nb=Shq->nb; break; }
01669 case lineartet:{ nb=Ltet->nb; break; }
01670 case quadrtet:{ nb=Qtet->nb; break; }
01671 case linearhex:{ nb=Lhex->nb; break; }
01672 case quadrhex:{ nb=Qhex->nb; break; }
01673 case lineartetrot:{ nb=Ltetrot->nb; break; }
01674 case linearhexrot:{ nb=Lhexrot->nb; break; }
01675 case linearwed:{ nb=Lwed->nb; break; }
01676 case quadrwed:{ nb=Qwed->nb; break; }
01677 case particleelem:{ nb=Pelem->nb; break; }
01678 default:{
01679 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01680 }
01681 }
01682 return nb;
01683 }
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723 strastrestate mechtop::give_ssst (long eid,long bi)
01724 {
01725 strastrestate ssst;
01726 elemtype te;
01727
01728 te = give_elem_type (eid);
01729 switch (te){
01730 case bar2d:{ ssst = Bar2d->ssst; break; }
01731 case bar3d:{ ssst = Bar3d->ssst; break; }
01732 case barq2d:{ ssst = Barq2d->ssst; break; }
01733 case barq3d:{ ssst = Barq3d->ssst; break; }
01734 case beam2d:{ ssst = Beam2d->ssst; break; }
01735 case beam3d:{ ssst = Beam3d->ssst; break; }
01736 case beamg3d:{ ssst = Beam3dg->ssst; break; }
01737 case subsoilbeam:{ ssst = Sbeam->ssst; break; }
01738 case beam2dsp:{ ssst = Spbeam2d->ssst; break; }
01739 case spring_1:{ ssst = Spring->ssst; break; }
01740 case spring_2:{ ssst = Spring->ssst; break; }
01741 case spring_3:{ ssst = Spring->ssst; break; }
01742 case spring_4:{ ssst = Spring->ssst; break; }
01743 case spring_5:{ ssst = Spring->ssst; break; }
01744 case spring_6:{ ssst = Spring->ssst; break; }
01745 case planeelementlt:{ ssst = elements[eid].ssst; break; }
01746 case planeelementqt:{ ssst = elements[eid].ssst; break; }
01747 case planeelementrotlt:{ ssst = elements[eid].ssst; break; }
01748 case planeelementlq:{ ssst = elements[eid].ssst; break; }
01749 case planeelementqq:{ ssst = elements[eid].ssst; break; }
01750 case planeelementrotlq:{ ssst = elements[eid].ssst; break; }
01751 case planeelementsubqt:{ ssst = elements[eid].ssst; break; }
01752 case planequadcontact:{ ssst = Pqcon->ssst; break; }
01753 case cctel:{ ssst = Cct->ssst; break; }
01754 case dktel:{ ssst = Dkt->ssst; break; }
01755 case dstel:{ ssst = Dst->ssst; break; }
01756 case q4plateel:{ ssst = Q4pl->ssst; break; }
01757
01758 case argyristr:{ ssst = Argtrpl->ssst; break; }
01759 case subsoilplatetr:{ ssst = Spltr->ssst; break; }
01760 case subsoilplateq:{ ssst = Splq->ssst; break; }
01761 case axisymmlt:{ ssst = Asymlt->ssst; break; }
01762 case axisymmlq:{ ssst = Asymlq->ssst; break; }
01763 case axisymmqq:{ ssst = Asymqq->ssst; break; }
01764 case shelltrelem:{
01765 if (bi==0 || bi==1) ssst = (strastrestate) planestress;
01766 if (bi==2) ssst = (strastrestate) platek;
01767 if (bi==3) ssst = (strastrestate) shell;
01768 break; }
01769 case shellqelem:{
01770 if (bi==0 || bi==1) ssst = (strastrestate) planestress;
01771 if (bi==2 || bi==3) ssst = (strastrestate) plates;
01772 if (bi==4) ssst = (strastrestate) shell;
01773 break; }
01774 case lineartet:{ ssst = Ltet->ssst; break; }
01775 case quadrtet:{ ssst = Qtet->ssst; break; }
01776 case linearhex:{ ssst = Lhex->ssst; break; }
01777 case quadrhex:{ ssst = Qhex->ssst; break; }
01778 case lineartetrot:{ ssst = Ltetrot->ssst; break; }
01779 case linearhexrot:{ ssst = Lhexrot->ssst; break; }
01780 case linearwed:{ ssst = Lwed->ssst; break; }
01781 case quadrwed:{ ssst = Qwed->ssst; break; }
01782 case particleelem:{ break; }
01783 default:{
01784 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01785 }
01786 }
01787 return ssst;
01788 }
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801 long mechtop::give_degree (long eid)
01802 {
01803 long deg=0;
01804 elemtype te;
01805
01806 te = give_elem_type (eid);
01807
01808 switch (te){
01809 case bar2d:{ deg=1; break; }
01810 case bar3d:{ deg=1; break; }
01811 case barq2d:{ deg=2; break; }
01812 case barq3d:{ deg=2; break; }
01813 case beam2d:{ deg=1; break; }
01814 case beam3d:{ deg=1; break; }
01815 case beamg3d:{ deg=3; break; }
01816 case spring_1:{ deg=1; break; }
01817 case spring_2:{ deg=1; break; }
01818 case spring_3:{ deg=1; break; }
01819 case spring_4:{ deg=1; break; }
01820 case spring_5:{ deg=1; break; }
01821 case spring_6:{ deg=1; break; }
01822 case subsoilbeam:{ deg=1; break; }
01823 case beam2dsp:{ deg=1; break; }
01824 case planeelementlt:{ deg=1; break; }
01825 case planeelementqt:{ deg=2; break; }
01826 case planeelementrotlt:{ deg=1; break; }
01827 case planeelementlq:{ deg=1; break; }
01828 case planeelementqq:{ deg=2; break; }
01829 case planeelementrotlq:{ deg=1; break; }
01830 case planeelementsubqt:{ deg=2; break; }
01831 case cctel:{ deg=1; break; }
01832 case dktel:{ deg=1; break; }
01833 case dstel:{ deg=1; break; }
01834 case q4plateel:{ deg=1; break; }
01835 case argyristr:{ deg=3; break; }
01836 case subsoilplatetr:{ deg=1; break; }
01837 case subsoilplateq:{ deg=1; break; }
01838 case axisymmlt:{ deg=1; break; }
01839 case axisymmlq:{ deg=1; break; }
01840 case axisymmqq:{ deg=1; break; }
01841 case shelltrelem:{ deg=1; break; }
01842 case shellqelem:{ deg=1; break; }
01843 case lineartet:{ deg=1; break; }
01844 case quadrtet:{ deg=2; break; }
01845 case linearhex:{ deg=1; break; }
01846 case quadrhex:{ deg=2; break; }
01847 case lineartetrot:{ deg=1; break; }
01848 case linearhexrot:{ deg=1; break; }
01849 case linearwed:{ deg=1; break; }
01850 case quadrwed:{ deg=2; break; }
01851 case particleelem:{ break; }
01852 default:{
01853 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01854 }
01855 }
01856 return deg;
01857 }
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872 void mechtop::give_edge_nodes (long eid,long edid,long *nodes)
01873 {
01874 long nne;
01875 elemtype te;
01876 ivector nod;
01877
01878 te = give_elem_type (eid);
01879 nne = give_nne (eid);
01880 allocv (nne,nod);
01881 give_elemnodes (eid,nod);
01882
01883 switch (te){
01884 case planeelementlt:{
01885 if (edid==0){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01886 if (edid==1){ nodes[0]=nod[1]; nodes[1]=nod[2]; }
01887 if (edid==2){ nodes[0]=nod[2]; nodes[1]=nod[0]; }
01888 break;
01889 }
01890 case planeelementqt:{
01891 if (edid==0){ nodes[0]=nod[0]; nodes[1]=nod[1]; nodes[2]=nod[3]; }
01892 if (edid==1){ nodes[0]=nod[1]; nodes[1]=nod[2]; nodes[2]=nod[4]; }
01893 if (edid==2){ nodes[0]=nod[2]; nodes[1]=nod[0]; nodes[2]=nod[5]; }
01894 break;
01895 }
01896 case planeelementrotlt:{
01897 if (edid==0){ nodes[0]=nod[1]; nodes[1]=nod[2]; }
01898 if (edid==1){ nodes[0]=nod[2]; nodes[1]=nod[0]; }
01899 if (edid==2){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01900 break;
01901 }
01902 case planeelementlq:{
01903 if (edid==0){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01904 if (edid==1){ nodes[0]=nod[1]; nodes[1]=nod[2]; }
01905 if (edid==2){ nodes[0]=nod[2]; nodes[1]=nod[3]; }
01906 if (edid==3){ nodes[0]=nod[3]; nodes[1]=nod[0]; }
01907 break;
01908 }
01909 case planeelementqq:{
01910 if (edid==0){ nodes[0]=nod[0]; nodes[1]=nod[4]; nodes[2]=nod[1]; }
01911 if (edid==1){ nodes[0]=nod[1]; nodes[1]=nod[5]; nodes[2]=nod[2]; }
01912 if (edid==2){ nodes[0]=nod[2]; nodes[1]=nod[6]; nodes[2]=nod[3]; }
01913 if (edid==3){ nodes[0]=nod[3]; nodes[1]=nod[7]; nodes[2]=nod[0]; }
01914 break;
01915 }
01916 case planeelementrotlq:{
01917 if (edid==0){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01918 if (edid==1){ nodes[0]=nod[1]; nodes[1]=nod[2]; }
01919 if (edid==2){ nodes[0]=nod[2]; nodes[1]=nod[3]; }
01920 if (edid==3){ nodes[0]=nod[3]; nodes[1]=nod[0]; }
01921 break;
01922 }
01923 case planeelementsubqt:{
01924 if (edid==0){ nodes[0]=nod[0]; nodes[1]=nod[1]; nodes[2]=nod[3]; }
01925 if (edid==1){ nodes[0]=nod[1]; nodes[1]=nod[2]; nodes[2]=nod[4]; }
01926 if (edid==2){ nodes[0]=nod[2]; nodes[1]=nod[0]; nodes[2]=nod[5]; }
01927 break;
01928 }
01929 case axisymmlt:{
01930 if (edid==0){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01931 if (edid==1){ nodes[0]=nod[1]; nodes[1]=nod[2]; }
01932 if (edid==2){ nodes[0]=nod[2]; nodes[1]=nod[0]; }
01933 break;
01934 }
01935 case linearhex:{
01936 if (edid==0){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01937 if (edid==1){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01938 if (edid==2){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01939 if (edid==3){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01940 if (edid==4){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01941 if (edid==5){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01942 if (edid==6){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01943 if (edid==7){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01944 if (edid==8){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01945 if (edid==9){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01946 if (edid==10){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01947 if (edid==11){ nodes[0]=nod[0]; nodes[1]=nod[1]; }
01948 break;
01949 }
01950 default:{
01951 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01952 }
01953 }
01954 destrv (nod);
01955 }
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980 void mechtop::give_elemnodes (long eid,ivector &nodes)
01981 {
01982 Gtm->give_original_nodes (eid,nodes);
01983 }
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997 void mechtop::give_code_numbers (long eid,long *cn)
01998 {
01999 Gtm->give_code_numbers (eid,cn);
02000 }
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014 void mechtop::give_node_code_numbers (long nid,long *cn)
02015 {
02016 Gtm->give_node_code_numbers (nid,cn);
02017 }
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032 void mechtop::give_mult_code_numbers (long nid,long lid,long *cn)
02033 {
02034 Gtm->give_mult_code_numbers (nid,lid,cn);
02035 }
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049 void mechtop::give_node_coord1d (vector &x,long eid)
02050 {
02051 Gtm->give_node_coord1d (x,eid);
02052 }
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066 void mechtop::give_node_coord2d (vector &x,vector &y,long eid)
02067 {
02068 Gtm->give_node_coord2d (x,y,eid);
02069 }
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083 void mechtop::give_node_coord2dxz (vector &x,vector &z,long eid)
02084 {
02085 Gtm->give_node_coord2dxz (x,z,eid);
02086 }
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098
02099
02100 void mechtop::give_node_coord3d (vector &x,vector &y,vector &z,long eid)
02101 {
02102 Gtm->give_node_coord3d (x,y,z,eid);
02103 }
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119 void mechtop::give_ipcoord_elem(long eid, matrix &coord)
02120 {
02121 long ii, jj, ipp, i, k, l;
02122 long nb = give_nb(eid);
02123 vector aux(3);
02124
02125 k = 0;
02126 for(ii=0; ii<nb; ii++)
02127 {
02128 for(jj=0; jj<nb; jj++)
02129 {
02130 ipp=elements[eid].ipp[ii][jj];
02131 for (i=0; i<give_nip(eid, ii, jj); i++)
02132 {
02133 fillv(0.0, aux);
02134 ipcoord(eid, ipp+i, ii, jj, aux);
02135 for (l=0; l<3; l++)
02136 coord[k][l] = aux[l];
02137 k++;
02138 }
02139
02140 }
02141 }
02142 }
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155 double mechtop::give_length (long eid)
02156 {
02157 long nn = give_nne(eid);
02158 double l;
02159 if (nn > 2)
02160 print_err("invalid number of nodes on element is detected", __FILE__, __LINE__, __func__);
02161
02162 vector x(nn), y(nn), z(nn);
02163 Gtm->give_node_coord3d (x, y, z, eid);
02164 l = sqrt(sqr(x[1]-x[0]) + sqr(y[1]-y[0]) + sqr(z[1]-z[0]));
02165 return l;
02166 }
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179 double mechtop::give_area (long eid)
02180 {
02181 long nn = give_nne(eid);
02182 double a = 0.0;
02183 if (nn < 2)
02184 print_err("invalid number of nodes on element is detected", __FILE__, __LINE__, __func__);
02185
02186 vector x(nn), y(nn), z(nn);
02187 Gtm->give_node_coord3d (x, y, z, eid);
02188 switch (Mt->give_elem_type(eid)){
02189 case planeelementlt:
02190 case planeelementqt:
02191 case planeelementrotlt:
02192 case axisymmlt:
02193 a = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
02194 break;
02195 case planeelementlq:
02196 case planeelementqq:
02197 case planeelementrotlq:
02198 case axisymmlq:
02199 case axisymmqq:
02200 a = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]) + (x[2]-x[0])*(y[3]-y[0])-(x[3]-x[0])*(y[2]-y[0]) ) / 2.0;
02201 break;
02202 case cctel:
02203 case dktel:
02204 case dstel:
02205 a = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
02206 break;
02207 case q4plateel:
02208 a = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]) + (x[2]-x[0])*(y[3]-y[0])-(x[3]-x[0])*(y[2]-y[0]) ) / 2.0;
02209 break;
02210 default:
02211 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
02212 }
02213 return a;
02214 }
02215
02216
02217
02218
02219
02220
02221
02222
02223
02224
02225
02226
02227 double mechtop::give_volume (long eid)
02228 {
02229 long nn = give_nne(eid);
02230 double v = 0.0;
02231 if (nn < 2)
02232 print_err("invalid number of nodes on element is detected", __FILE__, __LINE__, __func__);
02233
02234 vector x(nn), y(nn), z(nn);
02235 Gtm->give_node_coord3d (x, y, z, eid);
02236 switch (Mt->give_elem_type(eid)){
02237 case lineartet:
02238 case quadrtet:
02239 jac_3d (v, x, y, z, 0.0, 0.0, 0.0);
02240 v/=6.0;
02241 break;
02242 case linearhex:
02243 case quadrhex:
02244 case lineartetrot:
02245 case linearhexrot:
02246 jac_3d (v, x, y, z, 1.0, 1.0, 1.0);
02247 v*=8.0;
02248 break;
02249 default:
02250 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
02251 }
02252
02253 if(v < 0.0){
02254 print_err("wrong numbering of nodes on 3D element number %ld, negative volume! v = %e",__FILE__,__LINE__,__func__,eid,v);
02255 abort();
02256 }
02257
02258 return v;
02259 }
02260
02261
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276 void mechtop::give_maxncompstr (long &max_ncompstrn, long &max_ncompstre)
02277 {
02278 long i;
02279 max_ncompstrn = max_ncompstre = 0;
02280
02281 for (i=0;i<ne;i++)
02282 {
02283 if (max_ncompstre < Mm->ip[elements[i].ipp[0][0]].ncompstr)
02284 max_ncompstre = Mm->ip[elements[i].ipp[0][0]].ncompstr;
02285 }
02286 for (i=0;i<nn;i++)
02287 {
02288 if (max_ncompstrn < nodes[i].ncompstr)
02289 max_ncompstrn = nodes[i].ncompstr;
02290 }
02291 }
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308 void mechtop::give_maxncompo (long &max_nncompo, long &max_encompo)
02309 {
02310 long i;
02311 max_nncompo = max_encompo = 0;
02312
02313 for (i=0;i<ne;i++)
02314 {
02315 if (max_encompo < Mm->ip[elements[i].ipp[0][0]].ncompeqother)
02316 max_encompo = Mm->ip[elements[i].ipp[0][0]].ncompeqother;
02317 }
02318 for (i=0;i<nn;i++)
02319 {
02320 if (max_nncompo < nodes[i].ncompother)
02321 max_nncompo = nodes[i].ncompother;
02322 }
02323 }
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337 long mechtop::give_maxndofn()
02338 {
02339 long i, ndofn;
02340 if (max_ndofn < 0)
02341 {
02342 for(i=0; i<nn; i++)
02343 {
02344 ndofn = give_ndofn(i);
02345 if (ndofn > max_ndofn)
02346 max_ndofn = ndofn;
02347 }
02348 }
02349 return max_ndofn;
02350 }
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363 long mechtop::locsystems (ivector &nod)
02364 {
02365 long i,j,n;
02366 n=nod.n;
02367 j=0;
02368 for (i=0;i<n;i++){
02369 j+=nodes[nod[i]].transf;
02370 }
02371 return j;
02372 }
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383 void mechtop::comreacnod (void)
02384 {
02385 long i,j,k,ndofn,nmn,nid;
02386
02387
02388 for (i=0;i<nn;i++){
02389
02390
02391 ndofn = Gtm->give_ndofn (i);
02392
02393 if (ndofn>0){
02394
02395
02396
02397 for (j=0;j<ndofn;j++){
02398 if (Gtm->give_dof (i,j)<1){
02399 nodes[i].react=1;
02400 if (nodes[i].r==NULL)
02401 nodes[i].r = new double [ndofn];
02402 memset (nodes[i].r,0,ndofn*sizeof(double));
02403 break;
02404 }
02405 }
02406 }
02407 if (ndofn<0){
02408
02409
02410
02411 nmn=0-ndofn;
02412
02413
02414 for (j=0;j<nmn;j++){
02415
02416 nid = Gtm->gnodes[i].mnodes[j];
02417
02418 ndofn = Gtm->give_ndofn (nid);
02419
02420
02421 for (k=0;k<ndofn;k++){
02422 if (Gtm->give_dof (nid,k)<1){
02423 nodes[i].react=1;
02424 if (nodes[i].r==NULL)
02425 nodes[i].r = new double [ndofn];
02426 memset (nodes[i].r,0,ndofn*sizeof(double));
02427 break;
02428 }
02429 }
02430
02431 }
02432
02433 }
02434
02435 }
02436 }
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447 void mechtop::comreacelem (void)
02448 {
02449 long i,j,k,nne;
02450 ivector nod;
02451
02452
02453 for (i=0;i<ne;i++){
02454 if (Gtm->leso[i]==1){
02455
02456
02457 nne=give_nne (i);
02458 reallocv (nne,nod);
02459
02460 give_elemnodes (i,nod);
02461
02462
02463 for (j=0;j<nne;j++){
02464 k=nod[j];
02465 if (nodes[k].react==1){
02466 elements[i].react=1;
02467 break;
02468 }
02469 }
02470 }
02471 }
02472 }
02473
02474
02475
02476
02477
02478
02479
02480
02481
02482
02483
02484 void mechtop::comreac (void)
02485 {
02486 comreacnod ();
02487 comreacelem ();
02488 }
02489
02490
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507 void mechtop::strain_nodal_values (ivector &nod,vector &nx,vector &ny,vector &nz,
02508 double *lhs,long dim,long fi,long ncomp,long lcid)
02509 {
02510 long i,j;
02511 vector eps;
02512
02513 allocv (ncomp,eps);
02514
02515 if (dim==1){
02516 for (i=0;i<nod.n;i++){
02517 for (j=0;j<ncomp;j++){
02518 eps[j] = lhs[j*2]+lhs[j*2+1]*nx[i];
02519 }
02520 nodes[nod[i]].storestrain (lcid,fi,ncomp,eps);
02521 }
02522 }
02523 if (dim==2){
02524 for (i=0;i<nod.n;i++){
02525 for (j=0;j<ncomp;j++){
02526 eps[j] = lhs[j*3]+lhs[j*3+1]*nx[i]+lhs[j*3+2]*ny[i];
02527 }
02528 nodes[nod[i]].storestrain (lcid,fi,ncomp,eps);
02529 }
02530 }
02531 if (dim==3){
02532 for (i=0;i<nod.n;i++){
02533 for (j=0;j<ncomp;j++){
02534 eps[j] = lhs[j*4]+lhs[j*4+1]*nx[i]+lhs[j*4+2]*ny[i]+lhs[j*4+3]*nz[i];
02535 }
02536 nodes[nod[i]].storestrain (lcid,fi,ncomp,eps);
02537 }
02538 }
02539
02540 destrv (eps);
02541 }
02542
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560 void mechtop::stress_nodal_values (ivector &nod,vector &nx,vector &ny,vector &nz,
02561 double *lhs,long dim,long fi,long ncomp,long lcid)
02562 {
02563 long i,j;
02564 vector sig;
02565
02566 allocv (ncomp,sig);
02567 if (dim==2){
02568 for (i=0;i<nod.n;i++){
02569 for (j=0;j<ncomp;j++){
02570
02571 sig[j] = lhs[j*3]+lhs[j*3+1]*nx[i]+lhs[j*3+2]*ny[i];
02572 }
02573 nodes[nod[i]].storestress (lcid,fi,ncomp,sig);
02574
02575 }
02576 }
02577 if (dim==3){
02578 for (i=0;i<nod.n;i++){
02579 for (j=0;j<ncomp;j++){
02580
02581 sig[j] = lhs[j*4]+lhs[j*4+1]*nx[i]+lhs[j*4+2]*ny[i]+lhs[j*4+3]*nz[i];
02582 }
02583 nodes[nod[i]].storestress (lcid,fi,ncomp,sig);
02584 }
02585 }
02586 destrv (sig);
02587 }
02588
02589
02590
02591
02592
02593
02594
02595
02596
02597
02598
02599
02600
02601
02602
02603
02604
02605 void mechtop::other_nodal_values (ivector &nod,vector &nx,vector &ny,vector &nz,
02606 double *lhs,long dim,long fi,long ncomp)
02607 {
02608 long i,j;
02609 vector other;
02610
02611 allocv (ncomp,other);
02612 if (dim==1){
02613 for (i=0;i<nod.n;i++){
02614 for (j=0;j<ncomp;j++){
02615 other[j] = lhs[j*2]+lhs[j*2+1]*nx[i];
02616 }
02617 nodes[nod[i]].storeother (fi,ncomp,other);
02618 }
02619 }
02620 if (dim==2){
02621 for (i=0;i<nod.n;i++){
02622 for (j=0;j<ncomp;j++){
02623 other[j] = lhs[j*3]+lhs[j*3+1]*nx[i]+lhs[j*3+2]*ny[i];
02624 }
02625 nodes[nod[i]].storeother (fi,ncomp,other);
02626
02627 }
02628 }
02629 if (dim==3){
02630 for (i=0;i<nod.n;i++){
02631 for (j=0;j<ncomp;j++){
02632 other[j] = lhs[j*4]+lhs[j*4+1]*nx[i]+lhs[j*4+2]*ny[i]+lhs[j*4+3]*nz[i];
02633 }
02634 nodes[nod[i]].storeother (fi,ncomp,other);
02635 }
02636 }
02637 destrv (other);
02638 }
02639
02640
02641
02642
02643
02644
02645
02646
02647
02648
02649
02650
02651
02652 void mechtop::give_nodal_strain(long nid, vector &stra)
02653 {
02654 long i, j, k;
02655 matrix a(3,3);
02656 strastrestate ssst;
02657
02658 ssst = guess_ssst(nodes[nid].ncompstr);
02659 for (i=0; i<nodes[nid].ncompstr; i++)
02660 stra(i) = nodes[nid].strain[i];
02661
02662 vector_tensor(stra, a, ssst, strain);
02663 k=0;
02664 for (i=0; i<3; i++)
02665 {
02666 for(j=0; j<3; j++)
02667 {
02668 stra[k] = a[i][j];
02669 k++;
02670 }
02671 }
02672 }
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686 void mechtop::give_nodal_stress(long nid, vector &stre)
02687 {
02688 long i, j, k;
02689 matrix a(3,3);
02690 strastrestate ssst;
02691
02692 ssst = guess_ssst(nodes[nid].ncompstr);
02693 for (i=0; i<nodes[nid].ncompstr; i++)
02694 stre(i) = nodes[nid].stress[i];
02695
02696 vector_tensor(stre, a, ssst, stress);
02697 k=0;
02698 for (i=0; i<3; i++)
02699 {
02700 for(j=0; j<3; j++)
02701 {
02702 stre[k] = a[i][j];
02703 k++;
02704 }
02705 }
02706 }
02707
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720 void mechtop::give_nodal_other(long nid, vector &othv)
02721 {
02722 long i;
02723 for (i=0; i<nodes[nid].ncompother; i++)
02724 {
02725 othv(i) = nodes[nid].other[i];
02726 }
02727 }
02728
02729
02730
02731
02732
02733
02734
02735
02736
02737
02738
02739
02740 void mechtop::store_code_num_elem(void)
02741 {
02742 long i, j, nne, ndofe, ndofn;
02743 elemtype te;
02744 ivector nod;
02745 long *cn;
02746 for(i=0; i < ne; i++)
02747 {
02748 te = give_elem_type(i);
02749 nne = give_nne(i);
02750 ndofe = give_ndofe(i);
02751 allocv(nne, nod);
02752 give_elemnodes(i, nod);
02753 switch(te)
02754 {
02755 case beam2d:
02756 break;
02757 case spring_1:
02758 {
02759 ndofn = give_ndofn(nod[0]);
02760 if (ndofn < 1)
02761 {
02762 print_err("spring_1 is connected to node %ld with wrong number of dofs(%ld)", __FILE__, __LINE__, __func__, nod[0], ndofn);
02763 abort();
02764 }
02765 if (Gtm->gelements[i].cne == 0)
02766 {
02767 Gtm->gelements[i].cne = 1;
02768 Gtm->gelements[i].cn = new long[ndofn];
02769 cn = new long[ndofn];
02770 memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
02771 Mt->give_node_code_numbers(nod[0], cn);
02772 Gtm->gelements[i].cn[0] = cn[0];
02773 delete [] cn;
02774 }
02775 break;
02776 }
02777 case spring_2:
02778 {
02779 ndofn = give_ndofn(nod[0]);
02780 if (ndofn < 2)
02781 {
02782 print_err("spring_2 is connected to node %ld with wrong number of dofs(%ld)", __FILE__, __LINE__, __func__, nod[0], ndofn);
02783 abort();
02784 }
02785 if (Gtm->gelements[i].cne == 0)
02786 {
02787 Gtm->gelements[i].cne = 1;
02788 Gtm->gelements[i].cn = new long[ndofn];
02789 cn = new long[ndofn];
02790 memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
02791 Mt->give_node_code_numbers(nod[0], cn);
02792 Gtm->gelements[i].cn[1] = cn[1];
02793 delete [] cn;
02794 }
02795 break;
02796 }
02797 case spring_3:
02798 {
02799 ndofn = give_ndofn(nod[0]);
02800 if (ndofn < 3)
02801 {
02802 print_err("spring_3 is connected to node %ld with wrong number of dofs(%ld)", __FILE__, __LINE__, __func__, nod[0], ndofn);
02803 abort();
02804 }
02805 if (Gtm->gelements[i].cne == 0)
02806 {
02807 Gtm->gelements[i].cne = 1;
02808 Gtm->gelements[i].cn = new long[ndofn];
02809 cn = new long[ndofn];
02810 memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
02811 Mt->give_node_code_numbers(nod[0], cn);
02812 Gtm->gelements[i].cn[2] = cn[2];
02813 delete [] cn;
02814 }
02815 break;
02816 }
02817 case spring_4:
02818 {
02819 ndofn = give_ndofn(nod[0]);
02820 if (ndofn < 1)
02821 {
02822 print_err("spring_4 is connected to node %ld with wrong number of dofs(%ld)", __FILE__, __LINE__, __func__, nod[0], ndofn);
02823 abort();
02824 }
02825 if (Gtm->gelements[i].cne == 0)
02826 {
02827 Gtm->gelements[i].cne = 1;
02828 Gtm->gelements[i].cn = new long[ndofn];
02829 cn = new long[ndofn];
02830 memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
02831 Mt->give_node_code_numbers(nod[0], cn);
02832 Gtm->gelements[i].cn[3] = cn[3];
02833 delete [] cn;
02834 }
02835 break;
02836 }
02837 case spring_5:
02838 {
02839 ndofn = give_ndofn(nod[0]);
02840 if (ndofn < 1)
02841 {
02842 print_err("spring_5 is connected to node %ld with wrong number of dofs(%ld)", __FILE__, __LINE__, __func__, nod[0], ndofn);
02843 abort();
02844 }
02845 if (Gtm->gelements[i].cne == 0)
02846 {
02847 Gtm->gelements[i].cne = 1;
02848 Gtm->gelements[i].cn = new long[ndofn];
02849 cn = new long[ndofn];
02850 memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
02851 Mt->give_node_code_numbers(nod[0], cn);
02852 Gtm->gelements[i].cn[4] = cn[4];
02853 delete [] cn;
02854 }
02855 break;
02856 }
02857 case spring_6:
02858 {
02859 if (ndofn < 1)
02860 {
02861 print_err("spring_6 is connected to node %ld with wrong number of dofs(%ld)", __FILE__, __LINE__, __func__, nod[0], ndofn);
02862 abort();
02863 }
02864 if (Gtm->gelements[i].cne == 0)
02865 {
02866 Gtm->gelements[i].cne = 1;
02867 Gtm->gelements[i].cn = new long[ndofn];
02868 cn = new long[ndofn];
02869 memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
02870 Mt->give_node_code_numbers(nod[0], cn);
02871 Gtm->gelements[i].cn[5] = cn[5];
02872 delete [] cn;
02873 }
02874 break;
02875 }
02876 case bar2d:
02877 case barq2d:
02878 case planeelementlt:
02879 case planeelementqt:
02880 case planeelementlq:
02881 case planeelementqq:
02882 if (Gtm->gelements[i].cne == 0)
02883 {
02884 Gtm->gelements[i].cne = 1;
02885 Gtm->gelements[i].cn = new long[ndofe];
02886 memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofe);
02887 for (j=0; j<nne; j++)
02888 {
02889 ndofn = Gtm->give_ndofn(nod[j]);
02890 cn = new long[ndofn];
02891 Mt->give_node_code_numbers(nod[j], cn);
02892 Gtm->gelements[i].cn[2*j] = cn[0];
02893 Gtm->gelements[i].cn[2*j+1] = cn[1];
02894 delete [] cn;
02895 }
02896 }
02897 break;
02898 default:
02899 print_err("unknown type of element is required", __FILE__, __LINE__, __func__);
02900 }
02901 }
02902 }
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914 void mechtop::alloc_nodes_arrays (void)
02915 {
02916 long i,j,nne,ncomp, ncompo;
02917 ivector nod;
02918 strastrestate ssst;
02919 elemtype te;
02920
02921 if (Mp->strainpos==2 || Mp->strainpos==3 || Mp->stresspos==2 || Mp->stresspos==3 || Mp->otherpos==2){
02922 for (i=0;i<ne;i++){
02923 if (Gtm->leso[i]==1){
02924
02925 nne = give_nne (i);
02926 reallocv (nne,nod);
02927 give_elemnodes (i,nod);
02928
02929 ncomp=give_tncomp (i);
02930
02931 ncompo=Mm->ip[elements[i].ipp[0][0]].ncompeqother;
02932
02933
02934 te = give_elem_type (i);
02935 if (te != shelltrelem && te != shellqelem){
02936
02937 ssst=Mt->give_ssst (i,0);
02938 if ((ssst == planestrain) || (ssst == planestress)){
02939 ncomp=4;
02940 }
02941 }
02942
02943
02944 for (j=0;j<nne;j++){
02945 if (ncomp > nodes[nod[j]].ncompstr){
02946
02947 if (Mp->strainpos > 1)
02948 nodes[nod[j]].alloc_strain(ncomp,Mb->nlc);
02949 if (Mp->stresspos > 1)
02950 nodes[nod[j]].alloc_stress(ncomp,Mb->nlc);
02951 }
02952 if (ncompo > nodes[nod[j]].ncompother){
02953
02954 if (Mp->otherpos > 1)
02955 nodes[nod[j]].alloc_other(ncompo);
02956 }
02957 }
02958 }
02959 }
02960 }
02961 }
02962
02963
02964
02965
02966
02967
02968
02969
02970
02971
02972
02973 void mechtop::alloc_enodes ()
02974 {
02975 if (Gtm->nen>0){
02976 enodes = new endnodem [Gtm->nen];
02977 }
02978 }
02979
02980
02981
02982
02983
02984
02985
02986
02987
02988
02989
02990 void mechtop::alloc_edges ()
02991 {
02992 if (Gtm->nged>0){
02993 edges = new edgem [Gtm->nged];
02994 }
02995 }
02996
02997
02998
02999
03000
03001
03002
03003
03004
03005
03006
03007
03008
03009 void mechtop::alloc_prep (long nn,long ne,elemtype *el_type)
03010 {
03011 long i;
03012 nodes = new node [nn];
03013 elements = new element [ne];
03014 for(i=0; i<ne; i++)
03015 elements[i].te = el_type[i];
03016 }
03017
03018
03019
03020
03021
03022
03023
03024
03025
03026
03027 void mechtop::alloc_meaning ()
03028 {
03029 long i;
03030
03031 for (i=0;i<nn;i++){
03032 nodes[i].alloc_meaning (i);
03033 }
03034 }
03035
03036
03037
03038
03039
03040
03041
03042
03043
03044
03045
03046 void mechtop::alloc_growstr ()
03047 {
03048 long i, ndofn;
03049
03050
03051 nodeforce = new double* [nn];
03052 memset(nodeforce, 0, sizeof(*nodeforce)*nn);
03053
03054 for (i=0;i<nn;i++){
03055 nodes[i].alloc_growstr (i);
03056 ndofn = give_ndofn(i);
03057 nodeforce[i] = new double[ndofn];
03058 memset(nodeforce[i], 0, sizeof(*nodeforce[i])*ndofn);
03059 }
03060 for (i=0;i<ne;i++){
03061 elements[i].alloc_growstr (i);
03062 }
03063 }
03064
03065
03066
03067
03068
03069
03070
03071
03072
03073
03074 void mechtop::define_meaning ()
03075 {
03076 long i;
03077 elemtype te;
03078
03079 for (i=0;i<ne;i++){
03080 if (Gtm->leso[i]==1){
03081 te = give_elem_type (i);
03082
03083 switch (te){
03084 case beam2d:{
03085 Beam2d->define_meaning (i);
03086 break;
03087 }
03088 case beam3d:{
03089 Beam3d->define_meaning (i);
03090 break;
03091 }
03092 case planeelementlq:{
03093 Pelq->define_meaning (i);
03094 break;
03095 }
03096 default:{
03097 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
03098 }
03099 }
03100 }
03101 }
03102 }
03103
03104
03105
03106
03107
03108
03109
03110
03111
03112
03113
03114
03115
03116
03117 void mechtop::gencodnumlagrmult (long &n)
03118 {
03119 long i,j,ndofn,nl,nmult,*cnl,*cnu;
03120
03121 n++;
03122
03123 for (i=0;i<nln;i++){
03124 nl=Gtm->lgnodes[i].nl;
03125 ndofn=Gtm->gnodes[Gtm->lgnodes[i].nodes[0]].ndofn;
03126 Gtm->lgnodes[i].ndofn=ndofn;
03127 Gtm->lgnodes[i].cn = new long* [nl+1];
03128
03129 if (Mp->tlm==1){
03130 nmult=2;
03131 Gtm->lgnodes[i].nmult = nmult;
03132 cnl = new long [ndofn];
03133 cnu = new long [ndofn];
03134
03135 Gtm->lgnodes[i].cn[0] = new long [nmult];
03136 Gtm->lgnodes[i].cn[0][0]=0;
03137 Gtm->lgnodes[i].cn[0][1]=0;
03138
03139 for (j=1;j<nl;j++){
03140 give_node_code_numbers (Gtm->lgnodes[i].nodes[j-1],cnl);
03141 give_node_code_numbers (Gtm->lgnodes[i].nodes[j],cnu);
03142
03143 Gtm->lgnodes[i].cn[j] = new long [nmult];
03144 if (cnl[0]>0 || cnl[2]>0 || cnu[0]>0 || cnu[2]>0){
03145 Gtm->lgnodes[i].cn[j][0]=n; n++;
03146 }
03147 if (cnl[1]>0 || cnu[1]>0){
03148 Gtm->lgnodes[i].cn[j][1]=n; n++;
03149 }
03150
03151 }
03152 Gtm->lgnodes[i].cn[nl] = new long [nmult];
03153 Gtm->lgnodes[i].cn[nl][0]=0;
03154 Gtm->lgnodes[i].cn[nl][1]=0;
03155
03156 delete [] cnl; delete [] cnu;
03157 }
03158
03159 if (Mp->tlm==2){
03160 nmult=4;
03161 Gtm->lgnodes[i].nmult = nmult;
03162 cnl = new long [ndofn];
03163 cnu = new long [ndofn];
03164
03165 Gtm->lgnodes[i].cn[0] = new long [nmult];
03166 Gtm->lgnodes[i].cn[0][0]=0;
03167 Gtm->lgnodes[i].cn[0][1]=0;
03168 Gtm->lgnodes[i].cn[0][2]=0;
03169 Gtm->lgnodes[i].cn[0][3]=0;
03170
03171 for (j=1;j<nl;j++){
03172 give_node_code_numbers (Gtm->lgnodes[i].nodes[j-1],cnl);
03173 give_node_code_numbers (Gtm->lgnodes[i].nodes[j],cnu);
03174
03175 Gtm->lgnodes[i].cn[j] = new long [nmult];
03176 if (cnl[0]>0 || cnl[4]>0 || cnu[0]>0 || cnu[4]>0){
03177 Gtm->lgnodes[i].cn[j][0]=n; n++;
03178 }
03179 if (cnl[1]>0 || cnl[3]>0 || cnu[1]>0 || cnu[3]>0){
03180 Gtm->lgnodes[i].cn[j][1]=n; n++;
03181 }
03182 if (cnl[2]>0 || cnu[2]>0){
03183 Gtm->lgnodes[i].cn[j][2]=n; n++;
03184 }
03185 if (cnl[5]>0 || cnu[5]>0){
03186 Gtm->lgnodes[i].cn[j][3]=n; n++;
03187 }
03188
03189 }
03190 Gtm->lgnodes[i].cn[nl] = new long [nmult];
03191 Gtm->lgnodes[i].cn[nl][0]=0;
03192 Gtm->lgnodes[i].cn[nl][1]=0;
03193 Gtm->lgnodes[i].cn[nl][2]=0;
03194 Gtm->lgnodes[i].cn[nl][3]=0;
03195
03196 delete [] cnl; delete [] cnu;
03197 }
03198 }
03199 n--;
03200 }
03201
03202
03203
03204
03205
03206
03207
03208
03209
03210
03211
03212
03213
03214
03215
03216
03217
03218
03219
03220
03221
03222
03223
03224
03225
03226
03227
03228
03229
03230
03231
03232
03233
03234
03235
03236
03237
03238
03239
03240
03241
03242
03243
03244
03245
03246
03247
03248
03249
03250
03251
03252
03253
03254
03255
03256
03257
03258
03259
03260
03261
03262
03263
03264
03265
03266
03267
03268
03269
03270
03271
03272
03273
03274
03275
03276
03277
03278
03279
03280
03281
03282
03283
03284
03285
03286
03287
03288
03289
03290
03291
03292
03293
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305
03306
03307
03308
03309
03310
03311
03312 long mechtop::compare_edges (long *enod, long *enodc, long nned)
03313 {
03314 long i, j, nsn = 0;
03315
03316 for (i = 0; i < nned; i++)
03317 {
03318 for (j = 0; j < nned; j++)
03319 {
03320 if (enod[i] == enodc[j])
03321 {
03322 nsn++;
03323 break;
03324 }
03325 }
03326 }
03327 if (nsn == nned)
03328 return 1;
03329
03330 return 0;
03331 }
03332
03333
03334
03335
03336
03337
03338
03339
03340
03341
03342
03343
03344
03345
03346
03347
03348
03349
03350
03351
03352
03353
03354
03355
03356
03357
03358
03359
03360
03361
03362
03363
03364
03365
03366
03367
03368
03369
03370
03371
03372
03373
03374
03375
03376
03377
03378
03379
03380
03381
03382
03383
03384
03385
03386
03387
03388
03389
03390
03391
03392
03393
03394
03395
03396
03397
03398
03399
03400
03401
03402
03403
03404
03405
03406
03407
03408
03409
03410
03411
03412
03413
03414
03415
03416
03417
03418
03419
03420
03421
03422
03423
03424
03425
03426
03427
03428
03429
03430
03431
03432
03433
03434
03435 void mechtop::save_elem_inidispl(long lcid, long *ifn)
03436 {
03437 long i, j, k, l, nne, ndofe, ndofn;
03438 ivector enod;
03439 vector r;
03440
03441
03442 for (i=0;i<ne;i++)
03443 {
03444 if (Gtm->leso[i] == 1)
03445
03446
03447 {
03448 ndofe = give_ndofe (i);
03449
03450 reallocv(ndofe, r);
03451
03452
03453 eldispl(lcid, i, r.a);
03454
03455
03456
03457
03458
03459 nne = give_nne(i);
03460 reallocv(nne, enod);
03461 give_elemnodes(i, enod);
03462 l = 0;
03463 for (j=0; j<nne; j++)
03464 {
03465 ndofn = give_ndofn(enod[j]);
03466 if (ifn[enod[j]])
03467 {
03468 for (k=0; k<ndofn; k++)
03469 {
03470 r[l] = nodes[enod[j]].nodval[k];
03471 l++;
03472 }
03473 }
03474 else
03475 l += ndofn;
03476 }
03477
03478
03479 elements[i].initdisplacement(r.a, ndofe);
03480 }
03481 }
03482 }
03483
03484
03485
03486
03487
03488
03489
03490
03491
03492
03493
03494
03495
03496
03497
03498
03499
03500
03501
03502 void mechtop::save_node_inidispl(long lcid, double time, double prev_time)
03503 {
03504 long i,j,k,l,m,ndofn;
03505 vector r;
03506
03507 for (i=0; i<nn; i++)
03508 {
03509 if (Gtm->lnso[i])
03510 {
03511 if (nodedispl[i])
03512 {
03513 ndofn = give_ndofn(i);
03514 reallocv(ndofn, r);
03515 noddispl(lcid, r.a, i);
03516 for (j=0;j<ndofn;j++)
03517 {
03518 k=Gtm->gnodes[i].tgf[j];
03519
03520 if (k < 0)
03521 continue;
03522
03523 l=Gtm->gf[k].getval_long (time);
03524 m=Gtm->gf[k].getval_long (prev_time);
03525
03526 if ((l < 1) && (m > 0))
03527 {
03528
03529 nodedispl[i][j] = r[j];
03530 }
03531 if ((l > 0) && (m < 1) && nodedispl[i])
03532 {
03533
03534
03535
03536 nodedispl[i][j] = 0.0;
03537 }
03538 if ((l > 0) && (m > 0) && nodedispl[i])
03539 {
03540
03541
03542 nodedispl[i][j] = 0.0;
03543 }
03544 }
03545 }
03546 }
03547 else
03548 {
03549 ndofn = give_ndofn(i);
03550 memset(nodes[i].nodval, 0, sizeof(*nodes[i].nodval)*ndofn);
03551 }
03552 }
03553 }
03554
03555
03556
03557
03558
03559
03560
03561
03562
03563
03564
03565
03566 void mechtop::clean_ip_new_elem()
03567 {
03568 long i, j, ipp, tnip;
03569
03570 for (i=0; i<ne; i++)
03571 {
03572 if (Gtm->leso[i] && (Gtm->gelements[i].auxinf == 0))
03573 {
03574 ipp = give_nip(i, 0, 0);
03575 tnip = give_tnip(i);
03576 for (j=0; j<tnip; j++)
03577 {
03578 Mm->ip[ipp].clean_strains(0);
03579 Mm->ip[ipp].clean_stresses(0);
03580
03581 ipp++;
03582 }
03583 }
03584 }
03585 }
03586
03587
03588
03589
03590
03591
03592
03593
03594
03595
03596
03597
03598
03599
03600
03601 void mechtop::give_noddispl_1d (ivector &nodes,vector &u)
03602 {
03603 long i;
03604
03605 for (i=0;i<nodes.n;i++){
03606 u[i]=nodedispl[nodes[i]][0];
03607 }
03608
03609 }
03610
03611
03612
03613
03614
03615
03616
03617
03618
03619
03620
03621
03622
03623
03624
03625
03626 void mechtop::give_noddispl_2d (ivector &nodes,vector &u,vector &v)
03627 {
03628 long i;
03629
03630 for (i=0;i<nodes.n;i++){
03631 u[i]=nodedispl[nodes[i]][0];
03632 v[i]=nodedispl[nodes[i]][1];
03633 }
03634 }
03635
03636
03637
03638
03639
03640
03641
03642
03643
03644
03645
03646
03647
03648
03649
03650
03651
03652 void mechtop::give_noddispl_3d (ivector &nodes,vector &u,vector &v,vector &w)
03653 {
03654 long i;
03655
03656 for (i=0;i<nodes.n;i++){
03657 u[i]=nodedispl[nodes[i]][0];
03658 v[i]=nodedispl[nodes[i]][1];
03659 w[i]=nodedispl[nodes[i]][2];
03660 }
03661
03662 }
03663
03664
03665
03666
03667
03668
03669
03670
03671
03672
03673
03674
03675 void mechtop::init_from_siftop (siftop *top)
03676 {
03677 long i,j;
03678
03679
03680 nn = top->nn;
03681
03682 ne = top->ne;
03683
03684 top->exporttop(Gtm);
03685
03686 if (Mespr==1) fprintf (stdout,"\n number of nodes %ld",nn);
03687 if (Mespr==1) fprintf (stdout,"\n number of elements %ld",ne);
03688 nodes = new node [nn];
03689 elements = new element [ne];
03690
03691 for (i=0; i < ne; i++){
03692 if (Gtm->leso[i]==1){
03693 for(j=0; j<Gtm->gelements[i].give_nne(); j++)
03694 Gtm->gelements[i].nodes[j]--;
03695 }
03696 }
03697 for (i=0;i<ne;i++){
03698 if (Gtm->leso[i]==1){
03699 elements[i].tm = new mattype[1];
03700 elements[i].idm = new long[1];
03701 elements[i].idm[0] = top->elements[i].prop-1;
03702 switch (top->elements[i].type)
03703 {
03704 case ISOLinear1D:
03705 elements[i].te = bar2d;
03706 if (Bar2d == NULL)
03707 Bar2d = new barel2d;
03708 break;
03709 case ISOQuadratic1D:
03710 elements[i].te = barq2d;
03711 if (Barq2d == NULL)
03712 Barq2d = new barelq2d;
03713 break;
03714 case TriangleLinear:
03715 elements[i].te = planeelementlt;
03716 if (Pelt == NULL)
03717 Pelt = new planeelemlt;
03718 break;
03719 case TriangleQuadratic:
03720 elements[i].te =planeelementqt;
03721 if (Peqt == NULL)
03722 Peqt = new planeelemqt;
03723 break;
03724 case ISOLinear2D:
03725 elements[i].te = planeelementlq;
03726 if (Pelq == NULL)
03727 Pelq = new planeelemlq;
03728 break;
03729 case ISOQuadratic2D:
03730 elements[i].te = planeelementqq;
03731 if (Peqq == NULL)
03732 Peqq = new planeelemqq;
03733 break;
03734 case TetrahedronLinear:
03735 elements[i].te = lineartet;
03736 if (Ltet == NULL)
03737 Ltet = new lintet;
03738 break;
03739 case TetrahedronQuadratic:
03740 elements[i].te = quadrtet;
03741 if (Qtet == NULL)
03742 Qtet = new quadtet;
03743 break;
03744 case VedgeLinear:
03745 elements[i].te = linearwed;
03746 if (Lwed == NULL)
03747 Lwed = new linwedge;
03748 break;
03749 case VedgeQuadratic:
03750 elements[i].te = quadrwed;
03751 if (Qwed == NULL)
03752 Qwed = new quadwedge;
03753 break;
03754 case ISOLinear3D:
03755 elements[i].te = linearhex;
03756 if (Lhex == NULL)
03757 Lhex = new linhex;
03758 break;
03759 case ISOQuadratic3D:
03760 elements[i].te = quadrhex;
03761 if (Qhex == NULL)
03762 Qhex = new quadhex;
03763 break;
03764 default:
03765 print_err("unknown type of element is used", __FILE__, __LINE__, __func__);
03766 abort();
03767 }
03768 }
03769 }
03770 }
03771
03772
03773
03774
03775
03776
03777
03778
03779
03780
03781
03782 void mechtop::clean_nodes ()
03783 {
03784 long i;
03785
03786 for (i=0;i<nn;i++){
03787 nodes[i].clean (Mb->nlc);
03788 }
03789 }
03790
03791
03792
03793
03794
03795
03796
03797
03798
03799
03800
03801
03802 void mechtop::initial_displ (long lcid)
03803 {
03804 long i,j,k,ndofe;
03805 vector r;
03806
03807
03808
03809 for (i=0;i<ne;i++){
03810
03811 j=Gtm->leso[i];
03812 k=Gtm->gelements[i].auxinf;
03813
03814 if (j==1 && k==0){
03815
03816
03817 ndofe = give_ndofe (i);
03818
03819 reallocv(ndofe, r);
03820
03821
03822 eldispl (lcid,i,r.a);
03823
03824 elements[i].initdisplacement (r.a,ndofe);
03825
03826 Gtm->gelements[i].auxinf=1;
03827 }
03828 }
03829 }
03830
03831
03832
03833
03834
03835
03836
03837
03838
03839
03840
03841
03842 void mechtop::save_nodval(long lcid)
03843 {
03844 long i, ndofn;
03845
03846 for (i=0;i<nn;i++){
03847 if (Gtm->lnso[i])
03848 {
03849 ndofn = give_ndofn (i);
03850
03851
03852
03853 noddispl(lcid, nodes[i].nodval, i);
03854 }
03855 }
03856 }
03857
03858
03859
03860
03861
03862
03863
03864
03865
03866
03867
03868
03869
03870
03871
03872
03873 void mechtop::restore_nodval(double *lhs,long n)
03874 {
03875 long i,j,ndofn;
03876 ivector cn;
03877
03878 nullv(lhs, n);
03879
03880 for (i=0;i<nn;i++){
03881 ndofn = Gtm->give_ndofn (i);
03882 if (ndofn < 0)
03883
03884 continue;
03885 reallocv(ndofn, cn);
03886 give_node_code_numbers (i,cn.a);
03887
03888 for (j=0;j<ndofn;j++){
03889 if (cn[j]>0){
03890 lhs[cn[j]-1]=nodes[i].nodval[j];
03891 }
03892 }
03893 }
03894 }
03895
03896
03897
03898
03899
03900
03901
03902
03903
03904
03905
03906
03907
03908
03909 void mechtop::save_nodforce(double *f)
03910 {
03911 long i, ndofn;
03912 vector nf;
03913
03914 for (i=0;i<nn;i++)
03915 {
03916 if (Gtm->lnso[i])
03917 {
03918 ndofn = give_ndofn (i);
03919 if (ndofn < 0)
03920
03921 continue;
03922 nf.n = ndofn;
03923 nf.a = nodeforce[i];
03924 select_nodforce(f, i, nf);
03925 }
03926 }
03927 nf.n = 0L;
03928 nf.a = NULL;
03929 }
03930
03931
03932
03933
03934
03935
03936
03937
03938
03939
03940
03941
03942
03943
03944
03945 void mechtop::restore_nodforce(double *rhs,long n)
03946 {
03947 long i,j,ndofn;
03948 ivector cn;
03949
03950 nullv(rhs, n);
03951
03952 for (i=0;i<nn;i++){
03953 ndofn = Gtm->give_ndofn (i);
03954 if (ndofn < 0)
03955
03956 continue;
03957 reallocv(ndofn, cn);
03958 give_node_code_numbers (i,cn.a);
03959
03960 for (j=0;j<ndofn;j++){
03961 if (cn[j]>0){
03962 rhs[cn[j]-1]=nodeforce[i][j];
03963 }
03964 }
03965 }
03966 }
03967
03968
03969
03970
03971
03972
03973
03974
03975
03976
03977
03978 long mechtop::mesh_check(void)
03979 {
03980 long i,j,err=0,kk,jj;
03981
03982
03983 for(i=0;i<Mt->ne;i++){
03984 for(j=0;j<Gtm->gelements[i].nne;j++){
03985
03986 jj = Gtm->gelements[i].nodes[j];
03987 for(kk = 0; kk < Gtm->gnodes[jj].ndofn;kk++){
03988
03989 if((Gtm->gelements[i].tgf) < (Gtm->gnodes[jj].tgf[kk])){
03990 print_err("different element node(%ld) time function number=%ld than element(%ld) time function number=%ld\n"
03991 " (node number=%ld, element number=%ld)",
03992 __FILE__, __LINE__, __func__, j+1,Gtm->gnodes[jj].tgf[kk]+1,i+1,Gtm->gelements[i].tgf+1,Gtm->gelements[i].nodes[j]+1,i+1);
03993 err = 1;
03994 return err;
03995 }
03996 }
03997 }
03998 }
03999 return err;
04000 }
04001
04002
04003
04004
04005
04006
04007
04008
04009
04010
04011
04012 double mechtop::give_domain_vol()
04013 {
04014 double ret = 0.0;
04015 long i;
04016
04017 for (i=0; i < ne; i++)
04018 ret += give_volume(i);
04019
04020 return ret;
04021 }