Add a New UniaxialMaterial C++

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 - ElasticPPecpp

In the following section we will provide all necessary code to add a new elastic perfectly plastic material into an OpenSees interpreter.

Header

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

#ifndef ElasticPPcpp_h
#define ElasticPPcpp_h

#include <UniaxialMaterial.h>

class ElasticPPcpp : public UniaxialMaterial
{
  public:
    ElasticPPcpp(int tag, double E, double eyp);
    ElasticPPcpp();

    ~ElasticPPcpp();

    int setTrialStrain(double strain, double strainRate = 0.0);
    double getStrain(void);
    double getStress(void);
    double getTangent(void);


    double getInitialTangent(void) {return E;};

    int commitState(void);
    int revertToLastCommit(void);
    int revertToStart(void);

    UniaxialMaterial *getCopy(void);

    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:
    double fyp, fyn;    // positive and negative yield stress                                                                                  
    double ezero;       // initial strain                                                                                                      
    double E;           // elastic modulus                                                                                                     
    double ep;          // plastic strain at last commit                                                                                       

    double trialStrain; // trial strain                                                                                                        
    double trialStress;      // current trial stress                                                                                           
    double trialTangent;     // current trial tangent                                                                                          
    double commitStrain;     // last commited strain                                                                                           
    double commitStress;     // last commited stress                                                                                           
    double commitTangent;    // last committed  tangent

The header file defines the interface and variables for the class ElasticPPcpp. It defines the new class to be a sublass of the UniaxialMaterial class. In the public interface, are two constructors and a destructor in addition to minimal set of methods we showed for the UniaxialMaterial class. There are no protected data or methods as we do not expect this class to be subclassed. In the private section we define a number of private variables and a number of variables. Some of these are parameter variable which do not change with each commit, e.g. E, and some state variable which do change, e.g. ep, fyp, and fyn.

Implementation

It another file, ElasticPPcpp.cpp, we place the code that details what the constructors, destructor and methods do. In addition we provide one additional procedure OPS_ElasticPPcpp() (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 "ElasticPPcpp.h"

#include <elementAPI.h>
#include <Vector.h>
#include <Channel.h>
#include <math.h>
#include <float.h>

Static Variables

Next, we initialize the static variables. For this example we are using 1 static-variables to keep track of the number of times the external procedure to parse and create such an object is called.

static int numElasticPPcpp = 0;

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 if you use pointers in your class. We will use none here.

The first constructor is the one most typically used. The arguments provide the materials tag, youngs modulus and initial yield point strain values. material. The code in the constructor simply computes the positive and negative yield stress based on the input provided.

                                                               
ElasticPPcpp::ElasticPPcpp(int tag, double e, double eyp)
:UniaxialMaterial(tag, 0),
 ezero(0.0), E(e), ep(0.0),
 trialStrain(0.0), trialStress(0.0), trialTangent(E),
 commitStrain(0.0), commitStress(0.0), commitTangent(E)
{
  fyp = E*eyp;
  fyn = -fyp;
}

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

ElasticPPcpp::ElasticPPcpp()
:UniaxialMaterial(0, 0),
 fyp(0.0), fyn(0.0), ezero(0.0), E(0.0), ep(0.0),
 trialStrain(0.0), trialStress(0.0), trialTangent(E),
 commitStrain(0.0), commitStress(0.0), commitTangent(E)
{

}

Destructor

The we provide the destructor. In the destructor all memory that the the object created or was passed to it in the constructor must be destroyed. For this example we have no such memory. We could have left the destructor out entirely. Hoowever, it is good practice to leave it in your source code.

   
ElasticPPcpp::~ElasticPPcpp()
{
  // does nothing                                                                                                                              
}

===== setTrialStrain() Method

This, as mentioned, is the method called when the element has computed a nw strain for the element. The element will make subsequent calls to getStress() and getTangent() to obtain new values of these for the new strain. This is typically the most complicated method to write and to determine the theory for before you even write the code. ALl subsequent methods are trivial.

int
ElasticPPcpp::setTrialStrain(double strain, double strainRate)
{
    if (fabs(trialStrain - strain) < DBL_EPSILON)
      return 0;

    trialStrain = strain;

    double sigtrial;    // trial stress                                                                                                        
    double f;           // yield function                                                                                                      

    // compute trial stress                                                                                                                    
    sigtrial = E * ( trialStrain - ezero - ep );

    //sigtrial  = E * trialStrain;                                                                                                             
    //sigtrial -= E * ezero;                                                                                                                   
    //sigtrial -= E *  ep;                                                                                                                     

    // evaluate yield function                                                                                                                 
    if ( sigtrial >= 0.0 )
        f =  sigtrial - fyp;
    else
        f = -sigtrial + fyn;

    double fYieldSurface = - E * DBL_EPSILON;
    if ( f <= fYieldSurface ) {

      // elastic                                                                                                                               
      trialStress = sigtrial;
      trialTangent = E;

    } else {

      // plastic                                                                                                                               
      if ( sigtrial > 0.0 ) {
        trialStress = fyp;
      } else {
        trialStress = fyn;
      }

      trialTangent = 0.0;
    }

    return 0;
}

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.

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;
}


OPS_Export void *
OPS_ElasticPPcpp()
{
  // print out some KUDO's                                                                                                                     
  if (numElasticPPcpp == 0) {
    opserr << "ElasticPPcpp unaxial material - Written by fmk UC Berkeley Copyright 2008 - Use at your Own Peril\n";
    numElasticPPcpp =1;
  }

  // Pointer to a uniaxial material that will be returned                                                                                      
  UniaxialMaterial *theMaterial = 0;

  //                                                                                                                                           
  // parse the input line for the material parameters                                                                                          
  //                                                                                                                                           

  int    iData[1];
  double dData[2];
  int numData;
  numData = 1;
  if (OPS_GetIntInput(&numData, iData) != 0) {
    opserr << "WARNING invalid uniaxialMaterial ElasticPP tag" << endln;
    return 0;
  }

  numData = 2;
  if (OPS_GetDoubleInput(&numData, dData) != 0) {
    opserr << "WARNING invalid E & ep\n";
    return 0;
  }

  //                                                                                                                                           
  // create a new material                                                                                                                     
  //                                                                                                                                           

  theMaterial = new ElasticPPcpp(iData[0], dData[0], dData[1]);

  if (theMaterial == 0) {
    opserr << "WARNING could not create uniaxialMaterial of type ElasticPPCpp\n";
    return 0;
  }

  // return the material                                                                                                                       
  return theMaterial;
}