elmt_shell_4n.c

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

                                  p->stiff->uMatrix.daa[j-1][i-1]);
                           printf("elmtNo = %d Stiffness matrix IS NOT SYMMETRIC \n",
                                  p->elmt_no);
                           break;
                       }
                       else {
                           if( i == size && j == size)
                              printf("elmtNo = %d Stiffness matrix IS SYMMETRIC \n",
                                      p->elmt_no);
                       }
                    }
                }
#endif
                MatrixFreeIndirectDouble(T_matrix, size);
                MatrixFreeIndirectDouble(T_Transpose, size);
                MatrixFreeIndirectDouble(stiff, size);

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

       /* Initiation of Stiffness Units Buffer                      */

       switch( UNITS_SWITCH ) {
         case ON:
           if(UNITS_TYPE == SI || UNITS_TYPE == SI_US ) {
              dp_stress = DefaultUnits("Pa");
              dp_length = DefaultUnits("m");
           }
           else {
              dp_stress = DefaultUnits("psi");
              dp_length = DefaultUnits("in");
           }

          /* node 1 */
           UnitsMultRep( &(p->stiff->spColUnits[0]), dp_stress, dp_length );
           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]) , dp_length );
           UnitsCopy( &(p->stiff->spColUnits[4]), &(p->stiff->spColUnits[3]) );

           ZeroUnits( &(p->stiff->spRowUnits[0]) );
           ZeroUnits( &(p->stiff->spRowUnits[1]) );
           ZeroUnits( &(p->stiff->spRowUnits[2]) );
           UnitsCopy( &(p->stiff->spRowUnits[3]), dp_length );
           UnitsCopy( &(p->stiff->spRowUnits[4]), dp_length );

          /* 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 *) dp_stress->units_name);
           free((char *) dp_stress);
           free((char *) dp_length->units_name);
           free((char *) dp_length);

          break;
          case OFF:
          break;
          default:
          break;
        }
             break;
	case PRESSLD:
	case EQUIV_NODAL_LOAD:

#ifdef DEBUG
       printf("*** In elmt_shell() : enter case EQUIV_NODAL_LOAD\n");
#endif
         /* Form external nodal load vector */
         /* due to distributed loading      */
            
       if(p->elmt_load_ptr != (ELEMENT_LOADS *) NULL) {
           p = sld04(p, PRESSLD); /* Equivalent Load in local_coordinate */
            
       /***************************************************/
       /* Transform nodal_load vector from local lamina   */
       /* coordinate system to global coordinate system   */
       /***************************************************/

       T_matrix    = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
       T_Transpose = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
       nodal_load  = MatrixAllocIndirectDouble(p->size_of_stiff, 1);

      /* store the direction matrix in T_matrix */
      /* and transfer vector p->nodal_load into */
      /* nodal_load as a matrix form            */

         for (i = 1; i <= p->size_of_stiff; i++)  {
             for(j = 1; j <= p->nodes_per_elmt; j++) {
                for( k = 1; k <= p->dof_per_node; k++) {
                   n1  = p->dof_per_node*(j-1);
                   n2  = p->dof_per_node*j;
                   n   = n1 + k;
                   if(i > n1 && i <= n2) {
                      ii = i - n1;
                      T_matrix[i-1][n-1] = Direction_Matrix[ii-1][k-1];
                   }
                   else
                      T_matrix[i-1][n-1] = 0.0;
                }
            }
         }

      /* =====================================*/
      /* Calculate the inverse of T_matrix    */
      /*     T_matrix inverse [T]^-1 =        */
      /*        = T_matrix trnspose : [T]^t   */
      /* =====================================*/

         for( i = 1; i <= p->size_of_stiff; i++) {
             for( j = 1; j <= p->size_of_stiff; j++) 
                T_Transpose[i-1][j-1] = T_matrix[j-1][i-1];
         }
	
         nodal_load = dMatrixMultRep(nodal_load,
                      T_Transpose, p->size_of_stiff,p->size_of_stiff,
                      p->equiv_nodal_load->uMatrix.daa, p->size_of_stiff,1);

         for(i = 1; i <= p->size_of_stiff; i++)
            p->equiv_nodal_load->uMatrix.daa[i-1][0]
            = nodal_load[i-1][0];

         MatrixFreeIndirectDouble(T_matrix, p->size_of_stiff);
         MatrixFreeIndirectDouble(T_Transpose, p->size_of_stiff);
         MatrixFreeIndirectDouble(nodal_load, p->size_of_stiff);
       }

/* ------------------------ UNITS -----------------------------------*/

      UNITS_SWITCH = CheckUnits();
      switch( UNITS_SWITCH ) {
        case ON:
         if(UNITS_TYPE == SI) {
            dp_length = DefaultUnits("m");
            dp_force  = DefaultUnits("N");
         }
         if(UNITS_TYPE == US) {
            dp_length = DefaultUnits("in");
            dp_force  = DefaultUnits("lbf");
         }
         dp_moment = UnitsMult( dp_force, dp_length );
             
         for(i= 1; i<= p->dof_per_node; i++) {
            for(j = 1; j <= p->nodes_per_elmt; j++) {
                k  = p->dof_per_node*(j-1)+i;
                k1 =  p->dof_per_node*j-3;
                if ( k <= k1) {      /* force units */
                   UnitsCopy( &(p->equiv_nodal_load->spRowUnits[k-1]), dp_force );
                }
                else {  /* k > k1 moment units */
                   UnitsCopy( &(p->equiv_nodal_load->spRowUnits[k-1]), dp_moment );
                }
             }
          }

          free((char *) dp_force->units_name);
          free((char *) dp_force);
          free((char *) dp_length->units_name);
          free((char *) dp_length);
          free((char *) dp_moment->units_name);
          free((char *) dp_moment);
      break;
      case OFF:
      break;
      default:
      break;
   }

#ifdef DEBUG
       printf("*** In elmt_shell() : print equiv_nodal_load\n\n");
       MatrixPrintIndirectDouble(p->equiv_nodal_load);
       printf("*** In elmt_shell() : end case EQUIV_NODAL_LOAD\n");
#endif

             break;

        case STRESS_UPDATE:

        /* update stress for given displacement */
        /* and displacement incremental         */
             break;
        case STRESS:
        case STRESS_LOAD:
        case LOAD_MATRIX:
        /* ====================================== */
        /* Form internal nodal load vector due to */
        /* stress at previous load step           */
        /* ====================================== */

#ifdef DEBUG
       printf("****** In elmt_shell_4nodes_implicit() : \n");
       printf("       enter cases: STRESS_UPDATE, LOAD_MATRIX, STRESS_LOAD \n");     
       printf(" elemt_state = %d \n", p->elmt_state);
#endif
          nodal_load      = MatrixAllocIndirectDouble(p->size_of_stiff, 1); 
          nodal_load_temp = MatrixAllocIndirectDouble(p->size_of_stiff, 1); 
          T_matrix        = MatrixAllocIndirectDouble(p->size_of_stiff, 
                                                      p->size_of_stiff);
          T_Transpose     = MatrixAllocIndirectDouble(p->size_of_stiff,
                                                      p->size_of_stiff);

          h  = p->work_section[11].value;    /* thickness of the shell */
          for(i = 1; i <= p->nodes_per_elmt; i++) {
              elmt_length = ABS(p->coord[0][i-1].value - p->coord[0][i].value);

              if(elmt_length != 0)
                 break;
           }

           aspect_ratio = h/elmt_length;
           k_bar = 2.0*(1+nu)*aspect_ratio*aspect_ratio; /* reduced shear factor */

/*         EE  = (p->mater_matrix->uMatrix.daa[0][0]
                 +p->mater_matrix->uMatrix.daa[1][1])*0.5; */

           EE  = E.value;

        /* =================================================*/
        /* nodal load due to the inital stress or stress at */
        /* previous step                                    */
        /* =================================================*/

      /* Integration loop over thickness */

          integ_coord = dVectorAlloc(no_integ_pts+1);
          weight      = dVectorAlloc(no_integ_pts+1);

          gauss(integ_coord,weight,no_integ_pts);

          for(ii = 1; ii <= no_integ_pts; ii++) {
              nodal_load_temp 
              = Shell_Nodal_Load_Plane(nodal_load_temp, p, co_coord,
                             integ_coord[ii], ii, E.value, nu, isw);
              for(i = 1; i<= p->size_of_stiff; i++) {
                  nodal_load[i-1][0] += nodal_load_temp[i-1][0]*weight[ii]*0.5*h;  
              }

          } /* gaussian integration ends */

          free((char *) integ_coord);
          free((char *) weight);

          for(i= 1; i<= p->size_of_stiff; i++) {
              p->nodal_loads[i-1].value = nodal_load[i-1][0];  
          }

       /***************************************************/
       /* Transform nodal_load vector from local lamina   */
       /* coordinate system to global coordinate system   */
       /***************************************************/

      /* =====================================*/
      /* Calculate the inverse of T_matrix    */
      /*     T_matrix inverse [T]^-1 =        */
      /*        = T_matrix trnspose : [T]^t   */
      /* =====================================*/

      /* store the direction matrix : T_matrix */

   if( isw == LOAD_MATRIX ) {
      for (i = 1; i <= p->size_of_stiff; i++)  {
           for(j = 1; j <= p->nodes_per_elmt; j++) {
               for( k = 1; k <= p->dof_per_node; k++) {
                    n1  = p->dof_per_node*(j-1);
                    n2  = p->dof_per_node*j;
                    n   = n1 + k;
                    if(i > n1 && i <= n2) {
                       ii = i - n1;
                       T_matrix[i-1][n-1] = Direction_Matrix[ii-1][k-1];
                    }
                    else 
                       T_matrix[i-1][n-1] = 0.0;
                }
           }
       }

      for( i = 1; i <= p->size_of_stiff; i++) {
         for( j = 1; j <= p->size_of_stiff; j++) 
           T_Transpose[i-1][j-1] = T_matrix[j-1][i-1];
      }
     
      size = p->size_of_stiff;	
      nodal_load_temp = dMatrixMultRep(nodal_load_temp,
                        T_Transpose, size, size, nodal_load, size,1);

      for(i= 1; i<= p->size_of_stiff; i++) {
          p->nodal_loads[i-1].value = nodal_load_temp[i-1][0];  
      }
    }

      MatrixFreeIndirectDouble(nodal_load,p->size_of_stiff);
      MatrixFreeIndirectDouble(nodal_load_temp,p->size_of_stiff);
      MatrixFreeIndirectDouble(T_Transpose,p->size_of_stiff);
      MatrixFreeIndirectDouble(T_matrix,p->size_of_stiff);


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

     switch( UNITS_SWITCH ) {
       case ON:
          if(UNITS_TYPE == SI || UNITS_TYPE == SI_US ) {
              dp_force  = DefaultUnits("N");
              dp_length = DefaultUnits("m");
          }
          else {
              dp_force  = DefaultUnits("lbf");
              dp_length = DefaultUnits("in");
          }

          /* node no 1 */
          UnitsCopy( p->nodal_loads[0].dimen, dp_force );
          UnitsCopy( p->nodal_loads[1].dimen, dp_force );
          UnitsCopy( p->nodal_loads[2].dimen, dp_force );
          UnitsMultRep( p->nodal_loads[3].dimen, dp_force, dp_length );
          UnitsCopy( p->nodal_loads[4].dimen, p->nodal_loads[3].dimen );

          /* 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->nodal_loads[k-1].dimen, p->nodal_loads[0].dimen );
                  else 
                     UnitsCopy( p->nodal_loads[k-1].dimen, p->nodal_loads[3].dimen );
              }
          }
          free((char *) dp_length->units_name);
          free((char *) dp_length);
          free((char *) dp_force->units_name);
          free((char *) dp_force);
          break;
       case OFF:
          break;
       default:
          break;
     }

#ifdef DEBUG
printf("*** In elmt_shell_4nodes_implicit() : end case LOAD_MATRIX, STRESS_UPDATE, STRESS_LOAD \n");
#endif
	     break;

        case MASS_MATRIX:  /* form mass matrix */