code.c

有限元分析源代码

    d1 = Pop();  /* 1st quantity, d1.q */
    d.q        = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));

    if(d2.q->value != 0)
       d.q->value = d1.q->value/d2.q->value;
    else
       FatalError("In Quantity_Div() : Attempt to divide by zero",(char *)NULL);

    switch(CheckUnits()) {
      case ON:
         d.q->dimen = UnitsDiv( d1.q->dimen, d2.q->dimen, NO );
         free((char *) d1.q->dimen->units_name);
         free((char *) d1.q->dimen);
         free((char *) d2.q->dimen->units_name);
         free((char *) d2.q->dimen);
         break;
      case OFF:
         break;
      default:          
         FatalError("In Div() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d1.q);
    free((char *) d2.q);
    Push(d);
}

int Quantity_Negate()
{
DATUM d1;

    d1 = Pop();
    d1.q->value = -d1.q->value;
    Push(d1);
}

#ifdef __STDC__
int dComp_Le( double x, double y )
#else
int dComp_Le(x, y)
double x, y;
#endif
{
static int FLAG; 

   FLAG = (x < y );
   if(dComp_Eq(x,y)) FLAG = 1;
   return (FLAG);
}

#ifdef __STDC__
int dComp_Ge( double x, double y )
#else
int dComp_Ge(x, y)
double x, y;
#endif
{
static int FLAG; 

   FLAG = (x > y );
   if(dComp_Eq(x,y)) FLAG = 1;
   return (FLAG);
}

#ifdef __STDC__
int dComp_Eq( double x, double y )
#else
int dComp_Eq(x,y)
double x,y;
#endif
{
static int FLAG; 
double ERROR;
       
   ERROR = x - y;
   if((ABS(ERROR) == 0.0) || (ABS(ERROR) != 0.0 && ABS(ERROR) < 1E-15))
       FLAG = 1;
   else
       FLAG = 0;
   return (FLAG);
}

#ifdef __STDC__
dComp_Ne( double x, double y )
#else
dComp_Ne(x,y)
double x,y;
#endif
{
static int FLAG; 
double ERROR;
       
   ERROR = x - y;
   if(ABS(ERROR) != 0.0 && ABS(ERROR) > 1E-15)
       FLAG = 1;
   else
       FLAG = 0;
   return (FLAG);
}

int Quantity_Gt()
{
DATUM d1, d2;

    d2     = Pop();  /* 2nd quantity, d2.q */
    d1     = Pop();  /* 1st quantity, d1.q */

    switch( CheckUnits() ) {
      case ON:
         if(SameUnits(d1.q->dimen, d2.q->dimen) == TRUE)
            d1.q->value = (double) (d1.q->value > d2.q->value);
         else
            FatalError("In Quantity_Gt() : Inconsistent Dimensions",(char *)NULL);
         free((char *) d2.q->dimen->units_name);
         free((char *) d2.q->dimen);
         free((char *) d1.q->dimen->units_name);
         d1.q->dimen->units_name = (char *)NULL;
         break;
      case OFF:
         d1.q->value = (double) (d1.q->value > d2.q->value);
         break;
      default:          
         FatalError("In Quantity_Gt() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d2.q);
    Push(d1);
}

int Quantity_Lt()
{
DATUM d1, d2;

    d2 = Pop();  /* 2nd quantity, d2.q */
    d1 = Pop();  /* 1st quantity, d1.q */

    switch( CheckUnits() ) {
      case ON:
         if(SameUnits(d1.q->dimen, d2.q->dimen) == TRUE) {
            d1.q->value = (double) (d1.q->value < d2.q->value);
         }
         else
            FatalError("In Quantity_Lt() : Inconsistent Dimensions",(char *)NULL);

         free((char *) d2.q->dimen->units_name);
         free((char *) d2.q->dimen);
         free((char *) d1.q->dimen->units_name);
         d1.q->dimen->units_name = (char *)NULL;
         break;
      case OFF:
         d1.q->value = (double) (d1.q->value < d2.q->value);
         break;
      default:          
         FatalError("In Quantity_Lt() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d2.q);
    Push(d1);
}

int Quantity_Ge()
{
DATUM d1, d2;

    d2 = Pop();  /* 2nd quantity, d2.q */
    d1 = Pop();  /* 1st quantity, d1.q */

    switch( CheckUnits() ) {
      case ON:
         if(SameUnits(d1.q->dimen, d2.q->dimen) == TRUE)
            d1.q->value = (double) dComp_Ge(d1.q->value, d2.q->value);
         else
            FatalError("In Quantity_Ge() : Inconsistent Dimensions",(char *)NULL);
         free((char *) d2.q->dimen->units_name);
         free((char *) d2.q->dimen);
         free((char *) d1.q->dimen->units_name);
         d1.q->dimen->units_name = (char *)NULL;
         break;
      case OFF:
         d1.q->value = (double) dComp_Ge(d1.q->value, d2.q->value);
         break;
      default:          
         FatalError("In Quantity_Ge() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d2.q);
    Push(d1);
}

int Quantity_Le()
{
DATUM d1, d2;

    d2 = Pop();  /* 2nd quantity, d2.q */
    d1 = Pop();  /* 1st quantity, d1.q */

    switch( CheckUnits() ) {
      case ON:
         if(SameUnits(d1.q->dimen, d2.q->dimen) == TRUE)
            d1.q->value = (double) dComp_Le(d1.q->value, d2.q->value);
         else
            FatalError("In Quantity_Le() : Inconsistent Dimensions",(char *)NULL);
         free((char *) d2.q->dimen->units_name);
         free((char *) d2.q->dimen);
         free((char *) d1.q->dimen->units_name);
         d1.q->dimen->units_name = (char *)NULL;
         break;
      case OFF:
         d1.q->value = (double) dComp_Le(d1.q->value, d2.q->value);
         break;
      default:          
         FatalError("In Quantity_Le() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d2.q);
    Push(d1);
}

int Quantity_Eq()
{
DATUM d1, d2;

    d2 = Pop();  /* 2nd quantity, d2.q */
    d1 = Pop();  /* 1st quantity, d1.q */

    switch( CheckUnits() ) {
      case ON:
         if(SameUnits(d1.q->dimen, d2.q->dimen) == TRUE)
            d1.q->value = (double) dComp_Eq(d1.q->value, d2.q->value);
         else
            FatalError("In Quantity_Eq() : Inconsistent Dimensions",(char *)NULL);
         free((char *) d2.q->dimen->units_name);
         free((char *) d2.q->dimen);
         free((char *) d1.q->dimen->units_name);
         d1.q->dimen->units_name = (char *)NULL;
         break;
      case OFF:
         d1.q->value = (double) dComp_Eq(d1.q->value, d2.q->value);
         break;
      default:          
         FatalError("In Quantity_Eq() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d2.q);
    Push(d1);
}

int Quantity_Ne()
{
DATUM d1, d2;

    d2 = Pop();  /* 2nd quantity, d2.q */
    d1 = Pop();  /* 1st quantity, d1.q */

    switch( CheckUnits() ) {
      case ON:
         if(SameUnits(d1.q->dimen, d2.q->dimen) == TRUE)
            d1.q->value = (double) dComp_Ne(d1.q->value, d2.q->value);
         else
            FatalError("In Quantity_Ne() : Inconsistent Dimensions",(char *)NULL);
         free((char *) d2.q->dimen->units_name);
         free((char *) d2.q->dimen);
         free((char *) d1.q->dimen->units_name);
         d1.q->dimen->units_name = (char *)NULL;
         break;
      case OFF:
         d1.q->value = (double) dComp_Ne(d1.q->value, d2.q->value);
         break;
      default:          
         FatalError("In Quantity_Ne() : CheckUnits is not ON or OFF",(char *)NULL);
         break;
    }
    free((char *) d2.q);
    Push(d1);
}

int Quantity_And()
{
DATUM d1, d2;

    d2 = Pop();  /* 2nd quantity, d2.q */
    d1 = Pop();  /* 1st quantity, d1.q */

    d1.q->value = (double) (d1.q->value != 0 && d2.q->value != 0);
    if(CheckUnits()==ON) {
       free((char *) d2.q->dimen->units_name);
       free((char *) d2.q->dimen);
       free((char *) d1.q->dimen->units_name);
       d1.q->dimen->units_name = (char *)NULL;
    }
    free((char *) d2.q);
    Push(d1);
}

int Quantity_Or()
{
DATUM d1, d2;

    d2 = Pop();  /* 2nd quantity, d2.q */
    d1 = Pop();  /* 1st quantity, d1.q */

    d1.q->value = (double) (d1.q->value != 0 || d2.q->value != 0);
    if(CheckUnits()==ON) {
       free((char *)d2.q->dimen->units_name);
       free((char *)d2.q->dimen);
       free((char *)d1.q->dimen->units_name);
       d1.q->dimen->units_name = (char *)NULL;
    }
    free((char *) d2.q);
    Push(d1);
}

int Quantity_Yes()
{
DATUM d1;

     d1.q = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));
     if(CheckUnits() == ON){
        d1.q->dimen = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
        d1.q->dimen->units_name = (char *)NULL;
     }
     else
        d1.q->dimen = (DIMENSIONS *)NULL;
     d1.q->value = (double) (TRUE);
     Push(d1);
}

int Quantity_Not()
{
DATUM d1;

     d1 = Pop();
     d1.q->value = (double) (d1.q->value == 0);

     if( CheckUnits() == ON ) {
        free((char *) d1.q->dimen->units_name);
        d1.q->dimen->units_name = (char *)NULL;
     }
     Push(d1);
}

int Quantity_Power()
{
DATUM     d1, d2, d;
int          length;

  /* [a] : Pop quantities from stack and allocate memory for result */

     d2  = Pop();  /* 2nd quantity, d2.q */
     d1  = Pop();  /* 1st quantity, d1.q */
     d.q = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));

  /* [a] : Pop quantities from stack and allocate memory for result */

     switch( CheckUnits() ) {
       case ON:
            if(((d1.q->dimen->length_expnt != 0) ||
                (d1.q->dimen->mass_expnt   != 0) ||
                (d1.q->dimen->time_expnt   != 0) ||
                (d1.q->dimen->temp_expnt   != 0)) &&
               ((d2.q->dimen->length_expnt == 0) &&
                (d2.q->dimen->mass_expnt   == 0) &&
                (d2.q->dimen->time_expnt   == 0) &&
                (d2.q->dimen->temp_expnt   == 0))) { 
                 d.q->value = pow(d1.q->value, d2.q->value);
                 d.q->dimen = UnitsPower( d1.q->dimen, d2.q->value, NO );
            } else if (((d1.q->dimen->length_expnt == 0) &&
                (d1.q->dimen->mass_expnt   == 0) &&
                (d1.q->dimen->time_expnt   == 0) &&
                (d1.q->dimen->temp_expnt   == 0)) &&
               ((d2.q->dimen->length_expnt == 0) ||
                (d2.q->dimen->mass_expnt   == 0) ||
                (d2.q->dimen->time_expnt   == 0) ||
                (d2.q->dimen->temp_expnt   == 0))) { 
                 d.q->value = pow(d1.q->value, d2.q->value);
                 d.q->dimen = UnitsPower( d2.q->dimen, 1.0, NO );
            } else if (((d1.q->dimen->length_expnt != 0) ||
                (d1.q->dimen->mass_expnt   != 0) ||
                (d1.q->dimen->time_expnt   != 0) ||
                (d1.q->dimen->temp_expnt   != 0)) &&
               ((d2.q->dimen->length_expnt != 0) ||
                (d2.q->dimen->mass_expnt   != 0) ||
                (d2.q->dimen->time_expnt   != 0) ||
                (d2.q->dimen->temp_expnt   != 0))) { 
                 FatalError("In Quantity_Power() : Illegal Power Operation",
                           (char *) NULL);
            }
            free((char *) d2.q->dimen->units_name);
            free((char *) d2.q->dimen);
            free((char *) d1.q->dimen->units_name);
            free((char *) d1.q->dimen);
            break;
       case OFF:
            d.q->value = pow(d1.q->value, d2.q->value);
            break;
       default:          
            FatalError("In Quantity_Power() : CheckUnits is not ON or OFF",
                        (char *)NULL);
            break;
     }
     free((char *) d2.q);
     free((char *) d1.q);
     Push(d);
}

int Quantity_Affirm()
{
    Push(Pop());
}

int Quantity_Extract()
{
DATUM d1, d2, d3, d;
int i, j, length;
int iMin, iMax;

     d3 = Pop();  /* input matrix, d3.m */
     d2 = Pop();  /* index j , d2.q */
     d1 = Pop();  /* index i , d1.q */
     i   = (int) d1.q->value;
     j   = (int) d2.q->value;
     d.q = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));

     /* Check arguments are defined */

     if((i < 1) || (j<1) || (i > d3.m->iNoRows) || (j > d3.m->iNoColumns)) {
        printf("ERROR >> In Quantity_Extract() : Array cpMatrixName = %s\n", d3.m->cpMatrixName);
        printf("ERROR >> In Quantity_Extract() : Array dimensions = [%4d][%4d]\n",
                                                       d3.m->iNoRows, d3.m->iNoColumns);
        printf("ERROR >> In Quantity_Extract() : Arguments [i][j] = [%4d][%4d]\n", i,j);
        FatalError(" Array arguments a[i][j] out of bounds",(char *)NULL);
     }

     if(d3.m->eType == DOUBLE_ARRAY) {
        switch( (int) d3.m->eRep ) {
            case INDIRECT:
                 d.q->value = d3.m->uMatrix.daa[i-1][j-1];
                 break;
            case SKYLINE:
                 iMin = MIN(i,j);
                 iMax = MAX(i,j);
                 if((iMax - iMin + 1) <= d3.m->uMatrix.daa[iMax-1][0])
                    d.q->value = d3.m->uMatrix.daa[iMax-1][iMax-iMin+1];
                 else
                    d.q->value = 0.0;
                 break;
            default:          
                 FatalError("In Quantity_Extract() : Cannot extract matrix type",
                            (char *) NULL);
        }

        switch( CheckUnits() ) {
          case ON:
             d.q->dimen =  UnitsMult( &(d3.m->spColUnits[j-1]), &(d3.m->spRowUnits[i-1]) );
             free((char *)d2.q->dimen->units_name);
             free((char *)d2.q->dimen);
             free((char *)d1.q->dimen->units_name);
             free((char *)d1.q->dimen);
             break;
          case OFF:
             break;
          default:          
             FatalError("In Quantity_Extract() : CheckUnits is not ON or OFF",(char *)NULL);
             break;