fe_print.c

有限元程序

                          array->displ->uMatrix.daa[k-1][i-1]
                          = frame->node[node_no -1].disp[ii-1].value;
                          if( UNITS_SWITCH == ON ) {
                              UnitsCopy( &(array->displ->spRowUnits[k-1]),
                                         frame->node[node_no -1].disp[ii-1].dimen);
                              ZeroUnits( &(array->displ->spColUnits[i-1]) );
                          }
                        }
                        k = k + 1;
                        break;
                  default:
                        break;
               }
           }
        }

        array = elmlib(array, STRESS);

        for( i=1 ; i <= frame->no_dof ; i++ )
        for( j=1 ; j <= frame->no_nodes_per_elmt ; j++ )
            frame->element[elmt_no-1].rp->Forces->uMatrix.daa[i-1][j-1]
            = array->nodal_loads[frame->no_dof*(j-1)+i-1].value;

        if( UNITS_SWITCH == ON ) {
            for( i=1 ; i <= frame->no_dof ; i++ )
                 UnitsCopy( &(frame->element[elmt_no-1].rp->Forces->spRowUnits[i-1]),
                              array->nodal_loads[i-1].dimen );
            for( j=1 ; j <= frame->no_nodes_per_elmt ; j++ )
                 ZeroUnits( &(frame->element[elmt_no-1].rp->Forces->spColUnits[j-1]) );
        }
    }

    /* 
     *  ====================================================================
     *  Print stresses corresponding to "nonlinear" response of structure
     *  These stresses are stored internally in Aladdin database       
     *                                                                 
     *  Aladdin command is : PrintStress() with no matrix argument.
     *                                                                 
     *  Note : Code in this section still needs to be checked          
     *  ====================================================================
     */ 

    if(m1 == (MATRIX *) NULL) {

       slp = lookup("InPlaneIntegPts");  /* no of integration pts in plane/surface */
       if(slp == NULL)
          iInPlaneIntegPts = 2*2;        /* 2x2 as default */
       else
          iInPlaneIntegPts = (int) slp->u.q->value;

       slp = lookup("ThicknessIntegPts"); /* no of integration pts in thickness direction */
       if(slp == NULL)
          iThicknessIntegPts = 2;        /* 2 as default */
       else
          iThicknessIntegPts = (int) slp->u.q->value;
  
       iNO_INTEG_pts = iInPlaneIntegPts*iThicknessIntegPts;

       lint = (int) 0;
       l = (int) iInPlaneIntegPts;
       for(i = 1; i<= l*l; i++) {
           xi[i-1]  = 0.0;
           eta[i-1] = 0.0;
           wg[i-1]  = 0.0;
       }

       if(l*l != lint)
          pgauss(l,&lint,xi,eta,wg);

       gamma = dVectorAlloc(iThicknessIntegPts+1);
       wt    = dVectorAlloc(iThicknessIntegPts+1);
       gauss(gamma,wt,iThicknessIntegPts);

       for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
           ep           = &frame->element[elmt_no-1];
           elmt_attr_no = ep->elmt_attr_no;
           eap          = &frame->eattr[elmt_attr_no-1];
           if(elmt_no == 1) 
              printf(" Element : %s \n Material : %s \n\n", eap->elmt_type, eap->material);

           printf("\n STRESS in  Element No  %d \n",elmt_no);
           printf(" =============================================================================================================== \n");
           printf(" Gaussion    xi       eta        gamma    Stre-xx         Stre-yy         Stre-xy         Stre-yz        Stre-zx \n");
           if(UNITS_SWITCH == OFF)
              printf("  Points \n");
           if(UNITS_SWITCH == ON)   {
              printf("  Points                                    %s             %s             %s              %s             %s \n",
                          ep->rp->stress->spRowUnits[0].units_name,
                          ep->rp->stress->spRowUnits[1].units_name,
                          ep->rp->stress->spRowUnits[2].units_name,
                          ep->rp->stress->spRowUnits[3].units_name,
                          ep->rp->stress->spRowUnits[4].units_name);
              printf(" ---------------------------------------------------------------------------------------------------------------\n \n");

              for(i = 1; i <= iThicknessIntegPts; i++)
              for(j = 1; j <= lint; j++) {
                  k = lint*(i-1)+j;
                  printf("   %d  %10.4f %10.4f %10.4f", k, xi[j-1], eta[j-1], gamma[i]);
                  printf("\t%11.4e\t%11.4e\t%11.4e\t%11.4e\t%11.4e\n", 
                             ep->rp->stress->uMatrix.daa[0][k-1]/ep->rp->stress->spRowUnits[0].scale_factor,
                             ep->rp->stress->uMatrix.daa[1][k-1]/ep->rp->stress->spRowUnits[1].scale_factor,
                             ep->rp->stress->uMatrix.daa[2][k-1]/ep->rp->stress->spRowUnits[2].scale_factor,
                             ep->rp->stress->uMatrix.daa[3][k-1]/ep->rp->stress->spRowUnits[3].scale_factor,
                             ep->rp->stress->uMatrix.daa[4][k-1]/ep->rp->stress->spRowUnits[4].scale_factor);
              }
              }
              if(UNITS_SWITCH == OFF)
                for(i = 1; i <= iThicknessIntegPts; i++)
                    for(j = 1; j <= lint; j++) {
                        k = lint*(i-1)+j;
                        printf("   %d  %10.4f %10.4f %10.4f", k, xi[j-1], eta[j-1], gamma[i]);
                        printf("\t%11.4e\t%11.4e\t%11.4e\t%11.4e\t%11.4e\n",
                             ep->rp->stress->uMatrix.daa[0][k-1],
                             ep->rp->stress->uMatrix.daa[1][k-1],
                             ep->rp->stress->uMatrix.daa[2][k-1],
                             ep->rp->stress->uMatrix.daa[3][k-1],
                             ep->rp->stress->uMatrix.daa[4][k-1]);
                }
        }
    }
}

/* 
 *  =========================================================
 *  Print_Displ() : Print Displacements
 *
 *  Input  : Matrix * : Matrix of Displacements to be printed
 *  Output : void
 *  =========================================================
 */ 

#ifdef __STDC__
void Print_Displ(MATRIX *m1)
#else
void Print_Displ(m1)
MATRIX *m1;
#endif
{
int node_no, j, k, ii, jj;
int          UNITS_SWITCH;
double                 da;

    UNITS_SWITCH = CheckUnits();

    if(frame->no_dof == 6) {
    printf("\n");
    printf("==================================================================================================================\n");
    printf(" Node                                             Displacement\n");
    printf("   No            displ-x           displ-y           displ-z             rot-x             rot-y             rot-z\n");
    printf("==================================================================================================================\n");

    if( UNITS_SWITCH == ON ) {
        printf(" units");
        for( ii=1 ; ii<=6 ; ii++ ) {
           if( ii <= 3 ) {
               if((CheckUnitsType()==SI) || (CheckUnitsType()==SI_US) )
                   printf("                m ");
               else 
                   printf("               in ");
           } else
                   printf("              rad ");
        }
    }
    printf("\n");
    }

    if(frame->no_dof == 5) {
    printf("\n");
    printf("================================================================================================\n");
    printf(" Node                                   Displacement\n");
    printf("   No            displ-x           displ-y           displ-z             rot-x             rot-y\n");
    printf("================================================================================================\n");

    if( UNITS_SWITCH == ON ) {
          printf(" units");
          for( ii=1 ; ii<=5 ; ii++ ) {
             if( ii <= 3 ) {
                if( (CheckUnitsType()==SI) || (CheckUnitsType()==SI_US) )
                   printf("                m ");
                else 
                   printf("               in ");
             }
             else
                   printf("              rad ");
          }
    }
    printf("\n");
    }

    if(frame->no_dof == 3) {
    printf("\n");
    printf("============================================================\n");
    printf(" Node                           Displacement                \n");
    printf("   No            displ-x           displ-y             rot-z\n");
    printf("============================================================\n");

    if( UNITS_SWITCH == ON ) {
        printf(" units");
        for( ii=1 ; ii<=3 ; ii++ ) {
             if( ii <= 2 ) {
                if( (CheckUnitsType()==SI) || (CheckUnitsType()==SI_US) )
                   printf("                m ");
                else 
                   printf("               in ");
             }
             else
                   printf("              rad ");
          }
    }
    printf("\n");
    }

    if(frame->no_dof == 2) {
    printf("\n");
    printf("==========================================\n");
    printf(" Node                  Displacement       \n");
    printf("   No            displ-x           displ-y\n");
    printf("==========================================\n");

    if( UNITS_SWITCH == ON ) {
        printf(" units");
        for( ii=1; ii <= 2; ii++ ) {
             if( (CheckUnitsType()==SI) || (CheckUnitsType()==SI_US) )
                 printf("                m ");
             else 
                 printf("               in ");
        }
    }
    printf("\n");
    } 

 /*
  *  ==================================================
  *  Print displacements for UNITS_SWITCH == ON and OFF
  *  ==================================================
  */

  if( UNITS_SWITCH == ON ) {

      for(node_no = 1; node_no <= frame->no_nodes; node_no++) {

      printf("%4d  ", node_no);
      for(j  = 1; j <= frame->no_dof; j++){

          jj = frame->node[node_no - 1].bound_id[j-1];
          if(jj > 0) {
             if(m1->spRowUnits[jj-1].units_name != NULL ) {
                UnitsTypeConvert(&(m1->spRowUnits[jj-1]), CheckUnitsType());
             }
             RadUnitsSimplify( &(m1->spRowUnits[jj-1]) );
             frame->node[node_no-1].disp[j-1].value = m1->uMatrix.daa[jj-1][0];
             UnitsCopy( frame->node[node_no-1].disp[j-1].dimen, &(m1->spRowUnits[jj-1]) );
             da =  MatrixContentScaleIndirectDouble(m1,jj, 1);
             printf("     %13.5e", da);
          } else {

             switch(CheckUnitsType()) {
                case SI:
                     if( (frame->node[node_no-1].disp[j-1].dimen->units_name != NULL) &&
                        !strcmp( frame->node[node_no-1].disp[j-1].dimen->units_name, "deg_F"))

                     frame->node[node_no -1].disp[j-1].value
                     = ConvertTempUnits(frame->node[node_no-1].disp[j-1].dimen->units_name,
                                        frame->node[node_no -1].disp[j-1].value, US);
                     break;
                case US:
                     if((frame->node[node_no-1].disp[j-1].dimen->units_name != NULL) &&
                         !strcmp(frame->node[node_no-1].disp[j-1].dimen->units_name, "deg_C"))

                         frame->node[node_no -1].disp[j-1].value
                         = ConvertTempUnits(frame->node[node_no-1].disp[j-1].dimen->units_name,
                                            frame->node[node_no -1].disp[j-1].value, SI);
                     break;
             }

             RadUnitsSimplify( frame->node[node_no-1].disp[j-1].dimen );
             printf("     %13.5e", frame->node[node_no -1].disp[j-1].value/
                          frame->node[node_no-1].disp[j-1].dimen->scale_factor);
         }
      }
      printf("\n");
      }
  }

  /* Print displacements when units are OFF */

  if( UNITS_SWITCH == OFF ) {
   
      for(node_no = 1; node_no <= frame->no_nodes; node_no++) {
          printf("%4d  ", node_no);
          for(j  = 1; j <= frame->no_dof; j++){
              jj = frame->node[node_no - 1].bound_id[j-1];
              if(jj > 0) {
                 frame->node[node_no-1].disp[j-1].value = m1->uMatrix.daa[jj-1][0];
                 printf("     %13.5e", m1->uMatrix.daa[jj-1][0]);
              }
              else
                 printf("     %13.5e", frame->node[node_no -1].disp[j-1].value);
          }
          printf("\n");
      }
  }

}