rowcol.cxx

不错的国外的有限元程序代码,附带详细的manual,可以节省很多的底层工作.

//   file ROWCOL.CXX

#include "rowcol.hxx"
#include "flotarry.hxx"
#include "intarray.hxx"
#include "mathfem.h"
#include "freestor.h"
#include <stdio.h>
#include <stdlib.h>


RowColumn :: RowColumn (int n, int st)
   // Constructor. Creates a row-column with number n, starting at index st.
{
   int size ;

   number = n ;
   start  = st ;
   diag   = 0. ;
   size   = number - start ;
   if (size) {
      row    = allocDouble(size) ;
      column = allocDouble(size) ;}
   else {
      row    = NULL ;
      column = NULL ;}
}


RowColumn :: ~RowColumn ()
   // Destructor.
{
   int size ;

   size = number - start ;
   if (size) {
      freeDouble(row) ;
      freeDouble(column) ;}
}


#ifdef DEBUG
double&  RowColumn :: atU (int i)
   // Returns the i-th coefficient of the column of the receiver. Slow but
   // safe.
{
   this -> checkBounds(i) ;
   return column[i-start] ;
}
#endif


#ifdef DEBUG
double&  RowColumn :: atL (int i)
   // Returns the i-th coefficient of the row of the receiver. Slow but safe.
{
   this -> checkBounds(i) ;
   return row[i-start] ;
}
#endif


void  RowColumn :: checkBounds (int i)
   // Checks that the receiver (k) possesses coefficients (i,k) and (k,i).
{
   if (i<start || i>=number) {
      printf ("error in bounds of RowColumn %d (start=%d) : no entry %d\n",
	       number,start,i) ;
      exit(0) ;}
}


void  RowColumn :: checkSizeTowards (IntArray* loc)
   // If loc points to a coefficient which out of the receiver's range,
   // enlarges the receiver.
{
   int n,i,j,c,first ;

   n = loc->giveSize() ;
   if (! n) {
      printf ("error in RowCol::checkBounds : 0-size location array\n") ;
      exit(0) ;}

   // get first non-zero coefficient in loc
   for (i=1 ; i<=n ; i++)
      if (first=loc->at(i))
	 break ;

   // get min coefficient in loc
   for (j=i+1 ; j<=n ; j++)
      if (c=loc->at(j))
	 first = min(c,first) ;

   // enlarge the receiver if necessary
   if (first < start)
      this -> growTo(first) ;
}


double  RowColumn :: dot (FloatArray* b, char c, int first, int last)
   // Returns the dot product of the row (c=='R') or the column (c=='C') of
   // the receiver with array 'b', from indices 'first' to 'last'.
{
   double answer,*p1,*p2 ;
   int    i ;

   answer = 0. ;
   i      = last - first + 1 ;
   if (i > 0) {
      if (c=='R')
	 p1 = row + first-start ;
      else
	 p1 = column + first-start ;
      p2 = b->givePointer() + first-1 ;            // bad practice !
      while (i--)
	 answer += *p1++ * *p2++ ;}

   return answer ;
}


void  RowColumn :: growTo (int newStart)
   // Enlarges the receiver (the column upwards, the row leftwards) to
   // 'start' = newStart.
{
   int    newSize,size,i ;
   double *newRow,*newColumn,*p1,*p2 ;

#  ifdef DEBUG
      if (newStart<=0 || newStart>=start) {
	 printf ("cannot enlarge RowCol %d (start=%d) to %d\n",
		  number,start,newStart) ;
	 exit(0) ;}
#  endif

   newSize   = number - newStart ;
   newRow    = allocDouble(newSize) ;
   newColumn = allocDouble(newSize) ;
   size      = number - start ;

   if (size) {
      p1  = row ;
      p2  = newRow + (start-newStart) ;
      i   = size ;
      while (i--)
	 *p2++ = *p1++ ;
      freeDouble(row) ;

      p1  = column ;
      p2  = newColumn + (start-newStart) ;
      i   = size ;
      while (i--)
	 *p2++ = *p1++ ;
      freeDouble(column) ;}

   row    = newRow ;
   column = newColumn ;
   start  = newStart ;
}


void  RowColumn :: printYourself ()
   // Prints the receiver on screen.
{
   int i ;

   printf ("Row-column %d : start = %d, diag = %.5f\n   col : ",
	    number,start,diag) ;
   for (i=0 ; i<number-start ; i++)
      printf (" % .5f",column[i]) ;
   printf ("\n   row : ") ;
   for (i=0 ; i<number-start ; i++)
      printf (" % .5f",row[i]) ;
   printf ("\n") ;
}


RowColumn*  RowColumn :: reinitialized ()
   // Sets all coefficients of the receiver to 0.
{
   int size = number-start ;

   diag = 0. ;
   if (size) {
      freeDouble(row) ;
      freeDouble(column) ;
      row    = allocDouble(size) ;
      column = allocDouble(size) ;}
   return this ;
}