00001 #include <math.h>
00002 #include "elemparticle.h"
00003 #include "global.h"
00004 #include "globmat.h"
00005 #include "genfile.h"
00006 #include "element.h"
00007 #include "node.h"
00008 #include "intpoints.h"
00009 #include "loadcase.h"
00010
00011
00012 elemparticle::elemparticle (long cnne,long cdim)
00013 {
00014 long i;
00015
00016
00017 dim = cdim;
00018
00019
00020
00021 nne=cnne;
00022
00023
00024
00025 ndofe = nne*cdim;
00026
00027
00028
00029
00030
00031 nb=1;
00032
00033
00034
00035
00036
00037 nip = new long* [nb];
00038 for (i=0;i<nb;i++){
00039 nip[i] = new long [nb];
00040 }
00041
00042 nip[0][0]=1;
00043
00044 tnip=0;
00045 for (i=0;i<nb;i++){
00046 tnip+=nip[i][i];
00047 }
00048
00049 }
00050
00051 elemparticle::~elemparticle (void)
00052 {
00053 long i;
00054
00055 for (i=0;i<nb;i++){
00056 delete [] nip[i];
00057 }
00058 delete [] nip;
00059
00060 }
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077 void elemparticle::direction_vector_1d (long eid,long i,long j,vector &s,vector &x,vector &u)
00078 {
00079 s[0]=x[j]+u[j] - x[i]-u[i];
00080 }
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 void elemparticle::direction_vector_2d (long eid,long i,long j,vector &s,
00094 vector &x,vector &y,vector &u,vector &v)
00095 {
00096 s[0]=x[j]+u[j] - x[i]-u[i];
00097 s[1]=y[j]+v[j] - y[i]-v[i];
00098 }
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111 void elemparticle::direction_vector_3d (long eid,long i,long j,vector &s,
00112 vector &x,vector &y,vector &z,
00113 vector &u,vector &v,vector &w)
00114 {
00115 s[0]=x[j]+u[j] - x[i]-u[i];
00116 s[1]=y[j]+v[j] - y[i]-v[i];
00117 s[2]=z[j]+w[j] - z[i]-w[i];
00118 }
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129 void elemparticle::stiffmat_1d_kii (long ipp,matrix &k,vector &s)
00130 {
00131 double f,g,r;
00132
00133
00134 r = normv(s);
00135
00136
00137 f = Mm->give_first_derivative (ipp,r);
00138
00139
00140 g = Mm->give_second_derivative (ipp,r);
00141
00142
00143 k[0][0] = g*s[0]*s[0]/r/r + f*(1.0/r-s[0]*s[0]/r/r/r);
00144 }
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 void elemparticle::stiffmat_1d_kij (long ipp,matrix &k,vector &s)
00156 {
00157 double f,g,r;
00158
00159
00160 r = normv(s);
00161
00162
00163 f = Mm->give_first_derivative (ipp,r);
00164
00165
00166 g = Mm->give_second_derivative (ipp,r);
00167
00168
00169 k[0][0] = -1.0*g*s[0]*s[0]/r/r + f*(-1.0/r+s[0]*s[0]/r/r/r);
00170 }
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180 void elemparticle::stiffness_matrix_1d (long eid,matrix &sm)
00181 {
00182 long i,j,ipp;
00183 ivector nodes(nne);
00184 vector x(nne),u(nne),s(1);
00185 matrix k(1,1);
00186
00187
00188 Mt->give_elemnodes (eid,nodes);
00189
00190 Mt->give_node_coord1d (x,eid);
00191
00192
00193 Mt->give_noddispl_1d (nodes,u);
00194
00195
00196 ipp = Mt->elements[eid].ipp[0][0];
00197
00198 fillm (0.0, sm);
00199
00200 for (i=0;i<nne;i++){
00201 for (j=i+1;j<nne;j++){
00202 direction_vector_1d (eid,i,j,s,x,u);
00203
00204 stiffmat_1d_kii (ipp,k,s);
00205
00206
00207
00208 sm[i][i]+=k[0][0];
00209 sm[j][j]+=k[0][0];
00210
00211 stiffmat_1d_kij (ipp,k,s);
00212 sm[i][j]+=k[0][0];
00213 sm[j][i]+=k[0][0];
00214 }
00215 }
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237 fprintf (Out,"\n\n STIFFNESS MATRIX");
00238 for (i=0;i<nne;i++){
00239 fprintf (Out,"\n");
00240 for (j=0;j<nne;j++){
00241 fprintf (Out," %le",sm[i][j]);
00242 }
00243 }
00244 fprintf (Out,"\n");
00245
00246 if (sm[0][0]<1.0e-2){
00247 sm[0][0]+=1.0;
00248 sm[0][1]-=1.0;
00249 sm[1][0]-=1.0;
00250 sm[1][1]+=1.0;
00251 }
00252
00253 }
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 void elemparticle::stiffmat_2d_kii (long ipp,matrix &k,vector &s)
00270 {
00271 double f,g,r;
00272
00273
00274 r = normv(s);
00275
00276
00277 f = Mm->give_first_derivative (ipp,r);
00278
00279
00280 g = Mm->give_second_derivative (ipp,r);
00281
00282
00283 k[0][0] = g*s[0]*s[0]/r/r + f*(1.0/r-s[0]*s[0]/r/r/r);
00284 k[0][1] = g*s[0]*s[1]/r/r - f*s[0]*s[1]/r/r/r;
00285
00286 k[1][0] = k[0][1];
00287 k[1][1] = g*s[1]*s[1]/r/r + f*(1.0/r-s[1]*s[1]/r/r/r);
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299 void elemparticle::stiffmat_2d_kij (long ipp,matrix &k,vector &s)
00300 {
00301 double f,g,r;
00302
00303
00304 r = normv(s);
00305
00306
00307 f = Mm->give_first_derivative (ipp,r);
00308
00309
00310 g = Mm->give_second_derivative (ipp,r);
00311
00312
00313 k[0][0] = -1.0*g*s[0]*s[0]/r/r + f*(-1.0/r+s[0]*s[0]/r/r/r);
00314 k[0][1] = -1.0*g*s[0]*s[1]/r/r + f*s[0]*s[1]/r/r/r;
00315
00316 k[1][0] = k[0][1];
00317 k[1][1] = -1.0*g*s[1]*s[1]/r/r + f*(-1.0/r+s[1]*s[1]/r/r/r);
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 void elemparticle::stiffness_matrix_2d (long eid,matrix &sm)
00329 {
00330 long i,j,ipp;
00331 ivector nodes(nne);
00332 vector x(nne),y(nne),u(nne),v(nne),s(2);
00333 matrix k(2,2);
00334
00335
00336 Mt->give_elemnodes (eid,nodes);
00337
00338 Mt->give_node_coord2d (x,y,eid);
00339
00340
00341 Mt->give_noddispl_2d (nodes,u,v);
00342
00343
00344 ipp = Mt->elements[eid].ipp[0][0];
00345
00346 fillm (0.0,sm);
00347
00348 for (i=0;i<nne;i++){
00349 for (j=i+1;j<nne;j++){
00350 direction_vector_2d (eid,i,j,s,x,y,u,v);
00351
00352 stiffmat_2d_kii (ipp,k,s);
00353
00354
00355
00356 sm[i*2+0][i*2+0]+=k[0][0];
00357 sm[i*2+0][i*2+1]+=k[0][1];
00358 sm[i*2+1][i*2+0]+=k[1][0];
00359 sm[i*2+1][i*2+1]+=k[1][1];
00360
00361 sm[j*2+0][j*2+0]+=k[0][0];
00362 sm[j*2+0][j*2+1]+=k[0][1];
00363 sm[j*2+1][j*2+0]+=k[1][0];
00364 sm[j*2+1][j*2+1]+=k[1][1];
00365
00366 stiffmat_2d_kij (ipp,k,s);
00367
00368 sm[i*2+0][j*2+0]+=k[0][0];
00369 sm[i*2+0][j*2+1]+=k[0][1];
00370 sm[i*2+1][j*2+0]+=k[1][0];
00371 sm[i*2+1][j*2+1]+=k[1][1];
00372
00373 sm[j*2+0][i*2+0]+=k[0][0];
00374 sm[j*2+0][i*2+1]+=k[0][1];
00375 sm[j*2+1][i*2+0]+=k[1][0];
00376 sm[j*2+1][i*2+1]+=k[1][1];
00377 }
00378 }
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400 fprintf (Out,"\n\n STIFFNESS MATRIX");
00401 for (i=0;i<nne;i++){
00402 fprintf (Out,"\n");
00403 for (j=0;j<nne;j++){
00404 fprintf (Out," %le",sm[i][j]);
00405 }
00406 }
00407 fprintf (Out,"\n");
00408
00409 }
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423 void elemparticle::stiffmat_3d_kii (long ipp,matrix &k,vector &s)
00424 {
00425 double f,g,r;
00426
00427
00428 r = normv(s);
00429
00430
00431 f = Mm->give_first_derivative (ipp,r);
00432
00433
00434 g = Mm->give_second_derivative (ipp,r);
00435
00436
00437 k[0][0] = g*s[0]*s[0]/r/r + f*(1.0/r-s[0]*s[0]/r/r/r);
00438 k[0][1] = g*s[0]*s[1]/r/r - f*s[0]*s[1]/r/r/r;
00439 k[0][2] = g*s[0]*s[2]/r/r - f*s[0]*s[2]/r/r/r;
00440
00441 k[1][0] = k[0][1];
00442 k[1][1] = g*s[1]*s[1]/r/r + f*(1.0/r-s[1]*s[1]/r/r/r);
00443 k[1][2] = g*s[1]*s[2]/r/r - f*s[1]*s[2]/r/r/r;
00444
00445 k[2][0] = k[0][2];
00446 k[2][1] = k[1][2];
00447 k[2][2] = g*s[2]*s[2]/r/r + f*(1.0/r-s[2]*s[2]/r/r/r);
00448 }
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459 void elemparticle::stiffmat_3d_kij (long ipp,matrix &k,vector &s)
00460 {
00461 double f,g,r;
00462
00463
00464 r = normv (s);
00465
00466
00467 f = Mm->give_first_derivative (ipp,r);
00468
00469
00470 g = Mm->give_second_derivative (ipp,r);
00471
00472
00473 k[0][0] = -1.0*g*s[0]*s[0]/r/r + f*(-1.0/r+s[0]*s[0]/r/r/r);
00474 k[0][1] = -1.0*g*s[0]*s[1]/r/r + f*s[0]*s[1]/r/r/r;
00475 k[0][2] = -1.0*g*s[0]*s[2]/r/r + f*s[0]*s[2]/r/r/r;
00476
00477 k[1][0] = k[0][1];
00478 k[1][1] = -1.0*g*s[1]*s[1]/r/r + f*(-1.0/r+s[1]*s[1]/r/r/r);
00479 k[1][2] = -1.0*g*s[1]*s[2]/r/r + f*s[1]*s[2]/r/r/r;
00480
00481 k[2][0] = k[0][2];
00482 k[2][1] = k[1][2];
00483 k[2][2] = -1.0*g*s[2]*s[2]/r/r + f*(-1.0/r+s[2]*s[2]/r/r/r);
00484 }
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494 void elemparticle::stiffness_matrix_3d (long eid,matrix &sm)
00495 {
00496 long i,j,ipp,ii,jj;
00497 ivector nodes(nne);
00498 vector x(nne),y(nne),u(nne),v(nne),s(2);
00499 matrix k(3,3);
00500
00501
00502 Mt->give_elemnodes (eid,nodes);
00503
00504 Mt->give_node_coord2d (x,y,eid);
00505
00506
00507 Mt->give_noddispl_2d (nodes,u,v);
00508
00509
00510 ipp = Mt->elements[eid].ipp[0][0];
00511
00512 fillm (0.0,sm);
00513
00514 for (i=0;i<nne;i++){
00515 for (j=i+1;j<nne;j++){
00516 direction_vector_2d (eid,i,j,s,x,y,u,v);
00517
00518 stiffmat_2d_kii (ipp,k,s);
00519 for (ii=0; ii<3; ii++)
00520 {
00521 for(jj=0; jj<3; jj++)
00522 {
00523 sm[i*3+ii][i*3+jj] += k[ii][jj];
00524 sm[j*3+ii][j*3+jj] += k[ii][jj];
00525 }
00526 }
00527
00528 stiffmat_2d_kij (ipp,k,s);
00529 for (ii=0; ii<3; ii++)
00530 {
00531 for(jj=0; jj<3; jj++)
00532 {
00533 sm[i*3+ii][j*3+jj] += k[ii][jj];
00534 sm[j*3+ii][i*3+jj] += k[ii][jj];
00535 }
00536 }
00537 }
00538 }
00539 fprintf (Out,"\n\n STIFFNESS MATRIX");
00540 for (i=0;i<nne;i++){
00541 fprintf (Out,"\n");
00542 for (j=0;j<nne;j++){
00543 fprintf (Out," %le",sm[i][j]);
00544 }
00545 }
00546 fprintf (Out,"\n");
00547
00548 }
00549
00550
00551
00552
00553
00554
00555
00556
00557 void elemparticle::res_stiffness_matrix (long eid,matrix &sm)
00558 {
00559 switch (dim){
00560 case 1:{
00561 stiffness_matrix_1d (eid,sm);
00562 break;
00563 }
00564 case 2:{
00565 stiffness_matrix_2d (eid,sm);
00566 break;
00567 }
00568 case 3:{
00569 stiffness_matrix_3d (eid,sm);
00570 break;
00571 }
00572 default:{
00573 fprintf (stderr,"\n\n unknown dimension of problem is required in function res_stiffness_matrix (file %s, line %d).\n",__FILE__,__LINE__);
00574 }
00575 }
00576 }
00577
00578
00579 void elemparticle::forces_1d (long ipp,vector &fij,vector &s)
00580 {
00581 double r,d;
00582
00583
00584 r = normv(s);
00585
00586
00587 d = Mm->give_first_derivative (ipp,r);
00588
00589 fij[0] = -s[0]*d/r;
00590 }
00591
00592 void elemparticle::inter_forces_1d (long eid,vector &f)
00593 {
00594 long i,j,ipp;
00595 ivector nodes(nne);
00596 vector x(nne),u(nne),s(1),fij(1);
00597
00598
00599 Mt->give_elemnodes (eid,nodes);
00600
00601 Mt->give_node_coord1d (x,eid);
00602
00603 Mt->give_noddispl_1d (nodes,u);
00604
00605 ipp = Mt->elements[eid].ipp[0][0];
00606
00607 fprintf (Out,"\n u %le %le\n",u[0],u[1]);
00608
00609 fillv (0.0,f);
00610
00611 for (i=0;i<nne;i++){
00612 for (j=i;j<nne;j++){
00613 if (i != j){
00614 direction_vector_1d (eid,i,j,s,x,u);
00615 forces_1d (ipp,fij,s);
00616
00617 fprintf (Out,"\n s %le fij %le\n",s[0],fij[0]);
00618
00619 f[i]+=fij[0];
00620 f[j]-=fij[0];
00621 }
00622 }
00623 }
00624 copyv(f.a, Mm->ip[ipp].stress, f.n);
00625 }
00626
00627
00628 void elemparticle::forces_2d (long ipp,vector &fij,vector &s)
00629 {
00630 double r,d;
00631
00632
00633 r = normv(s);
00634
00635 d = Mm->give_first_derivative (ipp,r);
00636
00637 fij[0] = -s[0]*d/r;
00638 fij[1] = -s[1]*d/r;
00639 }
00640
00641 void elemparticle::inter_forces_2d (long eid,vector &f)
00642 {
00643 long i,j,ipp,ii;
00644 ivector nodes(nne);
00645 vector x(nne),y(nne),u(nne),v(nne),s(2),fij(2);
00646
00647
00648 Mt->give_elemnodes (eid,nodes);
00649
00650 Mt->give_node_coord2d (x,y,eid);
00651
00652 Mt->give_noddispl_2d (nodes,u, v);
00653
00654 ipp = Mt->elements[eid].ipp[0][0];
00655
00656 fillv (0.0,f);
00657 for (i=0;i<nne;i++)
00658 {
00659 for (j=i+1;j<nne;j++)
00660 {
00661 direction_vector_2d (eid,i,j,s,x,y,u,v);
00662 forces_2d (ipp,fij,s);
00663 for(ii=0; ii<2; ii++)
00664 {
00665 f[i*2+ii]+=fij[ii];
00666 f[j*2+ii]-=fij[ii];
00667 }
00668 }
00669 }
00670 }
00671
00672
00673 void elemparticle::forces_3d (long ipp,vector &fij,vector &s)
00674 {
00675 double r,d;
00676
00677
00678 r = normv(s);
00679
00680 d = Mm->give_first_derivative (ipp,r);
00681
00682 fij[0] = -s[0]*d/r;
00683 fij[1] = -s[1]*d/r;
00684 fij[2] = -s[2]*d/r;
00685 }
00686
00687 void elemparticle::inter_forces_3d (long eid,vector &f)
00688 {
00689 long i,j,ii,ipp;
00690 ivector nodes(nne);
00691 vector x(nne),y(nne),z(nne),u(nne),v(nne),w(nne),s(3),fij(3);
00692
00693
00694 Mt->give_elemnodes (eid,nodes);
00695
00696 Mt->give_node_coord3d (x,y,z,eid);
00697
00698 Mt->give_noddispl_3d (nodes,u,v,w);
00699
00700 ipp = Mt->elements[eid].ipp[0][0];
00701
00702 fillv (0.0,f);
00703 for (i=0;i<nne;i++)
00704 {
00705 for (j=i+1;j<nne;j++)
00706 {
00707 direction_vector_3d (eid,i,j,s,x,y,z,u,v,w);
00708 forces_3d (ipp,fij,s);
00709 for(ii=0; ii<3; ii++)
00710 {
00711 f[i*3+ii]+=fij[ii];
00712 f[j*3+ii]-=fij[ii];
00713 }
00714 }
00715 }
00716 }
00717
00718
00719 void elemparticle::res_internal_forces (long eid,vector &ifor)
00720 {
00721 switch (dim){
00722 case 1:{
00723 inter_forces_1d (eid,ifor);
00724 break;
00725 }
00726 case 2:{
00727 inter_forces_2d (eid,ifor);
00728 break;
00729 }
00730 case 3:{
00731 inter_forces_3d (eid,ifor);
00732 break;
00733 }
00734 default:{
00735 fprintf (stderr,"\n\n unknown dimension of problem is required in function res_stiffness_matrix (file %s, line %d).\n",__FILE__,__LINE__);
00736 }
00737 }
00738 }