element.cxx

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

      return resultant ;
}


FloatArray*  Element :: ComputeRhsAt (TimeStep* stepN)
   // Computes the contribution of the receiver to the right-hand side of the
   // linear system.
{
   TimeIntegrationScheme* scheme ;

   scheme = domain -> giveTimeIntegrationScheme() ;
   if (scheme -> isStatic())
      return  this -> ComputeStaticRhsAt (stepN) ;
   else if (scheme -> isNewmark())
      return  this -> ComputeNewmarkRhsAt(stepN) ;
   else {
      printf ("Error : unknown time integration scheme : %c\n",scheme) ;
      exit(0) ;}
}


FloatMatrix*  Element :: ComputeStaticLhsAt (TimeStep* stepN)
   // Computes the contribution of the receiver to the left-hand side of the
   // linear system, in a static analysis.
{
   return  this->giveStiffnessMatrix()->GiveCopy() ;
}


FloatArray*  Element :: ComputeStaticRhsAt (TimeStep* stepN)
   // Computes the contribution of the receiver to the right-hand side of the
   // linear system, in a static analysis.
{
   return  this->ComputeLoadVectorAt(stepN) ;
}


FloatMatrix*  Element :: computeStiffnessMatrix ()
   // Computes numerically the stiffness matrix of the receiver.
{
   int         i ;
   double      dV ;
   FloatMatrix *b,*db,*d ;
   GaussPoint  *gp ;

   stiffnessMatrix = new FloatMatrix() ;
   for (i=0 ; i<numberOfGaussPoints ; i++) {
      gp = gaussPointArray[i] ;
      b  = this -> ComputeBmatrixAt(gp) ;
      d  = this -> giveConstitutiveMatrix() ;
      dV = this -> computeVolumeAround(gp) ;
      db = d -> Times(b) ;
      stiffnessMatrix -> plusProduct(b,db,dV) ;

      delete b ;
      delete db ;}

   return stiffnessMatrix -> symmetrized() ;
}


FloatArray*  Element :: computeStrainVector (GaussPoint* gp, TimeStep* stepN)
   // Computes the vector containing the strains at the Gauss point gp of
   // the receiver, at time step stepN. The nature of these strains depends
   // on the element's type.
{
   FloatMatrix *b ;
   FloatArray  *u,*Epsilon ;

   b       = this -> ComputeBmatrixAt(gp) ;
   u       = this -> ComputeVectorOf('d',stepN) ;
   Epsilon = b -> Times(u) ;

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


FloatArray*  Element :: computeStressVector (GaussPoint* gp, TimeStep* stepN)
   // Computes the vector containing the stresses at the Gauss point gp of
   // the receiver, at time step stepN. The nature of these stresses depends
   // on the element's type.
{
   FloatArray *Epsilon,*Sigma ;

   Epsilon = gp -> giveStrainVector() ;
   if (! Epsilon)
      Epsilon = this->computeStrainVector(gp,stepN) ;
   Sigma = this->giveConstitutiveMatrix() -> Times(Epsilon) ;

   gp -> letStressVectorBe(Sigma) ;      // gp stores Sigma, not a copy
   return Sigma ;
}


FloatArray* Element :: ComputeVectorOf (char u, TimeStep* stepN)
   // Forms the vector containing the values of the unknown 'u' (e.g., the
   // displacement) of the dofs of the receiver's nodes.
{
   Node       *nodeI ;
   FloatArray *answer ;
   int        i,j,k,nDofs ;

   answer = new FloatArray(this->computeNumberOfDofs()) ;
   k      = 0 ;
   for (i=1 ; i<=numberOfNodes ; i++) {
      nodeI = this->giveNode(i) ;
      nDofs = nodeI->giveNumberOfDofs() ;
      for (j=1 ; j<=nDofs ; j++)
	 answer->at(++k) = nodeI->giveDof(j)->giveUnknown(u,stepN) ;}

   return answer ;
}


FloatArray* Element :: ComputeVectorOfPrescribed (char u, TimeStep* stepN)
   // Forms the vector containing the prescribed values of the unknown 'u'
   // (e.g., the prescribed displacement) of the dofs of the receiver's
   // nodes. Puts 0 at each free dof.
{
   Node       *nodeI ;
   Dof        *dofJ ;
   FloatArray *answer ;
   int        i,j,k,nDofs ;

   answer = new FloatArray(this->computeNumberOfDofs()) ;
   k      = 0 ;
   for (i=1 ; i<=numberOfNodes ; i++) {
      nodeI = this->giveNode(i) ;
      nDofs = nodeI->giveNumberOfDofs() ;
      for (j=1 ; j<=nDofs ; j++) {
	 dofJ = nodeI->giveDof(j) ;
	 if (dofJ -> hasBc())
	    answer->at(++k) = dofJ->giveUnknown(u,stepN) ;
	 else
	    answer->at(++k) = 0. ;}}

   return answer ;
}


IntArray*  Element :: giveBodyLoadArray ()
   // Returns the array which contains the number of every body load that act
   // on the receiver.
{
   int i,numberOfLoads ;

   if (! bodyLoadArray) {
      numberOfLoads = this -> readIfHas("bodyLoads") ;
      bodyLoadArray = new IntArray(numberOfLoads) ;
      for (i=1 ; i<=numberOfLoads ; i++)
	 bodyLoadArray->at(i) = this->readInteger("bodyLoads",i+1) ;}

   return bodyLoadArray ;
}


FloatMatrix*  Element :: giveConstitutiveMatrix ()
   // Returns the elasticity matrix {E} of the receiver.
{
   if (! constitutiveMatrix)
      this -> computeConstitutiveMatrix() ;

   return constitutiveMatrix ;
}


IntArray*  Element :: giveLocationArray ()
   // Returns the location array of the receiver. This array is obtained by
   // simply appending the location array of every node of the receiver.
{
   IntArray* nodalArray ;
   int       i ;

   if (! locationArray) {
      locationArray = new IntArray(0) ;
      for (i=1 ; i<=numberOfNodes ; i++) {
	 nodalArray    = this -> giveNode(i) -> giveLocationArray() ;
	 locationArray = locationArray -> followedBy(nodalArray) ;}}

   return locationArray ;
}


FloatMatrix*  Element :: giveMassMatrix ()
   // Returns the mass matrix of the receiver.
{
   if (! massMatrix)
      this -> computeMassMatrix() ;
   return massMatrix ;
}


Material*  Element :: giveMaterial ()
   // Returns the material of the receiver.
{
   if (! material)
      material = this -> readInteger("mat") ;

   return  domain -> giveMaterial(material) ;
}



Node*  Element :: giveNode (int i)
   // Returns the i-th node of the receiver.
{
   int  n ;

   if (! nodeArray)
      nodeArray = new IntArray(numberOfNodes) ;

   n = nodeArray->at(i) ;
   if (! n) {
      n = this -> readInteger("nodes",i) ;
      nodeArray->at(i) = n ;}

   return  domain -> giveNode(n) ;
}


FloatMatrix*  Element :: giveStiffnessMatrix ()
   // Returns the stiffness matrix of the receiver.
{
   if (! stiffnessMatrix)
      this -> computeStiffnessMatrix() ;
   return stiffnessMatrix ;
}


void  Element :: instanciateYourself ()
   // Gets from input file all data of the receiver.
{
   int i ;

#  ifdef VERBOSE
      printf ("instanciating element %d\n",number) ;
#  endif

   material = this -> readInteger("mat") ;
   nodeArray = new IntArray(numberOfNodes) ;
   for (i=1 ; i<=numberOfNodes ; i++)
      nodeArray->at(i) = this->readInteger("nodes",i) ;
   this -> giveBodyLoadArray() ;
}


Element*  Element :: ofType (char* aClass)
   // Returns a new element, which has the same number than the receiver,
   // but belongs to aClass (PlaneStrain, or Truss2D,..).
{
   Element* newElement ;

   if (! strncmp(aClass,"PN",2))
      newElement = new PlaneStrain(number,domain) ;
   else if (! strncmp(aClass,"T2",2))
      newElement = new Truss2D(number,domain) ;
   else if (! strncmp(aClass,"Beam2D",2))
      newElement = new Beam2D(number,domain) ;
   else {
      printf ("%s : unknown element type \n",aClass) ;
      exit(0) ;}

   return newElement ;
}


void  Element :: printOutputAt (TimeStep* stepN)
   // Performs end-of-step operations.
{
   int         i ;
   GaussPoint* gp ;
   FILE*       file ;

#  ifdef VERBOSE
      printf ("element %d printing output\n",number) ;
#  endif

   file = domain -> giveOutputStream() ;
   fprintf (file,"element %d :\n",number) ;

   for (i=1 ; i<=numberOfGaussPoints ; i++) {
      gp = gaussPointArray[i-1] ;
      this -> computeStrainVector(gp,stepN) ;
      this -> computeStressVector(gp,stepN) ;
      gp   -> printOutput() ;}
}


Element*  Element :: typed ()
   // Returns a new element, which has the same number than the receiver,
   // but is typed (PlaneStrain, or Truss2D,..).
{
   Element* newElement ;
   char     type[32] ;

   this -> readString("class",type) ;
   newElement = this -> ofType(type) ;

   return newElement ;
}


void  Element :: updateYourself ()
   // Updates the receiver at end of step.
{
   int i ;

#  ifdef VERBOSE
      printf ("updating element %d\n",number) ;
#  endif

   for (i=1 ; i<=numberOfGaussPoints ; i++)
      gaussPointArray[i-1] -> updateYourself() ;
}