00001 #include "argyristriangle.h"
00002
00003 #include "global.h"
00004 #include "gadaptivity.h"
00005 #include "globmat.h"
00006 #include "genfile.h"
00007 #include "element.h"
00008 #include "node.h"
00009 #include "intpoints.h"
00010 #include <stdlib.h>
00011 #include <math.h>
00012
00013
00014
00015
00016 ArgyrisTriangle::ArgyrisTriangle (void)
00017 {
00018 long i,j;
00019
00020
00021 nne=3;
00022
00023
00024 ndofe = 21;
00025
00026 tncomp=3;
00027
00028 napfun=21;
00029
00030 ned=3;
00031
00032 nned=2;
00033
00034
00035 intordmm = 6;
00036
00037 ssst=platek;
00038
00039
00040
00041 nb = 1;
00042
00043
00044 ncomp = new long [nb];
00045 ncomp[0] = 3;
00046
00047
00048 cncomp = new long [nb];
00049 cncomp[0] = 3;
00050
00051
00052
00053 nip = new long* [nb];
00054 intordsm = new long* [nb];
00055 for (i=0;i<nb;i++){
00056 nip[i] = new long [nb];
00057 intordsm[i] = new long [nb];
00058 }
00059
00060 nip[0][0]=12;
00061
00062
00063 tnip=0;
00064 for (i=0;i<nb;i++){
00065 for (j=0;j<nb;j++){
00066 tnip+=nip[i][j];
00067 }
00068 }
00069
00070 intordsm[0][0]=12;
00071 }
00072
00073 ArgyrisTriangle::~ArgyrisTriangle (void)
00074 {
00075 long i;
00076
00077 for (i=0;i<nb;i++){
00078 delete [] intordsm[i];
00079 delete [] nip[i];
00080 }
00081 delete [] intordsm;
00082 delete [] nip;
00083
00084 delete [] cncomp;
00085 delete [] ncomp;
00086 }
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 double ArgyrisTriangle::approx (vector &areacoord,vector &nodval)
00098 {
00099 double f;
00100 scprd (areacoord,nodval,f);
00101 return f;
00102 }
00103
00104
00105
00106
00107
00108
00109
00110
00111 double ArgyrisTriangle::approx_nat (double xi,double eta,vector &nodval)
00112 {
00113 double f;
00114 vector areacoord(3);
00115 areacoord[0]=xi;
00116 areacoord[1]=eta;
00117 areacoord[2]=1.0-areacoord[0]-areacoord[1];
00118 scprd (areacoord,nodval,f);
00119 return f;
00120 }
00121
00122 void ArgyrisTriangle::bf_matrix (matrix& n, vector& x, vector& y, vector& areacoord)
00123 {
00124 double pointX = 0.0, pointY = 0.0;
00125
00126 double* xCoords = x.a;
00127 double* yCoords = y.a;
00128 double* pointAreaCoords = areacoord.a;
00129
00130 convertAreaCoordinatesToNaturalCoordinates(xCoords, yCoords, pointAreaCoords, &pointX, &pointY);
00131
00132 for (int i = 0; i < 3; i++)
00133 {
00134 n[0][i*7 + 0] = get_Hi(i, xCoords, yCoords, pointX, pointY);
00135 n[0][i*7 + 1] = get_Hix(i, xCoords, yCoords, pointX, pointY);
00136 n[0][i*7 + 2] = get_Hiy(i, xCoords, yCoords, pointX, pointY);
00137 n[0][i*7 + 3] = get_Hixx(i, xCoords, yCoords, pointX, pointY);
00138 n[0][i*7 + 4] = get_Hiyy(i, xCoords, yCoords, pointX, pointY);
00139 n[0][i*7 + 5] = get_Hixy(i, xCoords, yCoords, pointX, pointY);
00140 n[0][i*7 + 6] = get_HjkNi(i, xCoords, yCoords, pointX, pointY);
00141 }
00142 }
00143
00144 void ArgyrisTriangle::geom_matrix (matrix& gm, vector& x, vector& y, vector& areacoord)
00145 {
00146 double pointX = 0.0, pointY = 0.0;
00147
00148 double* xCoords = x.a;
00149 double* yCoords = y.a;
00150 double* pointAreaCoords = areacoord.a;
00151
00152 convertAreaCoordinatesToNaturalCoordinates(xCoords, yCoords, pointAreaCoords, &pointX, &pointY);
00153
00154 for (int d = 0; d < 3; d++)
00155 {
00156 Derivatives deriv = (Derivatives)d;
00157 for (int i = 0; i < 3; i++)
00158 {
00159 gm[d][i*7 + 0] = get_Hi(i, deriv, xCoords, yCoords, pointX, pointY);
00160 gm[d][i*7 + 1] = get_Hix(i, deriv, xCoords, yCoords, pointX, pointY);
00161 gm[d][i*7 + 2] = get_Hiy(i, deriv, xCoords, yCoords, pointX, pointY);
00162 gm[d][i*7 + 3] = get_Hixx(i, deriv, xCoords, yCoords, pointX, pointY);
00163 gm[d][i*7 + 4] = get_Hiyy(i, deriv, xCoords, yCoords, pointX, pointY);
00164 gm[d][i*7 + 5] = get_Hixy(i, deriv, xCoords, yCoords, pointX, pointY);
00165 gm[d][i*7 + 6] = get_HjkNi(i, deriv, xCoords, yCoords, pointX, pointY);
00166 }
00167 }
00168 }
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 void ArgyrisTriangle::dmat (matrix& d, double t)
00281 {
00282 double c;
00283
00284 c = t*t*t;
00285
00286 d[0][0] = c*d[0][0]; d[0][1] = c*d[0][1]; d[0][2] = c*d[0][2];
00287 d[1][0] = c*d[1][0]; d[1][1] = c*d[1][1]; d[1][2] = c*d[1][2];
00288 d[2][0] = c*d[2][0]; d[2][1] = c*d[2][1]; d[2][2] = c*d[2][2];
00289 }
00290
00291
00292
00293
00294
00295
00296 void ArgyrisTriangle::transf_matrix (ivector &nodes,matrix &tmat)
00297 {
00298 long i,n,m;
00299 n=nodes.n;
00300 m=tmat.m;
00301 for (i=0;i<m;i++){
00302 tmat[i][i]=1.0;
00303 }
00304
00305 for (i=0;i<n;i++){
00306 if (Mt->nodes[nodes[i]].transf>0){
00307 tmat[i*3+1][i*3+1]=Mt->nodes[nodes[i]].e1[0]; tmat[i*3+1][i*3+2]=Mt->nodes[nodes[i]].e2[0];
00308 tmat[i*3+2][i*3+1]=Mt->nodes[nodes[i]].e1[1]; tmat[i*3+2][i*3+2]=Mt->nodes[nodes[i]].e2[1];
00309 }
00310 }
00311 }
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323 void ArgyrisTriangle::stiffness_matrix (long eid,long ri,long ci,matrix &sm,vector &x,vector &y)
00324 {
00325 long i,ii,jj,ipp;
00326 double jac,det,thick;
00327 ivector nodes(nne);
00328 vector l(3),t(nne),gp1,gp2,w,areacoord(3);
00329 matrix gm,d(tncomp,tncomp);
00330
00331 Mt->give_elemnodes (eid,nodes);
00332 Mc->give_thickness (eid,nodes,t);
00333
00334
00335 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00336
00337 fillm (0.0,sm);
00338
00339 ii=0;
00340 {
00341 allocm (ncomp[ii],ndofe,gm);
00342 jj=0;
00343 if (intordsm[ii][jj] != 0)
00344 {
00345 allocv (intordsm[ii][jj],w);
00346 allocv (intordsm[ii][jj],gp1);
00347 allocv (intordsm[ii][jj],gp2);
00348
00349 gauss_points_tr(gp1.a, gp2.a, w.a, intordsm[ii][jj]);
00350
00351 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00352
00353 for (i=0;i<intordsm[ii][jj];i++)
00354 {
00355 areacoord[0]=gp1[i];
00356 areacoord[1]=gp2[i];
00357 areacoord[2]=1.0-areacoord[0]-areacoord[1];
00358
00359 thick = approx (areacoord,t);
00360
00361
00362 geom_matrix (gm,x,y,areacoord);
00363
00364
00365 Mm->matstiff(d, ipp);
00366 dmat(d, thick);
00367
00368 jac = w[i]*det;
00369
00370
00371 bdbjac(sm, gm, d, gm, jac);
00372
00373 ipp++;
00374 }
00375
00376 destrv (gp2); destrv (gp1); destrv (w);
00377 }
00378 destrm (gm);
00379 }
00380
00381 }
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391 void ArgyrisTriangle::res_stiffness_matrix (long eid,matrix &sm)
00392 {
00393 long transf;
00394 ivector nodes(nne);
00395 vector x(nne),y(nne);
00396 Mt->give_node_coord2d (x,y,eid);
00397 Mt->give_elemnodes (eid,nodes);
00398
00399 stiffness_matrix (eid,0,0,sm,x,y);
00400
00401 fprintf (Out,"\n\n matice tuhosti \n");
00402 for (long i=0;i<21;i++){
00403 fprintf (Out,"\n");
00404 for (long j=0;j<21;j++){
00405 fprintf (Out,"\n r %2ld %2ld %le",i+1,j+1,sm[i][j]);
00406 }
00407 }
00408 fprintf (Out,"\n\n konec matice tuhosti \n");
00409
00410
00411 transf = Mt->locsystems (nodes);
00412 if (transf>0){
00413 matrix tmat (ndofe,ndofe);
00414 transf_matrix (nodes,tmat);
00415 glmatrixtransf (sm,tmat);
00416 }
00417 }
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428 void ArgyrisTriangle::mass_matrix (long eid,matrix &mm,vector &x,vector &y)
00429 {
00430 long i;
00431 double det,ww,thick,rho;
00432 ivector nodes(nne);
00433 vector l(3),gp1(intordmm),gp2(intordmm),w(intordmm),t(nne),dens(nne);
00434 matrix n(1,ndofe);
00435
00436 Mt->give_elemnodes (eid,nodes);
00437 Mc->give_thickness (eid,nodes,t);
00438 Mc->give_density (eid,nodes,dens);
00439
00440 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00441
00442 gauss_points_tr (gp1.a,gp2.a,w.a,intordmm);
00443
00444 fillm (0.0,mm);
00445
00446 for (i=0;i<intordmm;i++){
00447 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00448 ww=w[i];
00449
00450 bf_matrix (n,x,y,l);
00451
00452 thick = t[0]*l[0]+t[1]*l[1]+t[2]*l[2];
00453 rho = dens[0]*l[0]+dens[1]*l[1]+dens[2]*l[2];
00454
00455 ww*=det*thick*rho;
00456
00457 nnj (mm.a,n.a,ww,1,ndofe);
00458 }
00459 }
00460
00461 void ArgyrisTriangle::res_mass_matrix (long eid,matrix &mm)
00462 {
00463 vector x(nne),y(nne);
00464 Mt->give_node_coord2d (x,y,eid);
00465 mass_matrix (eid,mm,x,y);
00466 }
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 void ArgyrisTriangle::load_matrix (long eid,matrix &lm)
00478 {
00479 long i;
00480 double det,ww,thick;
00481 ivector nodes(nne);
00482 vector x(nne),y(nne),l(3),gp1(intordmm),gp2(intordmm),w(intordmm),t(nne);
00483 matrix n(1,ndofe);
00484
00485 Mt->give_elemnodes (eid,nodes);
00486 Mc->give_thickness (eid,nodes,t);
00487 Mt->give_node_coord2d (x,y,eid);
00488
00489 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00490
00491 gauss_points_tr (gp1.a,gp2.a,w.a,intordmm);
00492
00493 fillm (0.0,lm);
00494
00495 for (i=0;i<intordmm;i++){
00496 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00497 ww=w[i];
00498
00499 bf_matrix (n,x,y,l);
00500
00501 thick = t[0]*l[0]+t[1]*l[1]+t[2]*l[2];
00502
00503 ww*=det*thick;
00504
00505 nnj (lm.a,n.a,ww,1,ndofe);
00506 }
00507 }
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520 void ArgyrisTriangle::ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00521 {
00522 long i,ii,ipp;
00523 vector gp1,gp2,w,eps(tncomp),areacoord(3);
00524 matrix gm(tncomp,ndofe);
00525
00526
00527 for (ii=0;ii<nb;ii++){
00528
00529 allocv (intordsm[ii][ii],gp1);
00530 allocv (intordsm[ii][ii],gp2);
00531 allocv (intordsm[ii][ii],w);
00532
00533 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00534
00535 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00536
00537 for (i=0;i<intordsm[ii][ii];i++){
00538 areacoord[0]=gp1[i];
00539 areacoord[1]=gp2[i];
00540 areacoord[2]=1.0-areacoord[0]-areacoord[1];
00541
00542
00543 geom_matrix (gm,x,y,areacoord);
00544
00545 mxv (gm,r,eps);
00546
00547 Mm->storestrain (lcid,ipp,eps);
00548 ipp++;
00549 }
00550 destrv (w); destrv (gp1); destrv (gp2);
00551 }
00552 }
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562 void ArgyrisTriangle::res_ip_strains (long lcid,long eid)
00563 {
00564 vector aux,x(nne),y(nne),r(ndofe);
00565 ivector nodes(nne);
00566 matrix tmat;
00567
00568 Mt->give_node_coord2d (x,y,eid);
00569 Mt->give_elemnodes (eid,nodes);
00570 eldispl (lcid,eid,r.a);
00571
00572
00573 long transf = Mt->locsystems (nodes);
00574 if (transf>0){
00575 allocv (ndofe,aux);
00576 allocm (ndofe,ndofe,tmat);
00577 transf_matrix (nodes,tmat);
00578
00579 lgvectortransf (aux,r,tmat);
00580 copyv (aux,r);
00581 destrv (aux);
00582 destrm (tmat);
00583 }
00584
00585 ip_strains (lcid,eid,0,0,x,y,r);
00586
00587 }
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 void ArgyrisTriangle::nod_strains_ip (long lcid,long eid,long ri,long ci)
00599 {
00600 long i,j,ipp;
00601 ivector ipnum(nne),nod(nne);
00602 vector eps(tncomp);
00603
00604
00605
00606 ipp=Mt->elements[eid].ipp[ri][ci];
00607 nodip_planelt (ipp,intordsm[0][0],ipnum);
00608
00609
00610 Mt->give_elemnodes (eid,nod);
00611
00612 for (i=0;i<nne;i++){
00613
00614 Mm->givestrain (lcid,ipnum[i],eps);
00615
00616
00617 j=nod[i];
00618 Mt->nodes[j].storestrain (lcid,0,eps);
00619 }
00620
00621 }
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631 void ArgyrisTriangle::res_ip_stresses (long lcid,long eid)
00632 {
00633 long ri,ci;
00634 ri=0;
00635 ci=0;
00636
00637 compute_nlstress (lcid,eid,ri,ci);
00638 }
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649 void ArgyrisTriangle::nod_stresses_ip (long lcid,long eid,long ri,long ci)
00650 {
00651 long i,j,ipp;
00652 ivector ipnum(nne),nod(nne);
00653 vector sig(tncomp);
00654
00655
00656
00657 ipp=Mt->elements[eid].ipp[ri][ci];
00658 nodip_planelt (ipp,intordsm[0][0],ipnum);
00659
00660
00661 Mt->give_elemnodes (eid,nod);
00662
00663 for (i=0;i<nne;i++){
00664
00665 Mm->givestress (lcid,ipnum[i],sig);
00666
00667
00668 j=nod[i];
00669 Mt->nodes[j].storestress (lcid,0,sig);
00670 }
00671
00672 }
00673
00674
00675 void ArgyrisTriangle::strains (long lcid,long eid,long ri,long ci)
00676 {
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734 }
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744 void ArgyrisTriangle::nodecoord (vector &xi,vector &eta)
00745 {
00746 xi[0] = 0.0; eta[0] = 0.0;
00747 xi[1] = 1.0; eta[1] = 0.0;
00748 xi[2] = 0.0; eta[2] = 1.0;
00749 }
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760 void ArgyrisTriangle::appval (double xi,double eta,long fi,long nc,vector &eps,double **val)
00761 {
00762 long i,j,k;
00763 vector nodval(nne),areacoord(3);
00764
00765 areacoord[0]=xi; areacoord[1]=eta; areacoord[2]=1.0-areacoord[0]-areacoord[1];
00766 k=0;
00767 for (i=fi;i<fi+nc;i++){
00768 for (j=0;j<nne;j++){
00769 nodval[j]=val[j][i];
00770 }
00771 eps[k]=approx (areacoord,nodval);
00772 k++;
00773 }
00774 }
00775
00776
00777 void ArgyrisTriangle::stresses (long lcid,long eid,long ri,long ci)
00778 {
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841 }
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857 void ArgyrisTriangle::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00858 {
00859 integratedquant iq;
00860 iq=locstress;
00861
00862
00863 compute_nlstress (lcid,eid,ri,ci);
00864
00865 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00866 }
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881 void ArgyrisTriangle::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00882 {
00883 integratedquant iq;
00884 iq=nonlocstress;
00885
00886
00887 compute_nonloc_nlstress (lcid,eid,ri,ci);
00888
00889 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00890 }
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905 void ArgyrisTriangle::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00906 {
00907 integratedquant iq;
00908 iq=stressincr;
00909
00910
00911 compute_nlstressincr (lcid,eid,ri,ci);
00912
00913 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00914 }
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928 void ArgyrisTriangle::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y)
00929 {
00930 integratedquant iq;
00931 iq=eigstress;
00932
00933
00934 if (Mp->eigstrcomp)
00935 compute_eigstress (lcid,eid,ri,ci);
00936
00937
00938 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y);
00939 }
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952 void ArgyrisTriangle::res_internal_forces (long lcid,long eid,vector &ifor)
00953 {
00954 long transf;
00955 ivector nodes (nne);
00956 vector v(ndofe),x(nne),y(nne);
00957
00958 Mt->give_node_coord2d (x,y,eid);
00959 internal_forces (lcid,eid,0,0,ifor,x,y);
00960
00961
00962 Mt->give_elemnodes (eid,nodes);
00963 transf = Mt->locsystems (nodes);
00964 if (transf>0){
00965 matrix tmat (ndofe,ndofe);
00966 transf_matrix (nodes,tmat);
00967 glvectortransf (ifor,v,tmat);
00968 copyv (v,ifor);
00969 }
00970 }
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983 void ArgyrisTriangle::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
00984 {
00985 long transf;
00986 ivector nodes (nne);
00987 vector v(ndofe),x(nne),y(nne);
00988
00989 Mt->give_node_coord2d (x,y,eid);
00990 nonloc_internal_forces (lcid,eid,0,0,ifor,x,y);
00991
00992
00993 Mt->give_elemnodes (eid,nodes);
00994 transf = Mt->locsystems (nodes);
00995 if (transf>0){
00996 matrix tmat (ndofe,ndofe);
00997 transf_matrix (nodes,tmat);
00998 glvectortransf (ifor,v,tmat);
00999 copyv (v,ifor);
01000 }
01001 }
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014 void ArgyrisTriangle::res_incr_internal_forces (long lcid,long eid,vector &ifor)
01015 {
01016 long transf;
01017 ivector nodes (nne);
01018 vector v(ndofe),x(nne),y(nne);
01019
01020 Mt->give_node_coord2d (x,y,eid);
01021 incr_internal_forces (lcid,eid,0,0,ifor,x,y);
01022
01023
01024 Mt->give_elemnodes (eid,nodes);
01025 transf = Mt->locsystems (nodes);
01026 if (transf>0){
01027 matrix tmat (ndofe,ndofe);
01028 transf_matrix (nodes,tmat);
01029 glvectortransf (ifor,v,tmat);
01030 copyv (v,ifor);
01031 }
01032 }
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045 void ArgyrisTriangle::res_eigstrain_forces (long lcid,long eid,vector &nfor)
01046 {
01047 long transf;
01048 ivector nodes (nne);
01049 vector v(ndofe),x(nne),y(nne);
01050
01051 Mt->give_node_coord2d (x,y,eid);
01052 eigstrain_forces (lcid,eid,0,0,nfor,x,y);
01053
01054
01055 Mt->give_elemnodes (eid,nodes);
01056 transf = Mt->locsystems (nodes);
01057 if (transf>0){
01058 matrix tmat (ndofe,ndofe);
01059 transf_matrix (nodes,tmat);
01060 glvectortransf (nfor,v,tmat);
01061 copyv (v,nfor);
01062 }
01063 }
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076 void ArgyrisTriangle::compute_nlstress (long lcid,long eid,long ri,long ci)
01077 {
01078 long i,ii,ipp;
01079 double thick;
01080 ivector nodes(nne);
01081 vector sig(5),gp1,gp2,w,l(3),t(nne);
01082
01083 Mt->give_elemnodes (eid,nodes);
01084 Mc->give_thickness (eid,nodes,t);
01085
01086 for (ii=0;ii<nb;ii++){
01087 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01088 allocv (intordsm[ii][ii],w);
01089 allocv (intordsm[ii][ii],gp1);
01090 allocv (intordsm[ii][ii],gp2);
01091 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[0][0]);
01092
01093 for (i=0;i<intordsm[ii][ii];i++){
01094
01095 if (Mp->strcomp==1){
01096 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
01097 thick = t[0]*l[0]+t[1]*l[1]+t[2]*l[2];
01098
01099 Mm->computenlstresses (ipp);
01100
01101 Mm->givestress (lcid,ipp,sig);
01102 sig[0]*=thick*thick*thick;
01103 sig[1]*=thick*thick*thick;
01104 sig[2]*=thick*thick*thick;
01105 sig[3]*=thick*5.0/6.0;
01106 sig[4]*=thick*5.0/6.0;
01107 Mm->storestress (lcid,ipp,sig);
01108
01109 ipp++;
01110 }
01111 }
01112 destrv (gp1);
01113 destrv (gp2);
01114 destrv (w);
01115 }
01116 }
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130 void ArgyrisTriangle::local_values(long lcid,long eid,long ri,long ci)
01131 {
01132 long i,ii,ipp;
01133
01134 for (ii=0;ii<nb;ii++){
01135 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01136
01137 for (i=0;i<intordsm[ii][ii];i++){
01138
01139 if (Mp->strcomp==1)
01140 Mm->computenlstresses (ipp);
01141 ipp++;
01142 }
01143 }
01144 }
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157 void ArgyrisTriangle::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01158 {
01159 long i,ii,ipp;
01160
01161 for (ii=0;ii<nb;ii++){
01162 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01163
01164 for (i=0;i<intordsm[ii][ii];i++){
01165
01166 if (Mp->strcomp==1)
01167 Mm->compnonloc_nlstresses (ipp);
01168 ipp++;
01169 }
01170 }
01171 }
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184 void ArgyrisTriangle::compute_nlstressincr (long lcid,long eid,long ri,long ci)
01185 {
01186 long i,ii,ipp;
01187
01188 for (ii=0;ii<nb;ii++){
01189 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01190
01191 for (i=0;i<intordsm[ii][ii];i++){
01192
01193 if (Mp->strcomp==1)
01194 Mm->computenlstressesincr (ipp);
01195 ipp++;
01196 }
01197 }
01198 }
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211 void ArgyrisTriangle::compute_eigstress(long lcid,long eid,long ri,long ci)
01212 {
01213 long i,ii,ipp;
01214 vector eigstr(tncomp),sig(tncomp);
01215 matrix d(tncomp,tncomp);
01216
01217 for (ii=0;ii<nb;ii++){
01218 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01219
01220 for (i=0;i<intordsm[ii][ii];i++){
01221 Mm->giveeigstrain (ipp,eigstr);
01222
01223 Mm->matstiff (d,ipp);
01224 mxv (d,eigstr,sig);
01225 Mm->storeeigstress (ipp,sig);
01226 ipp++;
01227 }
01228 }
01229 }
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 void ArgyrisTriangle::elem_integration (integratedquant iq,long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y)
01247 {
01248 long i,ii,ipp;
01249 double jac;
01250 ivector nodes(nne);
01251 vector w,gp1,gp2,areacoord(3);
01252 vector ipv,contr(ndofe);
01253 matrix gm;
01254
01255 Mt->give_elemnodes (eid,nodes);
01256 fillv (0.0,nv);
01257
01258 ii=0;
01259 {
01260 allocv (intordsm[ii][ii],gp1);
01261 allocv (intordsm[ii][ii],gp2);
01262 allocv (intordsm[ii][ii],w);
01263 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01264 ipp=Mt->elements[eid].ipp[ii][ii];
01265
01266 for (i=0;i<intordsm[ii][ii];i++){
01267 areacoord[0]=gp1[i];
01268 areacoord[1]=gp2[i];
01269 areacoord[2]=1.0-areacoord[0]-areacoord[2];
01270
01271 jac_2d (jac,x,y,areacoord[0],areacoord[1]);
01272 allocv (ncomp[ii],ipv);
01273
01274
01275 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
01276
01277
01278 allocm (ncomp[ii],ndofe,gm);
01279 geom_matrix (gm,x,y,areacoord);
01280
01281 mtxv (gm,ipv,contr);
01282 destrm (gm);
01283 cmulv (jac*w[i],contr);
01284
01285
01286 addv(contr,nv,nv);
01287 ipp++;
01288 destrv (ipv);
01289 }
01290 destrv (w); destrv (gp2); destrv (gp1);
01291 }
01292 }
01293
01294
01295
01296 void ArgyrisTriangle::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
01297 {
01298 long i, j, k, ipp;
01299 long ii, jj, nv = nodval.n;
01300 long nstra;
01301 double xi, eta, ipval;
01302 vector w, gp1, gp2, anv(nne);
01303
01304 nstra = 0;
01305 for (j = 0; j < nv; j++)
01306 {
01307 for(i = 0; i < nne; i++)
01308 anv[i] = nodval[i][j];
01309 for (ii = 0; ii < nb; ii++)
01310 {
01311 for (jj = 0; jj < nb; jj++)
01312 {
01313 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01314 if (intordsm[ii][jj] == 0)
01315 continue;
01316 allocv (intordsm[ii][jj],gp1);
01317 allocv (intordsm[ii][jj],gp2);
01318 allocv (intordsm[ii][jj],w);
01319 gauss_points_tr (gp1.a, gp2.a, w.a, intordsm[ii][jj]);
01320 for (k = 0; k < intordsm[ii][jj]; k++)
01321 {
01322 xi=gp1[k];
01323 eta=gp2[k];
01324
01325 ipval = approx_nat (xi,eta,anv);
01326 if ((ictn[i] & inistrain) && (j < Mm->ip[ipp].ncompstr))
01327 {
01328 Mm->ip[ipp].strain[j] += ipval;
01329 ipp++;
01330 continue;
01331 }
01332 if ((ictn[i] & inistress) && (j < nstra + Mm->ip[ipp].ncompstr))
01333 {
01334 Mm->ip[ipp].stress[j] += ipval;
01335 ipp++;
01336 continue;
01337 }
01338 if ((ictn[i] & iniother) && (j < nv))
01339 {
01340 Mm->ip[ipp].other[j] += ipval;
01341 ipp++;
01342 continue;
01343 }
01344 ipp++;
01345 }
01346 destrv(gp1); destrv (gp2); destrv (w);
01347 }
01348 }
01349 if (ictn[i] & inistrain) nstra++;
01350 }
01351 }
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363 inline double cube(double x)
01364 {
01365 return x * x * x;
01366 }
01367
01368
01369
01370
01371 void ArgyrisTriangle::computeTriangleConstants( double* xCoords, double* yCoords, double* A, double* L, double* c1, double* c2)
01372 {
01373
01374 *A = xCoords[0]*yCoords[1] + xCoords[1]*yCoords[2] + xCoords[2]*yCoords[0] - xCoords[0]*yCoords[2] - xCoords[1]*yCoords[0] - xCoords[2]*yCoords[1];
01375
01376
01377 for (int i = 0; i < 3; i++)
01378 {
01379 int j = (i + 1) % 3;
01380 int k = (i + 2) % 3;
01381 L[i] = sqrt(sqr(xCoords[j] - xCoords[k]) + sqr(yCoords[j] - yCoords[k]));
01382 }
01383
01384
01385 c1[0] = (yCoords[1] - yCoords[2]) / *A;
01386 c1[1] = (yCoords[2] - yCoords[0]) / *A;
01387 c1[2] = (yCoords[0] - yCoords[1]) / *A;
01388
01389 c2[0] = (xCoords[2] - xCoords[1]) / *A;
01390 c2[1] = (xCoords[0] - xCoords[2]) / *A;
01391 c2[2] = (xCoords[1] - xCoords[0]) / *A;
01392 }
01393
01394 double ArgyrisTriangle::getEij(int i, double* L)
01395 {
01396 int j = (i + 1) % 3;
01397 int k = (i + 2) % 3;
01398 return (sqr(L[k]) - sqr(L[i]) - sqr(L[j])) / (2 * sqr(L[j]));
01399 }
01400
01401 double ArgyrisTriangle::getEik(int i, double* L)
01402 {
01403 int j = (i + 1) % 3;
01404 int k = (i + 2) % 3;
01405 return (sqr(L[j]) - sqr(L[i]) - sqr(L[k])) / (2 * sqr(L[k]));
01406 }
01407
01408 void ArgyrisTriangle::computeAreaCoordinatesForPoint( double x, double y, double* xCoords, double* yCoords, double* omega)
01409 {
01410 double A = xCoords[0]*yCoords[1] + xCoords[1]*yCoords[2] + xCoords[2]*yCoords[0] - xCoords[0]*yCoords[2] - xCoords[1]*yCoords[0] - xCoords[2]*yCoords[1];
01411
01412 omega[0] = ((yCoords[1] - yCoords[2]) / A) * x + ((xCoords[2] - xCoords[1]) / A) * y + (xCoords[1] * yCoords[2] / A) - (xCoords[2] * yCoords[1] / A);
01413 omega[1] = ((yCoords[2] - yCoords[0]) / A) * x + ((xCoords[0] - xCoords[2]) / A) * y + (xCoords[2] * yCoords[0] / A) - (xCoords[0] * yCoords[2] / A);
01414 omega[2] = ((yCoords[0] - yCoords[1]) / A) * x + ((xCoords[1] - xCoords[0]) / A) * y + (xCoords[0] * yCoords[1] / A) - (xCoords[1] * yCoords[0] / A);
01415 }
01416
01417 void ArgyrisTriangle::convertAreaCoordinatesToNaturalCoordinates( double* xCoords, double* yCoords, double* omega, double* pointX, double* pointY)
01418 {
01419 double x = 0.0;
01420 double y = 0.0;
01421 for (int i = 0; i < 3; i++)
01422 {
01423 x += xCoords[i] * omega[i];
01424 y += yCoords[i] * omega[i];
01425 }
01426 *pointX = x;
01427 *pointY = y;
01428 }
01429
01430
01431
01432
01433
01434 double ArgyrisTriangle::get_Hi(int i, double* xCoords, double* yCoords, double pointX, double pointY)
01435 {
01436 double A;
01437 double L[3];
01438 double omega[3];
01439 double c1[3];
01440 double c2[3];
01441
01442 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01443 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01444
01445 int j = (i + 1) % 3;
01446 int k = (i + 2) % 3;
01447
01448 double result = cube(omega[i]) * (10.0 - 15.0 * omega[i] + 6.0 * sqr(omega[i])) - 30.0 * (getEij(i, L) * omega[k] + getEik(i, L) * omega[j]) * sqr(omega[i]) * omega[j] * omega[k];
01449
01450 return result;
01451 }
01452
01453 double ArgyrisTriangle::get_Hix(int i, double* xCoords, double* yCoords, double pointX, double pointY)
01454 {
01455 double A;
01456 double L[3];
01457 double omega[3];
01458 double c1[3];
01459 double c2[3];
01460
01461 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01462 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01463
01464 int j = (i + 1) % 3;
01465 int k = (i + 2) % 3;
01466
01467 double alpha = 7 * (xCoords[i] - xCoords[k]) * getEij(i, L) - 5 * (xCoords[j] - xCoords[k]);
01468 double beta = 7 * (xCoords[i] - xCoords[j]) * getEik(i, L) - 5 * (xCoords[k] - xCoords[j]);
01469
01470 double result = cube(omega[i]) * (4.0 - 3.0 * omega[i]) * (pointX - xCoords[i])
01471 + alpha * sqr(omega[i]) * omega[j] * sqr(omega[k])
01472 + beta * sqr(omega[i]) * sqr(omega[j]) * omega[k];
01473
01474 return result;
01475 }
01476
01477 double ArgyrisTriangle::get_Hiy(int i, double* xCoords, double* yCoords, double pointX, double pointY)
01478 {
01479 double A;
01480 double L[3];
01481 double omega[3];
01482 double c1[3];
01483 double c2[3];
01484
01485 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01486 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01487
01488 int j = (i + 1) % 3;
01489 int k = (i + 2) % 3;
01490
01491 double gamma = 7 * (yCoords[i] - yCoords[k]) * getEij(i, L) - 5 * (yCoords[j] - yCoords[k]);
01492 double delta = 7 * (yCoords[i] - yCoords[j]) * getEik(i, L) - 5 * (yCoords[k] - yCoords[j]);
01493
01494 double result = cube(omega[i]) * (4.0 - 3.0 * omega[i]) * (pointY - yCoords[i])
01495 + gamma * sqr(omega[i]) * omega[j] * sqr(omega[k])
01496 + delta * sqr(omega[i]) * sqr(omega[j]) * omega[k];
01497
01498 return result;
01499 }
01500
01501 double ArgyrisTriangle::get_Hixx(int i, double* xCoords, double* yCoords, double pointX, double pointY)
01502 {
01503 double A;
01504 double L[3];
01505 double omega[3];
01506 double c1[3];
01507 double c2[3];
01508
01509 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01510 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01511
01512 int j = (i + 1) % 3;
01513 int k = (i + 2) % 3;
01514
01515 double epsilon = -(xCoords[i] - xCoords[k]) * ((xCoords[i] - xCoords[k]) * getEij(i, L) - 2 * (xCoords[j] - xCoords[k]));
01516 double phi = -(xCoords[i] - xCoords[j]) * ((xCoords[i] - xCoords[j]) * getEik(i, L) - 2 * (xCoords[k] - xCoords[j]));
01517
01518 double result = 0.5 * ( cube(omega[i])*sqr(pointX-xCoords[i]) + epsilon*sqr(omega[i])*omega[j]*sqr(omega[k]) + phi*sqr(omega[i])*sqr(omega[j])*omega[k] );
01519
01520 return result;
01521 }
01522
01523 double ArgyrisTriangle::get_Hiyy(int i, double* xCoords, double* yCoords, double pointX, double pointY)
01524 {
01525 double A;
01526 double L[3];
01527 double omega[3];
01528 double c1[3];
01529 double c2[3];
01530
01531 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01532 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01533
01534 int j = (i + 1) % 3;
01535 int k = (i + 2) % 3;
01536
01537 double psi = -(yCoords[i] - yCoords[k]) * ((yCoords[i] - yCoords[k]) * getEij(i, L) - 2 * (yCoords[j] - yCoords[k]));
01538 double rho = -(yCoords[i] - yCoords[j]) * ((yCoords[i] - yCoords[j]) * getEik(i, L) - 2 * (yCoords[k] - yCoords[j]));
01539
01540 double result = 0.5 * ( cube(omega[i])*sqr(pointY-yCoords[i]) + psi*sqr(omega[i])*omega[j]*sqr(omega[k]) + rho*sqr(omega[i])*sqr(omega[j])*omega[k] );
01541
01542 return result;
01543 }
01544
01545 double ArgyrisTriangle::get_Hixy(int i, double* xCoords, double* yCoords, double pointX, double pointY)
01546 {
01547 double A;
01548 double L[3];
01549 double omega[3];
01550 double c1[3];
01551 double c2[3];
01552
01553 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01554 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01555
01556 int j = (i + 1) % 3;
01557 int k = (i + 2) % 3;
01558
01559 double mju = -0.5 * ((yCoords[i]-yCoords[k])*((xCoords[i]-xCoords[k])*getEij(i, L) - 2*(xCoords[j]-xCoords[k])) + (xCoords[i]-xCoords[k])*((yCoords[i]-yCoords[k])*getEij(i, L) - 2*(yCoords[j]-yCoords[k])));
01560 double pi = -0.5 * ((yCoords[i]-yCoords[j])*((xCoords[i]-xCoords[j])*getEik(i, L) - 2*(xCoords[k]-xCoords[j])) + (xCoords[i]-xCoords[j])*((yCoords[i]-yCoords[j])*getEik(i, L) - 2*(yCoords[k]-yCoords[j])));
01561
01562 double result = cube(omega[i])*(pointX-xCoords[i])*(pointY-yCoords[i]) + mju*sqr(omega[i])*omega[j]*sqr(omega[k]) + pi*sqr(omega[i])*sqr(omega[j])*omega[k];
01563
01564 return result;
01565 }
01566
01567 double ArgyrisTriangle::get_HjkNi(int i, double* xCoords, double* yCoords, double pointX, double pointY)
01568 {
01569 double A;
01570 double L[3];
01571 double omega[3];
01572 double c1[3];
01573 double c2[3];
01574
01575 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01576 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01577
01578 int j = (i + 1) % 3;
01579 int k = (i + 2) % 3;
01580
01581 double result = (-16.0 * A / L[i]) * omega[i]*sqr(omega[j])*sqr(omega[k]);
01582
01583 return result;
01584 }
01585
01586
01587
01588
01589
01590 double ArgyrisTriangle::get_Hi(int i, Derivatives dType, double* xCoords, double* yCoords, double pointX, double pointY)
01591 {
01592 double A;
01593 double L[3];
01594 double omega[3];
01595 double c1[3];
01596 double c2[3];
01597
01598 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01599 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01600
01601 double result = -60 * getEij(i, L) * getFirstCommonExpressionValue(i, dType, omega, c1, c2) - 60 * getEik(i, L) * getSecondCommonExpressionValue(i, dType, omega, c1, c2);
01602 switch (dType)
01603 {
01604 case dX2:
01605 result += 60 * sqr(c1[i]) * (omega[i] - 3 * sqr(omega[i]) + 2 * cube(omega[i]));
01606 break;
01607 case dY2:
01608 result += 60 * sqr(c2[i]) * (omega[i] - 3 * sqr(omega[i]) + 2 * cube(omega[i]));
01609 break;
01610 case dXdY:
01611 result += 60 * c1[i] * c2[i] * (omega[i] - 3 * sqr(omega[i]) + 2 * cube(omega[i]));
01612 break;
01613 }
01614 return result;
01615 }
01616
01617 double ArgyrisTriangle::get_Hix(int i, Derivatives dType, double* xCoords, double* yCoords, double pointX, double pointY)
01618 {
01619 double A;
01620 double L[3];
01621 double omega[3];
01622 double c1[3];
01623 double c2[3];
01624
01625 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01626 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01627
01628 int j = (i + 1) % 3;
01629 int k = (i + 2) % 3;
01630 double a = c1[i];
01631 double f = c2[i];
01632
01633 double alpha = 7 * (xCoords[i] - xCoords[k]) * getEij(i, L) - 5 * (xCoords[j] - xCoords[k]);
01634 double beta = 7 * (xCoords[i] - xCoords[j]) * getEik(i, L) - 5 * (xCoords[k] - xCoords[j]);
01635 double result = 2 * alpha * getFirstCommonExpressionValue(i, dType, omega, c1, c2) + 2 * beta * getSecondCommonExpressionValue(i, dType, omega, c1, c2);
01636
01637 switch (dType)
01638 {
01639 case dX2:
01640 result += 24*sqr(a)*omega[i]*pointX + 24*a*sqr(omega[i]) - 36*sqr(a)*sqr(omega[i])*pointX - 24*a*cube(omega[i]) - 24*sqr(a)*xCoords[i]*omega[i] + 36*sqr(a)*xCoords[i]*sqr(omega[i]);
01641 break;
01642 case dY2:
01643 result += 12*sqr(f) * (2*pointX*omega[i] - 3*pointX*sqr(omega[i]) - 2*xCoords[i]*omega[i] + 3*xCoords[i]*sqr(omega[i]));
01644 break;
01645 case dXdY:
01646 result += 24*a*f*omega[i]*pointX + 12*f*sqr(omega[i]) - 36*a*f*sqr(omega[i])*pointX - 12*f*cube(omega[i]) - 24*a*f*xCoords[i]*omega[i] + 36*a*f*xCoords[i]*sqr(omega[i]);
01647 break;
01648 }
01649 return result;
01650 }
01651
01652 double ArgyrisTriangle::get_Hiy(int i, Derivatives dType, double* xCoords, double* yCoords, double pointX, double pointY)
01653 {
01654 double A;
01655 double L[3];
01656 double omega[3];
01657 double c1[3];
01658 double c2[3];
01659
01660 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01661 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01662
01663 int j = (i + 1) % 3;
01664 int k = (i + 2) % 3;
01665 double a = c1[i];
01666 double f = c2[i];
01667
01668 double gamma = 7 * (yCoords[i] - yCoords[k]) * getEij(i, L) - 5 * (yCoords[j] - yCoords[k]);
01669 double delta = 7 * (yCoords[i] - yCoords[j]) * getEik(i, L) - 5 * (yCoords[k] - yCoords[j]);
01670 double result = 2 * gamma * getFirstCommonExpressionValue(i, dType, omega, c1, c2) + 2 * delta * getSecondCommonExpressionValue(i, dType, omega, c1, c2);
01671
01672 switch (dType)
01673 {
01674 case dX2:
01675 result += 24*sqr(a)*omega[i]*pointY - 36*sqr(a)*sqr(omega[i])*pointY - 24*sqr(a)*yCoords[i]*omega[i] + 36*sqr(a)*yCoords[i]*sqr(omega[i]);
01676 break;
01677 case dY2:
01678 result += 24*sqr(f)*omega[i]*pointY + 24*f*sqr(omega[i]) - 36*sqr(f)*sqr(omega[i])*pointY - 24*f*cube(omega[i]) - 24*sqr(f)*yCoords[i]*omega[i] + 36*sqr(f)*yCoords[i]*sqr(omega[i]);
01679 break;
01680 case dXdY:
01681 result += 24*a*f*omega[i]*pointY + 12*a*sqr(omega[i]) - 36*a*f*sqr(omega[i])*pointY - 12*a*cube(omega[i]) - 24*a*f*yCoords[i]*omega[i] + 36*a*f*yCoords[i]*sqr(omega[i]);
01682 break;
01683 }
01684 return result;
01685 }
01686
01687 double ArgyrisTriangle::get_Hixx(int i, Derivatives dType, double* xCoords, double* yCoords, double pointX, double pointY)
01688 {
01689 double A;
01690 double L[3];
01691 double omega[3];
01692 double c1[3];
01693 double c2[3];
01694
01695 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01696 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01697
01698 int j = (i + 1) % 3;
01699 int k = (i + 2) % 3;
01700 double a = c1[i];
01701 double f = c2[i];
01702
01703 double epsilon = -(xCoords[i] - xCoords[k]) * ((xCoords[i] - xCoords[k]) * getEij(i, L) - 2 * (xCoords[j] - xCoords[k]));
01704 double phi = -(xCoords[i] - xCoords[j]) * ((xCoords[i] - xCoords[j]) * getEik(i, L) - 2 * (xCoords[k] - xCoords[j]));
01705 double result = 2 * epsilon * getFirstCommonExpressionValue(i, dType, omega, c1, c2) + 2 * phi * getSecondCommonExpressionValue(i, dType, omega, c1, c2);
01706
01707 switch (dType)
01708 {
01709 case dX2:
01710 result += 6*sqr(a)*omega[i]*sqr(pointX-xCoords[i]) + 6*a*sqr(omega[i])*(pointX-xCoords[i]) + 6*sqr(omega[i])*a*(pointX-xCoords[i]) + 2*cube(omega[i]);
01711 break;
01712 case dY2:
01713 result += 6*sqr(f)*omega[i]*sqr(pointX-xCoords[i]);
01714 break;
01715 case dXdY:
01716 result += 6*a*f*omega[i]*sqr(pointX-xCoords[i]) + 6*f*sqr(omega[i])*(pointX-xCoords[i]);
01717 break;
01718 }
01719 result *= 0.5;
01720 return result;
01721 }
01722
01723 double ArgyrisTriangle::get_Hiyy(int i, Derivatives dType, double* xCoords, double* yCoords, double pointX, double pointY)
01724 {
01725 double A;
01726 double L[3];
01727 double omega[3];
01728 double c1[3];
01729 double c2[3];
01730
01731 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01732 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01733
01734 int j = (i + 1) % 3;
01735 int k = (i + 2) % 3;
01736 double a = c1[i];
01737 double f = c2[i];
01738
01739 double psi = -(yCoords[i] - yCoords[k]) * ((yCoords[i] - yCoords[k]) * getEij(i, L) - 2 * (yCoords[j] - yCoords[k]));
01740 double rho = -(yCoords[i] - yCoords[j]) * ((yCoords[i] - yCoords[j]) * getEik(i, L) - 2 * (yCoords[k] - yCoords[j]));
01741 double result = 2 * psi * getFirstCommonExpressionValue(i, dType, omega, c1, c2) + 2 * rho * getSecondCommonExpressionValue(i, dType, omega, c1, c2);
01742
01743 switch (dType)
01744 {
01745 case dX2:
01746 result += 6*sqr(a)*omega[i]*sqr(pointY-yCoords[i]);
01747 break;
01748 case dY2:
01749 result += 6*sqr(f)*omega[i]*sqr(pointY-yCoords[i]) + 6*f*sqr(omega[i])*(pointY-yCoords[i]) + 6*f*sqr(omega[i])*(pointY-yCoords[i]) + 2*cube(omega[i]);
01750 break;
01751 case dXdY:
01752 result += 6*a*f*omega[i]*sqr(pointY-yCoords[i]) + 6*a*sqr(omega[i])*(pointY-yCoords[i]);
01753 break;
01754 }
01755 result *= 0.5;
01756 return result;
01757 }
01758
01759 double ArgyrisTriangle::get_Hixy(int i, Derivatives dType, double* xCoords, double* yCoords, double pointX, double pointY)
01760 {
01761 double A;
01762 double L[3];
01763 double omega[3];
01764 double c1[3];
01765 double c2[3];
01766
01767 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01768 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01769
01770 int j = (i + 1) % 3;
01771 int k = (i + 2) % 3;
01772 double a = c1[i];
01773 double f = c2[i];
01774
01775 double mju = -0.5 * ((yCoords[i]-yCoords[k])*((xCoords[i]-xCoords[k])*getEij(i, L) - 2*(xCoords[j]-xCoords[k])) + (xCoords[i]-xCoords[k])*((yCoords[i]-yCoords[k])*getEij(i, L) - 2*(yCoords[j]-yCoords[k])));
01776 double pi = -0.5 * ((yCoords[i]-yCoords[j])*((xCoords[i]-xCoords[j])*getEik(i, L) - 2*(xCoords[k]-xCoords[j])) + (xCoords[i]-xCoords[j])*((yCoords[i]-yCoords[j])*getEik(i, L) - 2*(yCoords[k]-yCoords[j])));
01777 double result = 2 * mju * getFirstCommonExpressionValue(i, dType, omega, c1, c2) + 2 * pi * getSecondCommonExpressionValue(i, dType, omega, c1, c2);
01778
01779 switch (dType)
01780 {
01781 case dX2:
01782 result += 6*sqr(a)*omega[i]*(pointX-xCoords[i])*(pointY-yCoords[i]) + 6*a*sqr(omega[i])*(pointY-yCoords[i]);
01783 break;
01784 case dY2:
01785 result += 6*sqr(f)*omega[i]*(pointX-xCoords[i])*(pointY-yCoords[i]) + 6*f*sqr(omega[i])*(pointX-xCoords[i]);
01786 break;
01787 case dXdY:
01788 result += 6*a*f*omega[i]*(pointX-xCoords[i])*(pointY-yCoords[i]) + 3*a*sqr(omega[i])*(pointX-xCoords[i]) + 3*f*sqr(omega[i])*(pointY-yCoords[i]) + cube(omega[i]);
01789 break;
01790 }
01791
01792 return result;
01793 }
01794
01795 double ArgyrisTriangle::get_HjkNi(int i, Derivatives dType, double* xCoords, double* yCoords, double pointX, double pointY)
01796 {
01797 double A;
01798 double L[3];
01799 double omega[3];
01800 double c1[3];
01801 double c2[3];
01802
01803 computeAreaCoordinatesForPoint( pointX, pointY, xCoords, yCoords, omega);
01804 computeTriangleConstants( xCoords, yCoords, &A, L, c1, c2);
01805
01806 double theta = -32 * A / L[i];
01807 return theta * getThirdCommonExpressionValue(i, dType, omega, c1, c2);
01808 }
01809
01810 double ArgyrisTriangle::getFirstCommonExpressionValue(int i, Derivatives dType, double* omega, double* c1, double* c2)
01811 {
01812 int j = (i + 1) % 3;
01813 int k = (i + 2) % 3;
01814
01815 double* ca;
01816 double* cb;
01817
01818 switch (dType)
01819 {
01820 case dX2:
01821 ca = cb = c1;
01822 break;
01823 case dY2:
01824 ca = cb = c2;
01825 break;
01826 case dXdY:
01827 ca = c1;
01828 cb = c2;
01829 break;
01830 }
01831
01832 return ca[i] * (cb[i] * omega[j] * sqr(omega[k]) + cb[j] * omega[i] * sqr(omega[k]) + 2 * cb[k] * omega[i] * omega[j] * omega[k])
01833 + ca[j] * (cb[i] * omega[i] * sqr(omega[k]) + cb[k] * sqr(omega[i]) * omega[k])
01834 + ca[k] * (2 * cb[i] * omega[i] * omega[j] * omega[k] + cb[j] * sqr(omega[i]) * omega[k] + cb[k] * sqr(omega[i]) * omega[j]);
01835 }
01836
01837 double ArgyrisTriangle::getSecondCommonExpressionValue(int i, Derivatives dType, double* omega, double* c1, double* c2)
01838 {
01839 int j = (i + 1) % 3;
01840 int k = (i + 2) % 3;
01841
01842 double* ca;
01843 double* cb;
01844
01845 switch (dType)
01846 {
01847 case dX2:
01848 ca = cb = c1;
01849 break;
01850 case dY2:
01851 ca = cb = c2;
01852 break;
01853 case dXdY:
01854 ca = c1;
01855 cb = c2;
01856 break;
01857 }
01858
01859 return ca[i] * (cb[i] * sqr(omega[j]) * omega[k] + 2 * cb[j] * omega[i] * omega[j] * omega[k] + cb[k] * omega[i] * sqr(omega[j]))
01860 + ca[j] * (2 * cb[i] * omega[i] * omega[j] * omega[k] + cb[j] * sqr(omega[i]) * omega[k] + cb[k] * sqr(omega[i]) * omega[j])
01861 + ca[k] * (cb[i] * omega[i] * sqr(omega[j]) + cb[j] * sqr(omega[i]) * omega[j]);
01862 }
01863
01864 double ArgyrisTriangle::getThirdCommonExpressionValue(int i, Derivatives dType, double* omega, double* c1, double* c2)
01865 {
01866 int j = (i + 1) % 3;
01867 int k = (i + 2) % 3;
01868
01869 double* ca;
01870 double* cb;
01871
01872 switch (dType)
01873 {
01874 case dX2:
01875 ca = cb = c1;
01876 break;
01877 case dY2:
01878 ca = cb = c2;
01879 break;
01880 case dXdY:
01881 ca = c1;
01882 cb = c2;
01883 break;
01884 }
01885
01886 return ca[i] * (cb[j]*omega[j]*sqr(omega[k]) + cb[k]*sqr(omega[j])*omega[k])
01887 + ca[j] * (cb[i]*omega[j]*sqr(omega[k]) + cb[j]*omega[i]*sqr(omega[k]) + 2*cb[k]*omega[i]*omega[j]*omega[k])
01888 + ca[k] * (cb[i]*sqr(omega[j])*omega[k] + 2*cb[j]*omega[i]*omega[j]*omega[k] + cb[k]*omega[i]*sqr(omega[j]));
01889 }
01890