fe_profile.c

有限元分析源代码

                      p->coord[i-1][j-1].value 
                      = frp->node[node_no -1].coord[i-1].value;
                      if(UNITS_SWITCH == ON )
                         UnitsCopy( p->coord[i-1][j-1].dimen,
                         frp->node[node_no -1].coord[i-1].dimen );
                  }
               }
            }

            if(task == STRESS_LOAD) {
/*
               Get_elmt_forces_and_strains(p,frp,elmt_no);
*/
               ;
            }
            else{                     /* PRESSLD */
               i = 0;
               for(j=1;j <= frp->no_element_loads; j++) {
                   if(elmt_no == frp->eforces[i].elmt_no ){ 
                      p->elmt_load_ptr =  &frp->eforces[i-1] ;
                      goto PRESS_END;
                   }
                   else p->elmt_load_ptr = (ELEMENT_LOADS *) NULL;
                   i = i+1;
               }
            }

            PRESS_END:
            break;
       case STRESS:

            /* For nonlinear problems, get element forces and strains */

            for(i = 1; i <= p->size_of_stiff; i++) {
                p->nodal_loads[i-1].value = 0.0;
                if( UNITS_SWITCH == ON )
                    ZeroUnits( p->nodal_loads[i-1].dimen );
            }

	    /* Element Displacements */

            if( UNITS_SWITCH == ON ) {
                for(i = 1; i <= p->dof_per_node; i++)
                    ZeroUnits( &(p->displ->spRowUnits[i-1]) );
                for(i = 1; i <= p->nodes_per_elmt; i++)
                    ZeroUnits( &(p->displ->spColUnits[i-1]) );
            }

            /* Zero and Transfer Element Displacements */

            for(i=1; i <= p->dof_per_node; i++) {
            for(j=1; j <= p->nodes_per_elmt; j++) {
                p->displ->uMatrix.daa[i-1][j-1]      = el->rp->displ->uMatrix.daa[i-1][j-1];
                p->displ_incr->uMatrix.daa[i-1][j-1] = 0.0;
            }
            }

            /* Element Coordinates */

            for(j = 1; j <= p->nodes_per_elmt; j++) {
                node_no = el->node_connect[j-1];
                if(node_no != 0){
                   for(i=1; i <= p->no_dimen; i++) {
                       p->coord[i-1][j-1].value 
                       = frp->node[node_no -1].coord[i-1].value;
                       if(UNITS_SWITCH == ON ) {
                          UnitsCopy( p->coord[i-1][j-1].dimen,
                          frp->node[node_no -1].coord[i-1].dimen );
                       }
                   }
                }
            }

            /* [g] : Get Elmt Load Array */

            i = 0;
            for(j = 1; j <= frp->no_element_loads; j++) {
                if(elmt_no == frp->eforces[i-1].elmt_no) { 
                   p->elmt_load_ptr = &frp->eforces[i-1];
                   goto  TRESS_END;
                }
                else 
                   p->elmt_load_ptr = (ELEMENT_LOADS *) NULL;
                i = i+1;
            }

            /* [h] : Element stress & strain  */

            for(j = 1; j <= iNO_INTEG_pts; j++) {

            p->effect_pl_strain[j-1]   = el->rp->effect_pl_strain[j-1];
            p->eff_pl_strain_incr[j-1] = 0.0;

            for(i = 1; i <= 9; i++) {
                p->stress->uMatrix.daa[i-1][j-1]         = el->rp->stress->uMatrix.daa[i-1][j-1];
                p->strain_pl->uMatrix.daa[i-1][j-1]      = el->rp->strain_pl->uMatrix.daa[i-1][j-1];
                p->strain_pl_incr->uMatrix.daa[i-1][j-1] = 0.0;
            }
            }
            break;

            TRESS_END:
            break;
       default:
            break; 
    }
    return(p);

#ifdef DEBUG
       printf("*** Leave Assign_p_Array()\n");
#endif
}


/* 
 *  ===================================================
 *  Alloc_p_Array() : Allocate memory for working array
 *  ===================================================
 */ 

ARRAY *Alloc_p_Array()
{
ARRAY               *pp;
int   i,j, no_integ_pts;
int                temp;
int        UNITS_SWITCH;

SYMBOL             *slp;
int    iInPlaneIntegPts;
int  iThicknessIntegPts;
int       iNO_INTEG_pts;

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

   slp = lookup("ThicknessIntegPts"); /* number 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;

   slp = lookup("NDimension");         /* No of dimensions : 2 or 3 */
   max_no_dimen = (int) slp->u.q->value;

   slp = lookup("NDofPerNode");        /* No gdof per node */
   max_no_dof = (int) slp->u.q->value;

   slp = lookup("MaxNodesPerElement"); /* Max no nodes per element */

   max_nodes_per_elmt = (int) slp->u.q->value;
   max_size_of_stiff  = max_nodes_per_elmt*max_no_dof;
   no_integ_pts       = iThicknessIntegPts;

   UNITS_SWITCH = CheckUnits();
   pp = (ARRAY *) MyCalloc(1,sizeof(ARRAY));

   pp->integ_ptr                = (INTEG_PTS *) MyCalloc(1,sizeof(INTEG_PTS));
   pp->integ_ptr->surface_pts   = (int) sqrt((double)iInPlaneIntegPts);  
   pp->integ_ptr->thickness_pts =  iThicknessIntegPts;
   pp->integ_ptr->integ_pts     =  no_integ_pts;

   pp->elmt_type      = (char *) NULL;
   pp->material_name  = (char *) NULL;
   pp->displ          = MatrixAllocIndirect("displ",DOUBLE_ARRAY,
                        max_no_dof,max_nodes_per_elmt);
   pp->work_section   = (QUANTITY *) MyCalloc(UNIT_SECTION_ATTR,sizeof(QUANTITY)); 
   pp->work_material  = (QUANTITY *) MyCalloc(UNIT_MATERIAL_ATTR, sizeof(QUANTITY));
   pp->d_array        = iVectorAlloc(max_size_of_stiff); 
   pp->coord          = (QUANTITY **) MyCalloc(max_no_dimen, sizeof(QUANTITY *));

   for (i = 1; i <= max_no_dimen; i++)
       pp->coord[i-1] = (QUANTITY *) MyCalloc(max_nodes_per_elmt, sizeof(QUANTITY));

   pp->node_connect   = iVectorAlloc(max_nodes_per_elmt);  
   pp->nodal_loads    = (QUANTITY *) MyCalloc(max_size_of_stiff, sizeof(QUANTITY));
   pp->nodal_temp     = (QUANTITY *) MyCalloc(max_no_dof, sizeof(QUANTITY)); 
   pp->nodal_body_force  = (QUANTITY **) MyCalloc(max_no_dof, sizeof(QUANTITY *));   

   for( i = 1; i <= max_no_dof; i++)
       pp->nodal_body_force[i-1] = (QUANTITY *) MyCalloc(max_nodes_per_elmt, sizeof(QUANTITY));

   pp->nodal_init_strain = (double **)MatrixAllocIndirectDouble(max_no_dof, max_nodes_per_elmt); 
   pp->nodal_init_stress = (QUANTITY **) MyCalloc(max_no_dof, sizeof(QUANTITY *));   

   for( i = 1; i <= max_no_dof; i++)
       pp->nodal_init_stress[i-1] = (QUANTITY *) MyCalloc(max_nodes_per_elmt, sizeof(QUANTITY));

   pp->nodal_traction    = (QUANTITY **) MyCalloc(max_no_dimen, sizeof(QUANTITY *));   

   for( i = 1; i <= max_no_dimen; i++)
        pp->nodal_traction [i-1] = (QUANTITY *)
                                   MyCalloc(max_nodes_per_elmt, sizeof(QUANTITY));

   pp->stiff            = MatrixAllocIndirect("stiff/mass",
                          DOUBLE_ARRAY, max_size_of_stiff, max_size_of_stiff);
   pp->equiv_nodal_load = MatrixAllocIndirect("equivalent_load",
                          DOUBLE_ARRAY,max_size_of_stiff, 1);
   pp->mater_matrix     = MatrixAllocIndirect("material_matrix",
                          DOUBLE_ARRAY,max_no_dof, max_no_dof);
   pp->displ_incr       = MatrixAllocIndirect("incremt_displ",
                          DOUBLE_ARRAY, max_no_dof, max_nodes_per_elmt);
   pp->strain_rate      = MatrixAllocIndirect("strain_rate", DOUBLE_ARRAY, 9, 1);
   pp->stress_rate      = MatrixAllocIndirect("stress_rate", DOUBLE_ARRAY, 9, 1);

   pp->effect_pl_strain   = (double *) MyCalloc(iNO_INTEG_pts, sizeof(double));
   pp->eff_pl_strain_incr = (double *) MyCalloc(iNO_INTEG_pts, sizeof(double));
   pp->direc_cos        = (double *) MyCalloc(2, sizeof(double));
   pp->direc_cos[0]     = 1.0;
   pp->direc_cos[1]     = 0.0;

   pp->LC_ptr           = (MATER_LOAD_CURVE *) MyCalloc(1,sizeof(MATER_LOAD_CURVE));
   pp->LC_ptr->name     = (char *)NULL;
   pp->LC_ptr->R        = (double *) MyCalloc(iNO_INTEG_pts, sizeof(double));
   pp->LC_ptr->H        = (double *) MyCalloc(iNO_INTEG_pts, sizeof(double));
   pp->LC_ptr->back_stress = MatrixAllocIndirectDouble(6, iNO_INTEG_pts);

   pp->stress           = MatrixAllocIndirect("stress_rate", DOUBLE_ARRAY, 9, iNO_INTEG_pts);
   pp->strain_pl        = MatrixAllocIndirect("strain", DOUBLE_ARRAY, 9, iNO_INTEG_pts);
   pp->strain_pl_incr   = MatrixAllocIndirect("strain", DOUBLE_ARRAY, 9, iNO_INTEG_pts);

/* element load pointer is not allocated until the load is present */

   if( UNITS_SWITCH == ON ) {

      for(i = 1; i <= UNIT_SECTION_ATTR; i++) {
          pp->work_section[i-1].dimen 
          = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS)); 
      }

      for(i = 1; i <= UNIT_MATERIAL_ATTR; i++) {
          pp->work_material[i-1].dimen
          = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS)); 
      }

      for (i = 1; i <= max_no_dimen; i++) {
          for (j = 1; j <= max_nodes_per_elmt; j++) {
             pp->coord[i-1][j-1].dimen 
             = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
          }
      }

      for (i = 1; i <= max_size_of_stiff; i++) {
         pp->nodal_loads[i-1].dimen 
         = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
      }

      for (i = 1; i <= max_no_dof; i++) {
         pp->nodal_temp[i-1].dimen 
         = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
      }

      for( i = 1; i <= max_no_dof; i++) {
           for (j = 1; j <= max_nodes_per_elmt; j++) {
              pp->nodal_body_force[i-1][j-1].dimen 
              = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
           }
      }

      for( i = 1; i <= max_no_dof; i++) {
          for (j = 1; j <= max_nodes_per_elmt; j++) {
             pp->nodal_init_stress[i-1][j-1].dimen 
             = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
          }
      }

      for( i = 1; i <= max_no_dimen; i++) {
          for (j = 1; j <= max_nodes_per_elmt; j++) {
             pp->nodal_traction[i-1][j-1].dimen 
             = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
          }
      }

      pp->eangle.dimen         = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
      pp->ealpha.dimen         = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
      pp->length.dimen         = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
   }
   return(pp);

#ifdef DEBUG
       printf("*** Leave fe_profile.c : Alloc_p_Array()\n\n");
#endif
}


/* ==================================================== */
/* Get_elmt_forces_and_strains(p,frp,elmt_no)           */
/* ==================================================== */

#ifdef __STDC__
Get_elmt_forces_and_strains(ARRAY *p, EFRAME *frp, int elmt_no)
#else
Get_elmt_forces_and_strains(p,frp,elmt_no)
ARRAY  *p;
EFRAME *frp;
int    elmt_no; 
#endif
{
 
int i, j, k, no_dof, no_nodes_per_elmt;
char                   *elmt_type_name;
ELEMENT_ATTR                      *eap;
int                       elmt_attr_no;
int                             length;
int                       UNITS_SWITCH;

    UNITS_SWITCH = CheckUnits();
    p->elmt_state = frp->element[elmt_no -1].esp->state;
/*
    p->eep        = frp->element[elmt_no -1].ep;
    p->ealpha     = frp->element[elmt_no -1].alpha;
*/
    
    elmt_attr_no = frp->element[elmt_no-1].elmt_attr_no;
    eap = &frp->eattr[elmt_attr_no-1];

    /* get forces from EFRAME *    */
    
    for (i = 1; i <= NO_ELEMENTS_IN_LIBRARY; i++) {
        if(!strcmp(elmt_library[i-1].name, frp->eattr[elmt_attr_no-1].elmt_type)) {
        elmt_type_name    = elmt_library[i-1].name;
        no_dof            = elmt_library[i-1].no_dof;
        no_nodes_per_elmt = elmt_library[i-1].no_node_per_elmt;
        }
    }
    
    for(i = 1; i <= no_dof+1; i++){ 
       for(j = 1; j <= no_nodes_per_elmt; j++) {

          k = no_dof*(j-1)+i;
          p->nodal_loads[k-1].value   = frp->element[elmt_no -1].rp->Forces->uMatrix.daa[i-1][j-1]; 
          if( UNITS_SWITCH == ON ) {
              UnitsCopy( p->nodal_loads[k-1].dimen,
              &(frp->element[elmt_no -1].rp->Forces->spRowUnits[i-1]) );
          }
       }
    }
}