elmt_shell_4n_q.c

利用语言编写的有限元分析软件

           dv1 = d[1-1] * dv;
           dv2 = d[2-1] * dv;
           dv3 = d[3-1] * dv;
           dv4 = d[4-1] * dv;
           dv5 = d[5-1] * dv;
           dv6 = d[6-1] * dv;
 
           i1 = 1;
           for( i=1; i<=4; i++ ) {

           /* Recover the 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 8 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 the 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 );

/*
#ifdef DEBUG
    dMatrixPrint("s -- Stiffness Matrix", s, nst, nst);
    printf("pp vector\n");
    for(i=1; i<=nst; i++)
       printf("%5d%15.8e\n", i, pp[i-1]);
    dMatrixPrint("p->displ->uMatrix.daa", p->displ->uMatrix.daa, p->displ->iNoRows, p->displ->iNoColumns);
#endif
*/

/* ---- 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];
        }
/*
    printf("pp vector\n");
    for(i=1; i<=nst; i++)
    printf("%5d%15.8e\n", i, pp[i-1]);
*/
            

/* Transfer stiff matrix and equiv. nodal load vector to ARRAY *p */

    switch(isw) {
      case STIFF:
        p->stiff->uMatrix.daa = s;

    /**************************************************/
    /* Assign Units to Stiffness Matrix               */
    /**************************************************/

       if( CheckUnits() == ON ) {
          if(UNITS_TYPE == SI || UNITS_TYPE == 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 LOAD_MATRIX:
         for( i=1; i<=nst; i++ ) 
           p->equiv_nodal_load->uMatrix.daa[i-1][0] = pp[i-1];
      break;
      default:
      break;
    }

/*
#ifdef DEBUG
    dMatrixPrint("p->equiv_nodal_load->uMatrix.daa", p->equiv_nodal_load->uMatrix.daa, 
                   p->equiv_nodal_load->iNoRows, p->equiv_nodal_load->iNoColumns);
#endif
*/

#ifdef DEBUG
    printf(" In elmt_shell_4nodes_q(): \n");
    printf("    : leaving case of STIFF\n");
#endif
        break;

        case STRESS:     /* Compute and output the element variables */

#ifdef DEBUG
        printf("*** elmt_shell_4nodes_q() : In case STRESS\n");
#endif

        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]; 
              }
           }
        }

        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( UNITS_TYPE == SI || UNITS_TYPE == 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("                       S h e l l   S t r e s s e s");
           printf("\n elmt,x-coord  xx-stress  xy-stress  yy-stress   1-stress   2-stress   angle");
           printf("\n matl,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");
           if( UNITS_TYPE == 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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -");
           printf("\n     %8.3lf%11.3e%11.3e%11.3e%11.3e%11.3e%8.2lf",
                   xx/coord,sigi[0]/stress,sigi[1]/stress,sigi[2]/stress,sigi[3]/stress,sigi[4]/stress,sigi[5]);
           printf("\n     %8.3lf%11.3e%11.3e%11.3e%11.3e%11.3e%8.2lf",
                   yy/coord,sigb[0]/moment,sigb[1]/moment,sigb[2]/moment,sigb[3]/moment,sigb[4]/moment,sigb[5]);
           printf("\n     %8.3lf%11.3e%11.3e%11.3e%11.3e%11.3e%11.3e\n\n",
                   zz/coord,eps[0],eps[1],eps[2],eps[3]/curvtr,eps[4]/curvtr,eps[5]/curvtr);
        }
        
#ifdef DEBUG
    printf(" In elmt_shell_4nodes_q(): \n");
    printf("    : leaving case of STRESS\n");
#endif
        break;

        case MASS_MATRIX:                   /* Compute element mass array */

#ifdef DEBUG
    printf(" In elmt_shell_4nodes_q(): \n");
    printf("    : entering case of MASS_MATRIX\n");
#endif
 
        for(i=1; i<=nst; i++)
           pp[i-1] = 0.0;

        for( k=1; k<=lint; k++ ) {
           dv1 = dvl[k-1]*d[12];
           i1 = 0;
           for( j=1; j<=4; j++ ) {
              for( i=1; i<=3; i++ ) {
                 pp[i1+i-1] += shp[2][j-1][k-1]*dv1;
                 s[i1+i-1][i1+i-1] = pp[i1+i-1];
              }
              i1 = i1 + ndf;
           }
        }
	p->stiff->uMatrix.daa = s;

 /* ------------ MASS UNITS ---------------------------- */
 /* The units type is determined by the SetUnitsType()   */
 /* ---------------------------------------------------- */

 /* Initiation of Mass Units Buffer                      */

      switch( CheckUnits() ) {
        case ON:
           if(UNITS_TYPE == SI || UNITS_TYPE == SI_US ) {
              dimen1 = DefaultUnits("Pa");
              dimen2 = DefaultUnits("m");
           }
           else {
              dimen1 = DefaultUnits("psi");
              dimen2 = DefaultUnits("in");
           }
           dimen3 = DefaultUnits("sec");

          /* node no 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]) );

           UnitsPowerRep( &(p->stiff->spRowUnits[0]), dimen3, 2.0, NO );
           UnitsCopy( &(p->stiff->spRowUnits[1]), &(p->stiff->spRowUnits[0]) );
           UnitsCopy( &(p->stiff->spRowUnits[2]), &(p->stiff->spRowUnits[0]) );
           UnitsMultRep( &(p->stiff->spRowUnits[3]), dimen2, &(p->stiff->spRowUnits[0]) );
           UnitsCopy( &(p->stiff->spRowUnits[4]), &(p->stiff->spRowUnits[3]) );
           UnitsCopy( &(p->stiff->spRowUnits[5]), &(p->stiff->spRowUnits[3]) );

           /* node no > 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;
                   if(j <= 3) {
                      UnitsCopy( &(p->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[0]) );
                      UnitsCopy( &(p->stiff->spRowUnits[k-1]), &(p->stiff->spRowUnits[0]) );
                   }
                   else {
                      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);
           free((char *) dimen3->units_name);
           free((char *) dimen3);

           break;
        case OFF:
           break;
        default:
           break;
     }