element.cxx

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

//   file ELEMENT.CXX
 
#include "element.hxx"
#include "planstrn.hxx"
#include "beam2d.hxx"
#include "truss2d.hxx"
#include "domain.hxx"
#include "timestep.hxx"
#include "timinteg.hxx"
#include "node.hxx"
#include "dof.hxx"
#include "material.hxx"
#include "bodyload.hxx"
#include "gausspnt.hxx"
#include "intarray.hxx"
#include "flotarry.hxx"
#include "flotmtrx.hxx"
#include "linsyst.hxx"
#include "skyline.hxx"
#include "debug.def"
#include "verbose.def"
#include <stdlib.h>
#include <stdio.h>


Element :: Element (int n, Domain* aDomain)
	 : FEMComponent (n, aDomain)
   // Constructor. Creates an element with number n, belonging to aDomain.
{
   material           = 0    ;
   numberOfNodes      = 0    ;
   nodeArray          = NULL ;
   locationArray      = NULL ;
   constitutiveMatrix = NULL ;
   massMatrix         = NULL ;
   stiffnessMatrix    = NULL ;
   bodyLoadArray      = NULL ;
   gaussPointArray    = NULL ;
}


Element :: ~Element ()
   // Destructor.
{
   int i ;

   delete nodeArray ;
   delete locationArray ;
   delete massMatrix ;
   delete stiffnessMatrix ;
   delete constitutiveMatrix ;
   if (gaussPointArray) {
      for (i=0 ; i<numberOfGaussPoints ; i++)
	 delete gaussPointArray[i] ;
      delete [numberOfGaussPoints] gaussPointArray ;}
   delete bodyLoadArray ;
}


void  Element :: assembleLhsAt (TimeStep* stepN)
   // Assembles the left-hand side (stiffness matrix) of the receiver to
   // the linear system' left-hand side, at stepN.
{
   FloatMatrix* elemLhs ;
   Skyline*     systLhs ;
   IntArray*    locArray ;

   elemLhs  = this -> ComputeLhsAt(stepN) ;
   systLhs  = domain -> giveLinearSystem() -> giveLhs() ;
   locArray = this -> giveLocationArray() ;
   systLhs -> assemble(elemLhs,locArray) ;

   delete elemLhs ;
}


void  Element :: assembleRhsAt (TimeStep* stepN)
   // Assembles the right-hand side (load vector) of the receiver to
   // the linear system' right-hand side, at stepN.
{
   FloatArray* elemRhs ;
   FloatArray* systRhs ;
   IntArray*   locArray ;

   elemRhs = this -> ComputeRhsAt(stepN) ;
   if (elemRhs) {
      systRhs  = domain -> giveLinearSystem() -> giveRhs() ;
      locArray = this -> giveLocationArray() ;
      systRhs -> assemble(elemRhs,locArray) ;
      delete elemRhs ;}
}


void  Element :: assembleYourselfAt (TimeStep* stepN)
   // Assembles the contributions of the receiver to the linear system, at
   // time step stepN. This may, or may not, require assembling the receiver's
   // left-hand side.
{
#  ifdef VERBOSE
      printf ("assembling element %d\n",number) ;
#  endif

   if (stepN -> requiresNewLhs())
      this -> assembleLhsAt(stepN) ;
   this -> assembleRhsAt(stepN) ;
}


FloatArray*  Element :: ComputeBcLoadVectorAt (TimeStep* stepN)
   // Computes the load vector due to the boundary conditions acting on the
   // receiver's nodes, at stepN. Returns NULL if this array contains only
   // zeroes.
{
   FloatArray *d,*answer ;

   d = this -> ComputeVectorOfPrescribed('d',stepN) ;
   if (d -> containsOnlyZeroes())
      answer = NULL ;
   else
      answer = this -> giveStiffnessMatrix() -> Times(d) -> negated() ;

   delete d ;
   return answer ;
}


FloatArray*  Element :: ComputeBodyLoadVectorAt (TimeStep* stepN)
   // Computes numerically the load vector of the receiver due to the body
   // loads, at stepN.
{
   int         i ;
   double      dV ;
   GaussPoint* gp ;
   FloatArray  *answer,*f,*ntf ;
   FloatMatrix *n,*nt ;

   if (this -> giveBodyLoadArray() -> isEmpty())         // no loads
      return NULL ;

   else {
      f = this -> ComputeResultingBodyForceAt(stepN) ;
      if (! f)                                           // nil resultant
	 return NULL ;
      else {
	 answer = new FloatArray(0) ;
	 for (i=0 ; i<numberOfGaussPoints ; i++) {
	    gp  = gaussPointArray[i] ;
	    n   = this -> ComputeNmatrixAt(gp) ;
	    dV  = this -> computeVolumeAround(gp) ;
	    nt  = n    -> GiveTransposition() ;
	    ntf = nt   -> Times(f) -> times(dV) ;
	    answer -> add(ntf) ;
	    delete n ;
	    delete nt ;
	    delete ntf ;}
	 delete f ;
	 return answer ;}}
}


FloatMatrix*  Element :: ComputeConsistentMassMatrix ()
   // Computes numerically the consistent (full) mass matrix of the receiver.
{
   int         i ;
   double      density,dV ;
   FloatMatrix *n,*answer ;
   GaussPoint  *gp ;

   answer  = new FloatMatrix() ;
   density = this -> giveMaterial() -> give('d') ;
   for (i=0 ; i<numberOfGaussPoints ; i++) {
      gp      = gaussPointArray[i] ;
      n       = this -> ComputeNmatrixAt(gp) ;
      dV      = this -> computeVolumeAround(gp) ;
      answer -> plusProduct(n,n,density*dV) ;
      delete n ;}

   return  answer->symmetrized() ;
}


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

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


FloatArray*  Element :: ComputeLoadVectorAt (TimeStep* stepN)
   // Computes the load vector of the receiver, at stepN.
{
   FloatArray* loadVector ;
   FloatArray* bodyLoadVector = NULL ;
   FloatArray* bcLoadVector   = NULL ;

   loadVector = new FloatArray(0) ;

   bodyLoadVector = this -> ComputeBodyLoadVectorAt(stepN) ;
   if (bodyLoadVector) {
      loadVector -> add(bodyLoadVector) ;
      delete bodyLoadVector ;}

   bcLoadVector = this -> ComputeBcLoadVectorAt(stepN) ;
   if (bcLoadVector) {
      loadVector -> add(bcLoadVector) ;
      delete bcLoadVector ;}

   if (loadVector -> isNotEmpty())
      return loadVector ;
   else {
      delete loadVector ;
      return NULL ;}
}


FloatMatrix*  Element :: computeMassMatrix ()
   // Returns the lumped mass matrix of the receiver.
{
   FloatMatrix* consistentMatrix ;

   consistentMatrix = this -> ComputeConsistentMassMatrix() ;
   massMatrix       = consistentMatrix -> Lumped() ;
   delete consistentMatrix ;

   return massMatrix ;
}


FloatMatrix*  Element :: ComputeNewmarkLhsAt (TimeStep* stepN)
   // Computes the contribution of the receiver to the left-hand side of the
   // linear system, using Newmark's formula.
{
   FloatMatrix *m,*k,*lhs ;
   double      beta,dt ;

   if (stepN->giveNumber() == 0)
      lhs = this -> giveMassMatrix() -> GiveCopy() ;

   else {
      beta = domain -> giveTimeIntegrationScheme() -> giveBeta() ;
      if (beta == 0.0)
	 lhs = this -> giveMassMatrix() -> GiveCopy() ;
      else {
	 m   = this -> giveMassMatrix() ;
	 k   = this -> giveStiffnessMatrix() ;
	 dt  = stepN -> giveTimeIncrement() ;
	 lhs = k -> Times(beta*dt*dt) -> plus(m) ;}}

   return lhs ;
}


FloatArray*  Element :: ComputeNewmarkRhsAt (TimeStep* stepN)
   // Computes the contribution of the receiver to the right-hand side of the
   // linear system, using Newmark's formula.
{
   FloatMatrix *K ;
   FloatArray  *f,*d,*dPred,*rhs ;

   f = this -> ComputeLoadVectorAt(stepN) ;
   K = this -> giveStiffnessMatrix () ;

   if (stepN->giveNumber() == 0) {
      d     = this -> ComputeVectorOf ('d',stepN) ;
      rhs   = K -> Times(d->negated()) -> add(f) ;
      delete d ;}
   else {
      dPred = this -> ComputeVectorOf ('D',stepN) ;
      rhs   = K -> Times(dPred->negated()) -> add(f) ;
      delete dPred ;}

   delete f ;
   return rhs ;
}


int  Element :: computeNumberOfDofs ()
   // Returns the total number of dofs of the receiver's nodes.
{
   int i,n ;

   n = 0 ;
   for (i=1 ; i<=numberOfNodes ; i++)
      n += this -> giveNode(i) -> giveNumberOfDofs() ;
   return n ;
}


FloatArray*  Element :: ComputeResultingBodyForceAt (TimeStep* stepN)
   // Computes at stepN the resulting force due to all body loads that act
   // on the receiver. This force is used by the element for computing its
   // body load vector.
{
   int         i,n,nLoads ;
   BodyLoad*   load ;
   FloatArray  *force,*resultant ;

   resultant = new FloatArray(0) ;
   nLoads    = this -> giveBodyLoadArray() -> giveSize() ;
   for (i=1 ; i<=nLoads ; i++) {
      n     = bodyLoadArray -> at(i) ;
      load  = (BodyLoad*) domain->giveLoad(n) ;
      force = load -> ComputeForceOn(this,stepN) ;
      resultant -> add(force) ;
      delete force ;}

   if (resultant->giveSize() == 0) {
      delete resultant ;
      return NULL ;}
   else