00001 #include <stdlib.h>
00002 #include <math.h>
00003 #include "quadwedge.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
00015
00016 quadwedge::quadwedge (void)
00017 {
00018 long i;
00019
00020
00021 nne=15;
00022
00023 ndofe=45;
00024
00025 tncomp=6;
00026
00027 napfun=3;
00028
00029 intordmm=1;
00030
00031 ned=9;
00032
00033 nned=3;
00034
00035 nsurf=5;
00036
00037 nnsurf=8;
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 intordsmt = new long* [nb];
00056 intordsmz = new long* [nb];
00057 for (i=0;i<nb;i++){
00058 nip[i] = new long [nb];
00059 intordsmt[i] = new long [nb];
00060 intordsmz[i] = new long [nb];
00061 }
00062
00063 nip[0][0]=9;
00064
00065
00066 tnip=9;
00067
00068
00069 intordsmt[0][0]=3;
00070 intordsmz[0][0]=3;
00071
00072 }
00073
00074 quadwedge::~quadwedge (void)
00075 {
00076 long i;
00077
00078 for (i=0;i<nb;i++){
00079 delete [] nip[i];
00080 delete [] intordsmt[i];
00081 delete [] intordsmz[i];
00082 }
00083 delete [] nip;
00084 delete [] intordsmt;
00085 delete [] intordsmz;
00086
00087 delete [] cncomp;
00088 delete [] ncomp;
00089 }
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 double quadwedge::approx (double xi,double eta,double zeta,vector &nodval)
00101 {
00102 double f;
00103 vector bf(nne);
00104
00105 bf_quad_wed_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
00120 void quadwedge::bf_matrix (matrix &n,double xi,double eta,double zeta)
00121 {
00122 long i,j,k,l;
00123 vector bf(nne);
00124
00125 fillm (0.0,n);
00126
00127 bf_quad_wed_3d (bf.a,xi,eta,zeta);
00128
00129 j=0; k=1; l=2;
00130 for (i=0;i<nne;i++){
00131 n[0][j]=bf[i]; j+=3;
00132 n[1][k]=bf[i]; k+=3;
00133 n[2][l]=bf[i]; l+=3;
00134 }
00135
00136 }
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 void quadwedge::geom_matrix (matrix &gm,vector &x,vector &y,vector &z,
00151 double xi,double eta,double zeta,double &jac)
00152 {
00153 long i,j,k,l;
00154 vector dx(nne),dy(nne),dz(nne);
00155
00156 dx_bf_quad_wed_3d (dx.a,xi,eta,zeta);
00157 dy_bf_quad_wed_3d (dy.a,xi,eta,zeta);
00158 dz_bf_quad_wed_3d (dz.a,xi,eta,zeta);
00159
00160 derivatives_3d (dx,dy,dz,jac,x,y,z,xi,eta,zeta);
00161
00162 fillm (0.0,gm);
00163
00164 j=0; k=1; l=2;
00165 for (i=0;i<nne;i++){
00166 gm[0][j]=dx[i];
00167 gm[1][k]=dy[i];
00168 gm[2][l]=dz[i];
00169
00170 gm[3][k]=dz[i];
00171 gm[3][l]=dy[i];
00172
00173 gm[4][j]=dz[i];
00174 gm[4][l]=dx[i];
00175
00176 gm[5][j]=dy[i];
00177 gm[5][k]=dx[i];
00178
00179 j+=3; k+=3; l+=3;
00180 }
00181
00182 }
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193 void quadwedge::transf_matrix (ivector &nodes,matrix &tmat)
00194 {
00195 long i,n,m;
00196
00197 fillm (0.0,tmat);
00198
00199 n=nodes.n;
00200 m=tmat.m;
00201 for (i=0;i<m;i++){
00202 tmat[i][i]=1.0;
00203 }
00204
00205 for (i=0;i<n;i++){
00206 if (Mt->nodes[nodes[i]].transf>0){
00207 tmat[i*3+0][i*3]=Mt->nodes[nodes[i]].e1[0];
00208 tmat[i*3+1][i*3]=Mt->nodes[nodes[i]].e1[1];
00209 tmat[i*3+2][i*3]=Mt->nodes[nodes[i]].e1[2];
00210
00211 tmat[i*3+0][i*3+1]=Mt->nodes[nodes[i]].e2[0];
00212 tmat[i*3+1][i*3+1]=Mt->nodes[nodes[i]].e2[1];
00213 tmat[i*3+2][i*3+1]=Mt->nodes[nodes[i]].e2[2];
00214
00215 tmat[i*3+0][i*3+2]=Mt->nodes[nodes[i]].e3[0];
00216 tmat[i*3+1][i*3+2]=Mt->nodes[nodes[i]].e3[1];
00217 tmat[i*3+2][i*3+2]=Mt->nodes[nodes[i]].e3[2];
00218 }
00219 }
00220 }
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232 void quadwedge::stiffness_matrix (long eid,long ri,long ci,matrix &sm)
00233 {
00234 long i,k,ii,jj,ipp,transf;
00235 double xi,eta,zeta,jac;
00236 vector x(nne),y(nne),z(nne),w,wt,gp,gp1,gp2;
00237 matrix gm,d(tncomp,tncomp);
00238
00239 Mt->give_node_coord3d (x,y,z,eid);
00240
00241 fillm (0.0,sm);
00242
00243 for (ii=0;ii<nb;ii++){
00244 allocm (ncomp[ii],ndofe,gm);
00245 for (jj=0;jj<nb;jj++){
00246 if (intordsmt[ii][jj]==0) continue;
00247
00248 allocv (intordsmz[ii][jj],w);
00249 allocv (intordsmz[ii][jj],gp);
00250 allocv (intordsmt[ii][jj],wt);
00251 allocv (intordsmt[ii][jj],gp1);
00252 allocv (intordsmt[ii][jj],gp2);
00253
00254 gauss_points (gp.a,w.a,intordsmz[ii][jj]);
00255 gauss_points_tr (gp1.a,gp2.a,wt.a,intordsmt[ii][jj]);
00256
00257 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00258
00259 for (i=0;i<intordsmt[ii][jj];i++){
00260 xi=gp1[i]; eta=gp2[i];
00261 for (k=0;k<intordsmz[ii][jj];k++){
00262 zeta=gp[k];
00263
00264
00265 geom_matrix (gm,x,y,z,xi,eta,zeta,jac);
00266
00267 Mm->matstiff (d,ipp); ipp++;
00268
00269 jac*=wt[i]*w[k];
00270
00271
00272 bdbjac (sm,gm,d,gm,jac);
00273
00274 }
00275 }
00276 destrv (gp); destrv (w);
00277 }
00278 destrm (gm);
00279 }
00280
00281
00282 ivector nodes (nne);
00283 Mt->give_elemnodes (eid,nodes);
00284 transf = Mt->locsystems (nodes);
00285 if (transf>0){
00286 matrix tmat (ndofe,ndofe);
00287 transf_matrix (nodes,tmat);
00288 glmatrixtransf (sm,tmat);
00289 }
00290 }
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 void quadwedge::res_stiffness_matrix (long eid,matrix &sm)
00301 {
00302 stiffness_matrix (eid,0,0,sm);
00303 }
00304