Add a New UniaxialMaterial C++

From OpenSeesWiki
Revision as of 18:56, 11 September 2013 by Fmk (talk | contribs) (Created page with 'To add a new Uniaxial Material module using the C++ language, the developer must: # provide a new C++ subclass of the UniaxialMaterial class # provide an interface function that ...')
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

To add a new Uniaxial Material module using the C++ language, the developer must:

  1. provide a new C++ subclass of the UniaxialMaterial class
  2. provide an interface function that will be used to parse the input and create the new material.

Unlike the C and Fortran modules, no information about the state of the model is passed as arguments to the material methods. It is the responsibility of the material to okeep whatever information it needs. This information will include both parameters (information needed to define the material) and state variables or history variables (information needed by the material to remember it's current state for the computation of the stress and tangent)

NOTE: This document assumes the reader is familiar with the C++ programming language.

UniaxialMaterial Class

The Uniaxial class itself is an abstract base class. It inherits from both the Material class, which is itself a subclass of TaggedObject class and the MovableObject class. The class has a large number of methods defined in the interface, not all these methods need to be included in a new UniaxialMaterial class. The following is the minimal interface that should be considered:

The UniaxialMaterial Class:

#ifndef ElasticPPcpp_h
#define ElasticPPcpp_h

// Written: fmk                                                                                                                                
//                                                                                                                                             
// Description: This file contains the class definition for                                                                                    
// ElasticPPcpp. ElasticPPcpp provides the abstraction                                                                                         
// of an elastic perfectly plastic uniaxial material,                                                                                          
//                                                                                                                                             
// What: "@(#) ElasticPPcpp.h, revA"                                                                                                           

#include <UniaxialMaterial.h>
class UniaxialMaterial : public Material
{
  public:
    UniaxialMaterial (int tag, int classTag);
    virtual ~UniaxialMaterial();

    virtual int setTrialStrain (double strain, double strainRate =0) =0;

    virtual double getStrain (void) = 0;
    virtual double getStrainRate (void);
    virtual double getStress (void) = 0;
    virtual double getTangent (void) = 0;
    virtual double getInitialTangent (void) = 0;

    virtual int commitState (void) = 0;
    virtual int revertToLastCommit (void) = 0;
    virtual int revertToStart (void) = 0;

    virtual UniaxialMaterial *getCopy (void) = 0;

    virtual Response *setResponse (const char **argv, int argc,
                                   OPS_Stream &theOutputStream);
    virtual int getResponse (int responseID, Information &matInformation);

    virtual int sendSelf(int commitTag, Channel &theChannel);
    virtual int recvSelf(int commitTag, Channel &theChannel,
                 FEM_ObjectBroker &theBroker);

    virtual void Print(OPS_Stream &s, int flag =0);

  protected:

  private:
}

The methods with =0; are methods that you must implement for the class to link successfully with OpenSees. The other classes are optional.

The setTriaStrain() is the method called by an element when a new strain in the material is to be set. Subsequent calls to getTangent() and getStress() are to return thecorresponding tangent and stress values for that stress.

The commit() method is invoked when a solution has been obtained on the solution path. It is the responsibility of the material to be able to back track to that solution if a revertToLastCOmmit() is invoked. This will happen if the algorithm fails to find a solution on the solution path.

The getCopy() method is invoked by an element in the elements constructor. The material is to return a unique copy of itself to the element. This way different elements can use the same material type with the same properties, with each element having it's own unique copy.

The setResponse()/getResponse() typically do not have to be provided. These are the methods called by a recorder after a commit(). If you are appy with the existing responses fro a UniaxialMaterial which responds to "stress", "strain", "tangent", "stressANDstrain" you do not have to implement these methods. The example below shows them just for those ew who want additional info out of their materials.

The sendSelf()/recvSelf() methods are used in parallel processing with OpenSeesSP and when using the database command. If you don't envision using the material in these situations you can again ignore these methods. Again I am only showing them in the code for those even fewer who would do this.


The commit() method is what is called

Example - Truss2D

In the following section we will provide all necessary code to add a new 2d planar truss element into an OpenSees interpreter. To demonstrate the power of object-oriented programming, the stress-strain relationship will be provided by a UniaxialMaterial object.

Header

The header for thew new class, which we will call Truss2D is as follows:

// include directives
#include <Element.h>
#include <Matrix.h>
#include <Vector.h>

// forward declarations
class UniaxialMaterial;

class Truss2D : public Element
{
  public:
    // constructors                                                                   
    Truss2D(int tag,
            int Nd1, int Nd2,
            UniaxialMaterial &theMaterial,
            double A);

    Truss2D();

    // destructor                                                                     
    ~Truss2D();

    // initialization
    int setDomain(Domain *theDomain);

    // public methods to obtain inforrmation about dof & connectivity                 
    int getNumExternalNodes(void) const;
    const ID &getExternalNodes(void);
    Node **getNodePtrs(void);
    int getNumDOF(void);

    // public methods to set the state of the element                                 
    int commitState(void);
    int revertToLastCommit(void);
    int revertToStart(void);
    int update(void);

    // public methods to obtain stiffness  
    const Matrix &getTangentStiff(void);
    const Matrix &getInitialStiff(void);

    // public method to obtain resisting force
    const Vector &getResistingForce(void);

    // method for obtaining information specific to an element 
    void Print(OPS_Stream &s, int flag =0);
    Response *setResponse(const char **argv, int argc, OPS_Stream &s);
    int getResponse(int responseID, Information &eleInformation);

    // public methods for database/parallel processing                                                      
    int sendSelf(int commitTag, Channel &theChannel);
    int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker);
    void Print(OPS_Stream &s, int flag =0);

  protected:

  private:
    // private member functions - only available to objects of the class              
    double computeCurrentStrain(void) const;

    // private attributes - a copy for each object of the class                       
    UniaxialMaterial *theMaterial;       // pointer to a material                     
    ID  externalNodes;                   // contains the id's of end nodes            
    Matrix trans;       // hold the transformation matrix                      
    double L;           // length of truss (undeformed configuration)                                                                              
    double A;           // area of truss                                              
    Node *theNodes[2];  // node pointers                                              

    // static data - single copy for all objects of the class                         
    static Matrix trussK;   // class wide matrix for returning stiffness                            
    static Vector trussR;   // class wide vector for returning residual               
};
#endif

The header file defines the interface and variables for the class Truss2D. It defines the new class to be a sublass of the Element class. In the public interface, are two constructors and a destructor in addition to minimal set of methods we showed for the Element class. There are no protected data or methods as we do not expect this class to be subclassed. In the private section, we define one private method, computeCurrentStrain(), and we define a number of private variables and a number of static variables.

The header has a number of #include directives, one is needed for the base class and every class used as a variable in the list of data (except those that are used as pointers). For those classes that only appear as pointers in the header file (Node, UniaxialMaterial) a forward declaration is all that is needed (the include could also have been used, but using the forward declaration simplifies dependencies and reduces the amount of code that ha to be recompiled later if changes are made).

Implementation

It another file, Truss2D.cpp, we place the code that details what the constructors, destructor and methods do. In addition we provide one additional procedure OPS_Truss2D() (NOTE it has the same name as the class with an OPS_ prefix). We will go through each part of the file.

Include Directives

The first part of the file contains the list of includes. It is necessary to have an #include directive for each class and api file that is used within the .cpp file and is not included in the header.

#include "Truss2D.h"

#include <elementAPI.h>
#include <G3Globals.h>
#include <Information.h>
#include <Domain.h>
#include <Node.h>
#include <Channel.h>
#include <Message.h>
#include <FEM_ObjectBroker.h>
#include <UniaxialMaterial.h>
#include <Renderer.h>
#include <ElementResponse.h>

#include <math.h>
#include <stdlib.h>
#include <string.h>


Static Variables

Next, we initialize the static variables. For this example we are using 2 static-variables (objects shared by each Truss2D object that is created), one to return the tangent matrix and the other to return the resisting force.

// initialise the class wide variables                                                
Matrix Truss2D::trussK(4,4);
Vector Truss2D::trussR(4);
Constructors

After the list of includes, we provide the 2 constructors. The constructors are rather simple. They just initialize all the data variables defined in the header. Note it is very important to set all pointer values to 0.

The first constructor is the one most typically used. The arguments provide the elements tag, the tags of the two end nodes, the element's area and a copy of the element's material. The code in the constructor does the following:

  1. The elements tag and a 0 are passed to the Element constructor.
  2. The matreial pointer, theMaterial, is set to a copy of the material obtained from the material that is passed in the arguments.
  3. The externalNodes array is set to be an array of size 2 and it's values are set to the nodal tags of the 2 nodes.
  4. The theNodes array components are set to be 0.

It should be noted that the static variables dealing with length, transformations, and nodes are set to 0 in the constructors. They will be filled in when the setDomain() method is invoked on the object.


                                                               
Truss2D::Truss2D(int tag,
                 int Nd1, int Nd2,
                 UniaxialMaterial &theMat,
                 double a)
:Element(tag, 0),
 externalNodes(2),
 trans(1,4), L(0.0), A(a)
{
  // get a copy of the material object for our own use                                
  theMaterial = theMat.getCopy();
  if (theMaterial == 0) {
   opserr << "FATAL TrussCPP::TrussCPP() - out of memory, could not get a copy of the Material\n";
   exit(-1);
  }

  // fill in the ID containing external node info with node id's                      
  if (externalNodes.Size() != 2) {
    opserr << "FATAL TrussCPP::TrussCPP() - out of memory, could not create an ID of size 2\n";
    exit(-1);
  }

  externalNodes(0) = Nd1;
  externalNodes(1) = Nd2;

  theNodes[0] = 0;
  theNodes[1] = 0;
}

The second constructor is called when paralell processing or the database feature of the OpenSees application is used. It's pupose is to create a blank Truss2D object, that will be filled in when the recvSelf() method is invoked on the object.

Truss2D::Truss2D()
:Element(0, 0),
 theMaterial(0),
 externalNodes(2),
 trans(1,4), L(0.0), A(0.0)
{
  theNodes[0] = 0;
  theNodes[1] = 0;
}
Destructor

The we provide the destructor. In the destructor all memory that the Truss2D created or was passed to it in the constructor must be destroyed. For our example, we need to invoke the destructor on the copy of the material object.

   
Truss2D::~Truss2D()
{                                       
    if (theMaterial != 0)
        delete theMaterial;
}

===== setDomain() Initialization Method

The setDomain() method is invoked when the truss element is being added to the Domain. It is in this method that most of the private variables of the object are determined. The method returns 0 if successfull, a negative number if not. In the method we obtain pointers to the end nodes, nodal coordinates are obtained and the elements length and transformation matrix is set once the coordinates have been obtained.

void
Truss2D::setDomain(Domain *theDomain)
{
    // check Domain is not null - invoked when object removed from a domain                                              
    if (theDomain == 0) {
        return;
    }

    // first ensure nodes exist in Domain and set the node pointers                                                      
    Node *end1Ptr, *end2Ptr;
    int Nd1 = externalNodes(0);
    int Nd2 = externalNodes(1);
    end1Ptr = theDomain->getNode(Nd1);
    end2Ptr = theDomain->getNode(Nd2);
    if (end1Ptr == 0) {
      opserr << "WARNING Truss2D::setDomain() - at truss " << this->getTag() << " node " <<
        Nd1 << "  does not exist in domain\n";
        return;  // don't go any further - otherwise segemntation fault                                                  
    }
    if (end2Ptr == 0) {
      opserr << "WARNING Truss2D::setDomain() - at truss " << this->getTag() << " node " <<
        Nd2 << "  does not exist in domain\n";
        return;  // don't go any further - otherwise segemntation fault                                                  
    }
    theNodes[0] = end1Ptr;
    theNodes[1] = end2Ptr;
    // call the DomainComponent class method THIS IS VERY IMPORTANT                                                      
    this->DomainComponent::setDomain(theDomain);

    // ensure connected nodes have correct number of dof's                                                               
    int dofNd1 = end1Ptr->getNumberDOF();
    int dofNd2 = end2Ptr->getNumberDOF();
    if ((dofNd1 != 2) || (dofNd2 != 2)) {
      opserr << "Truss2D::setDomain(): 2 dof required at nodes\n";
      return;
    }

    // now determine the length & transformation matrix                                                                  
    const Vector &end1Crd = end1Ptr->getCrds();
    const Vector &end2Crd = end2Ptr->getCrds();

    double dx = end2Crd(0)-end1Crd(0);
    double dy = end2Crd(1)-end1Crd(1);

    L = sqrt(dx*dx + dy*dy);

    if (L == 0.0) {
      opserr << "WARNING Truss2D::setDomain() - Truss2D " << this->getTag() <<
        " has zero length\n";
      return;  // don't go any further - otherwise divide by 0 error                                                     
    }

    double cs = dx/L;
    double sn = dy/L;

    trans(0,0) = -cs;
    trans(0,1) = -sn;
    trans(0,2) = cs;
    trans(0,3) = sn;
}
Methods Dealing With Nodes

Next comes 4 rather simple methods that return basic information about the elements nodes. These are one line methods that should not need any explanation!

int
Truss2D::getNumExternalNodes(void) const
{
    return 2;
}

const ID &
Truss2D::getExternalNodes(void)
{
  return externalNodes;
}

Node **
Truss2D::getNodePtrs(void)
{
  return theNodes;
}

int
Truss2D::getNumDOF(void) {
    return 4;
}


Methods Dealing With Current State
int
Truss2D::commitState()
{
    return theMaterial->commitState();
}

int
Truss2D::revertToLastCommit()
{
    return theMaterial->revertToLastCommit();
}

int
Truss2D::revertToStart()
{
    return theMaterial->revertToStart();
}

int
Truss2D::update()
{
  // determine the current strain given trial displacements at nodes                  
  double strain = this->computeCurrentStrain();

  // set the strain in the materials                                                  
  theMaterial->setTrialStrain(strain);

  return 0;
}
Methods To Return Tangent Matrix

In both methods, we obtain the appropriate tangent from the material and use this to return the transformed matrix.

const Matrix &
Truss2D::getTangentStiff(void)
{
    if (L == 0.0) { // length = zero - problem in setDomain() warning message already printed                                          
        trussK.Zero();
        return trussK;
    }

    // get the current E from the material for the last updated strain                                                                 
    double E = theMaterial->getTangent();

    // form the tangent stiffness matrix                                                                                               
    trussK = trans^trans;
    trussK *= A*E/L;

    // return the matrix                                                                                                               
    return trussK;
}

const Matrix &
Truss2D::getInitialStiff(void)
{
    if (L == 0.0) { // length = zero - problem in setDomain() warning message already printed                                          
        trussK.Zero();
        return trussK;
    }

    // get the current E from the material for the last updated strain                                                                 
    double E = theMaterial->getInitialTangent();

    // form the tangent stiffness matrix                                                                                               
    trussK = trans^trans;
    trussK *= A*E/L;

    // return the matrix                                                                                                               
    return trussK;
}
Methods To Return Resisting Force

In this method we obtain the stress from the material and use this to return the transformed force vector.

const Vector &
Truss2D::getResistingForce()
{
    if (L == 0.0) { // if length == 0, problem in setDomain()                                                                          
        trussR.Zero();
        return trussR;
    }

    // want: R = Ku - Pext                                                                                                             

    // force = F * transformation                                                                                                      
    double force = A*theMaterial->getStress();
    for (int i=0; i<4; i++)
        trussR(i) = trans(0,i)*force;

    return trussR;
}
Methods Dealing With Output

Information is obtained by the user when the print command is invoked by the user and also when the user issues the recorder command. When the print command is invoked the Print method is invoked. This method simply prints information about the element, and then asks the material to do likewise.

void
Truss2D::Print(OPS_Stream &s, int flag)
{
  s << "Element: " << this->getTag();
  s << " type: Truss2D  iNode: " << externalNodes(0);
  s << " jNode: " << externalNodes(1);
  s << " Area: " << A;
  s << " \t Material: " << *theMaterial;
}

There are two methods used by the element recorders.

  1. The first method, setResponse(), is called when the recorder is created. The element informs the recorder that it can respond to a request of that type, if it cannot respond to the request it returns a 0, otherwise it returns an Response object. The response object includes a pointer to the element, an integer flag used to id the response when the getResponse() method is called, and a Vector detailing the size of the response.
  2. The second method, getReponse(), is called by the recorder when it is actually recording the information.
Response *
Truss2D::setResponse(const char **argv, int argc, OPS_Stream &output)
{
  Response *theResponse = 0;

  output.tag("ElementOutput");
  output.attr("eleType",this->getClassType());
  output.attr("eleTag",this->getTag());
  int numNodes = this->getNumExternalNodes();
  const ID &nodes = this->getExternalNodes();
  static char nodeData[32];

  for (int i=0; i<numNodes; i++) {
    sprintf(nodeData,"node%d",i+1);
    output.attr(nodeData,nodes(i));
  }

  if (strcmp(argv[0],"force") == 0 || strcmp(argv[0],"forces") == 0 ||
      strcmp(argv[0],"globalForce") == 0 || strcmp(argv[0],"globalForces") == 0) {
    const Vector &force = this->getResistingForce();
    int size = force.Size();
    for (int i=0; i<size; i++) {
      sprintf(nodeData,"P%d",i+1);
      output.tag("ResponseType",nodeData);
    }
    theResponse = new ElementResponse(this, 1, this->getResistingForce());
  }

  else if (strcmp(argv[0],"dampingForce") == 0 || strcmp(argv[0],"dampingForces") == 0) {
    const Vector &force = this->getResistingForce();
    int size = force.Size();
    for (int i=0; i<size; i++) {
      sprintf(nodeData,"P%d",i+1);
      output.tag("ResponseType",nodeData);
    }
    theResponse = new ElementResponse(this, 2, this->getResistingForce());
  } else if (strcmp(argv[0],"axialForce") ==0)
      return new ElementResponse(this, 3, 0.0);

  output.endTag();
  return theResponse;
}


int
Truss2D::getResponse(int responseID, Information &eleInfo)
{
 double strain;

  switch (responseID) {
  case -1:
    return -1;
  case 1: // global forces                                                                                               
    return eleInfo.setVector(this->getResistingForce());
  case 2:
    return eleInfo.setVector(this->getRayleighDampingForces());
  case 3:
    theMaterial->setTrialStrain(strain);
    return eleInfo.setDouble(A * theMaterial->getStress());
  default:
    return 0;
  }
}
Methods Dealing With Databases/Parallel Processing

There are two methods provided which are required is the user uses to use the database or parallel procesing features of the OpenSees applications. If neither are to be used, the developer need simply return a negative value in both methods. The idea is that the element must pack up it's information using Vector and ID objects and send it off to a Channel object. On the flip side, the receiving blank element must receive the same Vector and ID data, unpack it and set the variables.


<source lang="cpp"> int Truss2D::sendSelf(int commitTag, Channel &theChannel) {

   int res;
   // note: we don't check for dataTag == 0 for Element                                                                 
   // objects as that is taken care of in a commit by the Domain                                                        
   // object - don't want to have to do the check if sending data                                                       
   int dataTag = this->getDbTag();
   // Truss2D packs it's data into a Vector and sends this to theChannel                                                
   // along with it's dbTag and the commitTag passed in the arguments                                                   
   Vector data(5);
   data(0) = this->getTag();
   data(1) = A;
   data(2) = theMaterial->getClassTag();
   int matDbTag = theMaterial->getDbTag();
   // NOTE: we do have to ensure that the material has a database                                                       
   // tag if we are sending to a database channel.                                                                      
   if (matDbTag == 0) {
       matDbTag = theChannel.getDbTag();
       if (matDbTag != 0)
           theMaterial->setDbTag(matDbTag);
   }
   data(3) = matDbTag;
   res = theChannel.sendVector(dataTag, commitTag, data);
   if (res < 0) {
     opserr << "WARNING Truss2D::sendSelf() - failed to send Vector\n";
     return -1;
   }
   // Truss2D then sends the tags of it's two end nodes                                                                 
   res = theChannel.sendID(dataTag, commitTag, externalNodes);
   if (res < 0) {
     opserr << "WARNING Truss2D::sendSelf() - failed to send ID\n";
     return -2;
   }
   // finally Truss2D asks it's material object to send itself                                                          
   res = theMaterial->sendSelf(commitTag, theChannel);
   if (res < 0) {
     opserr << "WARNING Truss2D::sendSelf() - failed to send the Material\n";
     return -3;
   }
  return 0;

}

int Truss2D::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker) {

   int res;
   int dataTag = this->getDbTag();
   // Truss2D creates a Vector, receives the Vector and then sets the                                                   
   // internal data with the data in the Vector                                                                         
   Vector data(5);
   res = theChannel.recvVector(dataTag, commitTag, data);
   if (res < 0) {
     opserr << "WARNING Truss2D::recvSelf() - failed to receive Vector\n";
     return -1;
   }
   this->setTag((int)data(0));
   A = data(1);
   // Truss2D now receives the tags of it's two external nodes                                                          
   res = theChannel.recvID(dataTag, commitTag, externalNodes);
   if (res < 0) {
     opserr << "WARNING Truss2D::recvSelf() - failed to receive ID\n";
     return -2;
   }
   // we create a material object of the correct type,                                                                  
   // sets its database tag and asks this new object to recveive itself.                                                
   int matClass = data(2);
   int matDb = data(3);
   theMaterial = theBroker.getNewUniaxialMaterial(matClass);
   if (theMaterial == 0) {
     opserr << "WARNING Truss2D::recvSelf() - failed to create a Material\n";
     return -3;
   }
   // we set the dbTag before we receive the material  - this is important                                              
   theMaterial->setDbTag(matDb);
   res = theMaterial->recvSelf(commitTag, theChannel, theBroker);
   if (res < 0) {
     opserr << "WARNING Truss2D::recvSelf() - failed to receive the Material\n";
     return -3;
   }
   return 0;

}