skylineu.cxx

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

//   file SKYLINEU.CXX

#include "skylineu.hxx"
#include "rowcol.hxx"
#include "domain.hxx"
#include "element.hxx"
#include "flotmtrx.hxx"
#include "flotarry.hxx"
#include "intarray.hxx"
#include "mathfem.h"
#include "debug.def"
#include <math.h>


SkylineUnsym :: SkylineUnsym ()
   // Constructor. Creates an empty skyline.
{
   size         = 0 ;
   rowColumns   = NULL ;
   isFactorized = FALSE ;
}


SkylineUnsym :: ~SkylineUnsym ()
   // Destructor.
{
   RowColumn** p ;
   int         i ;

   if (size) {
      i = size ;
      p = rowColumns ;
      while (i--)
	 delete *p++ ;
      delete [size] rowColumns ;}
}


FloatMatrix*  SkylineUnsym :: AsFloatMatrix ()
   // Returns a copy of the receiver, instance of FloatMatrix. Useful for
   // printing when debugging.
{
   FloatMatrix* answer ;
   int          i,j,start ;

   answer = new FloatMatrix(size,size) ;
   for (j=1 ; j<=size ; j++) {
      start = this->giveRowColumn(j)->giveStart() ;
      for (i=start ; i<=j ; i++) {
	 answer->at(i,j) = this->at(i,j) ;
	 answer->at(j,i) = this->at(j,i) ;}}
   return answer ;
}


void  SkylineUnsym :: assemble (FloatMatrix* k, IntArray* loc)
   // Assembles the elemental matrix 'k' to the receiver, using 'loc' as a
   // location array. The values in k corresponding to a zero coefficient
   // in loc are not assembled.
   // Warning : k is not supposed to be an instance of DiagonalMatrix.
{
   int        i,j,ii,jj,dim ;
   RowColumn* rowColumnJJ ;

#  ifdef DEBUG
      dim = k -> giveNumberOfRows() ;
      if (dim != loc->giveSize()) {
	 printf ("error : dimension of 'k' and 'loc' mismatch \n") ;
	 exit(0) ;}

      this -> checkSizeTowards(loc) ;
#  endif

   if (k->isDiagonal()) {
      printf ("error in SkylineUnsym::assemble ; k cannot be diagonal\n") ;
      exit(0) ;}

   dim = k->giveNumberOfRows() ;

   for (j=1; j<=dim; j++) {
	   jj = loc->at(j);
	   if(jj){
		   for (i=1; i<=dim; i++) {
			   ii = loc->at(i);
			   if(ii){
				   this->at(ii,jj) += k->at(i,j); 
			   }
		   }
	   }
   }

}


double&  SkylineUnsym :: at (int i, int j)
   // Returns the (i,j) coefficient of the receiver. Not very efficient.
{
   if (i < j)
      return  this->giveRowColumn(j)->atU(i) ;
   else if (i > j)
      return  this->giveRowColumn(i)->atL(j) ;
   else
      return  this->giveRowColumn(i)->atDiag() ;
}


FloatArray*  SkylineUnsym :: backSubstitutionWith (FloatArray* y)
   // Returns the solution x of the system U.x = y , where U is the upper
   // half of the receiver. Note : x overwrites y.
{
   int        i,k,start ;
   double     yK ;
   RowColumn* rowColumnK ;

   for (k=size ; k>0 ; k--) {
      rowColumnK = this->giveRowColumn(k) ;
      yK       = y->at(k) ;
      start    = rowColumnK->giveStart() ;
      for (i=start ; i<k ; i++)
	 y->at(i) -= rowColumnK->atU(i) * yK ;}

   return y ;
}


void  SkylineUnsym :: carveYourselfFor (Domain* aDomain)
   // Instanciates the profile of the receiver and initializes all coeffi-
   // cients to zero.
   // Warning : case diagonal (lumped) matrix to expected.
{
   IntArray* LocArray ;
   int       n,nElem,i,ii ;

   nElem = aDomain -> giveNumberOfElements() ;
   for (n=1 ; n<=nElem ; n++) {
      LocArray = aDomain -> giveElement(n) -> giveLocationArray() ;
      for (i=1 ; i <= LocArray->giveSize() ; i++) {
	 ii = LocArray->at(i) ;
	 if (ii)
	    this -> giveRowColumn(ii) -> checkSizeTowards(LocArray) ;}}
}


void  SkylineUnsym :: checkSizeTowards (IntArray*loc)
   // Increases the number of columns of the receiver if 'loc' points to
   // out-of-range columns.
{
   int i,dim,maxCol ;

   maxCol = 0 ;                           // the largest coefficient in 'loc'
   dim    = loc->giveSize() ;
   for (i=1 ; i<=dim ; i++)
      maxCol = max (maxCol,loc->at(i)) ;

   if (maxCol > size)                     // enlarge the matrix
      this -> growTo(maxCol) ;

   for (i=1 ; i<=size ; i++)
      this -> giveRowColumn(i) -> checkSizeTowards(loc) ;
}


Skyline*  SkylineUnsym :: diagonalScalingWith (FloatArray* y)
   // Scales y by the diagonal of the receiver. Returns the receiver.
{
   double       diag ;
   int          k ;
   const double precision = 1.e-10 ;

   for (k=1 ; k<=size ; k++) {
      diag = this->giveRowColumn(k)->atDiag() ;
#     ifdef DEBUG
	 if (fabs(diag) < precision) {
	   printf ("error in SkylineUnsym::diagScaling : pivot %d is zero \n",k);
	   exit(0) ;}
#     endif
      y->at(k) /= diag ;}

   return this ;
}


Skyline*  SkylineUnsym :: factorized ()
   // Returns the receiver in L.D.U factorized form. From Golub & van Loan,
   // 1rst edition, pp 83-84.
{
   RowColumn     *rowColumnK,*rowColumnI ;
   FloatArray    *r,*w ;
   double        diag ;
   int           k,i,p,start,startK,startI ;
   const double  PRECISION = 1.e-10 ;

   if (isFactorized)
      return this ;

   if (! size) {
      printf ("error : cannot factorize null-sized matrix\n") ;
      exit(0) ;}

   for (k=1 ; k<=size ; k++) {

      rowColumnK = this -> giveRowColumn(k) ;
      startK   = rowColumnK -> giveStart() ;

      // compute vectors r and w
      r = new FloatArray(k-1) ;
      w = new FloatArray(k-1) ;
      for (p=startK ; p<k ; p++) {
	 diag     = this->giveRowColumn(p)->atDiag() ;
	 r->at(p) = diag * rowColumnK->atU(p) ;
	 w->at(p) = diag * rowColumnK->atL(p) ;}

      // compute diagonal coefficient of rowColumn k
      rowColumnK->atDiag() -= rowColumnK->dot(r,'R',startK,k-1) ;
      diag = rowColumnK->atDiag() ;

      // test pivot not too small
      if (fabs(diag) < PRECISION) {
	 printf ("too small pivot %d\n",k) ;
	 exit(0) ;}

      // compute off-diagonal coefficients of rowColumns i>k
      for (i=k+1 ; i<=size ; i++) {
	 rowColumnI = this->giveRowColumn(i) ;
	 startI   = rowColumnI->giveStart() ;
	 if (startI <= k) {
	    start = max(startI,startK) ;
	    rowColumnI->atL(k) -= rowColumnI->dot(r,'R',start,k-1) ;
            rowColumnI->atL(k) /= diag ;
	    rowColumnI->atU(k) -= rowColumnI->dot(w,'C',start,k-1) ;
            rowColumnI->atU(k) /= diag ;}}

      delete r ;
      delete w ;}

   isFactorized = TRUE ;
   return this ;
}


Skyline*  SkylineUnsym :: forwardReductionWith (FloatArray* b)
   // Computes y, the solution of the system  L.y = b , where L is the lower
   // half of the receiver (assumed to be in a factorized form). Returns the 
   // receiver. Note: y overwrites b.
{
   int        k,start ;
   RowColumn* rowColumnK ;

#  ifdef DEBUG
      if (! isFactorized) {
	 printf ("error : a rowColumn must be factorized before reductions\n");
	 exit(0) ;}
#  endif

   if (! size)
      return this ;                                      // null size system
   if (b->giveSize() != size) 
      b -> growTo(size) ;

   for (k=1 ; k<=size ; k++) {
      rowColumnK  = this->giveRowColumn(k) ;
      start     = rowColumnK->giveStart() ;
      b->at(k) -= rowColumnK->dot(b,'R',start,k-1) ;}

   return this ;
}


RowColumn*  SkylineUnsym :: giveRowColumn (int j)
   // Returns the j-th rowColumn of the receiver. Creates it if necessary.
{
   if (size < j)
      this -> growTo(j) ;
   if (! rowColumns[j-1])
      rowColumns[j-1] = new RowColumn(j,j) ;

   return rowColumns[j-1] ;
}


void  SkylineUnsym :: growTo (int n)
   // Enlarges the receiver to n columns.
{
   RowColumn **newRowColumns,**p1,**p2 ;
   int       i ;

#  ifdef DEBUG
      if (n <= size) {
	 printf ("error : SkylineUnsym(%d) cannot grow to size %d\n",size,n) ;
	 exit(0) ;}
#  endif

   newRowColumns = new RowColumn* [n] ;
   p1          = rowColumns ;
   p2          = newRowColumns ;

   i = size ;
   while (i--)
      *p2++ = *p1++ ;

   i = n - size ;
   while (i--)
      *p2++ = NULL ;

   if (rowColumns)
      delete [n] rowColumns ;
   rowColumns = newRowColumns ;
   size     = n ;
}


void  SkylineUnsym :: printYourself ()
  // Prints the receiver on screen.
{
  FloatMatrix* copy ;

  copy = this -> AsFloatMatrix() ;
  copy -> printYourself() ;
  delete copy ;
}


Skyline*  SkylineUnsym :: reinitialized ()
   // Returns the receiver with all coefficients set to zero.
{
   int j ;

   for (j=1 ; j<=size ; j++)
      this -> giveRowColumn(j) -> reinitialized() ;
   isFactorized = FALSE ;

   return this ;
}