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