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. setTrialStrain() is invoked as the solution algorithm tries a number of trial solution steps as it goes from one commited solution to the next on the solution path.

The commitState() method is invoked when a trial solution has been determined to be 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                                                                                                                              
}

getCopy() Method

This is the method called by each element or section to get unique copies of a material.

  
UniaxialMaterial *
ElasticPPcpp::getCopy(void)
{
  ElasticPPcpp *theCopy =
    new ElasticPPcpp(this->getTag(),E,fyp/E);
  theCopy->ep = this->ep;

  return theCopy;
}

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

Trivial Methods

Next comes 3 rather simple methods that return basic information computed in the setTrialStrain(). You do of course have the option to ignore the setTrialStrain() method and compute the stress and tangent quantities again in the interests of saving memory.

double
ElasticPPcpp::getStrain(void)
{
  return trialStrain;
}

double
ElasticPPcpp::getStress(void)
{
  return trialStress;
}


double
ElasticPPcpp::getTangent(void)
{
  return trialTangent;
}

Methods Dealing With Current State

As mentioned, when the algorithm finds a solution state as it goes from one converged solution to the next. As it attempts to find these solutions it goes through a number of trial steps (each setTrialStrain() is invoked in each of these steps). Once it finds a trial step that is on the solution path it will stop and invoke commitState() on the material. Any state variables that the material uses needs to be updated at this time. Should the algorithm fail to find a solution it may return to the last converged step or indeed the start. You the developer must provide code so that your mateial can indeed go back to these states and report correct getTangent() and getStress() values for subsequent analysis atte,pts.

int
ElasticPPcpp::commitState(void)
{
    double sigtrial;    // trial stress                                                                                                        
    double f;           // yield function                                                                                                      

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

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

    double fYieldSurface = - E * DBL_EPSILON;
    if ( f > fYieldSurface ) {
      // plastic                                                                                                                               
      if ( sigtrial > 0.0 ) {
        ep += f / E;
      } else {
        ep -= f / E;
      }
    }

    commitStrain = trialStrain;
    commitTangent=trialTangent;
    commitStress = trialStress;

    return 0;
}

int
ElasticPPcpp::revertToLastCommit(void)
{
  trialStrain = commitStrain;
  trialTangent = commitTangent;
  trialStress = commitStress;

  return 0;
}


int
ElasticPPcpp::revertToStart(void)
{
  trialStrain = commitStrain = 0.0;
  trialTangent = commitTangent = E;
  trialStress = commitStress = 0.0;

  ep = 0.0;

  return 0;
}

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
ElasticPPcpp::Print(OPS_Stream &s, int flag)
{
  s << "ElasticPPcpp tag: " << this->getTag() << endln;
  s << "  E: " << E << endln;
  s << "  ep: " << ep << endln;
  s << "  stress: " << trialStress << " tangent: " << trialTangent << endln;
}

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.

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 material 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
ElasticPPcpp::sendSelf(int cTag, Channel &theChannel)
{
  int res = 0;
  static Vector data(9);
  data(0) = this->getTag();
  data(1) = ep;
  data(2) = E;
  data(3) = ezero;
  data(4) = fyp;
  data(5) = fyn;
  data(6) = commitStrain;
  data(7) = commitStress;
  data(8) = commitTangent;

  res = theChannel.sendVector(this->getDbTag(), cTag, data);
  if (res < 0)
    opserr << "ElasticPPcpp::sendSelf() - failed to send data\n";

  return res;
}

int
ElasticPPcpp::recvSelf(int cTag, Channel &theChannel,
                                 FEM_ObjectBroker &theBroker)
{
  int res = 0;
  static Vector data(9);
  res = theChannel.recvVector(this->getDbTag(), cTag, data);
  if (res < 0)
    opserr << "ElasticPPcpp::recvSelf() - failed to recv data\n";
  else {
    this->setTag(data(0));
    ep    = data(1);
    E     = data(2);
    ezero = data(3);
    fyp   = data(4);
    fyn   = data(5);
    commitStrain=data(6);
    commitStress=data(7);
    commitTangent=data(8);
    trialStrain = commitStrain;
    trialTangent = commitTangent;
    trialStress = commitStress;
  }

  return res;
}

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