fe_profile.c

有限元分析源代码

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *         fe_profile.c : Major fuctions for FE Profiler
 *                                                                     
 *  Copyright (C) 1995 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 for any 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 notice is to remain intact.
 *                                                                    
 *  Written by: Mark Austin, Xiaoguang Chen, and Wane-Jang Lin      December 1995
 *  ============================================================================= 
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

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

/*
#define DEBUG 
*/

static int max_size_of_stiff;
static int max_no_dof;
static int max_nodes_per_elmt;
static int max_no_dimen;

/* ================================== */
/* profile() :  Input: EFRAME *frp    */
/* ================================== */

#ifdef __STDC__
EFRAME *profile(EFRAME *frame, int prt, int iop)
#else
EFRAME *profile(frame, prt, iop)
EFRAME *frame;
int    prt, iop;
#endif
{
int n, i,j,ii,jj,m,neq,nad;
int idl_vector[30];
int nmin;

#ifdef DEBUG
       printf("*** Enter profile() : iop        = %4d \n", iop);
       printf("                    : frame->no_dof = %4d \n", frame->no_dof);
       printf("                    : frame->no_nodes = %4d \n", frame->no_nodes);
#endif

   prt = PRINT_PROFILE; /* PRINT_PROFILE = YES has been set in */
                        /* fe_setflags.c prt = PRINT_PROFILE;  */
                        /* option allows profile to be printed */
                        /* all the time                        */

   switch(iop){
       case EQNOS:
            nmin = 0; nad = 0;

            for(i = 1; i <= frame->no_nodes; i++) {
                for(j = 1; j <= frame->no_dof; j++) {
	            if(frame->node[i-1].bound_id[j-1] < nmin)
	               nmin = frame->node[i-1].bound_id[j-1];
             }
            }
         
            neq = -nmin;
            for(n = 1 ; n <= frame->no_nodes; n++) {
                for(i = 1; i <= frame->no_dof; i++)  {
                    j =  frame->node[n-1].bound_id[i-1];
                    if(j < 0) {    
                       frame->node[n-1].bound_id[i-1] 
                       = -frame->node[n-1].bound_id[i-1]; 
                    }             
                    else if(j == 0) {
                       neq = neq + 1;
                       frame->node[n-1].bound_id[i-1] = neq; 

                       /* unconstrained nodes have positive bound_id */

                       CheckJdiagSpace(frame, neq);
                       frame->jdiag[neq-1]      = 0; 
                    }
                    else {
                       nad = nad - 1;
                       frame->node[n-1].bound_id[i-1]  = nad; 
                        /* change restrained node ID from 1 to -1 */
                    }
                }
            }
            frame->no_eq = neq;

            /* compute diagonal pointers for profile */

            nad = 1;  
            frame->jdiag[0] = 1;
            if(frame->no_eq == 1)
               return(frame);

            CheckJdiagSpace(frame, frame->no_eq);
            for(n = 2; n <= frame->no_eq; n++)
                frame->jdiag[n-1] = frame->jdiag[n-1] + frame->jdiag[n-2] +1;

            nad = frame->jdiag[frame->no_eq-1];
  /**********    frame->nad = nad;           *******/
            break;
        default:
            break;
   } 

   /* 
    *  ===============
    *  Print profile()
    *  ===============
    */ 

   if(prt != YES)
      return (frame);

   printf("\n");
   printf("===================================\n");
   printf("           PROFILE                 \n");  
   printf("===================================\n");
   printf("\n");

   printf("node");
   for(i=1; i <= frame->no_dof; i++) {
       switch(i) {
           case 1:
                printf("    dof1");
                break;
           case 2:
                printf("    dof2");
                break;
           case 3:
                printf("    dof3");
                break;
           case 4:
                printf("    dof4");
                break;
           case 5:
                printf("    dof5");
                break;
           case 6:
                printf("    dof6");
                break;
           default:
                break;
         }
   }

   printf("\n");
   printf("----------------------------------------------------\n");

   for(i = 1 ; i <= frame->no_nodes; i++) {
       printf("\n%4d ",i);
       for(j = 1; j <= frame->no_dof; j++)
           printf("\t%4d", frame->node[i-1].bound_id[j-1]);
   }
   printf("\n");
   printf("----------------------------------------------------\n\n");

   return(frame);
}


/* ===========================================================*/
/*  rplink                                                    */
/*  Input: EFRAME *frp                                        */
/*  objective: This subroutine is used to link rigid body and */
/*  flexible columns connected by nodes                       */ 
/* ===========================================================*/

#ifdef __STDC__
EFRAME *rplink(EFRAME *frp, MATRIX *dp, int *idl)
#else
EFRAME *rplink(frp,dp,idl)
EFRAME *frp;
MATRIX *dp;
int    *idl;
#endif
{
MATRIX *sp;
int j,ii,i,nmax,ntemp,nelim, ip, dofl;
int prt = PRINT_PLINK;
int n1;

#ifdef DEBUG
       printf("Enter rplink() : idl[]\n");
       printf("               : ");
       for(i = 1; i<= frp->no_dof;i++)
           printf("%4d ", idl[i-1]);
       printf("\n");
#endif

    /* Find max -ve index in Nodal bound_id[] array */

    nmax = 0;
    for(i = 1; i <= frp->no_nodes; i++) {
        for(j = 1; j <= frp->no_dof; j++)
            if(frp->node[i-1].bound_id[j] < nmax)
               nmax = frp->node[i-1].bound_id[j];
    }

    return(frp);
}


/* =================================================== */
/*  plink                                              */
/*  Input: EFRAME *frp                                 */
/* =================================================== */

#ifdef __STDC__
EFRAME *plink(EFRAME *frp, MATRIX *dofp, int *idl)
#else
EFRAME *plink(frp, dofp, idl)
EFRAME *frp;
MATRIX *dofp;
int    *idl;
#endif
{
MATRIX *sp;
int j,ii,i,nmax,ntemp,nelim, ip, dofl;
int prt = PRINT_PLINK;
int n1;

#ifdef DEBUG
       printf("Enter plink() : idl[]\n");
       printf("              : ");
       for(i = 1; i <= frp->no_dof; i++)
           printf("%4d ", idl[i-1]);
       printf("\n");
#endif

    /* Find max -ve index in Nodal id[fcp->no_dof] array */

       nmax = 0;
       for(i = 1; i <= frp->no_nodes; i++) {
           for(j = 1; j <= frp->no_dof; j++)
	       if(frp->node[i-1].bound_id[j] < nmax)
	          nmax = frp->node[i-1].bound_id[j];
       }

    return(frp);
}


/* ===================================================== */
/*  pload  ;form load vector in compact form             */
/*  Input: EFRAME *frp,                                  */
/*         load vector f(1,TDOF)                         */
/*         multiplier factor p                           */
/*         nneq = no_nodes * ndf = TDOF                  */
/* ===================================================== */

#ifdef __STDC__
QUANTITY *pload(EFRAME *frp, QUANTITY *f, QUANTITY *b, int nneq, QUANTITY p)
#else
QUANTITY *pload(frp,f,b,nneq,p)
EFRAME   *frp;
QUANTITY *f, *b, p;
int      nneq;
#endif
{
int n, i,j,jj;

    for(n = 1; n <= nneq; n++)
        b[n].value = 0.0;
        for(n =1; n <=frp->no_nodes; n++){
        for(i =1; i <=frp->no_dof; i++){
            j = frp->node[n-1].bound_id[i-1];
            if(j >0){
               jj = (n-1) * frp->no_dof+ i;
               b[j].value = f[jj].value * p.value + b[j].value;
            }
        }
        }

    return(b);

}


/* 
 *  ----------------------------------------------------------------- 
 *  Function : ARRAY  *Assign_p_Array()                           
 *                                                                 
 *  Transfers info from EFRAME *frp to working array ARRAY *p
 *  for a particular element elmt_no.                          
 *  Depending on task, elmt arrays in ARRAY *p              
 *  are checked to suit for element arrays &                    
 *  
 *  Input:    EFRAME *frame                                     
 *            int    elmt_no                                 
 *            ARRAY  *p                                        
 *            int        task                                     
 *  
 *  Output:   ARRAY  *p                                         
 *  ----------------------------------------------------------------- 
 */ 

#ifdef __STDC__
ARRAY *Assign_p_Array(EFRAME *frp, int elmt_no, ARRAY *p, int task)
#else
ARRAY *Assign_p_Array(frp, elmt_no, p, task)
EFRAME *frp;
ARRAY  *p;
int    elmt_no, task;
#endif
{
DIMENSIONS  *dp_length_temp, *dp_force_temp;
NODE                                    *np;
ELEMENT                                 *el;
ELEMENT_ATTR                           *eap;
MATRIX                                   *m;
int      ma,nen,j,i,node_no, jj,num_node, k;
int                         n, elmt_attr_no;
int                            dof_per_elmt;
int                                  length;
int                            UNITS_SWITCH;

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

#ifdef DEBUG