📄 elmt_psps.c
字号:
if(no_integ_pts*no_integ_pts != lint)
pgauss(no_integ_pts,&lint,sg,tg,wg);
/* start gaussian integration */
for(l = 1; l <= lint; l++) {
shape(sg[l-1],tg[l-1],p->coord,shp,&jacobain,p->no_dimen,p->nodes_per_elmt, p->node_connect,0);
/* output: shp[0][i-1] = dN_i/dx */
/* shp[1][i-1] = dN_i/dy */
/* compute [B] matrix at each Gaussian */
/* integration point */
/*****************************************************/
/* Derivative matrix (B_i)3x2; */
/* B_i[0][0] = dNi/dx, B_i[0][1] = 0 */
/* B_i[1][0] = 0, B_i[1][1] = dNi/dy */
/* B_i[2][0] = dNi/dy, B_i[2][2] = dNi/dx */
/* [B] = [B_1, B_2, B_3, B_4, ..., B_n] */
/* n = no of node */
/*****************************************************/
/* material matrix */
Mater_matrix = MATER_MAT_PLANE(E.value, nu); /* need to be changed */
/* Form [B] matrix */
for(j = 1; j <= p->nodes_per_elmt; j++) {
k = 2*(j-1);
B_matrix[0][k] = shp[0][j-1];
B_matrix[1][k+1] = shp[1][j-1];
B_matrix[2][k] = shp[1][j-1];
B_matrix[2][k+1] = shp[0][j-1];
}
/* muti. Jacobain determinant with weight coefficents and mater matrix */
jacobain = jacobain* wg[l-1];
/* [a] CALCULATE EQUIVALENT NODAL FORCE DUE TO INITIAL STRAIN */
if(p->nodal_init_strain != NULL) {
for( i = 1; i <=3 ; i++ )
for (j = 1; j <= 3; j++)
Mater_matrix[i-1][j-1] *= jacobain;
/* Transpose [B] matrix */
for(i = 1; i <= 3; i++)
for(j = 1; j <= dof_per_elmt; j++)
B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];
/* calculate strain[] at gaussian points : strain = sum Ni*nodal_strain */
for (i = 1; i <= 3; i++)
strain[i-1][0] = 0.0;
for (i = 1; i <= 3; i++) {
for( j = 1; j <= p->nodes_per_elmt; j++) {
strain[i-1][0] += shp[2][j-1] * p->nodal_init_strain[i-1][j-1];
}
}
/* [Mater]_3x3 * [strain]_3x1 and save in [stress]_3x1 */
stress = dMatrixMultRep(stress, Mater_matrix, 3, 3, strain, 3, 1);
/* mutiply [B]^T * [Mater]* [strain] */
if(nodal_load == NULL ) {
nodal_load = dMatrixMultRep(nodal_load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
} else {
load = dMatrixMultRep(load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
for (i = 1; i<= dof_per_elmt; i++) {
nodal_load[i-1][0] += load[i-1][0];
}
}
}
/* [b] CALCULATE EQUIVALENT NODAL FORCE DUE TO INITIAL STRESS */
if(p->nodal_init_stress != NULL) {
/* Transpose [B] matrix */
for(i = 1; i <= 3; i++)
for(j = 1; j <= dof_per_elmt; j++)
B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];
/* calculate stress[] at gaussian points : stress = sum Ni*nodal_stress */
for (i = 1; i <= 3; i++)
stress[i-1][0] = 0.0;
for (i = 1; i <= 3; i++) {
for( j = 1; j <= p->nodes_per_elmt; j++) {
stress[i-1][0] += shp[2][j-1] * p->nodal_init_stress[i-1][j-1].value;
}
stress[i-1][0] *= jacobain;
}
/* mutiply [B]^T * [stress] */
if(nodal_load == NULL) {
nodal_load = dMatrixMultRep(nodal_load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
} else {
load = dMatrixMultRep(load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
for (i = 1; i<= dof_per_elmt; i++)
nodal_load[i-1][0] -= load[i-1][0];
}
}
/* [c] CALCULATE EQUIVALENT NODAL FORCE DUE TO BODY FORCE */
if(p->nodal_body_force != NULL) {
/* calculate body_force[] at gaussian points : body_force = sum Ni*body_force */
for (i = 1; i <= p->no_dimen; i++)
body_force[i-1][0] = 0.0;
for(i = 1; i <= p->no_dimen; i++) {
for( j = 1; j <= p->nodes_per_elmt; j++) {
body_force[i-1][0]
+= shp[2][j-1] * p->nodal_body_force[i-1][j-1].value*jacobain;
/* mutiply [N]^T * [body_force] */
k = (j-1)*p->no_dimen + i;
if(nodal_load == NULL) {
nodal_load[k-1][0] = shp[2][j-1]*body_force[i-1][0];
} else
nodal_load[k-1][0] += shp[2][j-1]*body_force[i-1][0];
}
}
}
/* [d] CALCULATE EQUIVALENT NODAL FORCE DUE TO TEMPERATURE CHANGE */
if(p->nodal_init_strain != NULL) {
for(i = 1; i <= 3 ; i++ )
for(j = 1; j <= 3; j++)
Mater_matrix[i-1][j-1] *= jacobain;
/* Transpose [B] matrix */
for(i = 1; i <= 3; i++)
for(j = 1; j <= dof_per_elmt; j++)
B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];
/* [B]^T * [Mater] and save in [B]^T */
temp_matrix =
dMatrixMultRep(temp_matrix, B_Transpose, dof_per_elmt, 3, Mater_matrix, 3, 3);
/* calculate strain[] at gaussian points : strain = sum Ni*nodal_strain */
for(i = 1; i <= 3; i++)
strain[i-1][0] = 0.0;
for(j = 1; j <= p->nodes_per_elmt; j++) {
temperature += shp[2][j-1] * p->nodal_temp[j-1].value;
}
for(i = 1; i <= 2; i++) {
strain[i-1][0] = temperature *alpha_thermal[i-1];
}
strain[2][0] = 0.0;
/* mutiply [B]^T * [Mater]* [strain] */
if(nodal_load == NULL) {
nodal_load = dMatrixMultRep(nodal_load, temp_matrix, dof_per_elmt, 3, strain, 3, 1);
} else {
load = dMatrixMultRep(load, temp_matrix, dof_per_elmt, 3, strain, 3, 1);
for(i = 1; i<= dof_per_elmt; i++)
nodal_load[i-1][0] += load[i-1][0];
}
}
/* [e] CALCULATE EQUIVALENT NODAL FORCE DUE TO SURFACE TRACTION */
/* add code later */
/* Transfer nodal_load to p->equiv_nodal_load */
for(i = 1; i <= dof_per_elmt; i++) {
p->equiv_nodal_load->uMatrix.daa[i-1][0] += nodal_load[i-1][0];
}
}
/* ------------ EQUIVALENT NODAL LOAD UNITS ------------*/
/* Young's Modulus E is in SI then Use SI, else use US */
/* ---------------------------------------------------- */
switch(UNITS_SWITCH) {
case ON:
/* Initiation of Equivalent nodal load Units Buffer */
if(UNITS_TYPE == SI)
d1 = DefaultUnits("N");
else
d1 = DefaultUnits("lbf");
/* node 1 */
UnitsCopy( &(p->equiv_nodal_load->spRowUnits[0]), d1 );
UnitsCopy( &(p->equiv_nodal_load->spRowUnits[1]), d1 );
/* node i i > 1*/
for ( i = 2; i <= p->nodes_per_elmt; i++) {
for( j = 1; j <= p->dof_per_node; j++) {
k = p->dof_per_node*(i-1) + j;
UnitsCopy( &(p->equiv_nodal_load->spRowUnits[k-1]), d1 );
}
}
ZeroUnits( &(p->equiv_nodal_load->spColUnits[0]) );
free((char *) d1->units_name);
free((char *) d1);
break;
case OFF:
break;
default:
break;
}
#ifdef DEBUG
printf("*** in elmt_psps() : leaving case EQUIV_NODAL_LOAD: \n");
#endif
break;
case STRESS_UPDATE:
break;
case STRESS:
case LOAD_MATRIX:
#ifdef DEBUG
printf("*** in elmt_psps() : start case STRESS: or case LOAD_MATRIX: \n");
#endif
lint = (int ) 0;
if(isw == STRESS) l= no_stress_pts; /* stress pts */
if(isw == LOAD_MATRIX) l= no_integ_pts; /* guassian integ pts */
/* initilize nodal_load */
for (i = 1; i <= dof_per_elmt; i++) {
nodal_load[i-1][0] = 0.0;
load[i-1][0]= 0.0;
}
for (i = 1; i<= no_integ_pts*no_integ_pts; i++) {
sg[i-1] = 0.0;
tg[i-1] = 0.0;
wg[i-1] = 0.0;
}
if(l*l != lint)
pgauss(l,&lint,sg,tg,wg);
/* element stress, strains and forces */
if(p->no_dimen == 2 && isw == STRESS) {
printf("\n STRESS in Element No %d \n", p->elmt_no);
printf(" ================================================================================================== \n");
printf(" Gaussion xi eta x y Stress-11 Stress-22 Stress-12 \n");
if(UNITS_SWITCH == OFF)
printf(" Points \n");
if(UNITS_SWITCH == ON){
if(UNITS_TYPE == SI)
dp_stress = DefaultUnits("Pa");
else
dp_stress = DefaultUnits("psi");
printf(" Points %s %s %s %s %s\n",
p->coord[0][0].dimen->units_name,
p->coord[1][0].dimen->units_name,
dp_stress->units_name,
dp_stress->units_name,
dp_stress->units_name);
free((char *) dp_stress->units_name);
free((char *) dp_stress);
}
printf(" --------------------------------------------------------------------------------------------------\n \n");
}
#ifdef DEBUG
if(p->no_dimen == 2 && isw == LOAD_MATRIX) {
printf("\n NODAL LOAD in Element No %d \n", p->elmt_no);
printf(" ===================================================================================\n");
printf(" Gaussion xi eta x y nodal_no nodal_load \n");
printf(" Points \n");
printf(" ----------------------------------------------------------------------------------- \n");
}
#endif
for( l= 1; l<= lint; l++) {
/* element shape function and their derivatives */
shape(sg[l-1], tg[l-1],p->coord,shp,&jacobain,p->no_dimen,p->nodes_per_elmt, p->node_connect,0);
/* Form [B] matrix */
for(j = 1; j <= p->nodes_per_elmt; j++) {
k = 2*(j-1);
B_matrix[0][k] = shp[0][j-1];
B_matrix[1][k+1] = shp[1][j-1];
B_matrix[2][k] = shp[1][j-1];
B_matrix[2][k+1] = shp[0][j-1];
}
Mater_matrix = MATER_MAT_PLANE(E.value, nu);
/* calculate strains at guassian integretion pts */
for(i = 1;i <= 3; i++)
strain[i-1][0] = 0.0;
xx = 0.0;
yy = 0.0;
for(j = 1; j <= p->nodes_per_elmt; j++) {
xx = xx + shp[2][j-1] * p->coord[0][j-1].value;
yy = yy + shp[2][j-1] * p->coord[1][j-1].value;
/* converting p->displ into a array */
for ( k = 1; k <= p->dof_per_node; k++) {
j1 = p->dof_per_node*(j-1) + k;
displ[j1-1][0] = p->displ->uMatrix.daa[k-1][j-1];
}
}
strain = dMatrixMultRep(strain, B_matrix, 3, dof_per_elmt, displ, dof_per_elmt, 1);
/* compute stress */
stress = dMatrixMultRep(stress, Mater_matrix, 3, 3, strain, 3, 1);
if(isw == LOAD_MATRIX) {
for(i = 1; i <= 3; i++)
for(j = 1; j <= dof_per_elmt; j++)
B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -