code.c

有限元程序

                     }
                     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) {
                           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 */
                           iMin = MIN(row, column);
                           iMax = MAX(row, column);
                           if((iMax-iMin+1) <= d1.sym->u.m->uMatrix.daa[iMax-1][0])
                               d1.sym->u.m->uMatrix.daa[iMax-1][iMax-iMin+1] = d2.q->value;
                           else {
                               d1.sym->u.m->uMatrix.daa[iMax-1]
                               = (double *) realloc(d1.sym->u.m->uMatrix.daa[iMax-1],
                                 (iMax-iMin+1));
                               for( j=(int) d1.sym->u.m->uMatrix.daa[iMax-1][0]+1; j < iMax-iMin+1 ; j++ )
                                   d1.sym->u.m->uMatrix.daa[iMax-1][j] = 0.0;
                               d1.sym->u.m->uMatrix.daa[iMax-1][0] = iMax-iMin+1;
                               d1.sym->u.m->uMatrix.daa[iMax-1][iMax-iMin+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:
                     iMin = MIN(row, column);
                     iMax = MAX(row, column);
                     if((iMax-iMin+1) <= d1.sym->u.m->uMatrix.daa[iMax-1][0])
                         d1.sym->u.m->uMatrix.daa[iMax-1][iMax-iMin+1] = d2.q->value;
                     else {
                         d1.sym->u.m->uMatrix.daa[iMax-1] 
                         = (double *) realloc(d1.sym->u.m->uMatrix.daa[iMax-1], (iMax-iMin+1));
                         for( j=(int) d1.sym->u.m->uMatrix.daa[iMax-1][0]+1; j < iMax-iMin+1 ; j++ )
                             d1.sym->u.m->uMatrix.daa[iMax-1][j] = 0.0;
                         d1.sym->u.m->uMatrix.daa[iMax-1][0] = iMax-iMin+1;
                         d1.sym->u.m->uMatrix.daa[iMax-1][iMax-iMin+1] = d2.q->value;
                     }
                     break;  /* end of case UNITS_SWITCH == OFF */
                  default:
                     FatalError("In Matrix_Assign_Item(): CheckUnits() is not ON or OFF",
                     (char *)NULL);
                     break;
                  }
                  break;  /* end of case SKYLINE, DOUBLE_ARRAY */
              default:
                  FatalError("In Assign_Matrix_Item() : Undefined m->eType",(char *)NULL);
                  break;
          }
          break;  /* end of case SKYLINE */
       default:
          FatalError("In Assign_Matrix_Item() : Undefined m->eRep",(char *)NULL);
          break;
    }  /* end of switch of m->eRep */

    free((char *) d2.q);

    Push(d1);
}

int Assign_Quantity()
{
DATUM     d1, d2;
int       length;
int UNITS_SWITCH;

    d1 = Pop();   /* variable name , d1.sym */
    d2 = Pop();   /* value, scale factor and dimensions, d2.q */

    if(d1.sym->type != VAR && d1.sym->type != QUAN) {
       FatalError("In Assign_Quantity() : Assignment to non-variable",(char *)NULL);
    }

    UNITS_SWITCH = CheckUnits();
    if(d1.sym->type != QUAN) {
       d1.sym->u.q   = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));
       d1.sym->type  = QUAN;
       if( UNITS_SWITCH == ON )
          d1.sym->u.q->dimen = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
    }

    switch( UNITS_SWITCH ) {
      case ON:
         d1.sym->u.q->value = d2.q->value;
         UnitsCopy( d1.sym->u.q->dimen,  d2.q->dimen );
         free((char *) d2.q->dimen->units_name);
         free((char *) d2.q->dimen);
         break;
      case OFF:
         d1.sym->u.q->value             = d2.q->value;
         break;
      default:          
         FatalError("In Assign_Quantity() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d2.q);
    Push(d1);
}


/*
 *  ------------------------------------------
 *  Functions to Print Quantities and Matrices
 *  ------------------------------------------
 */

int Print_Expr()
{
DATUM d;
int UnitsType;

    d = Pop();

    switch(CheckUnits()) {
      case ON:
         UnitsSimplify(d.q->dimen);
         if(d.q->dimen->units_name != NULL) {
           UnitsType = CheckUnitsType(); 
           switch(UnitsType) { 
             case SI:
                  if(!strcmp(d.q->dimen->units_name, "deg_F") )
                     d.q->value = ConvertTempUnits(d.q->dimen->units_name, d.q->value, US); 
                  break;
             case US:
                  if(!strcmp(d.q->dimen->units_name, "deg_C") ) 
                     d.q->value = ConvertTempUnits(d.q->dimen->units_name, d.q->value, SI); 
                  break;
           }
           printf("%10.4g ", d.q->value/d.q->dimen->scale_factor);
           printf("%s ", d.q->dimen->units_name);
           free((char *) d.q->dimen->units_name);
         } else {
           printf("%10.4g ", d.q->value);
         }
         fflush(stdout);
         free((char *) d.q->dimen);
         break;
      case OFF:
         printf("%10.4g ", d.q->value);
         break;
      default:          
         FatalError("In Print_Expr() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d.q);
}

/*
 *  =================================================
 *  Print_Dimen_Expr() : Print Dimensional Expression
 *  =================================================
 */

int Print_Dimen_Expr() {
DATUM d1, d2;
int UnitsType;

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

    d2 = Pop(); /* units    */
    d1 = Pop(); /* quantity */

#ifdef DEBUG
       printf("*** d1.q->value = %10.5f\n", d1.q->value );
#endif

    switch(CheckUnits()) {
      case ON:
           if(SameUnits(d1.q->dimen, d2.q->dimen) == TRUE) {
           if(d2.q->dimen->units_name != NULL) {
              UnitsType = CheckUnitsType(); 
              switch(UnitsType) {
                 case SI:
                      if(!strcmp(d2.q->dimen->units_name, "deg_F")) 
                         d1.q->value = ConvertTempUnits(d1.q->dimen->units_name, d1.q->value, US);
                      break;
                 case US:
                      if(!strcmp(d2.q->dimen->units_name, "deg_C")) 
                         d1.q->value = ConvertTempUnits(d1.q->dimen->units_name, d1.q->value, SI);
                      break;
              }
              printf("%10.4g ", d1.q->value/d2.q->dimen->scale_factor);
              printf("%s ", d2.q->dimen->units_name);
              free((char *) d2.q->dimen->units_name);
           } else {
              printf("%10.4g ", d1.q->value);
           }
           fflush(stdout);
           free((char *) d1.q->dimen->units_name);
           free((char *) d1.q->dimen);
           free((char *) d2.q->dimen);
           } else {
             printf(" Wanted units_name   = %s\n", d2.q->dimen->units_name);
             printf(" quantity units_name = %s\n", d1.q->dimen->units_name);
             FatalError(" In Print_Dimen_Expr(): ",
                        "Try to print quantity with inconsistent units",
                        (char *) NULL);
           }
           break;
      case OFF:
           printf("%10.4g ", d1.q->value);
           break;
      default:          
           FatalError("In Print_Dimen_Expr() : CheckUnits is not ON or OFF",(char *)NULL);
           break;
    }
    free((char *) d1.q);
    free((char *) d2.q);

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

}

int Print_String()
{
    printf("%s", (char *) *pc);
    pc = pc+1;
    fflush(stdout);
}


/* 
 *  ---------------------------------------------
 *  Call Built-in Expression and Matrix Functions 
 *  ---------------------------------------------
 */

int Bltin_Break()
{
     iBREAK_FLAG = (int) 1;
}

int Check_Break()
{
    return(iBREAK_FLAG);
}

int After_Break()
{
    iBREAK_FLAG = (int) 0;
}

int Bltin_Quantity()
{
DATUM  d;

  /* random number function, return a new quantity */
     d.q  = (*(QUANTITY * (*)()) *pc)();
     pc = pc+1;
     Push(d);
}

int Bltin1_Quantity()
{
DATUM  d;

  /* mathematic functions, return the old quantity */
     d = Pop();  /* parameter value, d.q */
     d.q  = (*(QUANTITY * (*)()) *pc)(d.q);
     pc = pc+1;
     Push(d);
}

int Bltin2_Quantity()
{
DATUM d2, d;

  /* mathematic function, Pow(quantity,value), return old quantity */
     d2 = Pop();  /* parameter 2 value, d2.q */
     d  = Pop();  /* parameter 1 value, d.q */
     d.q = (*(QUANTITY * (*)()) *pc)(d.q, (double) d2.q->value);

     if( CheckUnits() == ON ) {
       free((char *) d2.q->dimen->units_name);
       free((char *) d2.q->dimen);
     }
     free((char *) d2.q);
     pc = pc+1;
     Push(d);
}

int Bltin3_Quantity()
{
DATUM d1, d2;

     d1 = Pop();  /* parameter matrix, d1.m */
     d2.q = (*(QUANTITY * (*)()) *pc)(d1.m);
     MatrixFree(d1.m);
     pc = pc+1;
     Push(d2);
}

int Bltin_Matrix()
{
DATUM d;

     d.m = (*(MATRIX * (*)()) *pc)();
     pc = pc+1;
     Push(d);
}

int Bltin1_Matrix()
{
DATUM d1;
DATUM d2;

     d1 = Pop();  /* matrix, d.m */
     d2.m = (*(MATRIX * (*)()) *pc)(d1.m);
     MatrixFree(d1.m);
     pc = pc+1;
     Push(d2);
}

int Bltin2_Matrix()
{
DATUM d1, d2, d3;

     d2 = Pop();
     d1 = Pop();
     d3.m = (*(MATRIX * (*)()) *pc)(d1.m, d2.m);
     MatrixFree( d1.m );
     MatrixFree( d2.m );
     pc = pc+1;
     Push(d3);
}

int Bltin3_Matrix()
{
DATUM d1, d2, d3, d4;

     d3 = Pop();
     d2 = Pop();
     d1 = Pop();
     d4.m = (*(MATRIX * (*)()) *pc)(d1.m, d2.m, d3.m);
     MatrixFree( d1.m );
     MatrixFree( d2.m );
     MatrixFree( d3.m );
     pc = pc+1;
     Push(d4);
}

int BltinVar_Matrix()
{
DATUM d1, d[5], d3;
int no_args, i;

    d1 = Pop();  /* constant = address of d1.sym */
    no_args = (int) d1.sym;

    switch(no_args) {
      case 0:
        d3.m = (*(MATRIX * (*)()) *pc)((MATRIX *) NULL);
      break;
      case 1:
        d[0]  = Pop();
        d3.m = (*(MATRIX * (*)()) *pc)(d[0].m, (MATRIX *) NULL);
        MatrixFree(d[0].m);
      break;
      
      case 2:
        d[1]  = Pop();
        d[0]  = Pop();
        d3.m = (*(MATRIX * (*)()) *pc)(d[0].m, d[1].m, (MATRIX *) NULL);
        MatrixFree(d[1].m);
        MatrixFree(d[0].m);
      break;

      case 3:
        d[2]  = Pop();
        d[1]  = Pop();
        d[0]  = Pop();
        d3.m = (*(MATRIX * (*)()) *pc)(d[0].m, d[1].m, d[2].m, (MATRIX *) NULL);
        MatrixFree(d[2].m);
        MatrixFree(d[1].m);
        MatrixFree(d[0].m);
      break;

      case 4:
        d[3]  = Pop();
        d[2]  = Pop();
        d[1]  = Pop();
        d[0]  = Pop();
        d3.m = (*(MATRIX * (*)()) *pc)(d[0].m,d[1].m, d[2].m, d[3].m, (MATRIX *) NULL);
        MatrixFree(d[3].m);
        MatrixFree(d[2].m);
        MatrixFree(d[1].m);
        MatrixFree(d[0].m);
      break;