code.c

有限元分析源代码

        }
     }
     free((char *)d1.q);
     free((char *)d2.q);
     MatrixFree( d3.m );
     Push(d);
}

/*
 *  ======================================================
 *  Evaluation of Variables and Strings
 *  
 *  Purpose : Transfer contents from Symbol Table to Stack
 *  ======================================================
 */

int String_Eval()
{
DATUM  d1;

    d1 = Pop();
    Push(d1);
}

int Variable_Eval()
{
DATUM  d1, d2;
int         i;
int   length1;

   /* variable name, d1.sym, stored in symbol table, cannot be freed */
    d1 = Pop();

    if(d1.sym->type != VAR && d1.sym->type != QUAN)
       FatalError("Attempt to evaluate non-variable\n", (char *)NULL);

    d2.q = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));
    if( d1.sym->u.q != (QUANTITY *)NULL ) {
        switch( CheckUnits() ) {
          case ON:
             d2.q->value  = d1.sym->u.q->value;
             d2.q->dimen  = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
             UnitsCopy( d2.q->dimen, d1.sym->u.q->dimen );
             break;
          case OFF:
             d2.q->value  = d1.sym->u.q->value;
             d2.q->dimen  = (DIMENSIONS *)NULL;
             break;
          default:          
             FatalError("In Variable_Eval() : CheckUnits is not ON or OFF",(char *)NULL);
             break;
        }
    }
    else
      if( CheckUnits()==ON )
          d2.q->dimen  = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));

    Push(d2);
}

int Matrix_Eval()
{
DATUM  d1, d2;
int iNoRows;
int iNoColumns;
int i, j;

   /* matrix, stored in symbol table, cannot be free */
    d1 = Pop();
    if(d1.sym->type != MATX) {
       printf("FATAL ERROR >> Execution halted in Matrix_Eval()\n");
       printf("FATAL ERROR >> Problem : %s->type != MATX\n", d1.sym->cpSymName);
       printf("FATAL ERROR >> Recommendation : Check input file !! \n");
       FatalError("Assignment to non-matrix",(char *)NULL);
    }
    d2.m = MatrixCopy(d1.sym->u.m);
    Push(d2);
}

int Matrix_Build()
{
DATUM    d1, d2, d3;
int iNoRows, row_length;
int             i,j,r,s;
int              length;
int        UNITS_SWITCH;

    d1 = Pop();  /* constant= no of rows = address of d1.sym */
    iNoRows = (int) d1.sym;

    d2.m = (MATRIX *) MyCalloc(1,sizeof(MATRIX));
    d2.m->cpMatrixName    = (char *) NULL;         /* matrix name */
    d2.m->eType    = DOUBLE_ARRAY;
    d2.m->eRep     = INDIRECT;
    d2.m->iNoRows  = iNoRows;
    d2.m->uMatrix.daa = (double **) MyCalloc(iNoRows, sizeof(double *));

    UNITS_SWITCH = CheckUnits();
    for(i = 1; i <= iNoRows; i++) {
        d3  = Pop();    /* constant= no of columns = address of d3.sym */
        row_length = (int) d3.sym;
        r = iNoRows - i;
        d2.m->uMatrix.daa[r]    = (double *) MyCalloc(row_length, sizeof(double));

        if(i == 1) {
          d2.m->iNoColumns    = row_length;
          if( UNITS_SWITCH == ON ) {
              d2.m->spColUnits = (DIMENSIONS *) MyCalloc(row_length, sizeof(DIMENSIONS));
              d2.m->spRowUnits = (DIMENSIONS *) MyCalloc(iNoRows, sizeof(DIMENSIONS));
          }
          else {
              d2.m->spColUnits = (DIMENSIONS *)NULL;
              d2.m->spRowUnits = (DIMENSIONS *)NULL;
          }
        }
        
        /* check for no of columns in matrix */

        if(i != 1 && d2.m->iNoColumns != row_length ) {
           printf(" In Matrix_Build(): \n");
           FatalError("==> no of columns in each row of a matrix must be same",(char *)NULL);
        }

        for(j = 1; j <= row_length; j++) {
            d3  = Pop();    /* matrix items, d3.q */
            s = row_length - j;
            d2.m->uMatrix.daa[r][s] = (double) d3.q->value;

        /* Units of matrix elements are stored in column_units_buf as an default */

           switch( UNITS_SWITCH ) {
             case ON:
                if( (row_length != 1) || ((iNoRows == 1) && (row_length == 1))  ) {
                   UnitsCopy( &(d2.m->spColUnits[s]), d3.q->dimen );
                }
                else {
                   ZeroUnits( &(d2.m->spColUnits[s]) );
                   d2.m->spColUnits[s].units_type = SI_US;
                }

                if((row_length == 1) && (iNoRows > 1) ) { /* a vector */
                   break;
                }
                else {
                   free((char *) d3.q->dimen->units_name);
                   free((char *) d3.q->dimen);
                   free((char *) d3.q);
                }
                break;
             case OFF:
                if((row_length == 1) && (iNoRows > 1) )   /* a vector */
                   break;
                else
                   free((char *) d3.q);
                break;
             default:          
                FatalError("In Matrix_Build() : CheckUnits is not ON or OFF",(char *)NULL);
                break;
           }
        }
        if( UNITS_SWITCH == ON ) {
            if((row_length == 1) && (iNoRows > 1) ) { /* a vector */
               UnitsCopy( &(d2.m->spRowUnits[r]), d3.q->dimen );
               free((char *) d3.q->dimen->units_name);
               free((char *) d3.q->dimen);
               free((char *) d3.q);
            } 
            else {
               ZeroUnits( &(d2.m->spRowUnits[r]) );
               d2.m->spRowUnits[r].units_type = SI_US;
            }
        }
   }
   Push(d2);
}

int Assign_Matrix()
{
DATUM        d1, d2, d3;
int iNoRows, iNoColumns;
int                i, j;
int              length;
int                 *ld;
int          iColHeight;
int        UNITS_SWITCH;

    d1     = Pop();  /* matrix name , d1.sym */
    d2     = Pop();  /* matrix contents , d2.m */

    iNoRows    = d2.m->iNoRows;
    iNoColumns = d2.m->iNoColumns;

    /* Check for assignment to "variable" or "matrix" */

    if(d1.sym->type != VAR && d1.sym->type != MATX)
       EXECUTION_ERROR("Assignment to non-variable\n", d1.sym->cpSymName);
 
    /* Check for change in "dimensions" and "type" of matrix */

    if(d1.sym->type == MATX && d1.sym->u.m != (MATRIX *) NULL)  {
       if(d1.sym->u.m->iNoRows    != iNoRows     ||
          d1.sym->u.m->iNoColumns != iNoColumns  ||
          d1.sym->u.m->eType      != d2.m->eType ||
          d1.sym->u.m->eRep       != d2.m->eRep  ) {
             MatrixFree(d1.sym->u.m);
             d1.sym->u.m = (MATRIX *) NULL;
          }
    }

    /* Allocate new matrix if needed */

    if( d1.sym->u.m == (MATRIX *) NULL || d1.sym->type != MATX) {

       d1.sym->type = MATX;
       switch((int) d2.m->eRep ) {
           case INDIRECT:
           {
               switch((int) d2.m->eType) {
                   case DOUBLE_ARRAY:
                        d1.sym->u.m 
                        = MatrixAllocIndirect((char *) d1.sym->cpSymName,
                          DOUBLE_ARRAY, iNoRows, iNoColumns);
                        break;
                   default:
                        FatalError("In Assign_Matrix() : Undefined Matrix Type",
                                  (char *)NULL);
                        break;
               }
               break;
           }
           case SKYLINE:
           {
               ld = (int *)iVectorAlloc( d2.m->iNoRows );
               for( i=0 ; i<d2.m->iNoRows ; i++ )
                   ld[i] = (int) d2.m->uMatrix.daa[i][0];
               switch((int) d2.m->eType) {
                   case DOUBLE_ARRAY:
                        d1.sym->u.m 
                        = MatrixAllocSkyline((char *) d1.sym->cpSymName,
                          DOUBLE_ARRAY, iNoRows, iNoColumns, ld);
                        free((char *) ld);
                        break;
                   default:
                        FatalError("In Assign_Matrix() : Undefined Matrix Type",
                        (char *)NULL);
                        break;
               }
               break;
           }
           default:
               FatalError("In Assign_Matrix() : Undefined Matrix eRep",
               (char *)NULL);
               break;
       } /* end of switch( d2.m->eRep ) */

    } /* end of if */

    /* Copy contents of Matrix d2.m to d1.sym->u.m with no reallocation of memory */

    UNITS_SWITCH = CheckUnits();
    switch((int) d2.m->eRep ) {
      case INDIRECT:
      {
         switch((int) d2.m->eType) {
           case DOUBLE_ARRAY:
           {
             for(i = 1; i <= iNoRows; i++) {
                if( UNITS_SWITCH == ON )
                    UnitsCopy( &(d1.sym->u.m->spRowUnits[i-1]), &(d2.m->spRowUnits[i-1]) ); 

                for(j = 1; j <= iNoColumns; j++)
                   d1.sym->u.m->uMatrix.daa[i-1][j-1]
                   = d2.m->uMatrix.daa[i-1][j-1];
             }
             if( UNITS_SWITCH == ON ) {
                 for(j = 1; j <= iNoColumns; j++)
                    UnitsCopy( &(d1.sym->u.m->spColUnits[j-1]), &(d2.m->spColUnits[j-1]) ); 
             }
             break;
           }
           default:
                FatalError("In Assign_Matrix() : Undefined Matrix Type",
               (char *)NULL);
                break;
         } /* end of case INDIRECT, switch( eType ) */
         break;
      } /* end of case INDIRECT */
      case SKYLINE:
      {
         switch((int) d2.m->eType) {
           case DOUBLE_ARRAY:
           {
             for(i = 1; i <= iNoRows; i++) {
                if( UNITS_SWITCH == ON )
                    UnitsCopy(&(d1.sym->u.m->spRowUnits[i-1]),&(d2.m->spRowUnits[i-1])); 

                iColHeight = (int) MAX(d1.sym->u.m->uMatrix.daa[i-1][0],
                                       d2.m->uMatrix.daa[i-1][0]);
                if( iColHeight > d1.sym->u.m->uMatrix.daa[i-1][0] ) {
                   free( d1.sym->u.m->uMatrix.daa[i-1] );
                   d1.sym->u.m->uMatrix.daa[i-1] 
                   = (double *)MyCalloc((iColHeight+1), sizeof(double));
                   d1.sym->u.m->uMatrix.daa[i-1][0] = iColHeight;
                }
                for(j = 1; j <= d2.m->uMatrix.daa[i-1][0]; j++) 
                   d1.sym->u.m->uMatrix.daa[i-1][j]  = d2.m->uMatrix.daa[i-1][j];
             }
             if( UNITS_SWITCH == ON ) {
                 for(j = 1; j <= iNoColumns; j++)
                    UnitsCopy(&(d1.sym->u.m->spColUnits[j-1]), &(d2.m->spColUnits[j-1])); 
             }
             break;
           } /* end of case SKYLINE, DOUBLE_ARRAY */
           default:
                FatalError("In Assign_Matrix() : Undefined Matrix Type",
                          (char *)NULL);
                break;
         } /* end of case SKYLINE , switch( eType ) */
         break;
      } /* end of case SKYLINE */
      default:
           FatalError("In Assign_Matrix() : Undefined Matrix eRep",(char *)NULL);
           break;
    } /* end of switch(d2.m->eRep) */

    MatrixFree(d2.m);
    Push(d1);
}

int Assign_Matrix_Item()
{
DATUM  d1, d2, d3, d4, d5;
int row, column;
int i, j, length;
int iMin, iMax;
int UNITS_SWITCH;

DIMENSIONS *dm1;
DIMENSIONS *dm_ptr;

    d1 = Pop();         /* matrix name , d1.sym    */
    d2 = Pop();         /* matrix item , d2.q      */
    d3 = Pop();         /* matrix iNoColumns, d3.q */
    d4 = Pop();         /* matrix iNoRows   , d4.q */

    UNITS_SWITCH = CheckUnits();

    row    = (int) d4.q->value;
    column = (int) d3.q->value;
    if( UNITS_SWITCH == ON ) {
      free((char *) d3.q->dimen->units_name);
      free((char *) d3.q->dimen);
      free((char *) d4.q->dimen->units_name);
      free((char *) d4.q->dimen);
    }
    free((char *) d3.q);
    free((char *) d4.q);

    /* Check Variable type and Matrix dimensions */

    if(d1.sym->type != MATX)
       FatalError("In Assign_Matrix_Item : Assignment to non-matrix",(char *)NULL);

    if((row < 1) || (row > d1.sym->u.m->iNoRows))
       FatalError("In Assign_Matrix_Item : row out of array bounds",(char *)NULL);
     
    if((column < 1) || (column > d1.sym->u.m->iNoColumns))
       FatalError("In Assign_Matrix_Item : column out of array bounds",(char *)NULL);

    /* Update Matrix Item */

    switch( d1.sym->u.m->eRep ) {
       case INDIRECT :
          switch((int) d1.sym->u.m->eType) {
              case DOUBLE_ARRAY:

                switch( UNITS_SWITCH ) {
                  case  ON:
                  /* [a]  calculate the units of m[row][col] */


                     dm1 = UnitsMult( &(d1.sym->u.m->spRowUnits[row-1]),
                                      &(d1.sym->u.m->spColUnits[column-1]) );

                     if(dm1->units_name  != (char *)NULL  ||
                        dm1->length_expnt != 0     ||
                        dm1->mass_expnt != 0       ||
                        dm1->time_expnt != 0       ||
                        dm1->temp_expnt != 0 ) {
                     
                        /* check units of m[row][col]    */
                     
                        if(SameUnits(dm1, d2.q->dimen) != TRUE) {
                           printf("\n **** WARNING : in  Assign_Matrix_Item(): \n");
                           printf(" =====> %s[%d][%d] and %15.4e %s have inconsistent units \n\n",
                                  d1.sym->u.m->cpMatrixName, row, column, d2.q->value/d2.q->dimen->scale_factor,
                                  d2.q->dimen->units_name);
                           printf(" =====> This assignment will NOT overwrite the previous units\n");
                        }
                        else
                          /* Assign the value to matrix element */
                          d1.sym->u.m->uMatrix.daa[row-1][column-1] = d2.q->value;
                     }
                     else {

                       /* [x] m[row][column] has no unit                */
                       /* Save units in the col_units_buf as an default */

                       if(d2.q->dimen->units_name != NULL    ||
                          d2.q->dimen->length_expnt != 0     ||
                          d2.q->dimen->mass_expnt != 0       ||
                          d2.q->dimen->time_expnt != 0       ||
                          d2.q->dimen->temp_expnt != 0 ) {
/*
                          printf("\n ****** WARNING: Matrix %s element %s[%d][%d] \n",
                                 d1.sym->u.m->cpMatrixName,d1.sym->u.m->cpMatrixName, row, column);
                          printf("****WARNING: assign a units to a single element might affect the units of other elements\n\n");
*/
                          d1.sym->u.m->uMatrix.daa[row-1][column-1] = d2.q->value;
                          UnitsCopy( &(d1.sym->u.m->spColUnits[column-1]), d2.q->dimen );
                       }
                       else /* both matrix and d2.q are non-dimensional */
                          d1.sym->u.m->uMatrix.daa[row - 1][column - 1] = d2.q->value;
                     }

                     free((char *) d2.q->dimen->units_name);
                     free((char *) d2.q->dimen);
                     free((char *) dm1->units_name);
                     free((char *) dm1);

                     break;   /* end of case UNITS_SWITCH==ON */
                  case OFF:
                     d1.sym->u.m->uMatrix.daa[row-1][column-1] = d2.q->value;
                     break;
                  default:          
                     FatalError("In Assign_Matrix_Item() : CheckUnits is not ON or OFF",(char *)NULL);
                     break;
                }
                break;  /* end of case INDIRECT, DOUBLE_ARRAY */
              default:
                FatalError("In Assign_Matrix_Item() : Undefined m->eType",(char *)NULL);
                break;
          }
          break;  /* end of case INDIRECT */
       case SKYLINE:
          switch((int) d1.sym->u.m->eType) {
              case DOUBLE_ARRAY:

                switch( UNITS_SWITCH ) {
                   case ON:
                     /* [a]  calculate the units of m[row][col] */

                     dm1 = UnitsMult( &(d1.sym->u.m->spRowUnits[row-1]),
                                      &(d1.sym->u.m->spColUnits[column-1]) );

                     if(dm1->units_name  != NULL  ||
                        dm1->length_expnt != 0     ||
                        dm1->mass_expnt != 0       ||
                        dm1->time_expnt != 0       ||
                        dm1->temp_expnt != 0 ) {
                     
                        /* check units of m[row][col]    */
                     
                        if(SameUnits(dm1, d2.q->dimen) != TRUE) {