📄 elmt_shell_4n_q.c
字号:
dv5 = d[5-1] * dv; dv6 = d[6-1] * dv; i1 = 1; for( i=1; i<=4; i++ ) { /* Recover previously computed shape functions */ shp1i = shp[0][i-1][k-1]; shp2i = shp[1][i-1][k-1]; shp3i = shp[2][i-1][k-1]; if( qdflg ) { shp11 = shp1[0][i-1][k-1]; shp12 = shp1[1][i-1][k-1]; shp13 = shp1[2][i-1][k-1]; shp21 = shp2[0][i-1][k-1]; shp22 = shp2[1][i-1][k-1]; shp23 = shp2[2][i-1][k-1]; a13 = -dv1*shp11 - dv2*shp22; a23 = -dv2*shp11 - dv1*shp22; a31 = dv3*shp2i; a32 = dv3*shp1i; a33 = -dv3*(shp12+shp21); } i2 = i1-1; /* Compute the loading term */ pp[i1-1] += shp3i*bl[0] *dv; pp[i1 ] += shp3i*bl[1] *dv; pp[i1+1] += shp3i*bl[2] *dv; pp[i1+2] += shp13*bl[2] *dv*tc; pp[i1+3] += shp23*bl[2] *dv*tc; pp[i1+4] -= (shp13*bl[0]+shp23*bl[1])*dv*tc; /* Form the stress-displacement matrix (Bi-trans*D) */ a11 = dv1*shp1i; a12 = dv2*shp2i; a21 = dv2*shp1i; a22 = dv1*shp2i; /* Form the plate stress-displacement matrix */ for( ii=ii1; ii<=ii2; ii++ ) { btd[0][ii-1] = dv4*bm[0][ii-1][i-1] + dv5*bm[1][ii-1][i-1]; btd[1][ii-1] = dv5*bm[0][ii-1][i-1] + dv4*bm[1][ii-1][i-1]; btd[2][ii-1] = dv6*bm[2][ii-1][i-1]; } /* Loop on columns */ j1 = i1; for( j=i; j<=4; j++ ) { j2 = j1 - 1; xn = shp[0][j-1][k-1]; yn = shp[1][j-1][k-1]; sn = shp[2][j-1][k-1]; if( qdflg ) { x1n = - shp1[0][j-1][k-1]; y2n = - shp2[1][j-1][k-1]; x2n = - shp1[1][j-1][k-1]; y1n = - shp2[0][j-1][k-1]; xyn = x2n + y1n; } /* Compute the membrane part */ s[i1-1][j1-1] += a11*xn + a31*yn; s[i1 ][j1-1] += a12*xn + a32*yn; s[i1-1][j1 ] += a21*yn + a31*xn; s[i1 ][j1 ] += a22*yn + a32*xn; s[i1+4][j1-1] += a13*xn + a33*yn; s[i1+4][j1 ] += a23*yn + a33*xn; s[i1-1][j1+4] += a11*x1n + a21*y2n + a31*xyn; s[i1 ][j1+4] += a12*x1n + a22*y2n + a32*xyn; s[i1+4][j1+4] += a13*x1n + a23*y2n + a33*xyn; /* Compute the bending part */ for( ii=ii1; ii<=ii2; ii++ ) { for( jj=ii1; jj<=ii2; jj++ ) { for( kk=1; kk<=3; kk++ ) s[i2+ii-1][j2+jj-1] += btd[kk-1][ii-1]*bm[kk-1][jj-1][j-1]; } } j1 = j1 + ndf; } i1 = i1 + ndf; } /* Compute Hughes/Brezzi rotation matrix */ if( ialph != 0 && qdflg ) { j1 = 0; for( j=1; j<=4; j++ ) { gg[j1 ] -= shp[1][j-1][k-1]*dv; gg[j1+1] += shp[0][j-1][k-1]*dv; gg[j1+5] -= 2.0*shp[2][j-1][k-1]*dv; if( ialph > 1 ) gg[j1+5] += (shp1[1][j-1][k-1] - shp2[0][j-1][k-1])*dv; j1 = j1 + ndf; } } } /* Perform the rank one update for the hughes/brezzi term */ if( ialph > 0 && qdflg ) { for( i=1; i<=nst; i++ ) { ggi = ggv * gg[i-1]; for( j=1; j<=nst; j++ ) s[i-1][j-1] += ggi * gg[j-1]; } } i1 = 1; if( yl[2][0] != 0.0 ) { for( i=1; i<=4; i++ ) { pp[i1+2] += yl[2][i-1]*pp[i1 ]; pp[i1+3] -= yl[2][i-1]*pp[i1-1]; j1 = i1; for( j=i; j<=4; j++ ) { proj06( i1-1, j1-1, s, yl[2][i-1], yl[2][j-1], nst); j1 = j1 + ndf; } i1 = i1 + ndf; } } rots06( s, pp, tr, nst, ndf ); /* Form residual */ ul = p->displ->uMatrix.daa; for( i=1; i<=nst; i++ ) { for( j=1; j<=nst; j++ ) pp[i-1] -= s[i-1][j-1]*ul[(j-1)/4][(j-1)%4]; } /* Transfer stiffness matrix/Nodal Load vector to ARRAY *p */ switch(isw) { case STIFF: p->stiff->uMatrix.daa = s; case LOAD_MATRIX: for( i=1; i<=nst; i++ ) p->equiv_nodal_load->uMatrix.daa[i-1][0] = pp[i-1]; break; default: break; } /* Assign Units to Stiffness Matrix */ if( CheckUnits() == ON && isw == STRESS ) { if( CheckUnitsType() == SI || CheckUnitsType() == SI_US ) { dimen1 = DefaultUnits("Pa"); dimen2 = DefaultUnits("m"); } else { dimen1 = DefaultUnits("psi"); dimen2 = DefaultUnits("in"); } /* node 1 */ UnitsMultRep( &(p->stiff->spColUnits[0]), dimen1, dimen2 ); UnitsCopy( &(p->stiff->spColUnits[1]), &(p->stiff->spColUnits[0]) ); UnitsCopy( &(p->stiff->spColUnits[2]), &(p->stiff->spColUnits[0]) ); UnitsMultRep( &(p->stiff->spColUnits[3]), &(p->stiff->spColUnits[0]), dimen2 ); UnitsCopy( &(p->stiff->spColUnits[4]), &(p->stiff->spColUnits[3]) ); UnitsCopy( &(p->stiff->spColUnits[5]), &(p->stiff->spColUnits[3]) ); ZeroUnits( &(p->stiff->spRowUnits[0]) ); ZeroUnits( &(p->stiff->spRowUnits[1]) ); ZeroUnits( &(p->stiff->spRowUnits[2]) ); UnitsCopy( &(p->stiff->spRowUnits[3]), dimen2 ); UnitsCopy( &(p->stiff->spRowUnits[4]), dimen2 ); UnitsCopy( &(p->stiff->spRowUnits[5]), dimen2 ); /* node i i > 1*/ for ( i = 2; i <= p->nodes_per_elmt; i++) { kk = p->dof_per_node*(i-1) + 3; for( j = 1; j <= p->dof_per_node; j++) { k = p->dof_per_node*(i-1) + j; if( k <= kk) { UnitsCopy( &(p->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[0]) ); UnitsCopy( &(p->stiff->spRowUnits[k-1]), &(p->stiff->spRowUnits[0]) ); } if(k > kk) { UnitsCopy( &(p->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[3]) ); UnitsCopy( &(p->stiff->spRowUnits[k-1]), &(p->stiff->spRowUnits[3]) ); } } } free((char *) dimen1->units_name); free((char *) dimen1); free((char *) dimen2->units_name); free((char *) dimen2); } break; case STRESS: /* Compute and output the element variables */ /* Allocate working memory */ if(vl == NULL) vl = MatrixAllocIndirectDouble( 6, 4 ); if(eps == NULL) eps = dVectorAlloc(6); if(sigi == NULL) sigi = dVectorAlloc(6); if(sigb == NULL) sigb = dVectorAlloc(6); ul = p->displ->uMatrix.daa; for( i=1; i<=4; i++ ) { for( j=1; j<=3; j++ ) { vl[j-1][i-1] = 0.0; vl[j+2][i-1] = 0.0; for( k=1; k<=3; k++ ) { vl[j-1][i-1] += tr[j-1][k-1]*ul[k-1][i-1]; vl[j+2][i-1] += tr[j-1][k-1]*ul[k+2][i-1]; } } } /* Loop over Gauss Integration Points */ l = d[15]; if( l*l != lint ) pgauss ( l, &lint, sg, tg, wg ); for( k=1; k<=lint; k++ ) { rshp06( k-1, sg[k-1], tg[k-1], yl, &xsj, ndm); /* Modify the rotational shape functions */ for( i=1; i<=3; i++ ) { for( j=1; j<=4; j++ ) { shp1[i-1][j-1][0] -= gshp1[i-1][j-1]/dv; shp2[i-1][j-1][0] -= gshp2[i-1][j-1]/dv; } } stre06(p,d,vl, ndm, nel,1, &xx, &yy, &zz, eps, sigi, sigb); if( CheckUnitsType() == SI || CheckUnitsType() == SI_US ) { dimen1 = DefaultUnits("N"); dimen2 = DefaultUnits("m"); } else { dimen1 = DefaultUnits("lbf"); dimen2 = DefaultUnits("in"); } coord = dimen2->scale_factor; stress = dimen1->scale_factor/dimen2->scale_factor; moment = dimen1->scale_factor; curvtr = 1.0/dimen2->scale_factor; free((char *) dimen1->units_name); free((char *) dimen1); free((char *) dimen2->units_name); free((char *) dimen2); printf("\n\n STRESS in Element No. %d, Element: %s; Material: %s", p->elmt_no, p->elmt_type, p->material_name); printf("\n =================================================================================="); if( CheckUnitsType() == SI ) printf("\n coord = m, stress = N/m, moment = N*m/m, curvtr = 1/m, angle = rad"); else printf("\n coord = in, stress = lbf/in, moment = lbf*in/in, curvtr = 1/in, angle = rad"); printf("\n\n x-coord xx-stress xy-stress yy-stress 1-stress 2-stress angle"); printf("\n y-coord xx-moment xy-moment yy-moment 1-moment 2-moment angle"); printf("\n z-coord xx-strain xy-strain yy-strain xx-curvtr xy-curvtr yy-curvtr"); printf("\n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"); printf("\n %8.3lf%12.3e%12.3e%12.3e%12.3e%12.3e%10.2lf", xx/coord,sigi[0]/stress,sigi[1]/stress,sigi[2]/stress,sigi[3]/stress,sigi[4]/stress,sigi[5]); printf("\n %8.3lf%12.3e%12.3e%12.3e%12.3e%12.3e%10.2lf", yy/coord,sigb[0]/moment,sigb[1]/moment,sigb[2]/moment,sigb[3]/moment,sigb[4]/moment,sigb[5]); printf("\n %8.3lf%12.3e%12.3e%12.3e%12.3e%12.3e%12.3e\n\n", zz/coord,eps[0],eps[1],eps[2],eps[3]/curvtr,eps[4]/curvtr,eps[5]/curvtr); } break; case MASS_MATRIX: /* Compute element mass array */ break; default: break; } /* [a] : Free working memory *//* MatrixFreeIndirectDouble(yl, 3); MatrixFreeIndirectDouble(tr, 3); MatrixFreeIndirectDouble(gshp1, 3); MatrixFreeIndirectDouble(gshp2, 3); MatrixFreeIndirectDouble(vl, 6); MatrixFreeIndirectDouble(btd, 3); MatrixFreeIndirectDouble(s, nst); free ((char *) eps); free ((char *) sigi); free ((char *) sigb); free ((char *) sg); free ((char *) tg); free ((char *) wg); free ((char *) gg); free ((char *) dvl); free ((char *) pp); free ((char *) d); free ((char *) bl);*/ return(p);} /* * ====================================================================== * dktb06( ) : * * Input : * Output : * ====================================================================== */void dktb06( sg, tg,xsj )double sg, tg,xsj;{double *el, **shn, *bd, *cd;double shm1, shm2;double *a1, *b1, *c1, *d1, *e1;double *a2, *b2, *c2, *d2, *e2;int i, j, k; #ifdef DEBUG printf("*** Entering dktb06() ***\n");#endif /* [a] : Allocate memory for working matrices */ el = dVectorAlloc(3); bd = dVectorAlloc(3); cd = dVectorAlloc(3); a1 = dVectorAlloc(3); a2 = dVectorAlloc(3); b1 = dVectorAlloc(3); b2 = dVectorAlloc(3); c1 = dVectorAlloc(3); c2 = dVectorAlloc(3); d1 = dVectorAlloc(3); d2 = dVectorAlloc(3); e1 = dVectorAlloc(3); e2 = dVectorAlloc(3); shn = MatrixAllocIndirectDouble( 2, 3 ); /* [b] : Compute the area coordinates for the point */ el[0] = 0.25 * (1.0 - sg) * (1.0 - tg); el[1] = 0.25 * (1.0 + sg) * (1.0 - tg); el[2] = 0.5 * (1.0 + tg); /* [c] : Form the shape function terms for trains */ for( i=1; i<=3; i++ ) { bd[i-1] = b[i-1]/xsj; cd[i-1] = c[i-1]/xsj; } for( i=1; i<=3; i++ ) { j = i%3 + 1; k = j%3 + 1; shn[0][i-1] = bd[i-1] * (el[i-1] - 1.0); shn[1][i-1] = cd[i-1] * (el[i-1] - 1.0); shm1 = bd[j-1]*el[k-1] + bd[k-1]*el[j-1]; shm2 = cd[j-1]*el[k-1] + cd[k-1]*el[j-1]; a1[i-1] = aa[i-1] * shm1; a2[i-1] = aa[i-1] * shm2; b1[i-1] = bb[i-1] * shm1; b2[i-1] = bb[i-1] * shm2; c1[i-1] = cc[i-1] * shm1; c2[i-1] = cc[i-1] * shm2; d1[i-1] = dd[i-1] * shm1; d2[i-1] = dd[i-1] * shm2; e1[i-1] = ee[i-1] * shm1; e2[i-1] = ee[i-1] * shm2; } /* [d] : Plate train displacement matrix */ for( i=1; i<=3; i++ ) { j = i%3 + 1; k = j%3 + 1; bm[0][2][i-1] = a1[k-1] - a1[j-1]; bm[0][3][i-1] = b1[k-1] + b1[j-1]; bm[0][4][i-1] = c1[k-1] + c1[j-1] - shn[0][i-1]; bm[1][2][i-1] = d2[k-1] - d2[j-1]; bm[1][3][i-1] =-e2[k-1] + e2[j-1] + shn[1][i-1]; bm[1][4][i-1] =-b2[k-1] + b2[j-1]; bm[2][2][i-1] = a2[k-1] - a2[j-1] + d1[k-1] - d1[j-1]; bm[2][3][i-1] =-e1[k-1] - e1[j-1] + shn[0][i-1] - bm[1][4][i-1]; bm[2][4][i-1] = c2[k-1] + c2[j-1] - shn[1][i-1] - bm[0][3][i-1]; } /* [e] : Free working memory */ MatrixFreeIndirectDouble(shn,2); free ((char *) el); free ((char *) bd); free ((char *) cd); free ((char *) a1); free ((char *) a2); free ((char *) b1); free ((char *) b2);⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -