00001 #include <stdlib.h>
00002 #include <math.h>
00003 #include "linhexrot.h"
00004 #include "global.h"
00005 #include "globmat.h"
00006 #include "genfile.h"
00007 #include "intpoints.h"
00008 #include "node.h"
00009 #include "element.h"
00010 #include "loadcase.h"
00011
00012
00013
00014
00015 linhexrot::linhexrot ()
00016 {
00017 long i,j;
00018
00019
00020 nne=8;
00021
00022 ndofe=48;
00023
00024 tncomp=6;
00025
00026 napfun=3;
00027
00028 intordmm=2;
00029
00030 ned=12;
00031
00032 nned=2;
00033
00034 intordb=2;
00035
00036 nsurf=6;
00037
00038 nnsurf=4;
00039
00040 ssst=spacestress;
00041
00042
00043 nb=1;
00044
00045
00046 ncomp = new long [nb];
00047 ncomp[0]=6;
00048
00049
00050 cncomp = new long [nb];
00051 cncomp[0]=0;
00052
00053
00054
00055 nip = new long* [nb];
00056 intordsm = new long* [nb];
00057 for (i=0;i<nb;i++){
00058 nip[i] = new long [nb];
00059 intordsm[i] = new long [nb];
00060 }
00061
00062 nip[0][0]=8;
00063
00064
00065 tnip=0;
00066 for (i=0;i<nb;i++){
00067 for (j=0;j<nb;j++){
00068 tnip+=nip[i][j];
00069 }
00070 }
00071
00072 intordsm[0][0]=2;
00073 }
00074
00075 linhexrot::~linhexrot ()
00076 {
00077 long i;
00078
00079 for (i=0;i<nb;i++){
00080 delete [] nip[i];
00081 delete [] intordsm[i];
00082 }
00083 delete [] nip;
00084 delete [] intordsm;
00085
00086 delete [] cncomp;
00087 delete [] ncomp;
00088 }
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 double linhexrot::approx (double xi,double eta,double zeta,vector &nodval)
00101 {
00102 double f;
00103 vector bf(nne);
00104
00105 bf_lin_hex_3d (bf.a,xi,eta,zeta);
00106
00107 scprd (bf,nodval,f);
00108
00109 return f;
00110 }
00111
00112
00113
00114
00115
00116
00117
00118
00119 void linhexrot::tran_side (matrix &t12,matrix &t23,matrix &t34,matrix &t41, matrix &t15, matrix &t26,
00120 matrix &t37,matrix &t48, matrix &t56,matrix &t67,matrix &t78, matrix &t85, vector &x,vector &y,vector &z)
00121 {
00122 long i;
00123 double pom;
00124 vector an(3),am(3),dl(ned);
00125 matrix a(3,ned);
00126
00127 a[0][0]=x[1]-x[0]; a[0][1]=x[2]-x[1]; a[0][2] =x[3]-x[2]; a[0][3] =x[0]-x[3];
00128 a[0][4]=x[4]-x[0]; a[0][5]=x[5]-x[1]; a[0][6] =x[6]-x[2]; a[0][7] =x[7]-x[3];
00129 a[0][8]=x[5]-x[4]; a[0][9]=x[6]-x[5]; a[0][10]=x[7]-x[6]; a[0][11]=x[4]-x[7];
00130
00131 a[1][0]=y[1]-y[0]; a[1][1]=y[2]-y[1]; a[1][2] =y[3]-y[2]; a[1][3] =y[0]-y[3];
00132 a[1][4]=y[4]-y[0]; a[1][5]=y[5]-y[1]; a[1][6] =y[6]-y[2]; a[1][7] =y[7]-y[3];
00133 a[1][8]=y[5]-y[4]; a[1][9]=y[6]-y[5]; a[1][10]=y[7]-y[6]; a[1][11]=y[4]-y[7];
00134
00135 a[2][0]=z[1]-z[0]; a[2][1]=z[2]-z[1]; a[2][2] =z[3]-z[2]; a[2][3] =z[0]-z[3];
00136 a[2][4]=z[4]-z[0]; a[2][5]=z[5]-z[1]; a[2][6] =z[6]-z[2]; a[2][7] =z[7]-z[3];
00137 a[2][8]=z[5]-z[4]; a[2][9]=z[6]-z[5]; a[2][10]=z[7]-z[6]; a[2][11]=z[4]-z[7];
00138
00139
00140 for (i=0;i<ned;i++){
00141 dl[i]=sqrt(a[0][i]*a[0][i]+a[1][i]*a[1][i]+a[2][i]*a[2][i]);
00142 if(fabs(a[0][i])<=1.e-6 && fabs(a[1][i])<=1.e-6)
00143 {am[0]=1.; am[1]=0.; am[2]=0.;
00144 }
00145 else
00146 {pom=sqrt(a[0][i]*a[0][i]+a[1][i]*a[1][i]);
00147 am[0]=-a[1][i]/pom; am[1]= a[0][i]/pom; am[2]=0.;
00148 }
00149
00150 an[0]=-am[1]*a[2][i]/dl[i];
00151 an[1]= am[0]*a[2][i]/dl[i];
00152 an[2]= am[1]*a[0][i]/dl[i]-am[0]*a[1][i]/dl[i];
00153
00154 if (i==0){
00155 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.;
00156 t12[1][0]=-t12[0][1]; t12[1][1]=0.; t12[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00157 t12[2][0]=-t12[0][2]; t12[2][1]=-t12[1][2]; t12[2][2]=0.;
00158 }
00159 else if (i==1){
00160 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.;
00161 t23[1][0]=-t23[0][1]; t23[1][1]=0.; t23[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00162 t23[2][0]=-t23[0][2]; t23[2][1]=-t23[1][2]; t23[2][2]=0.;
00163 }
00164 else if (i==2){
00165 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.;
00166 t34[1][0]=-t34[0][1]; t34[1][1]=0.; t34[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00167 t34[2][0]=-t34[0][2]; t34[2][1]=-t34[1][2]; t34[2][2]=0.;
00168 }
00169 else if (i==3){
00170 t41[0][0]=0.; t41[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t41[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00171 t41[1][0]=-t41[0][1]; t41[1][1]=0.; t41[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00172 t41[2][0]=-t41[0][2]; t41[2][1]=-t41[1][2]; t41[2][2]=0.;
00173 }
00174 else if (i==4){
00175 t15[0][0]=0.; t15[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t15[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00176 t15[1][0]=-t15[0][1]; t15[1][1]=0.; t15[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00177 t15[2][0]=-t15[0][2]; t15[2][1]=-t15[1][2]; t15[2][2]=0.;
00178 }
00179 else if (i==5){
00180 t26[0][0]=0.; t26[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t26[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00181 t26[1][0]=-t26[0][1]; t26[1][1]=0.; t26[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00182 t26[2][0]=-t26[0][2]; t26[2][1]=-t26[1][2]; t26[2][2]=0.;
00183 }
00184 else if (i==6){
00185 t37[0][0]=0.; t37[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t37[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00186 t37[1][0]=-t37[0][1]; t37[1][1]=0.; t37[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00187 t37[2][0]=-t37[0][2]; t37[2][1]=-t37[1][2]; t37[2][2]=0.;
00188 }
00189 else if (i==7){
00190 t48[0][0]=0.; t48[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t48[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00191 t48[1][0]=-t48[0][1]; t48[1][1]=0.; t48[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00192 t48[2][0]=-t48[0][2]; t48[2][1]=-t48[1][2]; t48[2][2]=0.;
00193 }
00194 else if (i==8){
00195 t56[0][0]=0.; t56[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t56[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00196 t56[1][0]=-t56[0][1]; t56[1][1]=0.; t56[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00197 t56[2][0]=-t56[0][2]; t56[2][1]=-t56[1][2]; t56[2][2]=0.;
00198 }
00199 else if (i==9){
00200 t67[0][0]=0.; t67[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t67[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00201 t67[1][0]=-t67[0][1]; t67[1][1]=0.; t67[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00202 t67[2][0]=-t67[0][2]; t67[2][1]=-t67[1][2]; t67[2][2]=0.;
00203 }
00204 else if (i==10){
00205 t78[0][0]=0.; t78[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t78[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00206 t78[1][0]=-t78[0][1]; t78[1][1]=0.; t78[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00207 t78[2][0]=-t78[0][2]; t78[2][1]=-t78[1][2]; t78[2][2]=0.;
00208 }
00209 else if (i==11){
00210 t85[0][0]=0.; t85[0][1]=(an[0]*am[1]-am[0]*an[1])*dl[i]/2.; t85[0][2]=(an[0]*am[2]-am[0]*an[2])*dl[i]/2.;
00211 t85[1][0]=-t85[0][1]; t85[1][1]=0.; t85[1][2]=(an[1]*am[2]-am[1]*an[2])*dl[i]/2.;
00212 t85[2][0]=-t85[0][2]; t85[2][1]=-t85[1][2]; t85[2][2]=0.;
00213 }
00214
00215 }
00216
00217 }
00218
00219
00220
00221
00222
00223
00224
00225
00226 void linhexrot::bf_matrix (matrix &n,double xi,double eta,double zeta)
00227 {
00228 long i,j,k,l;
00229 vector bf(nne);
00230
00231 fillm (0.0,n);
00232
00233 bf_lin_hex_3d (bf.a,xi,eta,zeta);
00234
00235 j=0; k=1; l=2;
00236 for (i=0;i<nne;i++){
00237 n[0][j]=bf[i]; j+=3;
00238 n[1][k]=bf[i]; k+=3;
00239 n[2][l]=bf[i]; l+=3;
00240 }
00241 }
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 void linhexrot::geom_matrix (matrix &gm,vector &x,vector &y,vector &z,
00256 double xi,double eta,double zeta,double &jac)
00257 {
00258 long i,i1,i2,i3,i4,i5,i6,i7,i8;
00259 vector dx(nne),dy(nne),dz(nne), dxr(ned),dyr(ned),dzr(ned);
00260 matrix t12(3,3),t23(3,3),t34(3,3),t41(3,3), t15(3,3),t26(3,3),t37(3,3),t48(3,3), t56(3,3),t67(3,3),t78(3,3),t85(3,3);
00261
00262 dx_bf_lin_hex_3d (dx.a, eta,zeta);
00263 dy_bf_lin_hex_3d (dy.a, xi,zeta);
00264 dz_bf_lin_hex_3d (dz.a, xi,eta);
00265 derivatives_3d (dx,dy,dz,jac, x,y,z,xi,eta,zeta);
00266
00267
00268 dx_bf_quad_hexrot_3d (dxr.a, xi,eta,zeta);
00269 dy_bf_quad_hexrot_3d (dyr.a, xi,eta,zeta);
00270 dz_bf_quad_hexrot_3d (dzr.a, xi,eta,zeta);
00271 derivatives_3d (dxr,dyr,dzr,jac, x,y,z,xi,eta,zeta);
00272
00273 tran_side (t12,t23,t34,t41, t15,t26,t37,t48, t56,t67,t78,t85, x,y,z);
00274
00275 fillm (0.0,gm);
00276
00277 i1=0; i2=1; i3=2;
00278 for (i=0;i<nne;i++){
00279 gm[0][i1]=dx[i];
00280 gm[1][i2]=dy[i];
00281 gm[2][i3]=dz[i];
00282
00283 gm[3][i2]=dz[i];
00284 gm[3][i3]=dy[i];
00285 gm[4][i1]=dz[i];
00286 gm[4][i3]=dx[i];
00287 gm[5][i1]=dy[i];
00288 gm[5][i2]=dx[i];
00289 i1+=6; i2+=6; i3+=6;
00290 }
00291
00292 i1=3; i2=9; i3=15; i4=21; i5=27; i6=33; i7=39; i8=45;
00293 for (i=0;i<3;i++){
00294
00295 gm[0][i1]=-t12[0][i]*dxr[0]+t41[0][i]*dxr[3]-t15[0][i]*dxr[4];
00296 gm[1][i1]=-t12[1][i]*dyr[0]+t41[1][i]*dyr[3]-t15[1][i]*dyr[4];
00297 gm[2][i1]=-t12[2][i]*dzr[0]+t41[2][i]*dzr[3]-t15[2][i]*dzr[4];
00298
00299 gm[3][i1]=-t12[1][i]*dzr[0]+t41[1][i]*dzr[3]-t15[1][i]*dzr[4]
00300 -t12[2][i]*dyr[0]+t41[2][i]*dyr[3]-t15[2][i]*dyr[4];
00301 gm[4][i1]=-t12[2][i]*dxr[0]+t41[2][i]*dxr[3]-t15[2][i]*dxr[4]
00302 -t12[0][i]*dzr[0]+t41[0][i]*dzr[3]-t15[0][i]*dzr[4];
00303 gm[5][i1]=-t12[0][i]*dyr[0]+t41[0][i]*dyr[3]-t15[0][i]*dyr[4]
00304 -t12[1][i]*dxr[0]+t41[1][i]*dxr[3]-t15[1][i]*dxr[4];
00305
00306 gm[0][i2]=t12[0][i]*dxr[0]-t23[0][i]*dxr[1]-t26[0][i]*dxr[5];
00307 gm[1][i2]=t12[1][i]*dyr[0]-t23[1][i]*dyr[1]-t26[1][i]*dyr[5];
00308 gm[2][i2]=t12[2][i]*dzr[0]-t23[2][i]*dzr[1]-t26[2][i]*dzr[5];
00309 gm[3][i2]=t12[1][i]*dzr[0]-t23[1][i]*dzr[1]-t26[1][i]*dzr[5]
00310 +t12[2][i]*dyr[0]-t23[2][i]*dyr[1]-t26[2][i]*dyr[5];
00311 gm[4][i2]=t12[2][i]*dxr[0]-t23[2][i]*dxr[1]-t26[2][i]*dxr[5]
00312 +t12[0][i]*dzr[0]-t23[0][i]*dzr[1]-t26[0][i]*dzr[5];
00313 gm[5][i2]=t12[0][i]*dyr[0]-t23[0][i]*dyr[1]-t26[0][i]*dyr[5]
00314 +t12[1][i]*dxr[0]-t23[1][i]*dxr[1]-t26[1][i]*dxr[5];
00315
00316 gm[0][i3]=t23[0][i]*dxr[1]-t34[0][i]*dxr[2]-t37[0][i]*dxr[6];
00317 gm[1][i3]=t23[1][i]*dyr[1]-t34[1][i]*dyr[2]-t37[1][i]*dyr[6];
00318 gm[2][i3]=t23[2][i]*dzr[1]-t34[2][i]*dzr[2]-t37[2][i]*dzr[6];
00319 gm[3][i3]=t23[1][i]*dzr[1]-t34[1][i]*dzr[2]-t37[1][i]*dzr[6]
00320 +t23[2][i]*dyr[1]-t34[2][i]*dyr[2]-t37[2][i]*dyr[6];
00321 gm[4][i3]=t23[2][i]*dxr[1]-t34[2][i]*dxr[2]-t37[2][i]*dxr[6]
00322 +t23[0][i]*dzr[1]-t34[0][i]*dzr[2]-t37[0][i]*dzr[6];
00323 gm[5][i3]=t23[0][i]*dyr[1]-t34[0][i]*dyr[2]-t37[0][i]*dyr[6]
00324 +t23[1][i]*dxr[1]-t34[1][i]*dxr[2]-t37[1][i]*dxr[6];
00325
00326 gm[0][i4]=t34[0][i]*dxr[2]-t41[0][i]*dxr[3]-t48[0][i]*dxr[7];
00327 gm[1][i4]=t34[1][i]*dyr[2]-t41[1][i]*dyr[3]-t48[1][i]*dyr[7];
00328 gm[2][i4]=t34[2][i]*dzr[2]-t41[2][i]*dzr[3]-t48[2][i]*dzr[7];
00329 gm[3][i4]=t34[1][i]*dzr[2]-t41[1][i]*dzr[3]-t48[1][i]*dzr[7]
00330 +t34[2][i]*dyr[2]-t41[2][i]*dyr[3]-t48[2][i]*dyr[7];
00331 gm[4][i4]=t34[2][i]*dxr[2]-t41[2][i]*dxr[3]-t48[2][i]*dxr[7]
00332 +t34[0][i]*dzr[2]-t41[0][i]*dzr[3]-t48[0][i]*dzr[7];
00333 gm[5][i4]=t34[0][i]*dyr[2]-t41[0][i]*dyr[3]-t48[0][i]*dyr[7]
00334 +t34[1][i]*dxr[2]-t41[1][i]*dxr[3]-t48[1][i]*dxr[7];
00335
00336 gm[0][i5]=t15[0][i]*dxr[4]-t56[0][i]*dxr[8]+t85[0][i]*dxr[11];
00337 gm[1][i5]=t15[1][i]*dyr[4]-t56[1][i]*dyr[8]+t85[1][i]*dyr[11];
00338 gm[2][i5]=t15[2][i]*dzr[4]-t56[2][i]*dzr[8]+t85[2][i]*dzr[11];
00339 gm[3][i5]=t15[1][i]*dzr[4]-t56[1][i]*dzr[8]+t85[1][i]*dzr[11]
00340 +t15[2][i]*dyr[4]-t56[2][i]*dyr[8]+t85[2][i]*dyr[11];
00341 gm[4][i5]=t15[2][i]*dxr[4]-t56[2][i]*dxr[8]+t85[2][i]*dxr[11]
00342 +t15[0][i]*dzr[4]-t56[0][i]*dzr[8]+t85[0][i]*dzr[11];
00343 gm[5][i5]=t15[0][i]*dyr[4]-t56[0][i]*dyr[8]+t85[0][i]*dyr[11]
00344 +t15[1][i]*dxr[4]-t56[1][i]*dxr[8]+t85[1][i]*dxr[11];
00345
00346 gm[0][i6]=t26[0][i]*dxr[5]+t56[0][i]*dxr[8]-t67[0][i]*dxr[9];
00347 gm[1][i6]=t26[1][i]*dyr[5]+t56[1][i]*dyr[8]-t67[1][i]*dyr[9];
00348 gm[2][i6]=t26[2][i]*dzr[5]+t56[2][i]*dzr[8]-t67[2][i]*dzr[9];
00349 gm[3][i6]=t26[1][i]*dzr[5]+t56[1][i]*dzr[8]-t67[1][i]*dzr[9]
00350 +t26[2][i]*dyr[5]+t56[2][i]*dyr[8]-t67[2][i]*dyr[9];
00351 gm[4][i6]=t26[2][i]*dxr[5]+t56[2][i]*dxr[8]-t67[2][i]*dxr[9]
00352 +t26[0][i]*dzr[5]+t56[0][i]*dzr[8]-t67[0][i]*dzr[9];
00353 gm[5][i6]=t26[0][i]*dyr[5]+t56[0][i]*dyr[8]-t67[0][i]*dyr[9]
00354 +t26[1][i]*dxr[5]+t56[1][i]*dxr[8]-t67[1][i]*dxr[9];
00355
00356 gm[0][i7]=t37[0][i]*dxr[6]+t67[0][i]*dxr[9]-t78[0][i]*dxr[10];
00357 gm[1][i7]=t37[1][i]*dyr[6]+t67[1][i]*dyr[9]-t78[1][i]*dyr[10];
00358 gm[2][i7]=t37[2][i]*dzr[6]+t67[2][i]*dzr[9]-t78[2][i]*dzr[10];
00359 gm[3][i7]=t37[1][i]*dzr[6]+t67[1][i]*dzr[9]-t78[1][i]*dzr[10]
00360 +t37[2][i]*dyr[6]+t67[2][i]*dyr[9]-t78[2][i]*dyr[10];
00361 gm[4][i7]=t37[2][i]*dxr[6]+t67[2][i]*dxr[9]-t78[2][i]*dxr[10]
00362 +t37[0][i]*dzr[6]+t67[0][i]*dzr[9]-t78[0][i]*dzr[10];
00363 gm[5][i7]=t37[0][i]*dyr[6]+t67[0][i]*dyr[9]-t78[0][i]*dyr[10]
00364 +t37[1][i]*dxr[6]+t67[1][i]*dxr[9]-t78[1][i]*dxr[10];
00365
00366 gm[0][i8]=t48[0][i]*dxr[7]+t78[0][i]*dxr[10]-t85[0][i]*dxr[11];
00367 gm[1][i8]=t48[1][i]*dyr[7]+t78[1][i]*dyr[10]-t85[1][i]*dyr[11];
00368 gm[2][i8]=t48[2][i]*dzr[7]+t78[2][i]*dzr[10]-t85[2][i]*dzr[11];
00369 gm[3][i8]=t48[1][i]*dzr[7]+t78[1][i]*dzr[10]-t85[1][i]*dzr[11]
00370 +t48[2][i]*dyr[7]+t78[2][i]*dyr[10]-t85[2][i]*dyr[11];
00371 gm[4][i8]=t48[2][i]*dxr[7]+t78[2][i]*dxr[10]-t85[2][i]*dxr[11]
00372 +t48[0][i]*dzr[7]+t78[0][i]*dzr[10]-t85[0][i]*dzr[11];
00373 gm[5][i8]=t48[0][i]*dyr[7]+t78[0][i]*dyr[10]-t85[0][i]*dyr[11]
00374 +t48[1][i]*dxr[7]+t78[1][i]*dxr[10]-t85[1][i]*dxr[11];
00375 i1+=1; i2+=1; i3+=1; i4+=1; i5+=1; i6+=1; i7+=1; i8+=1;
00376 }
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 }
00389
00390 void linhexrot::geom_matrix_shear (matrix &gm,vector &x,vector &y,vector &z,
00391 double xi,double eta,double zeta,double &jac)
00392 {
00393 long i,i1,i2,i3,i4,i5,i6,i7,i8;
00394 vector dx(nne),dy(nne),dz(nne), dxr(ned),dyr(ned),dzr(ned), bf(nne);
00395 matrix t12(3,3),t23(3,3),t34(3,3),t41(3,3), t15(3,3),t26(3,3),t37(3,3),t48(3,3), t56(3,3),t67(3,3),t78(3,3),t85(3,3);
00396
00397 dx_bf_lin_hex_3d (dx.a,eta,zeta);
00398 dy_bf_lin_hex_3d (dy.a,xi,zeta);
00399 dz_bf_lin_hex_3d (dz.a,xi,eta);
00400 derivatives_3d (dx,dy,dz,jac,x,y,z,xi,eta,zeta);
00401
00402
00403 dx_bf_quad_hexrot_3d (dxr.a, xi,eta,zeta);
00404 dy_bf_quad_hexrot_3d (dyr.a, xi,eta,zeta);
00405 dz_bf_quad_hexrot_3d (dzr.a, xi,eta,zeta);
00406 derivatives_3d (dxr,dyr,dzr,jac, x,y,z,xi,eta,zeta);
00407
00408 tran_side (t12,t23,t34,t41, t15,t26,t37,t48, t56,t67,t78,t85, x,y,z);
00409
00410 fillm (0.0,gm);
00411
00412
00413 i1=3; i2=9; i3=15; i4=21; i5=27; i6=33; i7=39; i8=45;
00414 for (i=0;i<3;i++){
00415
00416 gm[0][i1]=(t12[1][i]*dzr[0]-t41[1][i]*dzr[3]+t15[1][i]*dzr[4]
00417 -t12[2][i]*dyr[0]+t41[2][i]*dyr[3]-t15[2][i]*dyr[4])/2.;
00418 gm[1][i1]=(t12[2][i]*dxr[0]-t41[2][i]*dxr[3]+t15[2][i]*dxr[4]
00419 -t12[0][i]*dzr[0]+t41[0][i]*dzr[3]-t15[0][i]*dzr[4])/2.;
00420 gm[2][i1]=(t12[0][i]*dyr[0]-t41[0][i]*dyr[3]+t15[0][i]*dyr[4]
00421 -t12[1][i]*dxr[0]+t41[1][i]*dxr[3]-t15[1][i]*dxr[4])/2.;
00422
00423 gm[0][i2]=(-t12[1][i]*dzr[0]+t23[1][i]*dzr[1]+t26[1][i]*dzr[5]
00424 +t12[2][i]*dyr[0]-t23[2][i]*dyr[1]-t26[2][i]*dyr[5])/2.;
00425 gm[1][i2]=(-t12[2][i]*dxr[0]+t23[2][i]*dxr[1]+t26[2][i]*dxr[5]
00426 +t12[0][i]*dzr[0]-t23[0][i]*dzr[1]-t26[0][i]*dzr[5])/2.;
00427 gm[2][i2]=(-t12[0][i]*dyr[0]+t23[0][i]*dyr[1]+t26[0][i]*dyr[5]
00428 +t12[1][i]*dxr[0]-t23[1][i]*dxr[1]-t26[1][i]*dxr[5])/2.;
00429
00430 gm[0][i3]=(-t23[1][i]*dzr[1]+t34[1][i]*dzr[2]+t37[1][i]*dzr[6]
00431 +t23[2][i]*dyr[1]-t34[2][i]*dyr[2]-t37[2][i]*dyr[6])/2.;
00432 gm[1][i3]=(-t23[2][i]*dxr[1]+t34[2][i]*dxr[2]+t37[2][i]*dxr[6]
00433 +t23[0][i]*dzr[1]-t34[0][i]*dzr[2]-t37[0][i]*dzr[6])/2.;
00434 gm[2][i3]=(-t23[0][i]*dyr[1]+t34[0][i]*dyr[2]+t37[0][i]*dyr[6]
00435 +t23[1][i]*dxr[1]-t34[1][i]*dxr[2]-t37[1][i]*dxr[6])/2.;
00436
00437 gm[0][i4]=(-t34[1][i]*dzr[2]+t41[1][i]*dzr[3]+t48[1][i]*dzr[7]
00438 +t34[2][i]*dyr[2]-t41[2][i]*dyr[3]-t48[2][i]*dyr[7])/2.;
00439 gm[1][i4]=(-t34[2][i]*dxr[2]+t41[2][i]*dxr[3]+t48[2][i]*dxr[7]
00440 +t34[0][i]*dzr[2]-t41[0][i]*dzr[3]-t48[0][i]*dzr[7])/2.;
00441 gm[2][i4]=(-t34[0][i]*dyr[2]+t41[0][i]*dyr[3]+t48[0][i]*dyr[7]
00442 +t34[1][i]*dxr[2]-t41[1][i]*dxr[3]-t48[1][i]*dxr[7])/2.;
00443
00444 gm[0][i5]=(-t15[1][i]*dzr[4]+t56[1][i]*dzr[8]-t85[1][i]*dzr[11]
00445 +t15[2][i]*dyr[4]-t56[2][i]*dyr[8]+t85[2][i]*dyr[11])/2.;
00446 gm[1][i5]=(-t15[2][i]*dxr[4]+t56[2][i]*dxr[8]-t85[2][i]*dxr[11]
00447 +t15[0][i]*dzr[4]-t56[0][i]*dzr[8]+t85[0][i]*dzr[11])/2.;
00448 gm[2][i5]=(-t15[0][i]*dyr[4]+t56[0][i]*dyr[8]-t85[0][i]*dyr[11]
00449 +t15[1][i]*dxr[4]-t56[1][i]*dxr[8]+t85[1][i]*dxr[11])/2.;
00450
00451 gm[0][i6]=(-t26[1][i]*dzr[5]-t56[1][i]*dzr[8]+t67[1][i]*dzr[9]
00452 +t26[2][i]*dyr[5]+t56[2][i]*dyr[8]-t67[2][i]*dyr[9])/2.;
00453 gm[1][i6]=(-t26[2][i]*dxr[5]-t56[2][i]*dxr[8]+t67[2][i]*dxr[9]
00454 +t26[0][i]*dzr[5]+t56[0][i]*dzr[8]-t67[0][i]*dzr[9])/2.;
00455 gm[2][i6]=(-t26[0][i]*dyr[5]-t56[0][i]*dyr[8]+t67[0][i]*dyr[9]
00456 +t26[1][i]*dxr[5]+t56[1][i]*dxr[8]-t67[1][i]*dxr[9])/2.;
00457
00458 gm[0][i7]=(-t37[1][i]*dzr[6]-t67[1][i]*dzr[9]+t78[1][i]*dzr[10]
00459 +t37[2][i]*dyr[6]+t67[2][i]*dyr[9]-t78[2][i]*dyr[10])/2.;
00460 gm[1][i7]=(-t37[2][i]*dxr[6]-t67[2][i]*dxr[9]+t78[2][i]*dxr[10]
00461 +t37[0][i]*dzr[6]+t67[0][i]*dzr[9]-t78[0][i]*dzr[10])/2.;
00462 gm[2][i7]=(-t37[0][i]*dyr[6]-t67[0][i]*dyr[9]+t78[0][i]*dyr[10]
00463 +t37[1][i]*dxr[6]+t67[1][i]*dxr[9]-t78[1][i]*dxr[10])/2.;
00464
00465 gm[0][i8]=(-t48[1][i]*dzr[7]-t78[1][i]*dzr[10]+t85[1][i]*dzr[11]
00466 +t48[2][i]*dyr[7]+t78[2][i]*dyr[10]-t85[2][i]*dyr[11])/2.;
00467 gm[1][i8]=(-t48[2][i]*dxr[7]-t78[2][i]*dxr[10]+t85[2][i]*dxr[11]
00468 +t48[0][i]*dzr[7]+t78[0][i]*dzr[10]-t85[0][i]*dzr[11])/2.;
00469 gm[2][i8]=(-t48[0][i]*dyr[7]-t78[0][i]*dyr[10]+t85[0][i]*dyr[11]
00470 +t48[1][i]*dxr[7]+t78[1][i]*dxr[10]-t85[1][i]*dxr[11])/2.;
00471 i1+=1; i2+=1; i3+=1; i4+=1; i5+=1; i6+=1; i7+=1; i8+=1;
00472 }
00473
00474
00475 bf_lin_hex_3d (bf.a,xi,eta,zeta);
00476 i1=0; i2=1; i3=2;
00477 for (i=0;i<nne;i++){
00478
00479
00480 gm[0][i1+3]-= bf.a[i];
00481
00482
00483
00484 gm[1][i1+4]-= bf.a[i];
00485
00486
00487
00488 gm[2][i1+5]-= bf.a[i];
00489 i1+=6; i2+=6; i3+=6;
00490 }
00491
00492
00493
00494 }
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507 void linhexrot::bvectors (vector &x,vector &y,vector &z,double xi,double eta,double zeta,double &jac,
00508 vector &b11,vector &b12,vector &b13,
00509 vector &b21,vector &b22,vector &b23,
00510 vector &b31,vector &b32,vector &b33)
00511 {
00512 vector dx(nne),dy(nne),dz(nne);
00513
00514 dx_bf_lin_hex_3d (dx.a,eta,zeta);
00515 dy_bf_lin_hex_3d (dy.a,xi,zeta);
00516 dz_bf_lin_hex_3d (dz.a,xi,eta);
00517
00518 derivatives_3d (dx,dy,dz,jac,x,y,z,xi,eta,zeta);
00519
00520 fillv (0.0,b11);
00521 fillv (0.0,b12);
00522 fillv (0.0,b13);
00523 fillv (0.0,b21);
00524 fillv (0.0,b22);
00525 fillv (0.0,b23);
00526 fillv (0.0,b31);
00527 fillv (0.0,b32);
00528 fillv (0.0,b33);
00529
00530
00531 b11[0]=dx[0]; b11[3]=dx[1]; b11[6]=dx[2]; b11[9]=dx[3]; b11[12]=dx[4]; b11[15]=dx[5]; b11[18]=dx[6]; b11[21]=dx[7];
00532
00533 b12[0]=dy[0]; b12[3]=dy[1]; b12[6]=dy[2]; b12[9]=dy[3]; b12[12]=dy[4]; b12[15]=dy[5]; b12[18]=dy[6]; b12[21]=dy[7];
00534
00535 b13[0]=dz[0]; b13[3]=dz[1]; b13[6]=dz[2]; b13[9]=dz[3]; b13[12]=dz[4]; b13[15]=dz[5]; b13[18]=dz[6]; b13[21]=dz[7];
00536
00537
00538 b21[1]=dx[0]; b21[4]=dx[1]; b21[7]=dx[2]; b21[10]=dx[3]; b21[13]=dx[4]; b21[16]=dx[5]; b21[19]=dx[6]; b21[22]=dx[7];
00539
00540 b22[1]=dy[0]; b22[4]=dy[1]; b22[7]=dy[2]; b22[10]=dy[3]; b22[13]=dy[4]; b22[16]=dy[5]; b22[19]=dy[6]; b22[22]=dy[7];
00541
00542 b23[1]=dz[0]; b23[4]=dz[1]; b23[7]=dz[2]; b23[10]=dz[3]; b23[13]=dz[4]; b23[16]=dz[5]; b23[19]=dz[6]; b23[22]=dz[7];
00543
00544
00545 b31[2]=dx[0]; b31[5]=dx[1]; b31[8]=dx[2]; b31[11]=dx[3]; b31[14]=dx[4]; b31[17]=dx[5]; b31[20]=dx[6]; b31[23]=dx[7];
00546
00547 b32[2]=dy[0]; b32[5]=dy[1]; b32[8]=dy[2]; b32[11]=dy[3]; b32[14]=dy[4]; b32[17]=dy[5]; b32[20]=dy[6]; b32[23]=dy[7];
00548
00549 b33[2]=dz[0]; b33[5]=dz[1]; b33[8]=dz[2]; b33[11]=dz[3]; b33[14]=dz[4]; b33[17]=dz[5]; b33[20]=dz[6]; b33[23]=dz[7];
00550
00551 }
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564 void linhexrot::gngeom_matrix (matrix &gm,vector &r,vector &x,vector &y,vector &z,double xi,double eta,double zeta,double &jac)
00565 {
00566 long i;
00567 double b11r,b12r,b13r,b21r,b22r,b23r,b31r,b32r,b33r;
00568 vector b11(ndofe),b12(ndofe),b13(ndofe),b21(ndofe),b22(ndofe),b23(ndofe),b31(ndofe),b32(ndofe),b33(ndofe),av(ndofe);
00569
00570 fillm (0.0,gm);
00571
00572 bvectors (x,y,z,xi,eta,zeta,jac,b11,b12,b13,b21,b22,b23,b31,b32,b33);
00573
00574 scprd (b11,r,b11r);
00575 scprd (b12,r,b12r);
00576 scprd (b13,r,b13r);
00577 scprd (b21,r,b21r);
00578 scprd (b22,r,b22r);
00579 scprd (b23,r,b23r);
00580 scprd (b31,r,b31r);
00581 scprd (b32,r,b32r);
00582 scprd (b33,r,b33r);
00583
00584
00585
00586
00587
00588
00589
00590 for (i=0;i<ndofe;i++){
00591 gm[0][i]+=b11[i];
00592 }
00593
00594
00595 cmulv(b11r,b11,av);
00596 for (i=0;i<ndofe;i++){
00597 gm[0][i]+=av[i];
00598 }
00599
00600
00601 cmulv(b21r,b21,av);
00602 for (i=0;i<ndofe;i++){
00603 gm[0][i]+=av[i];
00604 }
00605
00606
00607 cmulv(b31r,b31,av);
00608 for (i=0;i<ndofe;i++){
00609 gm[0][i]+=av[i];
00610 }
00611
00612
00613
00614
00615
00616
00617
00618 for (i=0;i<ndofe;i++){
00619 gm[1][i]+=b22[i];
00620 }
00621
00622
00623 cmulv(b12r,b12,av);
00624 for (i=0;i<ndofe;i++){
00625 gm[1][i]+=av[i];
00626 }
00627
00628
00629 cmulv(b22r,b22,av);
00630 for (i=0;i<ndofe;i++){
00631 gm[1][i]+=av[i];
00632 }
00633
00634
00635 cmulv(b32r,b32,av);
00636 for (i=0;i<ndofe;i++){
00637 gm[1][i]+=av[i];
00638 }
00639
00640
00641
00642
00643
00644
00645 for (i=0;i<ndofe;i++){
00646 gm[2][i]+=b22[i];
00647 }
00648
00649
00650 cmulv(b13r,b13,av);
00651 for (i=0;i<ndofe;i++){
00652 gm[2][i]+=av[i];
00653 }
00654
00655
00656 cmulv(b23r,b23,av);
00657 for (i=0;i<ndofe;i++){
00658 gm[2][i]+=av[i];
00659 }
00660
00661
00662 cmulv(b33r,b33,av);
00663 for (i=0;i<ndofe;i++){
00664 gm[2][i]+=av[i];
00665 }
00666
00667
00668
00669
00670
00671
00672
00673 for (i=0;i<ndofe;i++){
00674 gm[3][i]+=b23[i]+b32[i];
00675 }
00676
00677
00678 cmulv(b13r,b12,av);
00679 for (i=0;i<ndofe;i++){
00680 gm[3][i]+=av[i];
00681 }
00682
00683
00684 cmulv(b12r,b13,av);
00685 for (i=0;i<ndofe;i++){
00686 gm[3][i]+=av[i];
00687 }
00688
00689
00690 cmulv(b23r,b22,av);
00691 for (i=0;i<ndofe;i++){
00692 gm[3][i]+=av[i];
00693 }
00694
00695
00696 cmulv(b22r,b23,av);
00697 for (i=0;i<ndofe;i++){
00698 gm[3][i]+=av[i];
00699 }
00700
00701
00702 cmulv(b33r,b32,av);
00703 for (i=0;i<ndofe;i++){
00704 gm[3][i]+=av[i];
00705 }
00706
00707
00708 cmulv(b32r,b33,av);
00709 for (i=0;i<ndofe;i++){
00710 gm[3][i]+=av[i];
00711 }
00712
00713
00714
00715
00716
00717
00718 for (i=0;i<ndofe;i++){
00719 gm[4][i]+=b31[i]+b13[i];
00720 }
00721
00722
00723 cmulv(b11r,b13,av);
00724 for (i=0;i<ndofe;i++){
00725 gm[4][i]+=av[i];
00726 }
00727
00728
00729 cmulv(b13r,b11,av);
00730 for (i=0;i<ndofe;i++){
00731 gm[4][i]+=av[i];
00732 }
00733
00734
00735 cmulv(b21r,b23,av);
00736 for (i=0;i<ndofe;i++){
00737 gm[4][i]+=av[i];
00738 }
00739
00740
00741 cmulv(b23r,b21,av);
00742 for (i=0;i<ndofe;i++){
00743 gm[4][i]+=av[i];
00744 }
00745
00746
00747 cmulv(b31r,b33,av);
00748 for (i=0;i<ndofe;i++){
00749 gm[4][i]+=av[i];
00750 }
00751
00752
00753 cmulv(b33r,b31,av);
00754 for (i=0;i<ndofe;i++){
00755 gm[4][i]+=av[i];
00756 }
00757
00758
00759
00760
00761
00762
00763
00764 for (i=0;i<ndofe;i++){
00765 gm[5][i]+=b12[i]+b21[i];
00766 }
00767
00768
00769 cmulv(b12r,b11,av);
00770 for (i=0;i<ndofe;i++){
00771 gm[5][i]+=av[i];
00772 }
00773
00774
00775 cmulv(b11r,b12,av);
00776 for (i=0;i<ndofe;i++){
00777 gm[5][i]+=av[i];
00778 }
00779
00780
00781 cmulv(b22r,b21,av);
00782 for (i=0;i<ndofe;i++){
00783 gm[5][i]+=av[i];
00784 }
00785
00786
00787 cmulv(b21r,b22,av);
00788 for (i=0;i<ndofe;i++){
00789 gm[5][i]+=av[i];
00790 }
00791
00792
00793 cmulv(b32r,b31,av);
00794 for (i=0;i<ndofe;i++){
00795 gm[5][i]+=av[i];
00796 }
00797
00798
00799 cmulv(b31r,b32,av);
00800 for (i=0;i<ndofe;i++){
00801 gm[5][i]+=av[i];
00802 }
00803
00804
00805
00806 }
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819 void linhexrot::gnl_grmatrix (matrix &grm,vector &x,vector &y,vector &z,double xi,double eta,double zeta,double &jac)
00820 {
00821 long i;
00822 vector b11(ndofe),b12(ndofe),b13(ndofe),b21(ndofe),b22(ndofe),b23(ndofe),b31(ndofe),b32(ndofe),b33(ndofe);
00823
00824 bvectors (x,y,z,xi,eta,zeta,jac,b11,b12,b13,b21,b22,b23,b31,b32,b33);
00825
00826 for (i=0;i<ndofe;i++){
00827 grm[0][i]=b11[i];
00828 grm[1][i]=b12[i];
00829 grm[2][i]=b13[i];
00830 grm[3][i]=b21[i];
00831 grm[4][i]=b22[i];
00832 grm[5][i]=b23[i];
00833 grm[6][i]=b31[i];
00834 grm[7][i]=b32[i];
00835 grm[8][i]=b33[i];
00836 }
00837 }
00838
00839
00840 void linhexrot::dmatblock (matrix &dd, matrix &d)
00841 {
00842 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];
00843
00844 }
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856 void linhexrot::transf_matrix (ivector &nodes,matrix &tmat)
00857 {
00858 long i,n,m;
00859
00860 fillm (0.0,tmat);
00861
00862 n=nodes.n;
00863 m=tmat.m;
00864 for (i=0;i<m;i++){
00865 tmat[i][i]=1.0;
00866 }
00867
00868 for (i=0;i<n;i++){
00869 if (Mt->nodes[nodes[i]].transf>0){
00870 tmat[i*3+0][i*3]=Mt->nodes[nodes[i]].e1[0];
00871 tmat[i*3+1][i*3]=Mt->nodes[nodes[i]].e1[1];
00872 tmat[i*3+2][i*3]=Mt->nodes[nodes[i]].e1[2];
00873
00874 tmat[i*3+0][i*3+1]=Mt->nodes[nodes[i]].e2[0];
00875 tmat[i*3+1][i*3+1]=Mt->nodes[nodes[i]].e2[1];
00876 tmat[i*3+2][i*3+1]=Mt->nodes[nodes[i]].e2[2];
00877
00878 tmat[i*3+0][i*3+2]=Mt->nodes[nodes[i]].e3[0];
00879 tmat[i*3+1][i*3+2]=Mt->nodes[nodes[i]].e3[1];
00880 tmat[i*3+2][i*3+2]=Mt->nodes[nodes[i]].e3[2];
00881 }
00882 }
00883 }
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896 void linhexrot::gl_stiffness_matrix (long eid,long ri,long ci,matrix &sm)
00897 {
00898 long i,j,k,ipp, ipshear;
00899 double xi,eta,zeta,jac;
00900 vector x(nne),y(nne),z(nne),w,gp;
00901 matrix gm,d(tncomp,tncomp), dd(3,3);
00902
00903 Mt->give_node_coord3d (x,y,z,eid);
00904 fillm (0.0,sm);
00905 fillm (0.0,dd);
00906
00907 allocv (intordsm[0][0],w);
00908 allocv (intordsm[0][0],gp);
00909 allocm (ncomp[0],ndofe,gm);
00910
00911 ipp=Mt->elements[eid].ipp[ri][ci];
00912 gauss_points (gp.a,w.a,intordsm[0][0]);
00913
00914 for (i=0;i<intordsm[0][0];i++){
00915 xi=gp[i];
00916 for (j=0;j<intordsm[0][0];j++){
00917 eta=gp[j];
00918 for (k=0;k<intordsm[0][0];k++){
00919 zeta=gp[k];
00920
00921 Mm->matstiff (d,ipp);
00922 ipp++;
00923
00924 geom_matrix (gm,x,y,z,xi,eta,zeta,jac);
00925
00926 jac*=w[i]*w[j]*w[k];
00927 if (jac<0.0)
00928 fprintf (stderr,"\n negative jacobian on element %ld",eid);
00929 bdbjac (sm,gm,d,gm,jac);
00930
00931
00932 dmatblock ( dd, d);
00933 }
00934 }
00935 }
00936
00937
00938 ipshear=2;
00939 allocv (ipshear,w);
00940 allocv (ipshear,gp);
00941 allocm (3,ndofe,gm);
00942 gauss_points (gp.a,w.a,ipshear);
00943 for (i=0;i<ipshear;i++){
00944 xi=gp[i];
00945 for (j=0;j<ipshear;j++){
00946 eta=gp[j];
00947 for (k=0;k<ipshear;k++){
00948 zeta=gp[k];
00949
00950 geom_matrix_shear (gm,x,y,z,xi,eta,zeta,jac);
00951
00952 jac*=w[i]*w[j]*w[k];
00953 bdbjac (sm,gm,dd,gm,jac);
00954
00955 }
00956 }
00957 }
00958 destrm (gm);
00959 destrv (gp); destrv (w);
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971 }
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985 void linhexrot::gnl_stiffness_matrix (long lcid,long eid,long ri,long ci,matrix &sm)
00986 {
00987 long i,j,k,ipp;
00988 double xi,eta,zeta,jac,jac2;
00989 vector w,gp,sig(tncomp),r(ndofe),x(nne),y(nne),z(nne);
00990 matrix gm(tncomp,ndofe),grm(9,ndofe),d(tncomp,tncomp),s(9,9);
00991
00992
00993 Mt->give_node_coord3d (x,y,z,eid);
00994
00995 eldispl (lcid,eid,r.a);
00996
00997
00998 fillm (0.0,sm);
00999
01000
01001 allocv (intordsm[0][0],w);
01002
01003 allocv (intordsm[0][0],gp);
01004
01005
01006 gauss_points (gp.a,w.a,intordsm[0][0]);
01007
01008
01009 ipp=Mt->elements[eid].ipp[ri][ci];
01010
01011
01012 for (i=0;i<intordsm[0][0];i++){
01013 xi=gp[i];
01014 for (j=0;j<intordsm[0][0];j++){
01015 eta=gp[j];
01016 for (k=0;k<intordsm[0][0];k++){
01017 zeta=gp[k];
01018
01019
01020
01021
01022
01023
01024 gngeom_matrix (gm,r,x,y,z,xi,eta,zeta,jac);
01025
01026
01027 Mm->matstiff (d,ipp);
01028
01029 jac*=w[i]*w[j]*w[k];
01030
01031
01032 bdbjac (sm,gm,d,gm,jac);
01033
01034
01035
01036
01037
01038
01039
01040 gnl_grmatrix (grm,x,y,z,xi,eta,zeta,jac2);
01041
01042
01043 Mm->givestress (lcid,ipp,sig);
01044
01045 s[0][0]=sig[0]; s[0][1]=sig[5]; s[0][2]=sig[4];
01046 s[1][0]=sig[5]; s[1][1]=sig[1]; s[1][2]=sig[3];
01047 s[2][0]=sig[4]; s[2][1]=sig[3]; s[2][2]=sig[2];
01048
01049 s[3][3]=sig[0]; s[3][4]=sig[5]; s[3][5]=sig[4];
01050 s[4][3]=sig[5]; s[4][4]=sig[1]; s[4][5]=sig[3];
01051 s[5][3]=sig[4]; s[5][4]=sig[3]; s[5][5]=sig[2];
01052
01053 s[6][6]=sig[0]; s[6][7]=sig[5]; s[6][8]=sig[4];
01054 s[7][6]=sig[5]; s[7][7]=sig[1]; s[7][8]=sig[3];
01055 s[8][6]=sig[4]; s[8][7]=sig[3]; s[8][8]=sig[2];
01056
01057
01058
01059
01060 bdbjac (sm,grm,s,grm,jac);
01061
01062 ipp++;
01063 }
01064 }
01065 }
01066
01067 destrv (gp); destrv (w);
01068
01069 }
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081 void linhexrot::res_stiffness_matrix (long lcid,long eid,matrix &sm)
01082 {
01083 gl_stiffness_matrix (eid,0,0,sm);
01084
01085 ivector nodes (nne);
01086 Mt->give_elemnodes (eid,nodes);
01087 long transf = Mt->locsystems (nodes);
01088 if (transf>0){
01089 matrix tmat (ndofe,ndofe);
01090 transf_matrix (nodes,tmat);
01091 glmatrixtransf (sm,tmat);
01092 }
01093
01094 }
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104 void linhexrot::mass_matrix (long eid,matrix &mm)
01105 {
01106 long i,j,k;
01107 double jac,xi,eta,zeta,rho;
01108 ivector nodes (nne);
01109 vector x(nne),y(nne),z(nne),w(intordmm),gp(intordmm),dens(nne);
01110 matrix n(napfun,ndofe);
01111
01112 Mt->give_elemnodes (eid,nodes);
01113 Mc->give_density (eid,nodes,dens);
01114
01115 Mt->give_node_coord3d (x,y,z,eid);
01116 gauss_points (gp.a,w.a,intordmm);
01117
01118 fillm (0.0,mm);
01119
01120 for (i=0;i<intordmm;i++){
01121 xi=gp[i];
01122 for (j=0;j<intordmm;j++){
01123 eta=gp[j];
01124 for (k=0;k<intordmm;k++){
01125 zeta=gp[k];
01126
01127 jac_3d (jac,x,y,z,xi,eta,zeta);
01128
01129 jac=fabs(jac);
01130
01131 bf_matrix (n,xi,eta,zeta);
01132
01133 rho = approx (xi,eta,zeta,dens);
01134
01135 jac*=w[i]*w[j]*w[k]*rho;
01136
01137 nnj (mm.a,n.a,jac,n.m,n.n);
01138 }
01139 }
01140 }
01141
01142 }
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152 void linhexrot::res_mass_matrix (long eid,matrix &mm)
01153 {
01154 long transf;
01155
01156 mass_matrix (eid,mm);
01157
01158
01159
01160 ivector nodes (nne);
01161 Mt->give_elemnodes (eid,nodes);
01162 transf = Mt->locsystems (nodes);
01163 if (transf>0){
01164 matrix tmat (ndofe,ndofe);
01165 transf_matrix (nodes,tmat);
01166 glmatrixtransf (mm,tmat);
01167 }
01168 }
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178 void linhexrot::load_matrix (long eid,matrix &lm)
01179 {
01180 long i,j,k;
01181 double jac,xi,eta,zeta,w1,w2,w3;
01182 ivector nodes (nne);
01183 vector x(nne),y(nne),z(nne),w(intordmm),gp(intordmm);
01184 matrix n(napfun,ndofe);
01185
01186 Mt->give_elemnodes (eid,nodes);
01187 Mt->give_node_coord3d (x,y,z,eid);
01188 gauss_points (gp.a,w.a,intordmm);
01189
01190 fillm (0.0,lm);
01191
01192 for (i=0;i<intordmm;i++){
01193 xi=gp[i]; w1=w[i];
01194 for (j=0;j<intordmm;j++){
01195 eta=gp[j]; w2=w[j];
01196 for (k=0;k<intordmm;k++){
01197 zeta=gp[k]; w3=w[k];
01198
01199 jac_3d (jac,x,y,z,xi,eta,zeta);
01200
01201 bf_matrix (n,xi,eta,zeta);
01202
01203 jac*=w1*w2*w3;
01204
01205 nnj (lm.a,n.a,jac,n.m,n.n);
01206 }
01207 }
01208 }
01209
01210 }
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220 void linhexrot::res_load_matrix (long eid,matrix &lm)
01221 {
01222 long transf;
01223 ivector nodes(nne);
01224
01225 load_matrix (eid,lm);
01226
01227
01228
01229 Mt->give_elemnodes (eid,nodes);
01230 transf = Mt->locsystems (nodes);
01231 if (transf>0){
01232 matrix tmat (ndofe,ndofe);
01233 transf_matrix (nodes,tmat);
01234 glmatrixtransf (lm,tmat);
01235 }
01236 }
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249 void linhexrot::res_ip_strains (long lcid,long eid)
01250 {
01251 ivector nodes(nne);
01252 vector x(nne),y(nne),z(nne),r(ndofe),aux;
01253 matrix tmat;
01254
01255 Mt->give_elemnodes (eid,nodes);
01256 Mt->give_node_coord3d (x,y,z,eid);
01257 eldispl (lcid,eid,r.a);
01258
01259
01260
01261 long transf = Mt->locsystems (nodes);
01262 if (transf>0){
01263 allocv (ndofe,aux);
01264 allocm (ndofe,ndofe,tmat);
01265 transf_matrix (nodes,tmat);
01266
01267 lgvectortransf (aux,r,tmat);
01268 copyv (aux,r);
01269 destrv (aux);
01270 destrm (tmat);
01271 }
01272
01273 gl_ip_strains (lcid,eid,0,0,x,y,z,r);
01274
01275
01276 }
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291 void linhexrot::gl_ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &z,vector &r)
01292 {
01293 long i,j,k,ii,ipp;
01294 double xi,eta,zeta,jac;
01295 vector gp,w,eps;
01296 matrix gm;
01297
01298 for (ii=0;ii<nb;ii++){
01299
01300 allocv (intordsm[ii][ii],gp);
01301 allocv (intordsm[ii][ii],w);
01302 allocv (ncomp[ii],eps);
01303 allocm (ncomp[ii],ndofe,gm);
01304
01305 gauss_points (gp.a,w.a,intordsm[ii][ii]);
01306
01307 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01308 for (i=0;i<intordsm[ii][ii];i++){
01309 xi=gp[i];
01310 for (j=0;j<intordsm[ii][ii];j++){
01311 eta=gp[j];
01312 for (k=0;k<intordsm[ii][ii];k++){
01313 zeta=gp[k];
01314
01315 geom_matrix (gm,x,y,z,xi,eta,zeta,jac);
01316
01317 mxv (gm,r,eps);
01318
01319 Mm->storestrain (lcid,ipp,cncomp[ii],eps);
01320
01321 ipp++;
01322 }
01323 }
01324 }
01325
01326 destrm (gm); destrv (w); destrv (gp);
01327 destrv (eps);
01328 }
01329 }
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344 void linhexrot::gnl_ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &z,vector &r)
01345 {
01346 long i,j,k,ipp;
01347 double xi,eta,zeta,jac,b11r,b12r,b13r,b21r,b22r,b23r,b31r,b32r,b33r;
01348 vector gp,w,eps;
01349 vector b11(ndofe),b12(ndofe),b13(ndofe),b21(ndofe),b22(ndofe),b23(ndofe),b31(ndofe),b32(ndofe),b33(ndofe);
01350
01351 allocv (intordsm[0][0],gp);
01352 allocv (intordsm[0][0],w);
01353 allocv (tncomp,eps);
01354
01355 gauss_points (gp.a,w.a,intordsm[0][0]);
01356
01357 ipp=Mt->elements[eid].ipp[ri][ci];
01358 for (i=0;i<intordsm[0][0];i++){
01359 xi=gp[i];
01360 for (j=0;j<intordsm[0][0];j++){
01361 eta=gp[j];
01362 for (k=0;k<intordsm[0][0];k++){
01363 zeta=gp[k];
01364
01365 bvectors (x,y,z,xi,eta,zeta,jac,b11,b12,b13,b21,b22,b23,b31,b32,b33);
01366
01367 scprd (b11,r,b11r);
01368 scprd (b12,r,b12r);
01369 scprd (b13,r,b13r);
01370 scprd (b21,r,b21r);
01371 scprd (b22,r,b22r);
01372 scprd (b23,r,b23r);
01373 scprd (b31,r,b31r);
01374 scprd (b32,r,b32r);
01375 scprd (b33,r,b33r);
01376
01377 eps[0] = b11r + 0.5*b11r*b11r + 0.5*b21r*b21r + 0.5*b31r*b31r;
01378 eps[1] = b22r + 0.5*b12r*b12r + 0.5*b22r*b22r + 0.5*b32r*b32r;
01379 eps[2] = b33r + 0.5*b13r*b13r + 0.5*b23r*b23r + 0.5*b33r*b33r;
01380
01381 eps[3] = b23r+b32r + b12r*b13r + b22r*b23r + b32r*b33r;
01382 eps[4] = b31r+b13r + b13r*b11r + b23r*b21r + b33r*b31r;
01383 eps[5] = b12r+b21r + b11r*b12r + b21r*b22r + b31r*b32r;
01384
01385
01386 Mm->storestrain (lcid,ipp,eps);
01387 ipp++;
01388 }
01389 }
01390 }
01391
01392 destrv (eps); destrv (w); destrv (gp);
01393
01394 }
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406 void linhexrot::nod_strains_ip (long lcid,long eid,long ri,long ci)
01407 {
01408 long i,j,ipp;
01409 ivector ipnum(nne),nod(nne);
01410 vector eps(tncomp);
01411
01412
01413
01414 ipp=Mt->elements[eid].ipp[ri][ci];
01415
01416
01417
01418 Mt->give_elemnodes (eid,nod);
01419
01420 for (i=0;i<nne;i++){
01421
01422 Mm->givestrain (lcid,ipnum[i],eps);
01423
01424
01425 j=nod[i];
01426 Mt->nodes[j].storestrain (lcid,0,eps);
01427 }
01428 }
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441 void linhexrot::nod_strains_comp (long lcid,long eid,double **stra)
01442 {
01443 long i,j;
01444 double jac;
01445 ivector nodes(nne);
01446 vector x(nne),y(nne),z(nne),nxi(nne),neta(nne),nzeta(nne),r(ndofe),eps(tncomp),aux;
01447 matrix tmat,gm(tncomp,ndofe);
01448
01449
01450 Mt->give_node_coord3d (x,y,z,eid);
01451
01452 Mt->give_elemnodes (eid,nodes);
01453
01454 eldispl (lcid,eid,r.a);
01455
01456
01457 long transf = Mt->locsystems (nodes);
01458 if (transf>0){
01459 allocv (ndofe,aux);
01460 allocm (ndofe,ndofe,tmat);
01461 transf_matrix (nodes,tmat);
01462
01463 lgvectortransf (aux,r,tmat);
01464 copyv (aux,r);
01465 destrv (aux);
01466 destrm (tmat);
01467 }
01468
01469
01470
01471
01472
01473
01474 for (i=0;i<nne;i++){
01475
01476 geom_matrix (gm,x,y,z,nxi[i],neta[i],nzeta[i],jac);
01477
01478 mxv (gm,r,eps);
01479
01480 for (j=0;j<eps.n;j++){
01481 stra[i][j]=eps[j];
01482 }
01483 }
01484
01485 }
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501 void linhexrot::appstrain (long lcid,long eid,double xi,double eta,double zeta,long fi,long ncomp,vector &eps)
01502 {
01503 long i,j,k;
01504 ivector nod(nne);
01505 vector nodval(nne);
01506
01507 if (ncomp != eps.n){
01508 fprintf (stderr,"\n\n wrong interval of indices in function strain (%s, line %d).\n",__FILE__,__LINE__);
01509 abort ();
01510 }
01511
01512 Mt->give_elemnodes (eid,nod);
01513 k=0;
01514 for (i=fi;i<fi+ncomp;i++){
01515 for (j=0;j<nne;j++){
01516 nodval[j]=Mt->nodes[nod[j]].strain[lcid*tncomp+i];
01517 }
01518 eps[k]=approx (xi,eta,zeta,nodval);
01519 k++;
01520 }
01521 }
01522
01523
01524
01525 void linhexrot::strains (long lcid,long eid,long ri,long ci)
01526 {
01527 long i,naep,ncp,sid;
01528 double **stra;
01529 vector coord,eps;
01530
01531 if (Mp->strainaver==0){
01532 stra = new double* [nne];
01533 for (i=0;i<nne;i++){
01534 stra[i] = new double [tncomp];
01535 }
01536
01537 }
01538
01539 switch (Mm->stra.tape[eid]){
01540 case nowhere:{
01541 break;
01542 }
01543 case intpts:{
01544
01545 break;
01546 }
01547 case enodes:{
01548 break;
01549 }
01550 case userdefined:{
01551
01552 naep = Mm->stra.give_naep (eid);
01553 ncp = Mm->stra.give_ncomp (eid);
01554 sid = Mm->stra.give_sid (eid);
01555 allocv (ncp,eps);
01556 allocv (3,coord);
01557 for (i=0;i<naep;i++){
01558 Mm->stra.give_aepcoord (sid,i,coord);
01559
01560 if (Mp->strainaver==0)
01561
01562 if (Mp->strainaver==1)
01563 appstrain (lcid,eid,coord[0],coord[1],coord[2],0,ncp,eps);
01564
01565 Mm->stra.storevalues(lcid,eid,i,eps);
01566 }
01567 destrv (eps);
01568 destrv (coord);
01569 break;
01570 }
01571 default:{
01572 fprintf (stderr,"\n\n unknown strain point is required in function planeelemlq::strains (%s, line %d).\n",__FILE__,__LINE__);
01573 }
01574 }
01575
01576 if (Mp->strainaver==0){
01577 for (i=0;i<nne;i++){
01578 delete [] stra[i];
01579 }
01580 delete [] stra;
01581 }
01582 }
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592 void linhexrot::res_ip_stresses (long lcid,long eid)
01593 {
01594 ip_stresses (lcid,eid,0,0);
01595 }
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608 void linhexrot::ip_stresses (long lcid,long eid,long ri,long ci)
01609 {
01610 long i,j,k,ipp;
01611
01612 ipp=Mt->elements[eid].ipp[ri][ci];
01613
01614 for (i=0;i<intordsm[0][0];i++){
01615 for (j=0;j<intordsm[0][0];j++){
01616 for (k=0;k<intordsm[0][0];k++){
01617
01618
01619 if (Mp->strcomp==1)
01620 Mm->computenlstresses (ipp);
01621
01622 ipp++;
01623 }
01624 }
01625 }
01626 }
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639 void linhexrot::ip_elast_stresses (long lcid,long eid,long ri,long ci)
01640 {
01641 long i,j,k,ipp;
01642 vector gp,w,eps(tncomp),sig(tncomp);
01643 matrix d(tncomp,tncomp);
01644
01645 allocv (intordsm[0][0],gp);
01646 allocv (intordsm[0][0],w);
01647
01648 gauss_points (gp.a,w.a,intordsm[0][0]);
01649 ipp=Mt->elements[eid].ipp[ri][ci];
01650
01651 for (i=0;i<intordsm[0][0];i++){
01652 for (j=0;j<intordsm[0][0];j++){
01653 for (k=0;k<intordsm[0][0];k++){
01654
01655 Mm->matstiff (d,ipp);
01656
01657 Mm->givestrain (lcid,ipp,eps);
01658
01659 mxv (d,eps,sig);
01660
01661 Mm->storestress (lcid,ipp,sig);
01662
01663 ipp++;
01664 }
01665 }
01666 }
01667
01668 destrv (w); destrv (gp);
01669 }
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680 void linhexrot::nod_stresses_ip (long lcid,long eid,long ri,long ci)
01681 {
01682 long i,j,ipp;
01683 ivector ipnum(nne),nod(nne);
01684 vector sig(tncomp);
01685
01686
01687
01688 ipp=Mt->elements[eid].ipp[ri][ci];
01689
01690
01691
01692 Mt->give_elemnodes (eid,nod);
01693
01694 for (i=0;i<nne;i++){
01695
01696 Mm->givestress (lcid,ipnum[i],sig);
01697
01698
01699 j=nod[i];
01700 Mt->nodes[j].storestress (lcid,0,sig);
01701 }
01702 }
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716 void linhexrot::appstress (long lcid,long eid,double xi,double eta,double zeta,long fi,long ncomp,vector &sig)
01717 {
01718 long i,j,k;
01719 ivector nodes(nne);
01720 vector nodval(nne);
01721
01722 if (ncomp != sig.n){
01723 fprintf (stderr,"\n\n wrong interval of indices in function stress (%s, line %d).\n",__FILE__,__LINE__);
01724 abort ();
01725 }
01726
01727 Mt->give_elemnodes (eid,nodes);
01728 k=0;
01729 for (i=fi;i<fi+ncomp;i++){
01730 for (j=0;j<nne;j++){
01731 nodval[j]=Mt->nodes[nodes[j]].stress[lcid*tncomp+i];
01732 }
01733 sig[k]=approx (xi,eta,zeta,nodval);
01734 k++;
01735 }
01736 }
01737
01738
01739
01740 void linhexrot::stresses (long lcid,long eid,long ri,long ci)
01741 {
01742 long i,naep,ncp,sid;
01743
01744 vector coord,sig;
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760 switch (Mm->stre.tape[eid]){
01761 case nowhere:{
01762 break;
01763 }
01764 case intpts:{
01765
01766 break;
01767 }
01768 case enodes:{
01769 break;
01770 }
01771 case userdefined:{
01772
01773 naep = Mm->stre.give_naep (eid);
01774 ncp = Mm->stre.give_ncomp (eid);
01775 sid = Mm->stre.give_sid (eid);
01776 allocv (ncp,sig);
01777 allocv (3,coord);
01778 for (i=0;i<naep;i++){
01779 Mm->stre.give_aepcoord (sid,i,coord);
01780
01781 if (Mp->stressaver==0)
01782
01783 if (Mp->stressaver==1)
01784 appstress (lcid,eid,coord[0],coord[1],coord[2],0,ncp,sig);
01785
01786 Mm->stre.storevalues(lcid,eid,i,sig);
01787 }
01788 destrv (sig);
01789 destrv (coord);
01790 break;
01791 }
01792 default:{
01793 fprintf (stderr,"\n\n unknown stress point is required in function planeelemlq::stresses (%s, line %d).\n",__FILE__,__LINE__);
01794 }
01795 }
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810 }
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821 void linhexrot::nod_eqother_ip (long lcid,long eid,long ri,long ci)
01822 {
01823 long i,j,ipp,ncompo;
01824 ivector ipnum(nne),nod(nne);
01825 vector eqother;
01826
01827
01828
01829 ipp=Mt->elements[eid].ipp[ri][ci];
01830
01831
01832
01833 Mt->give_elemnodes (eid,nod);
01834
01835 for (i=0;i<nne;i++){
01836
01837 ncompo = Mm->givencompeqother (ipnum[i],0);
01838 allocv (ncompo,eqother);
01839
01840 Mm->giveeqother (ipnum[i],0,ncompo,eqother.a);
01841
01842
01843 j=nod[i];
01844 Mt->nodes[j].storeother (lcid,0,ncompo,eqother);
01845
01846 destrv (eqother);
01847 }
01848 }
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875 void linhexrot::gl_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x, vector &y, vector &z)
01876 {
01877 integratedquant iq;
01878 iq=locstress;
01879
01880
01881 compute_nlstress (lcid,eid,ri,ci);
01882
01883
01884 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
01885 }
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898 void linhexrot::gnl_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
01899 {
01900 long i,j,k,ipp;
01901 double xi,eta,zeta,jac;
01902 vector w,gp,x(nne),y(nne),z(nne),sig(tncomp),contr(ndofe),r(ndofe);
01903 matrix gm(tncomp,ndofe);
01904
01905
01906 Mt->give_node_coord3d (x,y,z,eid);
01907
01908 eldispl (lcid,eid,r.a);
01909
01910
01911 fillv (0.0,ifor);
01912
01913
01914 allocv (intordsm[0][0],gp);
01915
01916 allocv (intordsm[0][0],w);
01917
01918
01919 gauss_points (gp.a,w.a,intordsm[0][0]);
01920
01921
01922 ipp=Mt->elements[eid].ipp[ri][ci];
01923
01924
01925 for (i=0;i<intordsm[0][0];i++){
01926 xi=gp[i];
01927 for (j=0;j<intordsm[0][0];j++){
01928 eta=gp[j];
01929 for (k=0;k<intordsm[0][0];k++){
01930 zeta=gp[k];
01931
01932
01933 if (Mp->strcomp==1)
01934 Mm->computenlstresses (ipp);
01935
01936 if (Mp->nodeintfor==1){
01937 Mm->givestress (lcid,ipp,sig);
01938
01939
01940 gngeom_matrix (gm,r,x,y,z,xi,eta,zeta,jac);
01941
01942 mtxv (gm,sig,contr);
01943
01944 cmulv (jac*w[i]*w[j]*w[k],contr);
01945
01946 addv (contr,ifor,ifor);
01947 }
01948
01949 ipp++;
01950 }
01951 }
01952 }
01953 destrv (w); destrv (gp);
01954 }
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969 void linhexrot::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y,vector &z)
01970 {
01971 integratedquant iq;
01972 iq=nonlocstress;
01973
01974
01975 compute_nonloc_nlstress (lcid,eid,ri,ci);
01976
01977
01978 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
01979
01980 }
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995 void linhexrot::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y,vector &z)
01996 {
01997 integratedquant iq;
01998 iq=stressincr;
01999
02000
02001 compute_nlstressincr (lcid,eid,ri,ci);
02002
02003 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
02004 }
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019 void linhexrot::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y,vector &z)
02020 {
02021 integratedquant iq;
02022 iq=eigstress;
02023
02024
02025 if (Mp->eigstrcomp)
02026 compute_eigstress (lcid,eid,ri,ci);
02027
02028
02029 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y,z);
02030 }
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043 void linhexrot::res_internal_forces (long lcid,long eid,vector &ifor)
02044 {
02045 long transf;
02046 ivector nodes (nne);
02047 vector v(ndofe),x(nne),y(nne),z(nne);
02048
02049 Mt->give_node_coord3d (x,y,z,eid);
02050
02051 gl_internal_forces (lcid,eid,0,0,ifor,x,y,z);
02052
02053
02054
02055
02056
02057 Mt->give_elemnodes (eid,nodes);
02058 transf = Mt->locsystems (nodes);
02059 if (transf>0){
02060 matrix tmat (ndofe,ndofe);
02061 transf_matrix (nodes,tmat);
02062
02063 glvectortransf (ifor,v,tmat);
02064 copyv (v,ifor);
02065 }
02066 }
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079 void linhexrot::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
02080 {
02081 long transf;
02082 ivector nodes (nne);
02083 vector v(ndofe),x(nne),y(nne),z(nne);
02084
02085 Mt->give_node_coord3d (x,y,z,eid);
02086
02087 nonloc_internal_forces (lcid,eid,0,0,ifor,x,y,z);
02088
02089
02090
02091 Mt->give_elemnodes (eid,nodes);
02092 transf = Mt->locsystems (nodes);
02093 if (transf>0){
02094 matrix tmat (ndofe,ndofe);
02095 transf_matrix (nodes,tmat);
02096
02097 glvectortransf (ifor,v,tmat);
02098 copyv (v,ifor);
02099 }
02100 }
02101
02102
02103
02104 void linhexrot::res_incr_internal_forces (long lcid,long eid,vector &ifor)
02105 {
02106 long transf;
02107 ivector nodes (nne);
02108 vector v(ndofe),x(nne),y(nne),z(nne);
02109
02110 Mt->give_node_coord3d (x,y,z,eid);
02111
02112 incr_internal_forces (lcid,eid,0,0,ifor,x,y,z);
02113
02114
02115
02116 Mt->give_elemnodes (eid,nodes);
02117 transf = Mt->locsystems (nodes);
02118 if (transf>0){
02119 matrix tmat (ndofe,ndofe);
02120 transf_matrix (nodes,tmat);
02121
02122 glvectortransf (ifor,v,tmat);
02123 copyv (v,ifor);
02124 }
02125 }
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137 void linhexrot::res_eigstrain_forces (long lcid,long eid,vector &nfor)
02138 {
02139 long transf;
02140 ivector nodes (nne);
02141 vector v(ndofe),x(nne),y(nne),z(nne);
02142
02143 Mt->give_node_coord3d (x,y,z,eid);
02144
02145 eigstrain_forces (lcid,eid,0,0,nfor,x,y,z);
02146
02147
02148
02149 Mt->give_elemnodes (eid,nodes);
02150 transf = Mt->locsystems (nodes);
02151 if (transf>0){
02152 matrix tmat (ndofe,ndofe);
02153 transf_matrix (nodes,tmat);
02154
02155 glvectortransf (nfor,v,tmat);
02156 copyv (v,nfor);
02157 }
02158 }
02159
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169
02170
02171 void linhexrot::compute_nlstress (long lcid,long eid,long ri,long ci)
02172 {
02173 long i,j,k,ipp;
02174
02175 ipp=Mt->elements[eid].ipp[ri][ci];
02176
02177 for (i=0;i<intordsm[0][0];i++){
02178 for (j=0;j<intordsm[0][0];j++){
02179 for (k=0;k<intordsm[0][0];k++){
02180
02181
02182 if (Mp->strcomp==1)
02183 Mm->computenlstresses (ipp);
02184
02185 ipp++;
02186 }
02187 }
02188 }
02189 }
02190
02191
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201
02202 void linhexrot::compute_nlstressincr (long lcid,long eid,long ri,long ci)
02203 {
02204 long i,j,k,ipp;
02205
02206 ipp=Mt->elements[eid].ipp[ri][ci];
02207
02208 for (i=0;i<intordsm[0][0];i++){
02209 for (j=0;j<intordsm[0][0];j++){
02210 for (k=0;k<intordsm[0][0];k++){
02211
02212
02213 if (Mp->strcomp==1)
02214 Mm->computenlstressesincr (ipp);
02215
02216 ipp++;
02217 }
02218 }
02219 }
02220 }
02221
02222
02223
02224
02225
02226
02227
02228
02229
02230
02231
02232
02233
02234 void linhexrot::local_values (long lcid,long eid,long ri,long ci)
02235 {
02236 long i,j,k,ipp;
02237
02238 ipp=Mt->elements[eid].ipp[ri][ci];
02239
02240 for (i=0;i<intordsm[0][0];i++){
02241 for (j=0;j<intordsm[0][0];j++){
02242 for (k=0;k<intordsm[0][0];k++){
02243
02244
02245 if (Mp->strcomp==1)
02246 Mm->computenlstresses (ipp);
02247
02248 ipp++;
02249 }
02250 }
02251 }
02252 }
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264
02265 void linhexrot::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
02266 {
02267 long i,j,k,ipp;
02268
02269 ipp=Mt->elements[eid].ipp[ri][ci];
02270 for (i=0;i<intordsm[0][0];i++){
02271 for (j=0;j<intordsm[0][0];j++){
02272 for (k=0;k<intordsm[0][0];k++){
02273
02274 if (Mp->strcomp==1)
02275 Mm->compnonloc_nlstresses (ipp);
02276 ipp++;
02277 }
02278 }
02279 }
02280 }
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293 void linhexrot::compute_eigstress (long lcid,long eid,long ri,long ci)
02294 {
02295 long i,j,k,ipp;
02296 vector eigstr(tncomp),sig(tncomp);
02297 matrix d(tncomp,tncomp);
02298
02299 ipp=Mt->elements[eid].ipp[ri][ci];
02300
02301 for (i=0;i<intordsm[0][0];i++){
02302 for (j=0;j<intordsm[0][0];j++){
02303 for (k=0;k<intordsm[0][0];k++){
02304 Mm->giveeigstrain (ipp,eigstr);
02305
02306 Mm->matstiff (d,ipp);
02307 mxv (d,eigstr,sig);
02308 Mm->storeeigstress (ipp,sig);
02309 ipp++;
02310 }
02311 }
02312 }
02313 }
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331 void linhexrot::elem_integration (integratedquant iq,long lcid,long eid,long ri,long ci,vector &nv,vector &x, vector &y, vector &z)
02332 {
02333 long i,j,k,ipp;
02334 double xi,eta,zeta,jac;
02335 vector w,gp,ipv(tncomp),contr(ndofe);
02336 matrix gm(tncomp,ndofe);
02337
02338 fillv (0.0,nv);
02339
02340 allocv (intordsm[0][0],gp);
02341 allocv (intordsm[0][0],w);
02342
02343 gauss_points (gp.a,w.a,intordsm[0][0]);
02344 ipp=Mt->elements[eid].ipp[ri][ci];
02345
02346 for (i=0;i<intordsm[0][0];i++){
02347 xi=gp[i];
02348 for (j=0;j<intordsm[0][0];j++){
02349 eta=gp[j];
02350 for (k=0;k<intordsm[0][0];k++){
02351 zeta=gp[k];
02352
02353 Mm->givequantity (iq,lcid,ipp,cncomp[0],ipv);
02354
02355 geom_matrix (gm,x,y,z,xi,eta,zeta,jac);
02356
02357 mtxv (gm,ipv,contr);
02358 cmulv (jac*w[i]*w[j]*w[k],contr);
02359
02360 addv (contr,nv,nv);
02361 ipp++;
02362 }
02363 }
02364 }
02365 destrv (w); destrv (gp);
02366 }
02367
02368
02369
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385 void linhexrot::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
02386 {
02387 long i,j,k,ii;
02388 double xi,eta,zeta;
02389 vector x(nne),y(nne),z(nne),w(intordsm[ri][ci]),gp(intordsm[ri][ci]);
02390
02391 gauss_points (gp.a,w.a,intordsm[ri][ci]);
02392 Mt->give_node_coord3d (x,y,z,eid);
02393 ii=Mt->elements[eid].ipp[ri][ci];
02394
02395 for (i=0;i<intordsm[ri][ci];i++){
02396 xi=gp[i];
02397 for (j=0;j<intordsm[ri][ci];j++){
02398 eta=gp[j];
02399 for (k=0;k<intordsm[ri][ci];k++){
02400 zeta=gp[k];
02401 if (ii==ipp){
02402 coord[0]=approx (xi,eta,zeta,x);
02403 coord[1]=approx (xi,eta,zeta,y);
02404 coord[2]=approx (xi,eta,zeta,z);
02405 }
02406 ii++;
02407 }
02408 }
02409 }
02410 }
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433
02434
02435
02436 void linhexrot::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
02437 {
02438 long i, j, k, l, m, ipp;
02439 long ii, jj, nv = nodval.n;
02440 double xi, eta, zeta, ipval;
02441 vector w, gp, anv(nne);
02442 long nstra, nstre, ncompstr, ncompeqother;
02443 long idstra, idstre, idoth, idic;
02444 inictype ict;
02445
02446 nstra = idstra = nstre = idstre = idoth = idic = 0;
02447
02448 ict = ictn[0];
02449 for (i=0; i<nne; i++)
02450 {
02451 if (ictn[i] != ict)
02452 {
02453 print_err("Incompatible types of initial conditions on element %ld\n"
02454 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
02455 abort();
02456 }
02457 }
02458 for (j = 0; j < nv; j++)
02459 {
02460 for(i = 0; i < nne; i++)
02461 anv[i] = nodval[i][j];
02462 for (ii = 0; ii < nb; ii++)
02463 {
02464 for (jj = 0; jj < nb; jj++)
02465 {
02466 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
02467 if (intordsm[ii][jj] == 0)
02468 continue;
02469 allocv (intordsm[ii][jj],gp);
02470 allocv (intordsm[ii][jj],w);
02471 gauss_points (gp.a,w.a,intordsm[ii][jj]);
02472 for (k = 0; k < intordsm[ii][jj]; k++)
02473 {
02474 xi=gp[k];
02475 for (l = 0; l < intordsm[ii][jj]; l++)
02476 {
02477 eta=gp[l];
02478 for (m = 0; m < intordsm[ii][jj]; m++)
02479 {
02480 zeta=gp[m];
02481
02482 ipval = approx (xi,eta,zeta,anv);
02483 ncompstr = Mm->ip[ipp].ncompstr;
02484 ncompeqother = Mm->ip[ipp].ncompeqother;
02485 if ((ictn[0] & inistrain) && (j < ncompstr))
02486 {
02487 Mm->ip[ipp].strain[idstra] += ipval;
02488 ipp++;
02489 continue;
02490 }
02491 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
02492 {
02493 Mm->ip[ipp].stress[idstre] += ipval;
02494 ipp++;
02495 continue;
02496 }
02497 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
02498 {
02499 Mm->ip[ipp].eqother[idoth] += ipval;
02500 ipp++;
02501 continue;
02502 }
02503 if ((ictn[0] & inicond) && (j < nv))
02504 {
02505 if (Mm->ic[ipp] == NULL)
02506 {
02507 Mm->ic[ipp] = new double[nv-j];
02508 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
02509 }
02510 Mm->ic[ipp][idic] += ipval;
02511 ipp++;
02512 continue;
02513 }
02514 ipp++;
02515 }
02516 }
02517 }
02518 destrv (gp); destrv (w);
02519 }
02520 }
02521 ipp=Mt->elements[eid].ipp[ri][ci];
02522 ncompstr = Mm->ip[ipp].ncompstr;
02523 ncompeqother = Mm->ip[ipp].ncompeqother;
02524 if ((ictn[0] & inistrain) && (j < ncompstr))
02525 {
02526 nstra++;
02527 idstra++;
02528 continue;
02529 }
02530 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
02531 {
02532 nstre++;
02533 idstre++;
02534 continue;
02535 }
02536 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
02537 {
02538 idoth++;
02539 continue;
02540 }
02541 if ((ictn[0] & inicond) && (j < nv))
02542 {
02543 idic++;
02544 continue;
02545 }
02546 }
02547 }
02548
02549
02550
02551
02552
02553
02554
02555
02556 void linhexrot::ipvolume (long eid,long ri,long ci)
02557 {
02558 long i,j,k,ii,jj,ipp;
02559 double xi,eta,zeta,jac;
02560 vector x(nne),y(nne),z(nne),w,gp;
02561
02562 Mt->give_node_coord3d (x,y,z,eid);
02563
02564 for (ii=0;ii<nb;ii++){
02565 for (jj=0;jj<nb;jj++){
02566 if (intordsm[ii][jj]==0) continue;
02567
02568 allocv (intordsm[ii][jj],w);
02569 allocv (intordsm[ii][jj],gp);
02570
02571 gauss_points (gp.a,w.a,intordsm[ii][jj]);
02572
02573 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
02574
02575 for (i=0;i<intordsm[ii][jj];i++){
02576 xi=gp[i];
02577 for (j=0;j<intordsm[ii][jj];j++){
02578 eta=gp[j];
02579 for (k=0;k<intordsm[ii][jj];k++){
02580 zeta=gp[k];
02581
02582 jac_3d (jac,x,y,z,xi,eta,zeta);
02583 jac=fabs(jac);
02584
02585 jac*=w[i]*w[j]*w[k];
02586
02587 Mm->storeipvol (ipp,jac);
02588 ipp++;
02589 }
02590 }
02591 }
02592 destrv (gp); destrv (w);
02593 }
02594 }
02595
02596 }
02597
02598
02599
02600
02601
02602
02603
02604
02605
02606
02607
02608 void linhexrot::node_forces_surf (long lcid,long eid,long *is,double *nv,vector &nf)
02609 {
02610 long i,j;
02611 double xi,eta,zeta,jac;
02612 double *tnv;
02613 vector x(nne),y(nne),z(nne),gp,w,av(ndofe),v(ndofe);
02614 matrix n(napfun,ndofe),an(napfun,ndofe),am(ndofe,ndofe),tran(3,3);
02615
02616 tnv = new double [12];
02617
02618
02619 Mt->give_node_coord3d (x,y,z,eid);
02620
02621 allocv (intordb,w);
02622 allocv (intordb,gp);
02623 gauss_points (gp.a,w.a,intordb);
02624
02625
02626
02627 if (is[0]>0 ){
02628 xi=1.0;
02629 for (i=0;i<intordb;i++){
02630 eta=gp[i];
02631 for (j=0;j<intordb;j++){
02632 zeta=gp[j];
02633
02634 jac2d_3d (jac,x,y,z,eta,zeta,0);
02635 bf_matrix (n,xi,eta,zeta);
02636 jac = jac*w[i]*w[j];
02637 nnj (am.a,n.a,jac,n.m,n.n);
02638 }
02639 }
02640
02641 if (is[0]==1){
02642 av[0] = nv[0*3+0];
02643 av[1] = nv[0*3+1];
02644 av[2] = nv[0*3+2];
02645
02646 av[9] = nv[1*3+0];
02647 av[10] = nv[1*3+1];
02648 av[11] = nv[1*3+2];
02649
02650 av[21] = nv[2*3+0];
02651 av[22] = nv[2*3+1];
02652 av[23] = nv[2*3+2];
02653
02654 av[12] = nv[3*3+0];
02655 av[13] = nv[3*3+1];
02656 av[14] = nv[3*3+2];
02657 }
02658
02659
02660 if (is[0]==2){
02661
02662 av[0] = nv[0*3+0];
02663 av[1] = nv[0*3+1];
02664 av[2] = nv[0*3+2];
02665
02666 av[3] = nv[1*3+0];
02667 av[4] = nv[1*3+1];
02668 av[5] = nv[1*3+2];
02669
02670 av[6] = nv[2*3+0];
02671 av[7] = nv[2*3+1];
02672 av[8] = nv[2*3+2];
02673
02674 av[9] = nv[3*3+0];
02675 av[10] = nv[3*3+1];
02676 av[11] = nv[3*3+2];
02677
02678 locglob_nodeval (0,av,tnv,x,y,z);
02679
02680 av[0] = tnv[0*3+0];
02681 av[1] = tnv[0*3+1];
02682 av[2] = tnv[0*3+2];
02683
02684 av[9] = tnv[1*3+0];
02685 av[10] = tnv[1*3+1];
02686 av[11] = tnv[1*3+2];
02687
02688 av[21] = tnv[2*3+0];
02689 av[22] = tnv[2*3+1];
02690 av[23] = tnv[2*3+2];
02691
02692 av[12] = tnv[3*3+0];
02693 av[13] = tnv[3*3+1];
02694 av[14] = tnv[3*3+2];
02695 }
02696
02697
02698 mxv (am,av,v);
02699 addv (v,nf,nf);
02700 }
02701
02702
02703
02704
02705 if (is[1]>0 ){
02706 eta=1.0;
02707 for (i=0;i<intordb;i++){
02708 xi=gp[i];
02709 for (j=0;j<intordb;j++){
02710 zeta=gp[j];
02711
02712 jac2d_3d (jac,x,y,z,xi,zeta,1);
02713 bf_matrix (n,xi,eta,zeta);
02714 jac = jac*w[i]*w[j];
02715 nnj (am.a,n.a,jac,n.m,n.n);
02716 }
02717 }
02718
02719 if (is[1]==1){
02720 av[3] = nv[12+0*3+0];
02721 av[4] = nv[12+0*3+1];
02722 av[5] = nv[12+0*3+2];
02723
02724 av[0] = nv[12+1*3+0];
02725 av[1] = nv[12+1*3+1];
02726 av[2] = nv[12+1*3+2];
02727
02728 av[12] = nv[12+2*3+0];
02729 av[13] = nv[12+2*3+1];
02730 av[14] = nv[12+2*3+2];
02731
02732 av[15] = nv[12+3*3+0];
02733 av[16] = nv[12+3*3+1];
02734 av[17] = nv[12+3*3+2];
02735 }
02736
02737
02738 if (is[1]==2){
02739
02740 av[0] = nv[12+0*3+0];
02741 av[1] = nv[12+0*3+1];
02742 av[2] = nv[12+0*3+2];
02743
02744 av[3] = nv[12+1*3+0];
02745 av[4] = nv[12+1*3+1];
02746 av[5] = nv[12+1*3+2];
02747
02748 av[6] = nv[12+2*3+0];
02749 av[7] = nv[12+2*3+1];
02750 av[8] = nv[12+2*3+2];
02751
02752 av[9] = nv[12+3*3+0];
02753 av[10] = nv[12+3*3+1];
02754 av[11] = nv[12+3*3+2];
02755
02756 locglob_nodeval (1,av,tnv,x,y,z);
02757
02758 av[3] = tnv[0*3+0];
02759 av[4] = tnv[0*3+1];
02760 av[5] = tnv[0*3+2];
02761
02762 av[0] = tnv[1*3+0];
02763 av[1] = tnv[1*3+1];
02764 av[2] = tnv[1*3+2];
02765
02766 av[12] = tnv[2*3+0];
02767 av[13] = tnv[2*3+1];
02768 av[14] = tnv[2*3+2];
02769
02770 av[15] = tnv[3*3+0];
02771 av[16] = tnv[3*3+1];
02772 av[17] = tnv[3*3+2];
02773
02774 }
02775
02776 mxv (am,av,v);
02777 addv (v,nf,nf);
02778 }
02779
02780
02781
02782 if (is[2]>0 ){
02783 xi=-1.0;
02784 for (i=0;i<intordb;i++){
02785 eta=gp[i];
02786 for (j=0;j<intordb;j++){
02787 zeta=gp[j];
02788
02789 jac2d_3d (jac,x,y,z,eta,zeta,2);
02790 bf_matrix (n,xi,eta,zeta);
02791 jac = jac*w[i]*w[j];
02792 nnj (am.a,n.a,jac,n.m,n.n);
02793 }
02794 }
02795
02796 if (is[2]==1){
02797 av[6] = nv[24+0*3+0];
02798 av[7] = nv[24+0*3+1];
02799 av[8] = nv[24+0*3+2];
02800
02801 av[3] = nv[24+1*3+0];
02802 av[4] = nv[24+1*3+1];
02803 av[5] = nv[24+1*3+2];
02804
02805 av[15] = nv[24+2*3+0];
02806 av[16] = nv[24+2*3+1];
02807 av[17] = nv[24+2*3+2];
02808
02809 av[18] = nv[24+3*3+0];
02810 av[19] = nv[24+3*3+1];
02811 av[20] = nv[24+3*3+2];
02812 }
02813
02814
02815 if (is[2]==2){
02816
02817 av[0] = nv[24+0*3+0];
02818 av[1] = nv[24+0*3+1];
02819 av[2] = nv[24+0*3+2];
02820
02821 av[3] = nv[24+1*3+0];
02822 av[4] = nv[24+1*3+1];
02823 av[5] = nv[24+1*3+2];
02824
02825 av[6] = nv[24+2*3+0];
02826 av[7] = nv[24+2*3+1];
02827 av[8] = nv[24+2*3+2];
02828
02829 av[9] = nv[24+3*3+0];
02830 av[10] = nv[24+3*3+1];
02831 av[11] = nv[24+3*3+2];
02832
02833 locglob_nodeval (2,av,tnv,x,y,z);
02834
02835 av[6] = tnv[0*3+0];
02836 av[7] = tnv[0*3+1];
02837 av[8] = tnv[0*3+2];
02838
02839 av[3] = tnv[1*3+0];
02840 av[4] = tnv[1*3+1];
02841 av[5] = tnv[1*3+2];
02842
02843 av[15] = tnv[2*3+0];
02844 av[16] = tnv[2*3+1];
02845 av[17] = tnv[2*3+2];
02846
02847 av[18] = tnv[3*3+0];
02848 av[19] = tnv[3*3+1];
02849 av[20] = tnv[3*3+2];
02850
02851 }
02852
02853 mxv (am,av,v);
02854 addv (v,nf,nf);
02855 }
02856
02857
02858
02859
02860 if (is[3]>0 ){
02861 eta=-1.0;
02862 for (i=0;i<intordb;i++){
02863 xi=gp[i];
02864 for (j=0;j<intordb;j++){
02865 zeta=gp[j];
02866
02867 jac2d_3d (jac,x,y,z,xi,zeta,3);
02868 bf_matrix (n,xi,eta,zeta);
02869 jac = jac*w[i]*w[j];
02870 nnj (am.a,n.a,jac,n.m,n.n);
02871 }
02872 }
02873
02874 if (is[3]==1){
02875 av[9] = nv[36+0*3+0];
02876 av[10] = nv[36+0*3+1];
02877 av[11] = nv[36+0*3+2];
02878
02879 av[6] = nv[36+1*3+0];
02880 av[7] = nv[36+1*3+1];
02881 av[8] = nv[36+1*3+2];
02882
02883 av[18] = nv[36+2*3+0];
02884 av[19] = nv[36+2*3+1];
02885 av[20] = nv[36+2*3+2];
02886
02887 av[21] = nv[36+3*3+0];
02888 av[22] = nv[36+3*3+1];
02889 av[23] = nv[36+3*3+2];
02890 }
02891
02892
02893 if (is[3]==2){
02894
02895 av[0] = nv[36+0*3+0];
02896 av[1] = nv[36+0*3+1];
02897 av[2] = nv[36+0*3+2];
02898
02899 av[3] = nv[36+1*3+0];
02900 av[4] = nv[36+1*3+1];
02901 av[5] = nv[36+1*3+2];
02902
02903 av[6] = nv[36+2*3+0];
02904 av[7] = nv[36+2*3+1];
02905 av[8] = nv[36+2*3+2];
02906
02907 av[9] = nv[36+3*3+0];
02908 av[10] = nv[36+3*3+1];
02909 av[11] = nv[36+3*3+2];
02910
02911 locglob_nodeval (3,av,tnv,x,y,z);
02912
02913 av[9] = tnv[0*3+0];
02914 av[10]= tnv[0*3+1];
02915 av[11]= tnv[0*3+2];
02916
02917 av[6] = tnv[1*3+0];
02918 av[7] = tnv[1*3+1];
02919 av[8] = tnv[1*3+2];
02920
02921 av[18] = tnv[2*3+0];
02922 av[19] = tnv[2*3+1];
02923 av[20] = tnv[2*3+2];
02924
02925 av[21] = tnv[3*3+0];
02926 av[22] = tnv[3*3+1];
02927 av[23] = tnv[3*3+2];
02928
02929 }
02930
02931 mxv (am,av,v);
02932 addv (v,nf,nf);
02933 }
02934
02935
02936
02937 if (is[4]>0 ){
02938 zeta=1.0;
02939 for (i=0;i<intordb;i++){
02940 xi=gp[i];
02941 for (j=0;j<intordb;j++){
02942 eta=gp[j];
02943
02944 jac2d_3d (jac,x,y,z,xi,eta,4);
02945 bf_matrix (n,xi,eta,zeta);
02946 jac = jac*w[i]*w[j];
02947 nnj (am.a,n.a,jac,n.m,n.n);
02948 }
02949 }
02950
02951 if (is[4]==1){
02952 av[0] = nv[48+0*3+0];
02953 av[1] = nv[48+0*3+1];
02954 av[2] = nv[48+0*3+2];
02955
02956 av[3] = nv[48+1*3+0];
02957 av[4] = nv[48+1*3+1];
02958 av[5] = nv[48+1*3+2];
02959
02960 av[6] = nv[48+2*3+0];
02961 av[7] = nv[48+2*3+1];
02962 av[8] = nv[48+2*3+2];
02963
02964 av[9] = nv[48+3*3+0];
02965 av[10] = nv[48+3*3+1];
02966 av[11] = nv[48+3*3+2];
02967 }
02968
02969
02970 if (is[4]==2){
02971
02972 av[0] = nv[48+0*3+0];
02973 av[1] = nv[48+0*3+1];
02974 av[2] = nv[48+0*3+2];
02975
02976 av[3] = nv[48+1*3+0];
02977 av[4] = nv[48+1*3+1];
02978 av[5] = nv[48+1*3+2];
02979
02980 av[6] = nv[48+2*3+0];
02981 av[7] = nv[48+2*3+1];
02982 av[8] = nv[48+2*3+2];
02983
02984 av[9] = nv[48+3*3+0];
02985 av[10] = nv[48+3*3+1];
02986 av[11] = nv[48+3*3+2];
02987
02988 locglob_nodeval (4,av,tnv,x,y,z);
02989
02990 av[0] = tnv[0*3+0];
02991 av[1] = tnv[0*3+1];
02992 av[2] = tnv[0*3+2];
02993
02994 av[3] = tnv[1*3+0];
02995 av[4] = tnv[1*3+1];
02996 av[5] = tnv[1*3+2];
02997
02998 av[6] = tnv[2*3+0];
02999 av[7] = tnv[2*3+1];
03000 av[8] = tnv[2*3+2];
03001
03002 av[9] = tnv[3*3+0];
03003 av[10] = tnv[3*3+1];
03004 av[11] = tnv[3*3+2];
03005
03006 }
03007
03008 mxv (am,av,v);
03009 addv (v,nf,nf);
03010 }
03011
03012
03013
03014
03015 if (is[5]>0 ){
03016 zeta=-1.0;
03017 for (i=0;i<intordb;i++){
03018 xi=gp[i];
03019 for (j=0;j<intordb;j++){
03020 eta=gp[j];
03021
03022 jac2d_3d (jac,x,y,z,xi,eta,5);
03023 bf_matrix (n,xi,eta,zeta);
03024 jac = jac*w[i]*w[j];
03025 nnj (am.a,n.a,jac,n.m,n.n);
03026 }
03027 }
03028
03029 if (is[5]==1){
03030 av[12] = nv[60+0*3+0];
03031 av[13] = nv[60+0*3+1];
03032 av[14] = nv[60+0*3+2];
03033
03034 av[21] = nv[60+1*3+0];
03035 av[22] = nv[60+1*3+1];
03036 av[23] = nv[60+1*3+2];
03037
03038 av[18] = nv[60+2*3+0];
03039 av[19] = nv[60+2*3+1];
03040 av[20] = nv[60+2*3+2];
03041
03042 av[15] = nv[60+3*3+0];
03043 av[16] = nv[60+3*3+1];
03044 av[17] = nv[60+3*3+2];
03045 }
03046
03047
03048 if (is[5]==2){
03049
03050 av[0] = nv[60+0*3+0];
03051 av[1] = nv[60+0*3+1];
03052 av[2] = nv[60+0*3+2];
03053
03054 av[3] = nv[60+1*3+0];
03055 av[4] = nv[60+1*3+1];
03056 av[5] = nv[60+1*3+2];
03057
03058 av[6] = nv[60+2*3+0];
03059 av[7] = nv[60+2*3+1];
03060 av[8] = nv[60+2*3+2];
03061
03062 av[9] = nv[60+3*3+0];
03063 av[10] = nv[60+3*3+1];
03064 av[11] = nv[60+3*3+2];
03065
03066 locglob_nodeval (5,av,tnv,x,y,z);
03067
03068 av[12] = tnv[0*3+0];
03069 av[13] = tnv[0*3+1];
03070 av[14] = tnv[0*3+2];
03071
03072 av[21] = tnv[1*3+0];
03073 av[22] = tnv[1*3+1];
03074 av[23] = tnv[1*3+2];
03075
03076 av[18] = tnv[2*3+0];
03077 av[19] = tnv[2*3+1];
03078 av[20] = tnv[2*3+2];
03079
03080 av[15] = tnv[3*3+0];
03081 av[16] = tnv[3*3+1];
03082 av[17] = tnv[3*3+2];
03083
03084 }
03085
03086 mxv (am,av,v);
03087 addv (v,nf,nf);
03088 }
03089
03090
03091
03092
03093
03094
03095 destrv (gp); destrv (w);
03096 delete [] tnv;
03097 }
03098
03099 void linhexrot::locglob_nodeval (long is,vector &nv,double *tnv,vector &x,vector &y,vector &z)
03100 {
03101 double norm;
03102 vector g1(3),g2(3),g3(3),lv(3),gv(3);
03103 matrix t(3,3);
03104
03105 if (is==0){
03106 g1[0]=x[4]-x[7];
03107 g1[1]=y[4]-y[7];
03108 g1[2]=z[4]-z[7];
03109
03110 g2[0]=x[3]-x[7];
03111 g2[1]=y[3]-y[7];
03112 g2[2]=z[3]-z[7];
03113 }
03114 if (is==1){
03115 g1[0]=x[5]-x[4];
03116 g1[1]=y[5]-y[4];
03117 g1[2]=z[5]-z[4];
03118
03119 g2[0]=x[0]-x[4];
03120 g2[1]=y[0]-y[4];
03121 g2[2]=z[0]-z[4];
03122 }
03123 if (is==2){
03124 g1[0]=x[6]-x[5];
03125 g1[1]=y[6]-y[5];
03126 g1[2]=z[6]-z[5];
03127
03128 g2[0]=x[1]-x[5];
03129 g2[1]=y[1]-y[5];
03130 g2[2]=z[1]-z[5];
03131 }
03132 if (is==3){
03133 g1[0]=x[7]-x[6];
03134 g1[1]=y[7]-y[6];
03135 g1[2]=z[7]-z[6];
03136
03137 g2[0]=x[2]-x[6];
03138 g2[1]=y[2]-y[6];
03139 g2[2]=z[2]-z[6];
03140 }
03141 if (is==4){
03142 g1[0]=x[3]-x[2];
03143 g1[1]=y[3]-y[2];
03144 g1[2]=z[3]-z[2];
03145
03146 g2[0]=x[1]-x[2];
03147 g2[1]=y[1]-y[2];
03148 g2[2]=z[1]-z[2];
03149 }
03150 if (is==5){
03151 g1[0]=x[5]-x[6];
03152 g1[1]=y[5]-y[6];
03153 g1[2]=z[5]-z[6];
03154
03155 g2[0]=x[7]-x[6];
03156 g2[1]=y[7]-y[6];
03157 g2[2]=z[7]-z[6];
03158 }
03159
03160
03161 scprd (g1,g1,norm);
03162 norm=1.0/sqrt(norm);
03163 cmulv (norm,g1,g1);
03164
03165 scprd (g1,g2,norm);
03166 g2[0]=g2[0]-norm*g1[0];
03167 g2[1]=g2[1]-norm*g1[1];
03168 g2[2]=g2[2]-norm*g1[2];
03169
03170 scprd (g2,g2,norm);
03171 norm=1.0/sqrt(norm);
03172 cmulv (norm,g2,g2);
03173
03174 g3[0]=g1[1]*g2[2]-g1[2]*g2[1];
03175 g3[1]=g1[2]*g2[0]-g1[0]*g2[2];
03176 g3[2]=g1[0]*g2[1]-g1[1]*g2[0];
03177
03178 t[0][0]=g1[0];
03179 t[1][0]=g1[1];
03180 t[2][0]=g1[2];
03181
03182 t[0][1]=g2[0];
03183 t[1][1]=g2[1];
03184 t[2][1]=g2[2];
03185
03186 t[0][2]=g3[0];
03187 t[1][2]=g3[1];
03188 t[2][2]=g3[2];
03189
03190 mxv (t.a,nv.a,tnv,3,3);
03191 mxv (t.a,nv.a+3,tnv+3,3,3);
03192 mxv (t.a,nv.a+6,tnv+6,3,3);
03193 mxv (t.a,nv.a+9,tnv+9,3,3);
03194
03195 }
03196
03197
03198
03199
03200
03201
03202
03203
03204
03205
03206
03207 void linhexrot::node_forces_surf_old (long lcid,long eid,long *is,double *nv,vector &nf)
03208 {
03209 long i,j,i1,ii;
03210 double xi=0.0,eta=0.0,zeta=0.0,jac, w1,w2;
03211 ivector nodes(nne);
03212 vector gx(nne),gy(nne),gz(nne),x(nnsurf),y(nnsurf),z(nnsurf),gp(intordb),w(intordb),av(ndofe),v(ndofe),pom(3);
03213 matrix n(napfun,ndofe),an(napfun,ndofe),am(ndofe,ndofe), tran(3,3);
03214
03215 Mt->give_elemnodes (eid,nodes);
03216 Mt->give_node_coord3d (gx,gy,gz,eid);
03217 gauss_points (gp.a,w.a,intordb);
03218 fillv (0.0,nf);
03219 fillm (0.0,an);
03220 for (ii=0;ii<6;ii++){
03221 fillv (0.0,av);
03222
03223 if (is[ii] ==0){}
03224 else {
03225 tran_mat( tran, gx, gy, gz, ii);
03226
03227 if (ii ==0){
03228 xi=1.0;
03229 x[0]=gx[0];x[1]=gx[3];x[2]=gx[7];x[3]=gx[4];
03230 y[0]=gy[0];y[1]=gy[3];y[2]=gy[7];y[3]=gy[4];
03231 z[0]=gz[0];z[1]=gz[3];z[2]=gz[7];z[3]=gz[4];
03232 for (i=0;i<intordb;i++){
03233 eta=gp[i]; w1=w[i];
03234 for (j=0;j<intordb;j++){
03235 zeta=gp[j]; w2=w[j];
03236 bf_matrix (n,xi,eta,zeta);
03237 jac2d3d (jac,x,y,z,eta,zeta);
03238 jac = jac*w1*w2;
03239 nnj (am.a,n.a,jac,n.m,n.n);
03240 }
03241 }
03242 for (i=0;i<3;i++){
03243 av[i]=nv[i]; av[9+i]=nv[3+i]; av[21+i]=nv[6+i]; av[12+i]=nv[9+i];
03244 }
03245
03246 if (is[ii] ==2){
03247 lgvectortransfblock (av, tran);
03248 }
03249 mxv (am,av,v);
03250 }
03251
03252 else if (ii ==1){
03253 eta=1.0;
03254 x[0]=gx[1];x[1]=gx[0];x[2]=gx[4];x[3]=gx[5];
03255 y[0]=gy[1];y[1]=gy[0];y[2]=gy[4];y[3]=gy[5];
03256 z[0]=gz[1];z[1]=gz[0];z[2]=gz[4];z[3]=gz[5];
03257 for (i=0;i<intordb;i++){
03258 xi=gp[i]; w1=w[i];
03259 for (j=0;j<intordb;j++){
03260 zeta=gp[j]; w2=w[j];
03261 bf_matrix (n,xi,eta,zeta);
03262 jac2d3d (jac,x,y,z,xi,zeta);
03263 jac = jac*w1*w2;
03264 nnj (am.a,n.a,jac,n.m,n.n);
03265 }
03266 }
03267 for (i=0;i<3;i++){
03268 av[3+i]=nv[12+i]; av[i]=nv[15+i]; av[12+i]=nv[18+i]; av[15+i]=nv[21+i];
03269 }
03270
03271 if (is[ii] ==2){
03272 lgvectortransfblock (av, tran);
03273 }
03274 mxv (am,av,v);
03275 }
03276
03277 else if (ii ==2){
03278 xi=-1.0;
03279 x[0]=gx[2];x[1]=gx[1];x[2]=gx[5];x[3]=gx[6];
03280 y[0]=gy[2];y[1]=gy[1];y[2]=gy[5];y[3]=gy[6];
03281 z[0]=gz[2];z[1]=gz[1];z[2]=gz[5];z[3]=gz[6];
03282 for (i=0;i<intordb;i++){
03283 eta=gp[i]; w1=w[i];
03284 for (j=0;j<intordb;j++){
03285 zeta=gp[j]; w2=w[j];
03286 bf_matrix (n,xi,eta,zeta);
03287 jac2d3d (jac,x,y,z,eta,zeta);
03288 jac = jac*w1*w2;
03289 nnj (am.a,n.a,jac,n.m,n.n);
03290 }
03291 }
03292 for (i=0;i<3;i++){
03293 av[6+i]=nv[24+i]; av[3+i]=nv[27+i]; av[15+i]=nv[30+i]; av[18+i]=nv[33+i];
03294 }
03295
03296 if (is[ii] ==2){
03297 lgvectortransfblock (av, tran);
03298 }
03299 mxv (am,av,v);
03300 }
03301
03302 else if (ii ==3){
03303 eta=-1.0;
03304 x[0]=gx[3];x[1]=gx[2];x[2]=gx[6];x[3]=gx[7];
03305 y[0]=gy[3];y[1]=gy[2];y[2]=gy[6];y[3]=gy[7];
03306 z[0]=gz[3];z[1]=gz[2];z[2]=gz[6];z[3]=gz[7];
03307 for (i=0;i<intordb;i++){
03308 xi=gp[i]; w1=w[i];
03309 for (j=0;j<intordb;j++){
03310 zeta=gp[j]; w2=w[j];
03311 bf_matrix (n,xi,eta,zeta);
03312 jac2d3d (jac,x,y,z,xi,zeta);
03313 jac = jac*w1*w2;
03314
03315 nnj (am.a,n.a,jac,n.m,n.n);
03316 }
03317 }
03318 for (i=0;i<3;i++){
03319 av[9+i]=nv[36+i]; av[6+i]=nv[39+i]; av[18+i]=nv[42+i]; av[21+i]=nv[45+i];
03320 }
03321
03322 if (is[ii] ==2){
03323 lgvectortransfblock (av, tran);
03324 }
03325 mxv (am,av,v);
03326 }
03327
03328 else if (ii ==4) {
03329 zeta=1.0;
03330 x[0]=gx[0];x[1]=gx[1];x[2]=gx[2];x[3]=gx[3];
03331 y[0]=gy[0];y[1]=gy[1];y[2]=gy[2];y[3]=gy[3];
03332 z[0]=gz[0];z[1]=gz[1];z[2]=gz[2];z[3]=gz[3];
03333 for (i=0;i<intordb;i++){
03334 xi=gp[i]; w1=w[i];
03335 for (j=0;j<intordb;j++){
03336 eta=gp[j]; w2=w[j];
03337 bf_matrix (n,xi,eta,zeta);
03338 jac2d3d (jac,x,y,z,xi,eta);
03339 jac = jac*w1*w2;
03340
03341
03342
03343 nnj (am.a,n.a,jac,n.m,n.n);
03344 }
03345 }
03346 for (i=0;i<3;i++){
03347 av[i]=nv[48+i]; av[3+i]=nv[51+i]; av[6+i]=nv[54+i]; av[9+i]=nv[57+i];
03348 }
03349
03350
03351
03352
03353
03354
03355
03356
03357 if (is[ii] ==2){
03358 lgvectortransfblock (av, tran);
03359 }
03360 mxv (am,av,v);
03361 }
03362
03363 else if (ii ==5) {
03364 zeta=-1.0;
03365 x[0]=gx[5];x[1]=gx[4];x[2]=gx[7];x[3]=gx[6];
03366 y[0]=gy[5];y[1]=gy[4];y[2]=gy[7];y[3]=gy[6];
03367 z[0]=gz[5];z[1]=gz[4];z[2]=gz[7];z[3]=gz[6];
03368 for (i=0;i<intordb;i++){
03369 xi=gp[i]; w1=w[i];
03370 for (j=0;j<intordb;j++){
03371 eta=gp[j]; w2=w[j];
03372 bf_matrix (n,xi,eta,zeta);
03373 jac2d3d (jac,x,y,z,xi,eta);
03374 jac = jac*w1*w2;
03375 nnj (am.a,n.a,jac,n.m,n.n);
03376 }
03377 }
03378 for (i=0;i<3;i++){
03379 av[12+i]=nv[60+i]; av[15+i]=nv[63+i]; av[18+i]=nv[66+i]; av[21+i]=nv[69+i];
03380 }
03381
03382 if (is[ii] ==2){
03383 lgvectortransfblock (av, tran);
03384 }
03385 mxv (am,av,v);
03386 }
03387
03388 fprintf (Out,"\n\n zatiz");
03389 for (i=0;i<nne;i++){
03390 i1=3*i;
03391 fprintf (Out,"\n %4ld %20.10e %20.10e %20.10e",i,av[i1],av[i1+1],av[i1+2]);}
03392 fprintf (Out,"\n\n vloc");
03393 for (i=0;i<nne;i++){
03394 i1=3*i;
03395 fprintf (Out,"\n %4ld %20.10e %20.10e %20.10e",i,v[i1],v[i1+1],v[i1+2]);}
03396
03397 fprintf (Out,"\n\n vglob");
03398 for (i=0;i<nne;i++){
03399 i1=3*i;
03400 fprintf (Out,"\n %4ld %20.10e %20.10e %20.10e",i,v[i1],v[i1+1],v[i1+2]);}
03401 addv (nf,v,nf);
03402 }
03403 }
03404
03405
03406
03407
03408
03409
03410
03411
03412
03413
03414
03415
03416 }
03417
03418
03419
03420
03421
03422
03423
03424
03425
03426
03427
03428 void linhexrot::tran_mat(matrix &tran, vector &gx, vector &gy, vector &gz, long is)
03429 {
03430 long i,i1;
03431 double dl;
03432 matrix a(3,3);
03433
03434 if (is ==4) {
03435 for (i=0; i<3; i++) {
03436 i1=i+1; if(i1>2)i1=i1-3;
03437 a[0][i]=gx[i1]-gx[i];
03438 a[1][i]=gy[i1]-gy[i];
03439 a[2][i]=gz[i1]-gz[i];
03440 }
03441 }
03442
03443 else if (is ==5) {
03444 for (i=4; i<7; i++) {
03445 i1=i+1; if(i1>6)i1=i1-3;
03446 a[0][i-4]=gx[i1]-gx[i];
03447 a[1][i-4]=gy[i1]-gy[i];
03448 a[2][i-4]=gz[i1]-gz[i];
03449 }
03450 }
03451
03452
03453
03454
03455 else {
03456 if(is>0)i1=is-4;else i1=is;
03457 a[0][0]=gx[i1+3]-gx[is];
03458 a[1][0]=gy[i1+3]-gy[is];
03459 a[2][0]=gz[i1+3]-gz[is];
03460 a[0][1]=gx[i1+7]-gx[i1+3];
03461 a[1][1]=gy[i1+7]-gy[i1+3];
03462 a[2][1]=gz[i1+7]-gz[i1+3];
03463 a[0][2]=gx[is]-gx[i1+7];
03464 a[1][2]=gy[is]-gy[i1+7];
03465 a[2][2]=gz[is]-gz[i1+7];
03466 }
03467
03468 dl=sqrt(a[0][0]*a[0][0]+a[1][0]*a[1][0]+a[2][0]*a[2][0]);
03469
03470 tran[0][0]=a[0][0]/dl;
03471 tran[1][0]=a[1][0]/dl;
03472 tran[2][0]=a[2][0]/dl;
03473
03474 tran[0][2]=a[1][0]*a[2][1]-a[2][0]*a[1][1];
03475 tran[1][2]=a[2][0]*a[0][1]-a[0][0]*a[2][1];
03476 tran[2][2]=a[0][0]*a[1][1]-a[1][0]*a[0][1];
03477
03478 dl=sqrt(tran[0][2]*tran[0][2]+tran[1][2]*tran[1][2]+tran[2][2]*tran[2][2]);
03479
03480 tran[0][2]=tran[0][2]/dl;
03481 tran[1][2]=tran[1][2]/dl;
03482 tran[2][2]=tran[2][2]/dl;
03483
03484 tran[0][1]=tran[1][2]*tran[2][0]-tran[2][2]*tran[1][0];
03485 tran[1][1]=tran[2][2]*tran[0][0]-tran[0][2]*tran[2][0];
03486 tran[2][1]=tran[0][2]*tran[1][0]-tran[1][2]*tran[0][0];
03487
03488 dl=sqrt(tran[0][1]*tran[0][1]+tran[1][1]*tran[1][1]+tran[2][1]*tran[2][1]);
03489
03490 tran[0][1]=tran[0][1]/dl;
03491 tran[1][1]=tran[1][1]/dl;
03492 tran[2][1]=tran[2][1]/dl;
03493
03494
03495
03496
03497
03498
03499
03500
03501
03502
03503
03504
03505
03506
03507
03508
03509
03510
03511 }
03512
03513
03514
03515
03516
03517
03518
03519
03520
03521 void linhexrot::intpointval (long eid,vector &nodval,vector &ipval)
03522 {
03523 long i,j,ii,jj,k,l;
03524 double xi,eta,zeta;
03525 vector w,gp;
03526
03527 l=0;
03528 for (ii=0;ii<nb;ii++){
03529 for (jj=0;jj<nb;jj++){
03530 if (intordsm[ii][jj]==0) continue;
03531
03532 allocv (intordsm[ii][jj],w);
03533 allocv (intordsm[ii][jj],gp);
03534
03535 gauss_points (gp.a,w.a,intordsm[ii][jj]);
03536
03537 for (i=0;i<intordsm[ii][jj];i++){
03538 xi=gp[i];
03539 for (j=0;j<intordsm[ii][jj];j++){
03540 eta=gp[j];
03541 for (k=0;k<intordsm[ii][jj];k++){
03542 zeta=gp[k];
03543
03544 ipval[l]=approx (xi,eta,zeta,nodval);
03545 l++;
03546 }
03547 }
03548 }
03549
03550 destrv (w);
03551 destrv (gp);
03552 }
03553 }
03554 }
03555
03556
03557
03558 void linhexrot::aver_strains (long lcid,long eid,long ri,long ci,vector &averstra,double &volume)
03559 {
03560 long i,j,k,l,ii,ipp;
03561 double xi,eta,zeta,jac;
03562 vector x(nne),y(nne),z(nne),gp,w,eps;
03563
03564 Mt->give_node_coord3d (x,y,z,eid);
03565
03566 for (ii=0;ii<nb;ii++){
03567 if (intordsm[ii][ii]==0) continue;
03568
03569 allocv (intordsm[ii][ii],gp);
03570 allocv (intordsm[ii][ii],w);
03571 allocv (ncomp[ii],eps);
03572
03573 gauss_points (gp.a,w.a,intordsm[ii][ii]);
03574
03575 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
03576 for (i=0;i<intordsm[ii][ii];i++){
03577 xi=gp[i];
03578 for (j=0;j<intordsm[ii][ii];j++){
03579 eta=gp[j];
03580 for (k=0;k<intordsm[ii][ii];k++){
03581 zeta=gp[k];
03582
03583
03584 Mm->givestress (lcid,ipp,eps);
03585
03586 jac_3d (jac,x,y,z,xi,eta,zeta);
03587 jac=fabs(jac);
03588
03589 volume+=w[i]*w[j]*w[k]*jac;
03590
03591 for (l=0;l<averstra.n;l++){
03592 averstra[l]+=eps[l]*w[i]*w[j]*w[k]*jac;
03593 }
03594
03595 ipp++;
03596 }
03597 }
03598 }
03599
03600
03601 destrv (w); destrv (gp); destrv (eps);
03602 }
03603
03604 }
03605
03606
03607
03608
03609
03610
03611
03612
03613
03614
03615
03616
03617
03618
03619 void linhexrot::surfnodeval (long surf,vector &nodval,double *list)
03620 {
03621 long i,j,k;
03622 ivector surfnod(nnsurf);
03623
03624 fillv (0.0,nodval);
03625 quadhexahedral_surfnod (surfnod.a,surf);
03626
03627 k=0;
03628 for (i=0;i<nnsurf;i++){
03629 for (j=0;j<napfun;j++){
03630 nodval[k]=list[surf*nnsurf*napfun+k];
03631 k++;
03632 }
03633 }
03634 }