matrix.c

有限元程序

      m2->uMatrix.daa[0][0]  = m1->iNoRows;
      m2->uMatrix.daa[0][1]  = m1->iNoColumns;
      if(CheckUnits()==ON) {
         ZeroUnits(&(m2->spRowUnits[0]));
         ZeroUnits(&(m2->spColUnits[0]));
         ZeroUnits(&(m2->spColUnits[1]));
      }

      return (m2);
}

/* 
 *  ================================================= 
 *  MatrixScale() : Scale a matrix by a double factor
 *  ================================================= 
 */ 

#ifdef __STDC__
MATRIX *MatrixScale(MATRIX *spA, double factor )
#else
MATRIX *MatrixScale(spA, factor )
MATRIX  *spA;
double  factor;
#endif
{
MATRIX *spB;

#ifdef DEBUG
       printf("*** Enter MatrixScale() : spA->iNoRows    = %4d\n", spA->iNoRows);
       printf("                        : spA->iNoColumns = %4d\n", spA->iNoColumns);
#endif
 
       switch((int) spA->eRep) {
           case INDIRECT:
                switch((int) spA->eType) {
                    case DOUBLE_ARRAY:
                         spB = MatrixScaleIndirectDouble( spA, factor );
                         break;
                    default:
                         FatalError("In MatrixScale() : Undefined spA->eType",
                                   (char *) NULL);
                         break;
                }
                break;
           case SKYLINE:
                spB = MatrixScaleSkyline( spA, factor );
                break;
          default:
                FatalError("In MatrixScale() : Undefined spA->eRep",
                          (char *) NULL);
                break;
       }

#ifdef DEBUG
       printf("*** Leave MatrixScale()\n");
#endif
       return ( spB );
}

#ifdef __STDC__
double MatrixContentScale(MATRIX *spA, int row_no, int col_no )
#else
double MatrixContentScale(spA, row_no, col_no )
MATRIX  *spA;
int row_no;        /* row number    */
int col_no;        /* column number */
#endif
{
double  da;

#ifdef DEBUG
       printf("*** Enter MatrixContentScale() : spA->iNoRows    = %4d\n", spA->iNoRows);
       printf("                        : spA->iNoColumns = %4d\n", spA->iNoColumns);
#endif
      if(CheckUnits()==OFF) 
         FatalError("You have to set units ON to use this function","In MatrixScale",(char *)NULL );
 
     /* [a] Scale for Column of matrix spA */

       switch((int) spA->eRep) {
           case INDIRECT:
                switch((int) spA->eType) {
                    case DOUBLE_ARRAY:
                         da = MatrixContentScaleIndirectDouble( spA, row_no, col_no );
                         break;
                    default:
                         FatalError("In MatrixContentScale() : Undefined spA->eType",
                                   (char *) NULL);
                         break;
                }
                break;
           case SKYLINE:
                da = MatrixContentScaleSkyline( spA, row_no, col_no );
                break;
          default:
                FatalError("In MatrixContentScale() : Undefined spA->eRep",
                          (char *) NULL);
                break;
       }

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

       return ( da );
}

/* ======================================
   Operations between QUANTITY and MATRIX
   ====================================== */

#ifdef __STDC__
MATRIX *MatrixQuanMult(QUANTITY *q1, MATRIX *m2)
#else
MATRIX *MatrixQuanMult(q1, m2)
QUANTITY *q1;
MATRIX   *m2;
#endif
{
MATRIX     *m3;
int      i,j,k;
int     length;
DIMENSIONS *d1;

   /* [a] : Multiply a quantity to a matrix */

   m3 = MatrixScale( m2, q1->value );
       
   if( CheckUnits() == ON ) {
       for(i = 1; i <= m2->iNoRows; i++) {
           UnitsCopy( &(m3->spRowUnits[i-1]), &(m2->spRowUnits[i-1]) );
           for(j = 1; j <= m2->iNoColumns; j++)
               UnitsMultRep( &(m3->spColUnits[j-1]), q1->dimen, &(m2->spColUnits[j-1]) );
       }
   }

   return (m3);
}

#ifdef __STDC__
MATRIX *MatrixQuanDiv(MATRIX *m2, QUANTITY *q1)
#else
MATRIX *MatrixQuanDiv(m2, q1)
QUANTITY *q1;
MATRIX   *m2;
#endif
{
MATRIX         *m3;
DIMENSIONS     *d1;
int          i,j,k;
int         length;

    /* [a] : matrix is divdied by a quantity */

    m3 = MatrixScale( m2, 1.0/q1->value );
       
    if( CheckUnits()==ON ) {
        for(i = 1; i <= m2->iNoRows; i++) {
            UnitsCopy( &(m3->spRowUnits[i-1]), &(m2->spRowUnits[i-1]) );
            for(j = 1; j <= m2->iNoColumns; j++)
                UnitsDivRep( &(m3->spColUnits[j-1]), &(m2->spColUnits[j-1]), q1->dimen, YES);
        }
    }

    return (m3);
}

/* 
 *  =================================================
 *  QuantityCast() : convert 1x1 matrix into quantity 
 *  =================================================
 */ 

#ifdef __STDC__
QUANTITY *QuantityCast( MATRIX *m )
#else
QUANTITY *QuantityCast(m)
MATRIX   *m;
#endif
{
QUANTITY         *q;
int          length;

   /* [a] : Check that matrix is either a column (or row) vector */

   if((m->iNoRows != 1) || (m->iNoColumns != 1)) {
       FatalError("In Quantity_Cast() : Matrix is not a 1x1 Matrix",(char *)NULL);
   }

   q            = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));
   q->value     = m->uMatrix.daa[0][0];

   switch(CheckUnits()) {
      case ON:
           q->dimen = UnitsMult( &(m->spRowUnits[0]), &(m->spColUnits[0]) );
           break;
      case OFF:
           q->dimen = (DIMENSIONS *)NULL;
           break;
   }
     
   return (q);
}


/* 
 *  ==========================================================
 *  MatrixColumnUnits() : Specify Column Units in a Matrix
 *  
 *  Usage:  1st Argument is pointer to matrix concerned.
 *          2nd Argument is matrix of units to be applied.
 *          3rd Argument column number for units.
 *  
 *  Examples : matrix = ColumnUnits( matrix, [ kg ] );
 *             matrix = ColumnUnits( matrix, [ kg ], [2] );
 *  ==========================================================
 */ 

#ifdef  __STDC__
MATRIX *MatrixColumnUnits(MATRIX *first_matrix, ...) {
#else
MATRIX *MatrixColumnUnits(va_alist)
va_dcl
{
MATRIX     *first_matrix;
#endif

va_list          arg_ptr;
MATRIX                *p;
MATRIX *m, *m1, *units_m;
int     i, j, k, counter;
int    NO_COL_WITH_UNITS;
int NO_COL_WITH_NO_UNITS;
int               *COUNT;
int           ID, COL_NO;
int               length;
int       iFLAG  = FALSE;
int       iFLAG1 = FALSE;

   /* [a] : Check units flags and get appropriate units */

   if( CheckUnits() == OFF) {
       FatalError("In MatrixColumnUnits () : To call this function, units must be set to ON ",
                 (char *)NULL );
   }

   /* [b] : Get appropriate units */

#ifdef  __STDC__
    va_start(arg_ptr, first_matrix);
#else
    va_start(arg_ptr);
    first_matrix = va_arg(arg_ptr, MATRIX *);
#endif

    m1 = first_matrix;
    units_m = va_arg(arg_ptr, MATRIX *);
    p       = va_arg(arg_ptr, MATRIX *);
    m       = MatrixCopyIndirectDouble(m1);
    va_end(arg_ptr);

   /* [c] : Retrieve and check appropriate units */

   if( p != NULL) {
       ID = 1; 
       COL_NO  = p->uMatrix.daa[0][0]; 
   } else
       ID = 2;

   if(units_m->iNoRows != 1) 
      FatalError("Fatal Error in MatrixColumnUnits(): Units_matrix should be a 1xn vector",
                (char *) NULL); 

   /* [d] : Case 1 : Update "column units" in one column only.              */
   /*                                                                       */
   /*       Two cases exist -- (1.1) Column i of matrix m has no units      */
   /*                          (1.2) Column i of matrix m has units.        */
   /*                                Check compatibility with units_m.      */
   /*                                Copy unmits_m to m.                    */
   /*                                                                       */
   /*       Examples : stiff = Matrix([2,2]);                               */
   /*                  stiff = ColumnUnits ( stiff, [  sec], [1] );         */
   /*                  stiff = ColumnUnits ( stiff, [sec^2], [2] );         */
   /*                                                                       */

   if (ID == 1) {

      if(units_m->iNoColumns != 1) 
         FatalError("In MatrixColumnUnits(): too many units for one column ",(char *)NULL);

      i = COL_NO; /* column no i is to be assigned an unit  */

      if(m->spColUnits[i-1].length_expnt == 0 &&
         m->spColUnits[i-1].time_expnt == 0   &&
         m->spColUnits[i-1].mass_expnt == 0   &&
         m->spColUnits[i-1].temp_expnt == 0) {

         /* Case (1.1) : Assign units from units_m to m */
          
         UnitsCopy( &(m->spColUnits[i-1]), &(units_m->spColUnits[0]) );
         for(j = 1; j <= m->iNoRows; j++)
             m->uMatrix.daa[j-1][i-1] = 
                units_m->spColUnits[0].scale_factor*m->uMatrix.daa[j-1][i-1];
      } else {

         /* Case (1.2) : Check with units_m, and copy units_m to units of m */
           
         if( SameUnits( &(m->spColUnits[i-1]), &(units_m->spColUnits[0])) == TRUE ){
             UnitsCopy( &(m->spColUnits[i-1]), &(units_m->spColUnits[0]) );
         } else {
             FatalError("In MatrixColumnUnits(): Inconsistent units",
                       (char *) NULL );
         }
      }
   }

   /* [e] : Case 2 : Update "units" in multiple matrix columns.                      */
   /*                                                                                */
   /*       Cases -- (2.1) Matrix m and units_m have the same number of columns.     */
   /*                (2.1.1) : Matrix m has no units. Assign units from units_m to m */
   /*                (2.1.2) : Matrix m has units. Copy units_m to units of m.       */
   /*                (2.2) Matrices m and units_m have different number of columns.  */
   /*                (2.2.1) : matrix m rows have no units. Assign all columns of m  */
   /*                          with units of units_m                                 */
   /*                (2.2.2) : Assign unit_m to columns of m where they match.       */
   /*                                                                                */
   /*       Examples : stiff = Matrix([3,3]);                                        */
   /*                  stiff = ColumnUnits ( stiff, [cm/sec^2], [1] );               */
   /*                  stiff = ColumnUnits ( stiff, [sec] );                         */
   /*                  stiff = ColumnUnits ( stiff, [in/sec^2] );                    */
   /*                                                                                */

   if (ID == 2) {

      /* Case (2.1) : Matrices m and units_m have same number of columns */

      if(units_m->iNoColumns == m->iNoColumns) {
         for (i = 1; i <= m->iNoColumns; i++) { 

            if(m->spColUnits[i-1].length_expnt == 0 && 
               m->spColUnits[i-1].time_expnt == 0   &&
               m->spColUnits[i-1].mass_expnt == 0   &&
               m->spColUnits[i-1].temp_expnt == 0) {

               /* Case (2.1.1) : Matrix m has no units -- Assign units from units_m to m */

               UnitsCopy( &(m->spColUnits[i-1]), &(units_m->spColUnits[i-1]) );
               for(j = 1; j <= m->iNoRows; j++)
                   m->uMatrix.daa[j-1][i-1] = 
                      units_m->spColUnits[i-1].scale_factor*m->uMatrix.daa[j-1][i-1];

            } else {

               /* Case (2.1.2) : Matrix m has units. Check, and copy units_m to units of m */
           
               if(SameUnits(&(m->spColUnits[i-1]), &(units_m->spColUnits[i-1])) == TRUE ){
                  UnitsCopy( &(m->spColUnits[i-1]), &(units_m->spColUnits[i-1]) );
               } else {
                  FatalError("Fatal Error in MatrixColumnUnits(): inconsistent units ",
                            (char *)NULL);
               }
            }
         }
      }

      /* Case (2.2) : Matrices m and units_m have different number of columns */

      if(units_m->iNoColumns != m->iNoColumns) {

         k = 0;
         COUNT = (int *) MyCalloc(m->iNoColumns, sizeof(int));
         counter = 0;

         /* Find and record column numbers with matching units */

         for (i = 1; i <= m->iNoColumns; i++) {
            if(m->spColUnits[i-1].length_expnt != 0 ||
               m->spColUnits[i-1].time_expnt   != 0 ||
               m->spColUnits[i-1].mass_expnt   != 0 ||
               m->spColUnits[i-1].temp_expnt   != 0) {
                k = k + 1;
                COUNT[k-1] = i; 
            }
         }

         NO_COL_WITH_UNITS = k;
         NO_COL_WITH_NO_UNITS = m->iNoColumns - NO_COL_WITH_UNITS;

         if( NO_COL_WITH_NO_UNITS == m->iNoColumns) {

            /* Case (2.2.1) : matrix m rows have no units. Assign all columns of m */
            /*                with units of units_m                                */

            if(units_m->iNoColumns != 1)
               FatalError(" in MatrixColumnUnits(): Too many/few units ",
                         (char *)NULL);

            if(units_m->iNoColumns == 1) { 
               for(i = 1; i <= m->iNoColumns; i++) {
                   UnitsCopy( &(m->spColUnits[i-1]), &(units_m->spColUnits[0]) );
                   for(j = 1; j <= m->iNoRows; j++)
                       m->uMatrix.daa[j-1][i-1] *= units_m->spColUnits[0].scale_factor;
               }
            }
         }

         if( NO_COL_WITH_NO_UNITS != m->iNoColumns) {

            /* Case (2.2.2) : Assign unit_m to columns of m where they match */

            for( k = 1; k <=  NO_COL_WITH_UNITS; k++) {
               i = COUNT[k-1];
               for(j = 1; j <= units_m->iNoColumns; j++) {
                   if(SameUnits(&(m->spColUnits[i-1]), &(units_m->spColUnits[j-1])) == TRUE) {
                      UnitsCopy( &(m->spColUnits[i-1]), &(units_m->spColUnits[j-1]) );
                      iFLAG1 = TRUE;
                   } else 
                      iFLAG1 = MAX(iFLAG1,