beam2d.cxx

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

//   file BEAM2D.CXX

#include "beam2d.hxx"
#include "node.hxx"
#include "material.hxx"
#include "gausspnt.hxx"
#include "flotmtrx.hxx"
#include <math.h>


Beam2D :: Beam2D (int n, Domain* aDomain)
       : Element (n,aDomain)
   // Constructor.
{
   numberOfNodes       = 2 ;
   rotationMatrix      = NULL ;
   length              = 0. ;
   pitch               = 10. ;          // a dummy value
   this -> computeGaussPoints() ;
}


FloatMatrix*  Beam2D :: ComputeBmatrixAt (GaussPoint* aGaussPoint)
   // Returns the strain matrix of the receiver.
{
   double       b,ksi,n1x,n4x,n2xx,n3xx,n5xx,n6xx ;
   FloatMatrix* answer ;

   b    = this->giveLength() / 2. ;
   ksi  = aGaussPoint->giveCoordinate(1) ;
   n1x  = -0.5 / b ;
   n4x  = -n1x ;
   n2xx = 1.5 * ksi / (b*b) ;
   n3xx = (-1.+3.*ksi) / (b+b) ;
   n5xx = -n2xx ;
   n6xx = (1.+3.*ksi) / (b+b) ;

   answer = new FloatMatrix(2,6) ;
   answer->at(1,1) =  n1x  ;
   answer->at(1,4) =  n4x  ;
   answer->at(2,2) = -n2xx ;
   answer->at(2,3) = -n3xx ;
   answer->at(2,5) = -n5xx ;
   answer->at(2,6) = -n6xx ;

   return answer ;
}


FloatMatrix*  Beam2D :: computeConstitutiveMatrix ()
   // Returns the 2-by-2 elasticity matrix {E} of the receiver.
{
   Material* mat ;
   double    e,a,i ;

   mat = this->giveMaterial() ;
   e = mat->give('E') ;
   a = mat->give('A') ;
   i = mat->give('I') ;

   constitutiveMatrix = new FloatMatrix(2,2) ;
   constitutiveMatrix->at(1,1) = (e*a) ;
   constitutiveMatrix->at(2,2) = (e*i) ;

   return constitutiveMatrix ;
}


void  Beam2D :: computeGaussPoints ()
   // Sets up the array of Gauss Points of the receiver. 
{
   double a = 1. / sqrt(3.) ;
   FloatArray *coord1, *coord2 ;

   numberOfGaussPoints = 2 ;
   gaussPointArray     = new GaussPoint* [numberOfGaussPoints] ;

   coord1 = new FloatArray(1) ;
   coord1->at(1) = -a ;
   gaussPointArray[0] = new GaussPoint(this,1,coord1,1.) ;

   coord2 = new FloatArray(1) ;
   coord2->at(1) = a ;
   gaussPointArray[1] = new GaussPoint(this,2,coord2,1.) ;
}


FloatMatrix*  Beam2D :: computeStiffnessMatrix ()
   // Returns the stiffness matrix of the receiver, expressed in the global
   // axes.
{
   this -> Element::computeStiffnessMatrix() ;
   this -> giveRotationMatrix() ;
   return  stiffnessMatrix -> rotatedWith(rotationMatrix) ;
}


FloatArray*  Beam2D :: computeStrainVector (GaussPoint* gp, TimeStep* stepN)
   // Returns the array with the strain and the curvature of the receiver.
{
   FloatMatrix *b,*rot ;
   FloatArray  *u,*epsilon ;

   b       = this -> ComputeBmatrixAt (gp) ;
   u       = this -> ComputeVectorOf ('d',stepN) ;
   rot     = this -> giveRotationMatrix () ;
   epsilon = b -> Times (u->rotatedWith(rot,'n')) ;

   gp -> letStrainVectorBe(epsilon) ;       // gp stores epsilon, not a copy
   delete b ;
   delete u ;
   return epsilon ;
}


double  Beam2D :: computeVolumeAround (GaussPoint* aGaussPoint)
   // Returns the half-length of the receiver.
{
   return  this->giveLength() / 2. ;
}


double  Beam2D :: giveLength ()
   // Returns the length of the receiver.
{
   double dx,dy ;
   Node   *nodeA,*nodeB ;

   if (length == 0.) {
      nodeA  = this->giveNode(1) ;
      nodeB  = this->giveNode(2) ;
      dx     = nodeB->giveCoordinate(1) - nodeA->giveCoordinate(1) ;
      dy     = nodeB->giveCoordinate(2) - nodeA->giveCoordinate(2) ;
      length = sqrt(dx*dx + dy*dy) ;}

   return length ;
}


double  Beam2D :: givePitch ()
   // Returns the angle in radians between the global X axis and the axis
   // of the receiver (the orientation of the receiver is node1->node2).
{
   Node   *nodeA, *nodeB ;
   double xA,xB,yA,yB ;

   if (pitch == 10.) {
      nodeA = this->giveNode(1) ;
      nodeB = this->giveNode(2) ;
      xA = nodeA->giveCoordinate(1) ;
      xB = nodeB->giveCoordinate(1) ;
      yA = nodeA->giveCoordinate(2) ;
      yB = nodeB->giveCoordinate(2) ;
      pitch = atan2(yB-yA,xB-xA) ;}

   return pitch ;
}


FloatMatrix*  Beam2D :: giveRotationMatrix ()
   // Returns the rotation matrix of the receiver. Computes it if it does
   // not exist yet.
{
   double sine,cosine ;

   if (! rotationMatrix) {
      sine           = sin(this->givePitch()) ;
      cosine         = cos(pitch) ;
      rotationMatrix = new FloatMatrix(6,6) ;
      rotationMatrix -> at(1,1) =  cosine ;
      rotationMatrix -> at(1,2) =  sine   ;
      rotationMatrix -> at(2,1) = -sine   ;
      rotationMatrix -> at(2,2) =  cosine ;
      rotationMatrix -> at(3,3) =  1.     ;
      rotationMatrix -> at(4,4) =  cosine ;
      rotationMatrix -> at(4,5) =  sine   ;
      rotationMatrix -> at(5,4) = -sine   ;
      rotationMatrix -> at(5,5) =  cosine ;
      rotationMatrix -> at(6,6) =  1.     ;}

   return rotationMatrix ;
}