fe_mesh.c

有限元程序

/*
 *  ============================================================================= 
 *  ALADDIN Version 2.1.
 *                                                                     
 *  fe_mesh.c : Generate Finite Element Mesh
 *                                                                     
 *  Copyright (C) 1995-2000 by Mark Austin, Xiaoguang Chen, and Wane-Jang Lin
 *  Institute for Systems Research,                                           
 *  University of Maryland, College Park, MD 20742                                   
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such and must not
 *     be misrepresented as being the original software.
 *  4. This software may not be sold or included in commercial software
 *     products without a license. 
 *  5. This notice is to remain intact.
 *                                                                    
 *  Written by: Mark Austin and Wane-Jang Lin                            May 1997
 *  ============================================================================= 
 */

#include <stdio.h>
#include <math.h>
#ifdef  __STDC__
#include <stdarg.h>
#else
#include <varargs.h>
#endif

#include "defs.h"
#include "miscellaneous.h"
#include "units.h"
#include "matrix.h"
#include "fe_database.h"
#include "symbol.h"
#include "vector.h"
#include "fe_functions.h"
#include "elmt.h"

/* External Declarations for global frame/working element data structures */

extern ARRAY     *array;
extern EFRAME    *frame;

int    TDOF;
int    TNEQ;
int    MDOF;

/*
#define DEBUG 
*/

/*
 *  -------------------------------------
 *  Add Node to Finite Element Mesh
 *  -------------------------------------
 */
#ifdef __STDC__
void Add_Node(double node_no, MATRIX *m)
#else
void Add_Node(node_no, m)
double node_no;
MATRIX      *m; /* coordinate vector/matrix */
#endif
{
int no = 0, num = 0;
int     i, j, length;
int     UNITS_SWITCH;
NODE             *np;

#ifdef DEBUG
       printf("*** Enter Add_Node() : node_no = %4d :\n ", (int) node_no);
#endif 

    UNITS_SWITCH = CheckUnits();
    frame = CheckNodeSpace(frame, (int) node_no);

    switch( UNITS_SWITCH ) {
      case ON:
         for (i =1; i <= frame->no_dimen; i++) {
            frame->node[(int) (node_no-1)].coord[i-1].value 
            = m->uMatrix.daa[0][i-1];
            UnitsCopy( frame->node[(int)(node_no-1)].coord[i-1].dimen, &(m->spColUnits[i-1]) );
         }
         break;
      case OFF:
         for(i =1; i <= frame->no_dimen; i++) 
             frame->node[(int) (node_no-1)].coord[i-1].value
             = m->uMatrix.daa[0][i-1];
         break;
      default:
         break;
    }
    frame->no_nodes = (int) MAX(node_no, frame->no_nodes); 

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

}


/*
 *  -------------------------------------
 *  Add Element to Finite Element Mesh
 *  -------------------------------------
 */

#ifdef __STDC__
void Add_Elmt(double elmt_no, MATRIX *nodes, char *elmt_attr_name)
#else
void Add_Elmt(elmt_no, nodes, elmt_attr_name)
double        elmt_no;
MATRIX         *nodes;
char  *elmt_attr_name;
#endif
{
int             i, no;
char            *type;

#ifdef DEBUG
       printf("\n*** Enter Add_Elmt(%3d) : nodal connectivity = [ ", (int) elmt_no);
       for(i = 1; i <= (int) nodes->iNoColumns; i++) 
           printf("%2d ", (int) nodes->uMatrix.daa[0][i-1]);
       printf(" ]\n");
       printf("          Add_Elmt(%3d) : elmt attribute  = \"%s\"\n", (int) elmt_no, (char *) elmt_attr_name);
       printf("\n");
#endif 


   frame = CheckElementSpace(frame, (int) elmt_no);

   frame->element[(int) (elmt_no-1)].elmt_attr_name
   = SaveString(elmt_attr_name);

   for(i = 1; i <= nodes->iNoColumns; i++)  {
       frame->element[(int)(elmt_no-1)].node_connect[i-1] 
       = (int) nodes->uMatrix.daa[0][i-1];
   }

   frame->no_elements = (int) MAX(elmt_no, frame->no_elements); 

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


/*
 *  ----------------------------------------------------
 *  Add Element and Section Attributes to Frame Database
 *  ----------------------------------------------------
 */

#ifdef __STDC__
int  Add_Element_Attr(char *name, ELEMENT_ATTR *ep)
#else
int  Add_Element_Attr(name, ep)
char          *name;
ELEMENT_ATTR    *ep;
#endif
{
SYMBOL          *hp;
int               i;

#ifdef DEBUG
       printf("*** Enter Add_Element_Attr() : name = \"%s\"\n", name);
#endif 

     /* [a] : Install elment attribute name into Symbol Table */

     hp = install(name);
     if(hp == NULL)
        FatalError("In Add_Element_Attr() : hp = NULL !!!\n",(char *)NULL);
     else
        hp->u.eap = ep;

     frame->no_element_attr = frame->no_element_attr + 1;

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

}

#ifdef __STDC__
int  Add_Fiber_Elmt_Attr(char *name, FIBER_ELMT *fep)
#else
int  Add_Fiber_Elmt_Attr(name, fep)
char          *name;
FIBER_ELMT    *fep;
#endif
{
SYMBOL          *hp;

#ifdef DEBUG
       printf("*** Enter Add_Fiber_Elmt_Attr() : name = \"%s\"\n", name);
#endif 

     /* [a] : Install elment attribute name into Symbol Table */

     hp = install(name);
     if(hp == NULL)
        FatalError("In Add_Fiber_Elmt_Attr() : hp = NULL !!!\n",(char *)NULL);
     else
        hp->u.fep = fep;

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

#ifdef __STDC__
int Add_Section_Attr(char *name, SECTION_ATTR *sp)
#else
int Add_Section_Attr(name, sp)
char          *name;
SECTION_ATTR    *sp;
#endif
{
SYMBOL          *hp;
int               i;

#ifdef DEBUG
       printf("*** Enter Add_Section_Attr() : name = \"%s\"\n", name);
#endif 

     /* [a] : Install Section into Symbol Table */

     hp = install(name);
     if(hp == NULL)
        FatalError("In Add_Section_Attr() : hp = NULL !!!\n",(char *)NULL);
     if(hp->u.sap == NULL)  hp->u.sap = sp;

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

}

#ifdef __STDC__
int Add_Material_Attr(char *name, MATERIAL_ATTR *mp)
#else
int Add_Material_Attr(name, mp)
char           *name;
MATERIAL_ATTR    *mp;
#endif
{
SYMBOL              *hp;
int    iInPlaneIntegPts;
int  iThicknessIntegPts;
int       iNO_INTEG_pts;
int                i, j;

#ifdef DEBUG
       printf("*** Enter Add_Material_Attr()\n");
#endif 

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

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

      iNO_INTEG_pts = iInPlaneIntegPts*iThicknessIntegPts;


     /* [0] calculate non-independent parameters */

     switch(((int) mp->LC_ptr->beta)){
        case 0:  /* Kinematic hardening */
        case 1:  /* Isotropic hardening */
          for(j = 1; j <= iNO_INTEG_pts; j++)
              for(i = 1; i <= 6; i++)
                  mp->LC_ptr->back_stress[i-1][j-1] = 0.0;
        break;
        default:
           printf(" In Add_Material_Attr(): UNDEFINED STRAIN_HARDENING \n");
           FatalError("***** Undefined strain hardening parameter found \n", (char *) NULL);
        break;
     }
     /* For inital yielding */
     for(j = 1; j <= iNO_INTEG_pts; j++)
         mp->LC_ptr->R[j-1] = sqrt(2.0/3.0)*mp->fy.value;
     
     if( mp->LC_ptr->name != (char *)NULL &&
         !strcmp(mp->LC_ptr->name, "Bi-Linear")) {
        for(j = 1; j <= iNO_INTEG_pts; j++) {
            mp->LC_ptr->H[j-1] = mp->ET.value/(1.0 - mp->ET.value/mp->E.value); 
        }
     }
     if(mp->LC_ptr->name != (char *)NULL &&